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

BMRat|XM_008846946.1

ATGAAGAACATCACAGAAGCCACCACCTTCATTCTCAAGGGACTCACAGACAATGTGGAACTACAGGTCA

TCCTCTTTTTTCTCTTTCTAGCGATTTATCTCTTCACTCTCATAGGAAATTTAGGACTTATTATTTTAGT

TATTGGGGATTCAAAACTCCACAACCCTATGTACTGTTTTCTGAGTGTATTGTCTTCTGTAGATGCCTGC

TATTCCTCAGACATCACCCCGAATATGTTAGTAGGCTTCCTGTCAAAAAACAAAGGCATTTCTCTCCATG

GATGTGCAACACAGTTGTTTCTCGCTGTTACTTTTGGAACCACAGAATGCTTTCTGTTGGCGGCAATGGC

TTATGACCGCTATGTAGCCATCCATGACCCACTTCTCTATGCAGTGAGCATGTCACCAAGGATCTATGTG

CCGCTCATCATTGCTTCCTATGCTGGTGGAATTCTGCATGCGATTATCCACACCGTGGCCACCTTCAGCC

TGTCCTTCTGTGGATCTAATGAAATCAGTCATATATTCTGTGACATCCCTCCTCTGCTGGCTATTTCTTG

TTCTGACACTTACATCAATGAGCTCCTGTTGTTCTTCTTTGTGAGCTCCATAGAAATAGTCACTATCCTC

ATCATCCTGGTCTCTTATGGTTTCATCCTTATGGCCATTCTGAAGATGAATTCAGCTGAAGGGAGGAGAA

AAGTCTTCTCTGCATGTGGGTCTCACCTAACTGGAGTGTCCATTTTCTATGGGACAAGCCTTTTCATGTA

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结束？

您是否有理由要忽略其他三个阅读框？ 您知道您的基因从哪条链转录？

您的问题表明所需的序列必须包含终止密码子。 您提供的代码在FASTA文件中每个序列的末尾都包含该段，该段以M开头且未终止。 您是否要在FASTA序列的末尾包含或排除以Met开头但未终止的片段？

是的，我也应该以相反的方式阅读它，谢谢。 @贝内特·布朗

我想要最长的段，以M开头，以终止密码子@BennettBrown结尾

相关脚本github.com/chris-rands/CR_bioinformatics_utilities/blob/master/

我认为遵循类似的原则编写新函数将是最容易的。 您的想法"跟踪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