脚本 | Shell | 提取每个转录本exon坐标

网上找了一些方法，但都不理想。

• 以拟南芥gff文件为例，查看$3，没有exon，所以没办法直接提取

$ grep -v "#" Arabidopsis_thaliana.gff | awk '{print $3}' | sort | uniq -c
 197160 CDS
  34621 five_prime_UTR
  27416 gene
  35386 mRNA
  30634 three_prime_UTR

• 所以，如何将UTR，CDS坐标整合成exon，就是关键
  经常使用gffread提取序列，同时我也注意到这个软件对exon也有一些参数，出于直觉...

$ gffread -h
...
 --force-exons: make sure that the lowest level GFF features are printed as
      "exon" features
...
# lowest ...巴拉巴拉，可能就是针对一个转录本多个five_prime_UTR、three_prime_UTR的，没有进行验证
$ gffread Arabidopsis_thaliana.gff --force-exons -o Arabidopsis_thaliana.exon.gff
$ grep -v "#" Arabidopsis_thaliana.exon.gff | awk '{print $3}' | sort | uniq -c
 197160 CDS
 207465 exon
  35386 mRNA

$ less -S Arabidopsis_thaliana.gff
Chr1    phytozomev10    mRNA    3631    5899    .       +       .       ID=AT1G01010.1.TAIR10;geneID=AT1G01010.TAIR10;gene_name=AT1G01010
Chr1    phytozomev10    exon    3631    3913    .       +       .       Parent=AT1G01010.1.TAIR10
Chr1    phytozomev10    exon    3996    4276    .       +       .       Parent=AT1G01010.1.TAIR10
Chr1    phytozomev10    exon    4486    4605    .       +       .       Parent=AT1G01010.1.TAIR10
Chr1    phytozomev10    exon    4706    5095    .       +       .       Parent=AT1G01010.1.TAIR10
Chr1    phytozomev10    exon    5174    5326    .       +       .       Parent=AT1G01010.1.TAIR10
Chr1    phytozomev10    exon    5439    5899    .       +       .       Parent=AT1G01010.1.TAIR10
Chr1    phytozomev10    CDS     3760    3913    .       +       0       Parent=AT1G01010.1.TAIR10
Chr1    phytozomev10    CDS     3996    4276    .       +       2       Parent=AT1G01010.1.TAIR10
Chr1    phytozomev10    CDS     4486    4605    .       +       0       Parent=AT1G01010.1.TAIR10
Chr1    phytozomev10    CDS     4706    5095    .       +       0       Parent=AT1G01010.1.TAIR10
Chr1    phytozomev10    CDS     5174    5326    .       +       0       Parent=AT1G01010.1.TAIR10
Chr1    phytozomev10    CDS     5439    5630    .       +       0       Parent=AT1G01010.1.TAIR10
...

• 那获取坐标就不用说了，直接长度统计

$ awk '$3 ~ /exon/{gsub(/.*Parent=/,"",$9);gsub(/;.*/,"",$9);if(a[$9]=="")a[$9]=$5-$4+1 ;else a[$9]+=$5-$4+1}END{for(i in a)print i"\t"a[i]}' Arabidopsis_thaliana.exon.gff | sort -k1 | les
AT1G01010.1.TAIR10      1688
AT1G01020.1.TAIR10      1623
AT1G01020.2.TAIR10      1085
AT1G01030.1.TAIR10      1905
AT1G01040.1.TAIR10      6251
AT1G01040.2.TAIR10      5877
AT1G01050.1.TAIR10      976
...

• 获得每个基因中最长exon的转录本
  需要结合 日常 | 2021-11-26 | 批量从*.gff 提取基因和转录本的对应关系 的结果geneid_transid.tsv

awk '$3 ~ /exon/{gsub(/.*Parent=/,"",$9);gsub(/;.*/,"",$9);if(a[$9]=="")a[$9]=$5-$4+1 ;else a[$9]+=$5-$4+1}END{for(i in a)print i"\t"a[i]}' Arabidopsis_thaliana.exon.gff  | \
awk 'ARGIND==1{a[$1]=$2}ARGIND==2{print $0"\t"a[$2]}' - geneid_transid.tsv | sort -k 1,1 -k 3,3nr | awk '{if(a[$1]==""){a[$1]=$2;print $1"\t"$2"\t"$3} else next}'

AT1G01010.TAIR10        AT1G01010.1.TAIR10      1688
AT1G01020.TAIR10        AT1G01020.1.TAIR10      1623
AT1G01030.TAIR10        AT1G01030.1.TAIR10      1905
AT1G01040.TAIR10        AT1G01040.1.TAIR10      6251
AT1G01050.TAIR10        AT1G01050.1.TAIR10      976
...

从基因AT1G01040可以观察到正确性

最后的代码，其实也可以用来“土方法”获取最长转录本。