在標識基因分析中，會以引子 (primer) 鎖定目標基因兩端的保守區，以截獲其中可提供分類資訊的變異區。例如由 Caporaso et al (2012) 設計，常用於研究人體腸道微生物的 515F/806R 引子，便能擴增 16S rRNA 的 V4 變異區，獲得長度約為 253 bp 的序列（順向擴增起點座標 - 逆向擴增起點座標 +1 - 順向引子長度 - 逆向引子長度 (806 - 515 + 1 - 19 -20 = 253 bp) ）。
僅憑印象判斷，可能會以為16S rRNA V4 區域兩端的保守區間距變化不大，所以引子擴增出的序列應具有相似長度，但實際上並非如此。下圖為 DIABIMMUNE 計畫，嬰兒腸道微生物 16S rRNA V4 區域的擴增結果。雖然多數的序列集中在 253 bp，但周圍仍有長度從 37 bp 到 258 bp 的序列。
 |
| 這張圖中呈現的資料已經事先經過許多處理，所以不是原始資料的長度分布。但我想以此圖說明，即使是從保守區擴增出來的片段，也不會只有一種長度。 |
這些長度不一的序列可能是具有生物含義的正確序列，也可能是定序和建庫過程中產生的偽序列，亦或是資料處理時引進的誤差。由於長度不一有許多成因，所以資料處理的方式也因成因而異。
此外，不同的標的基因分析流程對於序列有不同的要求，所以處理這些序列的方式也因研究議題和採用的工具而異。例如，使用 Deblur（將序列降噪為 ASVs 的工具）和 UPARSE（將序列聚類為 OTUs 的工具）等要求輸入序列長度一致的演算法時，便得在資料前處理階段統一裁切樣本中的序列；但使用 DADA2（另一種將序列降噪為 ASVs 的工具）則可省略此步。
然而，無論採取何種方法，由於多數序列的長度集中在特定範圍，序列長度在 16S rRNA 分析中可以作為判斷序列可信度的依據，透過移除極端長度者能排除偽序列或汙染序列的影響。
一、哪些原因會導致擴增序列的長度不同？
1. 變異區的長度變異
即使 16S rRNA 因為組成組成穩定而被選為分析物種親緣關係的標的基因，仍有可能因為插入、刪除或重複等遺傳變異，導致變異區出現長度變異，以至於擴增序列的長度不一。由於這些長度不一的序列屬於微生物的基因體，所以得設法在資料處理中保留。
 |
| 此圖是以 515F/806R 引子在 SILVA 資料庫模擬 PCR 擴增 V4 區域的結果，多數微生物的 V4 區域為 253 bp，但長度分布很廣泛。由此得知，雖然在前文曾提到可以用長度篩選潛在的偽序列，但序列長度實際上不是正確序列的絕對指標。 |
2. 建庫或定序過程中產生的偽序列
如果長度不一的序列皆屬樣本微生物的基因體，那麼長度便不是資料處理時最大的問題。然而，建庫和定序過程中有許多機會引進偽序列，例如引子因脫靶而擴增到非目標區域，或是順逆向引子配對成二聚體，又或在擴增階段因故中斷而形成的嵌合體等，這些偽序列的長度往往與目標基因的長度不同。
 |
| 建庫或定序過程中產生的偽序列 |
此外，定序過程也可能產生長度不同的偽序列。以 454 焦磷酸定序法為例，由於仰賴螢光強度來判斷鹼基添加的數量，當序列中出現連續鹼基時，儀器誤判的機率便會提升，導致插入或刪除的錯誤 (homopolymer misreads)。所以即使擴增同樣區域，也可能得出長度不一的序列。
3. 資料處理造成的長度差異
對於 illumina 旗下採用雙端定序的儀器而言，定序錯誤多源於誤判鹼基的種類，而不是誤判鹼基的數量，所以可以避免定序造成的長度差異。然而，合併雙端定序結果時，也可能因為順逆向序列重合處太短或是鹼基種類單調，導致合併方式錯誤，以至於結果的長度不符預期。
其次，由於 illumina 定序儀器的品質隨長度遞減，若要執行品質裁切，則會在裁切低品質的片段後留下長度不一的序列。
二、為什麼要處理長度不一的序列？
1. 序列長度不一是怎樣的問題？
處理定序資料的最終目的，是取得具有生物含義且能為分析工具所用的正確序列。序列的長度與數量、鹼基的定序品質和比例，皆為 FastQC 等品質查驗軟體常用的統計量，它們反映了樣本序列的不同特性，但本身並沒有正確與錯誤之分。
以 16S rRNA 分析為例，由於各變異區的長度變異不大，所以資料處理時會特別警戒極端長度的序列，以免偽序列膨脹了多樣性的估計值。然而長度變異常見於 18S/ITS rRNA，所以長度不一的序列可能不是問題。
因此，長度不一致的序列是怎樣的問題，端視解讀數據的假設與方法而定。
2. 區分具有生物意義的正確序列、偽序列和汙染序列
由於 16S rRNA 各變異區的長度變異不大，所以長度可以作為檢驗資料合理性的指標。極端長度的序列可能是建庫和定序過程引入的偽序列，或是未知來源的汙染序列。
以我目前在處理的 V4 變異區擴增定序資料為例，我發現其中一份樣本有五成序列的長度只有其他序列長的一半。比對核酸資料庫得知，這些序列很多都不屬於 16S rRNA 基因，隨後也在這些序列上發現用於其他建庫方式的 Nextera transposase sequences。
綜合這些線索，我猜測這份樣本受到其他文庫的汙染。由於長度差異很大，所以可以用明確的閾值篩選。然而16S rRNA 各變異區仍存在長度變異，所以長度並不是判斷序列真偽的最佳指標， 需要其他證據佐證才能判斷。
3. 奠定比較基礎以正確計算序列的數量
在標識基因分析中，會將序列聚類為 OTUs 或降噪為 ASVs，以校正定序時產生的誤差。而構建 OTUs/ASVs 的前提是 (1) 序列擴增自基因體的相同位置，(2) 正確序列的數量大於出錯序列的數量。如此一來，聚類或降噪演算法才有一致的比較基礎，得以從序列的數量推論其真偽。
然而，長度問題不利於計算序列數量和確認序列來源。假設有甲和乙兩條僅有一個鹼基差異的序列，此時有四種計算序列數量的策略。
ATCGATCGATCGCATCGATCGATCG
(1) 假設甲和乙源於不同物種的基因體，個別計算甲與乙的數量：雙端定序能確保源於相同變異區的序列有相同長度，所以得假設長度不一的序列源於不同的基因體。然而在單端定序中，源於相同區域的序列可能因品質裁切而得出長度不一的序列，導致此舉的假設不成立，亦即把相同的序列分到不同的 OTUs/ASVs。(2) 假設甲為乙的母序列，把乙的數量歸於甲：由於甲末端的鹼基沒有其它序列佐證，而乙末端遺失的序列是否與甲一致也未知。此舉意味著基於兩不可靠的資訊，得出參雜雜訊的序列。(3) 假設乙為甲的子序列，把甲的數量歸於乙：由於甲和乙前段的序列一致，可以互相印證。此舉意味著捨棄甲末端的資訊，以獲得較短但較可靠的序列。(4) 切除最後一個鹼基，泯除甲乙之別：策略 (1) 其實是將長度納為分類依據，所以能更精細地分類序列。然而，長度不一的序列往往會獨立出來聚類 OTUs ，導致不同 OTUs 生成標準不一致（部分序列依照長度聚類，部分序列依照鹼基差異聚類）。為了統一分類的標準，必須將序列裁至一定長度。雖然意義不同，但策略 (4) 的效果與策略 (3) 相似。
除了策略 (2) 明顯不合理以外，其它策略都各有優劣。然而，無論何種方法都沒有辦法完全解決此問題，因為這是使用局部基因體來區分物種的限制。
4. 符合聚類或降噪工具的要求
另外，因為處理序列的手段受限於現存的工具，所以除了品質控管和取得生物含意的序列以外，處理長度不一的序列也是為了滿足聚類或降噪工具的要求，例如 QIIME2 內的 Deblur 即要求輸入長度一致的序列。（當然這些要求也有各自的理論背景，需要閱讀原始文獻才能得知）
三、有什麼手段可以處理長度不一的序列？
1. 不處理
假設樣本中沒有任何偽序列，為了保留所有具生物含義的序列，可以選擇不處理長度不一的序列。然而如之前的討論，無論是建庫或定序過程都可能引入偽序列，而且有些軟體（例如 Deblur）要求輸入的序列長度一致，因此通常得處理長度不一的序列。
2. 長度裁切
由於單端定序往往產生長度不一的序列，必須經過裁切才能確保各序列源於標的基因上相同的起始點。而在雙端定序中，合併後的序列始於順向引子，終於逆向引子，所以能確保序列源於相同變異區，因此長度裁切反而是為了符合軟體的輸入格式要求或統一聚類的標準。
3. 長度篩選
長度篩選的重點是把偽序列或汙染序列從樣本中移除。以 16S V4 區域為例，由於這段區域的長度變異不大，所以極端長或極端短的序列很有可能是建庫擴成中引進的偽序列。但是貿然地把理論長度以外的序列全部刪除，會喪失部分序列並引入長度偏誤（參考變異區的長度分布圖）。因此，初次處理的時候可以先保留，再根據隨後的分析結果調整策略。
4. 資料庫比對篩選
由於最終的目標是取得變異區的擴增序列，所以可以透過註解序列，了解其是否具有生物意義來篩選偽序列。例如刪除相似度和 16S rRNA 基因低於閾值，或是無法註解到特定分類階層的序列。
四、怎麼決定處理序列的方式？
1. 合併雙端序列的時機為何者？
(1) 聚類/降噪前合併雙端序列，見「問題 2」
(2) 聚類/降噪後合併雙端序列（或僅使用單端序列，此時長度變異未必有生物含意），見「處理方式 A」
2. 使用的軟體是否要求輸入序列的長度一致？
(1) 是，見「處理方式 B」
(2) 否，見「問題 3」
3. 研究重點是什麼？
(1) 強調差異表達分析，犧牲聚類標準的一致性，見「處理方式 C」(2) 強調多樣性分析，犧牲分類精準度，見「處理方式 B」
處理方式 A：考量定序品質和合併雙端序列所需的重合長度，各將順逆股裁成一致長度（兩股長度不必一致）。合併雙端序列後依照擴增區域的預期長度分布，移除可能的嵌合體和偽序列（例：DADA2）。
處理方式 B：首先使用引子在資料庫中模擬擴增反應取得標區域的長度分布，並據此移除可能的嵌合體和偽序列，接著移除長度短於閾值的序列，裁切長度長於閾值的序列（例：Deblur）。
處理方式 C：不裁切序列，但仍依照擴增區域的預期長度分布，移除可能的嵌合體和偽序列。
五、結論
在 16S rRNA 標的基因分析中，處理長度不一的序列是為了取得具有生物含意的正確序列並確立比較的基礎。處理長度變異的手段取決於合併雙端序列的時機、聚類/降噪軟體的要求以及研究的重點。至於長度裁切或篩選的參數設定則沒有絕對標準，必須檢視研究目的並根據結果不斷調整，以下是我的策略：
(1) 取得標的基因的長度分布，這可以透過模擬擴增或是直接合併雙端序列得知。
(2) 依照長度分布決定一組寬鬆的長度篩選標準
(3) 參照核酸資料庫，判斷長度異常的序列是否具有生物含意，再據此調整篩選標準
(4) 參考稀釋曲線，檢查長度裁切是否太過寬鬆，誇大了樣本多樣性，導致曲線無法飽和（即樣本代表性不足）
(5) 參考多樣性分析結果，檢查長度裁切是否太過嚴格/寬鬆，以至於結果差異甚微/變異過大。
七、參考資料
Comparison of DADA2 with Deblur. Qiime2 forum.
Dada and Deblur, cutoffs and comparison. Qiime2 forum.
Defining unique sequence abundances. Usearch
Yet another question regarding OTU clustering. Qiime2 forum.
Read preparation: length trimming. usearch
To cluster or not to cluster? Qiime2 forum.
沒有留言:
張貼留言