mskcc的vcf2maf极简解决方案

Posted on by

为了写这个教程，我特意在唐医生的共享云服务器上面测试了，从头到尾运行过，验证过，你一定可以follow成功的哈！

首先是安装miniconda

https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/?C=M&O=A

wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
source ~/.bashrc
vi .condarc

查看 https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/ 内容，修改你的.condarc 文件，主要是方便中国大陆的朋友们下载软件。

然后使用conda进行软件安装（vcf2maf和VEP）

conda create -n vep -y 
conda activate vep

conda install -y -c bioconda vcf2maf
conda install -y -c bioconda ensembl-vep
conda remove samtools
conda install -y -c bioconda samtools

vcf2maf.pl --hlep
vep --help 
samtools -v
perl -e '{print join"\n",@INC}'

接着配置VEP的数据库资源文件夹

export VEP_PATH=$HOME/vep
export VEP_DATA=$HOME/.vep
mkdir -p $VEP_PATH $VEP_DATA; cd $VEP_PATH
# Ensembl Variant Effect Predictor (VEP) 
vep_install -a cf -s homo_sapiens -y GRCh38 -c $VEP_PATH --CONVERT
# downloading ftp://ftp.ensembl.org/pub/release-88/variation/VEP/homo_sapiens_vep_88_GRCh38.tar.gz
# 同时也会下载对应的参考基因组fa文件的压缩包哦！
ls -lh $HOME/vep 
# vep_install -a cf -s homo_sapiens -y GRCh37 -c $VEP_PATH --CONVERT
# vep_install -a cf -s Mus_musculus -y GRCm38 -c $VEP_PATH --CONVERT

一般来说，网络速度是一个限制！也可以手动选择下载最新版数据库文件，下面的代码无需运行，仅仅是举个例子给大家哈：

mkdir -p $HOME/.vep
cd $HOME/.vep
nohup wget ftp://ftp.ensembl.org/pub/release-101/variation/indexed_vep_cache/homo_sapiens_vep_101_GRCh38.tar.gz & 
tar xzf homo_sapiens_vep_101_GRCh38.tar.gz
nohup wget ftp://ftp.ensembl.org/pub/release-101/variation/indexed_vep_cache/homo_sapiens_vep_101_GRCh37.tar.gz & 
nohup wget ftp://ftp.ensembl.org/pub/release-101/variation/indexed_vep_cache/mus_musculus_refseq_vep_101_GRCm38.tar.gz & 
nohup wget ftp://ftp.ensembl.org/pub/release-101/variation/indexed_vep_cache/mus_musculus_vep_101_GRCm38.tar.gz &

注意一下，两个数据库文件夹目录不一样哦，调用VEP的时候注意区分它们。有意思的是，我下载的这几个最新版数据库文件居然会报错？？后来我还是使用的默认的 homo_sapiens_vep_88_GRCh38.tar.gz 版本文件。

单独运行VEP

VEP的全称是variant_effect_predictor，就是把vcf文件里面的每个变异位点的坐标，根据VEP软件自带的数据集，进行overlap后，就能给出每个变异位点的一些注释信息：

conda activate vep 
vep --help
# By default the VEP uses $HOME/.vep/ 
vep -i input_varscan.snp.Somatic.hc.vcf -o test.vcf \
--cache --dir_cache $HOME/vep --force_overwrite --assembly GRCh38 --vcf

注释前后的文件都是vcf格式，这样的注释，通常是针对germline的突变信息；

最后运行 mskcc的vcf2maf

因为mskcc的vcf2maf运行的时候也是会调用VEP，所以需要先测试VEP软件是否成功，然后使用下面的脚本：

conda activate vep 
vcf2maf.pl --help
# 注意下面的 --normal-id NORMAL --tumor-id TUMOR ，不同的somatic 软件不一样哦
# 我举例的这个来自于 varscan 哈！
ref=$HOME/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
vcf2maf.pl --input-vcf T1520021_varscan.snp.Somatic.hc.vcf \
--output-maf test.maf --normal-id NORMAL --tumor-id TUMOR \
--ref-fasta $ref \
--vep-data $HOME/vep \
--vep-path ~/miniconda3/envs/vep/bin/ \
--ncbi-build GRCh38 
# 然后我上面的vep因为是conda安装，所以路径是自定义。

这样的mskcc的vcf2maf流程通常是针对somatic的突变信息。

有一个小麻烦，就是conda 安装软件呢，容易造成版本不匹配，mskcc的vcf2maf的VEP如果不匹配，会报错，需要手动编辑文件 vcf2maf.pl ：

$which vcf2maf.pl
~/miniconda3/envs/vep/bin/vcf2maf.pl

参考 https://github.com/mskcc/vcf2maf/blob/master/vcf2maf.pl 定位到下面的代码，注释掉里面的专门针对newer VEP设置的参数即可。

if( $species eq "homo_sapiens" ) {
 # Slight change in options if in offline mode, or if using the newer VEP
 $vep_cmd .= " --polyphen b" . ( $vep_script =~ m/vep$/ ? " --af" : " --gmaf" );
 $vep_cmd .= ( $vep_script =~ m/vep$/ ? " --af_1kg --af_esp --af_gnomad" : " --maf_1kg --maf_esp" ) unless( $online );
}

这样你的vcf文件就全部转为maf格式的啦，后续只需要走maftools即可进行花式统计可视化啦。

如果是多个vcf文件批量转maf

写一个脚本，我的脚本如下：

针对varscan软件的somatic的snp：

ls *_varscan.snp.Somatic.hc.vcf |while read id;do 
 sample=$(basename "$id" _varscan.snp.Somatic.hc.vcf)
 ref=$HOME/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
 grep -v '_' $id > ${sample}_filter_snp.vcf
 vcf2maf.pl --input-vcf ${sample}_filter_snp.vcf \
 --output-maf ${sample}_snp.maf --normal-id NORMAL --tumor-id TUMOR \
 --ref-fasta $ref \
 --vep-data $HOME/vep \
 --vep-path ~/miniconda3/envs/vep/bin/ \
 --ncbi-build GRCh38 
done

如果是indel文件

ls *_varscan.indel.Somatic.hc.vcf |while read id;do 
 sample=$(basename "$id" _varscan.indel.Somatic.hc.vcf)
 ref=$HOME/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
 grep -v '_' $id > ${sample}_filter_indel.vcf
 vcf2maf.pl --input-vcf ${sample}_filter_indel.vcf \
 --output-maf ${sample}_indel.maf --normal-id NORMAL --tumor-id TUMOR \
 --ref-fasta $ref \
 --vep-data $HOME/vep \
 --vep-path ~/miniconda3/envs/vep/bin/ \
 --ncbi-build GRCh38 
done

多个maf文件可以合并，但是呢，对于这个varscan软件的somatic的vcf转为的maf文件，合并的时候需要注意了，小心 —normal-id 和 —tumor-id 信息的丢失。

如果是mutect2的vcf文件，代码稍微复杂一点，需要config文件，主要是3列信息：

  • 第一列是肿瘤命名
  • 第二列是肿瘤病人的normal组织的bam文件地址
  • 第三列是肿瘤病人的肿瘤组织的bam文件地址。
conda activate vep 
ref=$HOME/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
cat config |while read id
do
 arr=($id)
 normal_bam=${arr[1]}
 tumor_bam=${arr[2]}
 sample=${arr[0]}
 if [ ! -f ${sample}.maf ]; then
 vcf2maf.pl --input-vcf results/${sample}_pass.vcf --output-maf ${sample}.maf \
 --ref-fasta $ref \
 --vep-data $HOME/vep \
 --vep-path ~/miniconda3/envs/vep/bin/ \
 --ncbi-build GRCh38 \
 --tumor-id $(basename "$tumor_bam" _recal.bam) --normal-id $(basename "$normal_bam" _recal.bam) 
 fi
done
#grep '^Hugo_Symbol' *.maf |head -1|cut -d":" -f 2 > ../human_brca_all_mutect2.maf
#cat *.maf|grep -v '^Hugo_Symbol' |grep -v "^#" >> ../human_brca_all_mutect2.maf
​

Comments are closed.