我有一个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