基因注释软件-VEP

发表于 | 分类于 ，

简介

在TCGA等大型项目中，也推荐使用了VEP作为注释软件，同时，TCGA项目提供了由VEP向MAF转换的工具 vcf2maf，该工具由MSK长期跟进及维护。VEP 是Ensemble旗下开发的软件，有众多的专业人员进行着更新维护，同时他也符合 CLIA

Variant Effect Predictor（VEP）是功能强大的注释、分析工具。它可以对二代测试产生的不同类型变异进行注释，包含SNPs, insertions, deletions, copy number variants和structural variants。也可以依据各种数据库的内容，根据需要，对变异进行过滤和排序。

VEP 输入格式
VEP 结果格式
VEP Calculated variant consequences

软件问题

VEP 本身的对齐，只可以进行全碱基的对齐；

原始变异 ：
GAATCATCATCATCATG
在vep中 ：
GAATCATCATCATCATG
而不会拆分碱基进行对齐
GAATCATCATCATCATG

安装

环境依赖

1
2
3
4
5
6
conda install -c bioconda  perl-DBI
conda install -c bioconda  perl-Archive-Zip
conda install -c bioconda  perl-Bio-DB-HTS-Tabix  # 使用gff进行注释时，该包必须安装
perl -e 'use DBI ; use Archive::Zip ; use Bio::DB:HTS::Tabix'

conda install -c bioconda  perl-DBD-mysql # 使用离线模式，可以跳过该包

VEP的安装可以直接参考官方说明

通过conda安装

VEP本身可以直接通过conda进行安装，由于官方安装文档未提及，在此进行补充说明

1
conda install -c bioconda ensembl-vep==108 # 如果需要安装新版本需要指定，本文撰写阶段最新版本108，但是默认安装版本是v92.

库文件配置

在线版进行数据库下载的速度非常慢，而conda也只会进行vep注释软件的安装，不包含相关库文件，因此需要单独准备库文件，速度会更快。

数据集下载

  1. 运行VEP软件自带的 INSTALL.pl 进行下载，非常之慢，具体方式可以参考官方安装说法说明。
  2. 从官方ftp进行下载
    cache文件可以访问官方ftp进行下载
    https://ftp.ensembl.org/pub/grch37/release-108/variation/vep/homo_sapiens_refseq_vep_108_GRCh37.tar.gz
    参考基因组的fasta序列，
    https://ftp.ensembl.org/pub/grch37/release-108/fasta/homo_sapiens/dna/Homo_sapiens.GRCh37.dna_sm.primary_assembly.fa.gz

扩展功能配置

VEP本身提供了一个外部扩展功能模块 Plugins
可以使用INSTALL.pl 进行安装，也可以直接通过Github仓库进行下载

相关问题及处理记录

20221104 安装VEP108时会发发现存在版本兼容性错误因为bioconda中

1
2
3
# 报错信息
# Compress::Raw::Zlib version 2.201 required--this is only version 2.105
conda install -c conda-forge perl-compress-raw-zlib  #

使用

使用自己的转录本库（transcript models）

VEP本身自带对应的转录本库版本，但是如果遇到默认转录本中缺失了自己需要的转录本信息时，可以构建自己的转录本库，进行注释。参考官方文档：自建注释数据库的参考文档

VEP can use a variety of annotation sources to retrieve the transcript models used to predict consequence types.

  • Cache - a downloadable file containing all transcript models, regulatory features and variant data for a species
  • GFF or GTF - use transcript models defined in a tabix-indexed GFF or GTF file
  • Database - connect to a MySQL database server hosting Ensembl databases

    GFF/GTF files

    VEP can use transcript annotations defined in GFF or GTF files. The files must be bgzipped and indexed with tabix and a FASTA file containing the genomic sequence is required in order to generate transcript models. 
1
2
3
grep -v "#" data.gff | sort -k1,1 -k4,4n -k5,5n -t$'\t' | bgzip -c > data.gff.gz
tabix -p gff data.gff.gz
./vep -i input.vcf --gff data.gff.gz --fasta genome.fa.gz

GFF file

Example of command line with GFF, using of flag –gff :

1
./vep -i input.vcf --cache --gff data.gff.gz --fasta genome.fa.gz

This functionality uses VEP’s custom annotation feature, and the –gff flag is a shortcut to:

1
--custom data.gff.gz,,gff

NOTE: You should use the longer custom annotation form if you wish to customise the name of the GFF as it appears in the SOURCE field and VEP output header.

GTF file

Example of command line with GTF, using of flag –gtf :

1
./vep -i input.vcf --cache --gtf data.gtf.gz --fasta genome.fa.gz

This functionality uses VEP’s custom annotation feature, and the –gtf flag is a shortcut to:

1
--custom data.gtf.gz,,gtf

NOTE: You should use the longer custom annotation form if you wish to customise the name of the GTF as it appears in the SOURCE field and VEP output header.

定制注释内容（Custom annotations）

VEP 可以使用 --custom 标志将标准格式文件中的自定义注释集成到您的结果中。
对应的数据库可以是任意使用 tabix 构建过索引的文件（BED、GFF、GTF、VCF）或 bigWig (文件包含自己的索引)。

Format Type Description Notes
GFF/GTF Gene/transcript annotations Formats to describe genes and other genomic features — format specifications: GFF3 and GTF Requires a FASTA file in offline mode or if the desired species or assembly is not part of the Ensembl species list.
VCF Variant data A format used to describe genomic variants VEP uses the 3rd column as the identifier. INFO and FILTER fields from records may be added to the VEP output.
BED Basic/uninterpreted data A simple tab-delimited format containing 3-12 columns of data. The first 3 columns contain the coordinates of the feature. VEP uses the 4th column (if available) as the feature identifier.
bigWig Basic/uninterpreted data A format for storage of dense continuous data. VEP uses the value for the given position as the identifier. BigWig files contain their own indices, and do not need to be indexed by tabix. Requires Bio::DB::BigFile.

定制内容的使用参数说明

1
2
3
4
5
6
./vep [...] --custom Filename,Short_name,File_type,Annotation_type,Force_report_coordinates,VCF_fields
# Filename: 定制注释对应的库文件名称
# Short_name：注释后展示的名称
# File type： 定制库文件的文件类型
# Annotation type： "exact" or "overlap"  当使用“exact”时，仅报告其坐标与变体的坐标完全匹配的注释。这适用于位置特定信息，例如保守分数、等位基因频率或表型信息。使用“overlap”，任何与变体重叠甚至 1bp 的注释都会被报告。
# VCF_fields: 您可以指定自定义输入 VCF 的 INFO 字段中存在的任何信息类型（例如“AC”），或为 FILTER 字段指定“FILTER”，以将这些添加为自定义注释：

详细的参数说明和具体使用示例，参考 VEP官方文档

VEP的结果逻辑说明

Ensembl Variation - Calculated variant consequences

image
同时由于一个突变可能会存在多个不同类型的变异描述方式（Stop lost; Inframe insertion），在某些情况下需要形成一个唯一的功能描述，因此Ensemble提供了一套用于评估不同变异描述的优先级顺序，可以参考官方文档。

VEP使用过程中发现的问题/风险汇总

####1. 注释结果未见终止密码子，功能输出 stop_lost

1
2
3
变异结果： chr1_10342506_TCAGGTGGGCTTGACGTCTGTGACCAGTATTCAAGAGAGGATCATGTCTA/CCAGGTGTAGACATGATCCTCTCTTGAATACTGGTCACAGACGTCAAGCC
注释结果： NM_015074.3：c.1211_1260delinsCCAGGTGTAGACATGATCCTCTCTTGAATACTGGTCACAGACGTCAAGCC	
注释结果： NP_055889.2：p.I404_L420delinsTRCRHDPLLNTGHRRQA

在UCSC中核查发现读码框不一致。

1
2
3
4
5
6
7
8
9
VEP读码框为    |a T C|A G G|T G G|G C T|T G A|C G T|C T G|T G A|C ......
UCSC读码框为  a g T|C A G|G T G|G G C|T T G|A C G|T C T|G T G|A C ......

在VEP使用的密码框中，存在多个读码框对应的氨基酸序列为终止密码子（TGA），同时namechecker和Varsome中的突变后氨基酸均基于VEP的读码框翻译产生。

Ref / Alt 
aTCAGGTGGGCTTGACGTCTGTGACCAGTATTCAAGAGAGGATCATGTCTA/aCCAGGTGTAGACATGATCCTCTCTTGAATACTGGTCACAGACGTCAAGCC
I  R  W  A  *  R  L  *  P  V  F  K  R  G  S  C  L  /T  R  C  R  H  D  P  L  L  N  T  G  H  R  R  Q  A
04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20

经过核查是用于使用了 –use_given_ref 参数，当不使用该参数是测试发现问题可以消除。

2. start-lost /init loss 的写法不规范，

其实密码子缺失应该输出为 “p.0” 但是VEP输出结果格式为 “p.Leu1?”

3. VEP 注释结果中，存在部分变异结果没有提供外显子亚区

测试阶段发现如下情况时，变异注释结果未输出外显子，

  1. 变异位于内含子和外显子交界处，示例如下：
    1
    2
    3
                 Exon-1          | ins |        Intron-1
    Ref: TCAGGTGGGCTTGACGTCTGTGAC|     |CAGTATTCAAGAGAGGATCATGTCTA
    Alt: TCAGGTGGGCTTGACGTCTGTGAC|TTGAA|CAGTATTCAAGAGAGGATCATGTCTA

在使用中，可以进行了人工数据库的构建，但是需要注意一点，bed文件格式为0base的左闭右开，但是VEP的位置信息1base的

4. 跨越多个区域的变异，是否需要有氨基酸变化的描述

方法 cHGVS pHGVS
VEP c.111+1_117delinsCTGCGC p.A38_*39delinsLR
Varsome c.111+1_117delGCTGCGGinsCTGCGC .
  • 不同方法注释得到的氨基酸变化存在差异，现阶段无法确定正确答案。
Tag value
CHR chr10
Pos 43609928
Ref .
Alt GTCCACTGTGCGACGAGCTGTGCCGCACGGTGGATCTGTGGGTGGGGGTGGTGTGAGGCTTGGCACCGCCACCCACAG
Trans NM_020975.6
方法 cHGVS pHGVS
VEP c.1880-1_1880ins GTCCACTGTGCGACGAGCTGTGCCGCACGGTGGATCTGTGGGTGGGGGTGGTGTGAGGCTTGGCACCGCCACCCACAGA p.Q626_D627insGPLCDELCRTVDLWVGVV*
Varsome c.1880-1_1880insGTCCACTGTGCGACGAGCTGTGCCGCACGGTGGATCTGTGGGTGGGGGTGGTGTGAGGCTTGGCACCGCCACCCACAG p.D627_S1114delinsRVHCATSCAARWICGWGWCEAWHRHPQIHCATSCAAR
Namechecker c.1879_1880insGTCCACTGTGCGACGAGCTGTGCCGCACGGTGGATCTGTGGGTGGGGGTGGTGTGAGGCTTGGCACCGCCACCCACAG p.Asp627_Ser1114delinsGlyProLeuCysAspGluLeuCysArgThrValAspLeuTrpValGlyValVal
  • HGVS未明确描述导致氨基酸截断时描述规范
Tag value
CHR chr19
Pos 42795240
Ref CGCCCCCTGCTGTCCAGTT
Alt GGCAATGAACTGGACAGCA
Trans NM_015125.5
方法 cHGVS pHGVS
VEP c.2321_2339delinsGGCAATGAACTGGACAGCA p.Ala774_Phe780delinsGlyGlnTer
Varsome c.2321_2339delinsGGCAATGAACTGGACAGCA p.Ala774_Arg1608delinsGlyGln
Namechecker c.2321_2339delinsGGCAATGAACTGGACAGCA p.Ala774_Arg1608delinsGlyGln

变异导致多个氨基酸均为同一突变

-------------本文结束感谢您的阅读-------------