脚本 | Shell | 批量从*.gff 提取基因和转录本的对

从日常的体验，gff文件五花八门：

1. 第3列没有`gene`；
2. 有的`gene`没有`CDS`；# 见于转录组分析中更新了注释的`gff`文件
3. 第一列在`genome.fa`中没有对应的染色体号；# 多是 `.*_ERROP.*`，需要删除后进行`CDS`提取等操作
4. 最后一个字符`^M`；# windows下的换行符，可以用 `sed 's/\r$//'`或 `dos2unix`去除
5. 同一转录本多个`five_prime_UTR`信息，不知道算不算问题；
6. 多了去了，想起来补充

• 欣赏一下我今天碰见的.gff文件，精简！

$ grep -v "#" test1.gff | awk '{print $3}' | sort | uniq -c
 174760 CDS
  32886 mRNA

• 下面一行代码，就是根据$3 == mRNA行，第9列 ID=.* 和 Parent=.*；如果没有Parent，输出两列ID

$ grep -v "#" test1.gff | awk '$3 ~ /mRNA/{print $9}' | \
awk -F ";|=" '{for(i=1;i<=NF/2;i++)a[$(2*i-1)]=$(2*i);if(a["Parent"]!="")print a["Parent"]"\t"a["ID"]; else print a["ID"]"\t"a["ID"]}' 

Distr1S00001    Distr1S00001
Distr1S00002    Distr1S00002
Distr1S00003    Distr1S00003
Distr1S00004    Distr1S00004
Distr1S00005    Distr1S00005
Distr1S00006    Distr1S00006
...

• 换拟南芥的gff

AT1G01010.TAIR10        AT1G01010.1.TAIR10
AT1G01020.TAIR10        AT1G01020.1.TAIR10
AT1G01020.TAIR10        AT1G01020.2.TAIR10
AT1G01030.TAIR10        AT1G01030.1.TAIR10
AT1G01040.TAIR10        AT1G01040.2.TAIR10
AT1G01040.TAIR10        AT1G01040.1.TAIR10
AT1G01050.TAIR10        AT1G01050.1.TAIR10
...

其实可以用python去构建gene，transcript，CDS，甚至exon之间的对应关系，基于我的python水平和gff文件的“多样性”，要耗费的心力和时间......虽然解决问题的过程可以提升代码水平，但暂时还是算了吧。