最近在上生物信息学原理，打算记录一些课上的作业。第一次作业：如题。
基本思路：
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
转载请保留出处
#bioinformatics homework
import re
class CDS2AA():
pa = re.compile(r'\s+')
Pa = re.compile(r'[TCAG]TG') # 细菌起始密码子NTG
def __init__(self,fna,gff):
self.fna = fna
self.gff = gff
def N2M(self,sequence):
hash = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'}
sequence = ''.join([hash[i] for i in sequence]) #正负链转换
return sequence[::-1]
def Get_CDS_index(self,line): #获取CDS信息
line = self.pa.split(line)
CDS = (line[0], line[3], line[4], line[6], line[7]) #这里字符串分割有的文件是会出问题的，所以要看文件格式而定
return CDS
def Seq2AA(self,sequence,hash):
AA = hash[sequence[:3]]
if self.Pa.match(sequence[:3]):
AA = 'M' #起始密码子
for i in range(3, len(sequence) - 3, 3):
AA += hash[sequence[i:i + 3]]
return AA
def CDS2AA(self):
fn = open(self.fna, 'r')
gf = open(self.gff,'r')
r = open('AA_sequence.txt', 'w')
w = open('CDS.txt', 'w')
hash_AA = {} # AA hash,sequence2AA
with open('AA.txt', 'r') as f: #AA.txt 为密码子表
for line in f:
line = line.strip()
if line:
line = self.pa.split(line)
hash_AA[line[0]] = line[1] #AA hash
hash_CDS = {} # CDS hash,CDS2sequence
for line in fn:
line = line.strip()
if line.startswith('>'):
A = self.pa.split(line)[0].replace('>', '')
hash_CDS[A] = ''
else:
hash_CDS[A] += line
fn.close()
for line in gf:
line = line.strip()
if 'CDS' in line:
i = self.Get_CDS_index(line)
sequence = hash_CDS[i[0]][int(i[1]) - 1:int(i[2])]
sequence = sequence[int(i[4]):] # i[4] 表示密码子开始位置
if i[3] == '-':
sequence = self.N2M(sequence)
#w.write(i[0] + '\n' + sequence + '\n')
#后面是一堆正则，主要是对序列做注释的，看文件格式而定 s1 = self.pa.split(line)
p1 = re.compile(r'ID=(.*?);.*?Dbxref=(.*?);.*?Name=(.*?);.*?gbkey=(.*?);.*?product=(.*?);.*?protein_id=(.*?);')
p2 = re.compile(r'.*?product=(.*?);.*?protein_id=(.*?);')
m1 = re.findall(p1,line)
m2 = re.findall(p2,line)
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]
w.write(s + '\n' + sequence + '\n')
p = '>' + s1[0]+'\tproduct:'+m2[0][0]+'\tproduct_id:'+m2[0][1]
AA = self.Seq2AA(sequence, hash_AA)
r.write(p + '\n' + AA + '\n')
w.close()
r.close() if __name__=='__main__':
fna = 'GCA_000160075.2_ASM16007v2_genomic.fna'
gff = 'GCA_000160075.2_ASM16007v2_genomic.gff'
m = CDS2AA(fna,gff)
m.CDS2AA()
出现的一些问题我会慢慢完善。后面的有意思作业题目会陆续上传。