肿瘤外显子数据处理系列教程（八）不同注释软件的比较(下)：可视化比较maf文件

在前面两节中，我们比较了不同软件的注释方法以及如何将注释结果转化为maf格式，由于snpeff注释结果暂时无法转换成maf，所以后面我们就比较ANNOVAR、VEP、GATK的注释结果得到的合并后的maf文件。

首先将拿到的结果载入到R语言中，然后用maftool可视化，代码如下：

读入R中先了解一下数据

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)
## gatk
# gatk.laml = read.maf(maf = 'gatk/gatk4.1.4.0_merge.maf')

library(data.table)
tmp=fread('gatk/gatk4.1.4.0_merge.maf')
gatk.laml = read.maf(maf = tmp)

## vep
vep.laml = read.maf(maf = 'vep/vep_merge.maf')
## for vep.laml
library(stringr)
vep.laml@data$Protein_Change = paste0("p.",
                                      str_sub(vep.laml@data$Amino_acids,1,1),
                                      vep.laml@data$Protein_position,
                                      str_sub(vep.laml@data$Amino_acids,3,3))

## save Rdata
save(annovar.laml, gatk.laml, vep.laml, file = 'laml.Rdata')

## Summary
laml=annovar.laml
unique(laml@data$Tumor_Sample_Barcode)
getSampleSummary(laml)
getGeneSummary(laml)
getFields(laml)

laml=gatk.laml
unique(laml@data$Tumor_Sample_Barcode)
getSampleSummary(laml)
getGeneSummary(laml)
getFields(laml)

laml=vep.laml
unique(laml@data$Tumor_Sample_Barcode)
getSampleSummary(laml)
getGeneSummary(laml)
getFields(laml)

上述步骤是把拿到的maf文件读入R中，初步处理后保存成Rdata，并且对数据进行初步了解。接下来是绘图可视化mafSummary

rm(list = ls())
require(maftools)
options(stringsAsFactors = F)
load(file = 'laml.Rdata')
laml=c(annovar.laml,gatk.laml,vep.laml)

## mafsummary
anno=c('annovar','gatk','vep')
for (i in 1:3) {
  #i = 1
  png(paste0('plotmafSummary_', anno[i], '.png'),res = 150,width = 1080,height = 1080)
  plotmafSummary(maf = laml[[i]],
                 rmOutlier = TRUE,
                 showBarcodes = T,
                 textSize = 0.4,
                 addStat = 'median',
                 dashboard = TRUE,
                 titvRaw = FALSE)
  dev.off()
}

可以看到，3个软件的注释结果总体上非常接近的，但是还是有一些小细节不同，比如我们看一下，top30基因的oncoplot

## oncoplot_top30
for (i in 1:3) {
  #i = 1
  png(paste0('oncoplot_top30_', anno[i], '.png'),res = 150,width = 1080,height = 1080)
  oncoplot(maf = laml[[i]],
           top = 30,
           fontSize = 0.5,
           sampleOrder = laml[[i]]@clinical.data$Tumor_Sample_Barcode,
           showTumorSampleBarcodes = T)
  dev.off()
}

从上面的图可以看到，3个软件注释后的可视化结果非常接近，当然也不是完全一模一样，比如说SP140基因，在ANNOVAR和VEP的注释结果都是排在了第5位，而在GATK的注释结果就排到了第13位，少了一个突变位点。下面通过lollipopPlot函数来进一步比较一下：

## lollipopPlot for SP140
gene='SP140'
protein=c("AAChange.refGene","Protein_Change","Protein_Change")
for (i in 1:3) {
  #i=3
  png(paste0(gene,'_', anno[i], '.png'),res = 150,width = 1080,height = 1080)
  maftools::lollipopPlot(maf = laml[[i]],
                         gene = gene,
                         AACol = protein[i],
                         labelPos = 'all')
  dev.off()

}

可以看到GATK比另外两个软件的注释结果少了一个位点，且三个软件的氨基酸的坐标也不完全一致。我猜是注释到不同转录本的原因，结合IGV进行查看一下，果然如此，GATK选择的转录本是NM_001278452.1，在IGV中这个转录本该位点230238870为内含子，同样导致氨基酸坐标不一致也是因为注释到不同的转录本：

最后我们再来比较TP53基因的注释结果，也是因为转录本的不同导致氨基酸坐标不一致：

## lollipopPlot for TP53
gene='TP53'
protein=c("AAChange.refGene","Protein_Change","Protein_Change")
for (i in 1:3) {
  #i=3
  png(paste0(gene,'_', anno[i], '.png'),res = 150,width = 1080,height = 1080)
  maftools::lollipopPlot(maf = laml[[i]],
                         gene = gene,
                         AACol = protein[i],
                         labelPos = 'all')
  dev.off()
}

maftools还可以绘制很多图，完成更多的分析。感兴趣的朋友可以试试，下一节我们将开始拷贝数变异的分析。

最后友情宣传生信技能树

• 全国巡讲:R基础，Linux基础和RNA-seq实战演练 : 预告：12月28-30长沙站