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

Question

谁能告诉我如何计算 DNA 序列中最长的开放阅读框 (ORF) 的简单解决方案？ ATG 是起始密码子（即 ORF 的开头），TAG、TGA 和 TAA 是终止密码子（即 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]

是否可以用最少的外部依赖来编写这个过程（或者至少让上面的代码工作）？

我的输入如下所示：

ACCGCCGCGAACATCGCCGAGATCCTGCCGCCGCAGCCGAGCCGGCTGGTCGAGTATGCGCAACGACGCG
CGTCCGGCAGCATCCCGGCGATCATGGCGCGCTGGGATGCACGCGTACTGCAGGACAACGAACCATTCAC
CGCAGTCTATGGCGGCGCGTCGTACATCAACAACGACCTGTTCCTCGCCCGCCTCGCCGACTGGGGCGTG
TCGGCCGGCAACTACAGCGGCGAGATCGGCGGCGCGACACCGCCGCTGCGCTGGCGCCCGCTGCGGCTGC
TGCGTTCGCTGCCGGTGTTCTGGCGCATGCTGCGTGTCGCGCGCGGGCACCTGCCGACGCTCGAGCGCGG
CTTGCAGCGCTTCGACCAGGAACTCGCGACGCTCGTCGAGCGACGCGCCGACGGCCAGCAACTGGCCGAC
TGGTTCACGCGCTTCTACGTGTTCGTCGTGCAGGGCAACCTGTGCATCGCGTCGTCGCTGGCCAGCAGCG
GCGGCGCACTGTGGGGCCGTCCGCCGACCGCATACGGCCAGCTCGACGACAGCCCGCACCGGCTGCCGTG
GGAAACCGATCCGGGCACCGCACGGCCCGCGCCCACCCACCTGCCGCTGCAGGCGTTTCCCGCCTGGCCG
CTGCCGGTCCGCGTGCTCCACGCGCTCGGCGCGCCCGGCATGCGCGGCTGGTATCTGCAGGTGCGCGAGT
GGTATCGCGACAACCTGATGCGCGTGTTCTTCCGCCTGCATCATGCGATGCCGGCCGCCGATCGCGACAC
GTGGTTCGCGCCCCATCCCGATCGCCGCGAACGCAACGGCAGCTTCTGGCAGGACGGCGGCGAAGGCACC
GACGAGGCAGCCGGCTTCATGATCTATCCGGGCCACACGCAAGGCGTGCTCGGCCACGACATCCTGCTGG
AAGACACGCTCGACCCGGGCCGGCACGCGCAGTACCAGGCCGCGCGCGCCGTGATCGCGCGCATGGGCGG
CCGGCTGTCGCACGGCGCGACGCTGCTGCGCGAGCTGCGCAAGCCGTCGGCCGTGCTGCCGCGCGTCGAT
GCGGCGTGGATCGGGCGCGAGGTGCGGCTCAGCGACGGCCAGCTGACGCTGGTCGAATGAACGCGATGCG
GTTGCCGCGCACCCGAGCACGGGCCCGGGCCTGAACTGCCGATCAGCGTACCGGCGTGCGGACGACTCCG
TCGACCTTCAGCGTGCGCCGGTCGTGCGCGGCTTCGTATTCGACCGTCTGCGCAGGCGTGACGGCGCCGT
ATGAATGGCCGTTCACGTAGACGGTGCCGTCCCGCAGCTCGACCCGGTCGCCGTTGACCGTCGCTGTGGC
CCGTTCACCCTGCAGCACCGCGCCCGAACAACCTGCAGTCGAAAAACTGCGGACCGACGTGCCCGGCATC
GCGGCGATCCCGCCCTGGTCCGCCGCATGCGCCGCGCTGCACGGCGGCGCATCCATGCTGCCGGCAGCGT
GGACCGCGCCGGCGCTGATGCCGCATCCGGCAAGCAGCGCAATCGTCATCGGCTTCAGATGGTTCATGGT
GAGCTCCGTTGTCCGCCGCCGCGGATCGATGACCGGCCGACGCCCGTGCTCGCATGGCAGGCCGGCCGGC
CGGATGCATCCAGTATGCGTCCGGTTCGCGGCATTCCGCCATCGTCGCCGATACCGCTCATCGCCGCCCG
GTTCGCTCCCGCAGCGGCCTCTGGAAGCACCTCCCGCGGGGCAACCCGTCCCCATGAAAATCCACCTTGA
TCAAGTTGCGACTCGCAACTATTATTGATTGCGATCCGCAACCTTTCCGGACCCGCCATGGACCTCATCG
ACGCTCCCGCCAAGCCCCGCGAAGCCACGATCCTCGAGCTGCGCGACTTCTCCCGCAAACTGGTTCGCGA
GCTCGGCTTCATGCGCGCGACGCTGGCCGACAGCGACTGGGCGCCTT

我的输出应该是：

以ATG 开头（即ORF 的开头）并以TAG、TGA 或TAA 作为终止密码子（即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 本身的序列，包括终止密码子。

(338,
 93,
 'MARWDARVLQDNEPFTAVYGGASYINNDLFLARLADWGVSAGNYSGEIGGATPPLRWRPLRLLRSLPVFWRMLRVARGHLPTLERGLQRFDQELATLVERRADGQQLADWFTRFYVFVVQGNLCIASSLASSGGALWGRPPTAYGQLDDSPHRLPWETDPGTARPAPTHLPLQAFPAWPLPVRVLHALGAPGMRGWYLQVREWYRDNLMRVFFRLHHAMPAADRDTWFAPHPDRRERNGSFWQDGGEGTDEAAGFMIYPGHTQGVLGHDILLEDTLDPGRHAQYQAARAVIARMGGRLSHGATLLRELRKPSAVLPRVDAAWIGREVRLSDGQLTLVE',
 'ATGGCGCGCTGGGATGCACGCGTACTGCAGGACAACGAACCATTCACCGCAGTCTATGGCGGCGCGTCGTACATCAACAACGACCTGTTCCTCGCCCGCCTCGCCGACTGGGGCGTGTCGGCCGGCAACTACAGCGGCGAGATCGGCGGCGCGACACCGCCGCTGCGCTGGCGCCCGCTGCGGCTGCTGCGTTCGCTGCCGGTGTTCTGGCGCATGCTGCGTGTCGCGCGCGGGCACCTGCCGACGCTCGAGCGCGGCTTGCAGCGCTTCGACCAGGAACTCGCGACGCTCGTCGAGCGACGCGCCGACGGCCAGCAACTGGCCGACTGGTTCACGCGCTTCTACGTGTTCGTCGTGCAGGGCAACCTGTGCATCGCGTCGTCGCTGGCCAGCAGCGGCGGCGCACTGTGGGGCCGTCCGCCGACCGCATACGGCCAGCTCGACGACAGCCCGCACCGGCTGCCGTGGGAAACCGATCCGGGCACCGCACGGCCCGCGCCCACCCACCTGCCGCTGCAGGCGTTTCCCGCCTGGCCGCTGCCGGTCCGCGTGCTCCACGCGCTCGGCGCGCCCGGCATGCGCGGCTGGTATCTGCAGGTGCGCGAGTGGTATCGCGACAACCTGATGCGCGTGTTCTTCCGCCTGCATCATGCGATGCCGGCCGCCGATCGCGACACGTGGTTCGCGCCCCATCCCGATCGCCGCGAACGCAACGGCAGCTTCTGGCAGGACGGCGGCGAAGGCACCGACGAGGCAGCCGGCTTCATGATCTATCCGGGCCACACGCAAGGCGTGCTCGGCCACGACATCCTGCTGGAAGACACGCTCGACCCGGGCCGGCACGCGCAGTACCAGGCCGCGCGCGCCGTGATCGCGCGCATGGGCGGCCGGCTGTCGCACGGCGCGACGCTGCTGCGCGAGCTGCGCAAGCCGTCGGCCGTGCTGCCGCGCGTCGATGCGGCGTGGATCGGGCGCGAGGTGCGGCTCAGCGACGGCCAGCTGACGCTGGTCGAATGA')

Answer 3

看看这个：

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

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

请注意，我将 ORF 长度限制设置为 1000bp 以上，您可以根据需要进行调整。

第一个：

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))

第二个，使用正则表达式：

from Bio import SeqIO
import re
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):
            index = frame
            while index < len(record) - 6:
                match = re.match('(ATG(?:\S{3})*?T(?:AG|AA|GA))', str(seq[index:]))
                if match:
                    orf = match.group()
                    index += len(orf)
                    if len(orf) > 1300:
                        pos = str(record.seq).find(orf) + 1 
                        print("{}...{} - length {}, strand {}, frame {}, pos {}, name {}".format\
                           (orf[:6], orf[-3:], len(orf), strand, frame, pos, record.id))
                else: index += 3