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.shbash Miniconda3-latest-Linux-x86_64.shsource ~/.bashrcvi .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 vcf2mafconda install -y -c bioconda ensembl-vepconda remove samtoolsconda install -y -c bioconda samtools
vcf2maf.pl --hlepvep --help samtools -vperl -e '{print join"\n",@INC}'

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

export VEP_PATH=$HOME/vepexport VEP_DATA=$HOME/.vepmkdir -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/.vepcd $HOME/.vepnohup 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.gznohup 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.fastavcf2maf.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 信息的丢失。

(0)

相关推荐

  • 保姆级参考基因组及其注释下载教程(图文详解)

    目录一.什么是参考基因组和基因组注释?二.参考基因组版本命名1.常用人参考基因组对应表2.常用小鼠参考基因组对应表三.下载1.NCBI2.Ensemble3.GENCODE4.UCSC5.iGenom ...

  • 这25条极简Python代码,你还不知道

    GitChat 3月3日 以下文章来源于法纳斯特 ,作者小F 作者 | 小F 来源 | 法纳斯特 自从我用Python编写第一行代码以来,就被它的简单性.出色的可读性和特别流行的一行代码所吸引. 下面 ...

  • 30段极简Python代码:这些小技巧你都Get了么

    数学算法俱乐部 日期 : 2021年06月14日 正文共 :3367字 来源 : 机器之心 学 Python 怎样才最快?当然是实战各种小项目,只有自己去想与写,才记得住规则.本文是 30 个极简任务 ...

  • 30个极简Python代码,拿走就能用!

    最码农 2月28日 ‍ 作者丨Fatos Morina 来源丨Python 技术 编辑丨极市平台 学 Python 怎样才最快,当然是实战各种小项目,只有自己去想与写,才记得住规则. 本文是 30 个 ...

  • 30个极简Python代码,拿走即用

    小白带你学编程 今天 学 Python 怎样才最快,当然是实战各种小项目,只有自己去想与写,才记得住规则.本文是 30 个极简任务,初学者可以尝试着自己实现:本文同样也是 30 段代码,Python ...

  • ​这25条极简Python代码,你还不知道

    人工智能头条 1周前 以下文章来源于法纳斯特 ,作者小F 法纳斯特分享学习Python爬虫.数据分析.数据挖掘的点滴. 作者 | 小F 来源 | 法纳斯特 自从我用Python编写第一行代码以来,就被 ...

  • 30 个极简Python代码,拿走即用(真干货)

    学Python怎样才最快,当然是实战各种小项目,只有自己去想与写,才记得住规则.本文是30个极简任务,初学者可以尝试着自己实现:本文同样也是 30段代码,Python 开发者也可以看看是不是有没想到的 ...

  • 30个极简Python代码,拿走即用(真干货)

    重复元素判定 def all_unique(lst): return len(lst)==len(set(lst)) x = [1,1,2,2,3,2,3,4,5,6] y=[1,2,3,4,5] a ...

  • 家居好物分享,极简主义灯饰,好看百搭的吊...

    家居好物分享,极简主义灯饰,好看百搭的吊灯. 有一家是某东店,其余都是某宝店. 1.太阳湖 追求简单线条感,百看不厌. 2.纽思维尔家居 很有创意,也适用于办公室装饰,分分钟提高逼格. 3.星月灯饰设 ...

  • 30个Python常用极简代码,拿走就能用

    作者丨Fatos Morina 来源丨Python 技术 导读 本文是 30 个极简任务,初学者可以尝试着自己实现:本文同样也是 30 段代码,Python 开发者也可以看看是不是有没想到的用法. 学 ...