生物信息文件夹

接近肿瘤或WES的VEP注释的完美状态？

  • 2024年 04月02日

太久没有更新了，也不知道写点什么。最近继续读了下VEP的文档，基本上单独用VEP可以取代以前snpeff+annovar的方案了。

以前也写过一篇VEP的安装和使用，这篇是更新版本。

以下均是基于GRCh37来写的。

安装

建议使用docker安装

docker pull ensemblorg/ensembl-vep:release_111.0

在集群中，一般使用singularity（已经改名叫apptainer了）

singularity pull ensembl-vep.111.sif docker://ensemblorg/ensembl-vep:release_111.0

数据库

下载VEP官方提供的cache，一般采用refseq的转录本编号，因此我只下了refseq版本

curl -O https://ftp.ensembl.org/pub/release-111/variation/indexed_vep_cache/homo_sapiens_refseq_vep_111_GRCh37.tar.gz
tar xzf homo_sapiens_refseq_vep_111_GRCh37.tar.gz

VEP现在的最新版本是Version 111，GRCh37的cache（即提供下载的数据库）内置了以下内容，意味着，如果只需要注释下面的内容的话，就不用去下载其他数据库了。但是，容易看到，GRCh37内置的数据库均比较旧，因此对于部分有用的字段，就要自行下载新的数据库了。

Source Version (GRCh37)
Ensembl database version 111
Genome assembly GRCh37.p13
GENCODE 19
RefSeq 105.20220307
(GCF_000001405.25_GRCh37.p13_genomic.gff)
Regulatory build 1.0
PolyPhen 2.2.2
SIFT 5.2.2
dbSNP 156
COSMIC 98
HGMD-PUBLIC 2020.4
ClinVar 2023-06
1000 Genomes Phase 3
gnomAD exomes r2.1, exomes only
gnomAD genomes -

以下是一般都会用到的数据库。

Clinvar

VEP 111版本内置的clinvar是2023年6月的版本，可以从NCBI官方下载Clinvar数据库，然后作为一个custom库在注释时导入

curl -O https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/clinvar.vcf.gz
curl -O https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/clinvar.vcf.gz.tbi

dbNSFP

这个库中整合了非常多的有害性预测评分，基本上用了这个其他的有害性评估都不用了，需要同时下载插件和数据库。

# 插件
wget https://github.com/Ensembl/VEP_plugins/raw/release/111/dbNSFP.pm

# 数据库
version=4.7a
wget ftp://dbnsfp:dbnsfp@dbnsfp.softgenetics.com/dbNSFP${version}.zip
unzip dbNSFP${version}.zip
zcat dbNSFP${version}_variant.chr1.gz | head -n1 > h
zgrep -h -v ^#chr dbNSFP${version}_variant.chr* | awk '$8 != "." ' | sort -k8,8 -k9,9n - | cat h - | bgzip -c > dbNSFP${version}_grch37.gz
tabix -s 8 -b 9 -e 9 dbNSFP${version}_grch37.gz

spliceAI

illumina的有害性剪接预测软件，预先对大量突变进行了评估后获得的数据库，需要同时下载插件和数据库。

# 插件
wget https://github.com/Ensembl/VEP_plugins/raw/release/111/SpliceAI.pm

注意这个数据库需要在Illumina的basespace中注册账户，登录，然后才能下载。

# 登录此网站
https://basespace.illumina.com/s/otSPW8hnhaZR

# 下载数据后进行处理
tabix -p vcf spliceai_scores.raw.snv.hg37.vcf.gz
tabix -p vcf spliceai_scores.raw.indel.hg37.vcf.gz

cytoband

VEP中貌似没有提供cytoband注释，可以作为一个custom库导入。

# 下载数据库
wget https://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/cytoBand.txt.gz
gunzip cytoBand.txt.gz
sed -i 's/chr//' cytoBand.txt
grep -v "#" cytoBand.txt | sort -k1,1 -k2,2n -k3,3n -t$'\t' | bgzip -c > cytoBand.bed.gz
tabix -p bed cytoBand.bed.gz

intervar

没错，这个intervar是从annovar那里下载的intervar_20180118版本数据库，而不是重分析的方案，作为一个custom库导入（何尝不是一种NTR）。

# 下载后调整格式
perl annotate_variation.pl --downdb --webfrom annovar --buildver hg19 intervar_20180118 /path/to/vep/database
awk 'BEGIN{FS=OFS="\t"} {gsub(/ /,"_",$6); print $1,$2,$3,$4,$5,".",".","Intervar="$6}' hg19_intervar_20180118.txt > intervar.vcf
sed -i '1d' intervar.vcf

需要制作一个vcf的header.txt文件，可以从现成的vcf里创建，也可以自行构建，一般需要下面这些内容

##fileformat=VCFv4.2
##reference=GRCh37
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO

然后

cat header.txt intervar.vcf | bgzip -c > intervar.vcf.gz
tabix -p vcf intervar.vcf.gz

注释命令

最后是注释命令，把大部分需求的内容都注释上，形成一个新的vcf，只是结果保存在CSQ的tag中。也可以将–vcf参数替换为–tab，然后输出一个tab分割的文件，但是没有vcf灵活。

# 如果使用singularity运行，注意记得把数据库bind进去容器中
singularity exec --bind /path/to/vep_data:/path/to/vep_data ensembl-vep.111.sif <command>

# 单独的VEP命令则是，注意几个需要填路径的点
vep_database=/data/.vep
plugin_dir=/data/.vep/plugin
threads=16

cache_str="Uploaded_variation,Location,REF_ALLELE,Allele,Consequence,IMPACT,DOMAINS,Feature,EXON,INTRON,SYMBOL,HGNC_ID,HGVSc,HGVSp,HGVSg,MAX_AF,Protein_position,Amino_acids,Codons,PUBMED,Existing_variation"
custom_str="cytoBand,ClinVar_CLNSIG,ClinVar_CLNREVSTAT,ClinVar_CLNDN,InterVar_Intervar"
dbnsfp_str="rs_dbSNP,gnomAD_exomes_AF,gnomAD_exomes_EAS_AF,gnomAD_exomes_POPMAX_AF,gnomAD_exomes_AC,gnomAD_exomes_nhomalt,ExAC_AC,ExAC_EAS_AF,1000Gp3_AF,1000Gp3_EAS_AF,REVEL_score,REVEL_rankscore,M-CAP_score,M-CAP_rankscore,M-CAP_pred,GERP++_RS,GERP++_RS_rankscore,MVP_score,MVP_rankscore,phyloP100way_vertebrate,phyloP100way_vertebrate_rankscore,phyloP470way_mammalian,phyloP470way_mammalian_rankscore"
spliceai_str="SpliceAI_pred_DP_AG,SpliceAI_pred_DP_AL,SpliceAI_pred_DP_DG,SpliceAI_pred_DP_DL,SpliceAI_pred_DS_AG,SpliceAI_pred_DS_AL,SpliceAI_pred_DS_DG,SpliceAI_pred_DS_DL"
vep --offline --cache --dir_cache ${vep_database} --refseq \
    --dir_plugins ${plugin_dir} \
    --force_overwrite --fork ${threads} \
    -i input.vcf -o output.vcf \
    --format vcf --vcf --fa human_g1k_v37_decoy.fasta \
    --shift_3prime 1 --assembly GRCh37 --no_escape --check_existing -exclude_predicted --uploaded_allele --show_ref_allele --numbers --domains --pick_allele \
    --total_length --hgvs --hgvsg --symbol --ccds --uniprot --max_af --pubmed \
    --custom file=${vep_database}/clinvar/clinvar.vcf.gz,short_name=ClinVar,format=vcf,type=exact,coords=0,fields=CLNSIG%CLNREVSTAT%CLNDN \
    --custom file=${vep_database}/intervar/intervar.vcf.gz,short_name=InterVar,format=vcf,type=exact,coords=0,fields=Intervar \
    --custom file=${vep_database}/cytoband/cytoBand.bed.gz,short_name=cytoBand,format=bed,type=overlap,coords=0 \
    --plugin dbNSFP,${vep_database}/dbNSFP/dbNSFP4.7a_grch37.gz,"${dbnsfp_str}" \
    --plugin SpliceAI,snv=${vep_database}/spliceAI/spliceai_scores.raw.snv.hg19.vcf.gz,indel=${vep_database}/spliceAI/spliceai_scores.raw.indel.hg19.vcf.gz,cutoff=0.5 \
    --fields "${cache_str},${custom_str},${dbnsfp_str},${spliceai_str}"
搜索