高通量短讀比對工具
在過去的十幾年里，隨著高通量測序(HTS)成本降低，出現(xiàn)了各種測序概念, DNA-Seq, ChIP-Seq, RNA-Seq, BS-Seq覆蓋了研究領(lǐng)域的方方面面耕肩。隨之而來的問題是呀闻，如何把這些短片段快速且準(zhǔn)確地回貼到參考基因組上厨相。
解決這個問題不能直接使用傳統(tǒng)的比對工具炼七，比如說BLAST，因為它們的任務(wù)是找到最多的聯(lián)配虱朵，而短序列比對工具則是要快速從眾多潛在可選聯(lián)配中找到最優(yōu)的位置莉炉。也就是說BLAST和短讀工具的目標(biāo)其實不太一樣。
在將海量的reads回貼到參考基因組上的過程碴犬，大量短讀比對工具就需要面對準(zhǔn)確度(accuracy)和精確度(precision)的平衡絮宁，也就是盡可能保證每一次的分析結(jié)果是相近的，并且也是符合真實情況翅敌。
mapping and alignment
對于alignment和mapping羞福，其實我對他們之前的區(qū)別一直都不太清楚，并且也不知道它們到底該如何翻譯蚯涮，總感覺這兩個詞說的是同一件事情治专。這里看下Heng Li是如何進行定義
Mapping(映射)
- A mapping is a region where a read sequence is placed
- A mapping is regarded to be correct if it overlaps the ture region
Alignment(聯(lián)配)
- An alignment is the detailed placement of each base in a read.
- An alignment is regarded to be correct if each base is placed correctly.
也就是說mapping側(cè)重于把序列放到正確的位置，而不管這個序列的一致性遭顶，而聯(lián)配則是主要讓序列和參考序列盡可能的配對张峰，而不管位置。目前來看棒旗，大多數(shù)工具都是想既能找到正確的位置喘批，也保證有足夠多的聯(lián)配撩荣，不過明白這兩者的區(qū)別對于不同項目的分析非常重要。比如說變異檢測就要優(yōu)先保證聯(lián)配饶深，而RNA-Seq則要盡可能保證把reads放到正確的位置餐曹。
如何挑選合適的短讀比對工具
2012年 Bioinformatics 有一篇文章^[Tools for mapping high-throughput sequencing data ]綜述了目前高通量數(shù)據(jù)的比對軟件，并且建立主頁https://www.ebi.ac.uk/~nf/hts_mappers/羅列并追蹤目前的比對軟件敌厘。
盡管看起來有那么多軟件台猴，但是實際使用就那么幾種，BWA(傲視群雄), TopHat(盡管官方都建議用HISAT2俱两，還是那么堅挺), SOAP(架不住華大業(yè)務(wù)多)饱狂。 由于這些工具都挺成熟，所以選擇軟件更多靠的是信仰宪彩，比如說Broad Institute的科學(xué)家喜歡bwa（畢竟是自家的）休讳，華大（BGI）喜歡用novoalign(也是自家出品)，只不過novoalign是商業(yè)工具尿孔，不買就只能用單核俊柔，因此限制了它的傳播。
除了信仰之外纳猫，我們挑選短序列比對工具的時候還要看什么呢婆咸？
- 聯(lián)配算法： 全局竹捉，局部還是半全局
- 需要報道非線性重排(non-linear arrangements)嘛
- 比對工具如何處理InDels
- 比對工具支持可變剪切嘛
- 比對工具能夠過濾出符合需要的聯(lián)配嘛
- 比對工具能找到嵌合聯(lián)配(chimeric alignments)嘛
最后我們的選擇就落到兩個工具：BWA和Bowtie2.
BWA和Bowtie的使用簡介
大部分比對工具的使用都可以分為兩步芜辕，建立索引和比對索引。值得注意的是BWA有兩種算法块差，aln和mem分別處理低于100bp和大于70bp的短讀侵续。bowtie也有1和2兩代，處理50bp以下和50bp以上的短讀憨闰，注意選擇状蜗。
建立索引
需要先用efetch下載ebola參考基因組，如果網(wǎng)絡(luò)不佳鹉动，直接去NCBI查找到下載也可以
mkdir -p ~/biostar/refs/ebola
cd ~/biostar
# efetch下載
efetch -db=nuccore -format=fasta -id=AF086833 > ~/refs/ebola/1976.fa
# wget下載
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/848/505/GCF_000848505.1_ViralProj14703/GCF_000848505.1_ViralProj14703_genomic.fna.gz
由于基因組特別小轧坎，所以建立索引的速度也會特別快。
REF=~/biostar/refs/ebola/1976.fa
# bwa
bwa index $REF
# bowtie2
bowtie2-build $REF $REF
序列比對
為了比對序列泽示，首先得準(zhǔn)備數(shù)據(jù)文件缸血，可以從SRA上下載項目的Ebola項目的所有runs, 選擇其中一個作為demo數(shù)據(jù)。
esearch -db sra -query PRJNA257197 | efetch -format runinfo > runinfo.csv
mkdir raw_data
cd raw_data
fastq-dump -X 10000 --split-files SRR1972739
比對其實很簡單,如果只用默認(rèn)參數(shù)的話
R1=raw_data/SRR1972739_1.fastq
R2=raw_data/SRR1972739_2.fastq
# bwa-mem
bwa mem $REF $R1 $R2 > bwa_mem_out.sam
# bowite2
bowtie2 -x $REF -1 $R1 -2 $R2 > bowtie2_out.sam
結(jié)果是個SAM文件械筛，那什么是SAM呢捎泻，后面繼續(xù)討論。
bwa和bowtie2到底選誰
比較不同的比對軟件是一個比較麻煩的事情埋哟。最常見的比較方法是笆豁，先模擬出一些序列，然后檢查默認(rèn)參數(shù)下的比對率和運行速度
- 10w條read，錯誤率為1%闯狱，默認(rèn)參數(shù)
# dwgsim的安裝方法見biostar handbook
~/bin/dwgsim -N 100000 -e 0.01 -E 0.01 $REF data
R1=data.bwa.read1.fastq.gz
R2=data.bwa.read2.fastq.gz
time bwa mem $REF $R1 $R2 > bwa.sam
# 4s 95.04%
time bowtie2 -x $REF -1 $R1 -2 $R2 > bowtie2.sam
# 10秒煞赢， 94.82%
- 10w條read，錯誤率為10%哄孤，默認(rèn)參數(shù)
~/bin/dwgsim -N 100000 -e 0.1 -E 0.1 $REF data
R1=data.bwa.read1.fastq.gz
R2=data.bwa.read2.fastq.gz
time bwa mem $REF $R1 $R2 > bwa.sam
samtools flagstat bwa.sam
# 7s 83.16%
time bowtie2 -x $REF -1 $R1 -2 $R2 > bowtie2.sam
# 4秒耕驰，29.01%
在默認(rèn)參數(shù)下，bowtie2的運行結(jié)果真的是差異巨大录豺，尤其是10%的錯誤率下朦肘，幾乎沒有啥能夠比對上了，讓我們不禁懷疑bowtie2這個軟件是不是不好使双饥。
讓我們換其他參數(shù)試試看
bowtie2 --very-sensitive-local -x $REF -1 $R1 -2 $R2 > bowtie.sam
# 10s, 63.21%
time bowtie2 -D 20 -R 3 -N 1 -L 20 -x $REF -1 $R1 -2 $R2 > bowtie.sam
# 11s, 87.11%
bowtie2在我們更換參數(shù)后比對率有著明顯的提高媒抠，但是-D 2O -R 3 -N 1 -L 20如何得來呢？
也就是說bwa的默認(rèn)參數(shù)是經(jīng)過很好的優(yōu)化來保證在默認(rèn)參數(shù)下的結(jié)果咏花，是不是我們都要選擇bwa呢趴生？也不能如此絕對，畢竟bowtie2的SAM結(jié)果保留了更多的信息昏翰。
最后說一句苍匆，選擇比對軟件在初學(xué)者時期真的是全靠信仰。
