如何对比vcf文件
这篇文章将为大家详细讲解有关如何对比vcf文件，小编觉得挺实用的，因此分享给大家做个参考，希望大家阅读完这篇文章后可以有所收获。
坚守“ 做人真诚 · 做事靠谱 · 口碑至上 · 高效敬业 ”的价值观，专业网站建设服务10余年为成都湿喷机小微创业公司专业提供企业网站建设营销网站建设商城网站建设手机网站建设小程序网站建设网站改版，从内容策划、视觉设计、底层架构、网页布局、功能开发迭代于一体的高端网站建设服务。
如果我们要比较的两个vcf文件的参考基因组版本不一致，就需要使用CrossMap等软件进行参考基因组版本转换，然后里使用 SnpSift 软件的 Concordance 命令比较它们。其中CrossMap软件依赖pyBigWig，使用conda进行安装，代码如下：
conda create -n py3 python=3.6
conda activate py3
conda install -c bioconda pyBigWig
pip3 install CrossMap
进行参考基因组版本转换的命令如下：
# 需要自行下载 hg19ToHg38.over.chain.gz 文件，以及参考基因组 Homo_sapiens_assembly38.fasta
python ~/miniconda3/envs/py3/bin/CrossMap.py \
vcf ~/data/liftover/hg19ToHg38.over.chain.gz test.snp.hg19.vcf \
~/data/Homo_sapiens_assembly38.fasta test.snp.hg38.vcf
可以把snp和indel的vcf文件都转换一下，然后拿到的转换好的文件如下：
1.3M Jul 8 05:16 test.indel.hg38.vcf
23K Jul 8 05:16 test.indel.hg38.vcf.unmap
1003K Jun 19 11:10 test.indel.vcf
13M Jul 8 05:18 test.snp.hg38.vcf
245K Jul 8 05:18 test.snp.hg38.vcf.unmap
13M Jun 19 18:29 test.snp.vcf
可以看到转换的成功率是非常高的！unmap的文件很小，因为确实参考基因组有变化，总有一下基因组片段被修改了。
但是，有意思的是，之前我们的vcf文件是严格按照基因组坐标排好序的，但是转换过后，出现了部分坐标乱序情况，如下：
这个很容易理解，因为同一个物种的不同版本参考基因组肯定是有
chr1 119955031 . G A
chr1 148483282 rs7513869 C T
chr1 144995248 rs6600697 A G
chr1 144995236 rs6600696 A C
chr1 144995050 rs1884147 C T
chr1 144995033 rs1884146 A G
也就是说，人类的参考基因组在由hg19进化到hg38的时候，不仅仅是片段的自然扩充，还包括一些以前组装顺序弄错了的片段的纠正。
这样坐标乱序的vcf文件，在很多下游分析都是不友好的，所以可以使用下面的代码进行简单过滤。
input=test.snps.VQSR.vcf
cat $input | java -jar ~/biosoft/snpEff/SnpSift.jar filter "( DP > 20 & FILTER = 'PASS' )" | \
perl -alne '{print unless $F[0] =~ /_/}' | \
awk '$1 ~ /^#/ {print $0;next} {print $0 | "sort -k1,1 -k2,2n"}' | \
grep -v '1/2' > test.filter.sort.vcf
# 检查不同染色体分布情况：
cat new.filter.sort.vcf |grep -v '^#' |cut -f 1 |sort |uniq
# 接下来就可以对干净的VCF文件进行注释啦
java -jar ~/biosoft/snpEff/snpEff.jar GRCh48.86 \
test.filter.sort.vcf > test.filter.sort.eff.vcf
关于“如何对比vcf文件”这篇文章就分享到这里了，希望以上内容可以对大家有一定的帮助，使各位可以学到更多知识，如果觉得文章不错，请把它分享出去让更多的人看到。
网页名称：如何对比vcf文件
文章位置：http://6mz.cn/article/jdcpcj.html