Python从任意读取fram计算orf

Question

我有一个大的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，输入为字符串：

^{pr2}$

其中Seq='''CTGATGCGTGCGCGCGGCCGCCTGCAGCCAGCGCGTCAGTTCCGGCGCCGCCGCGCGGCTGTAGTTCAGCGCGCCGCCGCGATCGACGGGCAGGTAATGGCCTTCGATGTCGATGCCGTCCGGCGGCGTGTTCGAGTTCGCGATCGAGCCGAACTTGCCGGTCTTGCGCGCCTCGACGTACGTGCCGTCGTCGACGTACTGGATCTTCAGGTCGACGCCGAGCCGCTGCCGCGCCTGCGCCTGCAGCGCCTGCAGCAGCACGTCGCGCTGGTCGCGCACGGTC'''请注意，这是上面大fasta文件格式的第3个条目。在

我的示例输出是：set([])，所以我显然做错了什么。我的代码甚至不能扩展到多个条目（也就是说，它只需要一个称为Seq）的DNA字符串

谁能给我指一下正确的方向吗？在

编辑：

注意：ATG是起始密码子（即ORF的开始），而TAG、TGA、和{}是终止密码子（即ORF的结束）。在

Answer 1

你的要求我不清楚。在你的问题中，你说“我想知道出现在任何序列的第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)

注意事项：

1. 通常你会从文件中读取fasta集合。因为这里是字符串格式，所以我们使用斯金吉奥使其易于与来自BioPython的SeqIO.parse兼容
2. get_orfs函数查找atg并返回每个atg的ORG。如果您还对第4帧到第6帧感兴趣，则需要序列的^{}。在
3. 如果您只对每个fasta序列中最长的ORF感兴趣，那么可以使用get_orfs函数return (max(orfs, key=len))
4. 如果你只对一个特定帧（例如帧3）中的ATG开始的orf感兴趣，则会稍微困难一些。ATG可能会从最简单的方法开始，而不是从最简单的框架中找到所有ORF。在

Answer 2

编辑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('')