Bio.SeqIO.QualityIO 模組
Bio.SeqIO 支援 FASTQ 和 QUAL 檔案格式。
請注意,您應透過 Bio.SeqIO 介面使用此程式碼,如下所示。
FASTQ 檔案格式經常在 Wellcome Trust Sanger Institute 中使用,以捆綁 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) 以單個 ASCI 字元編碼,這意味著零對應於字元「!」,例如 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」表示使用 PHRED 分數和 33 的 ASCII 偏移量的 Sanger 風格 FASTQ 檔案(例如,來自 NCBI 短讀取封存和 Illumina 1.8+)。 這些檔案可能包含從 0 到 93 的 PHRED 分數。
「fastq-sanger」是「fastq」的別名。
「fastq-solexa」表示舊的 Solexa(以及非常早期的 Illumina)風格 FASTQ 檔案,使用 Solexa 分數,ASCII 偏移量為 64。 這些檔案可以包含 -5 到 62 的 Solexa 分數。
「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 慣例(具有 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 變體(Sanger 標準使用 PHRED 值時使用「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 的偏移量產生 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 品質是誤差機率的以 10 為底的對數的負十倍
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 分數表示誤差機率為 1(也就是說,鹼基呼叫絕對錯誤),這比隨機還糟! 實際上,零的 PHRED 品質通常表示預設值,或可能是隨機值 - 因此將其對應到 -5 的最低 Solexa 分數是合理的。
總之,我們遵循 EMBOSS,並採用此對數公式,但也為 Solexa 品質套用 -5.0 的最小值,並且也會將零的 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
請注意,小於 -5 的 solexa_quality 是不預期的,會觸發警告,但仍會按照對數映射進行轉換(回傳一個介於 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]]
以字串元組(而非 SeqRecord 物件)的形式迭代 Fastq 記錄。
- 參數
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(++(%
這是四個 PHRED 編碼的 FASTQ 條目,最初來自 NCBI 來源(給定讀取長度為 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)編碼 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]
若要修改解析器傳回的記錄,您可以使用產生器函式。例如,若要將平均 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
如果您想查看品質,它們會以簡單的整數列表形式記錄在每個紀錄的 per-letter-annotation 字典中
>>> 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/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
注意 - 真正的 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 檔案本身不包含序列字串。
品質分數本身可以整數列表形式在每個紀錄的 per-letter-annotation 中取得
>>> 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 分數 0。這將觸發警告,先前會引發 ValueError 例外。
- __abstractmethods__ = frozenset({})
- __parameters__ = ()
- class Bio.SeqIO.QualityIO.FastqPhredWriter(target: IO | PathLike | str | bytes, mode: str = 'w')
Bases:
SequenceWriter
用於寫入標準 FASTQ 格式檔案(使用 PHRED 品質分數)的類別(已過時)。
儘管您可以直接使用此類別,但強烈建議您改用
as_fastq
函式,或透過格式名稱 “fastq” 或別名 “fastq-sanger” 使用頂層的Bio.SeqIO.write()
函式。例如,這段程式碼會讀取一個標準 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,即波浪符號 (~)。如果您的品質分數被截斷以符合此上限,則會發出警告。
附註:為了避免您的工作目錄雜亂,您可以立即刪除這個暫存檔案。
>>> 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”) 方法,也會呼叫這段程式碼。
附註:如果您不再需要暫存檔案,請不要忘記清理它。
>>> 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 分數和 64 的 ASCII 偏移量。這些與標準 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,即波浪符號 (~)。如果您的品質分數必須截斷以符合此上限,則會發出警告。
附註:如果您不再需要暫存檔案,請不要忘記刪除它。
>>> 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]
以 SeqRecord 物件的形式迭代配對的 FASTA 和 QUAL 檔案。
例如,考慮以下包含 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 解析器一樣,如果您想查看品質,它們會以簡單的整數列表形式存在於每個記錄的 per-letter-annotation 字典中。
>>> 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")