2020年1月7日 星期二

序列品質篩選 (Read quality filtering)

Robert Edgar 是 UPARSE 和 UNOISE 等演算法的發明人,也是 Usearch 軟體的開發者,他的網站介紹了標識基因分析的重要概念。本文統整 Edgar 討論序列品質篩選的文章,介紹品質篩選方法和原則。

一、品質篩選的目的是移除無法校正的偽序列


Illumina 和 Roche 454 定序儀皆採邊合成編定序,並偵測添加鹼基發出的螢光訊號來辨識鹼基種類。假如螢光訊號因某些原因變得模糊,儀器便無法準確讀取定序結果,導致定序品質低落。

以 Illumina 為例,隨著定序時間增加,聚合酶降解和反應速率差會導致待測序列的延長速率不同調,以致於相同定序簇 (sequencing cluster) 的鹼基發出不同螢光,使定序的正確率下降。

由於標識基因分析仰賴局部的 16S rRNA 基因以辨識物種(例如 V4 變異區,長約 250 bp),定序錯誤產生的偽序列會混淆正確序列,誇大樣本內真實存在的序列種類數。

每個圓圈代表一種序列,圓圈大小則是該序列的讀數,圓圈間距代表序列相似度。定序錯誤產生的偽序列可能源於與之相似而且數量豐富的其他序列。
假如定序的品質很好,可以預期偽序列的數量稀少而且其與真實序列僅有少量差異。此時,採聚類 (clustering) 或是降噪 (denoising) 等方法,能將偽序列合併到或復原為真實序列。

然而,一旦偽序列過多 ( 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,如果合併後的序列長度偏離太多,代表引子可能接錯地方或是折疊為二級結構。剔除長度在意料外的序列可以移除這類偽序列。

品質過濾是資料量和資料品質的取捨,嚴格篩選序列會減少可用的資訊,但放寬標準又會增加樣本內的偽序列。由於定序條件差異,序列品質篩選方式應視研究目的調整,不能一味使用預設標準。

在品質篩選前,可以先檢查鹼基錯誤的位置和數量分布,預估品質篩選結果的最佳/最差狀況,再利用此資訊得到初步結果。接著依照初步結果不斷試誤,逐漸從上下限逼近最佳解。或是用模擬資料判斷品質篩選的影響。

由於雙端序列合併時能校正重疊處的鹼基品質分數,而且長度差異也會影響錯誤鹼基數的期望值,所以品質分數篩選得在雙端序列合併、長度裁減/篩選之後才能執行。因此,合理的序列品質篩選步驟為:雙端定序合併 → 長度裁剪/篩選 → 品質分數/錯誤鹼基數篩選 → 序列出現頻率篩選。


六、結論


序列品質篩選是為了移除無法校正的偽序列,依據序列的長度、錯誤鹼基數和序列出現頻率能剔除不同種類的錯誤。然而品質篩選是品質(特異度)和數量(靈敏度)的取捨,應視後續分析要求,反覆試誤以找出適當的標準。


沒有留言:

張貼留言