- 開(kāi)始之前特別感謝從美國(guó)“人肉背回”這兩個(gè)軟件的小伙伴！@Olivia阿儀_鴉雀
- 之前講過(guò)目前用的比較多的兩個(gè)二代測(cè)序基因組比對(duì)軟件茫藏，bwa 和 bowtie2缚态，也比較了一下效果，總的來(lái)說(shuō)兩者的效果相差不大，但是 bwa 可能更加的準(zhǔn)確一些粹庞，而很多文章以及一些組學(xué)計(jì)劃都是用的 bwa。
- 今天主要來(lái)比較一下一些新的 mapping 工具洽损，vg 和 Graph Genome庞溜，相對(duì)于傳統(tǒng)基于本文進(jìn)行比對(duì)的工具（bwa等），這類(lèi)工具基于 Graph 進(jìn)行比對(duì)碑定，賣(mài)點(diǎn)是速度更快流码，更加準(zhǔn)確。論文鏈接如下延刘，文章是開(kāi)放的漫试，直接可以下載：
- 文章都是發(fā)表在高影響因子的 nature 子刊上，基本都是踩著 bwa 上位的碘赖，那我就想測(cè)試一下這倆軟件（在我使用的時(shí)候我查了一下百度驾荣，發(fā)現(xiàn)國(guó)內(nèi)根本沒(méi)人用這玩意兒）有沒(méi)有文章說(shuō)的那么好，再加之 bwa 更新了新的 mem2 普泡，其速度又有大幅度的提升播掷，但是結(jié)果沒(méi)有什么改變。
- 首先是測(cè)試數(shù)據(jù)撼班，這里我用GIAB（瓶中基因組計(jì)劃歧匈，具體看鏈接文章，也是開(kāi)放論文砰嘁，總的這三篇論文篇幅都很短件炉，感興趣的不妨讀一讀。）中的 ChineseTrio 中的父母的 illumina 測(cè)序數(shù)據(jù)矮湘，隨機(jī)取了兩個(gè)個(gè)體的每個(gè)批次（共四個(gè)批次）中的一組數(shù)據(jù)斟冕，就是總共8對(duì) fastq 文件。
數(shù)據(jù)預(yù)處理
- 因?yàn)槲覀冎苯酉螺d的是 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
- 這里我是以其中一對(duì) fastq 為例貼的命令行券时，但沒(méi)有改變輸出后的名字孤里，只是修改了輸出路徑，也方便后續(xù)對(duì)照數(shù)據(jù)的名字橘洞。
- 這里不貼結(jié)果了捌袜，但是可以說(shuō)明一下，數(shù)據(jù)的質(zhì)量很高的炸枣，即使不進(jìn)行這一步質(zhì)控差別也不大虏等，而且用的是 pcr-free 的建庫(kù)方式弄唧，所以基本沒(méi)有什么duplicate。
- 下面就介紹一下三個(gè)軟件的使用方式以及測(cè)試結(jié)果霍衫。
bwa mem2
- 標(biāo)題鏈接即新版本 bwa mem 的 github 地址候引，使用非常之簡(jiǎn)單。
$ 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
- 直接下載解壓就可以用了敦跌，而且比對(duì)的參數(shù)和原來(lái)的基本沒(méi)什么變化澄干，當(dāng)然也可以不加什么參數(shù)，直接按照上面的來(lái)進(jìn)行柠傍。
- 值得注意的一點(diǎn)是麸俘，新的 bwa 需要重新建立索引，速度還是比較久的惧笛，但是這種建立一次从媚，終身管用的操作，也無(wú)需特別計(jì)較時(shí)間患整。
- 除此之外就是拜效，新的軟件建立的索引文件和原來(lái)的索引文件有所區(qū)別的（部分一樣），而且新的索引文件著實(shí)耗費(fèi)空間并级，我以 hs37d5 作為參考基因組進(jìn)行實(shí)驗(yàn)（不同參考基因組之間的區(qū)別可以見(jiàn)這里）：
reference 索引文件 - 可以看到以
.2bit.64和.8bit.32這倆文件很占空間拂檩，(這里漏截一個(gè)749M大小的.pac文件)，這里我之所以用 hs37d5 （b37_decoy）作為參考基因組也是和另外兩個(gè)軟件有關(guān)嘲碧，而且就主體序列而言，hs37d5 和 hg19 沒(méi)什么區(qū)別父阻，所以不用過(guò)度糾結(jié)愈涩。
# 兩個(gè)下載該參考基因組的地址，一個(gè)是 EBI 網(wǎng)站上的加矛，一個(gè)是 GATK 上的履婉，沒(méi)啥區(qū)別，除了文件名斟览，下載后請(qǐng)自行解壓毁腿。
$ 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對(duì)數(shù)據(jù)中的一對(duì)為例，把命令行貼一下：
$ 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ù)苛茂，三個(gè)軟件我都用12線程進(jìn)行運(yùn)行
# -R 是頭文件已烤，這里可以不設(shè)置
# time命令用來(lái)統(tǒng)計(jì)時(shí)間
vg
- 標(biāo)題即 github 鏈接，軟件提供了打包好的二進(jìn)制文件妓羊，也可以不安裝直接用胯究。
$ wget https://github.com/vgteam/vg/releases/download/v1.22.0/vg
$ chmod +x vg
- 值得一提的是，無(wú)論是 vg 還是 Graph Genome 都有自己的一套 pipeline 躁绸，其中包含了 calling 的功能裕循，尤其是 vgtools臣嚣，不僅能 call SNVs/INDELs 還有專(zhuān)門(mén)的call SVs 工具，這工具也發(fā)表在了高影響因子的 Genome Biology 上剥哑，Genotyping structural variants in
pangenome graphs using the vg toolkit硅则。而 Graph Genome 除了具備 calling 功能（一個(gè)命令同時(shí) call snv/indel/sv ）也提供了一套基于trio的分析工具VBT-TrioAnalysis。 -
但是vg真是一款令人火大的軟件株婴，在我測(cè)試開(kāi)始抢埋，到寫(xiě)下這篇文章，都是懷著一種悲憤的心情督暂。
vg pipeline - 這是官方提供的pipeline的上半部分揪垄，到 mapping 為止，看著還挺復(fù)雜的逻翁，最主要的是饥努，你會(huì)發(fā)現(xiàn)第一步有一個(gè) vcf 文件，這是 vg 所需要的先驗(yàn)文件八回，以你所要使用的 reference 得到的 vcf酷愧，于是就有了一個(gè)先有雞還是先有蛋的問(wèn)題……好了，開(kāi)玩笑的缠诅，但是有一點(diǎn)確實(shí)很讓人詬病溶浴，軟件也沒(méi)說(shuō)這個(gè) vcf 是任意一個(gè) vcf 還是需要越大的群體得到的結(jié)果越好，畢竟不同個(gè)體的變異肯定是相差很大的管引，如果是群體的 vcf士败，那么作為一個(gè)軟件，像 GATK 一樣給一個(gè)參考的數(shù)據(jù)庫(kù)豈不是更好褥伴，Graph Genome 也需要先驗(yàn) vcf 文件谅将，但是 Graph Genome 是提供了的。
- 官方的 README 文件的說(shuō)明簡(jiǎn)單的令人發(fā)指重慢，找了半天終于找到了一個(gè)通過(guò)公共數(shù)據(jù)庫(kù)建立索引文件的文章饥臂。
- 文章中每一步寫(xiě)的還是很詳細(xì)的，包括可能需要的運(yùn)行時(shí)間似踱，臨時(shí)空間大小隅熙，產(chǎn)生文件大小。但就是因?yàn)閷?xiě)的越細(xì)才越讓人生氣核芽，因?yàn)檫@部分的計(jì)算資源還是挺大的囚戚，官方既然自己做過(guò)一遍，為什么不直接把結(jié)果文件掛出來(lái)讓人下載作為使用參考狞洋，就像 GATK Bundle 一樣弯淘。
- 氣歸氣，接下來(lái)還是一步一步演示一下吉懊，不過(guò)這個(gè)流程還是很花費(fèi)時(shí)間的庐橙。
- 首先說(shuō)一下假勿，這里用的公共數(shù)據(jù)庫(kù)是千人基因組三期的結(jié)果，千人基因組提供的 vcf 結(jié)果是用 hs37d5 作為 reference 的态鳖，所以我后面的操作也是用 hs37d5 作為reference转培，這也是為什么前面用 hs37d5 的原因。
- 第一步：
- 下載 1000G VCF浆竭，原文提供的鏈接有問(wèn)題浸须，這是我自己找的，東西是一樣的邦泄，索引文件就不用下載了删窒，這個(gè)自帶的索引文件有問(wèn)題，我在進(jìn)行下一步的時(shí)候報(bào)錯(cuò)了顺囊，所以我們自己重建一個(gè)索引肌索。
- hs37d5在上面給過(guò)下載鏈接，下載后解壓就行特碳。
$ 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
- 第二步：
- 建立索引诚亚，這一步有五句命令，每一句都很費(fèi)時(shí)間午乓，請(qǐng)盡量放在后臺(tái)操作站宗，用
nohup …… &。 - parallel 請(qǐng)自行安裝一下益愈，我是Ubuntu的平臺(tái)梢灭，直接使用
sudo apt-get install parallele就行。 - 下面一共五步腕唧，注意vg腳本位置或辖，我是直接加到了環(huán)境變量，所以可以直接用枣接，請(qǐng)保留盡量大的空間用于輸出，最后的結(jié)果有80G+缺谴。
- 第一步注意 reference 的位置和 vcf 的位置但惶，vcf我直接放在了當(dāng)前目錄下。
-
第五步注意建立一個(gè) tmp 目錄（-b參數(shù)）湿蛔，要保證這個(gè)目錄的空間有2T膀曾，這一步挺勸退的，我也是刪了好多東西才在一個(gè)盤(pán)里湊出那么多空間阳啥，如果沒(méi)有那么大的空間就會(huì)報(bào)錯(cuò)添谊，如下：
報(bào)錯(cuò)
- 建立索引诚亚，這一步有五句命令，每一句都很費(fèi)時(shí)間午乓，請(qǐng)盡量放在后臺(tái)操作站宗，用
#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)
- 接下來(lái)就可以進(jìn)行比對(duì)操作了，也是以1對(duì)數(shù)據(jù)為例貼一下命令:
- 第一 你會(huì)發(fā)現(xiàn)這個(gè)命令沒(méi)有輸入reference察迟，我覺(jué)得那是因?yàn)樾畔⒍荚?index 文件里斩狱。
- 第二 你會(huì)發(fā)現(xiàn)輸出文件也不是bam文件耳高，之前也說(shuō)過(guò) vg 有一套自己的 pipeline，基本可以完成從 fastq 到 vcf所踊，甚至可以無(wú)需 samtools 參與泌枪，因此也定義了一套自己的文件格式，這個(gè)不打緊秕岛，后面會(huì)轉(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
- 專(zhuān)門(mén)把 gam 格式轉(zhuǎn)換為 bam 格式放到單獨(dú)的一步是考慮到 vg 軟件的pipeline中是不需要 bam 文件的修壕，如果計(jì)算進(jìn)時(shí)間確實(shí)也太欺負(fù)這個(gè)軟件了。但是統(tǒng)計(jì)必須得轉(zhuǎn)換到 bam 文件才可以遏考。
$ vg surject -x wg.xg -b in.gam > out.bam
Graph Genome
- 不同于前兩個(gè)軟件慈鸠，這個(gè)軟件是一個(gè)商業(yè)公司開(kāi)發(fā)的，也不是開(kāi)源的诈皿，但是可以通過(guò)鏈接進(jìn)入申請(qǐng)下載林束。下載下來(lái)的也不是軟件安裝包，只是一個(gè) docker 文件稽亏，不過(guò)相對(duì)來(lái)說(shuō)壶冒，操作都比 vg 友好一些，下面操作都假設(shè)已經(jīng)申請(qǐng)到軟件（限于版權(quán)截歉，我就不把軟件分享出來(lái)了胖腾，感興趣的自行申請(qǐng)）。
- 下載并解壓以后文件結(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
- 主要是三個(gè)文件：
docker-bpa-0.9.1.tar,docker-rasm-0.5.20.tar和SBG.Graph.B37.V6.rc6.vcf.gz咸作。前兩個(gè)用于安裝軟件，后面這個(gè)是官方先驗(yàn) vcf 文件宵睦，與 vg 類(lèi)似记罚。- 優(yōu)點(diǎn)：提供了文件，且使用方便壳嚎，無(wú)需建立索引等
- 缺點(diǎn)：只有 hs37d5 的對(duì)應(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é)果比較
- 首先說(shuō)明一下说庭，在測(cè)試的時(shí)候發(fā)生了一些問(wèn)題会通，用 fastp 對(duì)原始數(shù)據(jù)進(jìn)行質(zhì)控以后采蚀，Graph Genome 報(bào)錯(cuò)了欢顷，但是原始數(shù)據(jù)進(jìn)行 mapping 就沒(méi)有報(bào)錯(cuò)淆院，看來(lái)這個(gè)軟件存在著一些 bug舷礼，由于這個(gè)軟件沒(méi)有對(duì)應(yīng)的社區(qū)可以直接交流赴邻，我只能把問(wèn)題發(fā)給開(kāi)發(fā)者握恳。所以以下結(jié)果都是未對(duì)數(shù)據(jù)進(jìn)行質(zhì)控直接 mapping 的結(jié)果迟隅。
- 從以下幾個(gè)方面進(jìn)行比較：
- 首先是時(shí)間比較，我用
time來(lái)粗略比較了軟件的運(yùn)行時(shí)間攻礼，但是如果要嚴(yán)格的從軟件的性能角度來(lái)比較的話业踢，單單這樣是不夠的，還要考慮其他很多因素礁扮，包括占用的內(nèi)存等知举，這個(gè)我也不是很懂，還得請(qǐng)教計(jì)算機(jī)專(zhuān)業(yè)的同學(xué)太伊，但是我們從結(jié)果能大致比較軟件運(yùn)行的速度雇锡。下面這張表格是軟件的運(yùn)行時(shí)間
運(yùn)行時(shí)間
- 單從時(shí)間上來(lái)說(shuō)，vg 直接 out 出局僚焦，慢的不是一星半點(diǎn)锰提。這個(gè)時(shí)間不是絕對(duì)的，跟服務(wù)器當(dāng)時(shí)的性能也有關(guān)系芳悲，但是大體還是能反映一些事情的立肘。
- graph genome 也沒(méi)有文獻(xiàn)所說(shuō)的那么好，bwa-mem2 的速度一騎絕塵名扛，但是如果用 bwa mem 直接跑的話谅年，所花的時(shí)間大概是現(xiàn)在的兩倍，也就是說(shuō) graph genome 和未更新的 bwa 在不考慮內(nèi)存等其余條件下肮韧，速度是相差不大的融蹂。
-
其次我注意了一下結(jié)果占用的空間。
vg
bwa
graph genome
- 值得注意的是弄企，gam 在 vg 中作為一個(gè)替代了 bam 的文件超燃，且不說(shuō)這個(gè)格式是不是更加的科學(xué)，但這個(gè)占用的空間確實(shí)大了很大拘领，這絕對(duì)是一個(gè)敗筆意乓。
- 結(jié)果比較，這里用 RSeQC 和 samtools flagstat 進(jìn)行統(tǒng)計(jì)约素。
- 因?yàn)榘藗€(gè)樣的結(jié)果在軟件的差異上表現(xiàn)的都差不多洽瞬，所以我就只放一個(gè)結(jié)果的比較作為例子。
- 下面是 RSeQC 的結(jié)果：
- 首先是時(shí)間比較，我用
# 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)
- 從這個(gè)結(jié)果上似乎是 bwa-mem2 > graph genome > vg
- 可以注意到业汰，graph genome 和 vg 在某些統(tǒng)計(jì)上是異于 bwa 的，比如
Total records菩颖。 - 相對(duì)來(lái)說(shuō)样漆，bwa 的結(jié)果和 graph genome 的更接近，而 vg 的結(jié)果比較奇怪晦闰。我的理解是放祟，vg 的 pipeline 里本來(lái)是不需要 bam 的鳍怨，用的是 gam 格式，而我通過(guò)官方的腳本轉(zhuǎn)換成 bam 只是一種“兼容”轉(zhuǎn)換跪妥，gam 的信息和 bam 應(yīng)該是不完全一樣的鞋喇，所以通過(guò)“兼容”的方式轉(zhuǎn)換過(guò)來(lái)，信息也不完全一樣眉撵。
-
這里乍一眼看似乎都是 bwa 效果好侦香，但其實(shí)相差也不大。文獻(xiàn)也提到在某些區(qū)域纽疟，graph based 的軟件 mapping 效果更好罐韩，但是這些我也沒(méi)法查看，或者是我也不知道怎么去查看污朽，反正文獻(xiàn)就這么寫(xiě)了散吵。
來(lái)自 graph genome 文獻(xiàn) - 其實(shí)效果怎么樣，光從 bam 文件看也不準(zhǔn)確蟆肆，還是需要 call 出來(lái) variation 比較比較矾睦。那這個(gè)等以后有機(jī)會(huì)再試試吧。
??水平有限炎功，要是存在什么錯(cuò)誤請(qǐng)指出枚冗，可發(fā)送郵件至shiyuant@outlook.com！請(qǐng)大家多多批評(píng)指正亡问，相互交流官紫，共同成長(zhǎng)，謝謝Ｖ菖骸Ｊ馈！
