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

转载请保留出处

从细菌GFF文件提取CDS序列并转换为氨基酸序列

 

 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):               #获取CDS信息
14         line = self.pa.split(line)
15         CDS = (line[0], line[3], line[4], line[6], line[7])  #这里字符串分割有的文件是会出问题的，所以要看文件格式而定
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:                         #AA.txt 为密码子表
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                 sequence = sequence[int(i[4]):]                         # i[4] 表示密码子开始位置
51                 if i[3] == '-':
52                     sequence = self.N2M(sequence)
53                 #w.write(i[0] + '\n' + sequence + '\n')
54                 #后面是一堆正则，主要是对序列做注释的，看文件格式而定
55                 
56                 s1 = self.pa.split(line)
57                 p1 = re.compile(r'ID=(.*?);.*?Dbxref=(.*?);.*?Name=(.*?);.*?gbkey=(.*?);.*?product=(.*?);.*?protein_id=(.*?);')
58                 p2 = re.compile(r'.*?product=(.*?);.*?protein_id=(.*?);')
59                 m1 = re.findall(p1,line)
60                 m2 = re.findall(p2,line)
61                 s = '>'+s1[0]+'_'+m1[0][0]+'\tName='+m1[0][2]+'\tdbxref='+m1[0][1]+'\tprotein='+m1[0][4]+'\tprotein_id='+m1[0][5]+'\tgbkey='+m1[0][3]
62                 w.write(s + '\n' + sequence + '\n')
63                 p = '>' + s1[0]+'\tproduct:'+m2[0][0]+'\tproduct_id:'+m2[0][1]
64                 AA = self.Seq2AA(sequence, hash_AA)
65                 r.write(p + '\n' + AA + '\n')
66         w.close()
67         r.close()
68 
69 if __name__=='__main__':
70     fna = 'GCA_000160075.2_ASM16007v2_genomic.fna'
71     gff = 'GCA_000160075.2_ASM16007v2_genomic.gff'
72     m = CDS2AA(fna,gff)
73     m.CDS2AA()

 

 

 

出现的一些问题我会慢慢完善。后面的有意思作业题目会陆续上传。

posted @ 2017-11-26 22:34  orange1002  阅读(5032)  评论(1)  编辑  收藏  举报