生物信息百Jia軟體（一）——bwa

編者按

前面寫了專題《手把手教你生物信息分析平台搭建》，然後又介紹了很多《生物神奇網站》資源，也介紹了《生物信息之獨孤九劍》Linux操作。那麼萬事俱備，就開始學習生物信息吧。所以，我們開始新的篇章——《生物信息百jia軟體》。百Jia是什麼意思呢？可以是百佳，也可以是百家，還可以是百加。從100家中選擇100款優秀軟體，掌握這些軟體，就可以擴展出更多內容，這就是百Jia。

一、功能分類：短序列比對

二、軟體官網：

https://sourceforge.net/projects/bio-bwa/

三、軟體介紹：

隨著以illumina為主的高通量測序技術的流行，短序列比對成為二代測序分析中的核心分析，短序列比對也就是將測序得到的短片段在回帖到基因組上，這個過程也被稱為mapping。mapping的結果包含了所有的信息，後續所有的分析都是基於這個mapping的結果。像目前流行的RNAseq分析，外顯子分析，全基因組WGS等都需要利用短序列比對，其實序列拼接裡面也是首先進行短序列比對。可以說如果想分析二代測序的數據，就必須進行短序列比對。而在眾多的短序列比對軟體中，BWA幾乎已經成為默認的行業標準。隨著三代長片段測序技術的興起，bwa也在進步，最新版本的軟體改進了很多演算法，增加了很多新功能。之前的版本不支持長reads，不能進行split reads的比對，因此在RNAseq中存在可變剪切，不能使用bwa比對。不過目前最新版本的bwa都解決了這些功能。

四、軟體演算法介紹：

目前版本的BWA主要使用三種演算法，分別是BWA-backtrack， BWA-SW和BWA-MEM。其實之前短序列比對主要使用的是bwt演算法，(Burrows–Wheeler transform），bwa軟體的名字也來源於利用bwt演算法的aligner。

首先，BWA-backtrack，軟體中的aln，samse，sampe比對都使用bwa-backtack演算法中，主要適合比較短的reads，尤其是100bp一下的，不過目前很少了；BWA-SW中的SW 表示affine-gap Smith-Waterman)，BWA-MEM中的MEM表示maximal exact matches。

BWA-SW和BWA-MEM都支持長read是和split比對，reads長度可以從70bp到1M，BWA-MEM是更新的演算法，更快也更準確，是最推薦使用的算，也是使用範圍最廣泛的，Illumina, 454, Ion Torrent，Sanger，拼接的contigs 或者BAC 都可以使用這種演算法，對於包含很多gap的情況，bwa-sw會更好些。所以很明顯，目前的版本我們使用bwa-mem演算法就行了。

五、下載安裝

bwa安裝並不難，可以直接下載源代碼進行編譯，只需敲make編譯即可，注意系統環境中需要安裝zlib庫。

也可以直接利用apt工具下載，或者使用bioconda進行安裝。

六、軟體使用：

直接敲bwa，彈出軟體選項參數。比對還是分成兩大步，建立索引，也叫做建庫，然後是比對。首先，利用bwa index可以建立索引，輸入參考序列的fasta格式文件，

-a 指定建立索引的演算法，bwtsw，is或者rb2，以前沒有rb2，而是div。is還是適合小於2G的基因組，bwtwt是大於10M的基因組才行，這個地方不要選錯了，細菌基因組一般都小於10M，只能用is演算法，人基因組是3G，只能用bwtsw，處於中間的選擇哪個都可以。

-p 選項是輸出前綴，一般不用設置，就是在參考序列名字後面多出很多文件。

這樣就開始建立索引，根據文件大小和機器配置不同，需要的時間也不同，對於人的基因組序列，可能需要幾個小時，因此很多網站在提供參考序列的時候，連索引文件也一起提供，下載索引比自己建立索引文件要快很多。所以建立完成之後就可以進行比對了，一般是兩條pairend的reads。當然其他測序平台的單端reads也可以。現在使用bwa-mem演算法只需要兩步就行了。bwa mem有很多選項參數，這裡我們介紹幾個比較重要的。

-R 設置reads標頭，放到一對引號中，也就是sam文件中的RG部分，為什麼要設置RG表頭呢，因為同一樣品可能包括多個測序結果，來自不同lane，不同文庫，或者不同樣品的比對結果合併到同一個文件中進行處理，就需要通過RG進行標記區分。RG每個標記用冒號分割鍵和值，不同標記用 " " 分隔。例如"@RG ID:foo SM:bar LB:library1"

-M 這個選項也比較重要，因為bwa men在比對的時候，對於一條序列同時比對到基因組不同區域的情況，bwa認為都是最優匹配，但是會與Picard tools不兼容，影響後面GATK檢測，這個時候可以設置-M選項，將較短的比對標記為此優，與picard兼容。

-t 設置線程數，多線程可以顯著提高比對效率，因為數據比對大，影響還是非常顯著的，現在整個行業想進各種方法來提高比對效率，包括提升硬體，在CPU底層優化，改進軟體演算法等，因為現在一個完整的人基因組比對還需要幾個小時，GATK還需要一兩天，如果提高計算速度，是非常有意義的。

-T ，給定一個分值，當比對的分值比設置的小時，不輸出該比對結果；

-a 將所有的比對結果都輸出，包括 single-end 和 unpaired paired-end的 reads

-p 智能的pairing，如果不設置，輸入文件只有1個，則進行單端比對；若輸入文件有2個，則作為paired reads進行比對。如果設置-p，則僅以第1個文件作為輸入，輸入的文件若有2個，則被忽略之

下面在看一下bwasw比對，因為使用了 Smith-Waterman，特別適合長reads比對，因為能夠容許更多的gap情況。幫助文檔中默認提示的就是query.fa，一般是長序列。幫助中提示對於長的Illumina 數據, 454，Sanger拼接的contigs或者 fosmids，BACs等，使用默認參數就行，如果是2010年之前Pacbio，設置"-b5 -q2 -r1 -z10"。

七、使用案例：

現有參考序列ref.fna，illumina測序兩條reads，reads.1.fq.gz與reads.2.fq.gz。

八、注意事項：

1、安裝時需要依賴zlib；

2、建庫時注意選擇正確的演算法；

3、最新版本bwa，需要使用-R選項。

---------- END ----------

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！