從細菌GFF文件提取CDS序列並轉換為氨基酸序列

最近在上生物信息學原理，打算記錄一些課上的作業。第一次作業：如題。

基本思路：

　　　　　　1.從GFF中讀取CDS的起始終止位置以及正負鏈信息。GFF格式見 http://blog.sina.com.cn/s/blog_8a4f556e0102yd3l.html.

　　　　　　2.利用起始/終止位置等信息從FNA文件中提取CDS序列。FNA格式見 http://boyun.sh.cn/bio/?p=1192.

　　　　　　3.利用CDS序列及密碼子表得到FAA文件並輸出。

註意：最需要註意的一點是：當GFF中CDS位於負鏈時，需要進行鹼基互補配對，即反向互補(5'到3'配3'到5')。

下麵給出python代碼。python3.6

以下問題有待解決：得到的FNA文件和FAA文件的註釋信息不全，並且沒有提供命令行參數，交完作業再補充。

 1 #bioinformatics homework
 2 import re
 3 class CDS2AA():
 4     pa = re.compile(r'\s+')
 5     Pa = re.compile(r'[TCAG]TG')                 # 細菌起始密碼子NTG
 6     def __init__(self,fna,gff):
 7         self.fna = fna
 8         self.gff = gff
 9     def N2M(self,sequence):
10         hash = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'}
11         sequence = ''.join([hash[i] for i in sequence])     #正負鏈轉換
12         return sequence[::-1]
13     def Get_CDS_index(self,line):
14         line = self.pa.split(line)
15         CDS = (line[0], line[3], line[4], line[6])
16         return CDS
17     def Seq2AA(self,sequence,hash):
18         AA = hash[sequence[:3]]
19         if self.Pa.match(sequence[:3]):
20             AA = 'M'                                 #起始密碼子
21         for i in range(3, len(sequence) - 3, 3):
22             AA += hash[sequence[i:i + 3]]
23         return AA
24     def CDS2AA(self):
25         fn = open(self.fna, 'r')
26         gf = open(self.gff,'r')
27         r = open('AA_sequence.txt', 'w')
28         w = open('CDS.txt', 'w')
29         hash_AA = {}  # AA hash,sequence2AA
30         with open('AA.txt', 'r') as f:              
31             for line in f:
32                 line = line.strip()
33                 if line:
34                     line = self.pa.split(line)
35                     hash_AA[line[0]] = line[1]      #AA hash
36         hash_CDS = {}  # CDS hash,CDS2sequence
37         for line in fn:
38             line = line.strip()
39             if line.startswith('>'):
40                 A = self.pa.split(line)[0].replace('>', '')
41                 hash_CDS[A] = ''
42             else:
43                 hash_CDS[A] += line
44         fn.close()
45         for line in gf:
46             line = line.strip()
47             if 'CDS' in line:
48                 i = self.Get_CDS_index(line)
49                 sequence = hash_CDS[i[0]][int(i[1]) - 1:int(i[2])]
50                 if i[3] == '-':
51                     sequence = self.N2M(sequence)
52                 #w.write(i[0] + '\n' + sequence + '\n')
53                 s = self.pa.split(line)
54                 p = re.compile(r'ID=(.*?);.*?Dbxref=(.*?);.*?Name=(.*?);.*?gbkey=(.*?);.*?product=(.*?);.*?protein_id=(.*?);')
55                 m = re.findall(p,line)
56                 s = s[0]+'_'+m[0][0]+m[0][2]+'\tdbxref='+m[0][1]+'\tprotein='+m[0][4]+'\tprotein_id='+m[0][5]+'\tgbkey='+m[0][3]
57                 w.write(s + '\n' + sequence + '\n')
58                 AA = self.Seq2AA(sequence, hash_AA)
59                 r.write(i[0] + '\n' + AA + '\n')
60         w.close()
61         r.close()
62 
63 if __name__=='__main__':
64     fna = 'GCA_000160075.2_ASM16007v2_genomic.fna'
65     gff = 'GCA_000160075.2_ASM16007v2_genomic.gff'
66     m = CDS2AA(fna,gff)
67     m.CDS2AA()

一些其他的問題我會在交作業前完善。後面的有意思作業題目會陸續上傳。