VVP是一个用来给位点打分的软件。
文章PMID:29463208
从Vcf到VVP的注释,要经过下面的过程: Vcf -> VEP注释 -> VVP注释。
首先我们来安装VEP。 需要注意的是,在装VEP之前,要确保已经安装了perl。 然后需要安装下面的perl模组。 可以使用cpan安装。(最好在root用户下进行)
cpan install DBI
cpan install Archive::Zip
cpan install DBD::mysql
cpan install Set::IntervalTree
cpan install JSON
cpan install PerlIO::gzip
cpan install Bio::DB::BigFile
然后才开始装VEP,目前最新版本是92。创建cache目录可以选否(n)。
git clone https://github.com/Ensembl/ensembl-vep.git
cd ensembl-vep
perl INSTALL.pl
有个bug是安装的时候,BioPerl-1.6.924.tar.gz这个鬼东西会因为网络问题下载失败。 多试几次。。。我是装了十几次才装好的。。。
然后需要下载数据库。—> FTP
wget ftp://ftp.ensembl.org/pub/release-92/variation/VEP/homo_sapiens_merged_vep_92_GRCh37.tar.gz
tar -zxvf homo_sapiens_merged_vep_92_GRCh37.tar.gz
接下来,来装VVP。和VEP相比,VVP真的好装的多。 安装GNU Scientific Library和zlib。一般都有。
sudo apt-get install libgsl0ldbl
sudo apt-get install zlib1g
下载并安装VVP。
git clone https://github.com/Yandell-Lab/VVP-pub.git
cd VVP-pub
make
就装好了。 可以运行教程给的例子测试一下。
cd example
../build_background -i 1KG_cftr_background.recode.vep.vcf.gz -o 1KG.build -b 2500 -v CSQ,4,6,1,15
../VVP -i target_spiked_simple.vcf.gz -d 1KG.build -v CSQ,4,6,1,15 1> target.spiked.vvp.out
注意这里,要先使用build_background这个程序来建立background文件。target_spiked_simple.vcf.gz必须是VEP注释过的。 VVP提供了一个基于gnomAD的background文件,点击下载 —-> gnomAD
目前看应该是b37基因组的。下载后解压直接使用。
软件和数据库都下载完了,接下来就是使用了。
首先用VEP注释vcf。
# --offline 是用离线模式
# --cache 使用cache文件
# --dir_cache 指定cache文件的文件夹
# --fasta 参考基因
# --vcf 在vcf后添加注释信息
vep \
-i input.vcf \
-o output.vcf \
--offline \
--assembly GRCh37 \
--cache \
--dir_cache /PATH/to/VEP/ \
--fasta human_g1k_v37_decoy.fasta \
--vcf
要注意如果不加–vcf的话,只是出来就是个单纯的tab分割VEP注释文件。
注释出来的vcf文件,有些的Alt列有两个值(逗号分割),接下来运行VVP的时候会报错。 可以使用下面的python脚本处理一下vcf文件。
inputfile = open('output.vcf', 'r')
outputfile = open('output.final.vcf', 'w')
for line in inputfile:
if line.startswith('#'):
outputfile.write(line)
else:
lineAS = line.split('\t')
chrom = lineAS[0]
POS = lineAS[1]
ID = lineAS[2]
Ref = lineAS[3]
Alt = lineAS[4]
Qual = lineAS[5]
fil = lineAS[6]
info = lineAS[7]
forma = lineAS[8]
RM = lineAS[9]
if Alt.__contains__(','):
alts = Alt.split(',')
for i in alts:
Alt = i
l = [chrom, POS, ID, Ref, Alt, Qual, fil, info, forma, RM]
s = '\t'.join(l)
outputfile.write(s)
else:
l = [chrom, POS, ID, Ref, Alt, Qual, fil, info, forma, RM]
s = '\t'.join(l)
outputfile.write(s)
outputfile.close()
inputfile.close()
得到的文件,可以压缩一下节约空间。
gzip output.final.vcf
最后用VVP进行评价。
VVP -i out.final.vcf.gz -d /PATH/to/gnomad.062717.build -v CSQ,4,6,1,15 1> output.vvp.out
根据文章里的描述,57及以上的得分,都可以认为是可能致病的。 那么我们用下面的命令来输出57分及以上的数据。
awk '$8>56 {print $0}' output.vvp.out > hemi_likelypathogenic.txt
awk '$13>56 {print $0}' output.vvp.out > het_likelypathogenic.txt
awk '$18>56 {print $0}' output.vvp.out > hom_likelypathogenic.txt
完。