最近在讀 read-mapping 的 paper 順便寫寫
這篇再補充我之前寫的 Read-mapping 的部分，還有下篇的 graph-genome 做準備
Seed Searching and Extension
先介紹一下術語，通常 read-mapping 會包含兩個階段，Seed searching跟 extension。如下圖，
- seed searching: 先找到一個 (接近)perfectly match 的一段短序列
- extension: 然後再從那個區段我們再把它往左右延伸
你會問說為什麼要兩段式，既然 seeding 跟 extension 都有 alignment 的效果，而且 extension 似乎有容錯更高，那直接用 Extension 的方式把所有 reads 貼上去不是很好嗎?
沒錯，就只是速度有差而已，所以我們才會這麼做，如下表
通常 extension 的方法，大概都是 Dynamic programming 居多，因為他是找到最準的 alignment 的方式，但是他很慢，所以有很多方式盡量少使用他。而像 Hash 這種方式是非常快的，但不接受 mismatch ，不過當 k-mer 不長 且 Read quality 很高的時候，其實還蠻好用的，所以我們稱這種乍看不精確，但其實差不了多少的東西叫 heuristic method(並沒有
我們開始以下的介紹ㄅ
Seed Searching
我粗估分成兩大主流
Hashing
這個簡單，先把 reference 跟 query 切成 k-mer (就是長度是 k 的 sequence) 然後把他弄成一個數字 (你也可以說加密字串成一個數字)。比對的時候，只要數字一樣，就代表著 k-mer 整條是依樣的，我們可以直接找到他在 reference 的哪裡。
像 BLAST(概念上)(Altschul et al., 1990), kallisto(Bray et al., 2016), Graph Genome Pipeline(Rakocevic et al., 2019，下文簡稱 GGP) 都是這樣操作
BWT 系列
Suffix Array(SA)
本名 Burrows–Wheeler transform 原本是拿來做資料壓縮的，後來發現可以拿來做 read-alignment
故事從 suffix array 說起，我們把一段序列，每個位置切斷後排序(最好的方式複雜度只要 O(N))。當我們要搜尋字串時， binary search 就可以得到了
STAR(Dobin et al., 2013) 在 paper 上是說這樣做
BWT string
然而 BWT 有個更強的特性，可以再加速。BWT string 就是 suffix array 的最後一個字，(這裡的 suffix 是頭接尾的方式，下圖的 ipssm$pissii)。
我們準備建兩個表格(壓縮起來是所謂的 FM-index)，第一個是第一個字母的位置，比如說下圖 i 在這個新的 sorted suffix array 的 index 1 到 index 4。第二個表格是，occurrence array 就是 BWT string 每個字母的出現次數。
BWT 有一個最重要的特性是… 比如說在 Index 4 的結尾 s 在 BWT string 是第二個，那麼 s 的前一個 index (Index 3) 必在第二個開頭為 s 的地方，這個特性歸功於 suffix array 是 sorted 過的。
利用這種方式，我們會倒著比對，我們把可能對到的字串的上下限做夾擠，每多一個字，夾擠的範圍會越小，詳細就直接看圖ㄅ。
如果還不懂，就看 https://web.stanford.edu/class/cs262/presentations/lecture4.pdf
附贈我的 python example code https://gist.github.com/linnil1/f6ce6001257b41c5196cf7c65e6f237b
特性二是他可以吃一點點 mismatch，預設是 3 (bwa)
工具有 bwa(Li et al., 2013), bwa-mem(Li et al., 2013), bowtie2(Langmead et al. 2012), tophat2(Kim et al., 2013), minmap2(Li et at., 2018), hisat2(Kim et al., 2019)
最後 還有加強版 bidirectional FM-index 會把 forward 跟 reverse 一起 build 起來(從 bwa-sw 開始就有)，號稱會比原本快 80% ，這裡細節請去看 paper
Extension
Before extension
在 extension 前，都還會有一個步驟，aggregating ，其實有很多種說法，比如說，clustering(GGP), seed chaining(bwa-mem), stitching(STAR)
為什麼要這個步驟呢，比如先前講的 hash 的方式，同一個 query sequences 會切成很多段去 mapping，如果每個 seeds 都去 extend 是很浪費時間的，而且有時候 seeds 是有問題的，所以會先把同一個 query string 裡面的 seeds 挑選然後合併，再做一次 extend 會比較快。
還有一個優點是你可以在合併 seeds 的時候，中間的 indel 就可以順利被抓起來。(看你合併的時候要不要有 overlap)
因為 seeds 可能會 map 到很多地方，也有人在這個步驟刪掉次佳解(bwa-mem, 留最長的 chain)，或是保留主要解(gapped-blast, 留對角線上的 seeds)，或是減少長序列的 hash set 的比對 (minimap, 保留最小的 hash 值)，任何減少浪費時間的可能
Extension
最常見的是 Dynamic Programming(DP)，獲得的 alignment 是最佳解，可是他跑很慢。
概念跟 Edit Distance 差不多，維護目前最小的 penalty score(distance) $d_{ij}$，也就是 reference string[0..i] 跟 query string [0..j] 比對的時候的能得到的最小值。你會意外的發現 substring 的最小值，會是全域的最小值
我們是允許 indel(insertion, deletion)，還順便允許 affine gap 也就是 gap open 跟 gap extend 的 penalty 不一樣
最有名的兩種
- Needleman-Wunsch(Global alignment)
- Smith-Waterman(前後 gap 不用進行計算，所以叫 local alignment)
或是參考 wiki
Bitap Algorithm
如果覺得 extension 時沒有那麼多 indels 的話，我們在 DP 前可以多做這個步驟，以減少花在 DP 的時間。 Bitap 看似暴力解，但是因為它可以用 bit operation 來操作，所以是很快的。
而且他對 mismatch 是容易處理(把 mismatch -1 的 前一個比對結果跟現在 AND 就搞定了) 而且 複雜度是跟 mismatches 數量是線性
Bitap 的概念請自己去看 wiki
BAM: A format for storing alignment data
Bam 檔可以把 alignment 的資料儲存起來的檔案，在這個年代幾乎所有 read-mapping tools 甚至是下游的分析都支援這個格式
BAM 檔有他的 format 我就不多作介紹了
官方文件 https://samtools.github.io/hts-specs/SAMv1.pdf
你可以用 samtools view xx.bam 去看裡面有什麼
然後 Visualize 的話，可以推薦 IGV http://software.broadinstitute.org/software/igv/
你可以自由選擇 genome 的區段，然後用肉眼看到哪些 read-mapping 的不好，哪些 coverage 不好，哪些看似有 variants。
接下來我們可以進到 graph-based read-mapping 的進階作法，事實上我覺得這就是未來式。
要按 claps 鼓勵作者阿!!!
Reference
- James T. Robinson, Helga Thorvaldsdóttir, Wendy Winckler, Mitchell Guttman, Eric S. Lander, Gad Getz, Jill P. Mesirov. Integrative Genomics Viewer. Nature Biotechnology 29, 24–26 (2011)
- Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., … & Gingeras, T. R. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 29(1), 15–21.
- Bray, N. L., Pimentel, H., Melsted, P., & Pachter, L. (2016). Near-optimal probabilistic RNA-seq quantification. Nature biotechnology, 34(5), 525–527.
- Rakocevic, G., Semenyuk, V., Lee, W. P., Spencer, J., Browning, J., Johnson, I. J., … & Ji, S. G. (2019). Fast and accurate genomic analyses using genome graphs. Nature genetics, 51(2), 354–362.
- Altschul, S. F., Gish, W., Miller, W., Myers, E. W., & Lipman, D. J. (1990). Basic local alignment search tool. Journal of molecular biology, 215(3), 403–410.
- Langmead, B., & Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nature methods, 9(4), 357.
- Kim, D., Paggi, J. M., Park, C., Bennett, C., & Salzberg, S. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nature biotechnology, 37(8), 907–915.
- Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., & Salzberg, S. L. (2013). TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome biology, 14(4), R36.
- Li, H., & Durbin, R. (2010). Fast and accurate long-read alignment with Burrows–Wheeler transform. Bioinformatics, 26(5), 589–595.
- Li, H. (2012). Exploring single-sample SNP and INDEL calling with whole-genome de novo assembly. Bioinformatics, 28(14), 1838–1844.
- Li, H. (2018). Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics, 34(18), 3094–3100.
- Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., … & Durbin, R. (2009). The sequence alignment/map format and SAMtools. Bioinformatics, 25(16), 2078–2079.
- Li, H., & Durbin, R. (2010). Fast and accurate long-read alignment with Burrows–Wheeler transform. Bioinformatics, 26(5), 589–595.
- Li, H. (2013). Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:1303.3997.
