多款軟件進行vcf合并-gatk售淡、vcftools得滤、bcftools

vcf文件儲存的是樣本的變異信息文件勋锤，在同一批次分析中饭玲，如果不是采用joint calling的方式進行分析，最終會獲得單個樣本的變異數(shù)據(jù)叁执。這種文件很難對同組不同樣本進行差異SNP分析茄厘，此處就需要對文件進行合并。vcf文件的合并有很多的軟件可以做谈宛，主要的就是GATK次哈、vcftools和bcftools三種，但是具體的合并方法需要根據(jù)不同vcf文件中的信息來判斷入挣。

1. 樣本相同亿乳、位點獨立的vcf文件合并

樣本相同、位點獨立指的是在不同文件中包含的樣本一致径筏，但是位點批次質(zhì)檢沒有交集葛假，分染色體call變異的結(jié)果就是這一類的典型。這類最簡單的方法可以直接使用cat命令進行合并滋恬，但是未免不太專業(yè)聊训，以下是三款軟件的合并方法。

1.1 gatk：GatherVcfs|MergeVcfs

gatk4提供了兩種合并vcf文件的方法恢氯，分別是GatherVcfs和MergeVcfs带斑，兩個方法都是對相同樣本數(shù)據(jù)集的變異結(jié)果進行合并鼓寺，命令示例如下。

# GatherVcfs
gatk GatherVcfs -I concat-a.vcf -I concat-b.vcf -O combine_a_b_samesample_diffsites.vcf
# MergeVcfs
gatk MergeVcfs -I concat-a.vcf -I concat-b.vcf -O combine_a_b_diffsample_allsites_gatk.vcf

兩個命令執(zhí)行結(jié)果完全一致勋磕，結(jié)果如下妈候。

##fileformat=VCFv4.2
##FILTER=<ID=q10,Description="Quality below 10">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##contig=<ID=1,length=17540695>
##contig=<ID=2,length=14896646>
##contig=<ID=3,length=12399606>
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  A
# concat-a.vcf文件位點
1       100     .       GTTT    G       1806    q10     DP=35   GT:GQ:DP        0/1:409:35
1       110     .       C       T,G     1792    PASS    DP=32   GT:GQ:DP        0/1:245:32
1       120     .       GA      G       628     q10     DP=21   GT:GQ:DP        1/1:21:21
1       130     .       GAA     GG      1016    PASS    DP=22   GT:GQ:DP        0/1:212:22
1       140     .       GT      G       727     PASS    DP=30   GT:GQ:DP        0/1:150:30
1       150     .       TAAAA   TA,T    246     PASS    DP=10   GT:GQ:DP        1/2:12:10
2       100     .       GTTT    G       1806    q10     DP=35   GT:GQ:DP        0/1:409:35
2       110     .       CAAA    C       1792    PASS    DP=32   GT:GQ:DP        0/1:245:32
2       120     .       GA      G       628     q10     DP=21   GT:GQ:DP        1/1:21:21
2       130     .       GAA     G       1016    PASS    DP=22   GT:GQ:DP        0/1:212:22
2       140     .       GT      G       727     PASS    DP=30   GT:GQ:DP        0/1:150:30
2       150     .       TAAAA   TA,T    246     PASS    DP=10   GT:GQ:DP        1/2:12:10
2       160     .       TAAAA   TA,TC,T 246     PASS    DP=10   GT:GQ:DP        0/2:12:10
# concat-b.vcf文件位點
3       241     .       GTTT    G       1806    q10     DP=35   GT:GQ:DP        0/1:409:35
3       251     .       CAAA    C       1792    PASS    DP=32   GT:GQ:DP        0/1:245:32
3       261     .       GA      G       628     q10     DP=21   GT:GQ:DP        1/1:21:21
3       271     .       GAA     G       1016    PASS    DP=22   GT:GQ:DP        0/1:212:22

1.2 vcftools：vcf-concat

vcf-concat并不是vcftools的一個子命令，而是軟件安裝目錄下附帶的功能模塊挂滓，vcftools安裝完成后可以直接使用苦银，命令如下。

/vcftools/bin/vcf-concat concat-a.vcf concat-b.vcf > combine_a_b_samesample_diffsites_vcftools.vc

在進行較大文件合并時赶站，最好將vcf文件進行壓縮并創(chuàng)建索引幔虏，效率會快。如果待合并的vcf文件很多贝椿，可以將文件名寫入到一個文件想括，由參數(shù)-f進行操作。該命令合并完的vcf位點變異信息與gatk結(jié)果一致烙博，只不過表頭信息順序會發(fā)生改變瑟蜈，不影響數(shù)據(jù)使用。

1.3 bcftools：concat

bcftools軟件在處理vcf文件時渣窜，某些情況下會優(yōu)于vcftools踪栋，合并vcf文件的命令如下。

bcftools concat concat-a.vcf concat-b.vcf -o combine_a_b_samesample_diffsites_bcftools.vcf

合并之后的vcf文件位點信息與其余兩款軟件處理結(jié)果一致图毕，只不過在輸出結(jié)果的header中會出現(xiàn)運行的命令，示例如下眷唉。

##bcftools_concatVersion=1.3.1+htslib-1.3.1
##bcftools_concatCommand=concat -o combine_a_b_samesample_diffsites_bcftools.vcf concat-a.vcf concat-b.vcf

2. 樣本不同予颤，位點相同或不同

這種合并方式主要是對不同樣本的變異文件進行合并，合并時共有位點會進行合并統(tǒng)計冬阳，非共有位點若在某一個樣本中沒有變異蛤虐，則會自動記為缺失

2.1 vcftools:vcf-merge

模塊名稱為vcf-merge，在進行merge操作時肝陪，會對文件中的位點進行重排驳庭，耗時較長。輸入文件需要壓縮后創(chuàng)建索引氯窍，示例命令如下饲常。

bgzip merge-test-a.vcf.gz && tabix merge-test-a.vcf.gz
bgzip merge-test-b.vcf.gz && tabix merge-test-b.vcf.gz
/vcftools/bin/vcf-merge merge-test-a.vcf.gz merge-test-b.vcf.gz > combine_a_b_diffsamples_allsites_vcftools.vcf

合并之后會對不同文件中的數(shù)據(jù)集進行整合，沒有變異的位點會自動標(biāo)記為缺失狼讨。

2.2 bcftools:merge

使用的方法為merge贝淤，示例如下。

bcftools merge merge-test-a.vcf.gz merge-test-b.vcf.gz -o combine_a_b_diffsamples_allsites_bcftools.vcf

該方法也需要實現(xiàn)對所有vcf文件進行壓縮并創(chuàng)建索引政供，否則程序無法運行播聪。

需要注意的就是朽基，merge合并之后，不同軟件生成的結(jié)果會存在很大的差異离陶，主要是統(tǒng)計結(jié)果的重新計算上稼虎，示例如下。

# merge-test-a.vcf
1   3184885 .   TAAAA   TA,T    246 PASS    DP=10   GT:GQ:DP    1/2:12:10
# merge-test-b.vcf
1   3184885 .   TAAA    T   598 PASS    DP=16   GT:GQ:DP    0/1:435:16
# combine_a_b_diffsamples_allsites_vcftools.vcf
1   3184885 .   TAAAA   TA,T    422.00  PASS    AC=2,1;AN=4;DP=26;SF=0,1    GT:DP:GQ    1/2:10:12   0/1:16:435
# combine_a_b_diffsamples_allsites_bcftools.vcf
1   3184885 .   TAAAA   TA,T    598 PASS    DP=26   GT:GQ:DP    1/2:12:10   0/1:435:16

gatk:CombineVariants

gatk3提供了一個CombineVariants可以進行變異數(shù)據(jù)的合并招刨，而gatk4中并沒有找到功能相同的模塊霎俩。CombineVariants使用如下。

java -jar /GenomeAnalysisTK-3.8/GenomeAnalysisTK.jar -T CombineVariants -V merge-test-a.vcf -V merge-test-b.vcf -o combine_a_b_diffsample_allsites_gatk.vcf -R ref.fna

合并結(jié)果示例如下计济。

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  A       B
1       3062915 .       GTTT    G,GT    1806    q10;q20 AC=1,1;AF=0.250,0.250;AN=4;DP=49;set=FilteredInAll      GT:DP:GQ        0/1:35:99       0/2:14:99
1       3106154 .       CAAAA   CA,C    1792    PASS    AC=1,1;AF=0.250,0.250;AN=4;DP=47;set=Intersection       GT:DP:GQ        0/1:32:99       0/2:15:99
1       3157410 .       GA      G       628     PASS    AC=3;AF=0.750;AN=4;DP=32;set=filterInvariant-variant2   GT:DP:GQ        1/1:21:21       0/1:11:49
1       3162006 .       GAA     G       1016    PASS    AC=2;AF=0.500;AN=4;DP=41;set=Intersection       GT:DP:GQ        0/1:22:99       0/1:19:99
1       3177144 .       GT      G       727     PASS    AC=2;AF=0.500;AN=4;DP=54;set=Intersection       GT:DP:GQ        0/1:30:99       0/1:24:99
1       3184885 .       TAAAA   TA,T    246     PASS    AC=2,1;AF=0.500,0.250;AN=4;DP=26;set=Intersection       GT:DP:GQ        1/2:10:12       0/1:16:99
2       3188209 .       GA      G       162     .       AC=1;AF=0.500;AN=2;DP=15;set=variant2   GT:DP:GQ        ./.     0/1:15:99
2       3199812 .       G       GTT,GT  481     PASS    AC=1,1;AF=0.500,0.500;AN=2;DP=26;set=variant    GT:DP:GQ        1/2:26:99       ./.
3       3199812 .       G       GTT,GT  353     PASS    AC=1,1;AF=0.500,0.500;AN=2;DP=19;set=variant2   GT:DP:GQ        ./.     1/2:19:99
3       3199815 .       G       A       353     PASS    AC=1;AF=0.500;AN=2;DP=19;set=variant2   GT:DP:GQ        ./.     0/1:19:99
3       3212016 .       CTT     C,CT    565     PASS    AC=1,1;AF=0.500,0.500;AN=2;DP=26;set=variant    GT:DP:GQ        1/2:26:91       ./.
4       3212016 .       CTT     C       677     q20     AC=1;AF=0.500;AN=2;DP=15;set=FilteredInAll      GT:DP:GQ        ./.     0/1:15:99
4       3258448 .       TACACACAC       T       325     PASS    AC=1;AF=0.500;AN=2;DP=31;set=variant    GT:DP:GQ        0/1:31:99       ./.

從上面的實例可以看出茸苇，在合并計算分型質(zhì)量時，vcftools會對統(tǒng)計結(jié)果進行平均取值沦寂，bcftools則取其中的最大值輸出学密，而gatk會取最小值輸出，并且gatk和vcftools在輸出原有數(shù)據(jù)基礎(chǔ)之上传藏，還會重新計算AC腻暮、AN等指標(biāo)。

對于vcf數(shù)據(jù)合并毯侦，由于不同樣本發(fā)生的變異位置很難保證一致哭靖，所以在單獨合并不同樣本的數(shù)據(jù)時，合并后的結(jié)果往往具有很高的缺失率侈离，后續(xù)的差異分型實質(zhì)上只是對不同樣本的共有位點進行分析试幽，丟失了某些樣本或群體的特異位點。為了避免這種情況卦碾，多樣本差異分析的項目最好是用joint calling進行變異檢測铺坞，能獲得更多的位點。

3. 參考資料

[1] https://gatk.broadinstitute.org/hc/en-us/articles/360046221931-GatherVcfs-Picard-
[2] http://samtools.github.io/bcftools/bcftools.html#concat
[3] https://gatk.broadinstitute.org/hc/en-us/articles/360045800732-MergeVcfs-Picard-
[4] https://vcftools.github.io/man_latest.html

多款軟件進行vcf合并-gatk、vcftools漩怎、bcftools

多款軟件進行vcf合并-gatk勋颖、vcftools、bcftools