WGS分析笔记（7）—— 新的mapping工具

• 开始之前特别感谢从美国“人肉背回”这两个软件的小伙伴！@Olivia阿仪_鸦雀
  

---

• 之前讲过目前用的比较多的两个二代测序基因组比对软件，bwa 和 bowtie2，也比较了一下效果，总的来说两者的效果相差不大，但是 bwa 可能更加的准确一些，而很多文章以及一些组学计划都是用的 bwa。
• 今天主要来比较一下一些新的 mapping 工具，vg 和 Graph Genome，相对于传统基于本文进行比对的工具（bwa等），这类工具基于 Graph 进行比对，卖点是速度更快，更加准确。论文链接如下，文章是开放的，直接可以下载：
  • Variation graph toolkit improves read mapping by representing genetic variation in the reference
  • Fast and accurate genomic analyses using genome graphs
• 文章都是发表在高影响因子的 nature 子刊上，基本都是踩着 bwa 上位的，那我就想测试一下这俩软件（在我使用的时候我查了一下百度，发现国内根本没人用这玩意儿）有没有文章说的那么好，再加之 bwa 更新了新的 mem2 ，其速度又有大幅度的提升，但是结果没有什么改变。
• 首先是测试数据，这里我用GIAB（瓶中基因组计划，具体看链接文章，也是开放论文，总的这三篇论文篇幅都很短，感兴趣的不妨读一读。）中的 ChineseTrio 中的父母的 illumina 测序数据，随机取了两个个体的每个批次（共四个批次）中的一组数据，就是总共8对 fastq 文件。

数据预处理

• 因为我们直接下载的是 GIAB 中的 fastq 文件，常规步骤还是需要在 mapping 前进行一下质控处理的，这里的话可以参考我之前的文章。

$ 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 为例贴的命令行，但没有改变输出后的名字，只是修改了输出路径，也方便后续对照数据的名字。
• 这里不贴结果了，但是可以说明一下，数据的质量很高的，即使不进行这一步质控差别也不大，而且用的是 pcr-free 的建库方式，所以基本没有什么duplicate。
• 下面就介绍一下三个软件的使用方式以及测试结果。

bwa mem2

• 标题链接即新版本 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

• 直接下载解压就可以用了，而且比对的参数和原来的基本没什么变化，当然也可以不加什么参数，直接按照上面的来进行。
• 值得注意的一点是，新的 bwa 需要重新建立索引，速度还是比较久的，但是这种建立一次，终身管用的操作，也无需特别计较时间。
• 除此之外就是，新的软件建立的索引文件和原来的索引文件有所区别的（部分一样），而且新的索引文件着实耗费空间，我以 hs37d5 作为参考基因组进行实验（不同参考基因组之间的区别可以见这里）：
  reference 索引文件
• 可以看到以.2bit.64和.8bit.32这俩文件很占空间，(这里漏截一个749M大小的.pac文件)，这里我之所以用 hs37d5 （b37_decoy）作为参考基因组也是和另外两个软件有关，而且就主体序列而言，hs37d5 和 hg19 没什么区别，所以不用过度纠结。

# 两个下载该参考基因组的地址，一个是 EBI 网站上的，一个是 GATK 上的，没啥区别，除了文件名，下载后请自行解压。
$ 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对数据中的一对为例，把命令行贴一下：

$ 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 是线程数，三个软件我都用12线程进行运行
# -R 是头文件，这里可以不设置
# time命令用来统计时间

vg

• 标题即 github 链接，软件提供了打包好的二进制文件，也可以不安装直接用。

$ 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 工具，这工具也发表在了高影响因子的 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 为止，看着还挺复杂的，最主要的是，你会发现第一步有一个 vcf 文件，这是 vg 所需要的先验文件，以你所要使用的 reference 得到的 vcf，于是就有了一个先有鸡还是先有蛋的问题……好了，开玩笑的，但是有一点确实很让人诟病，软件也没说这个 vcf 是任意一个 vcf 还是需要越大的群体得到的结果越好，毕竟不同个体的变异肯定是相差很大的，如果是群体的 vcf，那么作为一个软件，像 GATK 一样给一个参考的数据库岂不是更好，Graph Genome 也需要先验 vcf 文件，但是 Graph Genome 是提供了的。
• 官方的 README 文件的说明简单的令人发指，找了半天终于找到了一个通过公共数据库建立索引文件的文章。
• 文章中每一步写的还是很详细的，包括可能需要的运行时间，临时空间大小，产生文件大小。但就是因为写的越细才越让人生气，因为这部分的计算资源还是挺大的，官方既然自己做过一遍，为什么不直接把结果文件挂出来让人下载作为使用参考，就像 GATK Bundle 一样。
• 气归气，接下来还是一步一步演示一下，不过这个流程还是很花费时间的。
• 首先说一下，这里用的公共数据库是千人基因组三期的结果，千人基因组提供的 vcf 结果是用 hs37d5 作为 reference 的，所以我后面的操作也是用 hs37d5 作为reference，这也是为什么前面用 hs37d5 的原因。
• 第一步：
  • 下载 1000G VCF，原文提供的链接有问题，这是我自己找的，东西是一样的，索引文件就不用下载了，这个自带的索引文件有问题，我在进行下一步的时候报错了，所以我们自己重建一个索引。
  • 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脚本位置，我是直接加到了环境变量，所以可以直接用，请保留尽量大的空间用于输出，最后的结果有80G+。
  • 第一步注意 reference 的位置和 vcf 的位置，vcf我直接放在了当前目录下。

  • 第五步注意建立一个 tmp 目录（-b参数），要保证这个目录的空间有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)

• 接下来就可以进行比对操作了，也是以1对数据为例贴一下命令:
  • 第一 你会发现这个命令没有输入reference，我觉得那是因为信息都在 index 文件里。
  • 第二 你会发现输出文件也不是bam文件，之前也说过 vg 有一套自己的 pipeline，基本可以完成从 fastq 到 vcf，甚至可以无需 samtools 参与，因此也定义了一套自己的文件格式，这个不打紧，后面会转换的。
  • -t是线程数，和前面保持一致，用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 格式转换为 bam 格式放到单独的一步是考虑到 vg 软件的pipeline中是不需要 bam 文件的，如果计算进时间确实也太欺负这个软件了。但是统计必须得转换到 bam 文件才可以。

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

Graph Genome

• 不同于前两个软件，这个软件是一个商业公司开发的，也不是开源的，但是可以通过链接进入申请下载。下载下来的也不是软件安装包，只是一个 docker 文件，不过相对来说，操作都比 vg 友好一些，下面操作都假设已经申请到软件（限于版权，我就不把软件分享出来了，感兴趣的自行申请）。
• 下载并解压以后文件结构是这样的。

/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 类似。
  • 优点：提供了文件，且使用方便，无需建立索引等
  • 缺点：只有 hs37d5 的对应 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

结果比较

• 首先说明一下，在测试的时候发生了一些问题，用 fastp 对原始数据进行质控以后，Graph Genome 报错了，但是原始数据进行 mapping 就没有报错，看来这个软件存在着一些 bug，由于这个软件没有对应的社区可以直接交流，我只能把问题发给开发者。所以以下结果都是未对数据进行质控直接 mapping 的结果。
• 从以下几个方面进行比较：
  1. 首先是时间比较，我用time来粗略比较了软件的运行时间，但是如果要严格的从软件的性能角度来比较的话，单单这样是不够的，还要考虑其他很多因素，包括占用的内存等，这个我也不是很懂，还得请教计算机专业的同学，但是我们从结果能大致比较软件运行的速度。下面这张表格是软件的运行时间
     运行时间
  • 单从时间上来说，vg 直接 out 出局，慢的不是一星半点。这个时间不是绝对的，跟服务器当时的性能也有关系，但是大体还是能反映一些事情的。
  • graph genome 也没有文献所说的那么好，bwa-mem2 的速度一骑绝尘，但是如果用 bwa mem 直接跑的话，所花的时间大概是现在的两倍，也就是说 graph genome 和未更新的 bwa 在不考虑内存等其余条件下，速度是相差不大的。

  2. 其次我注意了一下结果占用的空间。

     
     vg
     bwa
     graph genome
  • 值得注意的是，gam 在 vg 中作为一个替代了 bam 的文件，且不说这个格式是不是更加的科学，但这个占用的空间确实大了很大，这绝对是一个败笔。
  3. 结果比较，这里用 RSeQC 和 samtools flagstat 进行统计。
  • 因为八个样的结果在软件的差异上表现的都差不多，所以我就只放一个结果的比较作为例子。
  • 下面是 RSeQC 的结果：

# 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 的结果:

# 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)

• 从这个结果上似乎是 bwa-mem2 > graph genome > vg
• 可以注意到，graph genome 和 vg 在某些统计上是异于 bwa 的，比如Total records。
• 相对来说，bwa 的结果和 graph genome 的更接近，而 vg 的结果比较奇怪。我的理解是，vg 的 pipeline 里本来是不需要 bam 的，用的是 gam 格式，而我通过官方的脚本转换成 bam 只是一种“兼容”转换，gam 的信息和 bam 应该是不完全一样的，所以通过“兼容”的方式转换过来，信息也不完全一样。

• 这里乍一眼看似乎都是 bwa 效果好，但其实相差也不大。文献也提到在某些区域，graph based 的软件 mapping 效果更好，但是这些我也没法查看，或者是我也不知道怎么去查看，反正文献就这么写了。

  
  来自 graph genome 文献
• 其实效果怎么样，光从 bam 文件看也不准确，还是需要 call 出来 variation 比较比较。那这个等以后有机会再试试吧。

---

  水平有限，要是存在什么错误请指出，可发送邮件至shiyuant@outlook.com！请大家多多批评指正，相互交流，共同成长，谢谢！！！