2023年12月7日发(作者：)

你知道如何比较两个VCF文件吗？

最近工作中遇到需要比对两个不同的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大小非常方便测试）。测试文件下载：

wgetftp://**************************************.org/bundle/b37/

使用bedtools，简单来看看改VCF文件含有的SNPs信息:

#检查SNPs的总数

bcftools view -v snps | grep -v "^#" | wc -l

###SNPs总数为281347

#检查只含有单个alternate SNP的位点，换句话说多出来的就是含有两个或者多个alternate变异：

bcftools view -v snps | 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 |

perl -lane 'BEGIN {srand(31)} if (/^#/) { print } elsif (length($F[3]) == 1) { if (rand(1) > 0.5) {print} }' |

bgzip >

tabix -p vcf

###构建第二个VCF文件

bcftools view -v snps |

perl -lane 'BEGIN {srand(1984)} if (/^#/) { print } elsif (length($F[3]) == 1) { if (rand(1) > 0.5) {print} }' |

bgzip >

tabix -p vcf

###分别检查两个VCF文件的变异数目

zcat | grep -v "^#" | wc -l

# 第一个VCF文件含有140879个变异

zcat | grep -v "^#" | wc -l

# 第二个文件含有140367个变异

使用不同工具来执行VCF文件之间的比较:

使用bedtools进行比较:

###先使用intersect function进行查看

###首先以文件为基准，看看两个文件重叠的变异位点数目

bedtools intersect -u -a -b | wc -l

### 结果是70451

### 接着以为基准，看看两个文件重叠的变异位点数目

bedtools intersect -u -a -b | wc -l

### 结果是70453

###这里也可以使用bedtools jaccard function

bedtools jaccard -a -b

###输出结果为

intersection union-intersection jaccard n_intersections

70372 210680 0.334023 70157

接着可以使用vcf-compare进行比较

使用下面的命令行

vcf-compare

生成的结果如下:

# This file was generated by vcf-compare.

# The command line was: vcf-compare(v0.1.14-12-gcdb80b8)

#

#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 (49.8%)

VN 70372 (50.0%) (50.2%)

VN 70410 (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

SN Number of lost sites: 97 0.1% 140367 140270

然后bcftools isec也能够做同样的事情:

bcftools isec -p first_second

# records private to

cat first_second/ | grep -v "^#" | wc -l

70529

# records private to

cat first_second/ | grep -v "^#" | wc -l

70017

# records from shared by both

cat first_second/ | grep -v "^#" | wc -l

70350

# records from shared by both

cat first_second/ | grep -v "^#" | wc -l

70350

小结

在这篇文章中，我演示了使用三种不同的工具来比较VCF文件。以下是每种工具的简短总结：

BEDTools可用于比较VCF文件，但只能通过比较基因组坐标进行比较；这可以提供对两个文件中有多少个重叠变异位点的快速解答，并且可

以用来计算Jaccard索引，从而指示总体两个文件重叠位点的数量

vcf-compare提供了BEDTools的其他统计信息，包括重复位点的数量和Venn-Diagram数字，它们显示了每个相应的VCF文件中非它变体的

数量

bcftools isec还提供了Venn-Diagram数字，并根据这些交集另外创建了VCF文件。

值得注意的是，不同工具之间显示的重复位点数字略有不同。我怀疑这可能是由于每种工具处理重复位点时略微有所致。