BWA-MEM算法記錄

• BWA-MEM算法總體流程以下：

1 讀入 bwt、options、reads；

2 利用mem_chain生成chain；

3 利用mem_chain_flt過濾掉部分chain；

4 利用mem_chain2aln生成比對結果元數據。

BWA採用seed-and-extend策略。在seed階段，BWA取read的鹼基片斷在reference上進行精確匹配，並選擇知足必定匹配次數和長度要求read片斷做爲seed，這個階段算法的核心是基於FM-index的精確匹配；在extend階段，BWA利用Smith-Waterman算法將seed在read和reference上向兩邊延伸比對（容忍gap），進而找到整個read在reference上符合條件的全局匹配。

                                                                                                                                            

• BWT變換

BWT算法的主要流程是：首先將給定的字符串每次右移一個字符，獲得的全部字符串構成一個字符串矩陣，而後對該字符串矩陣中的全部字符串執行字典排序，取該字符串矩陣的最後一列即爲BWT變換字符串。下圖爲對字符串T=「acaacg」進行BWT變換（a）、BWT逆變換（b）、查找字符串（c）的過程。

• FM-index數據結構

FM-index是基於BWT的後綴查找算法，時間複雜度爲O（m），m爲待查找字符串長度。基於FM-index的序列比對算法內存佔用較小。

下圖爲字符串T（「acaacg」）構建SA（Suffix Array，後綴數組）和BWT變換的過程。

FM-index主要包含SA數組、OCC數組和BWT索引串三個輔助數據結構。SA數組保存了BW矩陣中後綴字符串在原始串中出現的位置座標。OCC[c,r]標誌BWT（T）中第r行以前字符c的出現次數。

sp和ep分別表明待匹配序列在SA數組中出現的起始和結束位置（BW矩陣的行號）。

查找的過程就是將待查找字符串從後向前逆序擴展，利用LF映射的性質進行匹配，延伸至待查找字符串的第一個字符時結束，此時區間值[sp，ep]就是在BW矩陣上精確匹配的區間。區間寬度ep-sp+1表明待查詢序列在T中的匹配個數。

LF映射有兩個特色：（1）同一行中 F的字母必定是L的後一個字母；（2）相同字母在排序後相對位置不變（例如：F列的第5個a必定是L列的第5個a）

• FM-index存儲優化

當reference比較長時，FM-index的內存佔用是至關可觀的，僅僅OCC數組可能達到幾個GB。BWA主要思想是：BWA只存儲SA後綴數組和OCC數組的一部分，其餘未存儲部分則在須要的時候即時計算出來。

BWA的OCC[a,i]只存儲i爲128的整數倍元素，知足這個條件的i稱爲一個checkpoint。若是i不是檢查點時，首先計算出離i最近的檢查點j，而後從OCC數組取出OCC[a,j],最後加上A從B[j]到B[i]出現的次數，結果就是a在BWT變換的位置i以前出現的次數。

SA[i]的空間壓縮：當i等於32的整數倍時，後綴數組的元素S[i]會被記錄下來，當待查找的序列結束位置不被SA存儲時，繼續作sp和ep的查找直至找到存儲在SA數組中的某一元素，此時用這個元素的值減去該查找步驟進行的步數，便可還原查找串的真實位置。

                                                                                                                                                            

• 基於SMEM的seed生成

MEM是一段不能繼續正向或逆向擴展的精確匹配，SMEM是一種MEM，不被其餘MEM所包含。BWA-MEM算法要求對於查詢序列的任一個位置，覆蓋該位置的最長精確匹配必須是一個SMEM。使用FMD-index和與之相關的正向-逆向擴展算法，能很快地找到全部的SMEM，SMEM查找算法以下圖：

BWA-MEM使用傳統的seed-and-extend範式，算法中以SMEM做爲seed來進行後續的比對。seed是read序列的一部分，用於在比對的初始階段衡量匹配程序，有的時候，SMEM並不出如今最終的比對結果中。爲了減小因爲SMEM選擇不當致使的誤比對，有時候須要從新生成種子（re-seeding）。好比，假設咱們找到一個長度l的SMEM，它在reference中出現了k次。當l太大時，須要從新生成種子，BWA-MEM會找出覆蓋原SMEM中間位置鹼基的最長精確匹配做爲種子，而且要求該種子在reference中至少出現k+1次，新種子的生成算法只須要對SMEM的生成算法進行少量修改。

生成的seed會被後續過程作鏈化處理。一組距離相近且共線的seed被稱爲一條鏈chain，全部的seed鏈會按照必定的規則作一次篩選。若是一條鏈的長度小於某一個閾值，或者該鏈被其餘更長的鏈所包含，則這條鏈會被拋棄。（該條鏈的seed也被捨棄？）

算法在seed階段對於FM-index的訪問和計算開銷超過總耗時的50%，主要瓶頸是對於FM-index數據結構中OCC矩陣和後綴數組SA的訪問。訪問OCC數組是爲了計算待查找字符串在BW矩陣的匹配區間，訪問SA後綴數組則爲了定位精確匹配的子序列（seed）在目標序列的起始位置。

• seed擴展

BWA-MEM使用帶狀動態規劃算法對全部seed進行擴展。在對seed進行擴展前，依據鏈的長度和seed的長度對全部的seed進行排序。若是一條鏈長度較大，則屬於這條鏈的全部seed獲得的排序就比較靠前，基於鏈長的排序結束後，再按照seed長度從長到端對鏈中全部的seed作一次排序。每個seed按照排序的前後進行擴展，首先判斷它是否被以前的比對結果所包含，若是沒有被包含，則使用帶狀動態規劃算法對它進行擴展。BWA-MEM使用的帶狀動態規劃算法經過對Smith-Waterman算法進行SSE2指令加速。

Banded DP的擴展過程當中，有可能會遇到這樣的一段序列：它的中間部分比對結果不好，而它的兩側部分的匹配質量很高。對於這樣的區域，BWA-MEM使用了Z-dropoff的啓發式思想，當比對的分數從最大值降低到低於某個閾值時，將會中止序列比對過程。

BWA-MEM算法會記錄比對到query末端的最大得分，當這個分數和局部對比的分數的差小於某一個閾值時，即便局部比對的分數較高，BWA-MEM也會將它丟棄而採用比對到末端的得分。BWA-MEM的總體算法流程以下圖。（？？？）

• OCC數組

OCC數組是按照必定的間隔存儲的，在BWA-MEM算法中，OCC數組每隔128bp存儲一個元素。在經典的FM-index結構中，OCC數組和BWT數組是分別存儲的，而BWA-MEM的FMD-index中，這兩個數組被緊湊地合併到了一個數組中，OCC和BWT組成的新數組的基本單元（一個cell）是由OCC數組和reference的BWT變換數組拼成的64個字節，其中前32個字節分別存儲在BWT數組的位置i*128以前ATCG四種鹼基出現的總次數，即四種鹼基對應的OCC數組元素，後32個字節用來存儲BWT數組在位置i*128以後128bp的鹼基。