用VVP分析一下致病位点看看

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

完。