Bio.SeqIO.QualityIO 模块
Bio.SeqIO 对 FASTQ 和 QUAL 文件格式的支持。
注意,您应该通过 Bio.SeqIO 接口使用此代码,如下所示。
FASTQ 文件格式在 Wellcome Trust Sanger 研究所经常用于捆绑 FASTA 序列及其 PHRED 质量数据(0 到 90 之间的整数)。通常使用成对的 FASTA 和 QUAL 文件分别包含序列和质量信息,而不是使用单个 FASTQ 文件。
PHRED 软件读取 DNA 测序轨迹文件,调用碱基,并使用对错误概率进行对数转换来为每个调用的碱基分配一个非负质量值,Q = -10 log10( Pe ),例如
Pe = 1.0, Q = 0
Pe = 0.1, Q = 10
Pe = 0.01, Q = 20
...
Pe = 0.00000001, Q = 80
Pe = 0.000000001, Q = 90
在典型的原始测序读取中,PHRED 质量值将在 0 到 40 之间。在 QUAL 格式中,这些质量值以空格分隔的文本形式保存在类似 FASTA 的文件格式中。在 FASTQ 格式中,每个质量值使用 chr(Q+33) 用单个 ASCII 字符编码,这意味着 0 映射到字符“!”,例如 80 映射到“q”。对于 Sanger FASTQ 标准,允许的 PHRED 分数范围为 0 到 93(包括)。然后,序列和质量以类似 FASTA 的格式成对存储。
不幸的是,没有官方文档描述 FASTQ 文件格式,更糟糕的是,存在几种相关但不同的变体。有关更多详细信息,请阅读此开放获取出版物
The Sanger FASTQ file format for sequences with quality scores, and the
Solexa/Illumina FASTQ variants.
P.J.A.Cock (Biopython), C.J.Fields (BioPerl), N.Goto (BioRuby),
M.L.Heuer (BioJava) and P.M. Rice (EMBOSS).
Nucleic Acids Research 2010 38(6):1767-1771
https://doi.org/10.1093/nar/gkp1137
好消息是 Roche 454 测序仪可以输出 QUAL 格式的文件,并且明智地它们使用与 Sanger 相同的 PHREP 风格分数。将一对 FASTA 和 QUAL 文件转换为 Sanger 风格的 FASTQ 文件很容易。要从 Roche 454 SFF 二进制文件中提取 QUAL 文件,请使用 Roche 仪器外命令行工具“sffinfo”,并使用 -q 或 -qual 参数。您可以使用 -s 或 -seq 参数提取匹配的 FASTA 文件。
坏消息是 Solexa/Illumina 做的事情不同 - 他们有自己的评分系统,并且还有自己的不兼容的 FASTQ 格式版本。Solexa/Illumina 质量分数使用 Q = - 10 log10 ( Pe / (1-Pe) ),它可以为负。PHRED 分数和 Solexa 分数不可互换(但可以在它们之间实现合理的映射,并且对于高质量读取它们大致相等)。
令人困惑的是,早期的 Solexa 管道生成类似于 FASTQ 的文件,但使用他们自己的分数映射和 64 的 ASCII 偏移量。更糟糕的是,对于 Solexa/Illumina 管道 1.3 及更高版本,他们引入了 FASTQ 文件格式的第三种变体,这次使用 PHRED 分数(这更一致),但 ASCII 偏移量为 64。
即,至少存在三种不同的、不兼容的 FASTQ 文件格式变体:原始的 Sanger PHRED 标准,以及 Solexa/Illumina 的两种变体。
好消息是,从 CASAVA 版本 1.8 开始,Illumina 测序仪将使用标准的 Sanger 编码生成 FASTQ 文件。
您应该通过 Bio.SeqIO 函数使用此模块,并使用以下格式名称
“qual”表示使用 PHRED 分数的简单质量文件(例如,来自 Roche 454)
“fastq”表示 Sanger 风格的 FASTQ 文件,使用 PHRED 分数和 33 的 ASCII 偏移量(例如,来自 NCBI 短读档案和 Illumina 1.8+)。这些文件可能包含 0 到 93 的 PHRED 分数。
“fastq-sanger”是“fastq”的别名。
“fastq-solexa”表示旧的 Solexa(以及非常早期的 Illumina)风格的 FASTQ 文件,使用 Solexa 分数和 64 的 ASCII 偏移量。这些文件可以包含 -5 到 62 的 Solexa 分数。
“fastq-illumina”表示较新的 Illumina 1.3 到 1.7 风格的 FASTQ 文件,使用 PHRED 分数,但 ASCII 偏移量为 64,允许 0 到 62 的 PHRED 分数。
我们可能会添加对“qual-solexa”的支持,这意味着包含 Solexa 分数的 QUAL 文件,但到目前为止,还没有任何理由使用此类文件。
例如,请考虑以下简短的 FASTQ 文件
@EAS54_6_R1_2_1_413_324
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;88
@EAS54_6_R1_2_1_540_792
TTGGCAGGCCAAGGCCGATGGATCA
+
;;;;;;;;;;;7;;;;;-;;;3;83
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGG
+
;;;;;;;;;;;9;7;;.7;393333
它包含三个长度为 25 的读取。从读取长度来看,这些读取可能最初来自早期的 Solexa/Illumina 测序仪,但此文件遵循 Sanger FASTQ 约定(使用 33 的 ASCII 偏移量的 PHRED 风格质量)。这意味着我们可以使用 Bio.SeqIO 使用“fastq”作为格式名称来解析此文件
>>> from Bio import SeqIO
>>> for record in SeqIO.parse("Quality/example.fastq", "fastq"):
... print("%s %s" % (record.id, record.seq))
EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC
EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA
EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG
质量作为整数列表存储在每个记录的注释中
>>> print(record)
ID: EAS54_6_R1_2_1_443_348
Name: EAS54_6_R1_2_1_443_348
Description: EAS54_6_R1_2_1_443_348
Number of features: 0
Per letter annotation for: phred_quality
Seq('GTTGCTTCTGGCGTGGGTGGGGGGG')
>>> print(record.letter_annotations["phred_quality"])
[26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]
您可以使用 SeqRecord 格式方法以 QUAL 格式显示它
>>> print(record.format("qual"))
>EAS54_6_R1_2_1_443_348
26 26 26 26 26 26 26 26 26 26 26 24 26 22 26 26 13 22 26 18
24 18 18 18 18
或返回 FASTQ 格式,使用“fastq”(或“fastq-sanger”)
>>> print(record.format("fastq"))
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGG
+
;;;;;;;;;;;9;7;;.7;393333
或者,使用 Illumina 1.3+ FASTQ 编码(带有 64 的 ASCII 偏移量的 PHRED 值)
>>> print(record.format("fastq-illumina"))
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGG
+
ZZZZZZZZZZZXZVZZMVZRXRRRR
您也可以让 Biopython 转换分数并显示 Solexa 风格的 FASTQ 文件
>>> print(record.format("fastq-solexa"))
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGG
+
ZZZZZZZZZZZXZVZZMVZRXRRRR
请注意,这实际上与上面使用“fastq-illumina”作为格式的输出相同!造成这种情况的原因是所有这些分数都足够高,以至于 PHRED 和 Solexa 分数几乎相等。对于低质量读取,差异会变得明显。有关更多详细信息,请参阅函数 solexa_quality_from_phred 和 phred_quality_from_solexa。
如果您想修剪您的序列(也许是为了去除低质量区域,或去除引物序列),请尝试对 SeqRecord 对象进行切片。例如
>>> sub_rec = record[5:15]
>>> print(sub_rec)
ID: EAS54_6_R1_2_1_443_348
Name: EAS54_6_R1_2_1_443_348
Description: EAS54_6_R1_2_1_443_348
Number of features: 0
Per letter annotation for: phred_quality
Seq('TTCTGGCGTG')
>>> print(sub_rec.letter_annotations["phred_quality"])
[26, 26, 26, 26, 26, 26, 24, 26, 22, 26]
>>> print(sub_rec.format("fastq"))
@EAS54_6_R1_2_1_443_348
TTCTGGCGTG
+
;;;;;;9;7;
如果您愿意,您可以读取此 FASTQ 文件,并将其另存为 QUAL 文件
>>> from Bio import SeqIO
>>> record_iterator = SeqIO.parse("Quality/example.fastq", "fastq")
>>> with open("Quality/temp.qual", "w") as out_handle:
... SeqIO.write(record_iterator, out_handle, "qual")
3
当然,您可以读取 QUAL 文件,例如我们刚刚创建的文件
>>> from Bio import SeqIO
>>> for record in SeqIO.parse("Quality/temp.qual", "qual"):
... print("%s read of length %d" % (record.id, len(record.seq)))
EAS54_6_R1_2_1_413_324 read of length 25
EAS54_6_R1_2_1_540_792 read of length 25
EAS54_6_R1_2_1_443_348 read of length 25
请注意,QUAL 文件没有真正的序列!但质量信息在那里
>>> print(record)
ID: EAS54_6_R1_2_1_443_348
Name: EAS54_6_R1_2_1_443_348
Description: EAS54_6_R1_2_1_443_348
Number of features: 0
Per letter annotation for: phred_quality
Undefined sequence of length 25
>>> print(record.letter_annotations["phred_quality"])
[26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]
为了保持整洁,如果您正在自己按照此示例操作,现在可以删除此临时文件
>>> import os
>>> os.remove("Quality/temp.qual")
有时您不会拥有 FASTQ 文件,而是只拥有一对 FASTA 和 QUAL 文件。由于 Bio.SeqIO 系统旨在读取单个文件,因此您必须分别读取这两个文件,然后将数据组合起来。但是,由于这是您想要执行的非常常见的事情,因此此模块中定义了一个帮助程序迭代器可以为您执行此操作 - PairedFastaQualIterator。
或者,如果您有足够的 RAM 将所有记录一次性保存在内存中,那么简单的字典方法就可以工作
>>> from Bio import SeqIO
>>> reads = SeqIO.to_dict(SeqIO.parse("Quality/example.fasta", "fasta"))
>>> for rec in SeqIO.parse("Quality/example.qual", "qual"):
... reads[rec.id].letter_annotations["phred_quality"]=rec.letter_annotations["phred_quality"]
然后,您可以通过其键访问任何记录,并获取序列和质量分数。
>>> print(reads["EAS54_6_R1_2_1_540_792"].format("fastq"))
@EAS54_6_R1_2_1_540_792
TTGGCAGGCCAAGGCCGATGGATCA
+
;;;;;;;;;;;7;;;;;-;;;3;83
重要的是,您明确告诉 Bio.SeqIO 您正在使用哪个 FASTQ 变体(对于使用 PHRED 值的 Sanger 标准,“fastq”或“fastq-sanger”,对于原始的 Solexa/Illumina 变体,“fastq-solexa”,或对于较新的变体,“fastq-illumina”),因为这无法通过自动方式可靠地检测到。
为了说明这个问题,让我们考虑一个人工示例
>>> from Bio.Seq import Seq
>>> from Bio.SeqRecord import SeqRecord
>>> test = SeqRecord(Seq("NACGTACGTA"), id="Test", description="Made up!")
>>> print(test.format("fasta"))
>Test Made up!
NACGTACGTA
>>> print(test.format("fastq"))
Traceback (most recent call last):
...
ValueError: No suitable quality scores found in letter_annotations of SeqRecord (id=Test).
我们创建了一个示例 SeqRecord,并且可以以 FASTA 格式显示它 - 但对于 QUAL 或 FASTQ 格式,我们需要提供一些质量分数。这些分数作为整数列表(每个碱基一个)保存在 letter_annotations 字典中
>>> test.letter_annotations["phred_quality"] = [0, 1, 2, 3, 4, 5, 10, 20, 30, 40]
>>> print(test.format("qual"))
>Test Made up!
0 1 2 3 4 5 10 20 30 40
>>> print(test.format("fastq"))
@Test Made up!
NACGTACGTA
+
!"#$%&+5?I
我们可以检查此 FASTQ 编码 - 第一个 PHRED 质量为零,它映射到感叹号,而最后一个分数为 40,它映射到字母“I”
>>> ord('!') - 33
0
>>> ord('I') - 33
40
>>> [ord(letter)-33 for letter in '!"#$%&+5?I']
[0, 1, 2, 3, 4, 5, 10, 20, 30, 40]
类似地,我们可以使用带有 64 偏移量的 PHRED 分数生成 Illumina 1.3 到 1.7 风格的 FASTQ 文件
>>> print(test.format("fastq-illumina"))
@Test Made up!
NACGTACGTA
+
@ABCDEJT^h
我们也可以检查它 - 第一个 PHRED 分数为零,它映射到“@”,而最后一个分数为 40,它映射到“h”
>>> ord("@") - 64
0
>>> ord("h") - 64
40
>>> [ord(letter)-64 for letter in "@ABCDEJT^h"]
[0, 1, 2, 3, 4, 5, 10, 20, 30, 40]
请注意,对于相同的数据,标准的 Sanger FASTQ 和 Illumina 1.3 到 1.7 风格的 FASTQ 文件看起来多么不同!然后我们需要考虑较旧的 Solexa/Illumina 格式,它对 Solexa 分数进行编码,而不是 PHRED 分数。
首先,让我们看看如果我们将 PHRED 分数转换为 Solexa 分数(舍入到小数点后一位),Biopython 会说什么
>>> for q in [0, 1, 2, 3, 4, 5, 10, 20, 30, 40]:
... print("PHRED %i maps to Solexa %0.1f" % (q, solexa_quality_from_phred(q)))
PHRED 0 maps to Solexa -5.0
PHRED 1 maps to Solexa -5.0
PHRED 2 maps to Solexa -2.3
PHRED 3 maps to Solexa -0.0
PHRED 4 maps to Solexa 1.8
PHRED 5 maps to Solexa 3.3
PHRED 10 maps to Solexa 9.5
PHRED 20 maps to Solexa 20.0
PHRED 30 maps to Solexa 30.0
PHRED 40 maps to Solexa 40.0
现在,这是使用旧的 Solexa 风格的 FASTQ 文件的记录
>>> print(test.format("fastq-solexa"))
@Test Made up!
NACGTACGTA
+
;;>@BCJT^h
同样,它使用 64 的 ASCII 偏移量,因此我们可以检查 Solexa 分数
>>> [ord(letter)-64 for letter in ";;>@BCJT^h"]
[-5, -5, -2, 0, 2, 3, 10, 20, 30, 40]
这解释了为什么此 FASTQ 输出的最后几个字母与使用 Illumina 1.3 到 1.7 格式的输出匹配 - 高质量的 PHRED 分数和 Solexa 分数大致相等。
- Bio.SeqIO.QualityIO.solexa_quality_from_phred(phred_quality: float) float
将 PHRED 质量(范围 0 到大约 90)转换为 Solexa 质量。
PHRED 和 Solexa 质量分数都是对错误概率的对数变换(高分数 = 低错误概率)。此函数获取 PHRED 分数,将其转换回错误概率,然后将其重新表示为 Solexa 分数。这假设错误估计是等效的。
这究竟是如何运作的?PHRED 质量是错误概率的十进制对数的负十倍
phred_quality = -10*log(error,10)
因此,将其反转
error = 10 ** (- phred_quality / 10)
现在,Solexa 质量使用不同的对数变换
solexa_quality = -10*log(error/(1-error),10)
代入后,经过一些运算,我们得到
solexa_quality = 10*log(10**(phred_quality/10.0) - 1, 10)
但是,真实的 Solexa 文件使用 -5 的最小质量。这有充分的理由 - 随机的碱基调用将有 25% 的时间是正确的,因此错误概率为 0.75,这将 PHRED 质量设为 1.25,或 Solexa 质量设为 -4.77。因此(舍入后),随机的核苷酸读取将具有 1 的 PHRED 质量,或 -5 的 Solexa 质量。
从字面上讲,此对数公式将 PHRED 质量为零映射到负无穷的 Solexa 质量。当然,从字面上讲,PHRED 分数为零意味着错误概率为一(即,碱基调用绝对错误),这比随机还差!实际上,PHRED 质量为零通常意味着默认值,或者可能是随机值 - 因此将其映射到 -5 的最小 Solexa 分数是合理的。
总之,我们遵循 EMBOSS,并采用此对数公式,但也应用 -5.0 的最小值作为 Solexa 质量,并将 PHRED 质量为零也映射到 -5.0。
注意,此函数将返回一个浮点数,如果需要,您需要将其四舍五入到最接近的整数。例如:
>>> print("%0.2f" % round(solexa_quality_from_phred(80), 2)) 80.00 >>> print("%0.2f" % round(solexa_quality_from_phred(50), 2)) 50.00 >>> print("%0.2f" % round(solexa_quality_from_phred(20), 2)) 19.96 >>> print("%0.2f" % round(solexa_quality_from_phred(10), 2)) 9.54 >>> print("%0.2f" % round(solexa_quality_from_phred(5), 2)) 3.35 >>> print("%0.2f" % round(solexa_quality_from_phred(4), 2)) 1.80 >>> print("%0.2f" % round(solexa_quality_from_phred(3), 2)) -0.02 >>> print("%0.2f" % round(solexa_quality_from_phred(2), 2)) -2.33 >>> print("%0.2f" % round(solexa_quality_from_phred(1), 2)) -5.00 >>> print("%0.2f" % round(solexa_quality_from_phred(0), 2)) -5.00
请注意,对于高质量读数,PHRED 和 Solexa 评分在数值上是相等的。对于低质量读数,差异很重要,其中 PHRED 的最小值为零,而 Solexa 评分可以为负。
最后,作为使用 None 表示“缺失值”的特殊情况,返回 None。
>>> print(solexa_quality_from_phred(None)) None
- Bio.SeqIO.QualityIO.phred_quality_from_solexa(solexa_quality: float) float
将 Solexa 质量(可以为负)转换为 PHRED 质量。
PHRED 和 Solexa 质量评分都是错误概率的对数变换(高评分 = 低错误概率)。此函数接受 Solexa 评分,将其转换回错误概率,然后将其重新表示为 PHRED 评分。这假设错误估计是等效的。
基础公式在姐妹函数 solexa_quality_from_phred 的文档中给出,在本例中,操作是
phred_quality = 10*log(10**(solexa_quality/10.0) + 1, 10)
这将返回一个浮点数,如果需要,您需要将其四舍五入到最接近的整数。例如:
>>> print("%0.2f" % round(phred_quality_from_solexa(80), 2)) 80.00 >>> print("%0.2f" % round(phred_quality_from_solexa(20), 2)) 20.04 >>> print("%0.2f" % round(phred_quality_from_solexa(10), 2)) 10.41 >>> print("%0.2f" % round(phred_quality_from_solexa(0), 2)) 3.01 >>> print("%0.2f" % round(phred_quality_from_solexa(-5), 2)) 1.19
注意,solexa_quality 小于 -5 是不期望的,将触发警告,但仍将根据对数映射进行转换(返回 0 到 1.19 之间的数字)。
作为使用 None 表示“缺失值”的特殊情况,返回 None。
>>> print(phred_quality_from_solexa(None)) None
- Bio.SeqIO.QualityIO.FastqGeneralIterator(source: IO[str] | PathLike | str | bytes) Iterator[tuple[str, str, str]]
以字符串元组的形式迭代 Fastq 记录(而不是 SeqRecord 对象)。
- 参数
source - 以文本模式打开的输入流,或指向文件的路径
此代码不会尝试以数值方式解释质量字符串。它只返回标题、序列和质量作为字符串的元组。对于序列和质量,任何空格(例如换行符)都将被删除。
我们基于 SeqRecord 的 FASTQ 迭代器在内部调用此函数,然后将字符串转换为 SeqRecord 对象,将质量字符串映射到数值评分列表。如果您想进行自定义质量映射,那么您可能需要考虑直接调用此函数。
对于解析 FASTQ 文件,记录开头“@”行中的标题字符串可以可选地在“+”行上省略。如果重复,则必须相同。
序列字符串和质量字符串可以可选地拆分到多行,尽管一些来源不鼓励这样做。相比之下,对于 FASTA 文件格式,60 到 80 个字符之间的换行符是规范。
警告 - 由于“@”字符可以出现在质量字符串中,这会导致问题,因为这也是新序列开头的标记。实际上,“+”符号也可以出现。一些来源建议在质量中不要有换行符以避免这种情况,但这还不够,请考虑以下示例
@071113_EAS56_0053:1:1:998:236 TTTCTTGCCCCCATAGACTGAGACCTTCCCTAAATA +071113_EAS56_0053:1:1:998:236 IIIIIIIIIIIIIIIIIIIIIIIIIIIIICII+III @071113_EAS56_0053:1:1:182:712 ACCCAGCTAATTTTTGTATTTTTGTTAGAGACAGTG + @IIIIIIIIIIIIIIICDIIIII<%<6&-*).(*%+ @071113_EAS56_0053:1:1:153:10 TGTTCTGAAGGAAGGTGTGCGTGCGTGTGTGTGTGT + IIIIIIIIIIIICIIGIIIII>IAIIIE65I=II:6 @071113_EAS56_0053:1:3:990:501 TGGGAGGTTTTATGTGGA AAGCAGCAATGTACAAGA + IIIIIII.IIIIII1@44 @-7.%<&+/$/%4(++(%
这是来自 NCBI 源的四个 PHRED 编码 FASTQ 条目(鉴于读长为 36,这些可能是 Solexa Illumina 读取,其中质量已映射到 PHRED 值上)。
此示例已编辑以说明 FASTQ 格式中允许的一些令人讨厌的事情。首先,在“+”行上,大多数但并非所有(冗余)标识符都被省略了。在实际文件中,这些额外的标识符很可能全部存在或不存在。
其次,虽然前三个序列已显示为没有换行符,但最后一个已拆分到多行。在实际文件中,任何换行符很可能都是一致的。
第三,一些质量字符串行以“@”字符开头。对于第二条记录,这是不可避免的。但是,对于第四个序列,这只会发生是因为它的质量字符串跨两行。一个简单的解析器可能会错误地将任何以“@”开头的行视为新序列的开始!此代码通过跟踪序列的长度来处理这种可能的歧义,这给出了质量字符串的预期长度。
使用此棘手的示例文件作为输入,这段简短的代码演示了此解析函数将返回什么
>>> with open("Quality/tricky.fastq") as handle: ... for (title, sequence, quality) in FastqGeneralIterator(handle): ... print(title) ... print("%s %s" % (sequence, quality)) ... 071113_EAS56_0053:1:1:998:236 TTTCTTGCCCCCATAGACTGAGACCTTCCCTAAATA IIIIIIIIIIIIIIIIIIIIIIIIIIIIICII+III 071113_EAS56_0053:1:1:182:712 ACCCAGCTAATTTTTGTATTTTTGTTAGAGACAGTG @IIIIIIIIIIIIIIICDIIIII<%<6&-*).(*%+ 071113_EAS56_0053:1:1:153:10 TGTTCTGAAGGAAGGTGTGCGTGCGTGTGTGTGTGT IIIIIIIIIIIICIIGIIIII>IAIIIE65I=II:6 071113_EAS56_0053:1:3:990:501 TGGGAGGTTTTATGTGGAAAGCAGCAATGTACAAGA IIIIIII.IIIIII1@44@-7.%<&+/$/%4(++(%
最后,我们注意到一些来源指出质量字符串应该以“!”开头(使用 PHRED 映射意味着第一个字母始终具有零的质量评分)。这个相当严格的规则没有得到广泛遵守,因此在这里被忽略了。关于“!”规则的一个优点是(只要质量序列中没有换行符),它会防止上述“@”字符问题。
- class Bio.SeqIO.QualityIO.FastqPhredIterator(source: IO[str] | PathLike | str | bytes, alphabet: None = None)
碱基:
SequenceIterator
[str
]用于 FASTQ 文件的解析器。
- __init__(source: IO[str] | PathLike | str | bytes, alphabet: None = None)
以 SeqRecord 对象的形式迭代 FASTQ 记录。
- 参数
source - 以文本模式打开的输入流,或指向文件的路径
alphabet - 可选的字母表,不再使用。保留为 None。
对于(Sanger 风格)FASTQ 文件中的每个序列,都有一个匹配的字符串,使用 ASCII 值(偏移量为 33)编码 PHRED 质量(0 到约 90 之间的整数)。
例如,考虑一个包含三个短读的
@EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC + ;;3;;;;;;;;;;;;7;;;;;;;88 @EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA + ;;;;;;;;;;;7;;;;;-;;;3;83 @EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG + ;;;;;;;;;;;9;7;;.7;393333
对于每个序列(例如“CCCTTCTTGTCTTCAGCGTTTCTCC”),都有一个匹配的字符串,使用 ASCII 值(偏移量为 33)(例如“;;3;;;;;;;;;;;;7;;;;;;;88”)编码 PHRED 质量。
直接使用此模块,您可能会运行
>>> with open("Quality/example.fastq") as handle: ... for record in FastqPhredIterator(handle): ... print("%s %s" % (record.id, record.seq)) EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG
但是,通常情况下,您会通过 Bio.SeqIO 而不是“fastq”(或“fastq-sanger”)作为格式来调用它
>>> from Bio import SeqIO >>> with open("Quality/example.fastq") as handle: ... for record in SeqIO.parse(handle, "fastq"): ... print("%s %s" % (record.id, record.seq)) EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG
如果您想查看质量,它们将记录在每个记录的每字母注释字典中,作为简单整数列表
>>> print(record.letter_annotations["phred_quality"]) [26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]
要修改解析器返回的记录,可以使用生成器函数。例如,要将平均 PHRED 质量存储在记录描述中,请使用
>>> from statistics import mean >>> def modify_records(records): ... for record in records: ... record.description = mean(record.letter_annotations['phred_quality']) ... yield record ... >>> with open('Quality/example.fastq') as handle: ... for record in modify_records(FastqPhredIterator(handle)): ... print(record.id, record.description) ... EAS54_6_R1_2_1_413_324 25.28 EAS54_6_R1_2_1_540_792 24.52 EAS54_6_R1_2_1_443_348 23.4
- __abstractmethods__ = frozenset({})
- __orig_bases__ = (Bio.SeqIO.Interfaces.SequenceIterator[str],)
- __parameters__ = ()
- Bio.SeqIO.QualityIO.FastqSolexaIterator(source: IO[str] | PathLike | str | bytes, alphabet: None = None) Iterator[SeqRecord]
解析旧版Solexa/Illumina FASTQ文件(质量映射方式不同)。
可选参数与FastqPhredIterator相同。
在Solexa/Illumina FASTQ文件中,每个序列对应一个字符串,该字符串使用ASCII值(偏移量为64)编码Solexa整数质量值。Solexa评分与PHRED评分的比例不同,Biopython在加载时不会进行任何自动转换。
注意 - 此文件格式由旧版Solexa/Illumina管道使用。另请参见FastqIlluminaIterator函数,用于新版。
例如,考虑一个包含以下五个记录的文件
@SLXA-B3_649_FC8437_R1_1_1_610_79 GATGTGCAATACCTTTGTAGAGGAA +SLXA-B3_649_FC8437_R1_1_1_610_79 YYYYYYYYYYYYYYYYYYWYWYYSU @SLXA-B3_649_FC8437_R1_1_1_397_389 GGTTTGAGAAAGAGAAATGAGATAA +SLXA-B3_649_FC8437_R1_1_1_397_389 YYYYYYYYYWYYYYWWYYYWYWYWW @SLXA-B3_649_FC8437_R1_1_1_850_123 GAGGGTGTTGATCATGATGATGGCG +SLXA-B3_649_FC8437_R1_1_1_850_123 YYYYYYYYYYYYYWYYWYYSYYYSY @SLXA-B3_649_FC8437_R1_1_1_362_549 GGAAACAAAGTTTTTCTCAACATAG +SLXA-B3_649_FC8437_R1_1_1_362_549 YYYYYYYYYYYYYYYYYYWWWWYWY @SLXA-B3_649_FC8437_R1_1_1_183_714 GTATTATTTAATGGCATACACTCAA +SLXA-B3_649_FC8437_R1_1_1_183_714 YYYYYYYYYYWYYYYWYWWUWWWQQ
直接使用此模块,您可能会运行
>>> with open("Quality/solexa_example.fastq") as handle: ... for record in FastqSolexaIterator(handle): ... print("%s %s" % (record.id, record.seq)) SLXA-B3_649_FC8437_R1_1_1_610_79 GATGTGCAATACCTTTGTAGAGGAA SLXA-B3_649_FC8437_R1_1_1_397_389 GGTTTGAGAAAGAGAAATGAGATAA SLXA-B3_649_FC8437_R1_1_1_850_123 GAGGGTGTTGATCATGATGATGGCG SLXA-B3_649_FC8437_R1_1_1_362_549 GGAAACAAAGTTTTTCTCAACATAG SLXA-B3_649_FC8437_R1_1_1_183_714 GTATTATTTAATGGCATACACTCAA
通常,您会通过Bio.SeqIO使用“fastq-solexa”作为格式来调用它
>>> from Bio import SeqIO >>> with open("Quality/solexa_example.fastq") as handle: ... for record in SeqIO.parse(handle, "fastq-solexa"): ... print("%s %s" % (record.id, record.seq)) SLXA-B3_649_FC8437_R1_1_1_610_79 GATGTGCAATACCTTTGTAGAGGAA SLXA-B3_649_FC8437_R1_1_1_397_389 GGTTTGAGAAAGAGAAATGAGATAA SLXA-B3_649_FC8437_R1_1_1_850_123 GAGGGTGTTGATCATGATGATGGCG SLXA-B3_649_FC8437_R1_1_1_362_549 GGAAACAAAGTTTTTCTCAACATAG SLXA-B3_649_FC8437_R1_1_1_183_714 GTATTATTTAATGGCATACACTCAA
如果要查看质量,它们在每个记录的每字母注释字典中作为一个简单的整数列表记录
>>> print(record.letter_annotations["solexa_quality"]) [25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 23, 25, 25, 25, 25, 23, 25, 23, 23, 21, 23, 23, 23, 17, 17]
这些分数不是很好,但它们足够高,以至于几乎完全映射到PHRED分数
>>> print("%0.2f" % phred_quality_from_solexa(25)) 25.01
让我们看一下更糟糕的伪造示例读取,其中Solexa分数和PHRED分数之间存在更明显的差异
@slxa_0001_1_0001_01 ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTNNNNNN +slxa_0001_1_0001_01 hgfedcba`_^]\[ZYXWVUTSRQPONMLKJIHGFEDCBA@?>=<;
同样,您通常会使用Bio.SeqIO读取此文件(而不是直接调用Bio.SeqIO.QualtityIO模块)。大多数FASTQ文件将包含数千个读取,因此您通常会使用Bio.SeqIO.parse()如上所示。此示例只有一个条目,因此我们可以使用Bio.SeqIO.read()函数
>>> from Bio import SeqIO >>> with open("Quality/solexa_faked.fastq") as handle: ... record = SeqIO.read(handle, "fastq-solexa") >>> print("%s %s" % (record.id, record.seq)) slxa_0001_1_0001_01 ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTNNNNNN >>> print(record.letter_annotations["solexa_quality"]) [40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, -1, -2, -3, -4, -5]
这些质量分数非常低,以至于从Solexa方案转换为PHRED分数时,它们看起来非常不同
>>> print("%0.2f" % phred_quality_from_solexa(-1)) 2.54 >>> print("%0.2f" % phred_quality_from_solexa(-5)) 1.19
注意,您可以使用Bio.SeqIO.write()函数或SeqRecord的format方法输出记录
>>> print(record.format("fastq-solexa")) @slxa_0001_1_0001_01 ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTNNNNNN + hgfedcba`_^]\[ZYXWVUTSRQPONMLKJIHGFEDCBA@?>=<;
注意,此输出与输入文件略有不同,因为Biopython省略了“+”行上可选的序列标识符重复。如果您想使用PHRED分数,请使用“fastq”或“qual”作为输出格式,Biopython将为您进行转换
>>> print(record.format("fastq")) @slxa_0001_1_0001_01 ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTNNNNNN + IHGFEDCBA@?>=<;:9876543210/.-,++*)('&&%%$$##""
>>> print(record.format("qual")) >slxa_0001_1_0001_01 40 39 38 37 36 35 34 33 32 31 30 29 28 27 26 25 24 23 22 21 20 19 18 17 16 15 14 13 12 11 10 10 9 8 7 6 5 5 4 4 3 3 2 2 1 1
如上所示,质量较低的Solexa读取已映射到等效的PHRED分数(例如,-5到1,如前所述)。
- Bio.SeqIO.QualityIO.FastqIlluminaIterator(source: IO[str] | PathLike | str | bytes, alphabet: None = None) Iterator[SeqRecord]
解析Illumina 1.3到1.7 FASTQ文件(质量映射方式不同)。
可选参数与FastqPhredIterator相同。
在Illumina 1.3+ FASTQ文件中,每个序列对应一个字符串,该字符串使用ASCII值(偏移量为64)编码PHRED整数质量值。
>>> from Bio import SeqIO >>> record = SeqIO.read("Quality/illumina_faked.fastq", "fastq-illumina") >>> print("%s %s" % (record.id, record.seq)) Test ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTN >>> max(record.letter_annotations["phred_quality"]) 40 >>> min(record.letter_annotations["phred_quality"]) 0
注意 - 旧版Solexa/Illumina管道使用ASCII偏移量为64编码Solexa分数。它们大约相等,但仅对于高质量读取。如果您有一个带有负Solexa分数的旧版Solexa/Illumina文件,并尝试将其读取为Illumina 1.3+文件,则会失败
>>> record2 = SeqIO.read("Quality/solexa_faked.fastq", "fastq-illumina") Traceback (most recent call last): ... ValueError: Invalid character in quality string
注意 - 真正的Sanger风格FASTQ文件使用PHRED分数,偏移量为33。
- class Bio.SeqIO.QualityIO.QualPhredIterator(source: IO[str] | PathLike | str | bytes, alphabet: None = None)
Bases:
SequenceIterator
解析包含PHRED质量分数但没有序列的QUAL文件。
- __init__(source: IO[str] | PathLike | str | bytes, alphabet: None = None) None
用于包含PHRED质量分数但没有序列的QUAL文件。
例如,考虑这个简短的QUAL文件
>EAS54_6_R1_2_1_413_324 26 26 18 26 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 26 26 26 23 23 >EAS54_6_R1_2_1_540_792 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 26 12 26 26 26 18 26 23 18 >EAS54_6_R1_2_1_443_348 26 26 26 26 26 26 26 26 26 26 26 24 26 22 26 26 13 22 26 18 24 18 18 18 18
直接使用此模块,您可能会运行
>>> with open("Quality/example.qual") as handle: ... for record in QualPhredIterator(handle): ... print("%s read of length %d" % (record.id, len(record.seq))) EAS54_6_R1_2_1_413_324 read of length 25 EAS54_6_R1_2_1_540_792 read of length 25 EAS54_6_R1_2_1_443_348 read of length 25
通常,您会通过Bio.SeqIO使用“qual”作为格式来调用它
>>> from Bio import SeqIO >>> with open("Quality/example.qual") as handle: ... for record in SeqIO.parse(handle, "qual"): ... print("%s read of length %d" % (record.id, len(record.seq))) EAS54_6_R1_2_1_413_324 read of length 25 EAS54_6_R1_2_1_540_792 read of length 25 EAS54_6_R1_2_1_443_348 read of length 25
只有序列长度是已知的,因为QUAL文件本身不包含序列字符串。
质量分数本身在每个记录的每字母注释中作为一个整数列表提供
>>> print(record.letter_annotations["phred_quality"]) [26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]
您仍然可以对这些SeqRecord对象进行切片
>>> sub_record = record[5:10] >>> print("%s %s" % (sub_record.id, sub_record.letter_annotations["phred_quality"])) EAS54_6_R1_2_1_443_348 [26, 26, 26, 26, 26]
从Biopython 1.59开始,此解析器将接受包含负质量分数的文件,但会将其替换为最低的可能PHRED分数零。这将触发警告,以前会引发ValueError异常。
- __abstractmethods__ = frozenset({})
- __parameters__ = ()
- class Bio.SeqIO.QualityIO.FastqPhredWriter(target: IO | PathLike | str | bytes, mode: str = 'w')
Bases:
SequenceWriter
用于写入标准FASTQ格式文件(使用PHRED质量分数)(已过时)。
虽然您可以直接使用此类,但强烈建议您使用
as_fastq
函数或顶级Bio.SeqIO.write()
函数,通过格式名称“fastq”或别名“fastq-sanger”。例如,此代码读取一个标准的Sanger风格的FASTQ文件(使用PHRED分数),并将其重新保存为另一个Sanger风格的FASTQ文件
>>> from Bio import SeqIO >>> record_iterator = SeqIO.parse("Quality/example.fastq", "fastq") >>> with open("Quality/temp.fastq", "w") as out_handle: ... SeqIO.write(record_iterator, out_handle, "fastq") 3
如果您想这样做,是因为原始文件包含额外的换行符,虽然有效,但可能不被所有工具支持。Biopython的输出文件将每个序列放在一行上,每个质量字符串放在一行上(这被认为是最大兼容性的理想选择)。
在下一个示例中,旧版Solexa/Illumina FASTQ文件(使用Solexa质量分数)被转换为使用PHRED质量的标准Sanger风格的FASTQ文件
>>> from Bio import SeqIO >>> record_iterator = SeqIO.parse("Quality/solexa_example.fastq", "fastq-solexa") >>> with open("Quality/temp.fastq", "w") as out_handle: ... SeqIO.write(record_iterator, out_handle, "fastq") 5
如果使用SeqRecord的.format(“fastq”)方法,或.format(“fastq-sanger”)方法(如果您喜欢该别名),也会调用此代码。
请注意,Sanger FASTQ 文件的 PHRED 质量值上限为 93,对应 ASCII 码 126,即波浪号。如果您的质量得分被截断以适应此限制,则会发出警告。
P.S. 为了避免工作目录混乱,您现在可以删除此临时文件。
>>> import os >>> os.remove("Quality/temp.fastq")
- Bio.SeqIO.QualityIO.as_fastq(record: SeqRecord) str
将 SeqRecord 转换为 Sanger FASTQ 格式的字符串。
此函数在 SeqRecord 的 .format(“fastq”) 方法和 SeqIO.write(…, …, “fastq”) 函数内部使用,以及在格式别名 “fastq-sanger” 中使用。
- class Bio.SeqIO.QualityIO.QualPhredWriter(handle: IO[str] | PathLike | str | bytes, wrap: int = 60, record2title: Callable[[SeqRecord], str] | None = None)
Bases:
SequenceWriter
用于写入 QUAL 格式文件(使用 PHRED 质量分数)的类(已过时)。
虽然您可以直接使用此类,但强烈建议您使用
as_qual
函数或顶层Bio.SeqIO.write()
函数。例如,以下代码读取 FASTQ 文件并将质量分数保存到 QUAL 文件中
>>> from Bio import SeqIO >>> record_iterator = SeqIO.parse("Quality/example.fastq", "fastq") >>> with open("Quality/temp.qual", "w") as out_handle: ... SeqIO.write(record_iterator, out_handle, "qual") 3
如果您使用 SeqRecord 的 .format(“qual”) 方法,也会调用此代码。
P.S. 如果不再需要临时文件,请务必清理它。
>>> import os >>> os.remove("Quality/temp.qual")
- __init__(handle: IO[str] | PathLike | str | bytes, wrap: int = 60, record2title: Callable[[SeqRecord], str] | None = None) None
创建 QUAL 写入器。
- 参数
handle - 输出文件的句柄,例如由 open(filename, “w”) 返回的句柄。
wrap - 用于包装序列行的可选行长。默认情况下,将序列包装为 60 个字符。使用零(或 None)表示不包装,从而为序列提供单个长行。
record2title - 用于返回每个记录的标题行中使用的文本的可选函数。默认情况下,使用 record.id 和 record.description 的组合。如果 record.description 以 record.id 开头,则仅使用 record.description。
record2title 参数是为了与 Bio.SeqIO.FastaIO 写入器类保持一致而提供的。
- Bio.SeqIO.QualityIO.as_qual(record: SeqRecord) str
将 SeqRecord 转换为 QUAL 格式的字符串。
此函数在 SeqRecord 的 .format(“qual”) 方法和 SeqIO.write(…, …, “qual”) 函数内部使用。
- class Bio.SeqIO.QualityIO.FastqSolexaWriter(target: IO | PathLike | str | bytes, mode: str = 'w')
Bases:
SequenceWriter
写入旧式 Solexa/Illumina FASTQ 格式文件(使用 Solexa 质量值)(已过时)。
此操作输出类似于早期 Solexa/Illumina 管道中的 FASTQ 文件,使用 Solexa 分数和 ASCII 偏移量 64。这些文件与标准 Sanger 式 PHRED FASTQ 文件不兼容。
如果您的记录在 letter_annotations 中包含 “solexa_quality” 条目,则使用此条目,否则,在使用 solexa_quality_from_phred 函数转换后使用任何 “phred_quality” 条目。如果两种类型的质量分数都不存在,则会引发异常。
虽然您可以直接使用此类,但强烈建议您使用
as_fastq_solexa
函数或顶层Bio.SeqIO.write()
函数。例如,以下代码读取 FASTQ 文件并将其重新保存为另一个 FASTQ 文件>>> from Bio import SeqIO >>> record_iterator = SeqIO.parse("Quality/solexa_example.fastq", "fastq-solexa") >>> with open("Quality/temp.fastq", "w") as out_handle: ... SeqIO.write(record_iterator, out_handle, "fastq-solexa") 5
如果您想这样做,可能是因为原始文件包含额外的换行符,尽管这些换行符是有效的,但可能不是所有工具都支持。Biopython 生成的输出文件将每个序列放在单行上,每个质量字符串也放在单行上(这被认为是为了最大程度的兼容性而理想的做法)。
如果您使用 SeqRecord 的 .format(“fastq-solexa”) 方法,也会调用此代码。例如:
>>> record = SeqIO.read("Quality/sanger_faked.fastq", "fastq-sanger") >>> print(record.format("fastq-solexa")) @Test PHRED qualities from 40 to 0 inclusive ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTN + hgfedcba`_^]\[ZYXWVUTSRQPONMLKJHGFECB@>;;
请注意,Solexa FASTQ 文件的 Solexa 质量值上限为 62,对应 ASCII 码 126,即波浪号。如果您的质量得分必须被截断以适应此限制,则会发出警告。
P.S. 如果不再需要临时文件,请务必删除它。
>>> import os >>> os.remove("Quality/temp.fastq")
- Bio.SeqIO.QualityIO.as_fastq_solexa(record: SeqRecord) str
将 SeqRecord 转换为 Solexa FASTQ 格式的字符串。
此函数在 SeqRecord 的 .format(“fastq-solexa”) 方法和 SeqIO.write(…, …, “fastq-solexa”) 函数内部使用。
- class Bio.SeqIO.QualityIO.FastqIlluminaWriter(target: IO | PathLike | str | bytes, mode: str = 'w')
Bases:
SequenceWriter
写入 Illumina 1.3+ FASTQ 格式文件(使用 PHRED 质量分数)(已过时)。
此操作输出类似于 Solexa/Illumina 1.3+ 管道中的 FASTQ 文件,使用 PHRED 分数和 ASCII 偏移量 64。请注意,这些文件与使用 ASCII 偏移量 32 的标准 Sanger 式 PHRED FASTQ 文件不兼容。
虽然您可以直接使用此类,但强烈建议您使用
as_fastq_illumina
或顶层Bio.SeqIO.write()
函数,并使用格式名称 “fastq-illumina”。如果您使用 SeqRecord 的 .format(“fastq-illumina”) 方法,也会调用此代码。例如:>>> from Bio import SeqIO >>> record = SeqIO.read("Quality/sanger_faked.fastq", "fastq-sanger") >>> print(record.format("fastq-illumina")) @Test PHRED qualities from 40 to 0 inclusive ACGTACGTACGTACGTACGTACGTACGTACGTACGTACGTN + hgfedcba`_^]\[ZYXWVUTSRQPONMLKJIHGFEDCBA@
请注意,Illumina FASTQ 文件的 PHRED 质量上限为 62,其 ASCII 编码为 126,即波浪号 (~)。如果您的质量分数被截断以适合文件,则会发出警告。
- Bio.SeqIO.QualityIO.as_fastq_illumina(record: SeqRecord) str
将一个 SeqRecord 转换为 Illumina FASTQ 格式的字符串。
这在 SeqRecord 的 .format(“fastq-illumina”) 方法和 SeqIO.write(…, …, “fastq-illumina”) 函数内部使用。
- Bio.SeqIO.QualityIO.PairedFastaQualIterator(fasta_source: IO[str] | PathLike | str | bytes, qual_source: IO[str] | PathLike | str | bytes, alphabet: None = None) Iterator[SeqRecord]
将匹配的 FASTA 和 QUAL 文件迭代为 SeqRecord 对象。
例如,考虑这个包含 PHRED 质量分数的简短 QUAL 文件
>EAS54_6_R1_2_1_413_324 26 26 18 26 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 26 26 26 23 23 >EAS54_6_R1_2_1_540_792 26 26 26 26 26 26 26 26 26 26 26 22 26 26 26 26 26 12 26 26 26 18 26 23 18 >EAS54_6_R1_2_1_443_348 26 26 26 26 26 26 26 26 26 26 26 24 26 22 26 26 13 22 26 18 24 18 18 18 18
以及匹配的 FASTA 文件
>EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC >EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA >EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG
您可以使用 Bio.SeqIO 的“qual”和“fasta”格式分别解析这两个文件,但您将得到一组没有序列的 SeqRecord 对象,以及一组具有序列但没有质量的匹配对象。由于 Bio.SeqIO 只能处理一个输入文件句柄,因此无法一起读取这两个文件——但这个函数可以!例如,
>>> with open("Quality/example.fasta") as f: ... with open("Quality/example.qual") as q: ... for record in PairedFastaQualIterator(f, q): ... print("%s %s" % (record.id, record.seq)) ... EAS54_6_R1_2_1_413_324 CCCTTCTTGTCTTCAGCGTTTCTCC EAS54_6_R1_2_1_540_792 TTGGCAGGCCAAGGCCGATGGATCA EAS54_6_R1_2_1_443_348 GTTGCTTCTGGCGTGGGTGGGGGGG
与 FASTQ 或 QUAL 解析器一样,如果您想查看质量,它们将作为简单的整数列表存储在每个记录的每字母注释字典中
>>> print(record.letter_annotations["phred_quality"]) [26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 26, 22, 26, 26, 13, 22, 26, 18, 24, 18, 18, 18, 18]
如果您有以 FASTQ 格式文件存储的数据,则直接使用它会更简单快捷。请注意,您可以轻松使用此函数将配对的 FASTA 和 QUAL 文件转换为 FASTQ 文件
>>> from Bio import SeqIO >>> with open("Quality/example.fasta") as f: ... with open("Quality/example.qual") as q: ... SeqIO.write(PairedFastaQualIterator(f, q), "Quality/temp.fastq", "fastq") ... 3
并且不要忘记在不再需要临时文件时将其清理掉
>>> import os >>> os.remove("Quality/temp.fastq")