我們將介紹在演算法上 DADA2 為何會比其他技術還要好
背景知識
Microbiota 就是微生物菌相，就是有很多不同種不同數量的微生物聚集在一起做事，比如說人體的皮膚上，人體的腸胃道菌。然後也衍伸許多相關的名詞: microbiome, metagenomics。
重要性在於微生物在人體裡面非常多，光是腸道的微生物數量跟人體細胞總數相當，約 1:1 (paper)。微生物們有許多人類沒有的基因跟功能，因此可能會吸收人體無法代謝的物質，然後代謝出人體能用的東西，甚至跟健康(paper)息息相關，最近的糞菌移植跟腸道微生物功能就跟這個有關。
我們想要知道那群為生物是誰，有多少量，才能勉強推斷他的功能。首先我們要做定序(sequencing)，我們只看微生物都有的 16S rRNA 的其中 V3-V4 區段，也就是約 420 bp 的長度來區分這是哪個物種，這段區域剛好可以用 Illumina 2x250bp (pair-end) 來做定序 (如下圖)，最後我們要把序列拼接起來。通常能看到 Genus(屬)，Species(種) 則比較難分辨，不過沒關係，同屬的功能通常相近，且從 Genus 以上就可以探討菌相上的差異了。
難題在於我們連最一開始序列的處理都有一些困難了，更何況其完整作用機制。遇到的困難如下:
- 微生物們的序列(sequences)有些長得很像，不好區分
- 即使 Illumina 的取樣的序列 quality 有 30，但是仍然會有 0.001 的 sequencing error，i.e. 期望值就是 每 1000 base pair(bp) 有一個 error
- 有些微生物的序列根本不在 Database 裡，或是我們組錯了
DADA2 introduction
我們先說明這些困難是怎麼回事，如下圖，將 raw sequences 投影在 2D 平面(e.g. 距離可以用 mismatch 的數量當距離)，並用顏色區分物種。
因為 sequencing error 的關係，來自紅色的序列可能會靠很近但不會重疊。相反的來自綠色跟來自紅色的序列長的很像，可是事實上不太一樣。所以如果使用這個距離來 cluster 的話，就是組成 OTUs (圖右)，分不太清楚；而 DADA2 是希望能抵銷 sequencing error 來做 cluster，希望得到圖左的結果。
結論就是 DADA2 比較好，跟 OTUs-based 方法: UPARSE, MED, mothur, QIIME 這四個來比較
下圖，y 軸是序列出現的次數(頻率)，越往上代表次數越多，出現越多次，所以正確判斷比較容易。x 軸是不同物種的序列距離，越往右代表距離越遠，序列長的不像，所以越容易分開。
所以 DADA2 (圖右上)，可以分出長得像的序列(往左的藍色)。而 OTUs-based methods(圖下三個) 都有許多 false positive
DADA2 Method
核心概念
建構 sequencing error model ，使得可以分辨出誰是 error、誰是源自不同 organism
從演算法的角度出發，如果認為源自相同的 organism ，那他們應該在歸在同一個 cluster。如何歸類取決於三個因素:
Abundance
判斷某條序列該歸類在某 cluster 中，我們使用 abundance p-value 來表示其機率，p-value 越小，代表該條序列越不可能來自該 cluster。如下圖，example1 有 1000 條 AAAAA 序列，定序難免有錯， 所以 AAAAT 來自這 1000 條的機率比較高。
Quality
同樣的，如下圖 example3，如果最後一個 base 的 quality 不高，那麼 A->T 的機率也會提高，AAAAT 來自 example3 的 cluster 的機率較高。
Base
跟 quality 一樣，可能隨著 base 不同，error 的機率可能不一樣(e.g. C->T)。如下圖，統計同一個 cluster 的序列們的 base 跟 quality ，就可以做出下下圖的黑色線，這條線也代表著機器的 sequencing error model 。
Model
有了這三個元素之後，我們可以開始建構我們的 model
Transition probability
使用 quality 的資訊 跟 base 的資訊，建構出 Transition probability
額這樣很難懂，看下圖，其實就是從 序列j 變成 序列i 的機率，這個機率必須對照剛剛那張 4x4 的表
這裡要注意的是計算前有先做 alignment，如果 align 完還太多 mismatch 或是產生太多 indel 的話，那段就不會拿來算。不然 indel 會大大影響 transition probability。
Abundance p-value
額 公式也很複雜，但其實可以寫成這樣
參考 abundance 那張圖，總之我們已知 序列j 變成 序列i 的機率ㄌ，然後我們也知道 序列j 的 n(abundance，可以想成出現次數)，所以總共的 transition probability 是 n * λ，然後觀察到 序列i n 次，所以可以套用 Poisson distribution 或是 negative binomial distribution 來做。舉例來說，你假設硬幣正面的機率是 0.5，那麼觀測到連續五十個正面的但符合假設的機率有多高就可以算出來。
要注意的點是，分母的 a 不會 = 0，因為觀察不到就無法估計。
Clustering Pipeline
如果你有發現
- 好的 cluster 需要靠 好的 sequencing error model 來建構，可是
- 好的 sequencing error model 需要靠 好的 cluster 來統計
那麼整個 clustering 的 pipeline 就顯而易見ㄌ：
0. 先假設只有一群
- 從這一群裡去統計 error
- 從 error 建構 error model 並計算 序列i 是否來自 序列j
- 如果 p-value 夠小(序列j 來自 序列i) 機率很低的話，就多一個 cluster
- 重新 clustering (用剛剛的 p-value 來當作指標)
- 重複 1–4
後處理
Merge
最後 每個 cluster 取出來的代表性序列就是 de-noise 的結果，也就是去除 sequencing error 的結果，然後把這個結果拿去做接下來的處理。
Merge 就是把剛剛 overlapping 的序列黏合回去，你可以設定要不要容許 mismatch，而要保留 12 bp 以上的原因就是怕會序列太短沒有辨識度會接錯。
如果該 cluster 只有一條，那稱作 singleton，通常直接移除就好了。
Chimera
Chimera 是一個奇妙的東西，序列們會自己重組成一個看似新的序列。Merge 完後，我們在這裡會移除 bimera(就是只有兩個序列重新排列組合的序列，chimera 的一種)，其實就普通的 global alignment 的應用而已。
Taxonomy
把最後已經是乾淨的序列拿來跟資料庫做比對，你就可以知道這個序列是哪種微生物，然後拿去跟 phenotype 來觀察，不過這裡不太算是 DADA2 的內容，之後的也無關ㄌ，通通跳過。
補充 前處理
Trim
跟 trimmomatic 不同，因為這個 pair-end 要 overlap 所以長度要控制得好，所以整理序列的概念有所不同。
- 直接把序列長度 truncate 掉，因為終究他們會 overlap 所以直接剪到適合的長度就好
- Minimal quality 可以設很低，因為 de-noise 會處理很多問題
- 使用 Expected error 來進行整條序列的 quality 檢測
其中 expected error 比較少看過，其實就是我們期待整條序列的 error < 某個值。如下圖有 10bp 的序列，他的期望錯誤個數是 0.028。
Deduplicate
為了減少 cluster 的複雜度，所以長的一模一樣的先合併，合併後 abundance 當然也就跟著增加。Deduplicate 完後，就可以使用剛剛介紹的 cluster 魔法了~
How to use
至於怎麼使用，官網有 step-by-step 的介紹
Reference
- Callahan, B. J., McMurdie, P. J., Rosen, M. J., Han, A. W., Johnson, A. J. A., & Holmes, S. P. (2016). DADA2: high-resolution sample inference from Illumina amplicon data. Nature methods, 13(7), 581–583.
覺得這篇文章有收穫的話，請幫我 clap 以表示支持。
