BWA使用詳解

一缠俺、簡介

BWA，即Burrows-Wheeler-Alignment Tool贷岸。BWA 是一種能夠?qū)⒉町惗容^小的序列比對到一個較大的參考基因組上的軟件包壹士。它由三個不同的算法：

BWA-backtrack: 是用來比對 Illumina 的序列的，reads 長度最長能到 100bp偿警。-
BWA-SW: 用于比對 long-read 躏救，支持的長度為 70bp-1Mbp；同時支持剪接性比對螟蒸。
BWA-MEM: 推薦使用的算法迹缀，支持較長的read長度肮帐，同時支持剪接性比對（split alignments)，但是BWA-MEM是更新的算法，也更快祟剔，更準(zhǔn)確，且 BWA-MEM 對于 70bp-100bp 的 Illumina 數(shù)據(jù)來說愤诱，效果也更好些顷蟀。

對于上述三種算法，首先需要使用索引命令構(gòu)建參考基因組的索引绍赛，用于后面的比對蔓纠。所以，使用BWA整個比對過程主要分為兩步吗蚌，第一步建索引腿倚，第二步使用BWA MEM進(jìn)行比對。

bwa的使用需要兩中輸入文件：

Reference genome data（fasta格式 .fa, .fasta, .fna）
Short reads data (fastaq格式 .fastaq, .fq)

二褪测、BWA的安裝

a.下載BWA http://bio-bwa.sourceforge.net/bwa.shtml

#解壓
$ tar -jxvf bwa-*.tar.bz2
$ cd bwa-*猴誊；

# 編譯BWA
$ make

$ echo 'PATH=$PATH:/path/bwa--*' >> ~/.bashrc
$ source ~/.bashrc

或

git clone https://github.com/lh3/bwa.git
cd bwa
make

安裝好之后潦刃，會出現(xiàn)一個名為bwa的可執(zhí)行文件，輸入下面命令可以查看幫助信息

./bwa

Program: bwa (alignment via Burrows-Wheeler transformation)
Version: 0.7.17-r1188
Contact: Heng Li <<a href="mailto:lh3@sanger.ac.uk" _href="mailto:lh3@sanger.ac.uk">lh3@sanger.ac.uk</a>>
Usage:   bwa <command> [options]
Command: index         index sequences in the FASTA format
         mem           BWA-MEM algorithm
         fastmap       identify super-maximal exact matches
         pemerge       merge overlapping paired ends (EXPERIMENTAL)
         aln           gapped/ungapped alignment
         samse         generate alignment (single ended)
         sampe         generate alignment (paired ended)
         bwasw         BWA-SW for long queries
         shm           manage indices in shared memory
         fa2pac        convert FASTA to PAC format
         pac2bwt       generate BWT from PAC
         pac2bwtgen    alternative algorithm for generating BWT
         bwtupdate     update .bwt to the new format
         bwt2sa        generate SA from BWT and Occ
Note: To use BWA, you need to first index the genome with `bwa index'.
      There are three alignment algorithms in BWA: `mem', `bwasw', and
      `aln/samse/sampe'. If you are not sure which to use, try `bwa mem'
      first. Please `man ./bwa.1' for the manual.

可以看到由很多的子命令構(gòu)成懈叹。

三乖杠、 BWA 的使用

bwa軟件的作用是將序列比對到參考基因組上，在比對之前澄成，首先需要對參考基因組建立索引胧洒。

建立索引 index

在進(jìn)行 reads 的比對前，需要對 fasta 文件構(gòu)建 FM-index墨状。

用法和參數(shù)如下：

Usage：bwa index [ –p prefix ] [ –a algoType ] <in.db.fasta>

Options:
　　-p STR 輸出數(shù)據(jù)庫的前綴卫漫；【默認(rèn)和輸入的文件名一致，輸出的數(shù)據(jù)庫在其輸入文件所在的文件夾肾砂，并以該文件名為前綴列赎。】
　　-a [is|bwtsw] 構(gòu)建index的算法镐确，有以下兩個選項(xiàng)：
　　　　-a is 是默認(rèn)的算法包吝，雖然相對較快，但是需要較大的內(nèi)存源葫，當(dāng)構(gòu)建的數(shù)據(jù)庫大于2GB的時候就不能正常工作了;
　　　　-a bwtsw 對于短的參考序列式不工作的诗越，必須要大于等于10MB, 但能用于較大的基因組數(shù)據(jù)，比如人的全基因組息堂。

bwa index 指令更多的用法及 options嚷狞，通過bwa index 命令來查看

# 根據(jù)reference genome data(e.g. ref.fa) 建立 Index File：
$ bwa index ref.fa -p genome        # 可以不加-p genome，這樣建立索引都是以ref.fa為前綴

構(gòu)建出參考基因組的 FM-index荣堰，建立好參考基因組之后床未，就可以進(jìn)行比對了。不同的比對算法持隧，命令不同即硼。

aln/samse/sampe ----> BWA-backtrack (samse 中的 se 是 single-end 的簡寫，而 sampe 中的 pe 是 paired-end 的簡寫）屡拨。
bwasw ----> BWA-SW
mem ----> BWA-MEM

１. BWA-MEM 算法

該算法先使用 MEM(maximal exact matches) 進(jìn)行 seeding alignments只酥，再使用 SW(affine-gap Smith-Waterman) 算法進(jìn)行 seeds 的延伸。BWA–MEM 算法執(zhí)行局部比對和剪接性呀狼×言剩可能會出現(xiàn) query 序列的多個不同的部位出現(xiàn)各自的最優(yōu)匹配，導(dǎo)致 reads 有多個最佳匹配位點(diǎn)哥艇。有些軟件如 Picard’s markDuplicates 跟 mem 的這種剪接性比對不兼容,在這種情況下绝编，可以使用 –M 選項(xiàng)來將 shorter split hits 標(biāo)記為次優(yōu)。

對應(yīng)的子命令為mem, 基本用法如下

Usage: bwa mem [options] ref.fa reads.fq [mates.fq]

Options:
　　-t INT 線程數(shù)，默認(rèn)是1十饥，增加線程數(shù)窟勃，會減少運(yùn)行時間。
　　-M 將 shorter split hits 標(biāo)記為次優(yōu)逗堵，以兼容 Picard’s markDuplicates 軟件秉氧。
　　-p 若無此參數(shù)：輸入文件只有1個，則進(jìn)行單端比對蜒秤；若輸入文件有2個汁咏，則作為paired reads進(jìn)行比對。若加入此參數(shù)：則僅以第1個文件作為輸入(會忽略第二個輸入序列文件作媚，把第一個文件當(dāng)做單端測序的數(shù)據(jù)進(jìn)行比對)攘滩，該文件必須是read1.fq和read2.fa進(jìn)行reads交叉的數(shù)據(jù)。
　　-R STR 完整的read group的頭部纸泡，可以用 '\t' 作為分隔符漂问，在輸出的SAM文件中被解釋為制表符TAB. read group 的ID，會被添加到輸出文件的每一個read的頭部女揭。
　　-T INT 當(dāng)比對的分值比 INT 小時级解，不輸出該比對結(jié)果，這個參數(shù)只影響輸出的結(jié)果田绑，不影響比對的過程。
　　-a 將所有的比對結(jié)果都輸出抡爹，包括 single-end 和 unpaired paired-end的 reads掩驱，但是這些比對的結(jié)果會被標(biāo)記為次優(yōu)。
　　-Y 對數(shù)據(jù)進(jìn)行soft clipping, 當(dāng)錯配或者gap數(shù)過多比對不上時冬竟，會對序列進(jìn)行切除欧穴，這里的切除并只是在比對時去掉這部分序列，最終輸出結(jié)果中序列還是存在的泵殴，所以稱為soft clipping涮帘。

特別說明：

如果 mates.fq 缺省，且參數(shù) –p 未設(shè)定笑诅，那么 reads.fq 被認(rèn)為是 single-end;
如果 mates.fq 存在调缨，且參數(shù) –p 未設(shè)定，那么 mem 命令會認(rèn)為 read.fq 和 mates.fq 中的 i-th reads 組成一個read對 (a read pair)吆你，這個模式是常用的 paired-end mode弦叶。
如果參數(shù) –p 被設(shè)定，那么妇多， mem 命令會認(rèn)為 read.fq 中的第 2i-th 和第 (2i + 1)-th 的 reads 組成一個 read 對（a read pair）伤哺，這種方式也被成為交錯式的（interleaved paired-end)。在這種情況下，即使有 mates.fq立莉，也會被忽略绢彤。

# SE （single end）
$ bwa mem ref.fa reads.fq > mem-se.sam

# PE（paired end）
$ bwa mem ref.fa read1.fq read2.fq > mem-pe.sam
$ bwa mem -t 4 -M -R "\@RG\tID:{library}\tLB:{library}\tPL:Illumina\tPU:{sample}\tSM:{sample}\" ref.fa read1.fastq read2.fastq > mem-pe.sam 2> ./mem-pe.log

對于超長讀長的reads,比如PacBio和Nanopore測序儀產(chǎn)生的reads, 用法如下

$ bwa mem -x pacbio ref.fa reads.fq > aln.sam
$ bwa mem -x ont2d ref.fa reads.fq > aln.sam

上述的代碼中, ref.fa指的是參考基因組索引的名字，對于上面提供的小鼠的例子而言蜓耻，參考基因組索引的名字為mm10.fasta, 注意是不包含后綴的茫舶。

２.BWA-backtrack 算法

對應(yīng)的子命令為aln/samse/sample

Usage: bwa aln [options] ref.fa read.fq > aln_sa.sai

Options:
　　-o int：允許出現(xiàn)的最大gap數(shù)。
　　-e int：每個gap允許的最大長度媒熊。
　　-d int：不允許在3’端出現(xiàn)大于多少bp的deletion奇适。
　　-i int：不允許在reads兩端出現(xiàn)大于多少bp的indel。
　　-l int：Read前多少個堿基作為seed芦鳍，如果設(shè)置的seed大于read長度嚷往，將無法繼續(xù)，最好設(shè)置在25-35柠衅，與-k 2 配合使用皮仁。
　　-k int：在seed中的最大編輯距離，使用默認(rèn)2菲宴，與-l配合使用贷祈。
　　-t int：要使用的線程數(shù)。
　　-R int：此參數(shù)只應(yīng)用于pair end中喝峦，當(dāng)沒有出現(xiàn)大于此值的最佳比對結(jié)果時势誊，將會降低標(biāo)準(zhǔn)再次進(jìn)行比對。增加這個值可以提高配對比對的準(zhǔn)確率谣蠢，但是同時會消耗更長的時間粟耻，默認(rèn)是32。
　　-I int：表示輸入的文件格式為Illumina 1.3+數(shù)據(jù)格式眉踱。
　　-B int：設(shè)置標(biāo)記序列挤忙。從5’端開始多少個堿基作為標(biāo)記序列，當(dāng)-B為正值時谈喳，在比對之前會將每個read的標(biāo)記序列剪切册烈，并將此標(biāo)記序列表示在BC SAM 標(biāo)簽里，對于pair end數(shù)據(jù)婿禽，兩端的標(biāo)記序列會被連接赏僧。
　　-b ：指定輸入格式為bam格式。

用法如下：
先使用 aln 命令將單獨(dú)的 reads 比對到參考序列扭倾，再使用 samse 或 sampe 生成 sam 文件次哈。常用例子：

# 對于single-read
$ bwa aln [options] ref.fa read.fq > aln_sa.sai                　#尋找 SA coordinates
$ bwa samse [options] ref.fa aln_sa.sai read.fq > aln-se.sam    # 轉(zhuǎn)換SA coordinates輸出為sam

# 對于pair-reads：兩個文件分別處理
$ bwa aln [options] ref.fa read1.fq > aln_sa1.sai
$ bwa aln [options] ref.fa read2.fq > aln_sa2.sai
$ bwa sampe [options] ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln-pe.sam

如果希望多線程運(yùn)行，在其中加入 -t這個參數(shù)吆录，另外-f這個參數(shù)可以指定結(jié)果輸出文件窑滞，如:

$ bwa aln -c -t 3 -f aln_sa1.sai ref.fa  read1.fq

３. BWA-SW 算法

對應(yīng)的子命令為bwasw, 基本用法如下

Usage: bwa bwasw [ options ] ref.fasta reads.fq [mate.fq] > aln.sam

Options:
　　-t INT 使用的線程數(shù)

例子：

$ bwa bwasw genome long_read.fq > aln.sam
$ bwa bwasw genome read1.fq read2.fq > aln-pe.sam

對輸入的第1個文件的所有序列進(jìn)行比對。如果輸如有 2 個文件，則進(jìn)行 paired-end 比對哀卫，此模式僅對 Illumina 的 short-insert 數(shù)據(jù)進(jìn)行比對巨坊。在 Paired-end 模式下，BWA-SW依然輸出剪接性比對結(jié)果此改，但是這些結(jié)果會標(biāo)記為 not properly paired; 同時如果有多個匹配位點(diǎn)趾撵，則不會寫入 mate 的匹配位置。

利用bwa進(jìn)行比對的例子：

1）對參考基因組（reference）構(gòu)建索引

echo "bwa index starts@"`date` && \
cd ref && \
bwa index -a bwtsw genome.fa && \
echo "bwa index ends@"`date`

生成的文件有以下幾種類型：

genome.fasta.ann
genome.fasta.pac
genome.fasta.amb
genome.fasta.bwt
genome.fasta.sa

2）用bwa mem 進(jìn)行比對

echo "bwa mem starts@"`date` && \
cd ref && \
bwa mem -t 4 -M -R "\@RG\tID:{library}\tLB:{library}\tPL:Illumina\tPU:{sample}\tSM:{sample}\" \
                /home/data/genome.fasta \
                R_1.fastq R_2.fastq \
                > bwa.mem.sam \
                2> ./bwa.log&& \
echo "bwa mem ends@"`date`

fai：是對ref基因組文件建的索引共啃，方便軟件快速隨機(jī)讀取基因組序列
sai：是將fastq比對后出來的文件占调，用于最后輸出比對結(jié)果sam文件的

REF:
http://bio-bwa.sourceforge.net/bwa.shtml
https://www.cnblogs.com/emanlee/p/4316587.html
https://blog.csdn.net/u013553061/article/details/53120973
http://www.reibang.com/p/3b86615d647b