高通量短讀比對(duì)工具

在過去的十幾年里，隨著高通量測(cè)序(HTS)成本降低，出現(xiàn)了各種測(cè)序概念, DNA-Seq, ChIP-Seq, RNA-Seq, BS-Seq覆蓋了研究領(lǐng)域的方方面面。隨之而來的問題是，如何把這些短片段快速且準(zhǔn)確地回貼到參考基因組上。

解決這個(gè)問題不能直接使用傳統(tǒng)的比對(duì)工具，比如說BLAST，因?yàn)樗鼈兊娜蝿?wù)是找到最多的聯(lián)配，而短序列比對(duì)工具則是要快速?gòu)谋姸酀撛诳蛇x聯(lián)配中找到最優(yōu)的位置。也就是說BLAST和短讀工具的目標(biāo)其實(shí)不太一樣。

在將海量的reads回貼到參考基因組上的過程，大量短讀比對(duì)工具就需要面對(duì)準(zhǔn)確度(accuracy)和精確度(precision)的平衡，也就是盡可能保證每一次的分析結(jié)果是相近的，并且也是符合真實(shí)情況。

mapping and alignment

對(duì)于alignment和mapping，其實(shí)我對(duì)他們之前的區(qū)別一直都不太清楚，并且也不知道它們到底該如何翻譯，總感覺這兩個(gè)詞說的是同一件事情。這里看下Heng Li是如何進(jìn)行定義

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è)重于把序列放到正確的位置，而不管這個(gè)序列的一致性，而聯(lián)配則是主要讓序列和參考序列盡可能的配對(duì)，而不管位置。目前來看，大多數(shù)工具都是想既能找到正確的位置，也保證有足夠多的聯(lián)配，不過明白這兩者的區(qū)別對(duì)于不同項(xiàng)目的分析非常重要。比如說變異檢測(cè)就要優(yōu)先保證聯(lián)配，而RNA-Seq則要盡可能保證把reads放到正確的位置。

如何挑選合適的短讀比對(duì)工具

2012年 Bioinformatics 有一篇文章^[Tools for mapping high-throughput sequencing data ]綜述了目前高通量數(shù)據(jù)的比對(duì)軟件，并且建立主頁https://www.ebi.ac.uk/~nf/hts_mappers/羅列并追蹤目前的比對(duì)軟件。

比對(duì)工具年譜

盡管看起來有那么多軟件，但是實(shí)際使用就那么幾種，BWA(傲視群雄), TopHat(盡管官方都建議用HISAT2，還是那么堅(jiān)挺), SOAP(架不住華大業(yè)務(wù)多)。 由于這些工具都挺成熟，所以選擇軟件更多靠的是信仰，比如說Broad Institute的科學(xué)家喜歡bwa（畢竟是自家的），華大（BGI）喜歡用novoalign(也是自家出品)，只不過novoalign是商業(yè)工具，不買就只能用單核，因此限制了它的傳播。

除了信仰之外，我們挑選短序列比對(duì)工具的時(shí)候還要看什么呢？

• 聯(lián)配算法： 全局，局部還是半全局
• 需要報(bào)道非線性重排(non-linear arrangements)嘛
• 比對(duì)工具如何處理InDels
• 比對(duì)工具支持可變剪切嘛
• 比對(duì)工具能夠過濾出符合需要的聯(lián)配嘛
• 比對(duì)工具能找到嵌合聯(lián)配(chimeric alignments)嘛

最后我們的選擇就落到兩個(gè)工具：BWA和Bowtie2.

BWA和Bowtie的使用簡(jiǎn)介

大部分比對(duì)工具的使用都可以分為兩步，建立索引和比對(duì)索引。值得注意的是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

由于基因組特別小，所以建立索引的速度也會(huì)特別快。

REF=~/biostar/refs/ebola/1976.fa
# bwa
bwa index $REF
# bowtie2
bowtie2-build $REF $REF

索引一覽

序列比對(duì)

為了比對(duì)序列，首先得準(zhǔn)備數(shù)據(jù)文件，可以從SRA上下載項(xiàng)目的Ebola項(xiàng)目的所有runs, 選擇其中一個(gè)作為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

比對(duì)其實(shí)很簡(jiǎn)單,如果只用默認(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é)果是個(gè)SAM文件，那什么是SAM呢，后面繼續(xù)討論。

bwa和bowtie2到底選誰

比較不同的比對(duì)軟件是一個(gè)比較麻煩的事情。最常見的比較方法是，先模擬出一些序列，然后檢查默認(rèn)參數(shù)下的比對(duì)率和運(yùn)行速度

• 10w條read，錯(cuò)誤率為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，錯(cuò)誤率為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的運(yùn)行結(jié)果真的是差異巨大，尤其是10%的錯(cuò)誤率下，幾乎沒有啥能夠比對(duì)上了，讓我們不禁懷疑bowtie2這個(gè)軟件是不是不好使。

讓我們換其他參數(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ù)后比對(duì)率有著明顯的提高，但是-D 2O -R 3 -N 1 -L 20如何得來呢？

也就是說bwa的默認(rèn)參數(shù)是經(jīng)過很好的優(yōu)化來保證在默認(rèn)參數(shù)下的結(jié)果，是不是我們都要選擇bwa呢？也不能如此絕對(duì)，畢竟bowtie2的SAM結(jié)果保留了更多的信息。

最后說一句，選擇比對(duì)軟件在初學(xué)者時(shí)期真的是全靠信仰。