食谱-很酷的事情与它有关
Biopython现在有两个“食谱”示例集--本章（已包含在本教程中多年，并且逐渐增长）和http://biopython.org/wiki/Category:Cookbook，这是我们维基上的用户贡献集。
我们试图鼓励Biopython用户向维基贡献自己的示例。除了帮助社区之外，分享这样的示例的一个直接好处是，您还可以从其他Biopython用户和开发人员那里获得一些关于代码的反馈-这可以帮助您改进所有Python代码。
从长远来看，我们最终可能会将本章中的所有示例移至维基或教程中的其他地方。
使用序列文件
本节使用 Bio.SeqIO 第章中描述的模块 序列输入/输出 .
过滤序列文件
通常，您会有一个包含许多序列的大文件（例如FASTA文件或基因，或读段的FASTQ或SFF文件）、感兴趣序列子集的单独较短的ID列表，并且想要为此子集创建新的序列文件。
假设ID列表位于一个简单的文本文件中，作为每一行的第一个单词。这可能是一个表格文件，其中第一列是ID。尝试类似这样的操作：
from Bio import SeqIO
input_file = "big_file.sff"
id_file = "short_list.txt"
output_file = "short_list.sff"
with open(id_file) as id_handle:
wanted = set(line.rstrip("\n").split(None, 1)[0] for line in id_handle)
print("Found %i unique identifiers in %s" % (len(wanted), id_file))
records = (r for r in SeqIO.parse(input_file, "sff") if r.id in wanted)
count = SeqIO.write(records, output_file, "sff")
print("Saved %i records from %s to %s" % (count, input_file, output_file))
if count < len(wanted):
print("Warning %i IDs not found in %s" % (len(wanted) - count, input_file))
请注意，我们使用Python set 而不是 list ，这使得测试会员资格的速度更快。
按照章节所述 低级FASTA和FASTQ解析器 ，对于大型FASTA或FASTQ文件，为了速度，您最好不使用高级 SeqIO 接口，但直接使用字符串。下一个示例展示了如何使用FASTQ文件执行此操作-它更复杂：
from Bio.SeqIO.QualityIO import FastqGeneralIterator
input_file = "big_file.fastq"
id_file = "short_list.txt"
output_file = "short_list.fastq"
with open(id_file) as id_handle:
# Taking first word on each line as an identifier
wanted = set(line.rstrip("\n").split(None, 1)[0] for line in id_handle)
print("Found %i unique identifiers in %s" % (len(wanted), id_file))
with open(input_file) as in_handle:
with open(output_file, "w") as out_handle:
for title, seq, qual in FastqGeneralIterator(in_handle):
# The ID is the first word in the title line (after the @ sign):
if title.split(None, 1)[0] in wanted:
# This produces a standard 4-line FASTQ entry:
out_handle.write("@%s\n%s\n+\n%s\n" % (title, seq, qual))
count += 1
print("Saved %i records from %s to %s" % (count, input_file, output_file))
if count < len(wanted):
print("Warning %i IDs not found in %s" % (len(wanted) - count, input_file))
产生随机基因组
假设您正在查看基因组序列，寻找一些序列特征-可能是极端的局部GC%偏差，或可能的限制性消化位点。一旦您的Python代码在真实的基因组上运行，尝试对同一基因组的随机版本运行相同的搜索以进行统计分析可能是明智的（毕竟，您发现的任何“特征”都可能只是偶然存在）。
在本讨论中，我们将使用来自 Yersinia pestis biovar Microtus .该文件包含在Gene文件夹下的Biopython单元测试中，或者您可以从我们的网站获取它， NC_005816.gb .该文件包含且仅包含一条记录，因此我们可以将其作为 SeqRecord 使用 Bio.SeqIO.read() 功能：
>>> from Bio import SeqIO
>>> original_rec = SeqIO.read("NC_005816.gb", "genbank")
那么，我们如何生成原始序列的洗牌版本呢？我会使用内置Python random 为此模块，特别是功能 random.shuffle - 但这适用于Python列表。我们的序列是 Seq 对象，因此为了对它进行洗牌，我们需要将它转换成一个列表：
>>> import random
>>> nuc_list = list(original_rec.seq)
>>> random.shuffle(nuc_list) # acts in situ!
现在，为了使用 Bio.SeqIO 要输出混洗序列，我们需要构建一个新的 SeqRecord 用新的 Seq 使用此洗牌列表的对象。为了做到这一点，我们需要将核苷酸列表（单字母字符串）变成一个长字符串-Python的标准方法是使用字符串对象的join方法。
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> shuffled_rec = SeqRecord(
... Seq("".join(nuc_list)), id="Shuffled", description="Based on %s" % original_rec.id
... )
让我们将所有这些部分放在一起，创建一个完整的Python脚本，该脚本生成一个FASTA文件，其中包含原始序列的30个随机洗牌版本。
这个第一个版本只是使用一个大的for循环并逐个写出记录（使用 SeqRecord 部分中描述的格式方法 将SeqRecord对象作为格式化字符串 ):
import random
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO
original_rec = SeqIO.read("NC_005816.gb", "genbank")
with open("shuffled.fasta", "w") as output_handle:
for i in range(30):
nuc_list = list(original_rec.seq)
random.shuffle(nuc_list)
shuffled_rec = SeqRecord(
Seq("".join(nuc_list)),
id="Shuffled%i" % (i + 1),
description="Based on %s" % original_rec.id,
)
output_handle.write(shuffled_rec.format("fasta"))
就我个人而言，我更喜欢下面的版本，它使用函数来洗牌记录和生成器表达，而不是for循环：
import random
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO
def make_shuffle_record(record, new_id):
nuc_list = list(record.seq)
random.shuffle(nuc_list)
return SeqRecord(
Seq("".join(nuc_list)),
id=new_id,
description="Based on %s" % original_rec.id,
)
original_rec = SeqIO.read("NC_005816.gb", "genbank")
shuffled_recs = (
make_shuffle_record(original_rec, "Shuffled%i" % (i + 1)) for i in range(30)
)
SeqIO.write(shuffled_recs, "shuffled.fasta", "fasta")
翻译CDS条目的FASTA文件
假设您有一个包含某些生物体CDS条目的输入文件，并且您想要生成一个包含其蛋白质序列的新FASTA文件。即从原始文件中取出每个核苷序列，并翻译它。回到部分 翻译 我们了解了如何使用 Seq 对象的 translate method ，以及可选的 cds 能够正确翻译替代开始密码子的参数。
我们可以将其与 Bio.SeqIO 如第2.1节中的反码示例所示 将序列文件转换为反向互补序列 .关键是对于每个核苷酸 SeqRecord ，我们需要创造一种蛋白质 SeqRecord - 并负责命名它。
您可以编写自己的函数来做到这一点，为您的序列选择合适的蛋白质标识符以及适当的遗传密码。在本例中，我们只需使用默认表并向标识符添加一个前置符：
from Bio.SeqRecord import SeqRecord
def make_protein_record(nuc_record):
"""Returns a new SeqRecord with the translated sequence (default table)."""
return SeqRecord(
seq=nuc_record.seq.translate(cds=True),
id="trans_" + nuc_record.id,
description="translation of CDS, using default table",
)
然后我们可以使用这个函数将输入的核苷酸记录转化为准备输出的蛋白质记录。一种优雅且内存高效的方法是使用生成器表达：
from Bio import SeqIO
proteins = (
make_protein_record(nuc_rec)
for nuc_rec in SeqIO.parse("coding_sequences.fasta", "fasta")
)
SeqIO.write(proteins, "translations.fasta", "fasta")
这应该适用于任何完整编码序列的FASTA文件。如果您正在研究部分编码序列，您可能更愿意使用 nuc_record.seq.translate(to_stop=True) 在上面的例子中，因为这不会检查有效的起始密码子等。
在FASTA文件大写字母中制作序列
通常，您会以FASTA文件的形式从合作者处获取数据，有时序列可能是大写和大写的混合形式。在某些情况下，这是故意的（例如，质量较差的地区的情况较低），但通常这并不重要。您可能想要编辑文件以使所有内容一致（例如全大写），并且您可以使用 upper() 方法 SeqRecord 对象（添加到Biopython 1.55中）：
from Bio import SeqIO
records = (rec.upper() for rec in SeqIO.parse("mixed.fas", "fasta"))
count = SeqIO.write(records, "upper.fas", "fasta")
print("Converted %i records to upper case" % count)
这是怎么回事？第一行只是导入 Bio.SeqIO module.第二行是有趣的一点-这是一个Python生成器表达，它给出了从输入文件解析的每个记录的大写版本 (mixed.fas ).在第三行中，我们将这个生成器表达赋予 Bio.SeqIO.write() 函数并将新的大写记录保存到我们的输出文件 (upper.fas ).
我们使用生成器表达（而不是列表或列表理解）的原因是，这意味着内存中一次只保留一条记录。如果您正在处理具有数百万个条目的大型文件，这一点非常重要。
排序序列文件
假设您想按长度对序列文件进行排序（例如来自程序集的一组重叠群），并且您正在使用FASTA或FASTQ等文件格式， Bio.SeqIO 可以读取、写入（和索引）。
如果文件足够小，您可以将其作为列表一次将其全部加载到内存中 SeqRecord 对象，对列表进行排序并保存：
from Bio import SeqIO
records = list(SeqIO.parse("ls_orchid.fasta", "fasta"))
records.sort(key=lambda r: len(r))
SeqIO.write(records, "sorted_orchids.fasta", "fasta")
唯一聪明的地方是指定如何对记录进行排序的比较方法（这里我们按长度对它们进行排序）。如果您想要首先获得最长的记录，则可以翻转比较或使用相反的参数：
from Bio import SeqIO
records = list(SeqIO.parse("ls_orchid.fasta", "fasta"))
records.sort(key=lambda r: -len(r))
SeqIO.write(records, "sorted_orchids.fasta", "fasta")
现在这是非常直接的-但是如果你有一个非常大的文件，你不能像这样把它全部加载到内存中，会发生什么呢？例如，您可能有一些下一代测序读取要按长度排序。这可以使用 Bio.SeqIO.index() 功能
from Bio import SeqIO
# Get the lengths and ids, and sort on length
len_and_ids = sorted(
(len(rec), rec.id) for rec in SeqIO.parse("ls_orchid.fasta", "fasta")
)
ids = reversed([id for (length, id) in len_and_ids])
del len_and_ids # free this memory
record_index = SeqIO.index("ls_orchid.fasta", "fasta")
records = (record_index[id] for id in ids)
SeqIO.write(records, "sorted.fasta", "fasta")
首先，我们扫描文件，使用 Bio.SeqIO.parse() ，将记录标识符及其长度记录在二元组列表中。然后我们对这个列表进行排序，以按照长度顺序获取它们，并丢弃长度。使用此排序的标识符列表 Bio.SeqIO.index() 允许我们一一检索记录，并将它们传递给 Bio.SeqIO.write() 用于输出。
这些例子都使用 Bio.SeqIO 将记录解析为 SeqRecord 使用输出的对象 Bio.SeqIO.write() .如果您想对以下文件格式进行排序该怎么办 Bio.SeqIO.write() 不支持纯文本SwissProt格式？这是使用 get_raw() 添加的方法 Bio.SeqIO.index() 在Biopython 1.54中（请参阅部分 获取原始数据以记录 ).
from Bio import SeqIO
# Get the lengths and ids, and sort on length
len_and_ids = sorted(
(len(rec), rec.id) for rec in SeqIO.parse("ls_orchid.fasta", "fasta")
)
ids = reversed([id for (length, id) in len_and_ids])
del len_and_ids # free this memory
record_index = SeqIO.index("ls_orchid.fasta", "fasta")
with open("sorted.fasta", "wb") as out_handle:
for id in ids:
out_handle.write(record_index.get_raw(id))
请注意，对于Python 3以后，我们必须打开文件才能以二进制模式进行写入，因为 get_raw() 方法返回 bytes 对象
作为奖励，因为它不会将数据解析为 SeqRecord 第二次对象应该更快。如果您只想将其与FASTA格式一起使用，我们可以通过使用低级FASTA解析器获取记录标识符和长度来进一步加快速度：
from Bio.SeqIO.FastaIO import SimpleFastaParser
from Bio import SeqIO
# Get the lengths and ids, and sort on length
with open("ls_orchid.fasta") as in_handle:
len_and_ids = sorted(
(len(seq), title.split(None, 1)[0])
for title, seq in SimpleFastaParser(in_handle)
)
ids = reversed([id for (length, id) in len_and_ids])
del len_and_ids # free this memory
record_index = SeqIO.index("ls_orchid.fasta", "fasta")
with open("sorted.fasta", "wb") as out_handle:
for id in ids:
out_handle.write(record_index.get_raw(id))
FASTQ文件的简单质量过滤
FASTQ文件格式是在Sanger引入的，现在被广泛用于保存核酸测序读数及其质量评分。FASTQ文件（以及相关的QUAL文件）是按字母注释的一个很好的例子，因为序列中的每个核苷酸都有相关的质量评分。任何按字母计算的注释都保存在 SeqRecord 在 letter_annotations 字典作为列表、tuple或字符串（元素数量与序列长度相同）。
一项常见的任务是获取一大组测序读段并根据其质量评分对其进行过滤（或裁剪）。下面的例子非常简单化，但应该说明在 SeqRecord object.我们在这里要做的就是读取FASTQ数据文件，并对其进行过滤，以仅选出PHRED质量分数均高于某个阈值（此处为20）的记录。
对于本示例，我们将使用从ENA序列读取存档ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz（2 MB）下载的一些真实数据，该存档解压为19 MB文件 SRR020192.fastq .这是来自感染病毒的加州海狮的Roche 454 GS FLX单端数据（详细信息，请参阅https://www.ebi.ac.uk/ena/data/view/SRS004476）。
首先，让我们计算阅读量：
from Bio import SeqIO
count = 0
for rec in SeqIO.parse("SRR020192.fastq", "fastq"):
count += 1
print("%i reads" % count)
现在让我们对最低PHRED质量为20进行简单的过滤：
from Bio import SeqIO
good_reads = (
rec
for rec in SeqIO.parse("SRR020192.fastq", "fastq")
if min(rec.letter_annotations["phred_quality"]) >= 20
)
count = SeqIO.write(good_reads, "good_quality.fastq", "fastq")
print("Saved %i reads" % count)
这只是拔出来的 \(14580\) 读出 \(41892\) 礼物更明智的做法是对读数进行质量修剪，但这仅作为示例。
FASTQ文件可以包含数百万个条目，因此最好避免将它们同时加载到内存中。此示例使用生成器表达，这意味着只有一个 SeqRecord 一次创建-避免任何内存限制。
请注意，使用低级会更快 FastqGeneralIterator 此处的解析器（请参阅部分 低级FASTA和FASTQ解析器 ），但这不会将质量字符串转化为整分。
修剪碱基序列
对于这个例子，我们将假装 GATGACGGTGT 是我们想要在一些FASTQ格式的读段数据中寻找的5 '碱基序列。与上面的例子一样，我们将使用 SRR020192.fastq 从ENA下载的文件（ftp：//ftp.sra.ebi.ac.uk/vol1/fastq/SRR 020/SRR 020192/SRR 020192.fastq.gz）。
通过使用主 Bio.SeqIO 界面时，相同的方法适用于任何其他支持的文件格式（例如FASTA文件）。然而，对于大型FASTQ文件，它会比低级文件更快 FastqGeneralIterator 此处的解析器（请参阅前面的示例和小节 低级FASTA和FASTQ解析器 ).
此代码使用 Bio.SeqIO 使用生成器表达（以避免将所有序列同时加载到内存中），以及 Seq 对象的 startswith 方法来查看读取是否以碱基序列开始：
from Bio import SeqIO
primer_reads = (
rec
for rec in SeqIO.parse("SRR020192.fastq", "fastq")
if rec.seq.startswith("GATGACGGTGT")
)
count = SeqIO.write(primer_reads, "with_primer.fastq", "fastq")
print("Saved %i reads" % count)
那应该找到 \(13819\) 读取 SRR014849.fastq 并将它们保存到新的FASTQ文件中， with_primer.fastq .
现在假设您想制作一个包含这些读段但删除了碱基序列的FASTQ文件？这只是一个小变化，因为我们可以切片 SeqRecord (see部分 切片SeqRecord ）删除前十一个字母（我们入门的长度）：
from Bio import SeqIO
trimmed_primer_reads = (
rec[11:]
for rec in SeqIO.parse("SRR020192.fastq", "fastq")
if rec.seq.startswith("GATGACGGTGT")
)
count = SeqIO.write(trimmed_primer_reads, "with_primer_trimmed.fastq", "fastq")
print("Saved %i reads" % count)
再说一次，这应该拉出 \(13819\) 读取 SRR020192.fastq ，但这次删除了前十个字符，并将它们保存到另一个新的FASTQ文件中， with_primer_trimmed.fastq .
现在，假设您想要创建一个新的FASTQ文件，其中这些读段的primer已被删除，但所有其他读段仍保持原样？如果我们想仍然使用生成器表达，最清楚的方法是定义我们自己的修剪函数：
from Bio import SeqIO
def trim_primer(record, primer):
if record.seq.startswith(primer):
return record[len(primer) :]
else:
return record
trimmed_reads = (
trim_primer(record, "GATGACGGTGT")
for record in SeqIO.parse("SRR020192.fastq", "fastq")
)
count = SeqIO.write(trimmed_reads, "trimmed.fastq", "fastq")
print("Saved %i reads" % count)
这需要更长的时间，因为这次输出文件包含所有 \(41892\) 阅读。同样，我们使用生成器表达来避免任何内存问题。您也可以使用生成器函数而不是生成器表达。
from Bio import SeqIO
def trim_primers(records, primer):
"""Removes perfect primer sequences at start of reads.
This is a generator function, the records argument should
be a list or iterator returning SeqRecord objects.
"""
len_primer = len(primer) # cache this for later
for record in records:
if record.seq.startswith(primer):
yield record[len_primer:]
else:
yield record
original_reads = SeqIO.parse("SRR020192.fastq", "fastq")
trimmed_reads = trim_primers(original_reads, "GATGACGGTGT")
count = SeqIO.write(trimmed_reads, "trimmed.fastq", "fastq")
print("Saved %i reads" % count)
如果您想执行更复杂的操作，其中仅保留部分记录，则此表格更加灵活-如下一个示例中所示。
修剪接头序列
这本质上是对前一个示例的简单扩展。我们要假装 GATGACGGTGT 是一些FASTQ格式的读数据中的适配器序列，同样是 SRR020192.fastq 来自NCBI的文件（ftp：//ftp.sra.ebi.ac.uk/vol1/fastq/SRR 020/SRR 020192/SRR 020192.fastq.gz）。
然而，这次我们将寻找序列 anywhere 在阅读中，不仅仅是在一开始：
from Bio import SeqIO
def trim_adaptors(records, adaptor):
"""Trims perfect adaptor sequences.
This is a generator function, the records argument should
be a list or iterator returning SeqRecord objects.
"""
len_adaptor = len(adaptor) # cache this for later
for record in records:
index = record.seq.find(adaptor)
if index == -1:
# adaptor not found, so won't trim
yield record
else:
# trim off the adaptor
yield record[index + len_adaptor :]
original_reads = SeqIO.parse("SRR020192.fastq", "fastq")
trimmed_reads = trim_adaptors(original_reads, "GATGACGGTGT")
count = SeqIO.write(trimmed_reads, "trimmed.fastq", "fastq")
print("Saved %i reads" % count)
由于我们在本例中使用了FASTQ输入文件，因此 SeqRecord 对象具有质量分数的按字母注释。进行切片加工而 SeqRecord 对象，在修剪的记录上使用适当的分数，因此我们也可以将它们输出为FASTQ文件。
与前一个例子的输出（我们只在每次读段开始时寻找一个碱基/接头）相比，您可能会发现一些修剪后的读段相当短（例如，如果在中间而不是在开始附近发现了接头）。因此，让我们也添加一个最小长度要求：
from Bio import SeqIO
def trim_adaptors(records, adaptor, min_len):
"""Trims perfect adaptor sequences, checks read length.
This is a generator function, the records argument should
be a list or iterator returning SeqRecord objects.
"""
len_adaptor = len(adaptor) # cache this for later
for record in records:
len_record = len(record) # cache this for later
if len(record) < min_len:
# Too short to keep
continue
index = record.seq.find(adaptor)
if index == -1:
# adaptor not found, so won't trim
yield record
elif len_record - index - len_adaptor >= min_len:
# after trimming this will still be long enough
yield record[index + len_adaptor :]
original_reads = SeqIO.parse("SRR020192.fastq", "fastq")
trimmed_reads = trim_adaptors(original_reads, "GATGACGGTGT", 100)
count = SeqIO.write(trimmed_reads, "trimmed.fastq", "fastq")
print("Saved %i reads" % count)
通过更改格式名称，您可以将其应用于FASTA文件。该代码还可以扩展为进行模糊匹配而不是精确匹配（可能使用成对对齐，或考虑读取质量分数），但那会慢得多。
转换FASTQ文件
回到部分 序列文件格式之间转换 我们展示了如何使用 Bio.SeqIO 在两种文件格式之间转换。在这里，我们将详细介绍有关第二代DNA测序中使用的FASTQ文件。请参考公鸡 et al. （2010年） [Cock2010] 以获取更长的描述。FASTQ文件存储DNA序列（作为字符串）和相关的读取质量。
PHRED分数（用于大多数FASTQ文件，也用于QUAL文件、ACE文件和SFF文件）已成为 de facto 代表测序错误概率的标准（这里表示为 \(P_e\) ）使用简单的以十为底的log转换在给定的基数：
这意味着读错 (\(P_e = 1\) ）获得PHRED质量 \(0\) ，虽然读起来很好， \(P_e = 0.00001\) 获得PHRED质量 \(50\) .虽然对于原始测序数据来说，质量高于此的情况很少见，但对于读段映射或组装等后处理，质量高达约 \(90\) 是可能的（事实上，MAQ工具允许PHRED分数在0到93范围内（含0和93）。
FASTQ格式有潜力成为 de facto 将测序读取的字母和质量分数存储在单个纯文本文件中的标准。唯一美中不足的是，FASTQ格式至少有三个版本不兼容且难以区分..
原始的Sanger FASTQ格式使用用ASC偏置33编码的PHRED质量。NCBI在其Short Read Archive中使用这种格式。我们称之为
fastq（或fastq-sanger）格式在Bio.SeqIO.Solexa（后来被Illumina收购）推出了他们自己的版本，该版本使用Solexa品质编码，以64的ASC偏差进行编码。我们称之为
fastq-solexa格式.Illumina管道1.3以后会生成具有PHRED质量（更一致）的FASTQ文件，但编码为64的ASC偏差。我们称之为
fastq-illumina格式.
Solexa质量分数使用不同的log转换定义：
鉴于Solexa/Illumina现已转向在其管道的1.3版本中使用PHRED分数，Solexa质量分数将逐渐不再使用。如果您等于误差估计 (\(P_e\) ）这两个方程允许在两个评分系统之间进行转换-Biopython包括在 Bio.SeqIO.QualityIO 模块，如果您使用，则会调用该模块 Bio.SeqIO 要将旧的Solexa/Illumina文件转换为标准的Sanger FASTQ文件：
from Bio import SeqIO
SeqIO.convert("solexa.fastq", "fastq-solexa", "standard.fastq", "fastq")
如果您想转换新的Illumina 1.3+ FASTQ文件，所有更改的只是ASC偏差，因为尽管编码不同，但分数都是PHRED质量：
from Bio import SeqIO
SeqIO.convert("illumina.fastq", "fastq-illumina", "standard.fastq", "fastq")
注意使用 Bio.SeqIO.convert() 像这样 much 比合并更快 Bio.SeqIO.parse() 和 Bio.SeqIO.write() 因为优化的代码用于FASTQ变体之间的转换（也用于FASTQ到FASTA的转换）。
对于高质量的阅读，PHRED和Solexa分数大致相等，这意味着由于两者 fasta-solexa 和 fastq-illumina 格式使用64的ASC偏差，文件几乎相同。这是Illumina深思熟虑的设计选择，这意味着应用程序期待旧的 fasta-solexa 使用较新的样式文件可能可以接受 fastq-illumina 文件（基于良好的数据）。当然，这两个变体与Sanger、NCBI和其他地方（格式名称）使用的原始FASTQ标准有很大不同 fastq 或 fastq-sanger ).
有关详细信息，请参阅内置帮助（也位于 Bio.SeqIO.QualityIO ):
>>> from Bio.SeqIO import QualityIO
>>> help(QualityIO)
将FASTA和QUAL文件转换为FASTQ文件
FASTQ文件保存 both 序列及其质量字符串。FASTA文件持有 just 序列，而QUAL文件保存 just 品质。因此，可以将单个FASTQ文件转换成或从其转换 paired FASTA和QUAL文件。
从FASTQ到FASTA很容易：
from Bio import SeqIO
SeqIO.convert("example.fastq", "fastq", "example.fasta", "fasta")
从FASTQ到QUAL也很容易：
from Bio import SeqIO
SeqIO.convert("example.fastq", "fastq", "example.qual", "qual")
然而，反之则有点棘手。您可以使用 Bio.SeqIO.parse() 在一个中重写记录 single 文件，但在本例中我们有两个输入文件。有多种可能的策略，但假设两个文件确实配对，那么最有效的内存方法就是将两者循环在一起。代码有点繁琐，所以我们提供了一个名为 PairedFastaQualIterator 在 Bio.SeqIO.QualityIO 模块来做到这一点。这需要两个手柄（FASTA文件和QUAL文件）并返回 SeqRecord 迭代器：
from Bio.SeqIO.QualityIO import PairedFastaQualIterator
for record in PairedFastaQualIterator(open("example.fasta"), open("example.qual")):
print(record)
该功能将检查FASTA和QUAL文件是否一致（例如，记录的顺序相同，并且具有相同的序列长度）。您可以将其与 Bio.SeqIO.write() 将一对FASTA和QUAL文件转换为单个FASTQ文件的函数：
from Bio import SeqIO
from Bio.SeqIO.QualityIO import PairedFastaQualIterator
with open("example.fasta") as f_handle, open("example.qual") as q_handle:
records = PairedFastaQualIterator(f_handle, q_handle)
count = SeqIO.write(records, "temp.fastq", "fastq")
print("Converted %i records" % count)
为FASTQ文件编制索引
FASTQ文件通常非常大，其中包含数百万次读取。由于数据量巨大，您无法一次将所有记录加载到内存中。这就是为什么上面的示例（过滤和修剪）会迭代查看一个文件 SeqRecord 一次。
然而，有时您不能使用大循环或迭代器-您可能需要随机访问读取。这里 Bio.SeqIO.index() 函数可能非常有帮助，因为它允许您通过其名称访问FASTQ文件中的任何读取内容（请参阅第节 将文件序列为词典-索引文件 ).
我们再次使用 SRR020192.fastq 来自ENA的文件（ftp：//ftp.sra.ebi.ac.uk/vol1/fastq/SRR 020/SRR 020192/SRR 020192.fastq.gz），尽管这实际上是一个相当小的FASTQ文件，少于 \(50,000\) 内容如下：
>>> from Bio import SeqIO
>>> fq_dict = SeqIO.index("SRR020192.fastq", "fastq")
>>> len(fq_dict)
41892
>>> list(fq_dict.keys())[:4]
['SRR020192.38240', 'SRR020192.23181', 'SRR020192.40568', 'SRR020192.23186']
>>> fq_dict["SRR020192.23186"].seq
Seq('GTCCCAGTATTCGGATTTGTCTGCCAAAACAATGAAATTGACACAGTTTACAAC...CCG')
当在具有700万次读取的FASTQ文件上测试这一点时，索引大约花了一分钟，但记录访问几乎是即时的。
姐妹功能 Bio.SeqIO.index_db() 允许您将索引保存到SQLite3数据库文件，以便几乎即时重复使用-请参阅第节 将文件序列为词典-索引文件 了解更多详细信息。
部分中的示例 排序序列文件 展示如何使用 Bio.SeqIO.index() 对大型FASTA文件进行排序的功能-这也可以用于FASTQ文件。
转换SFF文件
如果您使用454（Roche）序列数据，您可能可以访问作为标准流图格式（SFF）文件的原始数据。这包含具有质量评分和原始流程信息的序列读段（称为碱基）。
常见任务是从SFF转换为一对FASTA和QUAL文件，或转换为单个FASTQ文件。使用 Bio.SeqIO.convert() 功能（请参阅部分 序列文件格式之间转换 ):
>>> from Bio import SeqIO
>>> SeqIO.convert("E3MFGYR02_random_10_reads.sff", "sff", "reads.fasta", "fasta")
10
>>> SeqIO.convert("E3MFGYR02_random_10_reads.sff", "sff", "reads.qual", "qual")
10
>>> SeqIO.convert("E3MFGYR02_random_10_reads.sff", "sff", "reads.fastq", "fastq")
10
请记住，convert函数返回记录数，在本例中仅为十个。这会给你的 untrimmed 读取，其中引导和尾随的质量差序列或适配器将为大写字母。如果希望 trimmed 读取（使用SFF文件中记录的剪辑信息）使用以下内容：
>>> from Bio import SeqIO
>>> SeqIO.convert("E3MFGYR02_random_10_reads.sff", "sff-trim", "trimmed.fasta", "fasta")
10
>>> SeqIO.convert("E3MFGYR02_random_10_reads.sff", "sff-trim", "trimmed.qual", "qual")
10
>>> SeqIO.convert("E3MFGYR02_random_10_reads.sff", "sff-trim", "trimmed.fastq", "fastq")
10
如果您运行Linux，您可以向罗氏索要他们的“非仪器”工具（通常称为Newbler工具）的副本。这提供了在命令行上进行SFF到FASTA或QUAL转换的替代方法（但目前不支持FASTQ输出），例如
$ sffinfo -seq -notrim E3MFGYR02_random_10_reads.sff > reads.fasta
$ sffinfo -qual -notrim E3MFGYR02_random_10_reads.sff > reads.qual
$ sffinfo -seq -trim E3MFGYR02_random_10_reads.sff > trimmed.fasta
$ sffinfo -qual -trim E3MFGYR02_random_10_reads.sff > trimmed.qual
Biopython使用混合大小写序列字符串来表示修剪点的方式故意模仿了Roche工具的功能。
有关Biopython SFF支持的更多信息，请咨询内置帮助：
>>> from Bio.SeqIO import SffIO
>>> help(SffIO)
识别开放阅读框架
识别可能基因的一个非常简单化的第一步是寻找开放阅读框（ORF）。我们的意思是在所有六个框架中查看没有终止密码子的长区域--开放阅读框只是没有框内终止密码子的核酸区域。
当然，要找到基因，您还需要担心定位起始密码子、可能的启动子--在真核生物中，也有需要担心的插入子。然而，这种方法对病毒和原核生物仍然有用。
为了展示如何使用Biopython来处理这一问题，我们需要一个序列来搜索，作为示例，我们将再次使用细菌载体-尽管这次我们将从没有预先标记基因的普通FASTA文件开始： NC_005816.fna .这是一个细菌序列，因此我们想要使用NCBI密码子表11（请参阅第节 翻译 关于翻译）。
>>> from Bio import SeqIO
>>> record = SeqIO.read("NC_005816.fna", "fasta")
>>> table = 11
>>> min_pro_len = 100
这里有一个巧妙的技巧， Seq 对象的 split 方法获得六个阅读框中所有可能的ORF翻译列表：
>>> for strand, nuc in [(+1, record.seq), (-1, record.seq.reverse_complement())]:
... for frame in range(3):
... length = 3 * ((len(record) - frame) // 3) # Multiple of three
... for pro in nuc[frame : frame + length].translate(table).split("*"):
... if len(pro) >= min_pro_len:
... print(
... "%s...%s - length %i, strand %i, frame %i"
... % (pro[:30], pro[-3:], len(pro), strand, frame)
... )
...
GCLMKKSSIVATIITILSGSANAASSQLIP...YRF - length 315, strand 1, frame 0
KSGELRQTPPASSTLHLRLILQRSGVMMEL...NPE - length 285, strand 1, frame 1
GLNCSFFSICNWKFIDYINRLFQIIYLCKN...YYH - length 176, strand 1, frame 1
VKKILYIKALFLCTVIKLRRFIFSVNNMKF...DLP - length 165, strand 1, frame 1
NQIQGVICSPDSGEFMVTFETVMEIKILHK...GVA - length 355, strand 1, frame 2
RRKEHVSKKRRPQKRPRRRRFFHRLRPPDE...PTR - length 128, strand 1, frame 2
TGKQNSCQMSAIWQLRQNTATKTRQNRARI...AIK - length 100, strand 1, frame 2
QGSGYAFPHASILSGIAMSHFYFLVLHAVK...CSD - length 114, strand -1, frame 0
IYSTSEHTGEQVMRTLDEVIASRSPESQTR...FHV - length 111, strand -1, frame 0
WGKLQVIGLSMWMVLFSQRFDDWLNEQEDA...ESK - length 125, strand -1, frame 1
RGIFMSDTMVVNGSGGVPAFLFSGSTLSSY...LLK - length 361, strand -1, frame 1
WDVKTVTGVLHHPFHLTFSLCPEGATQSGR...VKR - length 111, strand -1, frame 1
LSHTVTDFTDQMAQVGLCQCVNVFLDEVTG...KAA - length 107, strand -1, frame 2
RALTGLSAPGIRSQTSCDRLRELRYVPVSL...PLQ - length 119, strand -1, frame 2
请注意，这里我们从5'结束（开始）开始计数帧， each 搁浅。有时总是从5 '结束（开始）开始计数会更容易 forward 搁浅。
您可以轻松地编辑上述基于循环的代码来建立候选蛋白质的列表，或将其转换为列表理解。现在，这个代码没有做的一件事就是跟踪蛋白质的位置。
您可以通过多种方式解决这个问题。例如，以下代码根据蛋白质计数跟踪位置，并通过乘以三转换回父序列，然后调整框架和链：
from Bio import SeqIO
record = SeqIO.read("NC_005816.gb", "genbank")
table = 11
min_pro_len = 100
def find_orfs_with_trans(seq, trans_table, min_protein_length):
answer = []
seq_len = len(seq)
for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]:
for frame in range(3):
trans = nuc[frame:].translate(trans_table)
trans_len = len(trans)
aa_start = 0
aa_end = 0
while aa_start < trans_len:
aa_end = trans.find("*", aa_start)
if aa_end == -1:
aa_end = trans_len
if aa_end - aa_start >= min_protein_length:
if strand == 1:
start = frame + aa_start * 3
end = min(seq_len, frame + aa_end * 3 + 3)
else:
start = seq_len - frame - aa_end * 3 - 3
end = seq_len - frame - aa_start * 3
answer.append((start, end, strand, trans[aa_start:aa_end]))
aa_start = aa_end + 1
answer.sort()
return answer
orf_list = find_orfs_with_trans(record.seq, table, min_pro_len)
for start, end, strand, pro in orf_list:
print(
"%s...%s - length %i, strand %i, %i:%i"
% (pro[:30], pro[-3:], len(pro), strand, start, end)
)
以及输出：
NQIQGVICSPDSGEFMVTFETVMEIKILHK...GVA - length 355, strand 1, 41:1109
WDVKTVTGVLHHPFHLTFSLCPEGATQSGR...VKR - length 111, strand -1, 491:827
KSGELRQTPPASSTLHLRLILQRSGVMMEL...NPE - length 285, strand 1, 1030:1888
RALTGLSAPGIRSQTSCDRLRELRYVPVSL...PLQ - length 119, strand -1, 2830:3190
RRKEHVSKKRRPQKRPRRRRFFHRLRPPDE...PTR - length 128, strand 1, 3470:3857
GLNCSFFSICNWKFIDYINRLFQIIYLCKN...YYH - length 176, strand 1, 4249:4780
RGIFMSDTMVVNGSGGVPAFLFSGSTLSSY...LLK - length 361, strand -1, 4814:5900
VKKILYIKALFLCTVIKLRRFIFSVNNMKF...DLP - length 165, strand 1, 5923:6421
LSHTVTDFTDQMAQVGLCQCVNVFLDEVTG...KAA - length 107, strand -1, 5974:6298
GCLMKKSSIVATIITILSGSANAASSQLIP...YRF - length 315, strand 1, 6654:7602
IYSTSEHTGEQVMRTLDEVIASRSPESQTR...FHV - length 111, strand -1, 7788:8124
WGKLQVIGLSMWMVLFSQRFDDWLNEQEDA...ESK - length 125, strand -1, 8087:8465
TGKQNSCQMSAIWQLRQNTATKTRQNRARI...AIK - length 100, strand 1, 8741:9044
QGSGYAFPHASILSGIAMSHFYFLVLHAVK...CSD - length 114, strand -1, 9264:9609
如果您注释掉排序声明，那么蛋白质序列将以与之前相同的顺序显示，因此您可以检查这是否正在做同样的事情。在这里，我们已按位置对它们进行了排序，以便更容易与ESB文件中的实际注释进行比较（如第节中所示 一个很好的例子 ).
然而，如果您只想找到开放阅读框的位置，那么翻译每个可能的密码子就是浪费时间，包括进行反向互补来搜索反向链。您需要做的就是搜索可能的停止密码子（及其反向互补）。在这里，使用正规表达是一种明显的方法（请参阅Python模块 re ).这些是描述搜索字符串的一种非常强大（但相当复杂）的方式，它受到许多编程语言以及命令行工具的支持，例如 grep 也是）。您可以找到有关这个主题的整本书！
序列解析加上简单的情节
本节使用 Bio.SeqIO module described in Chapter 序列输入/输出, plus the Python library matplotlib’s pyplot plotting interface (see the matplotlib website for a tutorial ).请注意，要遵循这些示例，您需要安装matplotlib-但如果没有它，您仍然可以尝试数据解析位。
序列长度柱状图
很多时候，您可能想要可视化数据集中序列长度的分布-例如基因组组装项目中的重叠群大小范围。在本示例中，我们将重复使用orch FASTA文件 ls_orchid.fasta 它只有94个序列。
首先，我们将使用 Bio.SeqIO 解析FASTA文件并编译所有序列长度的列表。您可以使用for循环来实现这一点，但我发现列表理解更令人满意：
>>> from Bio import SeqIO
>>> sizes = [len(rec) for rec in SeqIO.parse("ls_orchid.fasta", "fasta")]
>>> len(sizes), min(sizes), max(sizes)
(94, 572, 789)
>>> sizes
[740, 753, 748, 744, 733, 718, 730, 704, 740, 709, 700, 726, ..., 592]
既然我们已经拥有了所有基因的长度（作为一系列积分），我们就可以使用matplotlib history函数来显示它。
from Bio import SeqIO
sizes = [len(rec) for rec in SeqIO.parse("ls_orchid.fasta", "fasta")]
import matplotlib.pyplot as plt
plt.hist(sizes, bins=20)
plt.title(
"%i orchid sequences\nLengths %i to %i" % (len(sizes), min(sizes), max(sizes))
)
plt.xlabel("Sequence length (bp)")
plt.ylabel("Count")
plt.show()
图 22 兰花序列长度的柱状图。
这应该会弹出一个新窗口，其中包含中所示的图表 图 22 .请注意，大多数兰花序列都是关于 \(740\) 碱基长，这里可能有两类不同的序列，其中有一个较短序列的子集。
Tip: 而不是使用 plt.show() 要在窗口中显示情节，您还可以使用 plt.savefig(...) 将图形保存到文件（例如，作为PNG或PDF）。
序列GC%图
另一个容易计算的碱基序列量是GC%。例如，您可能想查看细菌基因组中所有基因的GC%，并调查最近可能通过水平基因转移获得的任何异常值。同样，对于这个示例，我们将重复使用orch FASTA文件 ls_orchid.fasta .
首先，我们将使用 Bio.SeqIO 解析FASTA文件并编制所有GC百分比的列表。同样，您可以使用for循环来完成这一任务，但我更喜欢这样：
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction
gc_values = sorted(
100 * gc_fraction(rec.seq) for rec in SeqIO.parse("ls_orchid.fasta", "fasta")
)
在读取每个序列并计算GC%后，我们然后将它们按递减顺序进行排序。现在我们将获取此浮点值列表并使用matplotlib绘制它们：
import matplotlib.pyplot as plt
plt.plot(gc_values)
plt.title(
"%i orchid sequences\nGC%% %0.1f to %0.1f"
% (len(gc_values), min(gc_values), max(gc_values))
)
plt.xlabel("Genes")
plt.ylabel("GC%")
plt.show()
图 23 兰花序列长度的柱状图。
与之前的示例一样，应该弹出一个新窗口，其中图形如中所示 图 23 .如果您在一种生物体的全套基因上尝试这种方法，您可能会得到比这更平滑的图。
核苷点图
点图是一种直观比较两个碱基序列彼此相似性的方法。滑动窗口用于将短子序列相互比较，通常具有不匹配阈值。为了简单起见，我们只寻找完美匹配（中以黑色显示 图 24 ).
图 24 使用图像显示的两个兰花序列的核点图。
首先，我们需要两个序列。为了方便讨论，我们只从Orange FASTA文件中取出前两个 ls_orchid.fasta :
from Bio import SeqIO
with open("ls_orchid.fasta") as in_handle:
record_iterator = SeqIO.parse(in_handle, "fasta")
rec_one = next(record_iterator)
rec_two = next(record_iterator)
我们将展示两种方法。首先，一个简单的朴素实现，将所有窗口大小的子序列相互比较以编制相似性矩阵。您可以构造矩阵或数组对象，但这里我们只使用通过嵌套列表理解创建的布尔列表列表：
window = 7
seq_one = rec_one.seq.upper()
seq_two = rec_two.seq.upper()
data = [
[
(seq_one[i : i + window] != seq_two[j : j + window])
for j in range(len(seq_one) - window)
]
for i in range(len(seq_two) - window)
]
请注意，我们有 not 在此处检查反向互补匹配。现在我们将使用matplotlib plt.imshow() 显示此数据的功能，首先请求灰色配色方案，因此以黑白方式完成：
import matplotlib.pyplot as plt
plt.gray()
plt.imshow(data)
plt.xlabel("%s (length %i bp)" % (rec_one.id, len(rec_one)))
plt.ylabel("%s (length %i bp)" % (rec_two.id, len(rec_two)))
plt.title("Dot plot using window size %i\n(allowing no mis-matches)" % window)
plt.show()
这应该会弹出一个新窗口，显示中的图表 图 24 .正如您可能预料的那样，这两个序列非常相似，沿着对角线有部分窗口大小的匹配线。不存在指示倒置或其他有趣事件的非对角线匹配。
上面的代码在小示例上运行良好，但将其应用于更大的序列会出现两个问题，我们将在下面解决这些问题。首先，这种针对所有比较的暴力方法非常缓慢。相反，我们将编译将窗口大小的子序列映射到它们的位置的字典，然后采用集合交集来查找在两个序列中找到的子序列。这会占用更多内存，但 much 快二是 plt.imshow() 该功能受其可以显示的矩阵大小限制。作为替代方案，我们将使用 plt.scatter() 功能
我们首先创建将窗口大小的子序列映射到位置的字典：
window = 7
dict_one = {}
dict_two = {}
for seq, section_dict in [
(rec_one.seq.upper(), dict_one),
(rec_two.seq.upper(), dict_two),
]:
for i in range(len(seq) - window):
section = seq[i : i + window]
try:
section_dict[section].append(i)
except KeyError:
section_dict[section] = [i]
# Now find any sub-sequences found in both sequences
matches = set(dict_one).intersection(dict_two)
print("%i unique matches" % len(matches))
为了使用 plt.scatter() 我们需要单独的列表 \(x\) 和 \(y\) 坐标：
# Create lists of x and y coordinates for scatter plot
x = []
y = []
for section in matches:
for i in dict_one[section]:
for j in dict_two[section]:
x.append(i)
y.append(j)
我们现在准备将修改后的点图绘制为散点图：
import matplotlib.pyplot as plt
plt.cla() # clear any prior graph
plt.gray()
plt.scatter(x, y)
plt.xlim(0, len(rec_one) - window)
plt.ylim(0, len(rec_two) - window)
plt.xlabel("%s (length %i bp)" % (rec_one.id, len(rec_one)))
plt.ylabel("%s (length %i bp)" % (rec_two.id, len(rec_two)))
plt.title("Dot plot using window size %i\n(allowing no mis-matches)" % window)
plt.show()
这应该会弹出一个新窗口，显示中的图表 图 25 .
图 25 使用散点图的两个兰花序列的核苷酸点图。
就我个人而言，我发现第二个情节更容易阅读！再次注意我们有 not 在这里检查反向补色匹配-您可以扩展这个示例来做到这一点，也许可以用一种颜色绘制正向匹配，用另一种颜色绘制反向匹配。
绘制测序读段数据的质量评分
如果您正在使用第二代测序数据，您可能需要尝试绘制质量数据。以下是一个使用两个包含成对末端读取的FASTQ文件的示例， SRR001666_1.fastq 对于向前读，以及 SRR001666_2.fastq 对于相反的阅读。这些内容是从ENA序列读取存档的RTP网站（ftp：//www.example.com和ftp：//ftp.sra.ebi.ac.uk/vol1/fastq/SRR 001/SRR 001666/SRR 00166_2.fastq.gz）下载的，来自ftp.sra.ebi.ac.uk/vol1/fastq/SRR001/SRR001666/SRR001666_1.fastq.gz E. coli - 详情请参阅https://www.ebi.ac.uk/ena/data/view/SRR001666。
在下面的代码中 plt.subplot(...) 使用该功能是为了并排显示两个子图上的向前和反向质量。还有一点代码仅绘制前五十次读取。
import matplotlib.pyplot as plt
from Bio import SeqIO
for subfigure in [1, 2]:
filename = "SRR001666_%i.fastq" % subfigure
plt.subplot(1, 2, subfigure)
for i, record in enumerate(SeqIO.parse(filename, "fastq")):
if i >= 50:
break # trick!
plt.plot(record.letter_annotations["phred_quality"])
plt.ylim(0, 45)
plt.ylabel("PHRED quality score")
plt.xlabel("Position")
plt.savefig("SRR001666.png")
print("Done")
您应该注意，我们正在使用 Bio.SeqIO 格式名 fastq 因为NCBI已经使用标准Sanger FASTQ格式和PHRED评分保存了这些读数。然而，正如您可能从读取长度猜测的那样，这些数据来自Illumina基因组分析仪，并且可能最初是两种Solexa/Illumina FASTQ变体文件格式之一。
此示例使用 plt.savefig(...) 函数而不是 plt.show(...) ，但如前所述，两者都是有用的。
图 26 一些配对末端读数的质量图。
中示出了结果 图 26 .
BioSQL -在关系数据库中存储序列
BioSQL 是双方的共同努力 OBF 项目（BioPerl、BioJava等）支持用于存储序列数据的共享数据库模式。从理论上讲，您可以使用BioPerl将SEN文件加载到数据库中，然后使用Biopython将其从数据库中提取为具有功能的记录对象-并获得或多或少相同的效果，就像您直接将SEN文件加载为SeqRecord一样使用 Bio.SeqIO （章 序列输入/输出 ).
Biopython的BioSQL模块目前记录在http://biopython.org/wiki/BioSQL上，该网站是我们维基页面的一部分。