• 開始之前特別感謝從美國“人肉背回”這兩個軟件的小伙伴！@Olivia阿儀_鴉雀
  

• 之前講過目前用的比較多的兩個二代測序基因組比對軟件，bwa 和 bowtie2，也比較了一下效果，總的來說兩者的效果相差不大，但是 bwa 可能更加的準(zhǔn)確一些，而很多文章以及一些組學(xué)計劃都是用的 bwa。
• 今天主要來比較一下一些新的 mapping 工具，vg 和 Graph Genome，相對于傳統(tǒng)基于本文進(jìn)行比對的工具（bwa等），這類工具基于 Graph 進(jìn)行比對，賣點是速度更快，更加準(zhǔn)確。論文鏈接如下，文章是開放的，直接可以下載：
  • Variation graph toolkit improves read mapping by representing genetic variation in the reference
  • Fast and accurate genomic analyses using genome graphs
• 文章都是發(fā)表在高影響因子的 nature 子刊上，基本都是踩著 bwa 上位的，那我就想測試一下這倆軟件（在我使用的時候我查了一下百度，發(fā)現(xiàn)國內(nèi)根本沒人用這玩意兒）有沒有文章說的那么好，再加之 bwa 更新了新的 mem2 ，其速度又有大幅度的提升，但是結(jié)果沒有什么改變。
• 首先是測試數(shù)據(jù)，這里我用GIAB（瓶中基因組計劃，具體看鏈接文章，也是開放論文，總的這三篇論文篇幅都很短，感興趣的不妨讀一讀。）中的 ChineseTrio 中的父母的 illumina 測序數(shù)據(jù)，隨機(jī)取了兩個個體的每個批次（共四個批次）中的一組數(shù)據(jù)，就是總共8對 fastq 文件。

數(shù)據(jù)預(yù)處理

• 因為我們直接下載的是 GIAB 中的 fastq 文件，常規(guī)步驟還是需要在 mapping 前進(jìn)行一下質(zhì)控處理的，這里的話可以參考我之前的文章。

$ fastp -i /your/path/of/data/male_data/141008_D00360_0058_AHB675ADXX/NA24694_GCCAAT_L001_R1_001.fastq.gz -I /your/path/of/data/male_data/141008_D00360_0058_AHB675ADXX/NA24694_GCCAAT_L001_R2_001.fastq.gz -o /bio-analysis1/mapper_test/data/N1/NA24694_GCCAAT_L001_R1_001.fastq.gz -O /bio-analysis1/mapper_test/data/N1/NA24694_GCCAAT_L001_R2_001.fastq.gz -h male.008.html -j male.008.json -c -q 20 -u 50 -n 15 -5 20 -3 20 -w 12

• 這里我是以其中一對 fastq 為例貼的命令行，但沒有改變輸出后的名字，只是修改了輸出路徑，也方便后續(xù)對照數(shù)據(jù)的名字。
• 這里不貼結(jié)果了，但是可以說明一下，數(shù)據(jù)的質(zhì)量很高的，即使不進(jìn)行這一步質(zhì)控差別也不大，而且用的是 pcr-free 的建庫方式，所以基本沒有什么duplicate。
• 下面就介紹一下三個軟件的使用方式以及測試結(jié)果。

bwa mem2

• 標(biāo)題鏈接即新版本 bwa mem 的 github 地址，使用非常之簡單。

$ curl -L https://github.com/bwa-mem2/bwa-mem2/releases/download/v2.0pre2/bwa-mem2-2.0pre2_x64-linux.tar.bz2 | tar jxf -
$ bwa-mem2-2.0pre2_x64-linux/bwa-mem2 index ref.fa
$ bwa-mem2-2.0pre2_x64-linux/bwa-mem2 mem ref.fa read1.fq read2.fq > out.sam

• 直接下載解壓就可以用了，而且比對的參數(shù)和原來的基本沒什么變化，當(dāng)然也可以不加什么參數(shù)，直接按照上面的來進(jìn)行。
• 值得注意的一點是，新的 bwa 需要重新建立索引，速度還是比較久的，但是這種建立一次，終身管用的操作，也無需特別計較時間。
• 除此之外就是，新的軟件建立的索引文件和原來的索引文件有所區(qū)別的（部分一樣），而且新的索引文件著實耗費空間，我以 hs37d5 作為參考基因組進(jìn)行實驗（不同參考基因組之間的區(qū)別可以見這里）：
  
  reference 索引文件
• 可以看到以.2bit.64和.8bit.32這倆文件很占空間，(這里漏截一個749M大小的.pac文件)，這里我之所以用 hs37d5 （b37_decoy）作為參考基因組也是和另外兩個軟件有關(guān)，而且就主體序列而言，hs37d5 和 hg19 沒什么區(qū)別，所以不用過度糾結(jié)。

# 兩個下載該參考基因組的地址，一個是 EBI 網(wǎng)站上的，一個是 GATK 上的，沒啥區(qū)別，除了文件名，下載后請自行解壓。
$ wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz
$ wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/human_g1k_v37_decoy.fasta.gz

• 下面我以8對數(shù)據(jù)中的一對為例，把命令行貼一下：

$ time bwa-mem2 mem -t 12 -R "@RG\tID:NA24694\tPL:ILLUMINA\tLB:L001\tSM:NA24694" /your/path/of/b37/human_g1k_v37_decoy.fasta /bio-analysis1/mapper_test/data/N1/NA24694_GCCAAT_L001_R1_001.fastq.gz /bio-analysis1/mapper_test/data/N1/NA24694_GCCAAT_L001_R2_001.fastq.gz | samtools view -Shb - > male008.bam
# -t 是線程數(shù)，三個軟件我都用12線程進(jìn)行運行
# -R 是頭文件，這里可以不設(shè)置
# time命令用來統(tǒng)計時間

vg

• 標(biāo)題即 github 鏈接，軟件提供了打包好的二進(jìn)制文件，也可以不安裝直接用。

$ wget https://github.com/vgteam/vg/releases/download/v1.22.0/vg
$ chmod +x vg

• 值得一提的是，無論是 vg 還是 Graph Genome 都有自己的一套 pipeline ，其中包含了 calling 的功能，尤其是 vgtools，不僅能 call SNVs/INDELs 還有專門的call SVs 工具，這工具也發(fā)表在了高影響因子的 Genome Biology 上，Genotyping structural variants in
  pangenome graphs using the vg toolkit。而 Graph Genome 除了具備 calling 功能（一個命令同時 call snv/indel/sv ）也提供了一套基于trio的分析工具VBT-TrioAnalysis。

• 但是vg真是一款令人火大的軟件，在我測試開始，到寫下這篇文章，都是懷著一種悲憤的心情。

  
  
  vg pipeline
• 這是官方提供的pipeline的上半部分，到 mapping 為止，看著還挺復(fù)雜的，最主要的是，你會發(fā)現(xiàn)第一步有一個 vcf 文件，這是 vg 所需要的先驗文件，以你所要使用的 reference 得到的 vcf，于是就有了一個先有雞還是先有蛋的問題……好了，開玩笑的，但是有一點確實很讓人詬病，軟件也沒說這個 vcf 是任意一個 vcf 還是需要越大的群體得到的結(jié)果越好，畢竟不同個體的變異肯定是相差很大的，如果是群體的 vcf，那么作為一個軟件，像 GATK 一樣給一個參考的數(shù)據(jù)庫豈不是更好，Graph Genome 也需要先驗 vcf 文件，但是 Graph Genome 是提供了的。
• 官方的 README 文件的說明簡單的令人發(fā)指，找了半天終于找到了一個通過公共數(shù)據(jù)庫建立索引文件的文章。
• 文章中每一步寫的還是很詳細(xì)的，包括可能需要的運行時間，臨時空間大小，產(chǎn)生文件大小。但就是因為寫的越細(xì)才越讓人生氣，因為這部分的計算資源還是挺大的，官方既然自己做過一遍，為什么不直接把結(jié)果文件掛出來讓人下載作為使用參考，就像 GATK Bundle 一樣。
• 氣歸氣，接下來還是一步一步演示一下，不過這個流程還是很花費時間的。
• 首先說一下，這里用的公共數(shù)據(jù)庫是千人基因組三期的結(jié)果，千人基因組提供的 vcf 結(jié)果是用 hs37d5 作為 reference 的，所以我后面的操作也是用 hs37d5 作為reference，這也是為什么前面用 hs37d5 的原因。
• 第一步：
  • 下載 1000G VCF，原文提供的鏈接有問題，這是我自己找的，東西是一樣的，索引文件就不用下載了，這個自帶的索引文件有問題，我在進(jìn)行下一步的時候報錯了，所以我們自己重建一個索引。
  • hs37d5在上面給過下載鏈接，下載后解壓就行。

$ wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20130502/ALL.wgs.phase3_shapeit2_mvncall_integrated_v5b.20130502.sites.vcf.gz
# 建立索引,tabix是bcftools下載完后自帶的
$ tabix ALL.wgs.phase3_shapeit2_mvncall_integrated_v5b.20130502.sites.vcf.gz

• 第二步：
  • 建立索引，這一步有五句命令，每一句都很費時間，請盡量放在后臺操作，用nohup …… &。
  • parallel 請自行安裝一下，我是Ubuntu的平臺，直接使用sudo apt-get install parallele就行。
  • 下面一共五步，注意vg腳本位置，我是直接加到了環(huán)境變量，所以可以直接用，請保留盡量大的空間用于輸出，最后的結(jié)果有80G+。
  • 第一步注意 reference 的位置和 vcf 的位置，vcf我直接放在了當(dāng)前目錄下。

  • 第五步注意建立一個 tmp 目錄（-b參數(shù)），要保證這個目錄的空間有2T，這一步挺勸退的，我也是刪了好多東西才在一個盤里湊出那么多空間，如果沒有那么大的空間就會報錯，如下：

    
    
    報錯

#1
(seq 1 22; echo X; echo Y) | parallel -j 24 "time vg construct -C -R {} -r /your/path/of/b37/human_g1k_v37_decoy.fasta -v ALL.wgs.phase3_shapeit2_mvncall_integrated_v5b.20130502.sites.vcf.gz -t 1 -m 32 >{}.vg"
#2
vg ids -j $(for i in $(seq 1 22; echo X; echo Y); do echo $i.vg; done)
#3
vg index -x wg.xg $(for i in $(seq 22; echo X; echo Y); do echo $i.vg; done)
#4
for chr in $(seq 1 22; echo X; echo Y);
do
    vg prune -r $chr.vg > $chr.pruned.vg
done
#5
vg index -b ./tmp -t 24 -g wg.gcsa $(for i in $(seq 22; echo X; echo Y); do echo $i.pruned.vg; done)

• 接下來就可以進(jìn)行比對操作了，也是以1對數(shù)據(jù)為例貼一下命令:
  • 第一 你會發(fā)現(xiàn)這個命令沒有輸入reference，我覺得那是因為信息都在 index 文件里。
  • 第二 你會發(fā)現(xiàn)輸出文件也不是bam文件，之前也說過 vg 有一套自己的 pipeline，基本可以完成從 fastq 到 vcf，甚至可以無需 samtools 參與，因此也定義了一套自己的文件格式，這個不打緊，后面會轉(zhuǎn)換的。
  • -t是線程數(shù)，和前面保持一致，用12線程

$ time vg map -t 12 -x wg.xg -g wg.gcsa -f /bio-analysis1/mapper_test/data/N1/NA24694_GCCAAT_L001_R1_001.fastq.gz -f /bio-analysis1/mapper_test/data/N1/NA24694_GCCAAT_L001_R2_001.fastq.gz > male_008.gam

• 專門把 gam 格式轉(zhuǎn)換為 bam 格式放到單獨的一步是考慮到 vg 軟件的pipeline中是不需要 bam 文件的，如果計算進(jìn)時間確實也太欺負(fù)這個軟件了。但是統(tǒng)計必須得轉(zhuǎn)換到 bam 文件才可以。

$ vg surject -x wg.xg -b in.gam > out.bam

Graph Genome

• 不同于前兩個軟件，這個軟件是一個商業(yè)公司開發(fā)的，也不是開源的，但是可以通過鏈接進(jìn)入申請下載。下載下來的也不是軟件安裝包，只是一個 docker 文件，不過相對來說，操作都比 vg 友好一些，下面操作都假設(shè)已經(jīng)申請到軟件（限于版權(quán)，我就不把軟件分享出來了，感興趣的自行申請）。
• 下載并解壓以后文件結(jié)構(gòu)是這樣的。

/biosoft/graph_genome
├── docker-bpa-0.9.1.tar
├── docker-rasm-0.5.20.tar
├── example_human_Illumina.pe_1.fastq
├── example_human_Illumina.pe_2.fastq
├── LICENSE.pdf
├── manifest.json
├── opensource-software-terms-and-copyrights.md
├── repositories
├── SAMPLE.sort.vcf
└── SBG.Graph.B37.V6.rc6.vcf.gz

• 主要是三個文件：docker-bpa-0.9.1.tar,docker-rasm-0.5.20.tar和SBG.Graph.B37.V6.rc6.vcf.gz。前兩個用于安裝軟件，后面這個是官方先驗 vcf 文件，與 vg 類似。
  • 優(yōu)點：提供了文件，且使用方便，無需建立索引等
  • 缺點：只有 hs37d5 的對應(yīng) vcf，如果有其余版本需求，需要自行建立。
• 首先是安裝：

$ docker load --input docker-bpa-0.9.1.tar
$ docker load --input docker-rasm-0.5.20.tar
# 查看安裝完的鏡像信息
$ docker images

• 然后就可以直接使用了

$ time docker run -v "/your/path/of/fastq/":"/input" -v "/your/path/of/vcf/":"/file" -v "/your/path/of/reference/":"/ref" gral-bpa:"0.9.1" /usr/local/bin/aligner --vcf /file/SBG.Graph.B37.V6.rc6.vcf.gz --reference /ref/human_g1k_v37_decoy.fasta -q /input/NA24694_GCCAAT_L001_R1_001.fastq.gz -Q /input/NA24694_GCCAAT_L001_R2_001.fastq.gz -o /file/out/male_008.bam --threads 12

結(jié)果比較

• 首先說明一下，在測試的時候發(fā)生了一些問題，用 fastp 對原始數(shù)據(jù)進(jìn)行質(zhì)控以后，Graph Genome 報錯了，但是原始數(shù)據(jù)進(jìn)行 mapping 就沒有報錯，看來這個軟件存在著一些 bug，由于這個軟件沒有對應(yīng)的社區(qū)可以直接交流，我只能把問題發(fā)給開發(fā)者。所以以下結(jié)果都是未對數(shù)據(jù)進(jìn)行質(zhì)控直接 mapping 的結(jié)果。
• 從以下幾個方面進(jìn)行比較：
  1. 首先是時間比較，我用time來粗略比較了軟件的運行時間，但是如果要嚴(yán)格的從軟件的性能角度來比較的話，單單這樣是不夠的，還要考慮其他很多因素，包括占用的內(nèi)存等，這個我也不是很懂，還得請教計算機(jī)專業(yè)的同學(xué)，但是我們從結(jié)果能大致比較軟件運行的速度。下面這張表格是軟件的運行時間
     
     運行時間
  • 單從時間上來說，vg 直接 out 出局，慢的不是一星半點。這個時間不是絕對的，跟服務(wù)器當(dāng)時的性能也有關(guān)系，但是大體還是能反映一些事情的。
  • graph genome 也沒有文獻(xiàn)所說的那么好，bwa-mem2 的速度一騎絕塵，但是如果用 bwa mem 直接跑的話，所花的時間大概是現(xiàn)在的兩倍，也就是說 graph genome 和未更新的 bwa 在不考慮內(nèi)存等其余條件下，速度是相差不大的。

  2. 其次我注意了一下結(jié)果占用的空間。

     
     
     vg
     
     
     bwa
     
     
     graph genome
  • 值得注意的是，gam 在 vg 中作為一個替代了 bam 的文件，且不說這個格式是不是更加的科學(xué)，但這個占用的空間確實大了很大，這絕對是一個敗筆。
  3. 結(jié)果比較，這里用 RSeQC 和 samtools flagstat 進(jìn)行統(tǒng)計。
  • 因為八個樣的結(jié)果在軟件的差異上表現(xiàn)的都差不多，所以我就只放一個結(jié)果的比較作為例子。
  • 下面是 RSeQC 的結(jié)果：

# bwa-mem2
#==================================================
#All numbers are READ count
#==================================================

Total records:                          8024912

QC failed:                              0
Optical/PCR duplicate:                  0
Non primary hits                        0
Unmapped reads:                         57103
mapq < mapq_cut (non-unique):           585944

mapq >= mapq_cut (unique):              7381865
Read-1:                                 3726141
Read-2:                                 3655724
Reads map to '+':                       3691175
Reads map to '-':                       3690690
Non-splice reads:                       7381865
Splice reads:                           0
Reads mapped in proper pairs:           7272853
Proper-paired reads map to different chrom:1110

# graph genome
#==================================================
#All numbers are READ count
#==================================================

Total records:                          8000000

QC failed:                              0
Optical/PCR duplicate:                  0
Non primary hits                        0
Unmapped reads:                         165747
mapq < mapq_cut (non-unique):           513535

mapq >= mapq_cut (unique):              7320718
Read-1:                                 3728036
Read-2:                                 3592682
Reads map to '+':                       3660744
Reads map to '-':                       3659974
Non-splice reads:                       7320718
Splice reads:                           0
Reads mapped in proper pairs:           7137762
Proper-paired reads map to different chrom:0

# vg
#==================================================
#All numbers are READ count
#==================================================

Total records:                          8000000

QC failed:                              0
Optical/PCR duplicate:                  0
Non primary hits                        0
Unmapped reads:                         170272
mapq < mapq_cut (non-unique):           738039

mapq >= mapq_cut (unique):              7091689
Read-1:                                 0
Read-2:                                 0
Reads map to '+':                       3545844
Reads map to '-':                       3545845
Non-splice reads:                       7091689
Splice reads:                           0
Reads mapped in proper pairs:           0
Proper-paired reads map to different chrom:0

• 下面是 samtools flagstat 的結(jié)果:

# bwa-mem2
8024912 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
24912 + 0 supplementary
0 + 0 duplicates
7967809 + 0 mapped (99.29% : N/A)
8000000 + 0 paired in sequencing
4000000 + 0 read1
4000000 + 0 read2
7746532 + 0 properly paired (96.83% : N/A)
7886064 + 0 with itself and mate mapped
56833 + 0 singletons (0.71% : N/A)
64324 + 0 with mate mapped to a different chr
45865 + 0 with mate mapped to a different chr (mapQ>=5)

# graph genome
8000000 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
7834253 + 0 mapped (97.93% : N/A)
8000000 + 0 paired in sequencing
4000000 + 0 read1
4000000 + 0 read2
7593524 + 0 properly paired (94.92% : N/A)
7671608 + 0 with itself and mate mapped
162645 + 0 singletons (2.03% : N/A)
23810 + 0 with mate mapped to a different chr
17071 + 0 with mate mapped to a different chr (mapQ>=5)

# vg
8000000 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
7829728 + 0 mapped (97.87% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

• 從這個結(jié)果上似乎是 bwa-mem2 > graph genome > vg
• 可以注意到，graph genome 和 vg 在某些統(tǒng)計上是異于 bwa 的，比如Total records。
• 相對來說，bwa 的結(jié)果和 graph genome 的更接近，而 vg 的結(jié)果比較奇怪。我的理解是，vg 的 pipeline 里本來是不需要 bam 的，用的是 gam 格式，而我通過官方的腳本轉(zhuǎn)換成 bam 只是一種“兼容”轉(zhuǎn)換，gam 的信息和 bam 應(yīng)該是不完全一樣的，所以通過“兼容”的方式轉(zhuǎn)換過來，信息也不完全一樣。

• 這里乍一眼看似乎都是 bwa 效果好，但其實相差也不大。文獻(xiàn)也提到在某些區(qū)域，graph based 的軟件 mapping 效果更好，但是這些我也沒法查看，或者是我也不知道怎么去查看，反正文獻(xiàn)就這么寫了。

  
  
  來自 graph genome 文獻(xiàn)
• 其實效果怎么樣，光從 bam 文件看也不準(zhǔn)確，還是需要 call 出來 variation 比較比較。那這個等以后有機(jī)會再試試吧。

??水平有限，要是存在什么錯誤請指出，可發(fā)送郵件至shiyuant@outlook.com！請大家多多批評指正，相互交流，共同成長，謝謝?。。?/p>