- 更新于2020.10.29
- 在進行正式的mapping記錄之前边苹，先記錄一下bwa與bowtie2在mapping這個環(huán)節(jié)的情況。
- 一般對于WGS結(jié)果的mapping杜顺，一般推薦的軟件有兩款财搁，分別是bwa和bowtie2，大多數(shù)的公司報告或者網(wǎng)上的教程躬络，我所看到的都是使用bwa進行比對的尖奔，這里，我來進行一下兩個軟件的對比穷当。
- 實驗對象還是之前文章提到的那個樣本的數(shù)據(jù)提茁，我只取用其中的一對數(shù)據(jù)進行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
- 關(guān)于安裝這里有一些需要注意的地方＼畋狻！；鸬恕５と酢！
- 首先是bowtie2铲咨，建議大家使用2.3.4的下載鏈接躲胳，我在下載的時候最新版是2.3.4.3，但是在使用的出現(xiàn)了報錯Ｏ死铡Ｅ髌弧！Ｒ√臁粹湃！
- 報錯的內(nèi)容如下（當時沒截圖）：
Segmentation fault (core dumped) (ERR): bowtie2-align exited with value 139 - 這個報錯只會出現(xiàn)在批量處理的腳本中，對單個樣本的處理并沒有影響泉坐，但是實際使用的時候为鳄，大家都是批量處理樣本，怎么可能一個樣本一個命令腕让，因此推薦2.3.4的版本孤钦，當然，下面的比較不會涉及這個問題纯丸。
- 還有就是BWA了偏形，這個軟件ubuntu用戶也可以直接使用
sudo apt-get install bwa命令進行安裝，我看了一下觉鼻，兩種方法的版本是一致的俊扭，都是0.7.17。
(注：bwa-mem2 已經(jīng)更新坠陈，可以直接下載編譯好的程序使用萨惑，在結(jié)果一致的前提下速度提升一倍左右，可以參考本文仇矾。) - 然后是參考序列咒钟，這里直接使用ucsc的hg19，下載與建立索引方式如下若未，根據(jù)自己的需要調(diào)整目錄
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
# 產(chǎn)生.bwt .pac .ann .amb .sa五個新文件
# -a：兩種構(gòu)建index算法朱嘴，bwtsw以及is，bwtsw適用大于10MB的參考基因組粗合，比如人萍嬉，is適用于小于2GB的數(shù)據(jù)庫，是默認的算法隙疚，速度較快壤追，需要較大的內(nèi)存，
# -p：輸出數(shù)據(jù)庫的前綴供屉，默認與輸入文件名一致行冰，這里我沒有加這個參數(shù)溺蕉，直接輸出到當前目錄
建立bowtie2索引：
$ bowtie2-build hg19.fa hg19.fa
# bowtie2-build命令在安裝bowtie2的目錄下找到
# 第一個hg19.fa代表輸入的參考序列
# 第二個hg19.fa代表輸出的索引文件前綴
#產(chǎn)生六個.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 &
-
這一步會比較久，我也是經(jīng)過漫長的半天等待終于迎來了結(jié)果梢褐，首先看一下速度吧徘层。先前的腳本使用了
time的命令，可以直接看到速度利职，在日志文件里趣效。
時間對比 日志文件的最后兩行就是
time命令輸出的結(jié)果，所以沒有必要用cat查看猪贪，而且bwa的日志文件跷敬，要是用cat怕是屏幕要炸。圖中可以看到兩個紅色的框热押，就是我標出來的時間西傀。（其實我原來用time命令，結(jié)果不長這樣的桶癣，這個結(jié)果不太利于觀看拥褂，但是也能說明問題了）很明顯的可以看到半套全基因組的數(shù)據(jù)（我只用了樣本一半的數(shù)據(jù)）bowtie2跑的更快一些，但其實大家不用糾結(jié)這個點牙寞。因為上一次我用24線程饺鹃，一樣的腳本一樣的數(shù)據(jù)，跑bowtie2花了六個多小時间雀，速度沒有bwa快悔详，同時以前在使用酵母的測序數(shù)據(jù)（數(shù)據(jù)量會比較小）的時候惹挟，明顯發(fā)現(xiàn)bwa速度比bowtie2快茄螃，甚至在說法上你也會發(fā)現(xiàn)不同的人給你的說法是不一樣的，有些人說bwa快有些人說bowtie2快连锯，網(wǎng)上看帖子也沒有一個十分明確的說法哪個速度快归苍。這里大家完全可以用自己的數(shù)據(jù)和腳本跑一下用狱，看看結(jié)果。
接下來我想看看比對效果拼弃，這里我先采用了samtools的flagstat分別進行統(tǒng)計夏伊，下面是安裝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其實我到現(xiàn)在為止裝的最崩潰的軟件之一了，因為在實際安裝的時候你會發(fā)現(xiàn)它需要各種各樣的庫的支持肴敛，對于使用新機器的我，我基本是安裝吗购，報錯缺什么庫医男，安裝缺的庫，重新安裝捻勉，這么折騰了一下午
- 接下來就是使用統(tǒng)計工具镀梭，其實很簡單。
$ samtools flagstat W2018001.bowtie2.sam >W2018001.bowtie2.flagstat
$ samtools flagstat W2018001.bwa.sam >W2018001.bwa.flagstat
-
這個也需要花一點時間踱启，但不會太長报账，看一下結(jié)果。
flagstat - 這個結(jié)果還是比較清楚的埠偿，bwa的結(jié)果比bowtie2稍稍好一點透罢。但相差不是很大，所以對于這兩個軟件冠蒋，一直是公說公有理羽圃，婆說婆有理，這里我用另一個軟件對結(jié)果進行統(tǒng)計抖剿，再進行對比試試朽寞。
- RSeQC是一個功能強大的軟件，里面有很多實用的小工具斩郎，其中的bam_stat就是一個實用的bam/sam結(jié)果統(tǒng)計工具脑融，安裝方式也是相當簡單了，就是一個python的包缩宜，支持python2.x和python3.x肘迎，這里我選用python3的pip來安裝，因為本人習慣使用python3锻煌。
$ pip3 install RSeQC
-
使用和結(jié)果如下膜宋，由于我這個sam文件比較大，運行起來比較慢炼幔，所以我開了倆終端秋茫。
bam_stat - 其實我最想看的unique mapping的reads，因為后期為了降低假陽性乃秀，在處理bam文件的時候會選擇unique mapped的reads肛著，但是在查看說明書無果后圆兵，找遍論壇沒有找到一個能夠說服我的篩選unique mapped的方式。
- 有這么幾個方式枢贿，一個說是看tag殉农，但是bwa的結(jié)果，你仔細看說明書和結(jié)果局荚，會發(fā)現(xiàn)超凳，這個tag并沒有什么用，bowtie2倒是還可以耀态。
-
第二個也是說的比較多的一個轮傍，看MAPQ。那么mapq是啥呢首装，我來貼幾張圖创夜。
bowtie2的說明
- 分別是sam格式官網(wǎng)的說明，bowtie2官網(wǎng)的說明仙逻，這兩個說明的公式是一樣的驰吓，都指向最后一個官網(wǎng)的說明∠捣睿看到這個官網(wǎng)的公式檬贰，我直接就傻掉了，反正到現(xiàn)在也沒推出個所以然來缺亮。但是前兩張圖就很好理解了偎蘸。但是和很多人說的MAPQ>=1就是unique mapping，我覺得是對不上的瞬内。對于這點我不多說迷雪，這里的解釋也是目前為止我最能接受的。
-
那上面那個結(jié)果虫蝶，bam_stat章咧，我在閱讀源碼后，發(fā)現(xiàn)是以MAPQ>=30作為閾值來挑選是否unique的能真。由于bwa和bowtie2的mapq的計算方式不一樣赁严，這個結(jié)果其實并不可信。于是我寫了一個腳本粉铐，看了一下mapq的分布情況疼约。
bowtie2
這個圖能說明什么呢，有待商榷蝙泼。
-
這個時候再回過來看一下bowtie2的輸出結(jié)果和大家說的bowtie2的篩選unique mapping的方法以及結(jié)果程剥。
bowtie2.log
bowtie2_tag 其實到這里我也不知道該怎么辦了，到現(xiàn)在還是不知道汤踏，bwa結(jié)果用mapq>=1是否能得到unique mapping织鲸。這個結(jié)果對于后續(xù)分析影響有多大我不好說舔腾，至于怎么選擇，我也不發(fā)表意見搂擦。
-
2019-1-9補充bwa結(jié)果稳诚，按mapq分布，分別計算>=1,10,20,30的比例瀑踢。為大家選擇MAPQ作為篩選提供一個參考扳还。
MAPQ比例
- 但是回歸主題，bwa和bowtie2橱夭，我決定選擇bwa氨距。(注：其實從結(jié)果上來看bwa和bowtie2的效果是差不多的，但是真的是不是差不多我覺得從 mapping 結(jié)果上看也是片面徘钥，要看看最終 call variants 的時候不同的 bam 文件是不是會導致不同的效果衔蹲，去比較 call 出來的變異的各項指標（recall肢娘，true positive rate等）呈础，最后，單個樣本上的結(jié)果也是片面的橱健，要從大量的樣本上去獲得統(tǒng)計結(jié)果才能有說服力而钞。之所以選擇 bwa 主要考慮的是在效果差不多的情況下大家的認可和使用程度，目前看來大部分的文獻中都是使用 bwa 作為上游分析工具拘荡，這不能說明 bwa 就是最好的臼节，但是至少能說明它是最容易被大家認可的。)
- 水平有限珊皿，要是存在什么錯誤請評論指出网缝！請大家多多批評指正，相互交流蟋定，共同成長粉臊，謝謝！Ｊ欢怠扼仲！
