Bio.SeqIO.QualityIO模块

Bio.SeqIO支持FASTQ和QUAL文件格式。

请注意,您需要通过Bio.SeqIO接口使用此代码,如下所示。

Wellcome Trust Sanger Institute经常使用FASTQ文件格式来捆绑FASTA序列及其Phred质量数据(介于0和90之间的整数)。与使用单个FASTQ文件不同,通常使用成对的FASTA和QUAL文件,分别包含序列和质量信息。

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质量值a将从0到40。在QUAL格式中,这些质量值以空格分隔的文本形式保存在类似FASTA的文件格式中。在FASTQ格式中,每个质量值都使用单个ASCI字符使用CHR(Q+33)进行编码,这意味着零映射到字符“!”例如,80映射到“q”。对于桑格FASTQ标准,Phred分数的允许范围是0到93(包括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

好消息是,罗氏454音序仪可以输出QUAL格式的文件,并且他们明智地使用了像Sanger这样的PHREP风格的分数。将一对FASTA和QUAL文件转换为Sanger样式的FASTQ文件很容易。要从Roche 454 SFF二进制文件中提取QUAL文件,请使用带有-q或-qual参数的Roche Off仪器命令行工具“sffinfo”。您可以改为使用-s或-seq参数提取匹配的FASTA文件。

坏消息是Solexa/Illumina的做法不同-他们有自己的评分系统和自己的FASTQ格式的不兼容版本。Solexa/Illumina质量分数使用q=-10log10(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分数的简单质量文件(例如,来自罗氏454)

  • “FASTQ”是指使用Phred分数和ASCII偏移量33的桑格式FASTQ文件(例如,来自NCBI Short Read Archive和Illumina 1.8+)。这些可能会将Phred分数保持在0到93之间。

  • “FASTQ-Sanger”是“FASTQ”的别名。

  • “FASTQ-Solexa”表示旧的Solexa(以及非常早期的Illumina)样式的FASTQ文件,使用带有ASCII偏移量64的Solexa分数。这些可以将Solexa分数保持在-5到62之间。

  • “FASTQ-Illumina”表示较新的Illumina 1.3到1.7样式的FASTQ文件,使用Phred分数,但具有ASCII偏移量64,允许Phred分数从0到62。

我们可以潜在地添加对“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约定(Phred样式质量,ASCII OFFET为33)。这意味着我们可以使用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 Format方法以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变体(“FASTQ”或“FASTQ-Sanger”用于使用Phred值的Sanger标准,“FASTQ-Solexa”用于原始Solexa/Illumina变体,或“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_Annotation字典中(每个基数一个):

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

将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,这给出了1.25作为Phred质量,或者-4.77作为Solexa质量。因此(舍入后),随机核苷酸读取的Phred质量为1,或Solexa质量为-5。

从字面上看,这个对数公式会将Phred质量为零映射到Solexa质量为负无穷大。当然,从字面上看,Phred分数为0意味着出错的概率为1(即基本调用肯定是错误的),这比随机调用更糟糕!在实践中,Phred质量为零通常意味着默认值,或者可能是随机的-因此将其映射到Solexa的最小分数-5是合理的。

总而言之,我们遵循EMBOSS,并采用这个对数公式,但也为Solexa质量应用了最小值-5.0,并且还将Phred质量从0映射到-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)

将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

注意,小于-5的SOLEXA_QUALITY不是预期的,它将触发警告,但仍将根据对数映射进行转换(返回一个介于0和1.19之间的数字)。

作为None用于“Missing Value”的特殊情况,将返回None:

>>> print(phred_quality_from_solexa(None))
None
Bio.SeqIO.QualityIO.FastqGeneralIterator(source)

作为字符串元组(而不是SeqRecord对象)迭代FASTQ记录。

参数:
  • source-以文本模式打开的输入流,或文件路径

此代码不会尝试用数字解释质量字符串。它只是以字符串的形式返回标题、序列和质量的元组。对于顺序和质量,所有空格(如新行)都将被删除。

我们基于SeqRecord的FASTQ迭代器在内部调用此函数,然后将字符串转换为SeqRecord对象,将质量字符串映射为数字分数列表。如果要进行自定义质量映射,则可以考虑直接调用此函数。

对于解析FASTQ文件,可以选择在“+”行中省略来自每个记录开头的“@”行的标题字符串。如果重复,则必须相同。

序列字符串和质量字符串可以任选地拆分在多行上,尽管几个来源不鼓励这样做。相比之下,对于FASTA文件格式,通常使用60到80个字符之间的换行符。

WARNING -因为“@”字符可能出现在质量字符串中,这可能会导致问题,因为这也是新序列开始的标记。事实上,“+”号也可以出现。一些来源建议在质量中不使用换行符来避免这种情况,但即使这样还不够,请考虑以下示例:

@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, alphabet=None, title2ids=None)

基类:SequenceIterator

FASTQ文件的解析器。

__init__(source, alphabet=None, title2ids=None)

将FASTQ记录作为SeqRecord对象迭代。

参数:
  • source-以文本模式打开的输入流,或文件路径

  • 字母表-可选的字母表,不再使用。保留为无。

  • title2ids-一个函数,当给定FASTQ文件的标题行(没有开头>)时,它将以字符串元组的形式返回记录的id、名称和描述(按该顺序)。如果未指定,则整个标题行将用作描述,第一个单词用作id和名称。

请注意,title 2id的使用与Bio.SeqIO.FastaIO的使用相匹配。

对于(Sanger Style)FASTQ文件中的每个序列,都有一个使用偏移量33的ASCII值编码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”)存在使用偏移量33的ASCII值编码Phred品质的匹配字符串(例如“;;3;7;88”)。

直接使用此模块,您可以运行:

>>> 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]
parse(handle)

开始解析文件,并返回SeqRecord生成器。

iterate(handle)

解析文件并生成SeqRecord对象。

__abstractmethods__ = frozenset({})
Bio.SeqIO.QualityIO.FastqSolexaIterator(source, alphabet=None, title2ids=None)

解析旧的类似Solexa/Illumina FASTQ的文件(它们在质量映射上有所不同)。

可选参数与FastqPhredIterator的参数相同。

对于Solexa/Illumina FASTQ文件中的每个序列,都有一个使用偏移量为64的ASCII值编码Solexa整数质量的匹配字符串。Solexa分数的缩放比例与Phred分数不同,并且Biopython在加载时不会执行任何自动转换。

注意-此文件格式由旧版本的Solexa/Illumina管道使用。另请参见新版本的FastqIlllightaIterator函数。

例如,考虑一个包含以下五条记录的文件:

@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

让我们看一下更糟糕的伪造示例Read,其中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, alphabet=None, title2ids=None)

解析Illumina 1.3到1.7类似FASTQ的文件(它们在质量映射上有所不同)。

可选参数与FastqPhredIterator的参数相同。

对于Illumina 1.3+FASTQ文件中的每个序列,都有一个使用偏移量为64的ASCII值编码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管道编码的Solexa分数的ASCII偏移量为64。它们大致相等,但仅限于高质量读取。如果您的旧Solexa/Illumina文件的Solexa分数为负值,并尝试将其作为Illumina 1.3+文件读取,则会失败:

>>> record2 = SeqIO.read("Quality/solexa_faked.fastq", "fastq-illumina")
Traceback (most recent call last):
   ...
ValueError: Invalid character in quality string

注意-True Sanger Style FASTQ文件使用偏移量为33的Phred分数。

class Bio.SeqIO.QualityIO.QualPhredIterator(source, alphabet=None, title2ids=None)

基类:SequenceIterator

具有Phred质量分数但没有序列的QUAL文件的解析器。

__init__(source, alphabet=None, title2ids=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]

从Biopython1.59开始,此解析器将接受质量分数为负值的文件,但会将其替换为Phred分数可能最低的0。这将触发警告,以前它引发了ValueError异常。

parse(handle)

开始解析文件,并返回SeqRecord生成器。

iterate(handle)

解析文件并生成SeqRecord对象。

__abstractmethods__ = frozenset({})
class Bio.SeqIO.QualityIO.FastqPhredWriter(target, mode='w')

基类:SequenceWriter

类编写标准FASTQ格式文件(使用Phred质量分数)(已过时)。

虽然您可以直接使用此类,但强烈建议您使用 as_fastq 函数或顶层 Bio.SeqIO.write() 而是通过格式名称“FASTQ”或别名“FASTQ-Sanger”来实现。

例如,此代码读入标准桑格样式的FASTQ文件(使用Phred分数),然后将其重新保存为另一个桑格样式的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的输出文件将使每个序列在一行上,每个质量字符串在一行上(这被认为是获得最大兼容性所必需的)。

在下一个示例中,使用Phred Quality将旧样式的Solexa/Illumina FASTQ文件(使用Solexa质量分数)转换为标准的桑格样式的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 Quality 93,其编码为ASCII126(波浪号)。如果您的质量分数被截断以适应,则会发出警告。

附注:为避免弄乱您的工作目录,您现在可以删除此临时文件:

>>> import os
>>> os.remove("Quality/temp.fastq")
write_record(record)

将单个FASTQ记录写入文件。

Bio.SeqIO.QualityIO.as_fastq(record)

将SeqRecord转换为Sanger FASTQ格式的字符串。

它由SeqRecord的.format(“FASTQ”)方法和SeqIO.write(.,.,“FASTQ”)函数在内部使用,也在格式别名“FASTQ-Sanger”下使用。

class Bio.SeqIO.QualityIO.QualPhredWriter(handle, wrap=60, record2title=None)

基类: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”)方法,也会调用此代码。

附注:如果您不再需要临时文件,请不要忘记清理它:

>>> import os
>>> os.remove("Quality/temp.qual")
__init__(handle, wrap=60, record2title=None)

创建一个QUAL编写器。

参数:
  • 句柄-输出文件的句柄,例如由open(filename,“w”)返回的句柄

  • 换行-用于换行序列行的可选行长。默认情况下,序列以60个字符换行。使用零(或无)表示不换行,为序列提供单个长行。

  • record2title-可选函数,用于返回用于每条记录标题行的文本。默认情况下,使用record.id和record.description的组合。如果record.description以record.id开头,则只使用record.description。

存在record2title参数是为了与Bio.SeqIO.FastaIO编写器类保持一致。

write_record(record)

将单个QUAL记录写入文件。

Bio.SeqIO.QualityIO.as_qual(record)

将SeqRecord转换为QUAL格式的字符串。

它由SeqRecord的.format(“qual”)方法和SeqIO.write(.,.,“qual”)函数在内部使用。

class Bio.SeqIO.QualityIO.FastqSolexaWriter(target, mode='w')

基类:SequenceWriter

编写旧式Solexa/Illumina FASTQ格式文件(具有Solexa品质)(过时)。

这将输出类似于早期Solexa/Illumina管道中的FASTQ文件,使用Solexa分数和ASCII偏移量64。这些文件与标准桑格样式Phred FASTQ文件不兼容。

如果您的记录在Letter_ANNOTIONS下包含“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 Quality 62,其编码为ASCII126(波浪号)。如果您的质量分数必须截断以适应,则会发出警告。

附注:如果您不再需要临时文件,请不要忘记将其删除:

>>> import os
>>> os.remove("Quality/temp.fastq")
write_record(record)

将单个FASTQ记录写入文件。

Bio.SeqIO.QualityIO.as_fastq_solexa(record)

将SeqRecord转换为Solexa FASTQ格式的字符串。

它由SeqRecord的.format(“FASTQ-Solexa”)方法和SeqIO.write(.,.,“FASTQ-Solexa”)函数在内部使用。

class Bio.SeqIO.QualityIO.FastqIlluminaWriter(target, mode='w')

基类:SequenceWriter

写入Illumina 1.3+FASTQ格式文件(带有Phred质量分数)(已过时)。

这将使用Phred分数和ASCII偏移量64输出类似于Solexa/Illumina 1.3+管道的FASTQ文件。注意:这些文件与使用ASCII偏移量32的标准桑格样式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 Quality 62的上限,其编码为ASCII126(波浪号)。如果您的质量分数被截断以适应,则会发出警告。

write_record(record)

将单个FASTQ记录写入文件。

Bio.SeqIO.QualityIO.as_fastq_illumina(record)

将SeqRecord转换为Illumina FASTQ格式的字符串。

它由SeqRecord的.format(“FASTQ-Illumina”)方法和SeqIO.write(.,.,“FASTQ-Illumina”)函数在内部使用。

Bio.SeqIO.QualityIO.PairedFastaQualIterator(fasta_source, qual_source, alphabet=None, title2ids=None)

作为SeqRecord对象迭代匹配的FASTA和QUAL文件。

例如,请考虑以下简短的QUAL文件,其中包含Phred质量分数:

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