Biopython 简明教程

Biopython - Sequence I/O Operations

Biopython 提供了一个模块 Bio.SeqIO, 分别用于从文件（任何流）读取和写入序列。它支持生物信息学中几乎所有可用的文件格式。大多数软件为不同的文件格式提供了不同的方法。但是, Biopython 认真地遵循一种方法, 通过其 SeqRecord 对象向用户显示解析的序列数据。

让我们在以下部分详细了解 SeqRecord。

SeqRecord

Bio.SeqRecord 模块提供 SeqRecord 来保存序列的元信息以及序列数据本身, 如下所示 −

seq − 这是实际序列。
id − 这是给定序列的主标识符。默认类型为字符串。
name − 这是序列的名称。默认类型为字符串。
description − 它显示关于该序列的人类可读信息。
annotations − 这是一个关于该序列的附加信息的字典。

以下是可以导入 SeqRecord 的指定方式

from Bio.SeqRecord import SeqRecord

让我们逐步使用真实的序列文件理解解析序列文件时的微妙差别。

Parsing Sequence File Formats

本部分阐述如何解析两种最流行的序列文件格式，即 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 格式不支持序列注释。

list comprehension

我们可以使用列表解析将可迭代对象转换为列表，如下所示

>>> 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 提供了一个函数来解析所有生物信息学格式，所以解析 GenBank 格式只需在解析方法中改变 format 选项即可。

以下是相对应的代码：

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