最近工作中遇到需要比對兩個不同的VCF文件，看看兩個VCF文件中哪些變異位點是相同的，哪些是不同的。為了尋找解決辦法，于是求救Google，發(fā)現(xiàn)了一篇非常不錯的教程。于是今天和大家整理學習一下如何比較兩個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的總數(shù)
bcftools view -v snps NA12878.knowledgebase.snapshot.20131119.b37.vcf.gz | grep -v "^#" | wc -l
###SNPs總數(shù)為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總數(shù)為280764
構(gòu)建測試VCF文件
講完基礎(chǔ)的部分后回到咱們這文章的主題，我們是要比對兩個VCF文件，所以首先要利用上面的測試數(shù)據(jù)，構(gòu)建兩個“相近的帶有不同子集”VCF文件。
這里通過perl為每個變體生成一個隨機數(shù)(使用種子）來創(chuàng)建可變的子集，并且僅在生成的隨機數(shù)大于0.5時才打印出變異的位點。這樣我們會生成兩個VCF文件原始文件SNP變異的一半，其SNP位點在50％內(nèi)重疊。
###構(gòu)建第一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
###構(gòu)建第二個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文件的變異數(shù)目
zcat first.vcf.gz| grep -v "^#" | wc -l
# 第一個VCF文件含有140879個變異
zcat second.vcf.gz| grep -v "^#" | wc -l
# 第二個文件含有140367個變異
使用不同工具來執(zhí)行VCF文件之間的比較:
使用bedtools進行比較:
###先使用intersect function進行查看
###首先以first.vcf.gz文件為基準，看看兩個文件重疊的變異位點數(shù)目
bedtools intersect -u -a first.vcf.gz -b second.vcf.gz | wc -l
### 結(jié)果是70451
### 接著以second.vcf.gz為基準，看看兩個文件重疊的變異位點數(shù)目
bedtools intersect -u -a second.vcf.gz -b first.vcf.gz | wc -l
### 結(jié)果是70453
###這里也可以使用bedtools jaccard function
bedtools jaccard -a first.vcf.gz -b second.vcf.gz
###輸出結(jié)果為
intersection union-intersection jaccard n_intersections
70372 210680 0.334023 70157
接著可以使用vcf-compare進行比較
使用下面的命令行
vcf-compare first.vcf.gz second.vcf.gz
生成的結(jié)果如下:
# 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
小結(jié)
在這篇文章中，我演示了使用三種不同的工具來比較VCF文件。以下是每種工具的簡短總結(jié)：
- BEDTools可用于比較VCF文件，但只能通過比較基因組坐標進行比較；這可以提供對兩個文件中有多少個重疊變異位點的快速解答，并且可以用來計算Jaccard索引，從而指示總體兩個文件重疊位點的數(shù)量
- vcf-compare提供了BEDTools的其他統(tǒng)計信息，包括重復(fù)位點的數(shù)量和Venn-Diagram數(shù)字，它們顯示了每個相應(yīng)的VCF文件中非它變體的數(shù)量
- bcftools isec還提供了Venn-Diagram數(shù)字，并根據(jù)這些交集另外創(chuàng)建了VCF文件。
值得注意的是，不同工具之間顯示的重復(fù)位點數(shù)字略有不同。我懷疑這可能是由于每種工具處理重復(fù)位點時略微有所致。
