从细菌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()

一些其他的问题我会在交作业前完善。后面的有意思作业题目会陆续上传。