一、品質篩選的目的是移除無法校正的偽序列
Illumina 和 Roche 454 定序儀皆採邊合成編定序，並偵測添加鹼基發出的螢光訊號來辨識鹼基種類。假如螢光訊號因某些原因變得模糊，儀器便無法準確讀取定序結果，導致定序品質低落。
以 Illumina 為例，隨著定序時間增加，聚合酶降解和反應速率差會導致待測序列的延長速率不同調，以致於相同定序簇 (sequencing cluster) 的鹼基發出不同螢光，使定序的正確率下降。
由於標識基因分析仰賴局部的 16S rRNA 基因以辨識物種（例如 V4 變異區，長約 250 bp），定序錯誤產生的偽序列會混淆正確序列，誇大樣本內真實存在的序列種類數。
 |
| 每個圓圈代表一種序列，圓圈大小則是該序列的讀數，圓圈間距代表序列相似度。定序錯誤產生的偽序列可能源於與之相似而且數量豐富的其他序列。 |
然而，一旦偽序列過多 ( A ) 或是錯誤太嚴重 ( B )，兩種方法皆無法校正定序錯誤導致的偽序列，只能移除低品質序列，以免影響往後的分析。
二、利用品質分數篩選偽序列
1. 品質分數代表定序結果出錯的機率
序列的品質可以用品質分數 ( quality score, Q ) 評估，不同儀器的品質分數含意不同。在 Illumina 定序平台裡，品質分數是定序儀誤判鹼基種類的機率。品質分數與錯誤率 ( error probability, P ) 的關係為 Q = -10 log( P )，值閾從 2( P = 63 %) 到 40 ( P = 0.01 %)。
Trimmomatic、Cutadapt 和 BBduk 等工具皆能依據平均品質分數篩選低品質序列，例如移除平均品質分數低於 15 (P = 3 %) 的序列，以免使用 97 % 相似度聚類 OTUs 時，誤認這些序列為新的 OTUs。
2. 平均品質分數不等同錯誤的鹼基數量
然而，平均品質分數較佳的序列可能因為鹼基錯誤率分布不均，反倒有較多誤判的鹼基（下表）。假設定序錯誤為隨機獨立事件，則長 150 bp 錯誤率為 0.0032 的序列預期出現 150 × 0.0032 = 0.48 個錯誤鹼基；而表中另一序列則預期出現 ( 140 × 0.00032 ) + ( 10 × 0.63 ) = 6.3 個錯誤鹼基。
 |
| adapted from Robert Edgar, Calculating average Q (Phred) scores is a bad idea |
3. 錯誤鹼基數量的期望值較直觀
相較於平均品質分數，直接估計誤判的鹼基數量較直觀。在品質分數可靠的情況下，可以從個別鹼基的品質分數推算一序列的錯誤鹼基數量期望值 (expected errors, E) ，其定義為所有鹼基的品質分數和。
例如一批鹼基品質分數皆為 10 ( P = 10 % ) 且長 250 bp 的序列，平均每條會有 250 × 0.1 = 25 個鹼基被定序儀器誤判。相較之下，品質分數皆為 30 ( P = 0.1 %) 者，平均只有 250 × 0.001 = 0.25 個錯誤鹼基，表示序列幾乎不會出錯。
是以，錯誤鹼基數愈多，序列品質愈差；錯誤鹼基數愈少，序列品質愈佳。選定可容忍的鹼基錯誤數以後 ( maximum number of expected errors, E_max)，移除 E > E_max 的序列即可確保剩餘序列的品質。
目前，DADA2、Vsearch、Usearch pipeline 皆能依據錯誤鹼基數量的期望值，篩選低品質的序列。
三、依品質分數篩選的限制
1. 理論計算往往低估序列的錯誤
由於使用聚類或降噪等方法能校正部分的定序錯誤，原則上 E_max 可以設定得較寬鬆。然而，計算鹼基錯誤數的前提是 (1) 定序錯誤為隨機獨立事件，(2) 品質分數忠實反映儀器誤判機率，但這些假設往往與實際情況不符，所以會設定較嚴格的標準確保資料的真確性（例如 E_max = 1，沒有錯誤）
2. 無法改善其他因素造成的誤差
此外，由於品質分數僅估計定序儀器誤判的機率，定序前引進的誤差沒辦法透過篩選低品質序列改善。以構建定序文庫的過程為例，引子 (primer) 脫靶產生的嵌合體 (chimera) 或非目標序列皆獨立於定序，所以樣本內的這些錯誤不會反映在品質分數，得藉由移除低頻率序列改善。
四、基於序列出現的頻率篩選
1. 定序錯誤產生的偽序列是 singletons
Singletons 是樣本中僅出現一次的序列，所以 singletons 都是獨特序列 (unique seqeunces) 。假設定序錯誤是隨機獨立事件，那麼便不容易在相同序列、相同位置發生相同的鹼基置換，所以定序錯誤產生的偽序列往往是 singletons（此處的前提是有足夠的定序品質）。
2. 獨特序列往往是偽序列
假設有一條長 250 bp 的序列，其鹼基品質分數皆為 30（P = 0.1 %，正確率 99.9 %），經 PCR 和定序後得 100 條序列：
- 一條序列出錯的機率 = 1 - 所有鹼基皆正確的機率 =1 - ( 0.999 ) ^ 250 = 0.22
- 錯誤序列的數量期望值 = 序列數 * 序列出錯機率 = 100 × 0.22 = 22
- 獨特序列的種類數 = 1（正確序列，n = 78） + 22（偽序列，n = 1 × 22）= 23 種
- 獨特序列的錯誤率 = 22（偽序列種類數） / 23（獨特序列種類）= 96 %
綜上所述，即使個別鹼基的判讀率達 99.9 %，仍然只有 78 % 的序列正確，而得出的獨特序列有 96 % 不可靠。
3. 多數 Singletons 是偽序列
由於 singletons 是獨特序列，而多數獨特序列屬於偽序列，所以 singletons 往往也是偽序列。在品質篩選後移除 singletons 能有效降低樣本內的偽序列。除了定序錯誤，移除 singletons 也能多少排除罕見但會影響 OTUs 或 ASVs 數量的偽序列（例如嵌合體或定序文庫汙染）。
4. 移除 singletons 利大於弊
雖然樣本裡部分 singletons 是（或近似）正確序列，移除 singletons 可能損失有用的資訊，但移除 singletons 能大幅增加特異度，僅微幅降低靈敏度。
這是因為多數的 singletons 是偽序列，而且實際比例應大於理論估算。其次，假如定序品質優良，定序錯誤產生 singletons 數量應遠低於與正確序列的數量，刪除這些 singletons 不會損害偵測到正確序列的機會。
 |
| If most bases are good, most unique sequences are bad, because good reads are all alike, but every bad read is bad in its own way. |
但如果樣本中含有大量真實的 singletons (loner tags) 呢？這可能表示沒有充分採集環境微生物的基因，得重新檢視樣本的收集與製備流程。（想像箱子裡有多顆色球。假如抽了一萬次卻只得到紅球，可以推測箱子裡的球應該沒幾種顏色；假如不斷抽到不同顏色的球，即 singletons，那應該會猜箱子裡還有別種顏色的球）
五、如何設定品質過濾的參數？
除了序列的品質和數量，雙端序列合併後的長度也能作為篩選低品質序列的依據。例如 16S rRNA 基因 V4 變異區長約 250 bp，如果合併後的序列長度偏離太多，代表引子可能接錯地方或是折疊為二級結構。剔除長度在意料外的序列可以移除這類偽序列。
品質過濾是資料量和資料品質的取捨，嚴格篩選序列會減少可用的資訊，但放寬標準又會增加樣本內的偽序列。由於定序條件差異，序列品質篩選方式應視研究目的調整，不能一味使用預設標準。
在品質篩選前，可以先檢查鹼基錯誤的位置和數量分布，預估品質篩選結果的最佳/最差狀況，再利用此資訊得到初步結果。接著依照初步結果不斷試誤，逐漸從上下限逼近最佳解。或是用模擬資料判斷品質篩選的影響。
由於雙端序列合併時能校正重疊處的鹼基品質分數，而且長度差異也會影響錯誤鹼基數的期望值，所以品質分數篩選得在雙端序列合併、長度裁減/篩選之後才能執行。因此，合理的序列品質篩選步驟為：雙端定序合併 → 長度裁剪/篩選 → 品質分數/錯誤鹼基數篩選 → 序列出現頻率篩選。
六、結論
序列品質篩選是為了移除無法校正的偽序列，依據序列的長度、錯誤鹼基數和序列出現頻率能剔除不同種類的錯誤。然而品質篩選是品質（特異度）和數量（靈敏度）的取捨，應視後續分析要求，反覆試誤以找出適當的標準。
沒有留言:
張貼留言