使用 Biopython 的翻译功能后，如何跟踪核苷酸序列中起始密码子 (ATG) 的位置？答案

【问题标题】：How to track the position of a start codon (ATG) in a nucleotide sequence after using the translate function of Biopython?使用 Biopython 的翻译功能后，如何跟踪核苷酸序列中起始密码子 (ATG) 的位置？
【发布时间】：2018-09-12 06:18:24
【问题描述】：

我有一个包含一堆序列的 FASTA 文件，格式如下：

BMRat|XM_008846946.1 ATGAAGAACATCACAGAAGCCACCACCTTCATTCTCAAGGGACTCACAGACAATGTGGAACTACAGGTCA TCCTCTTTTTTCTCTTTCTAGCGATTTATCTCTTCACTCTCATAGGAAATTTAGGACTTATTATTTTAGT TATTGGGGATTCAAAACTCCACAACCCTATGTACTGTTTTCTGAGTGTATTGTCTTCTGTAGATGCCTGC TATTCCTCAGACATCACCCCGAATATGTTAGTAGGCTTCCTGTCAAAAAACAAAGGCATTTCTCTCCATG GATGTGCAACACAGTTGTTTCTCGCTGTTACTTTTGGAACCACAGAATGCTTTCTGTTGGCGGCAATGGC TTATGACCGCTATGTAGCCATCCATGACCCACTTCTCTATGCAGTGAGCATGTCACCAAGGATCTATGTG CCGCTCATCATTGCTTCCTATGCTGGTGGAATTCTGCATGCGATTATCCACACCGTGGCCACCTTCAGCC TGTCCTTCTGTGGATCTAATGAAATCAGTCATATATTCTGTGACATCCCTCCTCTGCTGGCTATTTCTTG TTCTGACACTTACATCAATGAGCTCCTGTTGTTCTTCTTTGTGAGCTCCATAGAAATAGTCACTATCCTC ATCATCCTGGTCTCTTATGGTTTCATCCTTATGGCCATTCTGAAGATGAATTCAGCTGAAGGGAGGAGAA AAGCTTCTCTGCATGTGGGTCTCACCTAACTGGAGTGTCCATTTTCTATGGGACAAGCCTTTTCATGTA TGTGAGACCAAGCTCCAACTATTCCTTGGCACATGACATGGTAGTGTCGACATTTTATACCATTGTGATT CCCATGCTGAACCCTGTCATCTACAGTCTGAGGAACAAAGATGTGAAAGAGGCAATGAGAAGATTTTTGA AGAAAAATTTTCAGAAACTTTAA

使用 biopython http://biopython.org/wiki/Seq 实现的代码允许我找到 FASTA 文件中每个序列的最长氨基酸序列，该序列以甲硫氨酸开头并以终止密码子结尾。

函数是find_largest_polypeptide_in_DNA。基本上，它使用 3 个不同的正向阅读框将 DNA 序列翻译成氨基酸序列，并在变量allPossibilities 中保存以 M（特定氨基酸）开头并以终止密码子结尾的片段。然后它比较可能性的长度并选择最长的可能性，返回该片段的蛋白质序列。

def find_largest_polypeptide_in_DNA(seq, translationTable=1):
    allPossibilities = []
    for frame in range(3):
        trans = str(seq[frame:].translate(translationTable))
        framePossibilitiesF = [i[i.find("M"):] for i in trans.split("*") if "M" in i]
        allPossibilities += framePossibilitiesF
    allPossibilitiesLengths = [len(i) for i in allPossibilities]

    if len(allPossibilitiesLengths) == 0:
        raise Exception("no candidate ORFs")

    proteinAsString = allPossibilities[allPossibilitiesLengths.index(max(allPossibilitiesLengths))]

    return Seq(proteinAsString, alphabet=ProteinAlphabet)

它工作得很好，但现在我想获得与函数返回的蛋白质序列相对应的 DNA 序列。我需要在函数中添加一些行以获得两个序列，但我真的不知道如何。我不知道是否可以跟踪 i.find("M") 的每个蛋氨酸的位置，然后使用该位置在核苷酸序列中跟踪它。

谢谢。

【问题讨论】：

您想修改函数，使其返回 DNA 序列，而不是返回以 Met 开头并以 STOP 结尾的最长片段的氨基酸序列？
您是否有理由要忽略相反的 3 个阅读框？你知道你的基因是从哪条链转录而来的吗？
您的问题表明所需的序列必须包含终止密码子。您提供的代码包括 FASTA 文件中每个序列末尾的段，该段以 M 开头且未终止。是否要在 FASTA 序列的末尾包含或排除以 Met 开头但未终止的片段？
是的，我也应该以相反的方式阅读它，谢谢。 @BennettBrown
我想要以 M 开头并以终止密码子 @BennettBrown 结尾的最长片段

标签： python sequence biopython dna-sequence

【解决方案1】：

我认为按照类似的原则编写一个新函数是最简单的。您的想法“跟踪 i.find('M') 的每个蛋氨酸的位置”基本上是在下面完成的。使用您开始的代码执行此操作的困难在于序列被split('*') 切碎，因此DNA起始位置是阅读框偏移量加上序列之前的所有片段密码子的总和忧虑。根据您的说明，我添加了一个封闭循环来迭代向前和向后的方向。

def find_largest_polypeptide_in_DNA(seq, translationTable=1):
    # Set the record to start with, then try to beat it
    longest_DNA = ''
    longest_amino_acid_sequence = 0

    for direction in [-1, 1]:
        forward_DNA = seq[::direction]
        # Check all three reading frames in this direction.
        for frame in range(3):
            trans = str(forward_DNA[frame:].translate(translationTable))
            cut_codons = 0
            while 'M' in trans:
                codons_before_Met = trans.find('M')
                cut_codons += codons_before_Met
                trans = trans[codons_before_Met:]
                if '*' in trans:
                    length = trans.find('*') + 1 
                    if length > longest_amino_acid_sequence:
                        longest_amino_acid_sequence = length
                        first_bp = frame + 3*cut_codons
                        last_bp = frame + 3*cut_codons + 3*(length)
                        longest_DNA = str(forward_DNA[first_bp:last_bp+1])
                    trans = trans[length:]
                else:
                    # Ignore sequence M... if ORF extends beyond FASTA?
                    trans = ''
    return longest_DNA

【讨论】：