gvcf-vcf文件合并

转载于http://blog.sciencenet.cn/blog-1469385-829144.html，该帖子详细介绍了gatk2和gatk3的不同点。本篇文章对该帖子部分内容进行修改，只保留gvcf生成条件下变异检测的流程及vcf文件合并的相关说明。

gatk3中检测变异从原来的gatk2的Joint Variant Calling变成了现在的Joint Genotyping，一方面顺利的解决了变异检测过程中的“N+1”难题，另一方面则大大减少了运行所需要的时间和消耗的资源。

在这种模式下只需要对每个sample的bam文件运行HaplotypeCaller，生成每个sample的gVCF文件，之后再合并这些文件进行genotyping，这样原先最耗时间的步骤(HaplotypeCaller这一步)从与sample数量的指数级关系就变成了线性关系，而合并genotyping消耗的时间则相对很少，因此整个模式带来的将是一次巨大的革新。

1. Variant calling

针对每个sample的BAM文件运行HaplotypeCaller（每个sample可以有多个文件(应当包含同样的SM标签)，运行的时候通过-I参数设置一起提供给caller），然后每个sample会生成一个gVCFs文件，需要设置的参数有：

 

--emitRefConfidenceGVCF 通过这个选项来输出与reference一致的碱基的相关信息，这样最终的结果其实就相当于包含了所有的位点，而不再仅仅是原先的variants位点。

--variant_index_type LINEAR 这个选项是主程序里面的，具体可以不管

--variant_index_parameter 128000 跟上面的选项相关的一个参数，也不用管

其实这个思路还是很明显的，原先情况下单纯的输出variants位点再合并会导致多出来的reference位点的实际情况无法确定，所以把所有位点一起输出再合并就可以比较了。因此理论上这种做法在原来的框架下也是能做的，因为本身大部分的caller就可以输出所有位点的信息，但是原来的做法会有以下几个问题:

(1)生成文件的体积会非常大，很容易就超过原来的bam文件，即使压缩后也很大，一方面极度浪费资源，另一方面下游处理效率会很差；

(2)通过直接输出的这些reference位点的信息可能不准，或者不完全，因为当初的目的可能更多的是用于debug之类；

(3)很多情况下通过合并得到的multi-sample vcf文件中的结果部分值不准确，或者根本得不到更新。

这些也是造成所谓“N+1”问题的根本原因。

第三点其实跟第二点本质上是一样的，都是由于单个vcf文件中存储的每个sample的信息不全(loss)导致合并时很多值无法更新，尤其是碰到多variants位点时这个问题最明显。不过解决的思路也很简单，就是尝试将每个sample BAM文件中的有用信息尽可能无损(lossless)的转换到一个文件中，再通过这个文件提取得到最后的vcf文件。

这边我想顺便介绍一些平常可以用来进行这种vcf文件合并操作的一些工具，并且讨论一下在使用过程中可能需要注意的一些问题:

a).vcf-merge

http://vcftools./perl_module.html#vcf-merge

这个工具是vcftools里面的一个Perl脚本，具体用法如下:

vcf-mergeA.vcf.gz B.vcf.gz C.vcf.gz > out.vcf

首先最大的缺点是慢(所以后来重新开发了C版本，见后)；然后另外一个比较麻烦的就在于对多variants位点的处理，碰到这些位点的时候不会更新各个sample，例如我们可以尝试合并下面两个vcf文件：

vcf1:

##fileformat=VCFv4.1

#CHROM      POS         ID    REF ALT QUAL     FILTER   INFO       FORMAT        Sample01

1       6324        .        T C       30    PASS       .        GT:AD    1/1:1,22

1       16324      .      T A     40    PASS       .        GT:AD    0/1:10,11

4       134861    .      G C       50    PASS       .        GT:AD    0/1:7,13

vcf2:

##fileformat=VCFv4.1

#CHROM      POS         ID    REF ALT QUAL     FILTER   INFO       FORMAT        Sample02

1       6324        .        T      A      45    PASS       .        GT:AD    1/1:0,23

1       16324      .        T      C,G 35    PASS       .        GT:AD    1/2:0,12,12

4       134861    .        G      GC   15    LowQual        .        GT:AD    1/1:1,19

我们可以尝试用vcf-merge来合并这两个文件，具体代码如下:

## step1: compress vcf files and create an index (this process

## is required for most perl APIs in vcftools)

./biosoft/bgzip./test/test1.vcf &&

   ./biosoft/tabix -p vcf ./test/test1.vcf.gz

./biosoft/bgzip./test/test2.vcf &&

    ./biosoft/tabix -p vcf ./test/test2.vcf.gz

 

## step2: merge with vcf-merge

perl ./scripts/vcf-merge./test/test1.vcf.gz./test/test2.vcf.gz

    > ./test/test_merge.vcf

第一步是压缩以及生成index，vcftools的很多perl API都要求input的vcf文件是通过bgzip压缩，并且用tabix生成对应的index文件，这样可以实现对vcf文件的随机读写。

得到的结果如下:

##fileformat=VCFv4.1

##source_20140401.1=vcf-merge(r840) ./test/test1.vcf.gz ./test/test2.vcf.gz

##sourceFiles_20140401.1=0:./test/test1.vcf.gz,1:./test/test2.vcf.gz

#CHROM      POS         ID    REF ALT QUAL     FILTER   INFO       FORMAT        Sample01      Sample02

1       6324        .        T      A,C 37.50       PASS       AC=2,2;AN=4;SF=0,1 GT:AD    2/2:1,22 1/1:0,23

1       16324      .        T      C,G,A      37.50       PASS       AC=1,1,1;AN=4;SF=0,1        GT:AD    0/3:10,11        1/2:0,12,12

4       134861    .        G      GC,C        32.50      LowQual       AC=2,1;AN=4;SF=0,1f         GT:AD    0/2:7,13 1/1:1,19

从更新的结果来看，Quality值应该是取得两者的平均值，然后INFO部分增加了AC，AN以及SF的信息，FILTER的标签是以被FILTER掉的为准，然后在SF这边有显示(1f就表示在第二个sample中是LowQual的)，假设我们当时定义的LowQual的标准是30的话，这边根据情况可能就要更新成PASS；然后可以看到，这边AD值(这个值反应的是每种形态的allele的depth)的信息是没有更新的，因为每个caller生成的vcf文件中包含的INFO或者FORMAT部分的标签各种各样(这边的AD值就是HaplotypeCaller或者UnifiedGenotyper默认输出的)，除非是配套的下游操作软件，否则一般是很难考虑到所有情况的，而且这边即使尝试更新，得到的结果也不一定完全准确，例如即使是0/1的情况下，也是有存在第三种形态的可能的，只是当时不足以对genotyping产生影响所以被舍弃了而已。

b).GATKCombineVariants

http://www./gatk/gatkdocs/org_broadinstitute_sting_gatk_walkers_variantutils_CombineVariants.html

这个是GATK里面用来合并vcf文件的工具，用法如下:

## test of CombineVariants in GATK

java -jar ./biosoft/GenomeAnalysisTK-3.1-1/GenomeAnalysisTK.jar

    -R ./ref/TAIR9_chr1.random.fasta

    -T CombineVariants

    --variant ./test/test1.vcf.gz

    --variant ./test/test2.vcf.gz

    -o ./test/test_combine.vcf

    -genotypeMergeOptions UNIQUIFY

这边顺便提一下GATK的大部分工具也是直接支持通过bgzip压缩的并且通过tabix index过的vcf文件的。

然后结果如下：

#CHROM      POS         ID    REF ALT QUAL     FILTER   INFO       FORMAT        Sample01.variant         Sample02.variant2

1       6324        .        T      C,A 30    PASS       AC=2,2;AF=0.500,0.500;AN=4;set=Intersection GT   1/1   2/2

1       16324      .        T      A,C,G      40    PASS       AC=1,1,1;AF=0.250,0.250,0.250;AN=4;set=Intersection   GT   0/1         2/3

4       134861    .        G      C      50    PASS       AC=1;AF=0.500;AN=2;set=variant GT:AD    0/1:7,13   ./.

4       134861    .        G      GC   15    LowQual        AC=2;AF=1.00;AN=2;set=FilteredInAll GT:AD    ./.     1/1:1,19

然后……(一瞬间有种不想评论的感觉)，可以看到可能是也注意到AD这个值的问题了，所以前两个位点直接把AD值给删掉了(删掉了……)，后面两个因为是Indel和SNP位点重合，而Indel的实际位置本身就应该是134861+1，所以这边分成两个可能应该更准确一点，但这样vcf文件当中就出现了duplicate的位置。

c). bcftools merge

http://samtools./bcftools/bcftools.html#merge

这个工具是后来重新开发的，既包含了原来samtools里面附带的bcftools的所有功能，也包含htslib (https://github.com/samtools/htslib) 里面的所有工具，用来替代vcftools里面提供的Perl API，这边的bcftools merge就是用来替代之前提到的vcf-merge的，具体用法如下：

## test of bcftools merge

./biosoft/bcftools merge./test/test1.vcf.gz ./test/test2.vcf.gz

    > ./test/test_merge2.vcf

生成的结果如下：

##fileformat=VCFv4.2

##bcftools_mergeVersion=0.0.1+htslib-0.0.1

##bcftools_mergeCommand=merge ./test/test1.vcf.gz ./test/test2.vcf.gz

#CHROM      POS         ID    REF ALT QUAL     FILTER   INFO       FORMAT        Sample01      Sample02

1       6324        .        T      C,A 45    PASS       .        GT:AD    1/1:1,22   2/2:0,23

1       16324      .        T      A,C,G      40    PASS       .        GT:AD    0/1:10,11         2/3:0,12,12

4       134861    .        G      C      50    PASS       .        GT:AD    0/1:7,13   ./.:.

4       134861    .        G      GC   15    LowQual        .        GT:AD    ./.:.   1/1:1,19

可以看到最新的bcftools输出的vcf已经更新到4.2版本了，然后默认情况下对于SNP和Indel重合的位点的处理和CombineVariants是一样的(另外可以通过--merge参数来修改)，然后这边的Quality值变成了两个的最大值，而AD值也同样是没更新的，如果存在INFO部分的话，很多值的合并方式也可以通过--info-rules来修改。

最后补充两点：

1) 这边讲的合并指的是多个single-sample的vcf文件合并成单个multi-samples的vcf文件，当然其中有一些工具也支持合并同样sample(s)但是包含不同位置结果的vcf文件(例如开始Call variants的时候是以每条染色体为单位的，后面可以通过CombineVariants将所有结果合并到一起);

2) 这边举得很多例子实际数据中可能占得比例并不大，而且很多时候会直接去掉这些位点(例如计算LD的时候只用bi-allelic的位点)，或者有些信息根本用不上(大部分下游软件用到比较多的是GT值)，这边讲的比较细节，大部分情况下对实际数据的影响都是可以忽略的，但是如果涉及到相关方面的一些分析，或者是处理低覆盖的数据的时候，这种细节方面引起的误差甚至是错误可能就不容忽视了。

2. Optional data aggregation step

这一步是通过CombineGVCFs

(http://www./gatk/gatkdocs/org_broadinstitute_sting_gatk_walkers_variantutils_CombineGVCFs.html)

先把上一步生成的所有gVCF文件进行合并，不过这一步主要是针对sample数量比较多的情况(上百个sample)，先合并后面操作会比较方便，所以是可选的。

3. Joint genotyping

这一步就是再用GenotypeGVCFs

(http://www./gatk/gatkdocs/org_broadinstitute_sting_gatk_walkers_variantutils_CombineGVCFs.html)

来生成相应的vcf格式的variants结果。

4. Variant recalibration

最后一步还是跟之前一样的VQSR来对variants进行一系列的校正，得到最终可以使用的vcf文件。