其他
肿瘤外显子数据处理系列教程（八）不同注释软件的比较(下)：可视化比较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还可以绘制很多图，完成更多的分析。感兴趣的朋友可以试试，下一节我们将开始拷贝数变异的分析。