小技巧-获取gff3格式的内含子注释信息

用python从基因组注释文件（gff3格式）整理出基因编码区、起始密码子上游几KB和内含子信息的gff3文件。

前言

多数情况下，参考基因组的注释文件不会显示内含子的位置，因为一个基因内的两个CDS间就是一个内含子。有时候做一些个性化分析就是需要仅包含内含子或其它基因结构信息的bed或gff文件。比如我现在想知道我已鉴定到的转座子在基因结构上的分布，需要通过各个基因结构的gff或bed文件和转座子位置的交集获得我想要的信息，在这种情况下就需要自己写脚本或者找现成的轮子了...

我的需求

不管是基因的上游区域还是内含子，后面必须标注出这个区域的类型（内含子？编码区？）以及对应的基因或转录本ID，形如 chr start end type geneID。

解决方案

找现成的轮子

GitHub：https://github.com/irusri/Extract-intron-from-gff3 图1 作者给出了3种解决方案，第一种我不熟悉；第二种作者给出了他写的脚本，我感觉依赖的库太多所以没有用；第三种分享了一个在线工具，但我打不开这个链接...最后以失败告终。

后来在公众号上看到一个伙伴遇到了和我同样的问题，推送在https://mp.weixin.qq.com/s/GWlba4tfyfE9VBA5m5kVTg，号主分享的解决方案是

显然也挺麻烦....

自己写脚本

搜索现成的轮子花了我不少的时间，最后还是自己动手，丰衣足食吧。gff文件内容如下

1.编码区和起始密码子上游的gff文件

对于编码区而言，判断第三列为CDS的行就行了；mRNA的起始位置往前推几Kb就是起始密码子的上游区域（因为在该注释文件中mRNA的起始位置和第一个CDS起始位置一致）。代码如下

gff = open(r'D:\doc\GrapeStudy\grape-genome\Newversion\代表性转录本的gff3文件.gff3')
upstream = open(r'D:\doc\GrapeStudy\grape-genome\Newversion\upstream10kb.gff3', 'w+')
cds = open(r'D:\doc\GrapeStudy\grape-genome\Newversion\cds.gff3', 'w+')
intron = open(r'D:\doc\GrapeStudy\grape-genome\Newversion\intron.gff3', 'w+')
aList = []

for line in gff:
    line = line.strip().split('\t')
    aList.append(line)
    # 上游10kb
    if line[2] == 'mRNA':
        up_start = int(line[3]) - 10000
        up_end = line[3]
        line[3] = str(up_start)
        line[4] = str(up_end)
        line[2] = 'upstream'
        upresult = '\t'.join(line) + '\n'
        upstream.write(upresult)
    # CDS区
    elif line[2] == 'CDS':
        cdsresult = '\t'.join(line) + '\n'
        cds.write(cdsresult)

起始密码子上游

CDS区 

2.内含子的gff文件

实现思路：

1. 遍历文件的每一行。并将文本转为列表；
2. 若该行的第三列为gene和mRNA则pass；
3. 若该行的第三列为CDS，前一行为mRNA，则处理并写入文件；
4. 若该行第三列为CDS。下一行仍为CDS，则处理并写入文件；
5. 若该行第三列为CDS，下一行为gene，则跳过；

代码如下

# 获得内含子区
for i in range(0, len(aList)-1):
    if aList[i][2] == 'gene':
        pass
    elif aList[i][2] == 'mRNA':
        pass
    elif aList[i-1][2] == 'mRNA' and aList[i][2] == 'CDS':
        intron_start = aList[i][4]
        intron_end = aList[i+1][3]
        aList[i][2] = 'intron'
        aList[i][3] = str(intron_start)
        aList[i][4] = str(intron_end)
        intronresult = '\t'.join(aList[i]) + '\n'
        intron.write(intronresult)
    elif aList[i][2] == 'CDS' and aList[i+1][2] == 'CDS':
        intron_start = aList[i][4]
        intron_end = aList[i+1][3]
        aList[i][2] = 'intron'
        aList[i][3] = str(intron_start)
        aList[i][4] = str(intron_end)
        intronresult = '\t'.join(aList[i]) + '\n'
        intron.write(intronresult.replace('CDS', 'intron'))
    elif aList[i][2] == 'CDS' and aList[i+1][2] == 'gene':
        continue

gff.close()
cds.close()
upstream.close()
intron.close()

只需1秒，60+M的文件就被处理完毕

总结

比较简单的小功能如果找不到现成的比较好的轮子，可以试着自己造一个，顺便提升下编程能力和思维。如果脑子里只想着找现成的东西来用，自己曾经“引以为傲”的编程能力很可能会被自己的懒惰侵蚀掉。