基因突变瀑布图的绘制

本期主要介绍使用R的GenVisR包绘制突变瀑布图。

要求输入的数据格式：

  GenVisR包支持现成的数据格式包括从突变注释工具VEP生成的MAF格式(TCGA使用的数据格式)，和从基因组建模系统GMS生成的MGI格式，以及自定义的数据框。所以如果有MAF或者 MGI格式的数据可以直接导入进行作图，比如从TCGA下载的突变位点数据就可以直接使用。但我们一般得到的是VCF格式的数据，可以使用SnpSift从VCF文件中提取指定列进行作图。

比如我们从下面代码得到下图数据框：

library(GenVisR)
col_brcaMAF<-c("Tumor_Sample_Barcode","Hugo_Symbol","Variant_Classification","amino_acid_change_WU","Chromosome")
brcaMAF_sub<-brcaMAF[,col_brcaMAF]

第一列样品名称，第二列基因名，第三列突变位点类别，第四列是氨基酸突变HGVS标签，第五列是染色体名，其实只要有前三列就可以作图了，注意列名和MAF格式相同。

根据上述数据作图：

#设置该突变在样本中出现频率的阈值，低于此频率的将被过滤掉
waterfall(brcaMAF_sub, fileType = "MAF", mainRecurCutoff = 0.5)

该图最上面的panel表示mutations in sample/coverage space∗1,000,000，也被称为Mutation Burden，这个coverage space可以自定义（默认是44100000），这里规定每个样本的coverage space是相同的，如果需要的话可以自定义每个样本的coverage space进行计算，不想显示这部分的话设置plotMutBurden = FALSE。左侧panel表示每个突变位点在样本中出现的频率，该值为0.5表示有一般样本携带此突变。中间panel表示每个样本每个基因携带的突变信息，颜色表示该样本指定基因的突变类型。

只显示指定的基因：

waterfall(brcaMAF_sub, fileType = "MAF", plotGenes =c("CDKN1B", "PIK3CA"))

只显示指定的样本：

samples_to_plot <- c("TCGA-A1-A0SO-01A-22D-A099-09", "TCGA-A2-A0EU-01A-22W-A071-09", "TCGA-A2-A0ER-01A-21W-A050-09", "TCGA-A1-A0SI-01A-11D-A142-09", "TCGA-A2-A0D0-01A-11W-A019-09")
waterfall(brcaMAF_sub, plotSamples = samples_to_plot, mainRecurCutoff = 0.25)

添加临床信息显示在下方：

# Create clinical data
subtype <- c("lumA", "lumB", "her2", "basal", "normal")
subtype <- sample(subtype, 50, replace = TRUE)
age <- c("20-30", "31-50", "51-60", "61+")
age <- sample(age, 50, replace = TRUE)
sample <- as.character(unique(brcaMAF$Tumor_Sample_Barcode))
clinical <- as.data.frame(cbind(sample, subtype, age))

# Melt the clinical data into 'long' format.
library(reshape2)
clinical <- melt(clinical, id.vars = c("sample"))

临床信息数据如下：

waterfall(brcaMAF_sub, clinDat = clinical, clinVarCol = c(lumA = "blue4", lumB = "deepskyblue", 
    her2 = "hotpink2", basal = "firebrick2", normal = "green4", `20-30` = "#ddd1e7", 
    `31-50` = "#bba3d0", `51-60` = "#9975b9", `61+` = "#7647a2"), mainRecurCutoff = 0.05, 
    maxGenes = 25, clinLegCol = 2, clinVarOrder = c("lumA", "lumB", "her2", "basal", 
        "normal", "20-30", "31-50", "51-60", "61+"))

方格中显示染色体：

waterfall(brcaMAF, mainRecurCutoff = 0.05, maxGenes = 10, mainLabelCol = "Chromosome")

方格中显示氨基酸替换：

显示样本名：

更多原创精彩视频敬请关注生信杂谈: