bwa 0.7.17-r1188

下载安装

地址：https://sourceforge.net/projects/bio-bwa/files/

wget http://superb-sea2.dl.sourceforge.net/project/bio-bwa/bwa-0.7.17.tar.bz2
tar -jxf bwa-0.7.17.tar.bz2
cd bwa-0.7.17/
make
make install
echo 'export PATH=/home/li.han/Softwares/bwa-0.7.17:$PATH ' >> ~/.bashrc

运行

./bwa

作用

BWA 是一种能够将差异度较小的序列比对到一个较大的参考基因组上的软件包

算法

BWA-backtrack

比对小于100bp的序列

参数：aln / samse / sampe

BWA-SW

比对70bp-1Mbp的序列

参数：bwasw

BWA-MEM

比对70bp-1Mbp的序列，但BWA-MEM更快更准确，在比对70-100bp的reads时也比BWA-backtrack好

参数：mem

使用：

需要两种输入文件：

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

步骤：

1. 使用BWA构建参考基因组的index数据库
2. BWA-MEM / BWA-SW / BWA-backtrack比对

格式：

fai是对ref基因组文件建的索引，方便软件快速随机读取基因组序列
sai是将fastq比对后出来的文件，用于最后输出比对结果sam文件的

step 1: 对参考基因组构建索引

根据reference genome data 建立 Index File，最后生成文件：hg19.fa.amb、hg19.fa.ann、hg19.fa.bwt、hg19.fa.pac和hg19.fa.sa。

bwa index -a bwtsw reference.fa

-a参数：

bwa构建索引有三种算法，三种算法都是基于BWT的，这三种算法通过参数-a is 、－a div和-a bwtsw进行选择。

-a is（效果和-a div是一样的）是默认参数，虽然相对较快，但是需要较大的内存，适用于小于等于2G的参考基因组文件。

-a bwtsw适用于大于等于10Mb的参考基因组文件。

step 2: 与参考基因组比对

使用BWA-backtrack算法：

1. 寻找输入reads文件的SA坐标

aln功能：

将目标序列与参考序列进行比对，生成sai文件

aln用法：bwa aln [options] <prefix> <in.fq>

处理pair-end 数据（read_1.fq和read_2.fq）

bwa aln -l 30 -k 2 -t 4 -I hg19.fa read_1.fq.gz > read_1.sai

bwa aln -l 30 -k 2 -t 4 -I hg19.fa read_2.fq.gz > read_2.sai

处理single 数据（read.fq）

bwa aln -l 30 -k 2 -t 4 -I hg19.fa read.fq.gz > read.sai

如果希望多线程运行，在其中加入 -t这个参数，另外-f这个参数可以指定结果输出文件，如：

bwa aln -t 3 -f read.sai hg19.fa read.fq

aln参数：

-b	输入格式为bam格式
-B int	设置标记序列。从5’端开始多少个碱基作为标记序列，当-B为正值时，在比对之前会将每个read的标记序列剪切，并将此标记序列表示在BC SAM 标签里，对于pair end数据，两端的标记序列会被连接。
-d int	不允许在3’端出现大于多少bp的deletion
-e int	每个gap允许的最大长度
-f	输出文件
-i int	不允许在reads两端出现大于多少bp的indel
-I int	大写 i，表示输入的文件格式为Illumina 1.3+的fq格式。
-k int	在seed中的最大编辑距离，使用默认2，与-l配合使用。
-l int	小写 L，Read前多少个碱基作为seed，如果设置的seed大于read长度，将无法继续，最好设置在25-35，与-k 2 配合使用。
-n	允许最大错配数：mismatch+gap
-o int	允许出现的最大gap数
-R int	此参数只应用于pair end中，当没有出现大于此值的最佳比对结果时，将会降低标准再次进行比对。增加这个值可以提高配对比对的准确率，但是同时会消耗更长的时间，默认是32。
-t	使用的线程数

2: 生成sam格式的比对文件。如果一条read比对到多个位置，会随机选择一种。

如果是pair-end数据

bwa sampe -a 500 -f pair-end.sam hg19.fa read_1.sai read_2.sai read_1.fq read_2.fq

参数说明：

-a int	最大插入片段大小
-f	输出Sam文件
-n int	如果reads比对次数超过多少次，就不在XA标签显示
-o int	pair end两reads中其中之一所允许配对的最大次数，超过该次数，将被视为single end。降低这个参数，可以加快运算速度，对于少于30bp的read，建议降低-o值。
-r str	定义头文件。‘@RG\tID:foo\tSM:bar’，如果在此步骤不进行头文件定义，在后续GATK分析中还是需要重新增加头文件。

如果是single reads数据

bwa samse -f single.sam hg19.fa read.sai read.fq

参数说明：

-n int	如果reads比对次数超过多少次，就不在XA标签显示
-r str	定义头文件。同pair end。

使用bwa-mem算法：

原理：

mem 使用MEMs(maximal exact matches)进行seedling alignments，再使用 SW(affine-gap Smith-Waterman 算法)进行seeding extending。
mem 进行局部比对，因此，对于一条序列的不同区域可能会产生多种最优匹配结果， 这对于long reads 来说尤为重要。有些软件如 Picard’s markDuplicates 跟 mem 的这种剪接性比对不兼容，在这种情况下，可以使用 –M 选项来将 shorter split hits 标记为次优。

用法：

bwa mem -M -t 15 -R '@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name' <ref.fa> <read1.fq> <read2.fq>  >  lane.sam

参数：

-M	处理同一个reads比对到参考基因组上不同位置的情况，将 shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 软件
-p	若无此参数，输入文件只有1个，则进行单端比对；若输入文件有两个，则作为paired reads进行比对。若加入此参数，则仅以第一个文件作为输入，该文件必须是read1.fq和read2.fq进行reads交叉的数据
-R	设置Read Group信息，“\t”分割。 它是用来将比对的read进行分组的，这个信息对于我们后续对比对数据进行错误率分析和Mark duplicate时非常重要。不设置-R参数不会报错，但使用GATK时是会报错的。ID，PL和SM信息在正式的项目中是不能缺少的(如果样本包含多个测序文库的话，LB信息也不要省略)
-t	线程数，默认是1

在Read Group中，有如下几个信息非常重要：

• ID，这是Read Group的分组ID，一般设置为测序的lane ID（不同lane之间的测序过程认为是独立的），下机数据中我们都能看到这个信息的，一般都是包含在fastq的文件名中；
• PL，指的是所用的测序平台，在GATK中，PL只允许被设置为：ILLUMINA，SLX，SOLEXA，SOLID，454，LS454，COMPLETE，PACBIO，IONTORRENT，CAPILLARY，HELICOS或UNKNOWN这12个中的一个，名字方面不区分大小写。如果你在分析的时候这里没设置正确，那么在后续使用GATK过程中可能会碰到类似如下的错误：

ERROR MESSAGE: The platform (xx) associated with read group GATKSAMReadGroupRecord @RG:xx is not a recognized platform.

这个时候你需要对比对文件的header信息进行重写，就会稍微比较麻烦。

• SM，样本ID，同样非常重要，有时候我们测序的数据比较多的时候，那么可能会分成多个不同的lane分布测出来，这个时候SM名字就是可以用于区分这些样本；
• LB，测序文库的名字，这个重要性稍微低一些，主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的fq文件名中找到，如果上面的lane ID足够用于区分的话，也可以不用设置LB。