Biopython 序列 I/O 操作


Biopython 提供了一个模块 Bio.SeqIO 来分别从文件(任何流)读取和写入序列。它支持生物信息学中几乎所有可用的文件格式。大多数软件为不同的文件格式提供了不同的方法。但是,Biopython 有意识地遵循一种单一的方法,通过其 SeqRecord 对象将解析的序列数据呈现给用户。

让我们在下一节中了解有关 SeqRecord 的更多信息。

序列记录


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

  • seq : 这是一个实际的序列。

  • id :是给定序列的主标识符。默认类型是字符串。

  • name : 是序列的名称。默认类型是字符串。

  • 描述 : 显示关于序列的人类可读信息。

  • 注释 : 是关于序列的附加信息的字典。

SeqRecord 可以按如下指定导入

from Bio.SeqRecord import SeqRecord

让我们在接下来的部分中了解使用真实序列文件解析序列文件的细微差别。

解析序列文件格式


本节介绍如何解析两种最流行的序列文件格式, FASTA and GenBank .

FASTA


FASTA 是存储序列数据的最基本的文件格式。最初,FASTA 是在生物信息学早期发展过程中开发的用于 DNA 和蛋白质序列比对的软件包,主要用于搜索序列相似性。

Biopython 提供了一个示例 FASTA 文件,可以在以下位置访问它 https://github.com/biopython/biopython/blob/master/Doc/examples/ls_orchid.fasta。

下载此文件并将其保存到你的 Biopython 示例目录中 '兰花.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")

该方法可以有效地转换如下指定的格式:

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 示例目录中 “兰花.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', ...)
    ]
}