BWA用法笔记

BWA下载与安装

$ wget https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.17.tar.bz2 
$ tar xvfj bwa-0.7.17.tar.bz2 
$ cd bwa-0.7.17
$ make
$./bwa

BWA简单比对

# 建立索引
# 建立索引可以通过子命令调用两种算法
# -is IS线性时间算法，无法处理2GB以上文件
# -Bwtsw BWT-SW中使用的算法，对于短的参考序列不工作，必须大于10Mb
$ bwa index ref.fa

# aln 算法
# 建立 .sai 文件 寻找 SA corrdinates
# pair-end 两个文件分别处理；single-end 直接处理 -t 多线程
$ bwa aln ref.fa read.fq > aln_sa.sai
# 单末端测序结果比对 
$ bwa samse ref.fa aln_sa.sai read.fq > aln_se.sam
# 双末端测序结果比対
$ bwa sampe ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln_pe.sam

# mem 算法
$ bwa mem ref.fa reads.fq > mem_se.sam
$ bwa mem ref.fa read1.fq read2.fq > mem_pe.sam

# 建立索引
-o int：允许出现的最大gap数。

-e int：每个gap允许的最大长度。

-d int：不允许在3’端出现大于多少bp的deletion。

-i int：不允许在reads两端出现大于多少bp的indel。

-l int：Read前多少个碱基作为seed，如果设置的seed大于read长度，将无法继续，最好设置在25-35，与-k 2 配合使用。

-k int：在seed中的最大编辑距离，使用默认2，与-l配合使用。

-t int：要使用的线程数。

-R int：此参数只应用于pair end中，当没有出现大于此值的最佳比对结果时，将会降低标准再次进行比对。增加这个值可以提高配对比对的准确率，但是同时会消耗更长的时间，默认是32。

-I int：表示输入的文件格式为Illumina 1.3+数据格式。

-B int：设置标记序列。从5’端开始多少个碱基作为标记序列，当-B为正值时，在比对之前会将每个read的标记序列剪切，并将此标记序列表示在BC SAM 标签里，对于pair end数据，两端的标记序列会被连接。

-b ：指定输入格式为bam格式。

 

0.BWA简介
bwa（全称Burrows-Wheeler Aligner），  BWA是一款基于BWT的快速比对工具，  项目主页：http://bio-bwa.sourceforge.net/，主要功能是将差异度较小的序列比对到一个较大的参考基因组上。其中提供了三种算法

算法	应用场景
BWA-backtrack	illumina测序结果（reads长度不超过100bp）
BWA-SW	支持序列长度70bp-1Mbp，同时支持剪接性比对(split alignments)
BWA-MEM	最常用，最新，最准确，支持序列长度70bp-1Mbp，表现比BWA-backtrack好

1.构建index
官方文档说明：

bwa index [-p prefix] [-a algoType] <in.db.fasta>

其中需要的数据库序列为fasta格式。

选项介绍：

-p 输出数据库的前缀（默认和输入的文件名一致，输出的数据库在其输入文件所在的文件夹，并以该文件名为前缀）

-a 表示构造BWT索引的算法。可用选项有:

is 是默认的算法，虽然相对较快，但是需要较大的内存，当构建的数据库大于2GB的时候就不能正常工作了。

bwtsw 对于短的参考序列式不工作的，必须要大于等于10MB, 但能用于较大的基因组数据，比如人的全基因组。

2.数据比对
数据比对大多使用MEM(maximal exact matches) 进行 seeding alignments，再使用 SW(affine-gap Smith-Waterman) 算法进行 seeds 的延伸。

BWA–MEM 算法执行局部比对和剪接性。可能会出现 query 序列的多个不同的部位出现各自的最优匹配，导致 reads 有多个最佳匹配位点。这对 long reads 的比对是比较重要的结果。但是却会和 Picard 的 markDuplicates 程序不兼容。

使用方法：

bwa mem [potions] ref.fa reads.fq [mates.fq]

#mem 进行局部比对，因此，对于一条序列的不同区域可能会产生多种最优匹配结果， 这对于long reads 来说尤为重要。 有些软件如 Picard’s markDuplicates 跟 mem 的这种剪接性比对不兼容,在这种情况下，可以使用 –M 选项来将 shorter split hits 标记为次优。

常用选项：

-t  线程数，默认为1。
-M  将 shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 软件。
-p  若无此参数：输入文件只有1个，进行单端比对；若输入文件有2个，则作为paired reads进行比对。若加入此参数：则仅以第1个文件作为输入(输入的文件若有2个，则忽略之)，该文件必须是read1.fq和read2.fa进行reads交叉的数据。
-R STR  完整的read group的头部，可以用 '\t' 作为分隔符， 在输出的SAM文件中被解释为制表符TAB. read group 的ID，会被添加到输出文件的每一个read的头部。
-T INT  当比对的分值比 INT 小时，不输出该比对结果，这个参数只影响输出的结果，不影响比对的过程。
-a      将所有的比对结果都输出，包括 single-end 和 unpaired paired-end的 reads，但是这些比对的结果会被标记为次优。

运行例子：

 bwa mem ref.fa reads.fq > mem-se.sam
 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

参考资料：

1. BWA项目主页
2. BWA命令详解-寂寞先生