用于计算分子生物学的 Python 工具

Biopython 版本 1.84
© 2024。保留所有权利。
在 GitHub 上编辑此页面

SeqIO 简介

本页介绍 Bio.SeqIO，它是 Biopython 1.43 及更高版本中用于标准序列输入/输出的接口。有关实现细节，请参阅 SeqIO 开发页面。

Python 初学者可能会发现 Peter 的介绍性 Biopython 工作坊 有用，该工作坊从使用 SeqIO 处理序列文件开始。

在 教程 (PDF) 中有一整章内容介绍了 Bio.SeqIO，虽然有些内容重复，但除了本 WIKI 页面之外，非常值得一读。此外，还有 API 文档（您可以在线阅读，或者使用 help 命令在 Python 中阅读）。

目标

Bio.SeqIO 提供了一个简单统一的接口，用于输入和输出各种序列文件格式（包括多序列比对），但它 *只* 处理作为 SeqRecord 对象的序列。还有一个姊妹接口 Bio.AlignIO，用于直接处理作为比对对象的序列比对文件。

该设计部分受到 BioPerl 的 SeqIO 的简单性的启发。从长远来看，我们希望匹配 BioPerl 支持的 序列文件格式 和 多序列比对格式 的令人印象深刻的列表。

请注意，在 Biopython 中包含 Bio.SeqIO（和 Bio.AlignIO）会导致在处理某些文件格式时出现一些重复或选择。例如，Bio.Nexus 也会从 Nexus 文件中读取序列 - 但 Bio.Nexus 还可以做更多的事情，例如读取 Nexus 文件中的任何系统发育树。

我的愿景是，对于操作序列数据，您应该尝试将 Bio.SeqIO 作为您的首选。除非您有一些非常特殊的要求，我希望这应该足够了。

Peter

文件格式

此表列出了 Bio.SeqIO 可以读取、写入和索引的文件格式，以及首次支持该格式的 Biopython 版本（或 git 表示这在最新的开发代码中受支持）。格式名称是一个简单的 lowercase 字符串。在可能的情况下，我们使用与 BioPerl 的 SeqIO 和 EMBOSS 相同的名称。

格式名称 读取 写入 索引 备注
abi 1.58 否 N/A 读取 ABI“Sanger”毛细管序列跟踪文件，包括碱基调用 PHRED 质量分数。这允许 ABI 到 FASTQ 的转换。请注意，每个 ABI 文件只包含一个序列（因此对文件进行索引毫无意义）。
abi-trim 1.71 否 N/A 与“abi”相同，但使用 Mott 算法进行质量修剪。
ace 1.47 否 1.52 读取 ACE 组装文件中 contigs 的序列。在内部使用 Bio.Sequencing.Ace。
cif-atom 1.73 否 否 使用 Bio.PDB.MMCIFParser 根据原子坐标确定结构中出现的（部分）蛋白质序列。
cif-seqres 1.73 否 否 读取宏观分子晶体学信息文件 (mmCIF) 文件以确定由 _pdbx_poly_seq_scheme 记录定义的完整蛋白质序列。
clustal 1.43 1.43 否 Clustal X 和 Clustal W 的比对格式。
embl 1.43 1.54 1.52 EMBL 平面文件格式。在内部使用 Bio.GenBank 进行解析。Biopython 1.48 到 1.50 只编写具有最小注释的基本 GenBank 文件，而 1.51 及更高版本还将编写特征表。
fasta 1.43 1.43 1.52 这指的是为 Bill Pearson 的 FASTA 工具引入的 FASTA 文件格式 输入，其中每个记录以“>”行开头。
fasta-2line 1.71 1.71 否 FASTA 格式变体，没有换行符，每个记录正好有两行。
fastq-sanger 或 fastq 1.50 1.50 1.52 FASTQ 文件 类似于 FASTA 文件，但还包括测序质量。在 Biopython 中，“fastq”（或别名“fastq-sanger”）指的是 Sanger 风格的 FASTQ 文件，它们使用 33 的 ASCII 偏移量编码 PHRED 质量。另见早期 Solexa/Illumina 管道中使用的不兼容的“fastq-solexa”和“fastq-illumina”变体，Illumina 管道 1.8 生成了 Sanger FASTQ。
fastq-solexa 1.50 1.50 1.52 在 Biopython 中，“fastq-solexa”指的是原始的 Solexa/Illumina 风格的 FASTQ 文件，它们使用 64 的 ASCII 偏移量编码 Solexa 质量。另见我们称为“fastq-illumina”格式的内容。
fastq-illumina 1.51 1.51 1.52 在 Biopython 中，“fastq-illumina”指的是早期 Solexa/Illumina 风格的 FASTQ 文件（来自管道版本 1.3 到 1.7），它们使用 64 的 ASCII 偏移量编码 PHRED 质量。对于 *良好* 质量的读数，PHRED 和 Solexa 分数大致相等，因此“fastq-solexa”和“fastq-illumina”变体几乎等效。
gck 1.75 否 否 Gene Construction Kit 使用的本机格式。
genbank 或 gb 1.43 1.48 / 1.51 1.52 GenBank 或 GenPept 平面文件格式。在内部使用 Bio.GenBank 进行解析。Biopython 1.48 到 1.50 只编写具有最小注释的基本 GenBank 文件，而 1.51 及更高版本还将编写特征表。
ig 1.47 否 1.52 这指的是 IntelliGenetics 文件格式，显然与 MASE 比对格式 相同。
imgt 1.56 1.56 1.56 这指的是 IMGT 版本的 EMBL 纯文本文件格式。
nexus 1.43 1.48 否 NEXUS 多序列比对 格式，也称为 PAUP 格式。在内部使用 Bio.Nexus。
pdb-seqres 1.61 否 否 读取蛋白质数据银行 (PDB) 文件以确定在标题中出现的完整蛋白质序列（不依赖于 Bio.PDB 和 NumPy）。
pdb-atom 1.61 否 否 使用 Bio.PDB 根据文件的原子坐标部分确定结构中出现的（部分）蛋白质序列（需要 NumPy）。
phd 1.46 1.52 1.52 PHD 文件 是 PHRED 的输出，由 PHRAP 和 CONSED 用于输入。在内部使用 Bio.Sequencing.Phd。
phylip 1.43 1.43 否 PHYLIP 文件。将名称截断为 10 个字符。
pir 1.48 1.71 1.52 由国家生物医学研究基金会 (NBRF) 为 蛋白质信息资源 (PIR) 数据库引入的“类似 FASTA”格式，现在是 UniProt 的一部分。
seqxml 1.58 1.58 否 简单的 序列 XML 文件格式。
sff 1.54 1.54 1.54 由 Roche 454 和 IonTorrent/IonProton 测序仪生成的标准流图格式 (SFF) 二进制文件。
sff-trim 1.54 否 1.54 应用文件中列出的修剪的标准流图格式。
snapgene 1.75 否 否 SnapGene 使用的本机格式。
stockholm 1.43 1.43 否 斯德哥尔摩比对格式 也称为 PFAM 格式。
swiss 1.43 否 1.52 Swiss-Prot 也称为 UniProt 格式。在内部使用 Bio.SwissProt。另见 UniProt XML 格式。
tab 1.48 1.48 1.52 简单的两列制表符分隔的序列文件，其中每一行包含一个记录的标识符和序列。例如，这在 Aligent 的 eArray 软件将微阵列探针保存到最小的制表符分隔文本文件中时使用。
qual 1.50 1.50 1.52 Qual 文件 类似于 FASTA 文件，但不是序列，而是记录空格分隔的整数测序值作为 PHRED 质量分数。匹配的 FASTA 和 QUAL 文件对通常用作单个 FASTQ 文件的替代方案。
uniprot-xml 1.56 否 1.56 UniProt XML 格式，是纯文本 Swiss-Prot 格式的继任者。
xdna 1.75 1.75 否 Christian Marck 的 DNA Strider 和 Serial Cloner 使用的本机格式。
         

使用 Bio.SeqIO，您可以将序列比对文件格式视为任何其他序列文件，但新的 Bio.AlignIO 模块旨在直接处理此类比对文件。您还可以将来自任何文件格式的一组 SeqRecord 对象转换为比对 - 只要它们都具有相同的长度。请注意，当使用 Bio.SeqIO 将序列写入比对文件格式时，所有（间隙）序列都应该具有相同的长度。

序列输入

主要功能是 Bio.SeqIO.parse()，它接受一个文件句柄（或文件名）和格式名称，并返回一个 SeqRecord 迭代器。这使您可以执行以下操作

from Bio import SeqIO

for record in SeqIO.parse("example.fasta", "fasta"):
    print(record.id)

或者使用句柄

from Bio import SeqIO

with open("example.fasta") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        print(record.id)

在上面的示例中，我们使用内置的 python 函数 open 打开了文件。with 语句确保文件在读取后正确关闭。如果您只使用文件名，则所有这些都应该自动发生。

请注意，您 *必须* 显式指定文件格式，与 BioPerl 的 SeqIO 不同，后者可以尝试使用文件名扩展名和/或文件内容来猜测。参见 *显式优于隐式* (Python 之禅)。

如果您有不同类型的文件，例如 Clustalw 比对文件，例如 opuntia.aln，其中包含七个序列，唯一区别在于您指定 "clustal" 而不是 "fasta"

from Bio import SeqIO

with open("opuntia.aln") as handle:
    for record in SeqIO.parse(handle, "clustal"):
        print(record.id)

迭代器非常适合您只需要按文件中的顺序逐个获取记录的情况。对于某些任务，您可能需要以任何顺序随机访问记录。在这种情况下，使用内置的 python list() 函数将迭代器转换为列表

from Bio import SeqIO

records = list(SeqIO.parse("example.fasta", "fasta"))
print(records[0].id)  # first record
print(records[-1].id)  # last record

另一个常见任务是按某个标识符索引您的记录。对于小文件，我们有一个函数 Bio.SeqIO.to_dict() 将 SeqRecord 迭代器（或列表）转换为字典（在内存中）

from Bio import SeqIO

record_dict = SeqIO.to_dict(SeqIO.parse("example.fasta", "fasta"))
print(record_dict["gi:12345678"])  # use any record ID

函数 Bio.SeqIO.to_dict() 默认情况下将使用记录 ID 作为字典键，但您可以使用可选参数 key_function 指定您喜欢的任何映射。

对于较大的文件，无法将所有内容保存在内存中，因此 Bio.SeqIO.to_dict 不适用。Biopython 1.52 及更高版本包含 Bio.SeqIO.index 函数来处理这种情况，但您也可以考虑使用 BioSQL.

from Bio import SeqIO

record_dict = SeqIO.index("example.fasta", "fasta")
print(record_dict["gi:12345678"])  # use any record ID

Biopython 1.45 引入了另一个函数，Bio.SeqIO.read()，它与 Bio.SeqIO.parse() 一样，需要一个句柄和格式。它用于句柄中只包含一个记录的情况，该记录将作为单个 SeqRecord 对象返回。如果没有任何记录或超过一个记录，则会引发异常。

from Bio import SeqIO

record = SeqIO.read("single.fasta", "fasta")

对于您只想获取第一条记录（并且愿意忽略任何后续记录）的相关情况，您可以使用内置的 Python 函数 next。

from Bio import SeqIO

first_record = next(SeqIO.parse("example.fasta", "fasta"))

序列输出

要将记录写入文件，请使用函数 Bio.SeqIO.write()，它接受一个 SeqRecord 迭代器（或列表）、输出句柄（或文件名）和格式字符串。

from Bio import SeqIO

sequences = ...  # add code here
with open("example.fasta", "w") as output_handle:
    SeqIO.write(sequences, output_handle, "fasta")

或者

from Bio import SeqIO

sequences = ...  # add code here
SeqIO.write(sequences, "example.fasta", "fasta")

在以下关于文件格式转换的部分中，有更多示例。

请注意，如果您正在写入对齐文件格式，则所有序列必须具有相同的长度。

如果您将序列作为 SeqRecord 迭代器提供，则对于 Fasta 或 GenBank 等顺序文件格式，可以逐个写入记录。因为一次只创建一个记录，所以需要的内存很少。请参阅下面的示例，该示例过滤了一组记录。

另一方面，对于 Clustal 等交错或非顺序文件格式，Bio.SeqIO.write() 函数将被迫自动将迭代器转换为列表。这将破坏使用生成器/迭代器方法节省内存的任何可能性。

文件格式转换

假设您有一个 GenBank 文件，要将其转换为 Fasta 文件。例如，让我们考虑文件 cor6_6.gb，该文件包含在 Biopython 单元测试的 GenBank 目录下。

您可以使用 Bio.SeqIO.parse() 函数像这样读取文件

from Bio import SeqIO

with open("cor6_6.gb") as input_handle:
    for record in SeqIO.parse(input_handle, "genbank"):
        print(record)

请注意，此文件包含六条记录。现在，让我们将 SeqRecord 迭代器传递给 Bio.SeqIO.write() 函数，而不是打印记录，以便将此 GenBank 文件转换为 Fasta 文件

from Bio import SeqIO

with open("cor6_6.gb") as input_handle, open(
    "cor6_6.fasta", "w"
) as output_handle:
    sequences = SeqIO.parse(input_handle, "genbank")
    count = SeqIO.write(sequences, output_handle, "fasta")

print("Converted %i records" % count)

或者，更简洁地说，使用 Bio.SeqIO.convert() 函数（在 Biopython 1.52 或更高版本中），只需

from Bio import SeqIO

count = SeqIO.convert("cor6_6.gb", "genbank", "cor6_6.fasta", "fasta")
print("Converted %i records" % count)

在此示例中，GenBank 文件的开头是这样的

LOCUS       ATCOR66M      513 bp    mRNA            PLN       02-MAR-1992
DEFINITION  A.thaliana cor6.6 mRNA.
ACCESSION   X55053
VERSION     X55053.1  GI:16229
...

生成的 Fasta 文件如下所示

>X55053.1 A.thaliana cor6.6 mRNA.
AACAAAACACACATCAAAAACGATTTTACAAGAAAAAAATA...
...

请注意，所有 Fasta 文件可以存储的只是标识符、描述和序列。

通过更改格式字符串，该代码可用于在任何受支持的文件格式之间进行转换。

示例

输入/输出示例 - 按序列长度过滤

虽然您可能只是想转换一个文件（如上所示），但更现实的例子是对数据进行某种操作或过滤。

例如，让我们将所有长度小于 300 个核苷酸的“短”序列保存到 Fasta 文件中

from Bio import SeqIO

short_sequences = []  # Setup an empty list
for record in SeqIO.parse("cor6_6.gb", "genbank"):
    if len(record.seq) < 300:
        # Add this record to our list
        short_sequences.append(record)

print("Found %i short sequences" % len(short_sequences))

SeqIO.write(short_sequences, "short_seqs.fasta", "fasta")

如果您了解 **列表推导**，那么您也可以这样编写上面的示例

from Bio import SeqIO

input_seq_iterator = SeqIO.parse("cor6_6.gb", "genbank")

# Build a list of short sequences:
short_sequences = [record for record in input_seq_iterator if len(record.seq) < 300]

print("Found %i short sequences" % len(short_sequences))

SeqIO.write(short_sequences, "short_seqs.fasta", "fasta")

我不确定这是否真的更容易理解，但它更短。

但是，如果您正在处理具有数千条记录的非常大的文件，那么您可以从使用 **生成器表达式** 中受益。这避免在内存中创建整个所需记录列表

from Bio import SeqIO

input_seq_iterator = SeqIO.parse("cor6_6.gb", "genbank")
short_seq_iterator = (record for record in input_seq_iterator if len(record.seq) < 300)

SeqIO.write(short_seq_iterator, "short_seqs.fasta", "fasta")

请记住，对于 Fasta 或 GenBank 等顺序文件格式，Bio.SeqIO.write() 将接受一个 SeqRecord 迭代器。上面代码的优势在于，一次只会有一个记录在内存中。

但是，正如输出部分所述，对于 Clustal 等非顺序文件格式，Bio.SeqIO.write() 被迫自动将迭代器转换为列表，因此这种优势就消失了。

如果这一切都令人困惑，不要惊慌，只需忽略这些花哨的东西。对于中等大小的数据集，在内存中一次拥有太多记录（例如在列表中）可能不会出现问题。

使用 SEGUID 校验和

在此示例中，我们将使用 Bio.SeqIO 和 Bio.SeqUtils.CheckSum 模块（在 Biopython 1.44 或更高版本中）。首先，我们将只打印出 GenBank 文件 ls_orchid.gbk 中每个序列的校验和

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()，使用 SEGUID 作为键来构建一个 SeqRecord 字典。这里的诀窍是使用 Python lambda 语法创建一个临时函数来获取每个 SeqRecord 的 SEGUID - 我们不能直接使用 seguid() 函数，因为它只适用于 Seq 对象或字符串。

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)
print(record.description)

给出以下输出

Z78439.1
P.barbatum 5.8S rRNA gene and ITS1 and ITS2 DNA.

随机子序列

此脚本将读取一个包含整个线粒体基因组的 Genbank 文件（例如，烟草线粒体，Nicotiana tabacum mitochondrion NC_006581），创建 500 条记录，包含该基因组的随机片段，并将它们保存为 fasta 文件。这些子序列是使用随机起点和固定的 200 长度创建的。

from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
from random import randint

# There should be one and only one record, the entire genome:
mito_record = SeqIO.read("NC_006581.gbk", "genbank")

mito_frags = []
limit = len(mito_record.seq)
for i in range(0, 500):
    start = randint(0, limit - 200)
    end = start + 200
    mito_frag = mito_record.seq[start:end]
    record = SeqRecord(mito_frag, "fragment_%i" % (i + 1), "", "")
    mito_frags.append(record)

SeqIO.write(mito_frags, "mitofrags.fasta", "fasta")

这应该会将以下内容作为输出文件，

>fragment_1
TGGGCCTCATATTTATCCTATATACCATGTTCGTATGGTGGCGCGATGTTCTACGTGAAT
CCACGTTCGAAGGACATCATACCAAAGTCGTACAATTAGGACCTCGATATGGTTTTATTC
TGTTTATCGTATCGGAGGTTATGTTCTTTTTTGCTCTTTTTCGGGCTTCTTCTCATTCTT
CTTTGGCACCTACGGTAGAG
...
>fragment_500
ACCCAGTGCCGCTACCCACTTCTACTAAGGCTGAGCTTAATAGGAGCAAGAGACTTGGAG
GCAACAACCAGAATGAAATATTATTTAATCGTGGAAATGCCATGTCAGGCGCACCTATCA
GAATCGGAACAGACCAATTACCAGATCCACCTATCATCGCCGGCATAACCATAAAAAAGA
TCATTAAAAAAGCGTGAGCC

写入字符串

有时您可能不希望将 SeqRecord 对象写入文件，而是写入字符串。例如，您可能正在准备输出，以便将其显示为网页的一部分。如果要将多个记录写入单个字符串，请使用 StringIO 创建一个基于字符串的句柄。教程（PDF）在 SeqIO 章节中对此进行了示例。

对于您希望将单个记录作为给定文件格式的字符串的特殊情况，Biopython 1.48 添加了一个新的格式方法

from Bio import SeqIO

mito_record = SeqIO.read("NC_006581.gbk", "genbank")
print(mito_record.format("fasta"))

格式方法将接受 Bio.SeqIO 支持的任何输出格式，其中文件格式可用于单个记录（例如 "fasta"、"tab" 或 "genbank"）。

请注意，我们不建议您将其用于文件输出 - 使用 Bio.SeqIO.write() 速度更快且更通用。

帮助！

如果您在使用 Bio.SeqIO 时遇到问题，请加入讨论邮件列表（请参阅 邮件列表）。

如果您认为您发现了一个错误，请在项目的 GitHub 页面 上报告它。