gatk3中檢測變異從原來的gatk2的Joint Variant Calling變成了現(xiàn)在的Joint Genotyping，一方面順利的解決了變異檢測過程中的“N+1”難題，另一方面則大大減少了運(yùn)行所需要的時(shí)間和消耗的資源。

在這種模式下只需要對每個(gè)sample的bam文件運(yùn)行HaplotypeCaller，生成每個(gè)sample的gVCF文件，之后再合并這些文件進(jìn)行g(shù)enotyping，這樣原先最耗時(shí)間的步驟(HaplotypeCaller這一步)從與sample數(shù)量的指數(shù)級關(guān)系就變成了線性關(guān)系，而合并genotyping消耗的時(shí)間則相對很少，因此整個(gè)模式帶來的將是一次巨大的革新。

1. Variant calling

針對每個(gè)sample的BAM文件運(yùn)行HaplotypeCaller（每個(gè)sample可以有多個(gè)文件(應(yīng)當(dāng)包含同樣的SM標(biāo)簽)，運(yùn)行的時(shí)候通過-I參數(shù)設(shè)置一起提供給caller），然后每個(gè)sample會(huì)生成一個(gè)gVCFs文件，需要設(shè)置的參數(shù)有：

--emitRefConfidenceGVCF 通過這個(gè)選項(xiàng)來輸出與reference一致的堿基的相關(guān)信息，這樣最終的結(jié)果其實(shí)就相當(dāng)于包含了所有的位點(diǎn)，而不再僅僅是原先的variants位點(diǎn)。

--variant_index_type LINEAR 這個(gè)選項(xiàng)是主程序里面的，具體可以不管

--variant_index_parameter 128000 跟上面的選項(xiàng)相關(guān)的一個(gè)參數(shù)，也不用管

其實(shí)這個(gè)思路還是很明顯的，原先情況下單純的輸出variants位點(diǎn)再合并會(huì)導(dǎo)致多出來的reference位點(diǎn)的實(shí)際情況無法確定，所以把所有位點(diǎn)一起輸出再合并就可以比較了。因此理論上這種做法在原來的框架下也是能做的，因?yàn)楸旧泶蟛糠值腸aller就可以輸出所有位點(diǎn)的信息，但是原來的做法會(huì)有以下幾個(gè)問題:

(1)生成文件的體積會(huì)非常大，很容易就超過原來的bam文件，即使壓縮后也很大，一方面極度浪費(fèi)資源，另一方面下游處理效率會(huì)很差；

(2)通過直接輸出的這些reference位點(diǎn)的信息可能不準(zhǔn)，或者不完全，因?yàn)楫?dāng)初的目的可能更多的是用于debug之類；

(3)很多情況下通過合并得到的multi-sample vcf文件中的結(jié)果部分值不準(zhǔn)確，或者根本得不到更新。

這些也是造成所謂“N+1”問題的根本原因。

第三點(diǎn)其實(shí)跟第二點(diǎn)本質(zhì)上是一樣的，都是由于單個(gè)vcf文件中存儲(chǔ)的每個(gè)sample的信息不全(loss)導(dǎo)致合并時(shí)很多值無法更新，尤其是碰到多variants位點(diǎn)時(shí)這個(gè)問題最明顯。不過解決的思路也很簡單，就是嘗試將每個(gè)sample BAM文件中的有用信息盡可能無損(lossless)的轉(zhuǎn)換到一個(gè)文件中，再通過這個(gè)文件提取得到最后的vcf文件。

這邊我想順便介紹一些平?？梢杂脕磉M(jìn)行這種vcf文件合并操作的一些工具，并且討論一下在使用過程中可能需要注意的一些問題:

a).vcf-merge

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

這個(gè)工具是vcftools里面的一個(gè)Perl腳本，具體用法如下:

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

首先最大的缺點(diǎn)是慢(所以后來重新開發(fā)了C版本，見后)；然后另外一個(gè)比較麻煩的就在于對多variants位點(diǎn)的處理，碰到這些位點(diǎn)的時(shí)候不會(huì)更新各個(gè)sample，例如我們可以嘗試合并下面兩個(gè)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來合并這兩個(gè)文件，具體代碼如下:

## 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生成對應(yīng)的index文件，這樣可以實(shí)現(xiàn)對vcf文件的隨機(jī)讀寫。

得到的結(jié)果如下:

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

從更新的結(jié)果來看，Quality值應(yīng)該是取得兩者的平均值，然后INFO部分增加了AC，AN以及SF的信息，F(xiàn)ILTER的標(biāo)簽是以被FILTER掉的為準(zhǔn)，然后在SF這邊有顯示(1f就表示在第二個(gè)sample中是LowQual的)，假設(shè)我們當(dāng)時(shí)定義的LowQual的標(biāo)準(zhǔn)是30的話，這邊根據(jù)情況可能就要更新成PASS；然后可以看到，這邊AD值(這個(gè)值反應(yīng)的是每種形態(tài)的allele的depth)的信息是沒有更新的，因?yàn)槊總€(gè)caller生成的vcf文件中包含的INFO或者FORMAT部分的標(biāo)簽各種各樣(這邊的AD值就是HaplotypeCaller或者UnifiedGenotyper默認(rèn)輸出的)，除非是配套的下游操作軟件，否則一般是很難考慮到所有情況的，而且這邊即使嘗試更新，得到的結(jié)果也不一定完全準(zhǔn)確，例如即使是0/1的情況下，也是有存在第三種形態(tài)的可能的，只是當(dāng)時(shí)不足以對genotyping產(chǎn)生影響所以被舍棄了而已。

b).GATKCombineVariants

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

這個(gè)是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文件的。

然后結(jié)果如下：

#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這個(gè)值的問題了，所以前兩個(gè)位點(diǎn)直接把AD值給刪掉了(刪掉了……)，后面兩個(gè)因?yàn)槭荌ndel和SNP位點(diǎn)重合，而Indel的實(shí)際位置本身就應(yīng)該是134861+1，所以這邊分成兩個(gè)可能應(yīng)該更準(zhǔn)確一點(diǎn)，但這樣vcf文件當(dāng)中就出現(xiàn)了duplicate的位置。

c). bcftools merge

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

這個(gè)工具是后來重新開發(fā)的，既包含了原來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

生成的結(jié)果如下：

##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已經(jīng)更新到4.2版本了，然后默認(rèn)情況下對于SNP和Indel重合的位點(diǎn)的處理和CombineVariants是一樣的(另外可以通過--merge參數(shù)來修改)，然后這邊的Quality值變成了兩個(gè)的最大值，而AD值也同樣是沒更新的，如果存在INFO部分的話，很多值的合并方式也可以通過--info-rules來修改。

最后補(bǔ)充兩點(diǎn)：

1) 這邊講的合并指的是多個(gè)single-sample的vcf文件合并成單個(gè)multi-samples的vcf文件，當(dāng)然其中有一些工具也支持合并同樣sample(s)但是包含不同位置結(jié)果的vcf文件(例如開始Call variants的時(shí)候是以每條染色體為單位的，后面可以通過CombineVariants將所有結(jié)果合并到一起);

2) 這邊舉得很多例子實(shí)際數(shù)據(jù)中可能占得比例并不大，而且很多時(shí)候會(huì)直接去掉這些位點(diǎn)(例如計(jì)算LD的時(shí)候只用bi-allelic的位點(diǎn))，或者有些信息根本用不上(大部分下游軟件用到比較多的是GT值)，這邊講的比較細(xì)節(jié)，大部分情況下對實(shí)際數(shù)據(jù)的影響都是可以忽略的，但是如果涉及到相關(guān)方面的一些分析，或者是處理低覆蓋的數(shù)據(jù)的時(shí)候，這種細(xì)節(jié)方面引起的誤差甚至是錯(cuò)誤可能就不容忽視了。

2. Optional data aggregation step

這一步是通過CombineGVCFs

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

先把上一步生成的所有gVCF文件進(jìn)行合并，不過這一步主要是針對sample數(shù)量比較多的情況(上百個(gè)sample)，先合并后面操作會(huì)比較方便，所以是可選的。

3. Joint genotyping

這一步就是再用GenotypeGVCFs

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

來生成相應(yīng)的vcf格式的variants結(jié)果。

4. Variant recalibration

最后一步還是跟之前一樣的VQSR來對variants進(jìn)行一系列的校正，得到最終可以使用的vcf文件。