Python在DNA序列中找到最长的ORF

Question

有人能告诉我一个简单的方法来计算DNA序列中最长的开放阅读框（ORF）吗？ATG是起始密码子（即ORF的开始），而{}、TGA、和{}是终止密码子（即ORF的结束）。在

下面是一些产生错误的代码（并使用名为BioPython的外部模块）：

import sys
from Bio import SeqIO

currentCid = ''
buffer = []

for record in SeqIO.parse(open(sys.argv[1]),"fasta"):
    cid = str(record.description).split('.')[0][1:]

    if currentCid == '':
        currentCid = cid
    else:
        if cid != currentCid:
            buffer.sort(key = lambda x : len(x[1]))
            print '>' + buffer[-1][0]
            print buffer[-1][1]
            currentCid = cid
            buffer = [(str(record.description),str(record.seq))]
        else:
            buffer.append((str(record.description),str(record.seq)))

buffer.sort(key = lambda x : len(x[1]))
print '>' + buffer[-1][0]
print buffer[-1][1]

有没有可能用最少的外部依赖来编写这个过程（或者至少让上面的代码正常工作）？在

我的输入如下：

^{pr2}$

我的输出应该是：

以ATG（即ORF的开始）开始并以TAG、TGA或{}作为终止密码子（即ORF的结尾）的最长子串。在

Answer 1

您应该研究正则表达式：

import re

max(re.findall(r'ATG(?:(?!TAA|TAG|TGA)...)*(?:TAA|TAG|TGA)',s), key = len)

有一个很好的教程here，重点介绍了如何使用带有DNA字符串的正则表达式

Answer 2

因为BioPython是一个成熟的、广泛使用的模块，专门为这类问题而设计，所以几乎没有理由回避它并重新发明轮子。也就是说，使用正则表达式来识别起始密码子很有用：

from Bio import Seq
import regex as re
startP = re.compile('ATG')
nuc = input_seq.replace('\n','')
longest = (0,)
for m in startP.finditer(nuc, overlapped=True):
    if len(Seq.Seq(nuc)[m.start():].translate(to_stop=True)) > longest[0]:
        pro = Seq.Seq(nuc)[m.start():].translate(to_stop=True)
        longest = (len(pro), 
                   m.start(), 
                   str(pro),
                   nuc[m.start():m.start()+len(pro)*3+3])

注意，这使用了regex模块，而不是re模块；前者允许更容易地识别重叠匹配。我们可以让BioPython计算三胞胎并寻找终止密码子，而不是尝试通过正则表达式来实现这一点。在

在这里，longest产生由ORF编码的蛋白质的长度、起始位点（注意，使用基于0的编号）、由ORF编码的蛋白质序列以及ORF本身的序列，包括终止密码子。在

^{pr2}$

Answer 3

看看这个：

https://www.kaggle.com/xiangma/orf-finder?scriptVersionId=6709465

如上面的链接所示，有两种方法可以做到这一点：

请注意，我设置了超过1000bp的ORF长度限制，您可以根据需要进行调整。在

第一个：

from Bio import SeqIO
records = SeqIO.parse('dna2.fasta', 'fasta')
for record in records:
    for strand, seq in (1, record.seq), (-1, record.seq.reverse_complement()):
        for frame in range(3):
            length = 3 * ((len(seq)-frame) // 3)
            for pro in seq[frame:frame+length].translate(table = 1).split("*")[:-1]:
                if 'M' in pro:
                    orf = pro[pro.find('M'):]
                    pos = seq[frame:frame+length].translate(table=1).find(orf)*3 + frame +1
                    if len(orf)*3 +3 > 1300:
                        print("{}...{} - length {}, strand {}, frame {}, pos {}, name {}".format\
                           (orf[:3], orf[-3:], len(orf)*3+3, strand, frame, pos, record.id))

第二个使用正则表达式：

^{pr2}$