python 从基因预测输出中提取序列

Question

我用SoftBerry做了基因预测，结果输出是这样的：

预测的蛋白质：

>FGENESH:[mRNA]   1  12 exon (s)   4267  -   6782  1296 bp, chain +
ATGATACGCACTGCGCTTTCACGAGCAGCGGCCATCGTCGCCGCCCGCACCTCCGCCAAG
CTCCGCCCTTCCCTCCTCGCTCGATCTCCGCCGTCCAGACTCCTCCACGATGGGATTAAC
GCCAACCCAGTTGCTCTTCAGATGATCAACTACGCCGTCTCTCTCGCCAGGTCTCAGAAA
>FGENESH:   1  12 exon (s)   4267  -   6782   431 aa, chain +
MIRTALSRAAAIVAARTSAKLRPSLLARSPPSRLLHDGINANPVALQMINYAVSLARSQK
SDESYGQAQLVLEQCLSSQPSEGQDLATHNSRAMVLMAMSTLLSERGKLDEAIEKLQKVE

等等：输出内容很长，所以手动编辑起来不太容易。

我需要找出以'>FGENESH:[mRNA]'开头的序列。所以，我尝试了这个：

for line in infile:
    if line.startswith('>FGENESH:[mRNA]'):
        print(line)
        outfile.write(line)

结果只给我了标题行：

>FGENESH:[mRNA]   1  12 exon (s)   4267  -   6782  1296 bp, chain +

不过，我希望输出能像这样：

>FGENESH:[mRNA]   1  12 exon (s)   4267  -   6782  1296 bp, chain +
ATGATACGCACTGCGCTTTCACGAGCAGCGGCCATCGTCGCCGCCCGCACCTCCGCCAAG
CTCCGCCCTTCCCTCCTCGCTCGATCTCCGCCGTCCAGACTCCTCCACGATGGGATTAAC
GCCAACCCAGTTGCTCTTCAGATGATCAACTACGCCGTCTCTCTCGCCAGGTCTCAGAAA

有没有人能告诉我怎么才能得到这个结果？

我会非常感激的，因为我还是个新手。

谢谢。

jd

Answer 1

也许可以用正则表达式来解决这个问题：

import re

data = '''>FGENESH:[mRNA]   1  12 exon (s)   4267  -   6782  1296 bp, chain +
ATGATACGCACTGCGCTTTCACGAGCAGCGGCCATCGTCGCCGCCCGCACCTCCGCCAAG
CTCCGCCCTTCCCTCCTCGCTCGATCTCCGCCGTCCAGACTCCTCCACGATGGGATTAAC
GCCAACCCAGTTGCTCTTCAGATGATCAACTACGCCGTCTCTCTCGCCAGGTCTCAGAAA
>FGENESH:   1  12 exon (s)   4267  -   6782   431 aa, chain +
MIRTALSRAAAIVAARTSAKLRPSLLARSPPSRLLHDGINANPVALQMINYAVSLARSQK
SDESYGQAQLVLEQCLSSQPSEGQDLATHNSRAMVLMAMSTLLSERGKLDEAIEKLQKVE
>FGENESH:[mRNA]   1  12 exon (s)   4267  -   6782  1296 bp, chain +
ATGATACGCACTGCGCTTTCACGAGCAGCGGCCATCGTCGCCGCCCGCACCTCCGCCAAG
CTCCGCCCTTCCCTCCTCGCTCGATCTCCGCCGTCCAGACTCCTCCACGATGGGATTAAC
GCCAACCCAGTTGCTCTTCAGATGATCAACTACGCCGTCTCTCTCGCCAGGTCTCAGAAA'''

pattern = re.compile(r'^(>FGENESH:\[mRNA\].*?[ACTG\n]+)$', re.MULTILINE|re.DOTALL)
for match in pattern.findall(data):
    print match

输出结果：

>FGENESH:[mRNA]   1  12 exon (s)   4267  -   6782  1296 bp, chain +
ATGATACGCACTGCGCTTTCACGAGCAGCGGCCATCGTCGCCGCCCGCACCTCCGCCAAG
CTCCGCCCTTCCCTCCTCGCTCGATCTCCGCCGTCCAGACTCCTCCACGATGGGATTAAC
GCCAACCCAGTTGCTCTTCAGATGATCAACTACGCCGTCTCTCTCGCCAGGTCTCAGAAA
>FGENESH:[mRNA]   1  12 exon (s)   4267  -   6782  1296 bp, chain +
ATGATACGCACTGCGCTTTCACGAGCAGCGGCCATCGTCGCCGCCCGCACCTCCGCCAAG
CTCCGCCCTTCCCTCCTCGCTCGATCTCCGCCGTCCAGACTCCTCCACGATGGGATTAAC
GCCAACCCAGTTGCTCTTCAGATGATCAACTACGCCGTCTCTCTCGCCAGGTCTCAGAAA

Answer 2

你的基因输出是标准的FASTA格式，这种格式是任何生物信息学软件都能准确读取的；使用这样的模块可以让后续的工作更加灵活。例如，BioPython模块（http://biopython.org/wiki/Biopython）可以帮助你提取mRNA序列：

from Bio import SeqIO 
from StringIO import StringIO

output = []

input = """>FGENESH:[mRNA]   1  12 exon (s)   4267  -   6782  1296 bp, chain +
ATGATACGCACTGCGCTTTCACGAGCAGCGGCCATCGTCGCCGCCCGCACCTCCGCCAAG
CTCCGCCCTTCCCTCCTCGCTCGATCTCCGCCGTCCAGACTCCTCCACGATGGGATTAAC
GCCAACCCAGTTGCTCTTCAGATGATCAACTACGCCGTCTCTCTCGCCAGGTCTCAGAAA
>FGENESH:   1  12 exon (s)   4267  -   6782   431 aa, chain +
MIRTALSRAAAIVAARTSAKLRPSLLARSPPSRLLHDGINANPVALQMINYAVSLARSQK
SDESYGQAQLVLEQCLSSQPSEGQDLATHNSRAMVLMAMSTLLSERGKLDEAIEKLQKVE"""

for record in SeqIO.parse(StringIO(input), 'fasta'):
    if "FGENESH:[mRNA]" in record.description:
        output.append(record)

print output[0].format('fasta')

输出结果：

>FGENESH:[mRNA]  1 12 exon (s)  4267 -  6782 1296 bp, chain +
ATGATACGCACTGCGCTTTCACGAGCAGCGGCCATCGTCGCCGCCCGCACCTCCGCCAAG
CTCCGCCCTTCCCTCCTCGCTCGATCTCCGCCGTCCAGACTCCTCCACGATGGGATTAAC
GCCAACCCAGTTGCTCTTCAGATGATCAACTACGCCGTCTCTCTCGCCAGGTCTCAGAAA

之后你可以轻松地对这些序列做更多的操作，比如翻译它们：

print output[0].seq.translate() 

MIRTALSRAAAIVAARTSAKLRPSLLARSPPSRLLHDGINANPVALQMINYAVSLARSQK

或者高效地将它们写入文件：

with open("fgene.fasta", "w") as f:
    SeqIO.write(output, f, "fasta")

Answer 3

在编程中，有时候我们需要把一些代码放在特定的地方，这样它们才能正常工作。比如说，如果你在一个函数里定义了一个变量，这个变量就只能在这个函数内部使用，外面是不能访问的。这就像是在一个房间里放东西，只有在这个房间里你才能找到它。

另外，有些时候我们需要把一些代码放在全局范围内，这样在程序的任何地方都能使用。这就像是在一个大仓库里放东西，任何人都可以随时去拿。

总之，代码的作用范围就像是一个房间或仓库，决定了你能否在不同的地方使用这些代码。

flag = False

for line in infile:
    if flag is True:
        if line.startswith('>'):
            flag = False
        else:
            outfile.write(line)
    if line.startswith('>FGENESH:[mRNA]'):
        flag = True
        outfile.write(line)