网上找了一些方法，但都不理想。
- 以拟南芥
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可以观察到正确性
最后的代码，其实也可以用来“土方法”获取最长转录本。
