Biopython - 序列 I/O 操作
Biopython 提供了一个模块 Bio.SeqIO 来分别从文件(任何流)读取和写入序列。 它支持生物信息学中几乎所有可用的文件格式。 大多数软件针对不同的文件格式提供不同的方法。 但是,Biopython 有意识地遵循一种单一的方法,通过其 SeqRecord 对象将解析后的序列数据呈现给用户。
让我们在下一节中了解有关 SeqRecord 的更多信息。
SeqRecord
Bio.SeqRecord 模块提供 SeqRecord 来保存序列的元信息以及序列数据本身,如下所示 −
seq − 这是一个实际的序列。
id − 它是给定序列的主要标识符。 默认类型为字符串。
name − 它是序列的名称。 默认类型为字符串。
description − 它显示有关序列的人类可读信息。
annotations − 它是有关序列的附加信息的字典。
SeqRecord 可以按照下面指定的方式导入
from Bio.SeqRecord import SeqRecord
让我们在接下来的部分中了解使用真实序列文件解析序列文件的细微差别。
解析序列文件格式
本节介绍如何解析两种最流行的序列文件格式,FASTA 和 GenBank。
FASTA
FASTA 是存储序列数据的最基本的文件格式。 FASTA最初是生物信息学早期发展过程中开发的用于DNA和蛋白质序列比对的软件包,主要用于搜索序列相似性。
Biopython 提供了一个示例 FASTA 文件,可以在 https://github.com/biopython/biopython/blob/master/Doc/examples/ls_orchid.fasta 访问。
下载此文件并将其作为 'orchid.fasta' 保存到您的 Biopython 示例目录中。
Bio.SeqIO 模块提供了parse()方法来处理序列文件,可以按如下方式导入 −
from Bio.SeqIO import parse
parse() 方法包含两个参数,第一个是文件句柄,第二个是文件格式。
>>> file = open('path/to/biopython/sample/orchid.fasta') >>> for record in parse(file, "fasta"): ... print(record.id) ... gi|2765658|emb|Z78533.1|CIZ78533 gi|2765657|emb|Z78532.1|CCZ78532 .......... .......... gi|2765565|emb|Z78440.1|PPZ78440 gi|2765564|emb|Z78439.1|PBZ78439 >>>
这里,parse() 方法返回一个可迭代对象,该对象在每次迭代时返回 SeqRecord。 作为可迭代的,它提供了许多复杂而简单的方法,让我们看到了一些特性。
next()
next() 方法返回可迭代对象中可用的下一个项目,我们可以使用它来获取第一个序列,如下所示 −
>>> first_seq_record = next(SeqIO.parse(open('path/to/biopython/sample/orchid.fasta'),'fasta')) >>> first_seq_record.id 'gi|2765658|emb|Z78533.1|CIZ78533' >>> first_seq_record.name 'gi|2765658|emb|Z78533.1|CIZ78533' >>> first_seq_record.seq Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC', SingleLetterAlphabet()) >>> first_seq_record.description 'gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA' >>> first_seq_record.annotations {} >>>
这里seq_record.annotations是空的,因为FASTA格式不支持序列注释。
列表理解
我们可以使用列表理解将可迭代对象转换为列表,如下所示
>>> seq_iter = SeqIO.parse(open('path/to/biopython/sample/orchid.fasta'),'fasta') >>> all_seq = [seq_record for seq_record in seq_iter] >>> len(all_seq) 94 >>>
在这里,我们使用了 len 方法来获取总计数。 我们可以得到最大长度的序列如下 −
>>> seq_iter = SeqIO.parse(open('path/to/biopython/sample/orchid.fasta'),'fasta') >>> max_seq = max(len(seq_record.seq) for seq_record in seq_iter) >>> max_seq 789 >>>
我们也可以使用下面的代码过滤序列 −
>>> seq_iter = SeqIO.parse(open('path/to/biopython/sample/orchid.fasta'),'fasta') >>> seq_under_600 = [seq_record for seq_record in seq_iter if len(seq_record.seq) < 600] >>> for seq in seq_under_600: ... print(seq.id) ... gi|2765606|emb|Z78481.1|PIZ78481 gi|2765605|emb|Z78480.1|PGZ78480 gi|2765601|emb|Z78476.1|PGZ78476 gi|2765595|emb|Z78470.1|PPZ78470 gi|2765594|emb|Z78469.1|PHZ78469 gi|2765564|emb|Z78439.1|PBZ78439 >>>
将一组 SqlRecord 对象(解析后的数据)写入文件就像调用 SeqIO.write 方法一样简单,如下所示 −
file = open("converted.fasta", "w) SeqIO.write(seq_record, file, "fasta")
这个方法可以有效的用来转换下面指定的格式 −
file = open("converted.gbk", "w) SeqIO.write(seq_record, file, "genbank")
GenBank
它是一种更丰富的基因序列格式,包括各种注释的字段。 Biopython 提供了一个示例 GenBank 文件,可以在 https://github.com/biopython/biopython/blob/master/Doc/examples/ls_orchid.fasta 上访问它。
将文件下载并保存到您的 Biopython 示例目录中,作为 'orchid.gbk'
因为,Biopython提供了一个单一的功能,parse来解析所有的生物信息学格式。 解析 GenBank 格式就像更改 parse 方法中的格式选项一样简单。
下面给出了相同的代码 −
>>> from Bio import SeqIO >>> from Bio.SeqIO import parse >>> seq_record = next(parse(open('path/to/biopython/sample/orchid.gbk'),'genbank')) >>> seq_record.id 'Z78533.1' >>> seq_record.name 'Z78533' >>> seq_record.seq Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC', IUPACAmbiguousDNA()) >>> seq_record.description 'C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA' >>> seq_record.annotations { 'molecule_type': 'DNA', 'topology': 'linear', 'data_file_division': 'PLN', 'date': '30-NOV-2006', 'accessions': ['Z78533'], 'sequence_version': 1, 'gi': '2765658', 'keywords': ['5.8S ribosomal RNA', '5.8S rRNA gene', 'internal transcribed spacer', 'ITS1', 'ITS2'], 'source': 'Cypripedium irapeanum', 'organism': 'Cypripedium irapeanum', 'taxonomy': [ 'Eukaryota', 'Viridiplantae', 'Streptophyta', 'Embryophyta', 'Tracheophyta', 'Spermatophyta', 'Magnoliophyta', 'Liliopsida', 'Asparagales', 'Orchidaceae', 'Cypripedioideae', 'Cypripedium'], 'references': [ Reference(title = 'Phylogenetics of the slipper orchids (Cypripedioideae: Orchidaceae): nuclear rDNA ITS sequences', ...), Reference(title = 'Direct Submission', ...) ] }