从基因组文件中找出可读框

前言

通过输入基因组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个氨基酸

注意：该代码未经优化，对于较大基因组文件运行时间会较长，用户可根据个人需求进行多进程等优化