最近工作中遇到需要比对两个不同的VCF文件，看看两个VCF文件中哪些变异位点是相同的，哪些是不同的。为了寻找解决办法，于是求救Google，发现了一篇非常不错的教程。于是今天和大家整理学习一下如何比较两个VCF文件,并且检验一下不同工具的效果。
准备工作
这里会用到几种不同的工具,都可以通过conda直接来安装:
conda install -y -c bioconda vcftools
conda install -y -c bioconda bcftools
conda install -y -c bioconda bedtools
conda install -y -c bioconda perl-vcftools-vcf
conda install -y -c bioconda tabix
测试文件会用到由GATK提供的文件（只有4.1M大小非常方便测试）。测试文件下载：
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/NA12878.knowledgebase.snapshot.20131119.b37.vcf.gz
使用bedtools，简单来看看改VCF文件含有的SNPs信息:
#检查SNPs的总数
bcftools view -v snps NA12878.knowledgebase.snapshot.20131119.b37.vcf.gz | grep -v "^#" | wc -l
###SNPs总数为281347
#检查只含有单个alternate SNP的位点，换句话说多出来的就是含有两个或者多个alternate变异：
bcftools view -v snps NA12878.knowledgebase.snapshot.20131119.b37.vcf.gz | grep -v "^#" | cut -f2 | sort -u | wc -l
###特异只含单个altern SNPs总数为280764
构建测试VCF文件
讲完基础的部分后回到咱们这文章的主题，我们是要比对两个VCF文件，所以首先要利用上面的测试数据，构建两个“相近的带有不同子集”VCF文件。
这里通过perl为每个变体生成一个随机数(使用种子）来创建可变的子集，并且仅在生成的随机数大于0.5时才打印出变异的位点。这样我们会生成两个VCF文件原始文件SNP变异的一半，其SNP位点在50％内重叠。
###构建第一VCF文件
bcftools view -v snps NA12878.knowledgebase.snapshot.20131119.b37.vcf.gz |
perl -lane 'BEGIN {srand(31)} if (/^#/) { print } elsif (length($F[3]) == 1) { if (rand(1) > 0.5) {print} }' |
bgzip > first.vcf.gz
tabix -p vcf first.vcf.gz
###构建第二个VCF文件
bcftools view -v snps NA12878.knowledgebase.snapshot.20131119.b37.vcf.gz |
perl -lane 'BEGIN {srand(1984)} if (/^#/) { print } elsif (length($F[3]) == 1) { if (rand(1) > 0.5) {print} }' |
bgzip > second.vcf.gz
tabix -p vcf second.vcf.gz
###分别检查两个VCF文件的变异数目
zcat first.vcf.gz| grep -v "^#" | wc -l
# 第一个VCF文件含有140879个变异
zcat second.vcf.gz| grep -v "^#" | wc -l
# 第二个文件含有140367个变异
使用不同工具来执行VCF文件之间的比较:
使用bedtools进行比较:
###先使用intersect function进行查看
###首先以first.vcf.gz文件为基准，看看两个文件重叠的变异位点数目
bedtools intersect -u -a first.vcf.gz -b second.vcf.gz | wc -l
### 结果是70451
### 接着以second.vcf.gz为基准，看看两个文件重叠的变异位点数目
bedtools intersect -u -a second.vcf.gz -b first.vcf.gz | wc -l
### 结果是70453
###这里也可以使用bedtools jaccard function
bedtools jaccard -a first.vcf.gz -b second.vcf.gz
###输出结果为
intersection union-intersection jaccard n_intersections
70372 210680 0.334023 70157
接着可以使用vcf-compare进行比较
使用下面的命令行
vcf-compare first.vcf.gz second.vcf.gz
生成的结果如下:
# This file was generated by vcf-compare.
# The command line was: vcf-compare(v0.1.14-12-gcdb80b8) first.vcf.gz second.vcf.gz
#
#VN 'Venn-Diagram Numbers'. Use `grep ^VN | cut -f 2-` to extract this part.
#VN The columns are:
#VN 1 .. number of sites unique to this particular combination of files
#VN 2- .. combination of files and space-separated number, a fraction of sites in the file
VN 69898 second.vcf.gz (49.8%)
VN 70372 first.vcf.gz (50.0%) second.vcf.gz (50.2%)
VN 70410 first.vcf.gz (50.0%)
#SN Summary Numbers. Use `grep ^SN | cut -f 2-` to extract this part.
SN Number of REF matches: 70371
SN Number of ALT matches: 70213
SN Number of REF mismatches: 1
SN Number of ALT mismatches: 158
SN Number of samples in GT comparison: 0
# Number of sites lost due to grouping (e.g. duplicate sites): lost, %lost, read, reported, file
SN Number of lost sites: 97 0.1% 140879 140782 first.vcf.gz
SN Number of lost sites: 97 0.1% 140367 140270 second.vcf.gz
然后bcftools isec也能够做同样的事情:
bcftools isec first.vcf.gz second.vcf.gz -p first_second
# records private to first.vcf.gz
cat first_second/0000.vcf | grep -v "^#" | wc -l
70529
# records private to second.vcf.gz
cat first_second/0001.vcf | grep -v "^#" | wc -l
70017
# records from first.vcf.gz shared by both first.vcf.gz second.vcf.gz
cat first_second/0002.vcf | grep -v "^#" | wc -l
70350
# records from second.vcf.gz shared by both first.vcf.gz second.vcf.gz
cat first_second/0003.vcf | grep -v "^#" | wc -l
70350
小结
在这篇文章中，我演示了使用三种不同的工具来比较VCF文件。以下是每种工具的简短总结：
- BEDTools可用于比较VCF文件，但只能通过比较基因组坐标进行比较；这可以提供对两个文件中有多少个重叠变异位点的快速解答，并且可以用来计算Jaccard索引，从而指示总体两个文件重叠位点的数量
- vcf-compare提供了BEDTools的其他统计信息，包括重复位点的数量和Venn-Diagram数字，它们显示了每个相应的VCF文件中非它变体的数量
- bcftools isec还提供了Venn-Diagram数字，并根据这些交集另外创建了VCF文件。
值得注意的是，不同工具之间显示的重复位点数字略有不同。我怀疑这可能是由于每种工具处理重复位点时略微有所致。
网友评论