怎麼決定合併雙端序列時允許的錯配數量？

illumina 的 Miseq 和 Hiseq 平台常用於定序 16S rRNA 基因特定變異區的增幅產物。因為這兩平台皆屬雙端定序 (paired-end sequencing)，所以每一條 DNA 都會獲得順向和逆向的序列（即定序結果）。在資料處理時，須比對這對序列的重疊區域，再基於重疊區域合併雙方，才能重建出代表該 DNA 的完整序列。

不管使用的軟體為 FLASH 、PEAR 還是 Usearch，合併雙端序列時都有許多參數可以設定。其中一個參數是合併時允許的錯配數量。這參數的必要性在於，即使順向和逆向序列都源於同一條 DNA，兩者的重疊區域仍會因為定序錯誤而略有差異。合併時允許少量的錯配數量能避免捨棄掉這些正確的序列。

然而，我們要怎麼知道洽當的錯配數量上限？使用預設值、反覆試誤、調查文獻都是常見的方法，不過我更想要知道能依照資料選定適用參數的策略。畢竟，如果能先推算出理論值，那麼測試時便有參考的基準，除了能縮減試誤時間，在寫論文或向他人解釋方法時也比較有依據。

因此，雖然在資料處理流程已成熟的今日，糾結參數設定可能被視為無關宏旨，我還是想以本文紀錄，自己思考這個問題的方向與看法。

首先，順逆向序列的重疊區只要有任一列發生定序錯誤便會發生錯配。由於 illumina 定序平台的逆向序列品質往往較順向序列品質低，假設不計兩序列同時出錯的狀況，重疊區的錯配數量可由逆向序列的預期鹼基錯誤數量估計。

由於 illumina 定序平台的逆向序列品質往往較順向序列低，而且重疊區任何一條序列發生定序錯誤便會造成錯配。假設不計兩序列於相同位點同時出錯的狀況，那麼重疊區的錯配數量上限，可以由逆向序列的預期鹼基錯誤數量估計：錯配數量 < 逆向序列重疊區的預期鹼基錯誤數 * 2

舉例來說，若增幅和定序 16S rRNA 基因 V4 區域之後，獲得一對長度為 250 bp 的序列。由於 V4 區域大約只有 253 bp 長，所以這對序列的重疊區長達 247 bp，幾乎覆蓋了整個 V4 區域。

假設考慮最極端的狀況，即要在順向和逆向序列的每個鹼基都完全對應時，我們才肯相信合併後的序列源於同一條 DNA。那麼，在逆向序列平均品質分數為 20 的狀況下，重疊區的預期錯誤數量為 $0.99 * 247 = 2.47$。

這意味著，即使是這對序列得自於同一條 DNA，平均而言也會出現 $2.47 * 2 = 4.94$ 個錯配。因此至少要允許 5 個錯配，才能確保這些源於同一條 DNA 的正確序列能順利合併。

推廣來說，錯配數量上限取決於重疊區的長度與品質。隨著對重疊區長度的要求放寬，預期的鹼基錯誤數會跟著下降，允許的錯配數量也得調低，才能確保合併雙端序列的特異性和靈敏度平衡。

依照這個推論，在取得資料時，可以先思考「若這對序列沒有定序錯誤，要有多長的重疊序列，才肯相信它們得自於同一條DNA」。

接著，找到平均品質最差的那個樣本的逆向序列，計算其 3’ 端重疊區域的預期鹼基錯誤數（該段區域的定序錯誤率之和），得出數值的兩倍，應該可以當作測試雙端序列合併的參數的起始值。