生物資訊學/bwa

軟件介紹

新的DNA測序技術產生了大量的短read，需要開發快速準確的read比對程式。MAQ是一種基於雜湊; 哈希=>zh-mo:雜湊表的方法，用於從單個個體比對短read，但對於頻繁出現插入缺失的較長read和數百個個體的比對不合適。為了解決這些限制，開發了Burrows-Wheeler Alignment工具（BWA），它可以有效地將短測序read與大型參考序列（例如人類基因組）進行比對，允許不配對（misMactch）和間隙（gap）。BWA支援基於鹼基空間的read，例如Illumina測序機器和來自AB SOLiD機器的顏色空間read。模擬和真實數據的評估表明，BWA比MAQ快約10-20倍，同時實現相似的準確性。此外，BWA輸出新的標準SAM（Sequence Alignment/Map）格式的比對。在比對之後進行變異分析和其他下游分析可以使用開源的SAMtools軟件包。

軟件網址：http://maq.sourceforge.net。

對參考基因組進行短讀比對方面已經發展了很多演算法，但是對於長讀（>200 bp）的比對，現有的演算法大多針對低測序錯誤率的短讀進行調整，不適用或者效率不高。因此，需要設計一種能夠高效比對長序列的演算法。作者基於Burrows-Wheeler演算法提出了一種新演算法，稱為Burrows-Wheeler Aligner's Smith-Waterman Alignment（BWA-SW），能夠使用少量記憶體; 内存=>zh-mo:記憶體對長達1 Mb的序列進行比對。結果顯示，該演算法比BLAT更加準確，與SSAHA2的準確率相當，並且比兩者都快幾十倍。BWA-SW的原始碼已公開。

bwa-sw軟件網站：http://bio-bwa.sourceforge.net/

軟件發表時的文章：

使用方法

本頁面主要介紹BWA的原理，關於使用請參考BWA的使用。

背景介紹

bwa軟件開發的背景

將Illumina/Solexa測序技術產生的短read序列比對到大型基因組的比對軟件，分為三個類別：基於雜湊; 哈希=>zh-mo:雜湊的比對程式、基於基因組的比對程式和基於字串合併排序的比對程式。最近，使用Burrows-Wheeler Transform（BWT）的字串配對理論引起了幾個研究小組的關注，這導致了SOAPv2(http://soap.genomics.org.cn/)、Bowtie(Langmead et al., 2009)和BWA等比對軟件的開發。BWA是一種新的比對程式，使用BWT和反向搜尋（Ferragina 和 Manzini，2000 年；Lippert，2005 年）來進行精確和近似配對，具有相對較小的記憶體; 内存=>zh-mo:記憶體佔用（Lam 等人，2008 年），能夠快速且準確地將短read序列比對到大型基因組。BWA比對軟件被證明在模擬和真實數據上的表現優於其他比對程式，特別是在記憶體; 内存=>zh-mo:記憶體使用方面。

基於雜湊; 哈希=>zh-mo:雜湊的比對程式：Eland (Cox, 2007, unpublished material)、RMAP (Smith et al., 2008)、MAQ (Li et al., 2008a)、ZOOM (Lin et al., 2008)、SeqMap (Jiang and Wong, 2008)、CloudBurst (Schatz, 2009) 和 SHRiMP (http://compbio.cs.toronto.edu/shrimp)
基於基因組的比對程式：SOAPv1 (Li et al., 2008b)、PASS (Campagna et al., 2009)、MOM (Eaves and Gao, 2009)、ProbeMatch (Jung Kim et al., 2009)、NovoAlign (http://www.novocraft.com)、ReSEQ (http://code.google.com/p/re-seq)、Mosaik (http://bioinformatics.bc.edu/marthlab/Mosaik) 和 BFAST (http://genome.ucla.edu/bfast)
基於字串合併排序：Malhis 等人，2009 年

BWA-SW演算法

隨着 1990 年左右的局部比對軟件的發展，例如 FASTA（Pearson 和 Lipman，1988）和 BLAST（Altschul 等，1997），自 2000 年以來開發了新一代更快的 DNA 序列配對方法，包括 MegaBLAST (Morgulis et al., 2008; Zhang et al., 2000)、SSAHA2 (Ning et al., 2001)、BLAT (Kent, 2002) 和 PatternHunter (Ma et al., 2002)，大大加快了配對測序read針對大型參考基因組。當新的測序技術出現並產生數百萬個短read（<100 bp）時，各種新演算法被開發出來，速度提高了 10-1000 倍，包括 SOAP (Li,R. et al., 2008)、MAQ (Li, H. et al., 2008)、 Bowtie（Langmead 等人，2009 年）和 BWA（Li 和 Durbin，2009 年）。然而，Roche/454 測序技術已經開始產生了 >400 bp 的read讀數，Illumina 逐漸增加了 >100 bp 的read讀數，而 Pacific Bioscience 在早期測試中產生了 1000 bp 的讀數（Eid 等，2009）。來自新測序技術的read不再是短的序列，這有效地排除了許多專為不超過 100 bp 的read設計的比對軟件。

隨着產生長read讀數的新測序技術的出現，現有的短read讀數比對演算法對於長read讀數比對來說並不高效。長read比對的目的是尋找局部配對，因為長read更容易受到結構變異和參考基因組序列中錯配組裝的影響，但受靠近read末端的錯配影響較小。此外，由於長read中更頻繁地出現插入和缺失，所以長讀數比對演算法必須對比對缺口gap持寬容態度。

BWA-SW是一種新的長read排列演算法，在兩個FM-indices之間使用動態編程來尋找種子，並在它們在參考序列中很少出現的情況下進行擴充。該演算法旨在減少高度重複序列的不必要的擴充，從而提高速度。

BWA原理方法

字首樹和字串配對

字串 X 的字首樹是一棵樹，其中每條邊都用一個符號標記，從葉子到根的路徑上邊符號的字串連接給出了唯一的字首 X。在字首樹上，字串連接從節點到根的邊符號給出了 X 的唯一子串，稱為節點所代表的串。請注意，X 的字首樹與 X 的反向字尾樹相同，因此字尾樹理論也可以應用於字首樹。

使用字首 trie，測試查詢 W 是否是 X 的精確子串相當於找到表示 W 的節點，這可以通過將 W 中的每個符號與一條邊配對，在 O(|W|) 時間內完成，從根。為了允許不配對，可以徹底遍歷特里樹並將 W 配對到每個可能的路徑。稍後將展示如何通過使用 W 的字首資訊來加速此搜尋。圖1給出了「GOOGOL」的字首樹的範例。每個節點中的字尾陣列（SA）間隔在第 2.3.3 節中解釋。

Burrows–Wheeler轉換

設 Σ 是一個字母表。符號 $ 不存在於 Σ 中，並且按字典順序小於 Σ 中的所有符號。字串 X=a0a1…an−1 總是以符號 $ 結尾（即 an−1=$），並且這個符號只出現在末尾。設 X[i]=ai, i=0, 1,..., n−1, 是 X 的第 i 個符號，X[i, j]=ai ...aj 是一個子串，Xi=X[i, n− 1] X 的字尾。X 的字尾陣列 S 是整數 0…n−1 的排列，使得 S(i) 是第 i 個最小字尾的起始位置。 X 的 BWT 在 S(i)=0 時定義為 B[i]=$，否則定義為 B[i]=X[S(i)-1]。還將字串 X 的長度定義為 |X| 因此 |X|=|B|=n。圖2 給出了如何構造 BWT 和字尾陣列的範例。

圖 2 所示的演算法在時間和空間上是二次的。然而，這不是必需的。在實踐中，通常先構造字尾陣列，然後生成 BWT。大多數構建字尾陣列的演算法至少需要 n⌈log2n⌉ 位的工作空間，對於人類基因組來說，這相當於 12 GB。最近，Hon 等人。 (2007) 給出了一種新演算法，該演算法使用 n 位工作空間，並且在構建人類基因組 BWT 的高峰時間只需要 <1 GB 記憶體; 内存=>zh-mo:記憶體。該演算法在 BWT-SW 中實現（Lam 等人，2008 年）。修改了原始碼以使其與 BWA 一起使用。

字尾陣列間隔和序列比對

如果字串 W 是 X 的子字串，則 X 中 W 每次出現的位置將出現在字尾陣列中的一個區間內。這是因為所有以 W 作為字首的字尾都排序在一起。基於這一觀察，定義：

${\underline {R}}(W)=min\left\{k:W\ is\ the\ prefix\ of\ X_{S(k)}\right\}$ (1)

${\overline {R}}(W)=max\left\{k:W\ is\ the\ prefix\ of\ X_{S(k)}\right\}$ (2)

特別是，如果 W 是空字串，則 R(W)=1 和 ${\textstyle {\overline {R}}=n-1}$ 。

${\textstyle [{\underline {R}}(W),{\overline {R}}(W)]}$ 稱為 W 的 SA 區間，X 中所有出現 W 的位置集合是 ${\textstyle \left\{S(k):{\underline {R}}(W)\leq k\leq {\overline {R}}(W)\right\}}$ 。例如在圖 2 中，字串『go』的 SA 區間為 [1, 2]。此區間中的字尾陣列值為 3 和 0，給出了所有「go」出現的位置。

知道字尾陣列中的間隔可以得到位置。因此，序列比對相當於搜尋X的子串的SA區間與查詢配對。對於精確配對問題，只能找到一個這樣的區間；對於不精確配對的問題，可能有很多。

精確配對：向後搜尋

令 C(a) 是 C(a)=#{0≤j≤n−2:X[j]<a}和O(a,i)=#{0≤j≤i:B[j]<a}中按字典順序小於 a ∈ Σ 的符號數，O(a, i) 是 a 在 B[0, i] 中的出現次數。 Ferragina 和 Manzini (2000) 證明了如果 W 是 X 的子串：

${\underline {R}}(aW)=C(a)+O(a,{\underline {R}}(W)-1)+1$ (3)

${\overline {R}}(aW)=C(a)+O(a,{\overline {R}}(W))$ (4)

若且唯若 aW 是 X 的子串時， ${\underline {R}}(aW)\leq {\overline {R}}(aW)$ 。此結果使得可以通過迭代計算 ${\underline {R}}$ 和 ${\overline {R}}$ 來測試 W 是否是 X 的子串並計算 W 在 O(|W|) 時間內的出現次數從W結尾。這個過程叫做反向搜尋。

需要注意的是，等式（3）和（4）實際上實現了對X的字首trie的自頂向下遍歷，因為如果知道其父節點的間隔，可以在恆定時間內計算出子節點的SA間隔 . 從這個意義上說，向後搜尋相當於字首樹上的精確字串配對，但沒有明確地將樹放入記憶體; 内存=>zh-mo:記憶體中。

不精確配對：有界遍歷/回溯

給出了一種遞歸演算法，用於搜尋 X 的子串的 SA 區間，該區間與查詢字串 W 配對且差異不超過 z（不配對或間隙）。本質上，該演算法使用向後搜尋從基因組中採樣不同的子串。此過程受 D(·) 陣列的限制，其中 D(i) 是 W[0, i] 中差異數量的下限。 D 估計得越好，搜尋空間越小，演算法效率越高。通過為所有 i 設置 D(i)=0 來實現一個樸素的界限，但由此產生的演算法在差異數量上顯然是指數級的，並且效率較低。

Precalculation:
	Calculate BWT string B for reference string X
	Calculate array C(·) and Ο(·,·) from B
	Calculate BWT string B' for the reverse reference
	Calculate array Ο'(·) from B'
Procedures:
	INEXACTSEARCH(W,z)
		CalculateD(W) 
		return InexRecur(W, |W|-1,z,1,|X|-1)
	CalculateD(W)
		k ← 1
		z ← 0
		for i = 0 to |W|-1 do
			k ← C(W[i])+O'(W[i],k-1)+1
			l ← C(W[i])+O'(W[i],l)
			if k > l then
				k ← l
				l ← |X|-1 
				z ← z+1
				D(i) ← z
	InexRecur( W, z, k, l)
		if z < D(i) then
			return Φ 
		if z <0 then 
			return {[k,l]}
		I ← Φ
		I ← I ∪ InexRecur(W, i-1,z-1,k,l) 
		for each b∈{A,C,G,T} do
			k ← C(b) + O(b,k-1)+1
			l ← C(b)+O(b,l)
            if k ≤ l then
            	I ← I ∪ InexRecur(W,i,z-1,k,l) 
            	if b = W[i] then
            		I ← I ∪ InexRecur(W,i- 1 ,z,k,l) 
            	else
            		I ← I ∪ InexRecur(W,i-1 ,z-1,k,l) 
      return I

圖3：不精確搜尋配對 W 的子串的 SA 區間的演算法。參考 X 以 $ 終止，而 W 以 A/C/G/T 終止。過程 InexactSearch(W, z) 返回與 W 配對的子串的 SA 間隔，差異不超過 z（不配對或間隙）； InexRecur(W, i, z, k, l) 在字尾Wi+1配對區間[k, l]的條件下，遞歸計算配對W[0, i]且差異不超過z的子串的SA區間。以星號開頭的行分別用於插入和刪除 X。 D(i) 是字串 W[0, i] 中差異數量的下限。圖 3 中的 CalculateD 過程提供了一個更好但不是最佳的界限。它在概念上與圖 4 中描述的相同，後者更易於理解。使用反向（未二補碼）參考序列的 BWT 來測試 W 的子串是否也是 X 的子串。請注意，單獨使用 BWT 字串 B 進行此測試將使CalculateD 成為 ${\textstyle O(|W|^{2})}$ O(|W|2) 過程，而不是圖 3 中描述的 ${\textstyle O(|W|)}$ 。

CalculateD(W) 
	z ← 0
	j ← 0
	for i=0 to |W| — 1 do
		if W[j,i] is not a substring of X then 
			z ← z+1
			j ← i+l
			D(i) ← z

圖4：計算 D(i) 的等效演算法。

為了理解D的作用，回到在X=GOOGOL$中搜尋W=LOL的例子（圖1）。如果為所有 i 設置 D(i)=0 並禁止間隙（去除演算法中的兩條星線），則 InexRecur 的呼叫圖是一棵樹，有效地模仿了圖 1 中虛線所示的搜尋路線. 然而，通過CalculateD，知道D(0)=0 和D(1)=D(2)=1。然後可以避免下降到字首樹中的「G」和「O」子樹以獲得更小的搜尋空間。

圖 3 中的演算法保證找到所有允許最大 z 差異的區間。理論上是完整的，但在實踐中，也做了各種修改。首先，對不配對、間隙開放和間隙擴充支付不同的懲罰，這對生物數據更現實。其次，使用類似堆的數據結構來保持部分命中而不是使用遞歸。堆狀結構優先於部分命中的對齊分數，以使 BWA 始終首先找到最佳間隔。同時處理反向互補讀取序列。請注意，圖 3 中描述的遞歸有效地模擬了字首樹上的深度優先搜尋 (DFS)，而 BWA 使用這種類似堆的數據結構實現了廣度優先搜尋 (BFS)。第三，採用迭代策略：如果頂部區間是重複的，預設不搜尋次優區間；如果頂部間隔是唯一的並且具有 z 差異，則只搜尋最多具有 z + 1 差異的命中。這種迭代策略在保持生成對映質素的能力的同時加速了 BWA。然而，這也使得 BWA 的速度對讀取和參考之間的不配對率敏感，因為找到具有更多差異的命中通常更慢。第四，允許對讀取的前幾十個鹼基對的最大允許差異設置限制，稱之為種子序列。給定 70 bp 模擬讀數，與 32 bp 種子中最大兩個差異的比對比沒有種子快 2.5 倍。對齊錯誤率，即模擬中置信對映中錯誤對齊的比例（另見第 2.3.3 節），僅從 0.08% 增加到 0.11%。對於較短的讀取，播種效果較差。

降低記憶體; 内存=>zh-mo:記憶體

上述演算法需要在記憶體; 内存=>zh-mo:記憶體中載入出現陣列O和字尾陣列S。儲存; 保存=>zh-mo:儲存完整的 O 和 S 陣列需要巨大的記憶體; 内存=>zh-mo:記憶體。幸運的是，可以通過只儲存 O 和 S 陣列的一小部分並動態計算其餘部分來減少記憶體; 内存=>zh-mo:記憶體。 BWT-SW（Lam 等人，2008 年）和 Bowtie（Langmead 等人，2009 年）使用了由 Ferragina 和 Manzini（2000 年）首次引入的類似策略。

給定大小為 n 的基因組，出現陣列 O(·, ·) 需要 $4n[log_{2}n]$ 位，因為每個整數需要 $[log_{2}n]$ 位，並且陣列中有 4n 個。在實踐中，在記憶體; 内存=>zh-mo:記憶體中儲存 O(·, k) 的 k 為 128 的因子，並使用 BWT 字串 B 計算其餘元素。當使用兩位表示核苷酸時，B 需要 2n 位。因此，用於反向搜尋的記憶體; 内存=>zh-mo:記憶體為 2n+n⌈log2n⌉/32 $2n+n[log_{2}n]/32$ 位。由於還需要儲存反向基因組的 BWT 來計算邊界，因此計算間隔所需的記憶體; 内存=>zh-mo:記憶體增加了一倍，或者對於 3 Gb 基因組大約需要 2.3 GB。

列舉每個出現的位置需要字尾陣列 S。如果將整個 S 放在記憶體; 内存=>zh-mo:記憶體中，它將使用 n⌈log2n⌉ 位。然而，當知道其中的一部分時，也可以重建整個 S。事實上，S 和逆壓縮字尾陣列 (inverse CSA) $\psi ^{-1}$ (Grossi and Vitter, 2000) 滿足：

$S(k)=S((\psi ^{-1})^{(j)}(k))+j$ (5)

其中 $(\psi ^{-1})^{(j)}$ 表示將變換 $\psi ^{-1}$ 重複應用 j 次。逆 CSA $\psi ^{-1}$ 可以用出現陣列 O 計算：

${\psi ^{-1}}(i)=C(B[i]+O(B[i],i)$ (6)

在 BWA 中，只在記憶體; 内存=>zh-mo:記憶體 S(k) 中儲存可以被 32 整除的 k。對於不是 32 的因數的 k，重複應用 $\psi ^{-1}$ 直到對於某些 j， $(\psi ^{-1})^{(j)}(k)$ 是 32 的因數，然後可以尋找 $S((\psi ^{-1})^{(j)}(k))$ 並且可以使用等式 (5) 計算 S(k)。

總之，比對過程使用 4n+n⌈log2n⌉/8 $4n+n[log_{2}n]/8$ 位，或 n 個位元組用於小於 4 Gb 的基因組。這包括原始和反向基因組的 BWT 字串、部分出現陣列和部分字尾陣列的記憶體; 内存=>zh-mo:記憶體。此外，堆、快取和其他數據結構需要幾百兆位元組的記憶體; 内存=>zh-mo:記憶體。

Illumina read的其他實際問題

不明確的鹼基

read上的非 A/C/G/T 鹼基被簡單地視為不配對，這在演算法中是隱含的（圖 3）。參考基因組上的非 A/C/G/T 鹼基被轉換為隨機核苷酸。這樣做可能會導致錯誤命中充滿模糊鹼基的區域。幸運的是，鑑於相對較長的讀取時間，發生這種情況的可能性非常小。嘗試了 200 萬個 32 bp 的讀數，但沒有看到任何讀數偶然對映到 poly-N 區域。

雙端對映

BWA 支援雙端對映。它首先找到所有好的命中的位置，根據染色體坐標對它們進行排序，然後對所有潛在的命中進行線性掃描以配對兩端。計算所有染色體坐標需要經常尋找字尾陣列。這種配對過程非常耗時，因為使用上述方法即時生成完整的字尾陣列非常昂貴。為了加速配對，快取了大的間隔。此策略將花費在配對上的時間減半。

在配對中，BWA 批次處理 256K 讀取對。在每個批次中，BWA 將完整的 BWA 索引載入到記憶體; 内存=>zh-mo:記憶體中，為每次出現生成染色體坐標，從兩端對映的對映質素高於 20 的讀取對估計插入大小分布，然後將它們配對。之後，BWA 從記憶體; 内存=>zh-mo:記憶體中清除 BWT 索引，載入 2 位編碼參考序列，並對未對映的讀取執行 Smith-Waterman 對齊，其配偶可以可靠地對齊。 Smith-Waterman 對齊挽救了一些差異過大的讀取。

確定允許的最大差異數

給定長度為 m 的讀取，BWA 僅容忍最多具有 k 個差異（錯配或間隙）的命中，其中選擇 k 使得 <4% 的 m 長讀取和 2% 的統一鹼基錯誤率可能包含大於 k 的差異.使用此組態，對於 15–37 bp 讀取，k 等於 2；對於 38–63 bp，k=3；對於 64–92 bp，k=4；對於 93–123 bp，k=5；對於 124–156 bp 讀數，k=6。

生成對映質素分數

對於每個比對，BWA 計算一個對映質素分數，這是比對不正確的 Phred 標度機率。該演算法與 MAQ 相似，不同之處在於在 BWA 中假設總能找到真正的命中。進行此修改是因為知道 MAQ 的公式高估了錯過真正命中的機率，從而導致低估了對映質素。模擬表明，由於這種修改，BWA 可能會高估對映質素，但偏差相對較小。例如，BWA 在對映質素為 60 的 1 569 108 個模擬 70 bp 讀數中錯誤地對齊了 11 個讀數。這些 Q60 對映的錯誤率 7 × 10−6 (= 11/1 569 108) 高於理論預期 10^-6。

SOLiD read的比對

對於 SOLiD 讀取，BWA 將參考基因組轉換為二核苷酸「顏色」序列，並為顏色基因組構建 BWT 索引。讀取被對映在顏色空間中，其中序列的反向二補碼與反向相同，因為顏色的二補碼是它本身。對於SOLiD配對末端作圖，如果兩種情況中的任何一種為真，則稱讀數對處於正確的方向：（i）兩端都對映到基因組的正向鏈，其中R3讀數具有較小的坐標； (ii) 兩端對映到基因組的反向鏈，F3 讀數具有較小的坐標。 Smith-Waterman 對齊也在色彩空間中完成。

比對後，BWA 使用動態編程將顏色讀取序列解碼為核苷酸序列。給定核苷酸參考子序列 $b_{1}b_{2}...b_{l+1}$ 和對映到該子序列的顏色讀取序列 $c_{1}c_{2}...c_{l}$ ，BWA 推斷出核苷酸序列 ${\hat {b}}_{1}{\hat {b}}_{2}...{\hat {b}}_{l+1}$ 。

$\textstyle \sum _{i=1}^{l+1}\displaystyle q^{\prime }\cdot (1-\delta _{{\hat {b}}_{i},b_{i}})+\textstyle \sum _{i=1}^{l}\displaystyle q_{i}[1-\delta _{{\hat {c}}_{i},g({\hat {b}}_{i},{\hat {b}}_{i+1})}]$

其中 q' 是突變的 Phred 標度機率，q_i 是顏色 c_i 的 Phred 質素，函數 g(b, b')=g(b', b) 給出對應於兩個相鄰核苷酸 b 和 b' 的顏色。本質上，如果 $b_{i}\neq {\hat {b}}_{i}$ ，將支付懲罰 q′。和一個懲罰 qi 如果 $c_{i}\neq g({\hat {b}}_{i},{\hat {b}}_{i+1})$ 。

這種最佳化; 优化=>zh-mo:最佳化可以通過動態編程來完成，因為超出位置 i 的最佳解碼僅取決於選擇 ${\hat {b}}_{i}$ 。 Let $f_{i}({\hat {b}}_{i})$ 是達到 i 的最佳解碼分數。迭代方程為

$f_{1}({\hat {b}}_{1})=q^{\prime }\cdot (1-\delta _{{\hat {b}}_{i},b_{i}})$

$f_{i+1}({\hat {b}}_{i+1})=min_{{\hat {b}}_{i}}\left\{f_{i}({\hat {b}}_{i})+q^{\prime }\cdot (1-\delta _{{\hat {b}}_{i+1},b_{i+1}})+q_{i}\cdot [1-\delta _{c_{i},g({\hat {b}}_{i},{\hat {b}}_{i+1})}]\right\}$

BWA 近似基本質素如下。 Let 一個包含圖片、插圖等的外部檔案。

物件名稱是 btp324i11.jpg $c_{i}=g({\hat {b}}_{i},{\hat {b}}_{i+1})$ 。第 i 個基本質素 ${\hat {q}}_{i}$ ，i=2…l，計算如下：

$f(n)={\begin{cases}q_{i-1}+q_{i}&{\text{if }}c_{i-1}={\hat {c}}_{i-1}{\text{ and }}c_{i}={\hat {c}}_{i}\\q_{i-1}-q_{i}&{\text{if }}c_{i-1}={\hat {c}}_{i-1}{\text{ but }}c_{i}\neq {\hat {c}}_{i}\\q_{i}-q_{i-1}&{\text{if }}c_{i}={\hat {c}}_{i}{\text{ but }}c_{i-1}\neq {\hat {c}}_{i-1}\\0&{\text{otherwise}}\\\end{cases}}$

BWA 輸出序列 ${\hat {b}}_{2}....{\hat {b}}_{l}$ ，質素為 ${\hat {q}}_{2}....{\hat {q}}_{l}$ 作為 SOLiD 對映的最終結果。

結果

說明

基於參考基因組的 BWT 實現了 BWA 來進行短讀比對。為單端讀取執行缺口比對，支援雙端比對，生成比對質素並在需要時提供多次命中。預設的輸出對齊格式是 SAM（序列比對/對映格式）。可以使用 SAMtools (http://samtools.sourceforge.net) 來提取區域中的比對、合併/排序比對、獲得單核苷酸多型性 (SNP) 和 indel 呼叫並視覺化比對。

文件; 文档=>zh-mo:文件和原始碼可在 MAQ 網站免費獲得：http://maq.sourceforge.net。

評估

為了評估 BWA 的效能，測試了另外三個比對程式：MAQ（Li 等人，2008a）、SOAPv2（http://soap.genomics.org.cn）和 Bowtie（Langmead 等人，2009 年）。 MAQ 索引使用雜湊; 哈希=>zh-mo:雜湊表read並掃描基因組，是之前為大規模read比對開發的軟件包。 SOAPv2 和 Bowtie 是另外兩個基於 BWT 的短讀比對軟件。最新的 SOAP-2.1.7（Li 等人，未公開數據）使用 2way-BWT（Lam 等人，未公開數據）進行對齊。可以接收超過 35 bp 種子序列的更多錯配，並支援僅限於一個缺口開放的缺口比對。 Bowtie（版本 0.9.9.2）部署了與 BWA 類似的演算法，並沒有通過用 D(i) 限制搜尋來減少搜尋空間，而是通過巧妙地對原始和反向讀取序列進行對齊以繞過對字首樹的根進行不必要的搜尋。預設情況下，Bowtie 對字首 trie 執行 DFS，並在找到第一個符合條件的命中時停止。因此，即使其播種策略被禁用，它也可能會錯過最佳的不精確命中。可以通過在命令行應用「--best」來使 Bowtie 執行 BFS，但這會使 Bowtie 變慢。

包括 BWA 在內的所有四個程式都會在多個同樣最佳的位置上隨機放置重複讀取。由於主要對實踐中的置信比對感興趣，因此需要排除重複命中。 SOAPv2 給出讀取的同等最佳命中數。僅保留唯一對映。還要求 SOAPv2 將可能的間隙大小限制為最多 3 bp。使用命令行選項「--best -k 2」執行 Bowtie，使 Bowtie 輸出讀取的前兩個命中。如果第二個最佳命中包含與最佳命中相同數量的不配對，則丟棄讀取對齊。 MAQ 和 BWA 生成對映質素。使用 MAQ 的比對質素閾值 1 和 BWA 的 10 來確定可信對映。使用不同的閾值是因為知道 MAQ 的對映質素被低估了，而 BWA 的對映質素被高估了。

模擬數據評估

使用 SAMtools 包中包含的 wgsim 程式模擬來自人類基因組的讀數，並執行四個程式將讀數比對回人類基因組。由於知道每個read的確切坐標，能夠計算比對錯誤率。

表 1 顯示 BWA 和 MAQ 實現了相似的比對精度。在可信對映讀取的比例和可信對映的錯誤率方面，BWA 比 Bowtie 和 SOAPv2 更準確。請注意，SOAP-2.1.7 針對長度超過 35 bp 的讀取進行了最佳化; 优化=>zh-mo:最佳化。對於 32 bp 讀取，SOAP-2.0.1 優於最新版本。

表1：模擬數據的評估結果
	Single-end			Paired-end
Program	Time (s)	Conf (%)	Err (%)	Time (s)	Conf (%)	Err (%)
Bowtie-32	1271	79.0	0.76	1391	85.7	0.57
BWA-32	823	80.6	0.30	1224	89.6	0.32
MAQ-32	19797	81.0	0.14	21589	87.2	0.07
SOAP2-32	256	78.6	1.16	1909	86.8	0.78
Bowtie-70	1726	86.3	0.20	1580	90.7	0.43
BWA-70	1599	90.7	0.12	1619	96.2	0.11
MAQ-70	17928	91.0	0.13	19046	94.6	0.05
SOAP2-70	317	90.3	0.39	708	94.5	0.34
bowtie-125	1966	88.0	0.07	1701	91.0	0.37
BWA-125	3021	93.0	0.05	3059	97.6	0.04
MAQ-125	17506	92.7	0.08	19388	96.3	0.02
SOAP2-125	555	91.5	0.17	1187	90.8	0.14
分別從人類基因組模擬了 100 萬對 32、70 和 125 bp 讀數，SNP 突變率為 0.09%，indel 突變率為 0.01%，均勻測序鹼基錯誤率為 2%。 32 bp 讀數的插入大小來自常態分佈 N(170, 25)，70 和 125 bp 讀數來自 N(500, 50)。表中顯示了 2.5 GHz Xeon E5420 處理器單核上的 CPU 時間（時間）、置信對映讀取百分比 (Conf) 和置信對映錯誤對齊百分比 (Err)。

在速度上，SOAPv2 是最快的，如果禁用間隙對齊，雙端對映實際上會快 30-80%。帶有預設選項（數據未顯示）的 Bowtie 比單端對映上的當前設置「–best -k 2」快幾倍。但是，速度的提高是以犧牲準確性為代價的。例如，使用預設選項，Bowtie 可以在 151 秒內對映 200 萬個單端 32 bp 讀數，但 6.4% 的置信對映是錯誤的。這種高對齊錯誤率可能會使結構變異的檢測複雜化，並可能影響 SNP 的準確性。在 BWA 和 MAQ 之間，BWA 快 6–18 倍，具體取決於讀取長度。 MAQ 的速度不受讀取長度的影響，因為它在內部將所有讀取視為 128 bp。通過不檢查類似於 Bowtie 和 SOAPv2 所做的次優命中，可以加速 BWA。然而，在這種情況下，計算對映質素是不可能的，相信生成適當的對映質素對各種下游分析有用，例如檢測結構變化。

在記憶體; 内存=>zh-mo:記憶體上，SOAPv2 使用 5.4 GB。 Bowtie 和 BWA 都使用 2.3 GB 用於單端對映，大約 3 GB 用於雙端對映，比 MAQ 的記憶體; 内存=>zh-mo:記憶體佔用 1 GB 大。然而，所有三個基於 BWT 的對齊器的記憶體; 内存=>zh-mo:記憶體使用與要對齊的讀取數量無關，而 MAQ 在其中是線性的。此外，所有基於 BWT 的對齊器都支援多線程，這減少了多核電腦; 计算机=>zh-mo:電腦上每個 CPU 內核的記憶體; 内存=>zh-mo:記憶體。在現代電腦; 计算机=>zh-mo:電腦伺服器上，記憶體; 内存=>zh-mo:記憶體不是基於 BWT 的對齊器的實際問題。

真實數據評估

為了評估真實數據的效能，從歐洲讀取檔案 (AC:ERR000589) 下載了大約 1220 萬對 51 bp 讀取。這些讀數由 Illumina 為 NA12750 生成，NA12750 是 1000 基因組計劃 (http://www.1000genomes.org) 中的一個男性。讀數被對映到人類基因組 NCBI build 36。表 2 顯示幾乎所有來自 MAQ 和 BWA 的可信比對都以一致的對存在，儘管 MAQ 給出的可信比對較少。較慢的 BWA 模式（無種子；即使最高命中是重複也搜尋次優命中）效果更好。在該模式下，BWA 在 6.3 小時內自信地對映了所有讀數的 89.2%，並以一致對的方式準確對映了 99.2%。

表2：真實數據評估
Program	Time (h)	Conf (%)	Paired (%)
Bowtie	5.2	84.4	96.3
BWA	4.0	88.9	98.8
MAQ	94.9	86.1	98.7
SOAP2	3.4	88.3	97.5
1220 萬個讀取對被對到人類基因組。 2.5 GHz 至強 E5420 處理器單核上的 CPU 時間（時間）、置信度對映讀取百分比 (Conf) 和以正確方向對映且在 300 bp 範圍內的配對置信度對映（配對），顯示在桌子。

在這個實驗中，如果間隙對齊被禁用，SOAPv2 的速度將提高兩倍，置信度對映 (Conf) 和配對百分比 (Paired) 均下降 1%。相比之下，BWA 僅執行無間隙對齊時的速度是它的 1.4 倍。但即使禁用了基於 BWT 的缺口比對，BWA 仍然能夠通過 Smith-Waterman 比對恢復許多短插入缺失，給出配對末端讀數。

還獲得了雞基因組序列（2.1 版），並將這 1220 萬個讀取對與人雞雜交參考序列進行比對。與純人類對齊相比，置信度對映百分比幾乎沒有變化。至於對映到雞基因組的讀取數，Bowtie 將 2640、BWA 2942、MAQ 3005 和 SOAPv2 對映到錯誤的基因組的讀取數為 4531。如果考慮雞序列占人雞雜交參考的四分之一，則 BWA 的比對錯誤率約為 0.06%（=2942×4/12.2M/0.889）。請注意，這種比對錯誤率的估計可能會被低估，因為錯誤比對的人類read往往與人類的重複序列相關，並比對回人類序列。由於存在高度保守的序列和人類組裝不完整或雞基因組組裝錯誤，估計也可能被高估。

如果想要更少的錯誤，BWA 和 MAQ 生成的對映質素允許選擇精度更高的對齊。如果將確定 BWA 的置信命中的比對質素閾值提高到 25，則 86.4% 的讀數可以與比對到雞基因組的 1927 個讀數可靠地對齊，在百分比置信對映和比對到的讀數數量方面均優於 Bowtie錯誤的基因組。

討論

對於與人類參考基因組的短讀比對，BWA 比 MAQ 快一個數量級，同時實現相似的比對精度。它支援單端讀取的間隙對齊，當讀取變得更長並且傾向於包含插入缺失時，這一點變得越來越重要。 BWA 以 SAM 格式輸出對齊，以利用 SAMtools 中實施的下游分析。 BWA 加 SAMtools 提供了 MAQ 包的大部分功能和附加功能。

與速度、記憶體; 内存=>zh-mo:記憶體和比對讀取的數量相比，在真實數據上評估對齊精度要困難得。本文使用了三個標準來評估比對軟件的準確性。第一個標準只能用模擬數據進行評估，它是置信比對數量和置信比對中比對錯誤率的組合。請注意，僅可信比對的數量可能不是一個好的標準：可以以犧牲準確性為代價比對更多。第二個標準是比對read的數量和比對成一致對的read數量的組合，它適用於真實數據，假設來自實驗的交配資訊是正確的，並且結構變異很少。儘管此標準與比對軟件定義配對「一致性」的方式有關，但根據經驗，如果配對參數設置正確，則具有很高的資訊量。第三個標準是如果有意添加來自不同物種的參考序列，read比對到錯誤參考序列的比例。

雖然理論上 BWA 適用於任意長的讀取，但它的效能在長read上會下降，尤其是在測序錯誤率很高的情況下。此外，BWA 總是需要比對完整的read，從第一個鹼基到最後一個（即相對於read的全域比對），但較長的read更有可能被參考基因組中的結構變異或錯誤組裝中斷，這將BWA 失敗。對於長讀段，可能更好的解決方案是將讀段分成多個短片段，用上述演算法分別比對片段，然後加入部分比對以獲得讀段的完整比對。

BWA-SW方法

BWA-SW演算法概述

BWA-SW 為參考序列和查詢序列構建 FM 索引。它隱含地表示字首樹中的參考序列，並表示字首有向無環詞圖（字首 DAWG；Blumer 等人，1985）中的查詢序列，該圖是從查詢序列的字首樹轉換而來的（演算法圖）。通過分別遍歷參考字首 trie 和查詢 DAWG，可以在 trie 和 DAWG 之間應用動態編程。如果不應用啟發式演算法，這種動態編程將找到所有本地配對項，但不會比 BWT-SW 快。在 BWA-SW 中，應用了兩個啟發式規則來大大加速這個過程。首先，查詢 DAWG 上的遍歷是在外迴圈中進行的，因此在不完成動態規劃的情況下，知道參考字首樹中所有與查詢節點配對且得分為正的節點。基於對真實對齊往往具有高對齊分數的觀察，可以在每個節點上修剪低分數配對以限制僅圍繞良好配對的動態規劃。動態規劃的規模因此可以顯着減少。在這個過程中可能會修剪真正的比對，但實際上，這可以通過使用啟發式方法來控制，並且很少發生，因為給定長或高質素的查詢序列。其次，BWA-SW 只報告在查詢序列上基本上不重疊的比對，而不是給出所有重要的局部比對。它啟發式地辨識並丟棄包含在較長對齊中的種子，從而節省不成功種子擴充的計算時間。

符號和定義

字尾陣列和 BWT

設 Σ={A, C, G, T} 是核苷酸字母表，$ 是一個在字典上比 Σ 中的所有符號都小的符號。給定一個核苷酸序列 X=a1…an−1 和 an−1=$，讓 X[i]=ai 是第 i 個符號，X[i, j]=ai…aj 是 X 和 Xi=X 的子序列 [i, n−1] X 的字尾。 X 的字尾陣列 S 是整數 0,..., n−1 的排列，使得 S(i)=j 若且唯若 Xj 是第 i 個字典序最小的字尾。 X 的 BWT 是 X 的排列，如果 S(i)=0，則 B[i]=$，否則 B[i]=X[S(i)−1]。

字尾陣列間隔

給定一個序列 W，字尾陣列 interval 或 SA interval ${\textstyle [{\underline {R}}(W),{\overline {R}}(W)]}$ 的 W 定義為

${\underline {R}}(W)=min\left\{k:W\ is\ the\ prefix\ of\ X_{S(k)}\right\}$

${\overline {R}}(W)=max\left\{k:W\ is\ the\ prefix\ of\ X_{S(k)}\right\}$

字尾陣列間隔和序列比對

不精確配對：有界遍歷/回溯

調頻指數

字尾陣列 S、陣列 C 和 O 足以在 X 中精確搜尋模式。 FM-index (Ferragina and Manzini, 2000) 是三個陣列的壓縮表示，由壓縮的 BWT 字串 B、計算 O 和字尾陣列 S 的一部分。然而，BWA-SW 使用簡化的 FM 索引，其中不壓縮 B 並在沒有輔助數據結構的情況下儲存部分出現陣列 O。簡化版本對於字母表非常小的 DNA 序列更有效。之前的論文 (Li and Durbin, 2009) 中介紹了有關構造的詳細資訊。

比對

對齊是一個元組 (W1, W2, A)，其中 W1 和 W2 是兩個字串，A 是一系列將 W2 轉換為 W1 的複製、替換、插入和刪除操作。插入和刪除是間隙。差距和替代是差異。對齊的編輯距離等於 A 中差異的總數。

可以為給定評分系統的比對計算分數。如果 W1 和 W2 可以與正分數對齊，說 W1 配對 W2，在這種情況下，也說 (W1, W2) 是配對的。

如果 W1 是 W'1 的子串，則稱配對 (W1, W2) 包含在第一個序列的 (W'1, W'2) 中。同樣，可以定義比對（更強的條件）之間以及比對和配對之間的「包含」關係。

字首樹和字首 DAWG

字串 X 的字首 trie 是一棵樹，每條邊都標有一個符號，這樣從葉子到根的路徑上的符號串聯給出了唯一的 X 字首。從節點到根的邊符號串聯為總是 X 的子串，稱為節點所代表的串。節點的SA區間定義為該節點所代表的字串的SA區間。不同的節點可能有相同的區間，但是回顧SA區間的定義，知道這些節點所代表的字串一定是同一個字串的字首，長度不同。

X 的字首 DAWG 通過摺疊具有相同間隔的節點從字首 trie 轉換而來。因此，在字首 DAWG 中，節點和 SA 區間具有一一對應的關係，並且一個節點可能代表 X 的多個子串，落入序列中，其中每個子串都是下一個的字首，如前一段所述。圖 1 給出了一個例子。

將字首樹與字首 DAWG 對齊

為查詢序列 W 構建了一個字首 DAWG 𝒢(W)，為參考 X 構建了一個字首 trie 𝒯(X)。用於計算 W 和 X 之間的最佳分數的動態規劃如下。當u是𝒢(W)的根，v是𝒯(X)的根時，設Guv=Iuv=Duv=0。在 𝒢(W) 中的節點 u，對於它的每個父節點 u′，計算

$I_{uv|u^{\prime }}=max\left\{I_{u^{\prime }v},G_{u^{\prime }v}-q\right\}-r$

$D_{uv|u^{\prime }}=max\left\{D_{uv^{\prime }},G_{uv^{\prime }}-q\right\}-r$

$G_{uv|u^{\prime }}=max\left\{G_{u^{\prime }v^{\prime }}+S(u^{\prime },u;v^{\prime },v),I_{u^{\prime }v^{\prime }},D_{u^{\prime }v^{\prime }},0\right\}-r$

其中 v′ 是 𝒯(X) 中 v 的父級，函數 S(u′, u; v′, v) 給出邊上的符號 (u′, u) 和邊上的符號 (v′, v) 和 q 和 r 分別是間隙開放和間隙擴充懲罰。 Guv、Iuv 和 Duv 的計算公式為：

$u^{*}=argmax_{u^{\prime }\in {\text{pre}}(u)}G_{u}v|u^{\prime }$

$(G_{uv},I_{uv},D_{uv})={\begin{cases}(G_{uv|u^{*}},I_{uv|u^{*}},D_{uv|u^{*}})&G_{uv|u^{*}}>0\\(-\infty ,-\infty ,-\infty )&{\text{otherwise}}\end{cases}}$

其中 pre(u) 是 u 的父節點集。 Guv 等於 u 表示的（可能是多個）子串與 v 表示的（一個）子串之間的最佳分數。如果 Guv>0，說節點 v 與 u 配對。

動態規劃是通過以反向後序（即所有父節點在子節點之前存取）以巢狀方式遍歷𝒢(W) 和 𝒯(X) 來執行的。注意到一旦u不配對v，u不配對v的任何節點，只需要存取靠近𝒯(X)根節點的節點，無需遍歷整個trie，相比之下大大減少了迭代次數到標準的 Smith-Waterman 演算法，該演算法總是通過整個參考序列。

標準 Smith-Waterman 的加速

與時間複雜度為 O(|X|·|W|) 的標準 Smith-Waterman 對齊相比，BWA-SW 具有更好的時間複雜度，因為它並不比時間複雜度為 O(|X|^0.628|W|) 的 BWT-SW 慢。（Lam 等人，2008 年）。這個結論是因為對於短序列比對，正在考慮多個可能的配對與單個 uv 比較。然而，由於與遍歷字首 trie 和字首 DAWG 相關的複雜步驟，與每次迭代相關的常數要大得多，這使得當使用 BWA-SW 擴充唯一對齊時 BWA-SW 效率低下。更有效的策略是使用 BWA-SW 尋找部分配對並應用 Smith-Waterman 演算法進行擴充。在動態規劃中，知道在任何對中考慮的部分配對的數量，因為這可以從 SA 間隔的大小計算。當Guv足夠好並且v的SA區間大小低於某個閾值（預設為3）時，儲存; 保存=>zh-mo:儲存(u,v)對，稱為種子區間對，不要從𝒯中的v節點往更深處（X）。通過尋找 X 和 W 的字尾陣列，可以從種子間隔對中匯出種子配對，或簡稱為種子。這些種子隨後由 Smith-Waterman 演算法擴充。如果整個查詢是一個高度重複的序列，它將完全按照上一節中描述的演算法進行對齊，而不使用 Smith-Waterman 擴充。

因為提前停止動態編程以生成種子，全域最佳對齊可能包含多個種子，實際上這往往是長對齊的情況。通常對於 1 kb 比對，將有 10-20 個種子。下面將利用這一觀察結果啟發式地加快搜尋速度。

BWT-SW 部署了類似的策略，在序列和字首嘗試之間執行動態編程以尋找種子配對，然後是 Smith-Waterman 擴充。與的演算法的主要區別在於，無論 SA 間隔大小如何，一旦分數足夠高，BWT-SW 就會啟動 Smith-Waterman 對齊。有時，重複序列可能與人類基因組中的數千個位置配對，每次擴充部分配對可能很慢。

啟發式加速

Z-最佳策略

到目前為止描述的演算法是精確的，因為它能夠提供與 Smith-Waterman 演算法相同的結果。雖然在給定長參考序列的情況下，它比標準演算法快得多，但對於比對大規模測序數據還不夠快。更仔細的調查表明，即使對於唯一的 500 bp 查詢序列，𝒯(X) 中的幾百萬個節點也可能與具有正比對分數的查詢配對。這些節點中的大多數是隨機配對或在短的低複雜度區域中配對。參觀所有這些都是浪費。

為了加速對齊，在外迴圈中遍歷𝒢(W)，在內迴圈中遍歷𝒯(X)，並且在𝒢(W) 中的每個節點u 處，只保留𝒯(X) 中配對的前Z 個最佳得分節點u，而不是保留所有配對的節點。這種啟發式策略稱為 Z-best。當然，當應用 Z-best 策略時，當錯誤配對具有更高的分數時，可能會錯過包含在真實對齊中的種子。但是，如果查詢與參照幾乎相同，則這種情況發生的頻率較低。另外，如果真對齊很長並且包含很多種子，那麼所有種子都是假的可能性很小。在模擬和真實數據（結果部分）中，發現即使 Z=1 也適用於高質素的 200 bp 讀數（<5% 測序錯誤率）。將 Z 增加到 10 或更高會略微提高準確性，但會大大降低對齊速度。

為了減少對齊錯誤，除了前向對齊之外，還將反向查詢序列與反向參考序列對齊，即反向反向對齊。理想情況下，正向-正向和反向-反向對齊應該產生相同的結果，但如果真正對齊中的種子具有低分字尾（或字首），則正向-正向（或反向-反向）對齊可能會錯過它，同時結合兩輪對齊減少機會。此外，如果前向對齊中的最佳對齊包含許多種子配對，則其為假的可能性也很小。在實現中，如果最佳對齊預設包含 5 個或更多種子，則不應用反向-反向對齊。

在 Smith-Waterman 擴充之前過濾種子

與 BLAST 一樣，BLAT 和 SSAHA2 都報告了所有重要的比對或通常是數十個得分最高的比對，但這並不是讀取對映中最需要的輸出。通常對覆蓋查詢序列的每個區域的最佳比對或最好的少數比對更感興趣。例如，假設一個 1000 bp 的查詢序列由來自一條染色體的 900 bp 片段和來自另一條染色體的 100 bp 片段組成； 900 bp 片段中的 400 bp 是高度重複的序列。對於 BLAST，要知道這是一個嵌合讀取，需要要求它報告 400 bp 重複的所有比對，這既昂貴又浪費，因為通常對包含在更長的唯一序列中的短重複序列的比對不感興趣序列。在這個例子中，一個有用的輸出是分別報告 900 bp 和 100 bp 片段的一個比對，並指出這兩個片段是否有可能使最佳比對不可靠的次優比對。這樣的輸出簡化了下游分析並節省了重建重複序列的詳細比對的時間。

在 BWA-SW 中，如果查詢上重疊區域的長度小於較短查詢段長度的一半，就說兩個對齊是不同的。的目標是找到一組不同的對齊方式，使集合中每個對齊方式的分數總和最大化。這個問題可以通過動態編程解決，但在的例子中，讀取通常是完全對齊的，貪婪的近似會很好地工作。在實際實現中，根據對齊分數對局部對齊進行排序，從最好的一個掃描排序列表，如果它與所有保留的分數較大的對齊不同，則保留一個對齊；如果對齊 a2 被拒絕，因為它與 a1 沒有區別，認為 a2 是 a1 的次優對齊，並使用此資訊來近似對映質素。

因為只保留在查詢序列上基本不重疊的比對，所以不妨丟棄對最終比對沒有貢獻的種子。可以在 Smith-Waterman 擴充之前使用另一種啟發式方法檢測此類種子，從而可以節省花費在不必要擴充上的時間。為了辨識這些種子，將包含在條帶中的種子連結起來（預設條頻寬度為 50 bp）。如果在查詢序列上，一條短鏈完全包含在一條長鏈中，並且短鏈中的種子數低於長鏈中種子數的十分之一，將丟棄短鏈中的所有種子，基於觀察到在這種情況下，短鏈很少能比長鏈產生更好的對齊。與 Z-best 策略不同，這種啟發式對對齊精度沒有明顯影響。在 1000 個 10 kb 模擬數據上，它將執行時間減半，而不會降低精度。

近似對映質素

Li，H。等。 (2008) 引入了對映質素的概念來估計查詢序列被放置在錯誤位置的機率。如果對齊演算法保證找到所有局部對齊，則對映質素僅由這些局部對齊決定。然而，隨着 BWA-SW 部署啟發式規則，產生錯誤對齊的機會也與啟發式有關。為了估計 BWA-SW 比對的對映質素，擬合了一個經驗公式：250·c1·c2·(S1−S2)/S1，其中 S1 是最佳比對的得分，S2 是次優比對的得分，如果對齊覆蓋超過四個種子，則 c1 等於 1，否則為 0.5，如果通過正向-正向和反向-反向對齊找到最佳對齊，則 c2 等於 1，否則為 0.2。

結果

簡介

BWA-SW 演算法是作為 BWA 程式（Li 和 Durbin，2009）的一個組件實現的，該程式在 GNU 通用公共許可證 (GPL) 下分發。該實現將 BWA 索引和查詢 FASTA 或 FASTQ 檔案作為輸入，並以 SAM 格式輸出對齊（Li et al., 2009）。查詢檔案通常包含許多序列（讀取）。依次處理每個查詢序列，如果適用，使用多個線程。記憶體; 内存=>zh-mo:記憶體使用由 FM-index 主導，人類基因組約為 3.7 GB。每個查詢所需的記憶體; 内存=>zh-mo:記憶體大致與序列長度成正比。在典型的測序讀取中，總記憶體; 内存=>zh-mo:記憶體小於 4 GB；在一個具有 100 萬個鹼基對 (Mbp) 的查詢序列上，峰值記憶體; 内存=>zh-mo:記憶體總共為 6.4 GB。

在實現中，嘗試根據讀取長度和測序錯誤率自動調整參數，使預設設置適用於不同特徵的輸入。這種行為對於不熟悉演算法的用戶來說很方便，並且有助於提高讀取長度和錯誤率混合的效能。

模擬數據評價

在模擬數據上，從比對中知道正確的染色體坐標，並且評估很簡單。

整體表現

表 1 顯示了給定不同讀取長度和錯誤率的 BLAT（v34）、BWA-SW（版本 0.5.3）和 SSAHA2（版本 2.4）的 CPU 時間、可靠對齊讀取的分數和對齊錯誤率。除非必要，嘗試使用每個對齊器的預設命令行選項。根據輸入數據的特徵微調選項可能會產生更好的效能。

表1：模擬數據評價
Program	Metrics	100 bp			200 bp			500 bp			1000 bp			10 000 bp
Program	Metrics	2%	5%	10%	2%	5%	10%	2%	5%	10%	2%	5%	10%	2%	5%	10%
BLAT	CPU sec	685	577	559	819	538	486	1078	699	512	1315	862	599	2628	1742	710
	Q20%	68.7	25.5	3.0	92.0	52.9	7.8	97.1	86.3	21.4	97.7	96.4	39.0	98.4	99.0	94.0
	errAln%	0.99	2.48	5.47	0.55	1.72	4.55	0.17	1.12	4.41	0.01	0.52	3.98	0.00	0.00	1.28
BWA-SW	CPU sec	165	125	84	222	168	118	249	172	152	234	168	150	158	134	120
	Q20%	85.1	62.2	19.8	93.8	88.7	49.7	96.1	95.5	85.1	96.9	96.5	95.0	98.4	98.5	98.1
	errAln%	0.01	0.05	0.17	0.00	0.02	0.13	0.00	0.00	0.04	0.00	0.00	0.00	0.00	0.00	0.00
SSAHA2	CPU sec	4872	7962	9345	1932	2236	5252	3311	8213	6863	1554	1583	3113	–	–	–
	Q20%	85.5	83.8	78.2	93.4	93.1	91.9	96.6	96.5	96.1	97.7	97.6	97.4	–	–	–
	errAln%	0.00	0.01	0.19	0.01	0.00	0.01	0.00	0.01	0.04	0.00	0.00	0.00	–	–	–
從人類基因組模擬了大約 10 000 000 bp 不同讀長和錯誤率的數據。 20% 的錯誤是從幾何分布（密度：0.7·0.3l−1）繪製的插入缺失長度的插入錯誤。這些模擬讀數分別使用 BLAT（選項 -fastMap）、BWA-SW 和 SSAHA2（選項 -454，用於 100 和 200 bp 讀數）與人類基因組對齊。然後將對齊的坐標與模擬坐標進行比較以尋找對齊錯誤。在此表的每個儲存格中，三個數字是 Intel E5420 2.5 GHz CPU 單核上的 CPU 秒數、對映質素大於或等於 20 (Q20) 的對齊百分比以及 Q20 對齊中的錯誤對齊百分比 . SSAHA2 和 BWA-SW 報告對映質素； BLAT 對映質素估計為最佳和次佳對齊分數之差除以最佳對齊分數的 250 倍（與 BWA-SW 的計算基本相同）。

從表 1 可以看出，BWA-SW 顯然是最快的，在所有輸入上都比 BLAT 和 SSAHA2 快幾倍，並且其速度對讀取長度或錯誤率不敏感。當查詢較長或錯誤率較低時，BWA-SW 的準確率與 SSAHA2 相當。鑑於短且容易出錯的讀取，SSAHA2 更準確，儘管它必須花費更多時間來對齊此類讀取。 SSAHA2 未在 10 kb 讀取上進行測試，因為它最初不是為此任務設計的，因此效能不佳。帶有 -fastMap 選項的 BLAT 比 SSAHA2 快，但精度較低。在預設選項下，BLAT 比 SSAHA2 慢幾到幾十倍。與 -fastMap 模式相比，精度更高，但總體上仍低於 BWA-SW（數據未顯示）。

在記憶體; 内存=>zh-mo:記憶體上，BWA-SW 和 BLAT 都使用 ∼4 GB 記憶體; 内存=>zh-mo:記憶體。 SSAHA2 使用 2.4 GB 使用預設選項進行≥500 bp 的讀取，使用 5.3 GB 使用 -454 選項進行更短的讀取，這會增加雜湊; 哈希=>zh-mo:雜湊表中儲存的種子序列的數量並因此增加記憶體; 内存=>zh-mo:記憶體。此外，BWA-SW 支援多線程，因此如果在多核電腦; 计算机=>zh-mo:電腦上執行，每個 CPU 內核可能佔用更少的記憶體; 内存=>zh-mo:記憶體。 SSAHA2 和 BLAT 目前不支援多線程。

首先研究給定嵌合讀數的每個對準器的行為。為此，製作了兩個嵌合讀數，它們都由一個染色體位置的 1000 bp 片段和另一個位置的 300 bp 片段組成。兩次讀取之間的主要區別在於，第二次讀取中的 1000 bp 片段具有~750 bp 的重複序列，而第一次讀取是高度獨特的。當將兩個嵌合讀數與人類基因組進行比對時，BWA-SW 會根據需要報告四種比對，每個片段一個。最新的 SSAHA2 無法在兩個讀取中找到 300 bp 片段的對齊，儘管如果將 300 bp 片段作為單獨讀取對齊，它能夠找到對齊。預設情況下，舊版本 (1.0.9) 能夠在第一次讀取中對齊 300 bp 片段，但對於第二次讀取，需要切換到更徹底但速度更慢的組態，將所有命中報告到 750 bp重複。帶有 -fastMap 的 BLAT 沒有找到第二次讀取的 300 bp 片段的對齊。在這兩個例子中，只有 BWA-SW 有足夠的能力來檢測嵌合體。

此外，BWA-SW 很少產生高質素的假嵌合比對。例如，假設 10 000 1 kb 讀數有 10% 的錯誤但在模擬中沒有嵌合體，BWA-SW 預測 18 個嵌合讀數。這些讀取上錯誤對齊的片段的對映質素僅為 2.4（最大 11），這意味着 BWA-SW 意識到這些嵌合體是不可靠的。正如預期的那樣，由於鹼基錯誤較低，BWA-SW 產生的假嵌合讀數較少。

真實數據評估

由於缺乏基本事實，對真實數據的評估變得複雜。然而，仍然可以通過比較兩個對齊器的結果來評估相對準確性，使用真實對齊往往具有相當高的對齊分數的原則，因為大多數錯誤是由於未能找到種子而產生的。

假設使用兩個對齊器 A 和 B 對齊讀取並得到不同的結果。如果 A 和 B 都給出低對映質素，則對齊是不明確的，並且任何對齊是否錯誤都無關緊要。如果 A 的對映質素很高，而 A 的對齊分數比 B 差，則 A 對齊可能是錯誤的；即使A對齊分數比B好一點，A對齊也不可靠，A給出的高對映質素仍然值得懷疑。實際上，定義「稍微好一點」的對齊分數需要對分數差異設置任意閾值，因此這種評估方法是近似的。

表 2 匯總了 454 個讀數，這些讀數僅由一個比對器對映或對映到不同位置，並且被 BWA-SW 或 SSAHA2 指定的對映質素大於或等於 20。可以看到 BWA-SW 往往會錯過高錯誤率的短對齊（其中 946 個），這與對模擬數據的評估一致。 SSAHA2 由於不同的原因錯過了對齊。在 1188 次讀取中，SSAHA2 產生明顯錯誤的比對。通過分配低對映質素意識到這些對齊是錯誤的，但無論如何都錯過了真正的對齊。


Condition	Count	BWA-SW			SSAHA2
Condition	Count	avgLen	avgDiff	avgMapQ	avgLen	avgDiff	avgMapQ
BWA-SW≥20; SSAHA2 unmapped	0	–	–	–	–	–	–
BWA-SW≥20 plausible; SSAHA2<20	1188	398.2	1.3%	178.4	198.3	13.1%	3.9
BWA-SW≥20 questionable	40	183.0	7.8%	41.2	280.3	9.4%	2.4
SSAHA2≥20; BWA-SW unmapped	946	–	–	–	75.4	6.3%	51.2
SSAHA2≥20 plausible; BWA-SW<20	47	129.0	9.3%	2.5	200.5	8.8%	34.4
SSAHA2≥20 questionable	185	400.2	1.7%	13.4	399.2	2.9%	216.4
從 SRR003161 中統一選擇的總共 137 670 454 個讀數分別用 BWA-SW 和 SSAHA2 對映到人類基因組。如果 BWA-SW 和 SSAHA2 對齊的最左側坐標相差超過 355 bp（平均讀取長度），則稱讀取不一致對齊。為每個比對計算一個分數，該分數等於配對數減去三乘以比對區域中的差異（錯配和缺口）數。如果從 BWA-SW 比對得出的分數減去從相同讀取的 SSAHA2 比對得出的分數大於或等於 20（即 BWA-SW 比對足夠好），則認為 BWA-SW 比對是合理的;否則，據說 BWA-SW 對齊是有問題的。合理的和有問題的 SSAHA2 比對以類似的方式定義。表中，『BWA-SW≥20』表示對映質素高於20的BWA-SW比對。BWA-SW總共遺漏了993（=946+47）個比對，其中SSAHA2比對好，而SSAHA2比對好1188； 40 個 BWA-SW Q20 對齊和 185 個 SSAHA2 Q20 對齊可能是錯誤的。

對於這兩個比對器，大多數錯誤對齊是由於忽略了與最佳報告對齊具有相似分數的對齊。例如，SSAHA2 將讀取 SRR003161.1261578 與對映質素為 244 的 X 染色體對齊，而 BWA-SW 將其與具有相同對齊長度和編輯距離的 2 號染色體對齊。兩個最佳評分比對的存在意味着無法唯一放置讀取，並且高達 244 的對映質素是不准確的。 SSAHA2 提供這種高對映質素可能是因為它忽略了染色體 2 上的配對。在這個特定範例中，BWA-SW 正確地給出了一個對映質素為零，儘管它可能會忽略其他範例中的替代配對。

在模擬的 100 和 200 bp 讀數上，帶有 -454 選項的 SSAHA2 比 BWA-SW 提供更好的比對。在這個真實的數據集上，BWA-SW 更準確，可能是因為平均讀取長度相對較長（355 bp）。為了證實這一推測，比較了來自執行 SRR002644 的 99 958 個讀數的兩個對齊器，平均讀數長度為 206 bp。這次 BWA-SW 錯過了 1092 個 SSAHA2 Q20 比對，並產生了 39 個有問題的比對； SSAHA2 漏掉了 325 個，產生了 10 個有問題的。 SSAHA2 在這個較短的數據集上更準確，儘管它比 BWA-SW 慢 9 倍，並且使用的記憶體; 内存=>zh-mo:記憶體多 40%。

討論

BWA-SW 是一種有效的演算法，用於將幾百個或更多鹼基對的查詢序列與長參考基因組進行比對。對於長查詢或低錯誤率的查詢，其敏感性和特異性往往更高，並且在此類查詢序列上，BWA-SW 的準確性可與迄今為止最準確的對齊器相媲美。此外，BWA-SW 能夠檢測可能由結構變化或參考錯誤組裝引起的嵌合體，這可能對 BLAT 和 SSAHA2 構成挑戰。

BWA-SW、BLAT 和 SSAHA2 都遵循種子和擴充範式。主要區別來自播種策略。 BLAT 和 SSAHA2 將短精確配對辨識為種子，通常長度為 11 或 12 bp。對於分別在長度為 L 和 l 的兩個序列之間的 k-mer 種子，種子的預期數量為 L·l/4k，或者是 105 的數量級，用於與人類基因組的比對。用 Smith-Waterman 演算法擴充這些種子是很昂貴的。為了減少不必要的種子擴充，BLAT 和 SSAHA2 預設都使用不重疊的種子，並且需要多個種子配對，這對於隨機序列應該可以很好地工作，但仍然涉及高度重複區域中的許多種子擴充。 BWA-SW 通過在獨特區域使用一些長間隙種子解決了這個問題。在真實的生物數據上，它節省了許多不必要的種子擴充，並帶來了更好的整體效能。然而，為了減少辨識長種子的時間，BWA-SW 只保留了動態規劃矩陣的一小部分，這可能會錯過所有真正配對的種子。這種啟發式是對齊錯誤的主要來源，特別是對於短查詢，當要對齊的序列之間只有很少的有效唯一種子時。幸運的是，在長對齊上，遺失所有種子的可能性很小。已經證明 BWA-SW 與 SSAHA2 一樣有效。

BWA-SW 在幾個方面與 BWT-SW 不同。首先，BWT-SW 保證找到所有本地配對，而 BWA-SW 是一種啟發式演算法，它可能會錯過真正的命中，但速度要快得多。其次，BWA-SW 對齊兩個 FM 索引，而 BWT-SW 對齊一個序列和一個 FM 索引。為查詢序列構建字首 DAWG 可能有助於避免在查詢中重複對齊相同的子字串，從而提高理論時間複雜度。第三，BWA-SW在內迴圈中遍歷參考字首trie，而BWT-SW在內迴圈中遍歷查詢序列。如果沒有啟發式方法，BWA-SW 方法會損害效能，因為必須在遍歷參照字首樹時用速度換取記憶體; 内存=>zh-mo:記憶體，而在外迴圈中遍歷它會更有效。儘管如此，應用 Z-best 策略需要知道與查詢子串配對的得分最高的參考節點，而無需完成動態規劃，因此僅在內部迴圈中遍歷參考時才有效。第四，BWA-SW 只報告與查詢序列基本不重疊的比對，而 BWT-SW 與 BLAST 一樣，報告所有具有統計意義的比對。 BWA-SW 保留了對齊的關鍵資訊，並生成更小更方便的輸出。對於 BWT-SW，終端使用者通常需要對結果進行後處理，以過濾掉許多他們不感興趣的對齊方式。總而言之，鑑於大規模真實數據，BWA-SW 被調整為具有實際實用性。

BWA-SW 的高速主要來自兩種策略：使用 FM 索引和抑制包含在更好配對中的短重複配對。雖然第一種策略不適用於基於雜湊; 哈希=>zh-mo:雜湊表的演算法，例如 SSAHA2 和 BLAT，但第二種策略可以在此類程式中實現，並且可以通過節省大量構建重複對齊的時間來顯着加速它們。雖然 BWT 的使用減少了重複中不必要的對齊，但與雜湊; 哈希=>zh-mo:雜湊表尋找相比，每個 BWT 操作都帶有一個很大的常數。如果基於雜湊; 哈希=>zh-mo:雜湊表的演算法結合了其中的一些功能，那麼它們仍有可能比 BWA-SW 更快。

參考文獻

Altschul SF, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25:3389–3402. [PMC free article] [PubMed] [Google Scholar]
Blumer A, et al. The smallest automaton recognizing the subwords of a text. Theor. Comput. Sci. 1985;40:31–55. [Google Scholar]
Burrows M, Wheeler DJ. Technical report 124. Palo Alto, CA: Digital Equipment Corporation; 1994. A block-sorting lossless data compression algorithm. [Google Scholar]
Campagna D, et al. PASS: a program to align short sequences. Bioinformatics. 2009;25:967–968. [PubMed] [Google Scholar]
Eaves HL, Gao Y. MOM: maximum oligonucleotide mapping. Bioinformatics. 2009;25:969–970. [PubMed] [Google Scholar]
Eid J, et al. Real-time DNA sequencing from single polymerase molecules. Science. 2009;323:133–138. [PubMed] [Google Scholar]
Ferragina P, Manzini G. Proceedings of the 41st Symposium on Foundations of Computer Science (FOCS 2000) IEEE Computer Society; 2000. Opportunistic data structures with applications; pp. 390–398. [Google Scholar]
Ferragina P, Manzini G. Proceedings of the 41st Symposium on Foundations of Computer Science (FOCS 2000) Redondo Beach, CA, USA: 2000. Opportunistic data structures with applications; pp. 390–398. [Google Scholar]
Grossi R, Vitter JS. Proceedings on 32nd Annual ACM Symposium on Theory of Computing (STOC 2000) ACM; 2000. Compressed suffix arrays and suffix trees with applications to text indexing and string matching; pp. 397–406. [Google Scholar]
Hon W.-K, et al. A space and time efficient algorithm for constructing compressed suffix arrays. Algorithmica. 2007;48:23–36. [Google Scholar]
Jiang H, Wong WH. SeqMap: mapping massive amount of oligonucleotides to the genome. Bioinformatics. 2008;24:2395–2396. [PMC free article] [PubMed] [Google Scholar]
Jung Kim Y, et al. ProbeMatch: a tool for aligning oligonucleotide sequences. Bioinformatics. 2009;25:1424–1425. [PMC free article] [PubMed] [Google Scholar]
Kent WJ. BLAT–the BLAST-like alignment tool. Genome Res. 2002;12:656–664. [PMC free article] [PubMed] [Google Scholar]
Lam TW, et al. Compressed indexing and local alignment of DNA. Bioinformatics. 2008;24:791–797. [PubMed] [Google Scholar]
Langmead B, et al. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25. [PMC free article] [PubMed] [Google Scholar]
Lin H, et al. ZOOM! Zillions of oligos mapped. Bioinformatics. 2008;24:2431–2437. [PMC free article] [PubMed] [Google Scholar]
Lippert RA. Space-efficient whole genome comparisons with Burrows-Wheeler transforms. J. Comput. Biol. 2005;12:407–415. [PubMed] [Google Scholar]
Li H, et al. Mapping short DNA sequencing reads and calling variants using mapping quality scores. Genome Res. 2008a;18:1851–1858. [PMC free article] [PubMed] [Google Scholar]
Li R, et al. SOAP: short oligonucleotide alignment program. Bioinformatics. 2008b;24:713–714. [PubMed] [Google Scholar]
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25:1754–1760. [PMC free article] [PubMed] [Google Scholar]
Li H, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–2079. [PMC free article] [PubMed] [Google Scholar]
Ma B, et al. PatternHunter: faster and more sensitive homology search. Bioinformatics. 2002;18:440–445. [PubMed] [Google Scholar]
Malhis N, et al. Slider–maximum use of probability information for alignment of short sequence reads and SNP detection. Bioinformatics. 2009;25:6–13. [PMC free article] [PubMed] [Google Scholar]
Meek C, et al. Proceedings of 29th International Conference on Very Large Data Bases (VLDB 2003) Berlin, Germany: 2003. OASIS: an online and accurate technique for local-alignment searches on biological sequences; pp. 910–921. [Google Scholar]
Morgulis A, et al. Database indexing for production megablast searches. Bioinformatics. 2008;24:1757–1764. [PMC free article] [PubMed] [Google Scholar]
Ning Z, et al. SSAHA: a fast search method for large DNA databases. Genome Res. 2001;11:1725–1729. [PMC free article] [PubMed] [Google Scholar]
Pearson WR, Lipman DJ. Improved tools for biological sequence comparison. Proc. Natl Acad. Sci. USA. 1988;85:2444–2448. [PMC free article] [PubMed] [Google Scholar]
Rumble SM, et al. SHRiMP: accurate mapping of short color-space reads. PLoS Comput. Biol. 2009;5:e1000386. [PMC free article] [PubMed] [Google Scholar]
Schatz M. Cloudburst: highly sensitive read mapping with mapreduce. Bioinformatics. 2009;25:1363–1369. [PMC free article] [PubMed] [Google Scholar]
Smith AD, et al. Using quality scores and longer reads improves accuracy of Solexa read mapping. BMC Bioinformatics. 2008;9:128. [PMC free article] [PubMed] [Google Scholar]
Weese D, et al. RazerS–fast read mapping with sensitivity control. Genome Res. 2009;19:1646–1654. [PMC free article] [PubMed] [Google Scholar]
Zhang Z, et al. A greedy algorithm for aligning DNA sequences. J. Comput. Biol. 2000;7:203–214. [PubMed] [Google Scholar]