python解析gff文件中的转录本

1.下载基因组注释文件，选择对应的版本： ftp://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/ARCHIVE/BUILD.37.3/GFF/ 

2.GTF 为General Transfer Format ，熟悉格式 http://www.huoyunjn.com/wuliuxinwen/2/33709819.htm。

第三列 feature  - 后面start和end之间区域代表的特征，如果此区域是基因，则此处为gene，如果是外显子，则为exon，如果是转录本，则为transcript，如果是非编码RNA则为lncRNA，如果是重复序列，则为TE，等等，主要表明这一块区域的特征。

3.每一个transcript对应的exon，所有长度加起来就是这个转录本的长度。与这个transcript后面的两列相减是有差别的。

4.用python 字典来统计每个转录本的长度。

            
              import pandas as pd
import pdb
df = pd.read_table(r'C:\Users\guosheng\Desktop\out.gff',sep = '\t',header= None)
out=open('./out.txt','a')
df =df[df.iloc[:,2].str.contains('exon')]  #提取第三列为exon的行
df['diff'] =df.iloc[:,4]-df.iloc[:,3]+1 #每个外显子的长度
name =  list(df.iloc[:,2])  #把data.frame中的一列转换为list
des =list(df.iloc[:,8])
length = list(df['diff'])   
dic ={}      
for index,value in enumerate(name):
    key=des[index].split(';')[-1].split('=')[-1] #获取每个转录本的名字
    old=0
    new=length[index]
    if dic.has_key(key):       #判断这个key是否在原有的字典中
        old=dic[key]
        del(dic[key])
    dic[key]=int(old)+int(new)
#print dic
for tran in dic:
    out.write(tran+'\t'+str(dic[tran])+'\n')
out.flush()
out.close()

            
          

5.后续找出每个基因的所有转录本，用heapq库找出最长的一个。库用法https://blog.csdn.net/Cassiel60/article/details/88344137

同样是解析这个文件，可以看出文件中的id是根据第三列进行编号的，没有实际意义，只是可以看出共有多少个gene、exon、cds等。不过在第三列为gene时，Name=和Dbxref=GeneID: 与第四列为exon时，Dbxref=GeneID:和transcript_id=进行基因与转录本的正确匹配。可以在上面代码中的字典加入Dbxref=，字典中一键对应多个值。

6.得到的两个文件进行merge，就可以得到基因，转录本，长度的文件了。