Python 从任意阅读框计算 ORF
Python calculate ORFs from any arbitrary reading frame
我有一个这种格式的大 fasta 文件：
>gi|142022655|gb|EQ086233.1|522 marine metagenome JCVI_SCAF_1096627390048 genomic scaffold, whole genome shotgun sequence
AAGACGGGCACCGTGTCCTTCGCGACGTACTCCGACCAGTTGTACACGTTCAGGTTGGTGTCGCCGGCAT
GGGCCGACAGGCTGGCCGCGACGGCCAGCGCCGCCGACGTGACGCGCGCGGCGCGCAACGCCGATTGACG
ACGGATACGGATACGCATGGGGATTCTCCTTGTGATGGGGATCGGCCGTTGCGCCCGGTCCGGGTCCGGA
CTCGCGTCAACGCCGTCGAGCGGTGTTCAGCACAAGGGCCAATGTAGAGATCGCGGCCGGCAGCGTCAGT
CCCGAAAACCGGGACAAACGGCGACGTCGATTCCCGCCGTTTGGGTAGATTCCCGCGTAGGCAGTCGAAA
ATATTCGTGATACCTGTAGCGCCACCTGAAAATCTTCGATACACGACGCCATGAGCGCTGCGCTGCCCGC
CCCCGATCTTCCGCTGAGCCACGTCGCGTTCGTGACTGAAACGCTGGGCGACATCGCACAAGCCGTCGGA
ACGCCGCAGTTCATGCGCGCCGTCTACGACACGCTCGTGCGCTACGTCGATTTCGACGCCGTGCACCTCG
ACTACGAGCGCAGCGCGTCTTCCGGCCGGCGCAGCGTCGGCTGGATCGGCAGCTTCGGCCGCGAGCCCGA
GCTGGTCGCGCAGGTGATGCGCCACTACTACCGCAGCTACGCGAGCGACGATGCAACTTACGCGGCGATC
GAAACCGAAAACGACGTGCAATTGCTGCAGGTGTCCGCGCAACGCGTGTCGAGCGAGCTACGGCATCTGT
TCTTCGATGCCGGCGACATTCATGACGAATGCGTGATCGCCGGCGTGACGGGCGGCACGCGCTACTCGAT
CTCGATCGCGCGCTCACGGCGGCTGCCGCCGTTTTCGCTGAAGGAACTGAGCCTGCTGAAGCAGCTTTCG
CAAGTCGTGCTGCCGCTGGCGTCCGCGCACAAGCGCCTGCTCGGCGCGATCTCCGCCGACGACGCACCGC
GCGACGAACTCGATCTCGACCTCGTCGCGCAATGGCTGCCGGAATGGCAGGAACGGTTGACCGCGCGCGA
GATGCATGTGTGTGCGTCGTTCATCCAGGGCATGACGTCGGCGGCCATCGCCCAATCGATGGGGCTCAAG
ACCTCCACCGTCGATACCTACGCGAAGCGCGCCTTCGCGAAGCTCGGCGTCGATTCGCGAAGGCAACTGA
TGACCCTCGTGCTGAGAAACGCGTCGCGGCGGCATGACGCATAGCATCC
>gi|142022655|gb|EQ086233.1|598 marine metagenome JCVI_SCAF_1096627390048 genomic scaffold, whole genome shotgun sequence
TTGCCGCCGGCCGCAGCCGGCTTGGCACCACGCTGCGGCTGGTCGCCGGACTTCGGCTTCGCGCCGGTGT
CCGCCGGCGCTGCCGGCCGCTTCGCGTTGCGCTCCTGCTTGGCCTTCGCTGCGAGCTGCGCCCGCAATTC
GGCAAGTTGTTCAAAACCCATAAATTCAATCCACCAGGAATATAAGGTGTGGTTCGTGCGGCCATGCCGC
GCGGCGCACGAGCTTCGCCGCCATGCGTGCGACCCGTCTGCCGCCGATGCGGAATACTACGGGGCCGCAT
>gi|142022655|gb|EQ086233.1|143 marine metagenome JCVI_SCAF_1096627390048 genomic scaffold, whole genome shotgun sequence
CTGATGCGTGCGCGCGGCCGCCTGCAGCCAGCGCGTCAGTTCCGGCGCCGCCGCGCGGCTGTAGTTCAGCGCG
CCGCCGCGATCGACGGGCAGGTAATGGCCTTCGATGTCGATGCCGTCCGGCGGCGTGTTCGAGTTCGCGA
TCGAGCCGAACTTGCCGGTCTTGCGCGCCTCGACGTACGTGCCGTCGTCGACGTACTGGATCTTCAGGTC
GACGCCGAGCCGCTGCCGCGCCTGCGCCTGCAGCGCCTGCAGCAGCACGTCGCGCTGGTCGCGCACGGTC
我想找出出现在任何序列的阅读框 3 中的最长开放阅读框 (ORF) 的长度？
到目前为止，我已经尝试了一些列出一个序列的所有 ORF 的代码，以字符串形式输入：
import re
from string import maketrans
pattern = re.compile(r'(?=(ATG(?:...)*?)(?=TAG|TGA|TAA))')
def revcomp(dna_seq):
return dna_seq[::-1].translate(maketrans("ATGC","TACG"))
def orfs(dna):
return set(pattern.findall(dna) + pattern.findall(revcomp(dna)))
print orfs(Seq)
其中 Seq='''CTGATGCGTGCGCGCGGCCGCCTGCAGCCAGCGCGTCAGTTCCGGCGCCGCCGCGCGGCTGTAGTTCAGCGCGCCGCCGCGATCGACGGGCAGGTAATGGCCTTCGATGTCGATGCCGTCCGGCGGCGTGTTCGAGTTCGCGATCGAGCCGAACTTGCCGGTCTTGCGCGCCTCGACGTACGTGCCGTCGTCGACGTACTGGATCTTCAGGTCGACGCCGAGCCGCTGCCGCGCCTGCGCCTGCAGCGCCTGCAGCAGCACGTCGCGCTGGTCGCGCACGGTC''' 请注意，这是上述大型 fasta 文件格式中的第 3 个条目。
我的示例输出是：set([])，所以我显然做错了什么。我的代码甚至没有扩展到多个条目（即，它只需要一个 DNA 字符串，称为 Seq）
谁能给我指出正确的方向？
编辑:
N.B.: ATG 是起始密码子（即一个 ORF 的开始），TAG、TGA 和 TAA 是终止密码子（即 ORF 的结尾）。
编辑 1：完全重写以更好地匹配问题描述。
我不知道这里的确切文件格式，所以我假设它的运行方式与您显示的三个序列相同——一个接一个。
如果我没理解错的话，你在第三个序列中没有看到匹配项的原因是那里实际上没有匹配项。不过前两个有匹配项，如果你运行这个你会看到它们。
'''
import re
import string
with open('dna.txt', 'rb') as f:
data = f.read()
data = [x.split('\n', 1) for x in data.split('>')]
data = [(x[0], ''.join(x[1].split())) for x in data if len(x) == 2]
start, end = [re.compile(x) for x in 'ATG TAG|TGA|TAA'.split()]
revtrans = string.maketrans("ATGC","TACG")
def get_longest(starts, ends):
''' Simple brute-force for now. Optimize later...
Given a list of start locations and a list
of end locations, return the longest valid
string. Returns tuple (length, start position)
Assume starts and ends are sorted correctly
from beginning to end of string.
'''
results = {}
# Use smallest end that is bigger than each start
ends.reverse()
for start in starts:
for end in ends:
if end > start and (end - start) % 3 == 0:
results[start] = end + 3
results = [(end - start, start) for
start, end in results.iteritems()]
return max(results) if results else (0, 0)
def get_orfs(dna):
''' Returns length, header, forward/reverse indication,
and longest match (corrected if reversed)
'''
header, seqf = dna
seqr = seqf[::-1].translate(revtrans)
def readgroup(seq, group):
return list(x.start() for x in group.finditer(seq))
f = get_longest(readgroup(seqf, start), readgroup(seqf, end))
r = get_longest(readgroup(seqr, start), readgroup(seqr, end))
(length, index), s, direction = max((f, seqf, 'forward'), (r, seqr, 'reverse'))
return length, header, direction, s[index:index + length]
# Process entire file
all_orfs = [get_orfs(x) for x in data]
# Put in groups of 3
all_orfs = zip(all_orfs[::3], all_orfs[1::3], all_orfs[2::3])
# Process each group of 3
for x in all_orfs:
x = max(x) # Only pring longest in each group
print(x)
print('')
我不清楚你的要求。在您的问题中，您说“我希望能够找出出现在任何序列的阅读框 3 中的最长开放阅读框 (ORF) 的长度？” 在随后的评论中，你说“我只对任意阅读框的每个 ORF 的综合计算感兴趣。它们对我来说都很重要。”
假设后者是您所追求的，这里有一种从 fasta 格式的序列集合中获取所有 ORF 的简单方法，使用 BioPython 来处理大部分工作。
import io # Only needed because input is in string form
from Bio import Seq, SeqIO
import regex as re
startP = re.compile('ATG')
def get_orfs(nuc):
orfs = []
for m in startP.finditer(nuc, overlapped=True):
pro = Seq.Seq(nuc)[m.start():].translate(to_stop=True)
orfs.append(nuc[m.start():m.start()+len(pro)*3+3])
return orfs
for fasta in SeqIO.parse(io.StringIO(fasta_inputs), 'fasta'):
header = fasta.description
orfs = get_orfs(str(fasta.seq))
print(header, orfs)
备注：
- 通常您会从文件中读取一个 fasta 集合。由于这里是字符串格式，我们使用 io.StringIO 使其与 BioPython
中的 SeqIO.parse 轻松兼容
get_orfs 函数查找 ATG 和 returns 来自每个的 ORG。如果您还对第 4 帧到第 6 帧感兴趣，则需要序列的 reverse_complement。- 如果您只对每个 fasta 序列中最长的 ORF 感兴趣，您可以使用 get_orfs 函数
return (max(orfs, key=len))
- 如果您只对特定框架（例如框架 3）中以 ATG 开头的 ORF 感兴趣，这会稍微困难一些。最简单的方法可能是简单地从框架中找到所有 ORF，然后丢弃那些不是以 ATG 开头的。
我有一个这种格式的大 fasta 文件：
>gi|142022655|gb|EQ086233.1|522 marine metagenome JCVI_SCAF_1096627390048 genomic scaffold, whole genome shotgun sequence
AAGACGGGCACCGTGTCCTTCGCGACGTACTCCGACCAGTTGTACACGTTCAGGTTGGTGTCGCCGGCAT
GGGCCGACAGGCTGGCCGCGACGGCCAGCGCCGCCGACGTGACGCGCGCGGCGCGCAACGCCGATTGACG
ACGGATACGGATACGCATGGGGATTCTCCTTGTGATGGGGATCGGCCGTTGCGCCCGGTCCGGGTCCGGA
CTCGCGTCAACGCCGTCGAGCGGTGTTCAGCACAAGGGCCAATGTAGAGATCGCGGCCGGCAGCGTCAGT
CCCGAAAACCGGGACAAACGGCGACGTCGATTCCCGCCGTTTGGGTAGATTCCCGCGTAGGCAGTCGAAA
ATATTCGTGATACCTGTAGCGCCACCTGAAAATCTTCGATACACGACGCCATGAGCGCTGCGCTGCCCGC
CCCCGATCTTCCGCTGAGCCACGTCGCGTTCGTGACTGAAACGCTGGGCGACATCGCACAAGCCGTCGGA
ACGCCGCAGTTCATGCGCGCCGTCTACGACACGCTCGTGCGCTACGTCGATTTCGACGCCGTGCACCTCG
ACTACGAGCGCAGCGCGTCTTCCGGCCGGCGCAGCGTCGGCTGGATCGGCAGCTTCGGCCGCGAGCCCGA
GCTGGTCGCGCAGGTGATGCGCCACTACTACCGCAGCTACGCGAGCGACGATGCAACTTACGCGGCGATC
GAAACCGAAAACGACGTGCAATTGCTGCAGGTGTCCGCGCAACGCGTGTCGAGCGAGCTACGGCATCTGT
TCTTCGATGCCGGCGACATTCATGACGAATGCGTGATCGCCGGCGTGACGGGCGGCACGCGCTACTCGAT
CTCGATCGCGCGCTCACGGCGGCTGCCGCCGTTTTCGCTGAAGGAACTGAGCCTGCTGAAGCAGCTTTCG
CAAGTCGTGCTGCCGCTGGCGTCCGCGCACAAGCGCCTGCTCGGCGCGATCTCCGCCGACGACGCACCGC
GCGACGAACTCGATCTCGACCTCGTCGCGCAATGGCTGCCGGAATGGCAGGAACGGTTGACCGCGCGCGA
GATGCATGTGTGTGCGTCGTTCATCCAGGGCATGACGTCGGCGGCCATCGCCCAATCGATGGGGCTCAAG
ACCTCCACCGTCGATACCTACGCGAAGCGCGCCTTCGCGAAGCTCGGCGTCGATTCGCGAAGGCAACTGA
TGACCCTCGTGCTGAGAAACGCGTCGCGGCGGCATGACGCATAGCATCC
>gi|142022655|gb|EQ086233.1|598 marine metagenome JCVI_SCAF_1096627390048 genomic scaffold, whole genome shotgun sequence
TTGCCGCCGGCCGCAGCCGGCTTGGCACCACGCTGCGGCTGGTCGCCGGACTTCGGCTTCGCGCCGGTGT
CCGCCGGCGCTGCCGGCCGCTTCGCGTTGCGCTCCTGCTTGGCCTTCGCTGCGAGCTGCGCCCGCAATTC
GGCAAGTTGTTCAAAACCCATAAATTCAATCCACCAGGAATATAAGGTGTGGTTCGTGCGGCCATGCCGC
GCGGCGCACGAGCTTCGCCGCCATGCGTGCGACCCGTCTGCCGCCGATGCGGAATACTACGGGGCCGCAT
>gi|142022655|gb|EQ086233.1|143 marine metagenome JCVI_SCAF_1096627390048 genomic scaffold, whole genome shotgun sequence
CTGATGCGTGCGCGCGGCCGCCTGCAGCCAGCGCGTCAGTTCCGGCGCCGCCGCGCGGCTGTAGTTCAGCGCG
CCGCCGCGATCGACGGGCAGGTAATGGCCTTCGATGTCGATGCCGTCCGGCGGCGTGTTCGAGTTCGCGA
TCGAGCCGAACTTGCCGGTCTTGCGCGCCTCGACGTACGTGCCGTCGTCGACGTACTGGATCTTCAGGTC
GACGCCGAGCCGCTGCCGCGCCTGCGCCTGCAGCGCCTGCAGCAGCACGTCGCGCTGGTCGCGCACGGTC
我想找出出现在任何序列的阅读框 3 中的最长开放阅读框 (ORF) 的长度？
到目前为止，我已经尝试了一些列出一个序列的所有 ORF 的代码，以字符串形式输入：
import re
from string import maketrans
pattern = re.compile(r'(?=(ATG(?:...)*?)(?=TAG|TGA|TAA))')
def revcomp(dna_seq):
return dna_seq[::-1].translate(maketrans("ATGC","TACG"))
def orfs(dna):
return set(pattern.findall(dna) + pattern.findall(revcomp(dna)))
print orfs(Seq)
其中 Seq='''CTGATGCGTGCGCGCGGCCGCCTGCAGCCAGCGCGTCAGTTCCGGCGCCGCCGCGCGGCTGTAGTTCAGCGCGCCGCCGCGATCGACGGGCAGGTAATGGCCTTCGATGTCGATGCCGTCCGGCGGCGTGTTCGAGTTCGCGATCGAGCCGAACTTGCCGGTCTTGCGCGCCTCGACGTACGTGCCGTCGTCGACGTACTGGATCTTCAGGTCGACGCCGAGCCGCTGCCGCGCCTGCGCCTGCAGCGCCTGCAGCAGCACGTCGCGCTGGTCGCGCACGGTC''' 请注意，这是上述大型 fasta 文件格式中的第 3 个条目。
我的示例输出是：set([])，所以我显然做错了什么。我的代码甚至没有扩展到多个条目（即，它只需要一个 DNA 字符串，称为 Seq）
谁能给我指出正确的方向？
编辑:
N.B.: ATG 是起始密码子（即一个 ORF 的开始），TAG、TGA 和 TAA 是终止密码子（即 ORF 的结尾）。
编辑 1：完全重写以更好地匹配问题描述。
我不知道这里的确切文件格式，所以我假设它的运行方式与您显示的三个序列相同——一个接一个。
如果我没理解错的话，你在第三个序列中没有看到匹配项的原因是那里实际上没有匹配项。不过前两个有匹配项，如果你运行这个你会看到它们。
'''
import re
import string
with open('dna.txt', 'rb') as f:
data = f.read()
data = [x.split('\n', 1) for x in data.split('>')]
data = [(x[0], ''.join(x[1].split())) for x in data if len(x) == 2]
start, end = [re.compile(x) for x in 'ATG TAG|TGA|TAA'.split()]
revtrans = string.maketrans("ATGC","TACG")
def get_longest(starts, ends):
''' Simple brute-force for now. Optimize later...
Given a list of start locations and a list
of end locations, return the longest valid
string. Returns tuple (length, start position)
Assume starts and ends are sorted correctly
from beginning to end of string.
'''
results = {}
# Use smallest end that is bigger than each start
ends.reverse()
for start in starts:
for end in ends:
if end > start and (end - start) % 3 == 0:
results[start] = end + 3
results = [(end - start, start) for
start, end in results.iteritems()]
return max(results) if results else (0, 0)
def get_orfs(dna):
''' Returns length, header, forward/reverse indication,
and longest match (corrected if reversed)
'''
header, seqf = dna
seqr = seqf[::-1].translate(revtrans)
def readgroup(seq, group):
return list(x.start() for x in group.finditer(seq))
f = get_longest(readgroup(seqf, start), readgroup(seqf, end))
r = get_longest(readgroup(seqr, start), readgroup(seqr, end))
(length, index), s, direction = max((f, seqf, 'forward'), (r, seqr, 'reverse'))
return length, header, direction, s[index:index + length]
# Process entire file
all_orfs = [get_orfs(x) for x in data]
# Put in groups of 3
all_orfs = zip(all_orfs[::3], all_orfs[1::3], all_orfs[2::3])
# Process each group of 3
for x in all_orfs:
x = max(x) # Only pring longest in each group
print(x)
print('')
我不清楚你的要求。在您的问题中，您说“我希望能够找出出现在任何序列的阅读框 3 中的最长开放阅读框 (ORF) 的长度？” 在随后的评论中，你说“我只对任意阅读框的每个 ORF 的综合计算感兴趣。它们对我来说都很重要。”
假设后者是您所追求的，这里有一种从 fasta 格式的序列集合中获取所有 ORF 的简单方法，使用 BioPython 来处理大部分工作。
import io # Only needed because input is in string form
from Bio import Seq, SeqIO
import regex as re
startP = re.compile('ATG')
def get_orfs(nuc):
orfs = []
for m in startP.finditer(nuc, overlapped=True):
pro = Seq.Seq(nuc)[m.start():].translate(to_stop=True)
orfs.append(nuc[m.start():m.start()+len(pro)*3+3])
return orfs
for fasta in SeqIO.parse(io.StringIO(fasta_inputs), 'fasta'):
header = fasta.description
orfs = get_orfs(str(fasta.seq))
print(header, orfs)
备注：
- 通常您会从文件中读取一个 fasta 集合。由于这里是字符串格式，我们使用 io.StringIO 使其与 BioPython 中的
get_orfs函数查找 ATG 和 returns 来自每个的 ORG。如果您还对第 4 帧到第 6 帧感兴趣，则需要序列的reverse_complement。- 如果您只对每个 fasta 序列中最长的 ORF 感兴趣，您可以使用 get_orfs 函数
return (max(orfs, key=len)) - 如果您只对特定框架（例如框架 3）中以 ATG 开头的 ORF 感兴趣，这会稍微困难一些。最简单的方法可能是简单地从框架中找到所有 ORF，然后丢弃那些不是以 ATG 开头的。
SeqIO.parse 轻松兼容