序列输入/输出

本章将更详细地讨论 Bio.SeqIO 模块,该模块在第 快速入门 - 你可以用 Biopython 做什么? 章中进行了简要介绍,并且在第 序列注释对象 章中也有使用。该模块旨在提供一个简单的接口,以统一的方式处理各种序列文件格式。另请参见 Bio.SeqIO wiki 页面 (https://biopython.pythonlang.cn/wiki/SeqIO),以及内置文档 Bio.Seq

>>> from Bio import SeqIO
>>> help(SeqIO)

需要注意的是,你需要使用 SeqRecord 对象(参见第 序列注释对象 章),其中包含一个 Seq 对象(参见第 序列对象 章)以及注释,例如标识符和描述。请注意,在处理非常大的 FASTA 或 FASTQ 文件时,使用所有这些对象的开销会导致脚本运行速度过慢。在这种情况下,请考虑使用低级 SimpleFastaParserFastqGeneralIterator 解析器,它们只为每个记录返回一个字符串元组(参见第 低级 FASTA 和 FASTQ 解析器 节)。

解析或读取序列

工作马函数 Bio.SeqIO.parse() 用于将序列数据读入 SeqRecord 对象。此函数需要两个参数

  1. 第一个参数是用于读取数据的句柄,或文件名。句柄通常是为读取而打开的文件,但也可以是来自命令行程序的输出,或从互联网下载的数据(参见第 从网络解析序列 节)。有关句柄的更多信息,请参见第 句柄到底是什么? 节。

  2. 第二个参数是一个小写字符串,指定序列格式 - 我们不会尝试为你猜测文件格式!有关支持的格式的完整列表,请参见 https://biopython.pythonlang.cn/wiki/SeqIO

Bio.SeqIO.parse() 函数返回一个迭代器,该迭代器给出 SeqRecord 对象。迭代器通常在 for 循环中使用,如下所示。

有时你会发现自己处理只包含单个记录的文件。对于这种情况,请使用函数 Bio.SeqIO.read(),它接受相同的参数。如果文件中只有一条记录,则会返回一个 SeqRecord 对象。否则会抛出异常。

读取序列文件

一般来说,Bio.SeqIO.parse() 用于将序列文件读入 SeqRecord 对象,通常与 for 循环一起使用,如下所示

from Bio import SeqIO

for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

上面的例子是从第 解析序列文件格式 节的介绍中重复的,它将以 FASTA 格式文件 ls_orchid.fasta 中的兰花 DNA 序列。如果要加载 GenBank 格式的文件,例如 ls_orchid.gbk,那么你只需要更改文件名和格式字符串

from Bio import SeqIO

for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

同样,如果要读取其他文件格式的文件,那么假设 Bio.SeqIO.parse() 支持该格式,你只需要根据需要更改格式字符串,例如,对于 SwissProt 文件为 "swiss",对于 EMBL 文本文件为 "embl"。在 wiki 页面 (https://biopython.pythonlang.cn/wiki/SeqIO) 和内置文档 Bio.SeqIO 中有完整的列表

另一种使用 Python 迭代器的常见方法是在列表推导(或生成器表达式)中使用。例如,如果你只想从文件中提取记录标识符列表,我们可以使用以下列表推导轻松实现

>>> from Bio import SeqIO
>>> identifiers = [seq_record.id for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank")]
>>> identifiers
['Z78533.1', 'Z78532.1', 'Z78531.1', 'Z78530.1', 'Z78529.1', 'Z78527.1', ..., 'Z78439.1']

在第 序列解析加上简单的绘图 节中,还有更多使用 SeqIO.parse() 在列表推导中进行类似操作的例子(例如,用于绘制序列长度或 GC%)。

遍历序列文件中的记录

在前面的例子中,我们通常使用 for 循环逐个遍历所有记录。你可以使用 for 循环处理所有支持迭代接口的 Python 对象(包括列表、元组和字符串)。

Bio.SeqIO 返回的对象实际上是一个迭代器,它返回 SeqRecord 对象。你可以依次看到每个记录,但只能看到一次。优点是迭代器在处理大型文件时可以节省内存。

除了使用 for 循环,你还可以对迭代器使用 next() 函数来逐步遍历条目,如下所示

from Bio import SeqIO

record_iterator = SeqIO.parse("ls_orchid.fasta", "fasta")

first_record = next(record_iterator)
print(first_record.id)
print(first_record.description)

second_record = next(record_iterator)
print(second_record.id)
print(second_record.description)

请注意,如果你尝试使用 next() 并且没有更多结果,你会得到特殊的 StopIteration 异常。

需要考虑的一种特殊情况是,你的序列文件有多条记录,但你只需要第一条。在这种情况下,以下代码非常简洁

from Bio import SeqIO

first_record = next(SeqIO.parse("ls_orchid.gbk", "genbank"))

这里需要注意的是 - 使用 next() 函数的方式会默默地忽略文件中所有其他记录。如果你的文件只有一条记录,例如本章后面的一些在线示例,或者单个染色体的 GenBank 文件,那么请改用新的 Bio.SeqIO.read() 函数。这将检查是否存在意外的额外记录。

获取序列文件中的记录列表

在上一节中,我们讨论了 Bio.SeqIO.parse() 会返回一个 SeqRecord 迭代器,并且你可以逐个获取记录。很多时候你需要能够以任何顺序访问记录。Python 的 list 数据类型非常适合这种情况,我们可以使用 Python 的内置函数 list() 将记录迭代器转换为 SeqRecord 对象列表,如下所示

from Bio import SeqIO

records = list(SeqIO.parse("ls_orchid.gbk", "genbank"))

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

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

输出结果

Found 94 records
The last record
Z78439.1
Seq('CATTGTTGAGATCACATAATAATTGATCGAGTTAATCTGGAGGATCTGTTTACT...GCC')
592
The first record
Z78533.1
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC')
740

当然,你仍然可以使用 for 循环处理 SeqRecord 对象列表。使用列表比使用迭代器更加灵活(例如,你可以从列表的长度确定记录的数量),但它需要更多的内存,因为它需要将所有记录同时保存在内存中。

提取数据

SeqRecord 对象及其注释结构在第 序列注释对象 章中有更详细的描述。作为注释存储方式的示例,我们将查看解析 GenBank 文件 ls_orchid.gbk 中的第一条记录的结果。

from Bio import SeqIO

record_iterator = SeqIO.parse("ls_orchid.gbk", "genbank")
first_record = next(record_iterator)
print(first_record)

这应该会得到类似的结果

ID: Z78533.1
Name: Z78533
Description: C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA.
Number of features: 5
/sequence_version=1
/source=Cypripedium irapeanum
/taxonomy=['Eukaryota', 'Viridiplantae', 'Streptophyta', ..., 'Cypripedium']
/keywords=['5.8S ribosomal RNA', '5.8S rRNA gene', ..., 'ITS1', 'ITS2']
/references=[...]
/accessions=['Z78533']
/data_file_division=PLN
/date=30-NOV-2006
/organism=Cypripedium irapeanum
/gi=2765658
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC')

这提供了大多数 SeqRecord 注释数据的可读摘要。对于这个例子,我们将使用 .annotations 属性,它只是一个 Python 字典。我们在上面打印记录时显示了该注释字典的内容。你也可以直接打印它们

print(first_record.annotations)

像任何 Python 字典一样,你可以轻松地获取键

print(first_record.annotations.keys())

或值

print(first_record.annotations.values())

一般来说,注释值是字符串,或者字符串列表。一个特殊情况是,文件中任何引用都作为引用对象存储。

假设你想从 ls_orchid.gbk GenBank 文件中提取物种列表。我们想要的信息,Cypripedium irapeanum,保存在注释字典中的 "source" 和 "organism" 下,我们可以像这样访问它们

>>> print(first_record.annotations["source"])
Cypripedium irapeanum

或者

>>> print(first_record.annotations["organism"])
Cypripedium irapeanum

一般来说,“生物体”用于科学名称(拉丁语,例如拟南芥),而“来源”通常是通用名称(例如阿拉伯芥)。在这个例子中,就像通常一样,这两个字段是相同的。

现在让我们遍历所有记录,建立每个兰花序列所属物种的列表。

from Bio import SeqIO

all_species = []
for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank"):
    all_species.append(seq_record.annotations["organism"])
print(all_species)

另一种编写此代码的方法是使用列表推导。

from Bio import SeqIO

all_species = [
    seq_record.annotations["organism"]
    for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank")
]
print(all_species)

无论哪种情况,结果都是

['Cypripedium irapeanum', 'Cypripedium californicum', ..., 'Paphiopedilum barbatum']

太棒了。这很简单,因为 GenBank 文件以标准化方式进行注释。

现在,假设你想从 FASTA 文件中提取物种列表,而不是 GenBank 文件。坏消息是你必须编写一些代码来从记录的描述行中提取你想要的数据——如果信息在文件中!我们示例 FASTA 格式文件ls_orchid.fasta 如下所示

>gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA
CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAACGATCGAGTG
AATCCGGAGGACCGGTGTACTCAGCTCACCGGGGGCATTGCTCCCGTGGTGACCCTGATTTGTTGTTGGG
...

你可以手动检查,但对于每条记录,物种名称都在描述行中作为第二个词。这意味着如果我们将每条记录的.description在空格处拆分,那么物种就在那里作为字段编号 1(字段 0 是记录标识符)。这意味着我们可以这样做

>>> from Bio import SeqIO
>>> all_species = []
>>> for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
...     all_species.append(seq_record.description.split()[1])
...
>>> print(all_species)  
['C.irapeanum', 'C.californicum', 'C.fasciculatum', ..., 'P.barbatum']

使用列表推导的简洁替代方案是

>>> from Bio import SeqIO
>>> all_species = [
...     seq_record.description.split()[1]
...     for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta")
... ]
>>> print(all_species)  
['C.irapeanum', 'C.californicum', 'C.fasciculatum', ..., 'P.barbatum']

总的来说,从 FASTA 描述行中提取信息并不理想。如果你可以将你的序列以 GenBank 或 EMBL 等良好注释的文件格式获取,那么这种类型的注释信息就容易处理得多。

修改数据

在上一节中,我们演示了如何从SeqRecord 中提取数据。另一个常见的任务是修改此数据。SeqRecord 的属性可以直接修改,例如

>>> from Bio import SeqIO
>>> record_iterator = SeqIO.parse("ls_orchid.fasta", "fasta")
>>> first_record = next(record_iterator)
>>> first_record.id
'gi|2765658|emb|Z78533.1|CIZ78533'
>>> first_record.id = "new_id"
>>> first_record.id
'new_id'

注意,如果你想更改 FASTA 在写入文件时的输出方式(参见第Writing Sequence Files 节),那么你应该修改iddescription 属性。为了确保正确行为,最好将id 加上一个空格放在想要的description 的开头。

>>> from Bio import SeqIO
>>> record_iterator = SeqIO.parse("ls_orchid.fasta", "fasta")
>>> first_record = next(record_iterator)
>>> first_record.id = "new_id"
>>> first_record.description = first_record.id + " " + "desired new description"
>>> print(first_record.format("fasta")[:200])
>new_id desired new description
CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGGAATAAA
CGATCGAGTGAATCCGGAGGACCGGTGTACTCAGCTCACCGGGGGCATTGCTCCCGTGGT
GACCCTGATTTGTTGTTGGGCCGCCTCGGGAGCGTCCATGGCGGGT

从压缩文件中解析序列

在上一节中,我们查看了从文件中解析序列数据。你不用使用文件名,而是可以给Bio.SeqIO 一个句柄(参见第What the heck is a handle? 节),在本节中,我们将使用句柄从压缩文件中解析序列。

如你所见,我们可以使用Bio.SeqIO.read()Bio.SeqIO.parse() 以及文件名——例如,这个快速示例使用生成器表达式计算多个记录 GenBank 文件中序列的总长度

>>> from Bio import SeqIO
>>> print(sum(len(r) for r in SeqIO.parse("ls_orchid.gbk", "gb")))
67518

在这里,我们使用文件句柄,使用with 语句自动关闭句柄。

>>> from Bio import SeqIO
>>> with open("ls_orchid.gbk") as handle:
...     print(sum(len(r) for r in SeqIO.parse(handle, "gb")))
...
67518

或者,使用旧式方法手动关闭句柄

>>> from Bio import SeqIO
>>> handle = open("ls_orchid.gbk")
>>> print(sum(len(r) for r in SeqIO.parse(handle, "gb")))
67518
>>> handle.close()

现在,假设我们有一个 gzip 压缩文件?这些文件在 Linux 上非常常见。我们可以使用 Python 的gzip 模块打开压缩文件以供读取——这会给我们一个句柄对象。

>>> import gzip
>>> from Bio import SeqIO
>>> with gzip.open("ls_orchid.gbk.gz", "rt") as handle:
...     print(sum(len(r) for r in SeqIO.parse(handle, "gb")))
...
67518

类似地,如果我们有一个 bzip2 压缩文件

>>> import bz2
>>> from Bio import SeqIO
>>> with bz2.open("ls_orchid.gbk.bz2", "rt") as handle:
...     print(sum(len(r) for r in SeqIO.parse(handle, "gb")))
...
67518

有一个 gzip(GNU Zip)变体称为 BGZF(Blocked GNU Zip Format),它可以像普通 gzip 文件一样进行读取,但对于以后的随机访问具有优势,我们将在第Indexing compressed files 节中讨论。

从网络解析序列

在上一节中,我们查看了从文件(使用文件名或句柄)和从压缩文件(使用句柄)中解析序列数据。在这里,我们将使用另一种类型的句柄,网络连接,以及Bio.SeqIO 来从互联网下载和解析序列。

请注意,即使你可以下载序列数据并将其解析为SeqRecord 对象,但这并不意味着这是一个好主意。总的来说,你可能应该下载序列一次,并将它们保存到一个文件中以备将来使用。

从网络解析 GenBank 记录

EFetch: Downloading full records from Entrez 节详细介绍了 Entrez EFetch 接口,但现在让我们只连接到 NCBI,并使用它们的 GI 号从 GenBank 中获取一些仙人掌(仙人掌)序列。

首先,让我们只获取一条记录。如果你不关心注释和特征,下载 FASTA 文件是一个不错的选择,因为它们很紧凑。现在记住,当你期望句柄只包含一条记录时,请使用Bio.SeqIO.read() 函数

from Bio import Entrez
from Bio import SeqIO

Entrez.email = "[email protected]"
with Entrez.efetch(
    db="nucleotide", rettype="fasta", retmode="text", id="6273291"
) as handle:
    seq_record = SeqIO.read(handle, "fasta")
print("%s with %i features" % (seq_record.id, len(seq_record.features)))

预期输出

gi|6273291|gb|AF191665.1|AF191665 with 0 features

NCBI 还会让你以其他格式索取文件,特别是作为 GenBank 文件。在 2009 年复活节之前,Entrez EFetch API 允许你使用“genbank”作为返回类型,但 NCBI 现在坚持使用官方返回类型“gb”(或蛋白质的“gp”),如EFetch for Sequence and other Molecular Biology Databases 中所述。因此,在 Biopython 1.50 及更高版本中,我们支持“gb”作为Bio.SeqIO 中“genbank”的别名。

from Bio import Entrez
from Bio import SeqIO

Entrez.email = "[email protected]"
with Entrez.efetch(
    db="nucleotide", rettype="gb", retmode="text", id="6273291"
) as handle:
    seq_record = SeqIO.read(handle, "gb")  # using "gb" as an alias for "genbank"
print("%s with %i features" % (seq_record.id, len(seq_record.features)))

此示例的预期输出为

AF191665.1 with 3 features

请注意,这一次我们有三个特征。

现在让我们获取多条记录。这一次,句柄包含多条记录,因此我们必须使用Bio.SeqIO.parse() 函数

from Bio import Entrez
from Bio import SeqIO

Entrez.email = "[email protected]"
with Entrez.efetch(
    db="nucleotide", rettype="gb", retmode="text", id="6273291,6273290,6273289"
) as handle:
    for seq_record in SeqIO.parse(handle, "gb"):
        print("%s %s..." % (seq_record.id, seq_record.description[:50]))
        print(
            "Sequence length %i, %i features, from: %s"
            % (
                len(seq_record),
                len(seq_record.features),
                seq_record.annotations["source"],
            )
        )

这应该会给出以下输出

AF191665.1 Opuntia marenae rpl16 gene; chloroplast gene for c...
Sequence length 902, 3 features, from: chloroplast Opuntia marenae
AF191664.1 Opuntia clavata rpl16 gene; chloroplast gene for c...
Sequence length 899, 3 features, from: chloroplast Grusonia clavata
AF191663.1 Opuntia bradtiana rpl16 gene; chloroplast gene for...
Sequence length 899, 3 features, from: chloroplast Opuntia bradtianaa

请参阅第Accessing NCBI’s Entrez databases 章,详细了解Bio.Entrez 模块,并确保阅读有关使用 Entrez 的 NCBI 指南(第Entrez Guidelines 节)。

从网络解析 SwissProt 序列

现在让我们使用句柄从 ExPASy 下载 SwissProt 文件,这在第Swiss-Prot and ExPASy 章中有更深入的介绍。如上所述,当你期望句柄只包含一条记录时,请使用Bio.SeqIO.read() 函数

from Bio import ExPASy
from Bio import SeqIO

with ExPASy.get_sprot_raw("O23729") as handle:
    seq_record = SeqIO.read(handle, "swiss")
print(seq_record.id)
print(seq_record.name)
print(seq_record.description)
print(repr(seq_record.seq))
print("Length %i" % len(seq_record))
print(seq_record.annotations["keywords"])

假设你的网络连接正常,你应该会得到

O23729
CHS3_BROFI
RecName: Full=Chalcone synthase 3; EC=2.3.1.74; AltName: Full=Naringenin-chalcone synthase 3;
Seq('MAPAMEEIRQAQRAEGPAAVLAIGTSTPPNALYQADYPDYYFRITKSEHLTELK...GAE')
Length 394
['Acyltransferase', 'Flavonoid biosynthesis', 'Transferase']

序列文件作为字典

SeqIO.parse 返回的迭代器循环一次将耗尽文件。对于自索引文件,例如 twoBit 格式的文件,SeqIO.parse 的返回值也可以用作字典,允许随机访问序列内容。在这种情况下,解析是在需要时进行的,因此文件必须保持打开状态,只要序列数据被访问。

>>> from Bio import SeqIO
>>> handle = open("sequence.bigendian.2bit", "rb")
>>> records = SeqIO.parse(handle, "twobit")
>>> records.keys()
dict_keys(['seq11111', 'seq222', 'seq3333', 'seq4', 'seq555', 'seq6'])
>>> records["seq222"]
SeqRecord(seq=Seq('TTGATCGGTGACAAATTTTTTACAAAGAACTGTAGGACTTGCTACTTCTCCCTC...ACA'), id='seq222', name='<unknown name>', description='<unknown description>', dbxrefs=[])
>>> records["seq222"].seq
Seq('TTGATCGGTGACAAATTTTTTACAAAGAACTGTAGGACTTGCTACTTCTCCCTC...ACA')
>>> handle.close()
>>> records["seq222"].seq
Traceback (most recent call last):
...
ValueError: cannot retrieve sequence: file is closed

对于其他文件格式,Bio.SeqIO 提供三个相关的模块函数,允许字典式的随机访问多序列文件。这里在灵活性与内存使用之间存在权衡。总之

  • Bio.SeqIO.to_dict() 选项最灵活,但也最占用内存(参见第Sequence files as Dictionaries – In memory 节)。这基本上是一个辅助函数,用于构建一个普通的 Python dictionary,其中每个条目都作为SeqRecord 对象保存在内存中,允许你修改记录。

  • Bio.SeqIO.index() 选项是一个有用的折衷方案,充当只读字典,并在需要时将序列解析为SeqRecord 对象(参见第Sequence files as Dictionaries – Indexed files 节)。

  • Bio.SeqIO.index_db() 选项也充当只读字典,但将标识符和文件偏移量存储在磁盘上的文件中(作为 SQLite3 数据库),这意味着它对内存的需求非常低(参见第Sequence files as Dictionaries – Database indexed files 节),但速度会稍慢。

请参阅讨论以获取概述(第Discussion 节)。

序列文件作为字典——在内存中

我们对我们无所不在的兰花文件要做的下一件事是展示如何索引它们,并使用 Python dictionary 数据类型(类似于 Perl 中的哈希)像数据库一样访问它们。这对于中等大小的文件非常有用,你只需要访问文件的某些元素,并且可以作为一个快速且便捷的数据库。对于处理内存成为问题的大型文件,请参阅下面的第Sequence files as Dictionaries – Indexed files 节。

你可以使用Bio.SeqIO.to_dict() 函数来创建一个 SeqRecord 字典(在内存中)。默认情况下,这将使用每个记录的标识符(即.id 属性)作为键。让我们尝试使用我们的 GenBank 文件

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.to_dict(SeqIO.parse("ls_orchid.gbk", "genbank"))

Bio.SeqIO.to_dict() 只有一个必需参数,即一个包含SeqRecord 对象的列表或生成器。在这里,我们只是使用了SeqIO.parse 函数的输出。顾名思义,它返回一个 Python 字典。

由于这个变量orchid_dict 是一个普通的 Python 字典,我们可以查看所有可用的键

>>> len(orchid_dict)
94
>>> list(orchid_dict.keys())
['Z78484.1', 'Z78464.1', 'Z78455.1', 'Z78442.1', 'Z78532.1', 'Z78453.1', ..., 'Z78471.1']

在 Python 3 中,像“.keys()“ 和 “.values()“ 这样的字典方法是迭代器,而不是列表。

如果你真的想这样做,你甚至可以一次查看所有记录

>>> list(orchid_dict.values())  # lots of output!

我们可以通过键访问单个SeqRecord 对象,并像往常一样操作该对象

>>> seq_record = orchid_dict["Z78475.1"]
>>> print(seq_record.description)
P.supardii 5.8S rRNA gene and ITS1 and ITS2 DNA
>>> seq_record.seq
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGATCACAT...GGT')

因此,创建 GenBank 记录的内存“数据库”非常容易。接下来,我们将尝试为 FASTA 文件执行此操作。

请注意,对于有 Python 经验的用户来说,手动构建这样的字典应该不成问题。然而,传统的字典构建方法在处理重复键时可能表现不佳。使用 Bio.SeqIO.to_dict() 函数将明确检查重复键,并在发现重复键时抛出异常。

指定字典键

使用与上面相同的代码,但针对 FASTA 文件。

from Bio import SeqIO

orchid_dict = SeqIO.to_dict(SeqIO.parse("ls_orchid.fasta", "fasta"))
print(orchid_dict.keys())

这次的键是

['gi|2765596|emb|Z78471.1|PDZ78471', 'gi|2765646|emb|Z78521.1|CCZ78521', ...
 ..., 'gi|2765613|emb|Z78488.1|PTZ78488', 'gi|2765583|emb|Z78458.1|PHZ78458']

您应该在之前解析 FASTA 文件的章节 简单的 FASTA 解析示例 中识别出这些字符串。假设您更希望使用其他内容作为键,例如登录号。这将我们自然地引向 SeqIO.to_dict() 函数的可选参数 key_function,它允许您定义用于记录的字典键。

首先,您必须编写自己的函数,当给定 SeqRecord 对象时,该函数将返回您想要的键(作为字符串)。一般来说,函数的细节将取决于您处理的输入记录类型。但对于我们的兰花,我们可以使用“管道”字符(垂直线)将记录的标识符拆分,并返回第四个条目(第三个字段)。

def get_accession(record):
    """Given a SeqRecord, return the accession number as a string.

    e.g. "gi|2765613|emb|Z78488.1|PTZ78488" -> "Z78488.1"
    """
    parts = record.id.split("|")
    assert len(parts) == 5 and parts[0] == "gi" and parts[2] == "emb"
    return parts[3]

然后,我们可以将该函数传递给 SeqIO.to_dict() 函数,在构建字典时使用。

from Bio import SeqIO

orchid_dict = SeqIO.to_dict(
    SeqIO.parse("ls_orchid.fasta", "fasta"), key_function=get_accession
)
print(orchid_dict.keys())

最后,如预期,新的字典键

>>> print(orchid_dict.keys())
['Z78484.1', 'Z78464.1', 'Z78455.1', 'Z78442.1', 'Z78532.1', 'Z78453.1', ..., 'Z78471.1']

我希望这并不太复杂!

使用 SEGUID 校验和索引字典

为了提供另一个使用 SeqRecord 对象字典的示例,我们将使用 SEGUID 校验和函数。这是一个相对较新的校验和,并且冲突应该非常罕见(即两个具有相同校验和的不同序列),比 CRC64 校验和有所改进。

再次使用兰花 GenBank 文件。

from Bio import SeqIO
from Bio.SeqUtils.CheckSum import seguid

for record in SeqIO.parse("ls_orchid.gbk", "genbank"):
    print(record.id, seguid(record.seq))

这应该会给出

Z78533.1 JUEoWn6DPhgZ9nAyowsgtoD9TTo
Z78532.1 MN/s0q9zDoCVEEc+k/IFwCNF2pY
...
Z78439.1 H+JfaShya/4yyAj7IbMqgNkxdxQ

现在,回想一下 Bio.SeqIO.to_dict() 函数的 key_function 参数期望一个函数,该函数将 SeqRecord 转换为字符串。我们不能直接使用 seguid() 函数,因为它期望给定 Seq 对象(或字符串)。但是,我们可以使用 Python 的 lambda 特性来创建一个“一次性”函数,以代替传递给 Bio.SeqIO.to_dict()

>>> from Bio import SeqIO
>>> from Bio.SeqUtils.CheckSum import seguid
>>> seguid_dict = SeqIO.to_dict(
...     SeqIO.parse("ls_orchid.gbk", "genbank"), lambda rec: seguid(rec.seq)
... )
>>> record = seguid_dict["MN/s0q9zDoCVEEc+k/IFwCNF2pY"]
>>> print(record.id)
Z78532.1
>>> print(record.description)
C.californicum 5.8S rRNA gene and ITS1 and ITS2 DNA

这应该已经检索到了 Z78532.1 记录,它是文件中的第二个条目。

序列文件作为字典 - 索引文件

正如前面的几个例子试图说明的,使用 Bio.SeqIO.to_dict() 非常灵活。然而,因为它将所有内容都保存在内存中,所以您可以处理的文件大小受限于计算机的 RAM。通常情况下,这只适用于小到中等大小的文件。

对于更大的文件,您应该考虑使用 Bio.SeqIO.index(),它工作方式略有不同。虽然它仍然返回一个类似字典的对象,但它不会所有内容都保存在内存中。相反,它只记录每个记录在文件中的位置 - 当您请求特定的记录时,它会按需解析该记录。

例如,让我们使用之前相同的 GenBank 文件。

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index("ls_orchid.gbk", "genbank")
>>> len(orchid_dict)
94
>>> orchid_dict.keys()
['Z78484.1', 'Z78464.1', 'Z78455.1', 'Z78442.1', 'Z78532.1', 'Z78453.1', ..., 'Z78471.1']
>>> seq_record = orchid_dict["Z78475.1"]
>>> print(seq_record.description)
P.supardii 5.8S rRNA gene and ITS1 and ITS2 DNA
>>> seq_record.seq
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGTTGAGATCACAT...GGT')
>>> orchid_dict.close()

请注意,Bio.SeqIO.index() 不接受句柄,而只接受文件名。这有充分的理由,但有点技术性。第二个参数是文件格式(在其他 Bio.SeqIO 函数中使用的一个小写字符串)。您可以使用许多其他简单文件格式,包括 FASTA 和 FASTQ 文件(请参阅章节 索引 FASTQ 文件 中的示例)。但是,不支持 PHYLIP 或 Clustal 等比对格式。最后,您可以选择提供一个键函数作为可选参数。

以下是用 FASTA 文件的相同示例 - 我们唯一改变的是文件名和格式名称。

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index("ls_orchid.fasta", "fasta")
>>> len(orchid_dict)
94
>>> orchid_dict.keys()
['gi|2765596|emb|Z78471.1|PDZ78471', 'gi|2765646|emb|Z78521.1|CCZ78521', ...
 ..., 'gi|2765613|emb|Z78488.1|PTZ78488', 'gi|2765583|emb|Z78458.1|PHZ78458']

指定字典键

假设您想使用与之前相同的键?就像章节 指定字典键 中使用 Bio.SeqIO.to_dict() 的示例一样,您需要编写一个微小的函数来映射从 FASTA 标识符(作为字符串)到您想要的键。

def get_acc(identifier):
    """Given a SeqRecord identifier string, return the accession number as a string.

    e.g. "gi|2765613|emb|Z78488.1|PTZ78488" -> "Z78488.1"
    """
    parts = identifier.split("|")
    assert len(parts) == 5 and parts[0] == "gi" and parts[2] == "emb"
    return parts[3]

然后,我们可以将该函数传递给 Bio.SeqIO.index() 函数,在构建字典时使用。

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index("ls_orchid.fasta", "fasta", key_function=get_acc)
>>> print(orchid_dict.keys())
['Z78484.1', 'Z78464.1', 'Z78455.1', 'Z78442.1', 'Z78532.1', 'Z78453.1', ..., 'Z78471.1']

知道怎么做之后就很简单了,对吧?

获取记录的原始数据

来自 Bio.SeqIO.index() 的类似字典的对象会将每个条目作为 SeqRecord 对象提供给您。但是,有时能够直接从文件中获取原始数据会很有用。为此,可以使用 get_raw() 方法,该方法接受单个参数(记录标识符)并返回一个字节字符串(从文件中提取,未经修改)。

一个动机示例是从大型文件中提取记录子集,在这种情况下,要么 Bio.SeqIO.write() 不支持输出文件格式(例如纯文本 SwissProt 文件格式),要么您需要完全保留文本(例如,来自 Biopython 的 GenBank 或 EMBL 输出尚未保留所有注释)。

假设您已经从他们的 FTP 站点下载了完整的 UniProt,以纯文本 SwissPort 文件格式存储(ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.dat.gz)并解压缩为 uniprot_sprot.dat 文件,并且您想从中提取几个记录。

>>> from Bio import SeqIO
>>> uniprot = SeqIO.index("uniprot_sprot.dat", "swiss")
>>> with open("selected.dat", "wb") as out_handle:
...     for acc in ["P33487", "P19801", "P13689", "Q8JZQ5", "Q9TRC7"]:
...         out_handle.write(uniprot.get_raw(acc))
...

请注意,在 Python 3 及更高版本中,我们必须以二进制模式打开文件进行写入,因为 get_raw() 方法返回字节字符串。

章节 排序序列文件 中有一个更长的示例,使用 SeqIO.index() 函数对大型序列文件进行排序(无需一次将所有内容都加载到内存中)。

序列文件作为字典 - 数据库索引文件

Biopython 1.57 引入了一种替代方法,Bio.SeqIO.index_db(),它可以处理非常大的文件,因为它将记录信息存储为磁盘上的文件(使用 SQLite3 数据库)而不是存储在内存中。此外,您可以索引多个文件(前提是所有记录标识符都是唯一的)。

Bio.SeqIO.index() 函数接受三个必需参数。

  • 索引文件名,建议使用以 .idx 结尾的文件名。该索引文件实际上是一个 SQLite3 数据库。

  • 要索引的序列文件名列表(或单个文件名)。

  • 文件格式(在 SeqIO 模块的其他部分中使用的小写字符串)。

例如,考虑来自 NCBI FTP 站点 ftp://ftp.ncbi.nih.gov/genbank/ 的 GenBank 平面文件发布版,它们是 gzip 压缩的 GenBank 文件。

截至 GenBank 发布版 \(210\),共有 \(38\) 个文件构成病毒序列,即 gbvrl1.seq、...、gbvrl38.seq,解压缩后在磁盘上大约占 8GB 空间,总共包含近 200 万条记录。

如果您对病毒感兴趣,您可以使用 rsync 命令非常轻松地从命令行下载所有病毒文件,然后使用 gunzip 命令解压缩它们。

# For illustration only, see reduced example below
$ rsync -avP "ftp.ncbi.nih.gov::genbank/gbvrl*.seq.gz" .
$ gunzip gbvrl*.seq.gz

除非您对病毒感兴趣,否则仅为了这个示例而下载如此多的数据 - 因此,让我们只下载前四个部分(每个压缩后约 25MB),并解压缩它们(总共占约 1GB 空间)。

# Reduced example, download only the first four chunks
$ curl -O ftp://ftp.ncbi.nih.gov/genbank/gbvrl1.seq.gz
$ curl -O ftp://ftp.ncbi.nih.gov/genbank/gbvrl2.seq.gz
$ curl -O ftp://ftp.ncbi.nih.gov/genbank/gbvrl3.seq.gz
$ curl -O ftp://ftp.ncbi.nih.gov/genbank/gbvrl4.seq.gz
$ gunzip gbvrl*.seq.gz

现在,在 Python 中,按如下方式索引这些 GenBank 文件。

>>> import glob
>>> from Bio import SeqIO
>>> files = glob.glob("gbvrl*.seq")
>>> print("%i files to index" % len(files))
4
>>> gb_vrl = SeqIO.index_db("gbvrl.idx", files, "genbank")
>>> print("%i sequences indexed" % len(gb_vrl))
272960 sequences indexed

在我的机器上,索引完整的病毒 GenBank 文件集大约需要十分钟,而只索引前四个文件大约需要一分钟。

但是,一旦完成,重复此操作将在几分之一秒内重新加载索引文件 gbvrl.idx

您可以将索引用作只读的 Python 字典 - 无需担心序列来自哪个文件,例如。

>>> print(gb_vrl["AB811634.1"].description)
Equine encephalosis virus NS3 gene, complete cds, isolate: Kimron1.

获取记录的原始数据

就像上面章节 获取记录的原始数据 中讨论的 Bio.SeqIO.index() 函数一样,类似字典的对象还允许您获取每个记录的原始字节。

>>> print(gb_vrl.get_raw("AB811634.1"))
LOCUS       AB811634                 723 bp    RNA     linear   VRL 17-JUN-2015
DEFINITION  Equine encephalosis virus NS3 gene, complete cds, isolate: Kimron1.
ACCESSION   AB811634
...
//

索引压缩文件

当您索引一个序列文件时,它通常会非常大 - 所以您可能希望在磁盘上压缩它。不幸的是,使用 gzip 和 bzip2 等更常见的文件格式很难实现有效的随机访问。在这种情况下,BGZF(块状 GNU 压缩格式)非常有用。它是 gzip 的变体(可以使用标准 gzip 工具解压缩),由 BAM 文件格式、samtoolstabix 推广使用。

要创建 BGZF 压缩文件,您可以使用命令行工具 bgzip,该工具随 samtools 提供。在我们的示例中,我们使用文件名扩展名 *.bgz,以便它们可以与普通 gzip 文件(名为 *.gz)区分开来。您也可以使用 Bio.bgzf 模块从 Python 中读取和写入 BGZF 文件。

Bio.SeqIO.index()Bio.SeqIO.index_db() 都可以与 BGZF 压缩文件一起使用。例如,如果您从一个未压缩的 GenBank 文件开始

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index("ls_orchid.gbk", "genbank")
>>> len(orchid_dict)
94
>>> orchid_dict.close()

您可以在命令行使用以下命令压缩它(同时保留原始文件)——但不用担心,压缩文件已经包含在其他示例文件中

$ bgzip -c ls_orchid.gbk > ls_orchid.gbk.bgz

您可以以完全相同的方式使用压缩文件

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index("ls_orchid.gbk.bgz", "genbank")
>>> len(orchid_dict)
94
>>> orchid_dict.close()

或者

>>> from Bio import SeqIO
>>> orchid_dict = SeqIO.index_db("ls_orchid.gbk.bgz.idx", "ls_orchid.gbk.bgz", "genbank")
>>> len(orchid_dict)
94
>>> orchid_dict.close()

SeqIO 索引会自动检测 BGZF 压缩。请注意,您不能对未压缩文件和压缩文件使用相同的索引文件。

讨论

那么,您应该使用哪种方法以及为什么?这取决于您要做什么(以及您处理多少数据)。但是,总的来说,选择 Bio.SeqIO.index() 是一个好的起点。如果您处理数百万条记录、多个文件或重复分析,那么请查看 Bio.SeqIO.index_db()

选择 Bio.SeqIO.to_dict() 而不是 Bio.SeqIO.index()Bio.SeqIO.index_db() 的原因归结为对灵活性的需求,尽管它需要高内存需求。在内存中存储 SeqRecord 对象的优点是,可以随意更改、添加或删除它们。除了高内存消耗的缺点外,索引也可能需要更长的时间,因为必须完全解析所有记录。

Bio.SeqIO.index()Bio.SeqIO.index_db() 仅按需解析记录。在索引时,它们会扫描文件一次,查找每个记录的开头,并尽可能少地做一些工作来提取标识符。

选择 Bio.SeqIO.index() 而不是 Bio.SeqIO.index_db() 的原因包括

  • 建立索引速度更快(在简单文件格式中更明显)

  • 访问 SeqRecord 对象的速度稍快(但差异只在解析简单文件格式时才真正显着)。

  • 可以使用任何不可变的 Python 对象作为字典键(例如字符串元组或冻结集),而不仅仅是字符串。

  • 如果要索引的序列文件已更改,则无需担心索引数据库过时。

选择 Bio.SeqIO.index_db() 而不是 Bio.SeqIO.index() 的原因包括

  • 不受内存限制——对于二代测序产生的文件,这已经很重要,因为通常会有数千万条序列,使用 Bio.SeqIO.index() 可能需要超过 4GB 的 RAM,因此需要 64 位版本的 Python。

  • 因为索引保存在磁盘上,所以可以重复使用。虽然构建索引数据库文件需要更长的时间,但如果您有一个脚本将在将来对相同的数据文件重新运行,那么从长远来看这可以节省时间。

  • 将多个文件索引在一起

  • get_raw() 方法可以快得多,因为对于大多数文件格式,记录的长度以及偏移量也会被存储。

写入序列文件

我们已经讨论了使用 Bio.SeqIO.parse() 进行序列输入(读取文件),现在我们将看看 Bio.SeqIO.write(),用于序列输出(写入文件)。这是一个接收三个参数的函数:一些 SeqRecord 对象、要写入的句柄或文件名以及序列格式。

以下是一个示例,我们首先以困难的方式(手动,而不是从文件加载)创建一些 SeqRecord 对象

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

rec1 = SeqRecord(
    Seq(
        "MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD"
        "GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK"
        "NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM"
        "SSAC",
    ),
    id="gi|14150838|gb|AAK54648.1|AF376133_1",
    description="chalcone synthase [Cucumis sativus]",
)

rec2 = SeqRecord(
    Seq(
        "YPDYYFRITNREHKAELKEKFQRMCDKSMIKKRYMYLTEEILKENPSMCEYMAPSLDARQ"
        "DMVVVEIPKLGKEAAVKAIKEWGQ",
    ),
    id="gi|13919613|gb|AAK33142.1|",
    description="chalcone synthase [Fragaria vesca subsp. bracteata]",
)

rec3 = SeqRecord(
    Seq(
        "MVTVEEFRRAQCAEGPATVMAIGTATPSNCVDQSTYPDYYFRITNSEHKVELKEKFKRMC"
        "EKSMIKKRYMHLTEEILKENPNICAYMAPSLDARQDIVVVEVPKLGKEAAQKAIKEWGQP"
        "KSKITHLVFCTTSGVDMPGCDYQLTKLLGLRPSVKRFMMYQQGCFAGGTVLRMAKDLAEN"
        "NKGARVLVVCSEITAVTFRGPNDTHLDSLVGQALFGDGAAAVIIGSDPIPEVERPLFELV"
        "SAAQTLLPDSEGAIDGHLREVGLTFHLLKDVPGLISKNIEKSLVEAFQPLGISDWNSLFW"
        "IAHPGGPAILDQVELKLGLKQEKLKATRKVLSNYGNMSSACVLFILDEMRKASAKEGLGT"
        "TGEGLEWGVLFGFGPGLTVETVVLHSVAT",
    ),
    id="gi|13925890|gb|AAK49457.1|",
    description="chalcone synthase [Nicotiana tabacum]",
)

my_records = [rec1, rec2, rec3]

现在我们有了一个 SeqRecord 对象列表,我们将把它们写入 FASTA 格式文件

from Bio import SeqIO

SeqIO.write(my_records, "my_example.faa", "fasta")

如果您在您最喜欢的文本编辑器中打开此文件,它应该看起来像这样

>gi|14150838|gb|AAK54648.1|AF376133_1 chalcone synthase [Cucumis sativus]
MMYQQGCFAGGTVLRLAKDLAENNRGARVLVVCSEITAVTFRGPSETHLDSMVGQALFGD
GAGAVIVGSDPDLSVERPLYELVWTGATLLPDSEGAIDGHLREVGLTFHLLKDVPGLISK
NIEKSLKEAFTPLGISDWNSTFWIAHPGGPAILDQVEAKLGLKEEKMRATREVLSEYGNM
SSAC
>gi|13919613|gb|AAK33142.1| chalcone synthase [Fragaria vesca subsp. bracteata]
YPDYYFRITNREHKAELKEKFQRMCDKSMIKKRYMYLTEEILKENPSMCEYMAPSLDARQ
DMVVVEIPKLGKEAAVKAIKEWGQ
>gi|13925890|gb|AAK49457.1| chalcone synthase [Nicotiana tabacum]
MVTVEEFRRAQCAEGPATVMAIGTATPSNCVDQSTYPDYYFRITNSEHKVELKEKFKRMC
EKSMIKKRYMHLTEEILKENPNICAYMAPSLDARQDIVVVEVPKLGKEAAQKAIKEWGQP
KSKITHLVFCTTSGVDMPGCDYQLTKLLGLRPSVKRFMMYQQGCFAGGTVLRMAKDLAEN
NKGARVLVVCSEITAVTFRGPNDTHLDSLVGQALFGDGAAAVIIGSDPIPEVERPLFELV
SAAQTLLPDSEGAIDGHLREVGLTFHLLKDVPGLISKNIEKSLVEAFQPLGISDWNSLFW
IAHPGGPAILDQVELKLGLKQEKLKATRKVLSNYGNMSSACVLFILDEMRKASAKEGLGT
TGEGLEWGVLFGFGPGLTVETVVLHSVAT

假设您想知道 Bio.SeqIO.write() 函数向句柄写入多少条记录?如果您的记录在列表中,您只需使用 len(my_records),但是当您的记录来自生成器/迭代器时,您无法这样做。 Bio.SeqIO.write() 函数返回写入文件的 SeqRecord 对象的数量。

注意 - 如果您告诉 Bio.SeqIO.write() 函数写入一个已经存在的文件,则旧文件将被覆盖,没有任何警告。

往返

有些人喜欢他们的解析器是“可往返的”,意思是如果您读取一个文件并将其写回,它将保持不变。这要求解析器必须提取足够的信息来精确地重现原始文件。 Bio.SeqIO打算这样做。

作为一个简单的例子,FASTA 文件中序列数据的任何换行都是允许的。从解析以下两个示例(它们只在换行符方面不同)将得到相同的 SeqRecord

>YAL068C-7235.2170 Putative promoter sequence
TACGAGAATAATTTCTCATCATCCAGCTTTAACACAAAATTCGCACAGTTTTCGTTAAGA
GAACTTAACATTTTCTTATGACGTAAATGAAGTTTATATATAAATTTCCTTTTTATTGGA

>YAL068C-7235.2170 Putative promoter sequence
TACGAGAATAATTTCTCATCATCCAGCTTTAACACAAAATTCGCA
CAGTTTTCGTTAAGAGAACTTAACATTTTCTTATGACGTAAATGA
AGTTTATATATAAATTTCCTTTTTATTGGA

为了创建一个可往返的 FASTA 解析器,您需要跟踪序列换行符的位置,而这些额外信息通常毫无意义。相反,Biopython 在输出时使用默认换行符 \(60\) 个字符。许多其他文件格式也存在同样的空格问题。在某些情况下,另一个问题是 Biopython 还没有(尚未)保留最后一个注释(例如 GenBank 和 EMBL)。

有时保留原始布局(包括它可能具有的任何怪癖)很重要。请参见第 获取记录的原始数据 节关于 Bio.SeqIO.index() 字典类对象中 get_raw() 方法的一个潜在解决方案。

在序列文件格式之间转换

在前面的示例中,我们使用了一个 SeqRecord 对象列表作为 Bio.SeqIO.write() 函数的输入,但它也接受 SeqRecord 迭代器,就像我们从 Bio.SeqIO.parse() 获得的一样——这让我们可以通过组合这两个函数来进行文件转换。

对于此示例,我们将读取 GenBank 格式文件 ls_orchid.gbk 并将其写入 FASTA 格式

from Bio import SeqIO

records = SeqIO.parse("ls_orchid.gbk", "genbank")
count = SeqIO.write(records, "my_example.fasta", "fasta")
print("Converted %i records" % count)

不过,这有点复杂。因此,由于文件转换是一项非常常见的任务,因此有一个帮助函数可以让你用以下内容替换它

from Bio import SeqIO

count = SeqIO.convert("ls_orchid.gbk", "genbank", "my_example.fasta", "fasta")
print("Converted %i records" % count)

Bio.SeqIO.convert() 函数接受句柄文件名。注意,如果输出文件已经存在,它将被覆盖!要了解更多信息,请查看内置帮助

>>> from Bio import SeqIO
>>> help(SeqIO.convert)

原则上,只需更改文件名和格式名,这段代码就可以用于在 Biopython 中提供的任何文件格式之间进行转换。但是,写入某些格式需要其他文件格式不包含的信息(例如质量得分)。例如,虽然您可以将 FASTQ 文件转换为 FASTA 文件,但您无法反过来操作。另请参见第 转换 FASTQ 文件 和第 将 FASTA 和 QUAL 文件转换为 FASTQ 文件 节,它们介绍了在不同 FASTQ 格式之间进行转换的示例。

最后,除了您的代码将更短之外,使用 Bio.SeqIO.convert() 函数的另一个好处是,这样做可能更快!这样做的原因是,convert 函数可以利用一些特定于文件格式的优化和技巧。

将序列文件转换为它们的互补反向序列

假设您有一个核苷酸序列文件,您想将其转换为包含其互补反向序列的文件。这一次,需要做一些工作才能将从输入文件获得的 SeqRecord 对象转换为适合保存到输出文件的对象。

首先,我们将使用 Bio.SeqIO.parse() 从文件加载一些核苷酸序列,然后使用 Seq 对象内置的 .reverse_complement() 方法(请参见第 核苷酸序列和(反向)互补序列 节)打印它们的互补反向序列

>>> from Bio import SeqIO
>>> for record in SeqIO.parse("ls_orchid.gbk", "genbank"):
...     print(record.id)
...     print(record.seq.reverse_complement())
...

现在,如果我们要将这些互补反向序列保存到文件,我们需要创建 SeqRecord 对象。我们可以使用 SeqRecord 对象内置的 .reverse_complement() 方法(请参见第 反向互补 SeqRecord 对象 节),但我们必须决定如何命名新的记录。

这是一个展示列表推导功能的绝佳地方,它可以在内存中创建一个列表

>>> from Bio import SeqIO
>>> records = [
...     rec.reverse_complement(id="rc_" + rec.id, description="reverse complement")
...     for rec in SeqIO.parse("ls_orchid.fasta", "fasta")
... ]
>>> len(records)
94

现在列表推导有一个很好的技巧,您可以在其中添加一个条件语句

>>> records = [
...     rec.reverse_complement(id="rc_" + rec.id, description="reverse complement")
...     for rec in SeqIO.parse("ls_orchid.fasta", "fasta")
...     if len(rec) < 700
... ]
>>> len(records)
18

这将在内存中创建一个互补反向记录列表,其中序列长度小于 700 个碱基对。但是,我们可以用生成器表达式做完全相同的事情——但优势在于它不会一次在内存中创建所有记录的列表

>>> records = (
...     rec.reverse_complement(id="rc_" + rec.id, description="reverse complement")
...     for rec in SeqIO.parse("ls_orchid.fasta", "fasta")
...     if len(rec) < 700
... )

作为一个完整的例子

>>> from Bio import SeqIO
>>> records = (
...     rec.reverse_complement(id="rc_" + rec.id, description="reverse complement")
...     for rec in SeqIO.parse("ls_orchid.fasta", "fasta")
...     if len(rec) < 700
... )
>>> SeqIO.write(records, "rev_comp.fasta", "fasta")
18

翻译 CDS 条目的 FASTA 文件 节中有一个相关的示例,它将 FASTA 文件中的每个记录从核苷酸翻译成氨基酸。

获取您的 SeqRecord 对象作为格式化的字符串

假设您并不想将记录写入文件或句柄——而是想要一个包含以特定文件格式表示的记录的字符串。 Bio.SeqIO 接口基于句柄,但 Python 有一个有用的内置模块,它提供了一个基于字符串的句柄。

例如,我们可以从我们的兰花 GenBank 文件中加载一批 SeqRecord 对象,并创建一个包含 FASTA 格式记录的字符串。

from Bio import SeqIO
from io import StringIO

records = SeqIO.parse("ls_orchid.gbk", "genbank")
out_handle = StringIO()
SeqIO.write(records, out_handle, "fasta")
fasta_data = out_handle.getvalue()
print(fasta_data)

第一次看到这个可能不太直观!不过,对于想要在特定文件格式中获取包含单个记录的字符串的特殊情况,可以使用 SeqRecord 类的 format() 方法(参见第 格式方法 部分)。

需要注意的是,虽然我们不鼓励这样做,但你可以使用 format() 方法写入文件,例如:

from Bio import SeqIO

with open("ls_orchid_long.tab", "w") as out_handle:
    for record in SeqIO.parse("ls_orchid.gbk", "genbank"):
        if len(record) > 100:
            out_handle.write(record.format("tab"))

虽然这种代码风格适用于像 FASTA 这样的简单顺序文件格式,或者这里使用的简单制表符分隔格式,但它不适用于更复杂或交织的文件格式。这就是我们仍然推荐使用 Bio.SeqIO.write() 的原因,如下例所示:

from Bio import SeqIO

records = (rec for rec in SeqIO.parse("ls_orchid.gbk", "genbank") if len(rec) > 100)
SeqIO.write(records, "ls_orchid.tab", "tab")

与多次调用 SeqRecord.format(...) 方法相比,一次调用 SeqIO.write(...) 速度要快得多。

低级 FASTA 和 FASTQ 解析器

在处理大型高通量 FASTA 或 FASTQ 测序文件时,使用低级的 SimpleFastaParserFastqGeneralIterator 通常比使用 Bio.SeqIO.parse 更实用,因为速度很重要。如本章引言中所述,文件格式中立的 Bio.SeqIO 接口需要为像 FASTA 这样的简单格式创建很多对象,因此会带来一些开销。

在解析 FASTA 文件时,Bio.SeqIO.parse() 会在内部使用文件句柄调用低级的 SimpleFastaParser。你可以直接使用它,它会遍历文件句柄,并返回每个记录作为一个包含两个字符串的元组,即标题行(> 字符后的所有内容)和序列(作为纯字符串)。

>>> from Bio.SeqIO.FastaIO import SimpleFastaParser
>>> count = 0
>>> total_len = 0
>>> with open("ls_orchid.fasta") as in_handle:
...     for title, seq in SimpleFastaParser(in_handle):
...         count += 1
...         total_len += len(seq)
...
>>> print("%i records with total sequence length %i" % (count, total_len))
94 records with total sequence length 67518

只要你不关心行换行(对于短读高通量数据,你可能并不关心),那么从这些字符串中输出 FASTA 格式也会非常快。

...
out_handle.write(">%s\n%s\n" % (title, seq))
...

类似地,在解析 FASTQ 文件时,Bio.SeqIO.parse() 会在内部使用文件句柄调用低级的 FastqGeneralIterator。如果你不需要将质量分数转换为整数,或者可以将它们作为 ASCII 字符串处理,那么这种方法非常理想。

>>> from Bio.SeqIO.QualityIO import FastqGeneralIterator
>>> count = 0
>>> total_len = 0
>>> with open("example.fastq") as in_handle:
...     for title, seq, qual in FastqGeneralIterator(in_handle):
...         count += 1
...         total_len += len(seq)
...
>>> print("%i records with total sequence length %i" % (count, total_len))
3 records with total sequence length 75

食谱 (第 食谱 – 使用它的酷方法 章) 中有更多关于这方面的示例,包括如何使用这段代码片段从字符串中有效地输出 FASTQ 格式。

...
out_handle.write("@%s\n%s\n+\n%s\n" % (title, seq, qual))
...