前言
通过输入基因组fasta文件即可找出对应可读框并输出为新的fasta文件
代码
global number
number=1
def get_cds(seq):
scaffolddict={}
for i in range(len(seq)-2):
if seq[i]=='A' and seq[i:i+3]=='ATG':
k=i
orf=''
while k<len(seq)-2:
codon=seq[k:k+3]
orf+=codon
if codon=='TAG' or codon=='TGA' or codon=='TAA':
break
if len(orf)>60: #可读框小于20个氨基酸则不算入
scaffolddict[f'sequence{number}']=orf
return scaffolddict
def reverse_complementary(seq):
reverseseq=seq[::-1]
resultseq=''
basepair={"A":"T","T":"A","C":"G","G":"C"}
for s in reverseseq:
resultseq+=basepair[s]
return resultseq
genomedict={}
#请修改输入文件名
with open("input.fasta",'rt') as infile:
for line in infile:
if line[0]=='>':
key=line[1:].strip()
genomedict[key]=''
elif len(line.strip())>0:
genomedict[key]+=line.strip().upper()
seqdict={}
for seqname in genomedict:
seqdict.update(get_cds(genomedict[seqname]))
seqdict.update(get_cds(reverse_complementary(genomedict[seqname])))
with open("annotated_ORF.fasta",'wt') as outfile:
for seqname in seqdict:
outfile.write(f">{seqname}\n{seqdict[seqname]}\n")
包含功能:
- 读入和输入fasta文件
- 获取反向互补序列
- 识别可读框,默认最小长度为20个氨基酸
注意:该代码未经优化,对于较大基因组文件运行时间会较长,用户可根据个人需求进行多进程等优化