序列输入/输出

在本章中，我们将更详细地讨论 Bio.SeqIO 模块，在第七章中简要介绍了该模块  快速入门-您可以用Biopython做什么？ 也用于章节  顺序注释对象 .这旨在提供一个简单的界面，以便以统一的方式处理各种序列文件格式。另见 Bio.SeqIO 维基页面（http：//biopython.org/wiki/SeqIO）和内置文档 Bio.Seq :

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

关键是你必须 SeqRecord 对象（参见第章  顺序注释对象 ），其中包含 Seq 对象（请参阅第章  序列对象 ）加上注释，如标识符和描述。请注意，在处理非常大的FASTA或FASTQ文件时，处理所有这些对象的开销可能会使脚本太慢。在这种情况下，考虑低级别 SimpleFastaParser 和 FastqGeneralIterator 解析器仅返回每条记录的字符串数组（请参阅第节  低级FASTA和FASTQ解析器 ).

解析或读取序列

主力功能 Bio.SeqIO.parse() 用于将序列数据作为SeqRecord对象读取。此函数需要两个参数：

1. 第一个参数是一个 handle 读取数据或文件名。手柄通常是打开以读取的文件，但可以是命令行程序的输出，也可以是从互联网下载的数据（请参阅第节  从网上解析序列 ).见章节  把手到底是什么？ 有关手柄的更多信息。

2. 第二个参数是指定序列格式的大小写字符串-我们不会尝试为您猜测文件格式！有关支持格式的完整列表，请参阅http://biopython.org/wiki/SeqIO。

的 Bio.SeqIO.parse() 函数返回一个 iterator 这给 SeqRecord 对象迭代器通常用于for循环，如下所示。

有时您会发现自己处理的文件只包含一条记录。对于这种情况，使用该功能 Bio.SeqIO.read() 这也采取了同样的论点。如果文件中只有一条记录，则该记录将作为 SeqRecord object.否则将引发异常。

读取序列文件

一般来说 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))

上面的例子在第节的介绍中重复  解析序列文件格式 ，并将兰花DNA序列加载到FASTA格式文件中 ls_orchid.fasta .如果您想加载像这样的基因库格式文件 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”。维基页面（http：//biopython.org/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() 在第节中类似的列表理解中  序列解析加上简单的情节 (e.g.用于绘制序列长度或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() 这样的函数将默默忽略文件中的任何其他记录。如果您的文件有 one and only one 记录，如本章后面的一些在线示例，或单个染色体的基因库文件，然后使用新的 Bio.SeqIO.read() 相反，功能。这将检查是否存在额外的意外记录。

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

在上一节中，我们讨论了这样一个事实： Bio.SeqIO.parse() 给你一个 SeqRecord 迭代器，并且您逐个获取记录。通常您需要能够以任何顺序访问记录。蟒蛇 list 数据类型非常适合这一点，我们可以将记录迭代器转换为 SeqRecord 使用内置Python函数的对象 list() 就像这样：

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 object and its annotation structures are described more fully in Chapter 顺序注释对象. As an example of how annotations are stored, we’ll look at the output from parsing the first record in the GenBank file 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 ，保存在注释词典中的“源”和“组织”下，我们可以这样访问：

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

或者：

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

一般来说，“有机体”用于学名（拉丁语，例如 Arabidopsis thaliana ），而“source”通常是常用名称（例如thale cress）。在本例中，与通常的情况一样，两个字段是相同的。

现在让我们回顾所有记录，建立每个兰花序列所属物种的列表：

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']

太好了这非常容易，因为基因库文件是以标准化的方式注释的。

现在，假设您想从FASTA文件（而不是家谱文件）中提取物种列表。坏消息是，您必须编写一些代码才能从记录的描述行中提取您想要的数据-如果信息首先就在文件中的话！我们的示例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写入文件时的输出方式（请参阅第节  写入序列文件 ），那么您应该修改 id 和 description 美德.先知-愿为了确保正确的行为，最好包括 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 手柄（请参阅部分  把手到底是什么？ ），在本节中，我们将使用手柄来解析压缩文件中的序列。

正如您在上面所看到的，我们可以使用 Bio.SeqIO.read() 或 Bio.SeqIO.parse() 带有文件名-例如，这个快速示例使用生成器公式计算多记录ESB文件中序列的总长度：

>>> 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 模块用于打开压缩文件进行读取-这为我们提供了一个handle对象：

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

有一个名为BGZR（阻塞的NU Zip格式）的gZip（NU Zip）变体，它可以像普通gZip文件一样进行读取，但具有随机访问的优势，稍后我们将在第一节中讨论  索引压缩文件 .

从网上解析序列

在前面的部分中，我们研究了从文件（使用文件名或手柄）和从压缩文件（使用手柄）中解析序列数据。在这里我们将使用 Bio.SeqIO 使用另一种类型的手柄（网络连接）来从互联网下载和解析序列。

请注意，仅仅因为你 can 下载序列数据并将其解析为 SeqRecord 一次性对象并不意味着这是一个好主意。一般来说，您可能应该下载序列 once 并将它们保存到文件中以供重复使用。

从网上解析基因库记录

部分  EFetch：从Deliverz下载完整记录 更详细地讨论了Deliverz EFetch接口，但现在让我们连接到NCBI并获取一些 Opuntia （花椒）序列，使用其GI编号来自基因库。

首先，让我们只获取一条记录。如果您不关心注释和功能，下载FASTA文件是一个不错的选择，因为这些文件很紧凑。现在请记住，当您希望该手柄包含且仅包含一条记录时，请使用 Bio.SeqIO.read() 功能：

from Bio import Entrez
from Bio import SeqIO

Entrez.email = "A.N.Other@example.com"
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还允许您请求其他格式的文件，特别是作为SEN文件。直到2009年复活节，Deliverz EFetch API允许您使用“genbank”作为返回类型，但NCBI现在坚持使用官方返回类型“gi”（或蛋白质的“GP”），如所述 EFetch for Sequence and other Molecular Biology Databases .因此，在Biopython 1.50以后，我们支持将“GB”作为中“genbank”的别名 Bio.SeqIO .

from Bio import Entrez
from Bio import SeqIO

Entrez.email = "A.N.Other@example.com"
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 = "A.N.Other@example.com"
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

参见Chapter  NCBI的 欲了解更多关于 Bio.Entrez 模块，并确保阅读NCBI使用Deliverz的指南（第部分  埃斯珀兹准则 ).

从网络中解析SwissProt序列

现在，让我们使用手柄从ExPAS下载SwissProt文件，第章将更深入地介绍这一点  Swiss-Prot和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() 是最灵活但内存要求最高的选项（请参阅第节  将文件序列为词典-在内存中 ).这基本上是构建普通Python的助手函数 dictionary 每个条目都作为 SeqRecord 内存中的对象，允许您修改记录。

• Bio.SeqIO.index() 是一个有用的中间立场，就像只读字典一样，将序列解析为 SeqRecord 按需对象（请参阅部分  将文件序列为词典-索引文件 ).

• Bio.SeqIO.index_db() 也像只读字典一样，但将标识符和文件补偿存储在磁盘上的文件中（作为SQLite3数据库），这意味着它的内存需求非常低（请参阅第节  作为词典的序列文件-数据库索引文件 ），但会慢一点。

请参阅讨论以了解广泛的概述（第节  讨论 ).

将文件序列为词典-在内存中

我们要对无处不在的兰花文件做的下一步是展示如何对它们进行索引并像使用Python访问数据库一样访问它们 dictionary 数据类型（就像Perl中的哈希）。这对于中等大的文件非常有用，在这些文件中您只需要访问文件的某些元素，并且可以形成一个快速且肮脏的数据库。有关处理内存成为问题的较大文件，请参阅部分  将文件序列为词典-索引文件 下面

您可以使用函数 Bio.SeqIO.to_dict() 制作SeqRecord字典（在内存中）。默认情况下，这将使用每个记录的标识符（即 .id 属性）作为关键。让我们使用我们的基因库文件尝试一下：

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

因此，创建我们的基因库记录的内存“数据库”非常容易。接下来，我们将尝试对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 object.一般来说，功能的详细信息取决于您正在处理的输入记录的类型。但对于我们的兰花，我们可以使用“管道”字符（垂直线）拆分记录的标识符并返回第四个条目（字段三）：

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']

我希望不要太复杂！

使用SEARCH检查和对字典进行索引

再举一个使用词典的例子 SeqRecord 对象时，我们将使用SESYS检查和函数。这是一个相对较新的检验和，冲突应该非常罕见（即两个不同的序列具有相同的检验和），这是对CRC 64检验和的改进。

再次使用兰花基因库文件：

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() ，其工作方式略有不同。尽管它仍然返回类似对象的字典，但这确实 not 保持 everything 在记忆中。相反，它只是记录每条记录在文件中的位置-当您请求特定记录时，它会根据需要解析它。

例如，让我们使用与之前相同的SEN文件：

>>> 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或Custal等对齐格式。最后，作为可选参数，您可以提供一个关键功能。

这是使用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 object.然而，能够直接从文件中获取原始原始数据有时很有用。对于此用途， get_raw() 方法，该方法采用单个参数（记录标识符）并返回字节字符串（未经修改从文件中提取）。

一个激励的例子是从一个大文件中提取一个记录的子集，其中 Bio.SeqIO.write() 不支持输出文件格式（例如纯文本SwissProt文件格式）或您需要精确保留文本的地方（例如Biopython的GenBank或EMBL输出尚未保留注释的每一位）。

假设您已经从他们的RTP网站（ftp：//ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/compete/uniprot_sprot.dat.gz）下载了纯文本SwissPort文件格式的整个UniProt，并将其解压为文件 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 TP网站ftp：//ftp.ncbi.nih.gov/genbank/的SEN平面文件，它们是gZip压缩的ESB文件。

截至基因库发布 \(210\) ，有 \(38\) 构成病毒序列的文件， gbvrl1.seq , …, gbvrl38.seq 解压后，磁盘上约有8 GB，总共包含近200万条记录。

如果您对病毒感兴趣，可以使用 rsync 命令，然后将它们与 gunzip :

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

除非您关心病毒，否则仅在这个示例中就需要下载大量数据-所以让我们下载 just 前四个块（每个压缩约25 MB），并初始化它们（占用约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

在我的机器上对全套病毒基因库文件进行索引大约花了十分钟，仅前四个文件就花了大约一分钟左右。

但是，一旦完成，重复此操作将重新加载索引文件 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和bzip 2等更常见的文件格式，有效的随机访问很困难。在这种设置中，BGZR（阻止的NU Zip格式）可能非常有帮助。这是gZip的一个变体（可以使用标准gZip工具解压），由BAM文件格式普及， samtools ，而且 tabix .

要创建BGZR压缩文件，您可以使用命令行工具 bgzip 它附带samtools。在我们的示例中，我们使用文件名扩展名 *.bgz ，因此它们可以与普通的gzipped文件（名为 *.gz ).您还可以使用 Bio.bgzf 模块用于在Python中读写BGZR文件。

的 Bio.SeqIO.index() 和 Bio.SeqIO.index_db() 两者都可以与BGZR压缩文件一起使用。例如，如果您从未压缩的SEN文件开始：

>>> 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 索引自动检测BGZR压缩。请注意，您不能对未压缩和压缩的文件使用相同的索引文件。

讨论

那么，您应该使用这些方法中的哪一种以及为什么？这取决于您试图做什么（以及您正在处理多少数据）。不过，一般采摘 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 写入文件的对象。

Note - 如果你告诉 Bio.SeqIO.write() 函数写入已存在的文件，旧文件将被覆盖，而没有任何警告。

往返

有些人喜欢他们的解析器是“可舍入的”，这意味着如果您读取文件并再次将其写回，它不会改变。这要求解析器必须提取足够的信息来重现原始文件 exactly . Bio.SeqIO 并 not 旨在做到这一点。

作为一个简单的例子，允许对FASTA文件中的序列数据进行任何行包装。相同的 SeqRecord 将通过解析以下两个示例给出，这两个示例仅在换行符上有所不同：

>YAL068C-7235.2170 Putative promoter sequence
TACGAGAATAATTTCTCATCATCCAGCTTTAACACAAAATTCGCACAGTTTTCGTTAAGA
GAACTTAACATTTTCTTATGACGTAAATGAAGTTTATATATAAATTTCCTTTTTATTGGA

>YAL068C-7235.2170 Putative promoter sequence
TACGAGAATAATTTCTCATCATCCAGCTTTAACACAAAATTCGCA
CAGTTTTCGTTAAGAGAACTTAACATTTTCTTATGACGTAAATGA
AGTTTATATATAAATTTCCTTTTTATTGGA

要制作可舍入三重的FASTA解析器，您需要跟踪序列行中断发生的位置，而这些额外的信息通常毫无意义。相反，Biopython使用默认的行包装 \(60\) 输出上的字符。空白的相同问题也适用于许多其他文件格式。在某些情况下，另一个问题是Biopython（尚未）保留注释的每一个最后一点（例如，SEN和MBE）。

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

序列文件格式之间转换

在前面的例子中，我们使用了一个 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() 功能将采用手柄 or 文件名。但请注意-如果输出文件已经存在，它将覆盖它！要了解更多信息，请参阅内置帮助：

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

原则上，只需更改文件名和格式名称，即可使用此代码在Biopython中可用的任何文件格式之间进行转换。然而，编写某些格式需要其他文件格式不包含的信息（例如质量分数）。例如，虽然您可以将FASTQ文件转换为FASTA文件，但您不能进行相反的操作。另请参阅部分  转换FASTQ文件 和  将FASTA和QUAL文件转换为FASTQ文件 在食谱章节中，该章节探讨了不同FASTQ格式之间的相互转换。

最后，作为使用 Bio.SeqIO.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有一个有用的内置模块，提供基于字符串的手柄。

对于如何使用它的示例，让我们加载一堆 SeqRecord 从我们的orchids基因库文件中获取对象，并创建一个包含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)

当您第一次看到它时，这并不完全简单！从好的方面来说，对于您想要包含 single 记录，请使用 SeqRecord 类' format() 方法（请参阅部分  format方法 ).

请注意，虽然我们不鼓励这样做，但您 can 使用 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等简单的顺序文件格式或这里使用的简单制表符分隔格式，但它将 not 适用于更复杂或交错的文件格式。这就是为什么我们仍然建议使用 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")

打一个电话 SeqIO.write(...) 也比多次拨打 SeqRecord.format(...) 法

低级FASTA和FASTQ解析器

与底层合作 SimpleFastaParser 或 FastqGeneralIterator 往往比实际 Bio.SeqIO.parse 在处理速度很重要的大型高通量FASTA或FASTQ测序文件时。正如本章介绍中所指出的，文件格式中立 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 带有文件手柄。如果您不需要将质量分数转换为integer，或者可以将它们作为ASC字符串处理，这是理想的选择：

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