mskcc的vcf2maf极简解决方案代码分享
为了写这个教程,我特意在唐医生的共享云服务器上面测试了,从头到尾运行过,验证过,你一定可以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 信息的丢失。
- 基础篇章:React Native 之 View 和 Text 的讲解
- CentOs7.3 修改主机名
- 基础篇章:React Native之Flexbox的讲解(Height and Width)
- PDF.js专题
- CentOs7.3 编译安装 Nginx 1.9.9
- 基础篇章:关于 React Native 之 RefreshControl 组件的讲解
- CentOs7.3 安装 JDK1.8
- 基础篇章:关于 React Native 之 ListView 组件的讲解
- maven环境快速搭建
- CentOs7.3 搭建 RabbitMQ 3.6 Cluster 集群服务
- CentOs7.3 搭建 Redis-4.0.1 Cluster 集群服务
- CentOs7.3 搭建 Redis-4.0.1 单机服务
- Shodan新手入坑指南
- 我用过的——Spring定时任务的几种用法
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- 关于 z-index,你可能一直存在误区
- Java浮点数机制及所存在的问题
- SAP UI5和React的页面渲染性能比较
- SAP CDS view自学教程之一:如何测试基于SAP CDS view自动生成的OData服务
- SAP CDS view自学教程之二:当SAP CDS view被激活时,背后发生了什么
- SAP Fiori Elements原理介绍之类型为Value Help的Smart Field工作原理
- 万字详述 MySQL ProxySQL
- 在SAP WebClient UI里使用AJAX进行异步数据读取
- Angular Component TypeScript代码和最后转换生成的JavaScript代码比较
- 如何使用Angular FormBuilder
- Angular HTML template的解析位置
- Angular FormBuilder的工作原理
- Angular HTTPClient的使用方法
- nodejs错误:PayloadTooLargeError: request entity too large
- 富文本编辑器 tinymce 的使用