生物信息百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进行安装。

conda install -y bwa

六、软件使用：

直接敲bwa，弹出软件选项参数。比对还是分成两大步，建立索引，也叫做建库，然后是比对。首先，利用bwa index可以建立索引，输入参考序列的fasta格式文件，

-a 指定建立索引的算法，bwtsw，is或者rb2，以前没有rb2，而是div。is还是适合小于2G的基因组，bwtwt是大于10M的基因组才行，这个地方不要选错了，细菌基因组一般都小于10M，只能用is算法，人基因组是3G，只能用bwtsw，处于中间的选择哪个都可以。

-p 选项是输出前缀，一般不用设置，就是在参考序列名字后面多出很多文件。

bwa index human_g1k_v37.fasta

这样就开始建立索引，根据文件大小和机器配置不同，需要的时间也不同，对于人的基因组序列，可能需要几个小时，因此很多网站在提供参考序列的时候，连索引文件也一起提供，下载索引比自己建立索引文件要快很多。所以建立完成之后就可以进行比对了，一般是两条pairend的reads。当然其他测序平台的单端reads也可以。现在使用bwa-mem算法只需要两步就行了。bwa mem有很多选项参数，这里我们介绍几个比较重要的。

-R 设置reads标头，放到一对引号中，也就是sam文件中的RG部分，为什么要设置RG表头呢，因为同一样品可能包括多个测序结果，来自不同lane，不同文库，或者不同样品的比对结果合并到同一个文件中进行处理，就需要通过RG进行标记区分。RG每个标记用冒号分割键和值，不同标记用 '\t' 分隔。例如'@RG\tID:foo\tSM:bar\tLB: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'。

bwa bwasw genome long_read.fq > all.sam
bwa bwasw genome read1.fq read2.fq > all.sam

七、使用案例：
现有参考序列ref.fna，illumina测序两条reads，reads.1.fq.gz与reads.2.fq.gz。
bwa index ref.fna

bwa mem -R '@RG\tID:E1602061\tSM:bar\tLB:library1' -t 4 ref.fna 
   reads.1.fq.gz reads.2.fq.gz  >result.sam

八、注意事项：

1、安装时需要依赖zlib；
2、建库时注意选择正确的算法；
3、最新版本bwa，需要使用-R选项。

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

（更多精彩内容，欢迎关注微信公众号）

您可能还会感兴趣的