肿瘤外显子数据处理系列教程（八）不同注释软件的比较（中）：注释后转成maf文件
上一节我们开始比较不同注释软件的安装和使用方法，但是后来检查的时候发现代码有一些的问题，GATK的Funcatator工具注释方法应该还需要一个interval文件，interval格式和bed格式几乎一致，只不过interval文件的编号是从“1”开始，而bed是从“0”开始，可以把bed文件转成interval文件，用gatk的BedToIntervalList工具：
GATK=/home/hcguo/biosoft/gatk/gatk-4.1.4.0/gatk
bed=/home/llwu/reference/index/human/CCDS/GRCh38.bed
dict=/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.dict
$GATK BedToIntervalList -I ${bed} -O GRCh38.interval_list -SD ${dict}
还有一个问题就是用上面这个数据库文件所做的注释会有很多条目是__UNKNOWN__，可以考虑换一个更大的数据库文件，压缩包有14G，解压后是18G：
nohup wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/funcotator/funcotator_dataSources.v1.6.20190124s.tar.gz &
tar -zxvf funcotator_dataSources.v1.6.20190124s.tar.gz
然后再走funcatator的注释脚本：
GATK=/home/hcguo/biosoft/gatk/gatk-4.1.4.0/gatk
ref=/public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
interval=/home/hcguo/SRP070662/8.mutect2/SRP070662_filter_vcf/GRCh38.interval_list
source=/home/hcguo/database/gatk_funcotator/funcotator_dataSources.v1.6.20190124s
cat config | while read id
do
sample=$id
echo "start Funcotator for ${id} " `date`
$GATK Funcotator -R $ref \
-V ${sample}_filter.vcf \
-O gatk/${sample}_funcotator.tmp.maf \
--data-sources-path ${source} \
--intervals ${interval} \
--output-file-format MAF \
--ref-version hg38
echo "end Funcotator for ${id} " `date`
done
ok，纠正过了这个问题，接下来我们继续后面的分析。
vcf文件转成maf
ANNOVAR注释结果转maf
ANNOVAR注释后拿到的结果，还需要做一个处理，就是在每一个*.hg38_multianno.txt文件中增加两列，为Tumor_Sample_Barcode和Matched_Norm_Sample_Barcode，而最后的一列Otherinfo我们暂时丢弃掉，添加上这两列信息之后再把所有的vcf文件合并成一个annovar_merge.vcf。处理的方法如下：
cat config | while read id
do
grep -v '^Chr' annovar/${id}.hg38_multianno.txt | cut -f 1-20 | awk -v T=${id} -v N=${id:0:5}_germline '{print $0"\t"T"\t"N}' >annovar/${id}.annovar.vcf
done
head -1 annovar/case1_biorep_A_techrep.hg38_multianno.txt| sed 's/Otherinfo/Tumor_Sample_Barcode\tMatched_Norm_Sample_Barcode/' >annovar/header
cat annovar/header annovar/*annovar.vcf >annovar/annovar_merge.vcf
最后可以用R包maftools中一个函数annovarToMaf将ANNOVAR注释的结果转为maf
rm(list = ls())
require(maftools)
options(stringsAsFactors = F)
## annovar
annovar.laml <- annovarToMaf(annovar = "annovar/annovar_merge.vcf",
refBuild = 'hg38',
tsbCol = 'Tumor_Sample_Barcode',
table = 'refGene',
MAFobj = T)
## 统计一下case1_biorep_A_techrep样本的突变类型
table(annovar.laml@data[annovar.laml@data$Tumor_Sample_Barcode=='case1_biorep_A_techrep',]$Func.knownGene)
## exonic intergenic
## 172 1
table(annovar.laml@data[annovar.laml@data$Tumor_Sample_Barcode=='case1_biorep_A_techrep',]$Variant_Type)
## DEL INS SNP
## 5 0 168
table(annovar.laml@data[annovar.laml@data$Tumor_Sample_Barcode=='case1_biorep_A_techrep',]$Variant_Classification)
## Frame_Shift_Del
## 4
## Frame_Shift_Ins
## 0
## In_Frame_Del
## 1
## In_Frame_Ins
## 1
## Missense_Mutation
## 153
## Nonsense_Mutation
## 14
## Nonstop_Mutation
## 0
最后拿到的annovar.laml，就是一个maftools对象了，可以进行后续的画图操作，这部分我们留到下一节再讲。
GATK注释结果转maf
gatk的Funcotator虽然注释上了很多数据库信息，但是很多结果为__UNKNOWN__。注释后直接输出的maf文件也是有缺陷的，Tumor_Sample_Barcode和Matched_Norm_Sample_Barcode这两列在maf文件中也是__UNKNOWN__，可能是我的参数设置有问题，但是看了帮助文档：https://software.broadinstitute.org/gatk/documentation/tooldocs/4.1.4.0/org_broadinstitute_hellbender_tools_funcotator_Funcotator.php，并没有发现相关的参数，知道的朋友可以留言。在此，我只能暂时手动添加一下了，方法如下：
## 添加Barcode
cat config | while read id
do
grep -v '^#' gatk/${id}_funcotator.tmp.maf| grep -v '^Hugo_Symbol'| awk -v T=${id} -v N=${id:0:5}_germline 'BEGIN{FS="\t";OFS="\t"}{$16=T;$17=N;print $0}' >gatk/${id}_funcotator.maf
done
## 取出一个列名
grep 'Hugo_Symbol' gatk/case1_biorep_A_techrep_funcotator.tmp.maf >gatk/header
## 删除掉旧文件
rm gatk/*tmp.maf
## 合并所有的样本的maf
cat gatk/header gatk/*maf >gatk/gatk4.1.4.0_merge.maf
这样每个样本就得到一个新的maf文件，当然也拿到了最后的合并的gatk4.1.4.0_merge.maf，新的maf文件里面的每一条记录都增加了Barcode的信息，如case1_biorep_A_techrep_funcotator.maf：
head -1 gatk/case1_biorep_A_techrep_funcotator.maf
CHD5 26038 __UNKNOWN__ hg38 chr1 6146376 6146376 + Missense_Mutation SNP G G T __UNKNOWN__ __UNKNOWN__ case1_biorep_A_techrep case1_germline __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ __UNKNOWN____UNKNOWN__ NA NA __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ g.chr1:6146376G>T ENST00000262450.7 - 11 1738 c.1638C>A c.(1636-1638)aaC>aaA p.N546K NM_015557.2 NP_056372 O00258 WRB_HUMAN chromodomain helicase DNA binding protein 5 __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ tail-anchored membrane protein insertion into ER membrane (GO:0071816) endoplasmicreticulum (GO:0005783)|integral component of membrane (GO:0016021)|nucleus (GO:0005634) biliary_tract(2)|breast(48)|central_nervous_system(44)|kidney(92)|large_intestine(37)|lung(1)|pancreas(22)|upper_aerodigestive_tract(120) 366 __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ __UNKNOWN__CATCCATGTCGTTCTTTCTTT 0.5860349127182045 __UNKNOWN__ __UNKNOWN__ __UNKNOWN__ 0.058 3 97 0 75 AL117491 HGNC:16816 Approved gene with protein product protein-coding gene 1p36.31 2016-10-05 AF425231 ENSG00000116254 11889561, 12592387 NM_015557 1305|88 NuRcomplex|PHD finger proteins CCDS57 OTTHUMG00000000952 610771 NM_015557 Q8TDI0 ENSG00000116254 uc001amb.3 A8KAP8|A8MQ44|D3DSH9|O60740 false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false false 26038 false 184 1 22.56 -1.888e+00 1.000e-05 0.00 -4.514e+01 5.41
在第9列和第10列分别记录了Variant_Classification和Variant_Type，我们也可以适当统计一下：
$less gatk/case1_biorep_A_techrep_funcotator.maf | cut -f 9 | sort | uniq -c
1 3'UTR
4 Frame_Shift_Del
1 In_Frame_Del
12 Intron
142 Missense_Mutation
14 Nonsense_Mutation
60 Silent
1 Splice_Site
1 START_CODON_SNP
$ less gatk/case1_biorep_A_techrep_funcotator.maf| cut -f 10 | sort | uniq -c
5 DEL
1 DNP
230 SNP
这里有一个DNP，我们就顺便解释一下关于几种变异类型：
INS - The variant allele is some kind of insertion.
DEL - The variant allele is some kind of deletion.
SNP - The variant allele is a single nucleotide polymorphism.
DNP - The variant allele is a di-nucleotide polymorphism.
TNP - The variant allele is a tri-nucleotide polymorphism.
ONP - The variant allele is an oligo-nucleotide polymorphism (Synonymous with MNP).
MNP - The variant allele is a multi-nucleotide polymorphism (Synonymous with ONP).
NA - The variant allele type cannot be determined.
VEP注释结果转maf
vcf2maf
VEP的注释结果拿到了vcf文件，也可以转成maf文件，不过要用到一个软件，叫vcf2maf，安装和使用的方法见曾老师教程：把vcf文件转换为maf格式或者官网vcf2maf：https://github.com/mskcc/vcf2maf，也很简单：
export VCF2MAF_URL=`curl -sL https://api.github.com/repos/mskcc/vcf2maf/releases | grep -m1 tarball_url | cut -d\" -f4`
curl -L -o mskcc-vcf2maf.tar.gz $VCF2MAF_URL
tar -zxf mskcc-vcf2maf.tar.gz
cd mskcc-vcf2maf-*
转换格式主要用到的是一个perl脚本vcf2maf.pl，同样的我们还是写一个shell脚本调用perl程序，shell脚本vcf2maf.sh如下，输入文件为vep注释后的vcf，输出文件为maf文件：
cat config | while read id
do
echo $id
sample=${id}_vep.vcf
normal=${id%%_*}_germline
perl ~/biosoft/vcf2maf/vcf2maf.pl \
--input-vcf VEP_annotation/${sample} \
--output-maf VEP_annotation/${id}_vep.maf \
--ref-fasta /public/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta \
--tumor-id ${id} \
--normal-id ${normal} \
--ncbi-build GRCh38
done
在第9列和第10列分别记录了Variant_Classification和Variant_Type，同样查看一下样本case1_biorep_A_techrep转换后的结果：
less vep/case1_biorep_A_techrep_vep.maf |cut -f 9|sort |uniq -c
4 Frame_Shift_Del
1 In_Frame_Del
153 Missense_Mutation
15 Nonsense_Mutation
59 Silent
3 Splice_Region
1 Translation_Start_Site
1 Variant_Classification
less vep/case1_biorep_A_techrep_vep.maf |cut -f 10 |sort |uniq -c
5 DEL
1 DNP
230 SNP
注释后的maf文件的记录和原先的vcf文件记录的条目是一致的，都是236条记录。
最后也是把所有的maf文件合并到一起，拿到VEP_merge.maf
cat VEP_annotation/*maf | grep -v '^#'| grep -v '^Hugo_Symbol' >VEP_annotation/tmp
grep 'Hugo_Symbol' case1_biorep_A_techrep_vep.maf >VEP_annotation/header
cat VEP_annotation/header VEP_annotation/tmp >VEP_annotation/VEP_merge.maf
SNPEFF注释结果转maf
SNPEFF的注释结果和vep的注释结果非常相似，都是在原来的vcf文件的INFO列添加了许多注释信息。但是，vcf2maf似乎不支持snpeff的注释结果，无法转成maf，搜索了很久也都没有找到有效的方法，按理来说稍稍改一下vcf2maf.pl的脚本就可以了，但是我没学过perl，面对上千行的perl脚本，无从下手，不知道各位是否已经造了轮子，如果有的话，希望可以邮件分享一下Nickier0510@163.com，谢谢~