2023年12月7日发(作者：)
多款软件进行vcf合并-gatk、vcftools、bcftools
vcf文件储存的是样本的变异信息文件，在同一批次分析中，如果不是采用joint calling的方式进行分析，最终会获得单个样本的变异数据。这种
文件很难对同组不同样本进行差异SNP分析，此处就需要对文件进行合并。vcf文件的合并有很多的软件可以做，主要的就是GATK、vcftools和
bcftools三种，但是具体的合并方法需要根据不同vcf文件中的信息来判断。
1. 样本相同、位点独立的vcf文件合并
样本相同、位点独立指的是在不同文件中包含的样本一致，但是位点批次质检没有交集，分染色体call变异的结果就是这一类的典型。这类最简单
的方法可以直接使用cat命令进行合并，但是未免不太专业，以下是三款软件的合并方法。
1.1 gatk：GatherVcfs|MergeVcfs
gatk4提供了两种合并vcf文件的方法，分别是GatherVcfs和MergeVcfs，两个方法都是对相同样本数据集的变异结果进行合并，命令示例如
下。
# GatherVcfs
gatk GatherVcfs -I -I -O combine_a_b_samesample_
# MergeVcfs
gatk MergeVcfs -I -I -O combine_a_b_diffsample_allsites_
两个命令执行结果完全一致，结果如下。
##fileformat=VCFv4.2
##FILTER=
##FORMAT=
##FORMAT=
##FORMAT=
##INFO=
##contig=
##contig=
##contig=
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT A
# 文件位点
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
# 文件位点
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 > combine_a_b_samesample_diffsites_
在进行较大文件合并时，最好将vcf文件进行压缩并创建索引，效率会快。如果待合并的vcf文件很多，可以将文件名写入到一个文件，由参数
-
f
进行操作。该命令合并完的vcf位点变异信息与gatk结果一致，只不过表头信息顺序会发生改变，不影响数据使用。
1.3 bcftools：concat
bcftools软件在处理vcf文件时，某些情况下会优于vcftools，合并vcf文件的命令如下。
bcftools concat -o combine_a_b_samesample_diffsites_
合并之后的vcf文件位点信息与其余两款软件处理结果一致，只不过在输出结果的header中会出现运行的命令，示例如下。
##bcftools_concatVersion=1.3.1+htslib-1.3.1
##bcftools_concatCommand=concat -o combine_a_b_samesample_diffsites_
2. 样本不同，位点相同或不同
这种合并方式主要是对不同样本的变异文件进行合并，合并时共有位点会进行合并统计，非共有位点若在某一个样本中没有变异，则会自动记为缺
失
2.1 vcftools:vcf-merge
模块名称为vcf-merge，在进行merge操作时，会对文件中的位点进行重排，耗时较长。输入文件需要压缩后创建索引，示例命令如下。
bgzip && tabix
bgzip && tabix
/vcftools/bin/vcf-merge > combine_a_b_diffsamples_allsites_
合并之后会对不同文件中的数据集进行整合，没有变异的位点会自动标记为缺失。
2.2 bcftools:merge
使用的方法为merge，示例如下。
bcftools merge -o combine_a_b_diffsamples_allsites_
该方法也需要实现对所有vcf文件进行压缩并创建索引，否则程序无法运行。
需要注意的就是，merge合并之后，不同软件生成的结果会存在很大的差异，主要是统计结果的重新计算上，示例如下。
#
1 3184885 . TAAAA TA,T 246 PASS DP=10 GT:GQ:DP 1/2:12:10
#
1 3184885 . TAAA T 598 PASS DP=16 GT:GQ:DP 0/1:435:16
# combine_a_b_diffsamples_allsites_
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_
1 3184885 . TAAAA TA,T 598 PASS DP=26 GT:GQ:DP 1/2:12:10 0/1:435:16
gatk:CombineVariants
gatk3提供了一个CombineVariants可以进行变异数据的合并，而gatk4中并没有找到功能相同的模块。CombineVariants使用如下。
java -jar /GenomeAnalysisTK-3.8/ -T CombineVariants -V -V -o combine_a_b_diffsample_allsites_ -R r
合并结果示例如下。
#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 ./.
从上面的实例可以看出，在合并计算分型质量时，vcftools会对统计结果进行平均取值，bcftools则取其中的最大值输出，而gatk会取最小值输
出，并且gatk和vcftools在输出原有数据基础之上，还会重新计算AC、AN等指标。
对于vcf数据合并，由于不同样本发生的变异位置很难保证一致，所以在单独合并不同样本的数据时，合并后的结果往往具有很高的缺失率，后续
的差异分型实质上只是对不同样本的共有位点进行分析，丢失了某些样本或群体的特异位点。为了避免这种情况，多样本差异分析的项目最好是用
joint calling进行变异检测，能获得更多的位点。
3. 参考资料
发布评论