1.下載基因組注釋文件，選擇對(duì)應(yīng)的版本：?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之間區(qū)域代表的特征，如果此區(qū)域是基因，則此處為gene，如果是外顯子，則為exon，如果是轉(zhuǎn)錄本，則為transcript，如果是非編碼RNA則為lncRNA，如果是重復(fù)序列，則為TE，等等，主要表明這一塊區(qū)域的特征。

3.每一個(gè)transcript對(duì)應(yīng)的exon，所有長(zhǎng)度加起來就是這個(gè)轉(zhuǎn)錄本的長(zhǎng)度。與這個(gè)transcript后面的兩列相減是有差別的。

4.用python 字典來統(tǒng)計(jì)每個(gè)轉(zhuǎn)錄本的長(zhǎng)度。

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 #每個(gè)外顯子的長(zhǎng)度 name = list(df.iloc[:,2]) #把data.frame中的一列轉(zhuǎn)換為list des =list(df.iloc[:,8]) length = list(df['diff']) dic ={} for index,value in enumerate(name):key=des[index].split(';')[-1].split('=')[-1] #獲取每個(gè)轉(zhuǎn)錄本的名字old=0new=length[index]if dic.has_key(key): #判斷這個(gè)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.后續(xù)找出每個(gè)基因的所有轉(zhuǎn)錄本，用heapq庫找出最長(zhǎng)的一個(gè)。庫用法https://blog.csdn.net/Cassiel60/article/details/88344137

同樣是解析這個(gè)文件，可以看出文件中的id是根據(jù)第三列進(jìn)行編號(hào)的，沒有實(shí)際意義，只是可以看出共有多少個(gè)gene、exon、cds等。不過在第三列為gene時(shí)，Name=和Dbxref=GeneID: 與第四列為exon時(shí)，Dbxref=GeneID:和transcript_id=進(jìn)行基因與轉(zhuǎn)錄本的正確匹配。可以在上面代碼中的字典加入Dbxref=，字典中一鍵對(duì)應(yīng)多個(gè)值。

6.得到的兩個(gè)文件進(jìn)行merge，就可以得到基因，轉(zhuǎn)錄本，長(zhǎng)度的文件了。

?

總結(jié)

以上是生活随笔為你收集整理的python解析gff文件中的转录本的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。