如何在Python中查找开放阅读框架

Question

我正在使用Python和正则表达式来查找一个ORF（开放阅读框架）。

查找一个子字符串，它仅由字母ATGC（没有空格或换行符）组成，该子字符串应满足以下条件：

以ATG开头，以TAG、TAA或TGA结尾，并且应从第一个字符、第二个字符和第三个字符考虑序列：

Seq= "CCTCAGCGAGGACAGCAAGGGACTAGCCAGGAGGGAGAACAGAAACTCCAGAACATCTTGGAAATAGCTCCCAGAAAAGC
AAGCAGCCAACCAGGCAGGTTCTGTCCCTTTCACTCACTGGCCCAAGGCGCCACATCTCCCTCCAGAAAAGACACCATGA
GCACAGAAAGCATGATCCGCGACGTGGAACTGGCAGAAGAGGCACTCCCCCAAAAGATGGGGGGCTTCCAGAACTCCAGG
CGGTGCCTATGTCTCAGCCTCTTCTCATTCCTGCTTGTGGCAGGGGCCACCACGCTCTTCTGTCTACTGAACTTCGGGGT
GATCGGTCCCCAAAGGGATGAGAAGTTCCCAAATGGCCTCCCTCTCATCAGTTCTATGGCCCAGACCCTCACACTCAGAT
CATCTTCTCAAAATTCGAGTGACAAGCCTGTAGCCCACGTCGTAGCAAACCACCAAGTGGAGGAGCAGCTGGAGTGGCTG
AGCCAGCGCGCCAACGCCCTCCTGGCCAACGGCATGGATCTCAAAGACAACCAACTAGTGGTGCCAGCCGATGGGTTGTA
CCTTGTCTACTCCCAGGTTCTCTTCAAGGGACAAGGCTGCCCCGACTACGTGCTCCTCACCCACACCGTCAGCCGATTTG
CTATCTCATACCAGGAGAAAGTCAACCTCCTCTCTGCCGTCAAGAGCCCCTGCCCCAAGGACACCCCTGAGGGGGCTGAG
CTCAAACCCTGGTATGAGCCCATATACCTGGGAGGAGTCTTCCAGCTGGAGAAGGGGGACCAACTCAGCGCTGAGGTCAA
TCTGCCCAAGTACTTAGACTTTGCGGAGTCCGGGCAGGTCTACTTTGGAGTCATTGCTCTGTGAAGGGAATGGGTGTTCA
TCCATTCTCTACCCAGCCCCCACTCTGACCCCTTTACTCTGACCCCTTTATTGTCTACTCCTCAGAGCCCCCAGTCTGTA
TCCTTCTAACTTAGAAAGGGGATTATGGCTCAGGGTCCAACTCTGTGCTCAGAGCTTTCAACAACTACTCAGAAACACAA
GATGCTGGGACAGTGACCTGGACTGTGGGCCTCTCATGCACCACCATCAAGGACTCAAATGGGCTTTCCGAATTCACTGG
AGCCTCGAATGTCCATTCCTGAGTTCTGCAAAGGGAGAGTGGTCAGGTTGCCTCTGTCTCAGAATGAGGCTGGATAAGAT
CTCAGGCCTTCCTACCTTCAGACCTTTCCAGATTCTTCCCTGAGGTGCAATGCACAGCCTTCCTCACAGAGCCAGCCCCC
CTCTATTTATATTTGCACTTATTATTTATTATTTATTTATTATTTATTTATTTGCTTATGAATGTATTTATTTGGAAGGC
CGGGGTGTCCTGGAGGACCCAGTGTGGGAAGCTGTCTTCAGACAGACATGTTTTCTGTGAAAACGGAGCTGAGCTGTCCC
CACCTGGCCTCTCTACCTTGTTGCCTCCTCTTTTGCTTATGTTTAAAACAAAATATTTATCTAACCCAATTGTCTTAATA
ACGCTGATTTGGTGACCAGGCTGTCGCTACATCACTGAACCTCTGCTCCCCACGGGAGCCGTGACTGTAATCGCCCTACG
GGTCATTGAGAGAAATAA"

我尝试过的方法：

# finding the  stop codon here 

def stop_codon(seq_0):

        for i in range(0,len(seq_0),3):
            if (seq_0[i:i+3]== "TAA" and i%3==0) or (seq_0[i:i+3]== "TAG" and i%3==0) or (seq_0[i:i+3]== "TGA" and i%3==0) :
                a =i+3

                break

            else:
                a = None

# finding the start codon here 

startcodon_find =[m.start() for m in re.finditer('ATG', seq_0)]

我该如何找到检查起始密码子并找到第一个终止密码子的方法。然后找到下一个起始密码子和下一个终止密码子。

我希望将此运行三个框架。如前所述，三个框架将考虑序列的第一个、第二个和第三个字符作为起始点。

此外，序列需要被分成小的3部分。因此应该像这样：

ATG TTT AAA ACA AAA TAT TTA TCT AAC CCA ATT GTC TTA ATA ACG CTG ATT TGA

任何帮助都将不胜感激。 我最终的答案：

def orf_find(st0):

    seq_0=""
    for i in range(0,len(st0),3):
        if len(st0[i:i+3])==3:
            seq_0 = seq_0 + st0[i:i+3]+ " "

    ms_1 =[m.start() for m in re.finditer('ATG', seq_0)]
    ms_2 =[m.start() for m in re.finditer('(TAA)|(TAG)|(TGA)', seq_0)]

    def get_next(arr,value):
        for a in arr:
            if a > value:
                return a
        return -1




    codons = []
    start_codon=ms_1[0]
    while (True):
        stop_codon = get_next(ms_2,start_codon)
        if stop_codon == -1:
            break
        codons.append((start_codon,stop_codon))
        start_codon = get_next(ms_1,stop_codon)
        if start_codon==-1:
            break

    max_val = 0
    selected_tupple = ()
    for i in codons:
        k=i[1]-i[0]
        if k > max_val:
            max_val = k
            selected_tupple = i

    print "selected tupple is ", selected_tupple

    final_seq=seq_0[selected_tupple[0]:selected_tupple[1]+3]

    print final_seq
    print "The longest orf length is " + str(max_val)



output_file = open('Longorf.txt','w')
output_file.write(str(orf_find(st0)))

output_file.close()

上述的写入函数对我在将内容写入文本文件时没有帮助。我在那里得到的只是 NONE。为什么会出现这种错误？有人能帮忙吗？

- Nodnin

---

你必须从函数中返回输出。我刚刚编辑了我的答案。 - MoRe

Answer 1

如果您想手动编写它：

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)

- Stefan Gruenwald

---

你好，你有Python3版本吗？ - pippo1980

Answer 2

根据您的标签，我猜您应该知道Biopython。您是否已经查看了文档？http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc231 可能会有所帮助。

我稍微调整了上面链接中的代码以适用于您的序列：

from Bio.Seq import Seq

seq = Seq("CCTCAGCGAGGACAGCAAGGGACTAGCCAGGAGGGAGAACAGAAACTCCAGAACATCTTGGAAATAGCTCCCAGAAAAGCAAGCAGCCAACCAGGCAGGTTCTGTCCCTTTCACTCACTGGCCCAAGGCGCCACATCTCCCTCCAGAAAAGACACCATGAGCACAGAAAGCATGATCCGCGACGTGGAACTGGCAGAAGAGGCACTCCCCCAAAAGATGGGGGGCTTCCAGAACTCCAGGCGGTGCCTATGTCTCAGCCTCTTCTCATTCCTGCTTGTGGCAGGGGCCACCACGCTCTTCTGTCTACTGAACTTCGGGGTGATCGGTCCCCAAAGGGATGAGAAGTTCCCAAATGGCCTCCCTCTCATCAGTTCTATGGCCCAGACCCTCACACTCAGATCATCTTCTCAAAATTCGAGTGACAAGCCTGTAGCCCACGTCGTAGCAAACCACCAAGTGGAGGAGCAGCTGGAGTGGCTGAGCCAGCGCGCCAACGCCCTCCTGGCCAACGGCATGGATCTCAAAGACAACCAACTAGTGGTGCCAGCCGATGGGTTGTACCTTGTCTACTCCCAGGTTCTCTTCAAGGGACAAGGCTGCCCCGACTACGTGCTCCTCACCCACACCGTCAGCCGATTTGCTATCTCATACCAGGAGAAAGTCAACCTCCTCTCTGCCGTCAAGAGCCCCTGCCCCAAGGACACCCCTGAGGGGGCTGAGCTCAAACCCTGGTATGAGCCCATATACCTGGGAGGAGTCTTCCAGCTGGAGAAGGGGGACCAACTCAGCGCTGAGGTCAATCTGCCCAAGTACTTAGACTTTGCGGAGTCCGGGCAGGTCTACTTTGGAGTCATTGCTCTGTGAAGGGAATGGGTGTTCATCCATTCTCTACCCAGCCCCCACTCTGACCCCTTTACTCTGACCCCTTTATTGTCTACTCCTCAGAGCCCCCAGTCTGTATCCTTCTAACTTAGAAAGGGGATTATGGCTCAGGGTCCAACTCTGTGCTCAGAGCTTTCAACAACTACTCAGAAACACAAGATGCTGGGACAGTGACCTGGACTGTGGGCCTCTCATGCACCACCATCAAGGACTCAAATGGGCTTTCCGAATTCACTGGAGCCTCGAATGTCCATTCCTGAGTTCTGCAAAGGGAGAGTGGTCAGGTTGCCTCTGTCTCAGAATGAGGCTGGATAAGATCTCAGGCCTTCCTACCTTCAGACCTTTCCAGATTCTTCCCTGAGGTGCAATGCACAGCCTTCCTCACAGAGCCAGCCCCCCTCTATTTATATTTGCACTTATTATTTATTATTTATTTATTATTTATTTATTTGCTTATGAATGTATTTATTTGGAAGGCCGGGGTGTCCTGGAGGACCCAGTGTGGGAAGCTGTCTTCAGACAGACATGTTTTCTGTGAAAACGGAGCTGAGCTGTCCCCACCTGGCCTCTCTACCTTGTTGCCTCCTCTTTTGCTTATGTTTAAAACAAAATATTTATCTAACCCAATTGTCTTAATAACGCTGATTTGGTGACCAGGCTGTCGCTACATCACTGAACCTCTGCTCCCCACGGGAGCCGTGACTGTAATCGCCCTACGGGTCATTGAGAGAAATAA")


table = 1
min_pro_len = 100

for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]:
    for frame in range(3):
        for pro in nuc[frame:].translate(table).split("*"):
            if len(pro) >= min_pro_len:
                print "%s...%s - length %i, strand %i, frame %i" % (pro[:30], pro[-3:], len(pro), strand, frame)

ORF也被翻译了。您可以选择不同的翻译表。请查看http://biopython.org/DIST/docs/tutorial/Tutorial.html#sec:translation

编辑：代码解释：

代码的开始我会用你提供的字符串创建一个序列对象。注意 seq = Seq("ACGT")。 两个 for 循环创建了 6 种不同的框架。内部的 for 循环根据所选的翻译表翻译每个框架，并返回一个氨基酸链，其中每个终止密码子都编码为 *。 split 函数将删除这些占位符并拆分该字符串，从而得到可能的蛋白质序列列表。通过设置 min_pro_len，您可以定义检测到蛋白质的最小氨基酸链长度。1 是标准表。请参阅 http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi#SG1，在这里您可以看到起始密码子是 AUG（等于 ATG），终止密码子（核苷酸序列中的 *）是 TAA、TAG 和 TGA，就像您想要的那样。您也可以使用不同的翻译表。

当您添加

print nuc[frame:].translate(table)

在第二个for循环的内部，你会得到类似这样的东西：

PQRGQQGTSQEGEQKLQNILEIAPRKASSQPGRFCPFHSLAQGATSPSRKDTMSTESMIRDVELAEEALPQKMGGFQNSRRCLCLSLFSFLLVAGATTLFCLLNFGVIGPQRDEKFPNGLPLISSMAQTLTLRSSSQNSSDKPVAHVVANHQVEEQLEWLSQRANALLANGMDLKDNQLVVPADGLYLVYSQVLFKGQGCPDYVLLTHTVSRFAISYQEKVNLLSAVKSPCPKDTPEGAELKPWYEPIYLGGVFQLEKGDQLSAEVNLPKYLDFAESGQVYFGVIAL*REWVFIHSLPSPHSDPFTLTPLLSTPQSPQSVSF*LRKGIMAQGPTLCSELSTTTQKHKMLGQ*PGLWASHAPPSRTQMGFPNSLEPRMSIPEFCKGRVVRLPLSQNEAG*DLRPSYLQTFPDSSLRCNAQPSSQSQPPSIYICTYYLLFIYYLFICL*MYLFGRPGCPGGPSVGSCLQTDMFSVKTELSCPHLASLPCCLLFCLCLKQNIYLTQLS**R*FGDQAVATSLNLCSPREP*L*SPYGSLREI

（注意星号位于终止密码子位置）

编辑：回答您的第二个问题：

您必须返回一个字符串，您想要写入文件中。创建一个输出字符串，并在函数末尾返回它：

output = "selected tupple is " + str(selected_tupple) + "\n"
output += final_seq + "\n"
output += "The longest orf length is " + str(max_val) + "\n"
return output

- MoRe

---

给我一个错误：“STR”没有属性“reverse_complement”。 - Nodnin

你有注意到我发布的前3行代码吗？那里是将你的字符串转换为Sequence对象。 - MoRe

但我看不到，序列如何以“ATG”开头并以终止密码子结尾..我如何在一个阅读框架中找到多个子字符串。 - Nodnin

所以这里星号需要是终止密码子。我们需要编写不同的代码来查找起始密码子吗？这里提供的文档很好，但解释非常少。我尝试进行更改，但出现了错误，最重要的是在Python中是否需要导入nio模块... - Nodnin

我已经实现了这段代码，但我想知道如何显示整个蛋白质序列。目前在最终的打印语句中，代码似乎是有意截断它。 - olliepower

Answer 3

以类似的方式工作：

#!/usr/bin/python3

import re

pattern = re.compile(r'(?=(ATG(?:...)*?)(?=TAG|TGA|TAA))')

Seq= """CATGCTCAGCGAGGACAGCAAGGGCCCATTTACAGGAGCATAGTAA"""

revcompseq = Seq[::-1].maketrans("ATGC", "TACG") #reverse complement

print (pattern.findall(Seq)) #forward search
print (pattern.findall(Seq[::-1].translate(revcompseq))) #backward search

使用这个小演示序列，你将得到以下结果:

['ATGCTCAGCGAGGACAGCAAGGGCCCATTTACAGGAGCA']
['ATGCTCCTG', 'ATGGGCCCTTGCTGTCCTCGC']

- Stefan Gruenwald