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。

下载并保存这个文件到你的Biopython样本目录，命名为 “orchid.fasta “。

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")

This method can be effectively used to convert the format as specified below −

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', ...)
   ]
}