太久没有更新了，也不知道写点什么。最近继续读了下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版本，当然，也可以下载merged版本（不能指定–refseq），在注释时如果需要优先注释refseq的转录本，则需要通过–pick_order命令调整优先度。
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,DISTANCE,EXON,INTRON,SYMBOL,STRAND,HGNC_ID,HGVSc,HGVSp,HGVSg,MAX_AF,Protein_position,Amino_acids,Codons,PUBMED,Existing_variation"
custom_str="cytoBand,ClinVar_CLNSIG,ClinVar_CLNREVSTAT,ClinVar_CLNDN,ClinVar_CLNHGVS,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 --pick_order mane_plus_clinical,mane_select,refseq,canonical,ensembl --pick_allele \ # 注释所有的转录本需删除这行
--transcript_filter "stable_id match N[MR]_" \ # 这行可以在注释merged库时只保留refseq的转录本
--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%CLNHGVS \
--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}"
上述命令使用了pick_order来根据优先度选择转录本，但在实际测试中，部分基因的转录本仍然未能将Mane Select等提前，具体原因未明（个人感觉是没有对refseq做优先）。这里建议是注释时不要加入与pick相关的参数，把所有转录本都注释出来，然后在下游自行构建一个过滤脚本，传递进需求的转录本列表用于过滤。