Python 找到 DNA 序列中最长的 ORF
Python find longest ORF in DNA sequence
有人可以告诉我如何计算 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 的结尾)。
你应该看看正则表达式:
import re
max(re.findall(r'ATG(?:(?!TAA|TAG|TGA)...)*(?:TAA|TAG|TGA)',s), key = len)
有一个很好的教程here,它着重于将正则表达式与 DNA 字符串一起使用
由于 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')
看看这个:
https://www.kaggle.com/xiangma/orf-finder?scriptVersionId=6709465
如上link所示,有两种方法:
请注意,我将 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
有人可以告诉我如何计算 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 的结尾)。
你应该看看正则表达式:
import re
max(re.findall(r'ATG(?:(?!TAA|TAG|TGA)...)*(?:TAA|TAG|TGA)',s), key = len)
有一个很好的教程here,它着重于将正则表达式与 DNA 字符串一起使用
由于 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')
看看这个:
https://www.kaggle.com/xiangma/orf-finder?scriptVersionId=6709465
如上link所示,有两种方法:
请注意,我将 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