序列输入// 输出

1 分析或读取序列

         Bio.SeqIO.parse()将序列数据读成 SeqRecord 对象。它有两个参数:第一个参数是要读取的 handle。handle 是一个要读的打开文件,但是可以从命令行输出,或者可以从网上下载。第二个参数是一个小写的序列特征格式--我们不会为你猜想文件格式。到

http://biopython.org/wiki/SeqIO 查 看 支 持 的 文 件 格 式 。 这 会 返 回 一 个 给 出SeqRecord 对象的迭代器,它主要用于 for 循环中。有时你会发现自己在处理仅仅包含单一记录的文件。对于此种情况,Biopython 1.45 引入了 Bio.SeqIO.read()函数。这也使用 handle 和格式作为参数。如果有一个或仅仅一个记录,将作为一个 SeqRecord 对象返回。

1.1 读取序列文件

总体上来将 Bio.SeqIO.parse() 是用来读取序列文件并作为 SeqRecord 对 象 ,往往用在 for 循环中:

from Bio import SeqIO
handle = open("ls_orchid.fasta")
for seq_record in SeqIO.parse(handle, "fasta") :
print seq_record.id
print repr(seq_record.seq)
print len(seq_record.seq)

handle.close()

这个例子是 2.4 节的,它会载入包含 orchid DNA 序列的 FASTA 格式文件
19Biopython 中文指南
生物信息学论坛 http://www.bioxxx.cn 倾力奉献
ls_orchid.fasta. 如果想载入一个 GenBank 格式文件,例如 ls_orchid.gbk,改下文
件名和格式就可以了:
from Bio import SeqIO
handle = open("ls_orchid.gbk")
for seq_record in SeqIO.parse(handle, "genbank") :
print seq_record.id
print seq_record.seq
print len(seq_record.seq)
handle.close()
类似的,如果你想读取另一种文件格式的文件,同时假定 Bio.SeqIO.parse()支持该格式,你需要做的就是将格式字符串改为合适的,例如 "swiss"代表SwissProt 文件,或者"embl"代表 EMBL 文本文件。具体文件格式列表可参看 wiki页面 (http://biopython.org/wiki/SeqIO).

1.2 在序列文件中重复记录

在以上的例子中,我们往往使用 for 循环来遍历记录。你可以使用 for 循环和其他支持遍历的 python 对象(列表,元组及字符串)。Bio.SeqIO 返回的对象往往是 SeqRecord 对象的遍历。你需要轮流查看每一个记录,但是仅仅一次。附加一点是当你处理大文件时,使用 iterator 会节约你的内存。除了使用 for 循 环 ,你也可以使用.next()方法来逐步调试条目:
from Bio import SeqIO
handle = open("ls_orchid.fasta")
record_iterator = SeqIO.parse(handle, "fasta")
first_record = record_iterator.next()
print first_record.id
print first_record.description
second_record = record_iterator.next()
print second_record.id
print second_record.description
handle.close()

如果你使用.next(),没有更多的结果,你将得到一个特殊的 python 对象 None或者一个 StopIteration 错误。我们也考虑了一个特殊的情况,当你的序列文件中包含多个记录,但是你仅仅想要第一个。这种情况下,以下代码就很简明:

from Bio import SeqIO
first_record = SeqIO.parse(open("ls_orchid.gbk"), "genbank").next()

一个注意事项--像这样使用.next()方法将会忽略文件中任何添加的记录。当你的文件中仅仅有一条记录时,就像这一章中的某些在线的例子,或者包含单个染色体的 GenBank 文件,那么使用新的 Bio.SeqIO.read()函数替代。这将会检查是否有其他非期待的记录存在。

1.3 在一个序列文件

在前一节中,我们讨论了 Bio.SeqIO.parse()给出 SeqRecord 遍历,你逐个的得到了记录。最常用的是使用两一个顺序来存取记录。 Python 列表数据类型很适合这个,我们可以使用 list()把记录遍历转变为一个包含 SeqRecord 对象的列表:
from Bio import SeqIO
handle = open("ls_orchid.gbk")
records = list(SeqIO.parse(handle, "genbank"))
handle.close()
print "Found %i records" % len(records)
print "The last record"
last_record = records[-1] #using Python's list tricks
print last_record.id
print repr(last_record.seq)
print len(last_record.seq)
print "The first record"
first_record = records[0] #remember, Python counts from zero
print first_record.id
print repr(first_record.seq)
print len(first_record.seq)
Giving:
Found 94 records

The last record
Z78439.1
Seq('CATTGTTGAGATCACATAATAATTGATCGAGTTAATCTGGAGGATCTGT
TTACT...GCC', IUPACAmbiguousDNA())
592
The first record
Z78533.1
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGA
CCGTGG...CGC', IUPACAmbiguousDNA())
740

当然,你一可以继续对 SeqRecord 对象列表使用 for 循环。使用一个 list 比一个 iterator 更灵活(例如,你可以通过 list 长度来确定记录的数目),但是会使用更多的内存,因为它把记录一次性全部包括在内存里。

1.4 提取数据

假定你想要从 ls_orchid.gbk 文件中提取物种列表。让我们先来看下文件的第一条记录,看下物种信息存储在那个部分。
from Bio import SeqIO
record_iterator = SeqIO.parse(open("ls_orchid.gbk"), "genbank")
first_record = record_iterator.next()
print first_record
会给出以下的信息:
ID: Z78533.1
Name: Z78533
Desription: C.irapeanum 5.8