全conda安装vep、vcf2maf

全conda安装vep、vcf2maf

vcf2maf的使用主要依赖vep，由于vep是由perl编译的，所以在配置和安装vep的时候会出现各种各样的问题。

我尝试了各种各样的方法后，摸索出如何只用conda安装使用vep和vcf2maf。

ps：看到有说可以使用vep docker的，但是我们的服务器上没有安装docker，所以只能一步一步用conda安装，如果想更省事点可以参考这篇帖子

vep docker下使用vcf2maf - 简书 (jianshu.com)

写在前面

其实如果在使用conda安装过程中遇到明明下载了依赖的perl模块，但在启动vep的时候仍然报错，那大概率是因为perl的版本和模块版本不兼容，默认下载的perl是5.34.0，而安装vep默认下载的perl模块是5.22.0，因此在安装perl和vep的时候必须指定版本。

创建环境

conda create -n vep
conda activate vep

安装perl

conda install perl=5.26.2

安装vep

基于perl安装的版本，可以使用conda search ensembl-vep，查看build列来查看是用哪个版本的perl编译的，下载对应的vep就好

 conda install -c "bioconda/label/cf201901" ensembl-vep=93.2
 vep

若vep安装成功，命令行输入vep会输出如下结果

这里可以看到我下载的vep版本是93.2，后面配置对应的数据库版本也要是93.2

Possible precedence issue with control flow operator at ~/.conda/envs/vep/lib/site_perl/5.26.2/Bio/DB/IndexedBase.pm line 845.
#----------------------------------#
# ENSEMBL VARIANT EFFECT PREDICTOR #
#----------------------------------#

Versions:
  ensembl              : 93.b4d45ee
  ensembl-funcgen      : 93.0c98373
  ensembl-io           : 93.cc29a66
  ensembl-variation    : 93.9b583cc
  ensembl-vep          : 93.2

Help: dev@ensembl.org , helpdesk@ensembl.org
Twitter: @ensembl

http://www.ensembl.org/info/docs/tools/vep/script/index.html

Usage:
./vep [--cache|--offline|--database] [arguments]

Basic options
=============

--help                 Display this message and quit

-i | --input_file      Input file
-o | --output_file     Output file
--force_overwrite      Force overwriting of output file
--species [species]    Species to use [default: "human"]

--everything           Shortcut switch to turn on commonly used options. See web
                       documentation for details [default: off]
--fork [num_forks]     Use forking to improve script runtime

For full option documentation see:
http://www.ensembl.org/info/docs/tools/vep/script/vep_options.html

安装vcf2maf

conda install -c "bioconda/label/cf201901" vcf2maf=1.6.16=pl526h18118b2_3
vcf2maf --help

会自动安装版本合适的samtools，若vcf2maf安装成功，命令行输入vcf2maf --help将会输出如下内容

Usage:
     perl vcf2maf.pl --help
     perl vcf2maf.pl --input-vcf WD4086.vcf --output-maf WD4086.maf --tumor-id WD4086 --normal-id NB4086

Options:
     --input-vcf      Path to input file in VCF format
     --output-maf     Path to output MAF file
     --tmp-dir        Folder to retain intermediate VCFs after runtime [Default: Folder containing input VCF]
     --tumor-id       Tumor_Sample_Barcode to report in the MAF [TUMOR]
     --normal-id      Matched_Norm_Sample_Barcode to report in the MAF [NORMAL]
     --vcf-tumor-id   Tumor sample ID used in VCF's genotype columns [--tumor-id]
     --vcf-normal-id  Matched normal ID used in VCF's genotype columns [--normal-id]
     --custom-enst    List of custom ENST IDs that override canonical selection
     --vep-path       Folder containing the vep script [~/vep]
     --vep-data       VEP's base cache/plugin directory [~/.vep]
     --vep-forks      Number of forked processes to use when running VEP [4]
     --buffer-size    Number of variants VEP loads at a time; Reduce this for low memory systems [5000]
     --any-allele     When reporting co-located variants, allow mismatched variant alleles too
     --ref-fasta      Reference FASTA file [~/.vep/homo_sapiens/91_GRCh37/Homo_sapiens.GRCh37.75.dna.primary_assembly.fa.gz]
     --filter-vcf     A VCF for FILTER tag common_variant. Set to 0 to disable [~/.vep/ExAC_nonTCGA.r0.3.1.sites.vep.vcf.gz]
     --max-filter-ac  Use tag common_variant if the filter-vcf reports a subpopulation AC higher than this [10]
     --species        Ensembl-friendly name of species (e.g. mus_musculus for mouse) [homo_sapiens]
     --ncbi-build     NCBI reference assembly of variants MAF (e.g. GRCm38 for mouse) [GRCh37]
     --cache-version  Version of offline cache to use with VEP (e.g. 75, 84, 91) [Default: Installed version]
     --maf-center     Variant calling center to report in MAF [.]
     --retain-info    Comma-delimited names of INFO fields to retain as extra columns in MAF []
     --min-hom-vaf    If GT undefined in VCF, minimum allele fraction to call a variant homozygous [0.7]
     --remap-chain    Chain file to remap variants to a different assembly before running VEP
     --help           Print a brief help message and quit
     --man            Print the detailed manual

配置vep

下载cache数据集

Index of /pub (ensembl.org) 在官网上进行对应vep版本的数据库下载，我安装的vep版本是93.2，所以我需要找到release-93，我选用的组装的版本是hg19，也就是GRCh37。

wget就可以了，也可以本地下载再上传至服务器。

6.8个G🙂，慢慢下吧。

下完后解压就行了

wget http://ftp.ensembl.org/pub/release-93/variation/indexed_vep_cache/homo_sapiens_vep_93_GRCh37.tar.gz
tar -zxvf homo_sapiens_vep_93_GRCh37.tar.gz

下载参考基因组

接下来下载参考基因组，Index of /pub (ensembl.org) 当然用自己的也可以。

wget -c ftp://ftp.ensembl.org/pub/grch37/current/fasta/homo_sapiens/dna/Homo_sapiens.GRCh37.dna.primary_assembly.fa.gz

如果基因组的fasta文件是单行序列，那么最好处理成多行，否则会报“单行长度超出上限”的错误

#指定每行序列的输出长度(为0的话，代表为一整行，默认的输出 长度是60个碱基)
seqkit seq test.fa -w 10 > test_10.fa (指定序列的长度为10)

测试

vcf2maf.pl --input-vcf input.vcf --output-maf test.maf --vep-path ~/.conda/envs/vep/bin/ --vep-data ~/vep_database/ --ref-fasta ~/database/hg19/ucsc_hg19_multirows.fa --filter-vcf 0 --species homo_sapiens

参数说明：

–input-vcf：输入的vcf文件

–output-maf：输出的maf文件名

–vep-path：可执行的vep程序，如果像我这样用conda安装，应该是在[conda_path]/envs/[env_name]/bin

–vep-data：下载的cache数据路径，只要到dataset目录即可

–ref-fasta：参考基因组文件，.fa或.fa.gz

–tumor-id ：输出的maf文件里面的tumor name

–normal-id ：输出的maf文件里面的normal name

–vcf-tumor-id ：输入的vcf文件里面的tumor name

–vcf-normal-id：输入的vcf文件里面的normal name

–filter-vcf ：是否进行vcf过滤，0就是不过滤

–species： 物种名称，这里要写到–vep-data 路径的下一级目录的名字。例如–species是 homo_sapiens，–vep-data写的是~/vep_dataset，那就意味着存在~/vep_dataset/homo_sapiens/，程序才能跑起来

报错解决

参考基因组一行碱基超过处理上限

STATUS: Running VEP and writing to: ./test.vep.vcf
Possible precedence issue with control flow operator at ~/.conda/envs/vep/lib/site_perl/5.26.2/Bio/DB/IndexedBase.pm line 845.

------------- EXCEPTION: Bio::Root::Exception -------------
MSG: Each line of the file must be less than 65,536 characters. Line 50 is 59373567 chars.
STACK: Error::throw

• 解决：

seqkit seq database/hg19/ucsc_hg19.fa > database/hg19/ucsc_hg19_multirows.fa    #默认每行60个碱基