【转】vcf文件合并操作
欢迎关注公众号：oddxix
最近在做PCA，需要将多个样本vcf合并起来，并且加入千人里中国人的数据，收集到下面的资料。
bcftools extract specified samples in vcf format
1、下载安装bcftools。
2、准备样本ID文件，这里命名为samplelistname.txt，一个样本一行，如下所示：
sample1
sample2
sample3
3、输入命令：
bcftools view -S samplelistname.txt /1000genomes/ALL.chr16.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz -Ov > samplelist_1000Genomes.vcf
参考链接：
https://www.biostars.org/p/184950/
https://samtools.github.io/bcftools/bcftools.html#view
原文：http://www.cnblogs.com/chenwenyan/p/9212170.html
单个vcf文件中存储的每个sample的信息不全(loss)导致合并时很多值无法更新，尤其是碰到多variants位点时这个问题最明显。通过尝试将每个sample BAM文件中的有用信息尽可能无损(lossless)的转换到一个文件中，再通过这个文件提取得到最后的vcf文件。
a).vcf-merge**
http://vcftools.sourceforge.net/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)
bgzip ./test/test1.vcf
tabix -p vcf ./test/test1.vcf.gz
bgzip ./test/test2.vcf
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).GATK CombineVariants**
这个是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.github.io/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值)，这边讲的比较细节，大部分情况下对实际数据的影响都是可以忽略的，但是如果涉及到相关方面的一些分析，或者是处理低覆盖的数据的时候，这种细节方面引起的误差甚至是错误可能就不容忽视了。
