一般找公司做的标准分析，会根据客户指定的要求，把基本的质控、比对、call SNP以及annovar注释等步骤做完，这些也比较流程化，难点是客户根据自己的项目需求做下游分析。如下图所示，就是公司交付的WES标准分析的结果文件，包括SNP以及INDEL分析
今天，我想用这个XXX.annovar.hg19.multianno.xls.gz文件，利用maftools去做瀑布图。那么首先我们需要将这个文件转换为maftools识别的文件格式，这里有一篇帖子给了我很大帮助利用maftools分析annovar的注释结果 - 知乎 (zhihu.com)，后面我有一些代码就是在这个帖子的基础上修改的。
1.首先，annovarToMaf()这个函数是在maftools包里,我们打开Rstudio,安装maftools,然后加载bao,查看函数annovarToMaf()的帮助文档，里面有一些示例：
library(maftools)
?annovarToMaf
这个函数的介绍文档详细介绍了各个参数的用法以及默认设置，往下接着看我们发现还有示例：
上面的示例，就是把一个名为"variants.hg19_multianno.txt"的annovar文件转换为maf格式的文件。那么我们要想利用此函数将我们自己的annovar文件转变为maf文件的话，首先格式上要和"variants.hg19_multianno.txt"示例文件相同。最直接的方式，找到这个示例文件，看一下格式。找到示例文件也不难，我们先运行一下示例
var.annovar <- system.file("extdata", "variants.hg19_multianno.txt", package = "maftools")
var.annovar.maf <- annovarToMaf(annovar = var.annovar, Center = 'CSI-NUS', refBuild = 'hg19',
tsbCol = 'Tumor_Sample_Barcode', table = 'ensGene')
然后在Rstudio右上方Environment模块我们可以看到示例文件"variants.hg19_multianno.txt"的具体存放位置：
找到后，用记事本打开，格式如下：
从上图可以看到，示例input文件一共包含11列，tab分隔，参照此文件格式，我开始改造我的annovar文件。
先看一下诺和致源给的XXX.annovar.hg19.multianno.xls.gz，有79个字段，包含了各个数据库的注释信息，但是我们并不需要那么多信息，下图显示的是前15个字段:
我先将所有的XXX.annovar.hg19.multianno.xls.gz解压缩，然后通过下面这个shell脚本，将所有的XXX.annovar.hg19.multianno.xls文件整合成一个类似"variants.hg19_multianno.txt"文件的格式：
#!/usr/bin/bash
for i in *.hg19_multianno.xls
do
sample=`echo $i|awk -F '.' '{print $1}'`
#过滤掉内含子区以及基因间区的SNP，以及千人基因组计划数据的中国人群中，该变异位点上突变碱基的等位基因频率<0.01的SNP位点
awk -F '\t' '$1!="L" && $11!="intergenic" && $11!="intronic" && $27<0.01 {print "Chr"$2"\t"$3"\t"$3"\t"$5"\t"$6"\t"$9"\t"$10"\t"$14"\t"$15"\t"$11}' $i|sed '1d'|sed "s/$/\t${sample}/">>all_sample4.txt
done
sed -i '1s/^/Chr\tStart\tEnd\tRef\tAlt\tGene.refGene\tGeneDetail.refGene\tExonicFunc.refGene\tAAChange.refGene\tFunc.refGene\tTumor_Sample_Barcode\n/' all_sample4.txt
完成啦，接下来在Rstudio中用all_sample4.txt就好啦
var.annovar.maf = annovarToMaf(annovar = "all_sample4.txt",
Center = 'NA',
refBuild = 'hg19',
tsbCol = 'Tumor_Sample_Barcode',
table = 'refGene',
sep = "\t")
write.table(var.annovar.maf,file="var_annovar.maf",quote= F,sep="\t",row.names=F)
var_maf = read.maf(maf ="var_annovar.maf")
plotmafSummary(maf = var_maf, rmOutlier = TRUE, addStat = 'median')
oncoplot(maf = var_maf, top = 30, fontSize = 12 ,showTumorSampleBarcodes = F )
somaticInteractions(maf = var_maf, top = 25, pvalue = c(0.05, 0.1))
