WGS分析笔记（2）—— bwa vs bowtie2

  在进行正式的mapping记录之前，先记录一下bwa与bowtie2在mapping这个环节的情况。
  一般对于WGS结果的mapping，一般推荐的软件有两款，分别是bwa和bowtie2，大多数的公司报告或者网上的教程，我所看到的都是使用bwa进行比对的，这里，我来进行一下两个软件的对比。
  实验对象还是之前文章提到的那个样本的数据，我只取用其中的一对数据进行mapping并比较。
  比较之前先进行一下软件安装、参考序列下载并建立索引

bowtie2
    $ wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip
    #https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4/bowtie2-2.3.4-linux-x86_64.zip
    $ unzip bowtie2-2.2.9-linux-x86_64.zip
    $ ln -sf /biosoft/bowtie2-2.3.4.3-linux-x86_64/bowtie2 /home/shiyuantong/bin/bowtie2
BWA:
    $ wget https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.17.tar.bz2
    $ tar -jxvf bwa-0.7.17.tar.bz2 # x extracts, v is verbose (details of what it is doing), f skips prompting for each individual file, and j tells it to unzip .bz2 files
    $ make

  关于安装这里有一些需要注意的地方！！！！！
  首先是bowtie2，建议大家使用2.3.4的下载链接，我在下载的时候最新版是2.3.4.3，但是在使用的出现了报错！！！！！
  报错的内容如下（当时没截图）：
    Segmentation fault (core dumped) (ERR): bowtie2-align exited with value 139
  这个报错只会出现在批量处理的脚本中，对单个样本的处理并没有影响，但是实际使用的时候，大家都是批量处理样本，怎么可能一个样本一个命令，因此推荐2.3.4的版本，当然，下面的比较不会涉及这个问题。
  还有就是BWA了，这个软件ubuntu用户也可以直接使用sudo apt-get install bwa命令进行安装，我看了一下，两种方法的版本是一致的，都是0.7.17。
  然后是参考序列，这里直接使用ucsc的hg19，下载与建立索引方式如下，根据自己的需要调整目录

hg19：
    $ cd /your/path/of/reference/
    $ wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
    $ tar zvfx chromFa.tar.gz
    $ cat *.fa > hg19.fa
    $ rm chr*.fa
建立bwa索引：
    $ bwa index -a bwtsw  hg19.fa
    # 产生.bwt .pac .ann .amb .sa五个新文件
    # -a：两种构建index算法，bwtsw以及is，bwtsw适用大于10MB的参考基因组，比如人，is适用于小于2GB的数据库，是默认的算法，速度较快，需要较大的内存，
    # -p：输出数据库的前缀，默认与输入文件名一致，这里我没有加这个参数，直接输出到当前目录
建立bowtie2索引：
    $ bowtie2-build hg19.fa hg19.fa
    #  bowtie2-build命令在安装bowtie2的目录下找到
    # 第一个hg19.fa代表输入的参考序列
    # 第二个hg19.fa代表输出的索引文件前缀
    #产生六个.bt2新文件

  上述程序建立索引速度较慢，尤其是bowtie2，但是一次建好，一劳永逸，请大家耐心等待，也可以放在后台，防止终端突然的中断。
  建立好索引就可以直接开始比对了，以下是我的比对程序，都开了24线程，用nohup …… &放在后台运行，用time记录时间。

$nohup time bowtie2 -p 24 -x /your/path/of/reference/ucsc.hg19.fasta --rg-id W2018001 --rg PL:ILLUMINA --rg LB:W2018001 --rg SM:W2018001 -1 W2018001_NZTD180602206_HCV5MDMXX_L1.cleaned.1.fq.gz -2 W2018001_NZTD180602206_HCV5MDMXX_L1.cleaned.2.fq.gz -S W2018001.bowtie2.sam > W2018001.bowtie2.log &
$nohup time bwa mem -t 24 -M -R "@RG\tID:W2018001\tPL:ILLUMINA\tLB:W2018001\tSM:W2018001" /your/path/of/reference/ucsc.hg19.fasta W2018001_NZTD180602206_HCV5MDMXX_L1.cleaned.1.fq.gz W2018001_NZTD180602206_HCV5MDMXX_L1.cleaned.2.fq.gz 1>W2018001.bwa.sam 2>W2018001.bwa.log &

  这一步会比较久，我也是经过漫长的半天等待终于迎来了结果，首先看一下速度吧。先前的脚本使用了time的命令，可以直接看到速度，在日志文件里。

时间对比

  日志文件的最后两行就是time命令输出的结果，所以没有必要用cat查看，而且bwa的日志文件，要是用cat怕是屏幕要炸。图中可以看到两个红色的框，就是我标出来的时间。（其实我原来用time命令，结果不长这样的，这个结果不太利于观看，但是也能说明问题了）
  很明显的可以看到半套全基因组的数据（我只用了样本一半的数据）bowtie2跑的更快一些，但其实大家不用纠结这个点。因为上一次我用24线程，一样的脚本一样的数据，跑bowtie2花了六个多小时，速度没有bwa快，同时以前在使用酵母的测序数据（数据量会比较小）的时候，明显发现bwa速度比bowtie2快，甚至在说法上你也会发现不同的人给你的说法是不一样的，有些人说bwa快有些人说bowtie2快，网上看帖子也没有一个十分明确的说法哪个速度快。这里大家完全可以用自己的数据和脚本跑一下，看看结果。
  接下来我想看看比对效果，这里我先采用了samtools的flagstat分别进行统计，下面是安装samtools的步骤：

samtools:
    $ wget https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
    $ tar xvfj samtools-1.9.tar.bz2
    $ cd samtools-1.9
    $ ./configure --prefix=/where/to/install
    $ make
    $ make install
#samtools其实我到现在为止装的最崩溃的软件之一了，因为在实际安装的时候你会发现它需要各种各样的库的支持，对于使用新机器的我，我基本是安装，报错缺什么库，安装缺的库，重新安装，这么折腾了一下午

  接下来就是使用统计工具，其实很简单。

$ samtools flagstat W2018001.bowtie2.sam >W2018001.bowtie2.flagstat
$ samtools flagstat W2018001.bwa.sam >W2018001.bwa.flagstat

  这个也需要花一点时间，但不会太长，看一下结果。

flagstat
  这个结果还是比较清楚的，bwa的结果比bowtie2稍稍好一点。但相差不是很大，所以对于这两个软件，一直是公说公有理，婆说婆有理，这里我用另一个软件对结果进行统计，再进行对比试试。
  RSeQC是一个功能强大的软件，里面有很多实用的小工具，其中的bam_stat就是一个实用的bam/sam结果统计工具，安装方式也是相当简单了，就是一个python的包，支持python2.x和python3.x，这里我选用python3的pip来安装，因为本人习惯使用python3。

$ pip3 install RSeQC

  使用和结果如下，由于我这个sam文件比较大，运行起来比较慢，所以我开了俩终端。

bam_stat

  其实我最想看的unique mapping的reads，因为后期为了降低假阳性，在处理bam文件的时候会选择unique mapped的reads，但是在查看说明书无果后，找遍论坛没有找到一个能够说服我的筛选unique mapped的方式。
  有这么几个方式，一个说是看tag，但是bwa的结果，你仔细看说明书和结果，会发现，这个tag并没有什么用，bowtie2倒是还可以。
  第二个也是说的比较多的一个，看MAPQ。那么mapq是啥呢，我来贴几张图。

bowtie2的说明 SAM的说明 官网说明
  分别是sam格式官网的说明，bowtie2官网的说明，这两个说明的公式是一样的，都指向最后一个官网的说明。看到这个官网的公式，我直接就傻掉了，反正到现在也没推出个所以然来。但是前两张图就很好理解了。但是和很多人说的MAPQ>=1就是unique mapping，我觉得是对不上的。对于这点我不多说，这里的解释也是目前为止我最能接受的。
  那上面那个结果，bam_stat，我在阅读源码后，发现是以MAPQ>=30作为阈值来挑选是否unique的。由于bwa和bowtie2的mapq的计算方式不一样，这个结果其实并不可信。于是我写了一个脚本，看了一下mapq的分布情况。
bowtie2 bwa

  这个图能说明什么呢，有待商榷。
  这个时候再回过来看一下bowtie2的输出结果和大家说的bowtie2的筛选unique mapping的方法以及结果。

bowtie2.log
bowtie2_tag

  其实到这里我也不知道该怎么办了，到现在还是不知道，bwa结果用mapq>=1是否能得到unique mapping。这个结果对于后续分析影响有多大我不好说，至于怎么选择，我也不发表意见。
  2019-1-9补充bwa结果，按mapq分布，分别计算>=1,10,20,30的比例。为大家选择MAPQ作为筛选提供一个参考。

MAPQ比例

  但是回归主题，bwa和bowtie2，我决定选择bwa。
  水平有限，要是存在什么错误请评论指出！请大家多多批评指正，相互交流，共同成长，谢谢！！！