生信代码：绘制基因组突变全景图

科研菌

发布于 2021-02-19 10:14:41

5.9K0

发布于 2021-02-19 10:14:41

文章被收录于专栏：科研菌

对于基因组突变全景图相信大家并不陌生，它是基因组学突变数据最基本的可视化展示方法之一。一张漂亮的，高大上的基因突变全景图不仅能展示出丰富的信息，还能为你的文章增色不少，其绘制方法也多种多样。今天我们则来看看最常用的两个包maftools和ComplexHeatmap在绘制基因组突变全景图上的异同。首先让我们来简单的了解下这两个包：

1.maftools包

maftools是一个特别强大的包，其基于MAF格式可以做很多分析以及可视化工作。

1.1 安装、加载包

if (!require("BiocManager"))
  install.packages("BiocManager")
BiocManager::install("maftools")
library(maftools)

1.2 数据格式要求

maftools包需要的数据格式为突变数据的maf文件（突变注释格式）以及临床信息文件，可以直接从官网下载：https://portal.gdc.cancer.gov/，也可以用TCGAbiolinks包来下载（具体下载方式及处理详见公众号笔记：TCGAbiolinks 学习笔记——汇总贴），当然从cBioProtal或者Xena等网站下载也可以，方法多种，任凭大家选择。本次以“KICH“数据示例，以下为下载好的数据： #突变的maf文件：

#下载并处理好的临床信息：

1.3 读取数据--read.maf

laml<-read.maf("TCGA.KICH.mutect.somatic.maf")
#看下maf数据
maf<-laml@data

#这里用到了read.maf()函数，我们来看下它的详细用法及参数意义吧

用法：

read.maf(maf,clinicalData = NULL,
removeDuplicatedVariants = TRUE,
useAll = TRUE,
gisticAllLesionsFile = NULL,
gisticAmpGenesFile = NULL,
gisticDelGenesFile = NULL,
gisticScoresFile = NULL,
cnLevel = "all",
cnTable = NULL,
isTCGA = FALSE,
vc_nonSyn = NULL,
verbose = TRUE)

参数及用法

主要参数	用法
maf	制表符分隔的MAF文件。文件也可以是gz压缩的，你也可以将已读取的MAF文件作为数据框提供
clinicalData	与MAF中每个样本/ Tumor_Sample_Barcode相关的临床数据。可以是文本文件或data.frame。默认为NULL
isTCGA	是来自TCGA源的输入MAF文件。如果TRUE仅使用Tumor_Sample_Barcode中的前12个字符
其余参数	详见官网说明文档

1.4 匹配临床信息（便于注释）

laml <-read.maf(maf=maf,isTCGA=TRUE,
                clinicalData = 'KICH-clinical.txt')

1.5 画突变全景图--oncoplot（）

画瀑布图的函数主要为oncoplot，先来了解下该函数的用法和参数,该函数参数很多，我们主要讲解其主要参数。

参数	用法
maf	则为你需要画图的数据
top	显示多少基因
genes	这个参数用来画指定的基因，默认为：NULL
altered	默认值为FALSE。根据突变状态选择最重要的基因。如果为TRUE，则选择基于顶级基因的变异（CNV或突变）
significant	基因和对应的q值作为侧边栏。默认为NULL
mutsig	Mutsig resuts，通常提供名sig_genes.txt的文件
includeColBarCN	是否在柱状图中包含CN
logColBar	在log10比例尺上绘制顶部条形图。默认值为FALSE。
ClinicalFeatures	要在图中绘制的MAF的“clinical.data”中的列名称， Dafault为NULL
annotationDat	如果读取的MAF文件没有临床数据，需提供一个自定义data.frame，该行的Tumor_Sample_Barcode列包含样品名称以及其余带有注释的列
genesToIgnore	不显示的基因，默认为：NULL
removeNonMutated	逻辑值。如果为TRUE，将删除在copcoplot中没有突变的样本，以实现更好的可视化。默认为TRUE
colors	每个Variant_Classification的颜色命名向量
sortByMutation	根据突变强制排序矩阵。默认值为FALSE
sortByAnnotation	通过提供的“clinicalFeatures”对comcomtrix（样本）进行逻辑排序。根据第一个“ clinicalFeatures”进行排序。默认为FALSE
draw_titv	如果为TRUE，则绘制TiTv plot，默认为：FALSE
writeMatrix	将用于生成图的突变矩阵矩阵写出
bgCol	野生型（未突变）样品的背景网格颜色。默认灰色-“＃CCCCCC”
borderCol	边界网格的颜色（未突变）样本。默认为“白色”
其他参数	详见官网说明文档

开始画图

#显示前20各基因，默认参数画图
oncoplot(maf = laml,top = 20)

#添加临床注释信息并根据注释排序
oncoplot(maf = laml,top = 20,clinicalFeatures
         =c("gender",'stage'),sortByAnnotation = T)

#添加titv图
oncoplot(maf = laml,top = 20,clinicalFeatures
         =c("gender",'stage'),sortByAnnotation = T,
         draw_titv = TRUE)

#设置边界不显示（建议大样本量时用）
oncoplot(maf = laml,top = 20,clinicalFeatures
         =c("gender",'stage'),sortByAnnotation = T,
         draw_titv = TRUE,borderCol=NULL)

以上为基本画图，想要调整字体大小，颜色，注释等可根据参数自由调整。
画基因突变全景图只是强大的mafools包的功能之一，如果有兴趣的小伙伴可阅读官网说明文档。

2.ComplexHeatmap包

maftools主要基于maf文件，但是有时候我们的文件并不是maf文件，那该如何画图呢？ComplexHeatmap包是个很不错的选择。（当然你也考虑将数据可以转换为maf后再用maftools包绘图）
ComplexHeatmap包也是一个超级强大的包，函数功能很多，今天则主要讲解该包如何绘制基因组突变全景图。

2.1 安装加载包

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("ComplexHeatmap")
library(ComplexHeatmap)

2.2 需要的数据格式

需要的数据格式为突变信息矩阵（这里我们依旧以KICH为例）---不管什么数据整理为该矩阵就可以画图了。

2.3 读入数据、匹配颜色

# 读入数据
library(ComplexHeatmap)
mat<-read.table(file="onco_matrix_KICH.txt",
                 header = T, stringsAsFactors = FALSE,sep = "\t")
mat[is.na(mat)]<-""
mat[1:3, 1:3]
#          TCGA.KO.8408      TCGA.KN.8428 TCGA.KL.8335
# TP53   Frame_Shift_Del Nonsense_Mutation  Splice_Site
# PTEN Nonsense_Mutation                             
# ZAN  Nonsense_Mutation Missense_Mutation

设置对应突变的颜色，背景设置为灰色

col <- c("Nonsense_Mutation" = "blue", "Missense_Mutation" = "red", 
        "Splice_Site" = "#008000","Frame_Shift_Del"="orange",
        "Frame_Shift_Ins"="black","Multi_Hit"="purple")
alter_fun <- list(
  background = function(x, y, w, h) {
    grid.rect(x, y, w-unit(0.5, "mm"), h-unit(0.5, "mm"), 
              gp = gpar(fill = "#CCCCCC", col = NA))
  },
  # big blue
  Nonsense_Mutation = function(x, y, w, h) {
    grid.rect(x, y, w-unit(0.5, "mm"), h-unit(0.5, "mm"), 
              gp = gpar(fill = col["Nonsense_Mutation"], col = NA))
  },
  # bug red
  Missense_Mutation = function(x, y, w, h) {
    grid.rect(x, y, w-unit(0.5, "mm"), h-unit(0.5, "mm"), 
              gp = gpar(fill = col["Missense_Mutation"], col = NA))
  },
  # small green
  Splice_Site = function(x, y, w, h) {
    grid.rect(x, y, w-unit(0.5, "mm"), h*0.33, 
              gp = gpar(fill = col["Splice_Site"], col = NA))
  },
  Frame_Shift_Del = function(x, y, w, h) {
    grid.rect(x, y, w-unit(0.5, "mm"), h*0.33, 
              gp = gpar(fill = col["Frame_Shift_Del"], col = NA))
  },
  Frame_Shift_Ins = function(x, y, w, h) {
    grid.rect(x, y, w-unit(0.5, "mm"), h*0.33, 
              gp = gpar(fill = col["Frame_Shift_Ins"], col = NA))
  },
  Multi_Hit = function(x, y, w, h) {
    grid.rect(x, y, w-unit(0.5, "mm"), h*0.33, 
              gp = gpar(fill = col["Multi_Hit"], col = NA))
  }
)

设置标签和标题,六个突变标签

column_title <- "OncoPrint for KICH"
heatmap_legend_param <- list(title = "Alternations", at = c("Nonsense_Mutation" , "Missense_Mutation", 
                                                            "Splice_Site" ,"Frame_Shift_Del",
                                                            "Frame_Shift_Ins","Multi_Hit"), 
                            labels = c("Nonsense_Mutation" , "Missense_Mutation", 
                                       "Splice_Site" ,"Frame_Shift_Del",
                                       "Frame_Shift_Ins","Multi_Hit"))

2.4画图

ComplexHeatmap包画图的主要函数为oncoprint(),先让我们来简单了解下这个函数

主要参数	用法
mat	用于画图的矩阵
get_type	如果在矩阵中将不同的突变编码为复杂的字符串，则此自定义函数将确定如何提取它们。仅当mat是矩阵时才有效，默认值为default_get_type
alter_fun	可以自定义不同的变异通过什么样子来进行显示。这个函数包括四个参数x,y,w,h分别代表变异的位置(x,y)以及高度(h)和宽度(w)
alter_fun_is_vectorized	是否将alter_fun实现矢量化
show_pct	是否在左侧显示百分比值
pct_gp	百分比值的图形参数
show_column_names	可以用来定义是否显示列名
row_names_side	定义行名的位置
pct_side	定义突变百分比的位置
anno_oncoprint_barplot	调整上面和有面barplot的具体参数
heatmap_legend_param	定义图例的变化
其他参数	详见官网说明

在初步了解各参数意义后，让我们一起画图吧

#画图并去除无突变的样本和基因
oncoPrint(mat,
          alter_fun = alter_fun, col = col, 
          column_title = column_title,
          remove_empty_columns = TRUE,
          remove_empty_rows = TRUE,
          heatmap_legend_param = heatmap_legend_param)

#对行显示位置进行调整
oncoPrint(mat,
          alter_fun = alter_fun, col = col, 
          remove_empty_columns = TRUE, remove_empty_rows = TRUE,
          pct_side = "right", row_names_side = "left",
          column_title = column_title, heatmap_legend_param = heatmap_legend_param)

颜色、字体、位置等修改可以参考参数自行设置

#添加临床信息的注释
pdata<-read.csv(file="KICH-clininc.csv",row.names = 1)
sOrder <- order(pdata$stage,
                pdata$gender,
                pdata$vital_status)
pdata <- pdata[sOrder,]
mat <- mat[, sOrder]
ha<-HeatmapAnnotation(stage=pdata$stage,
                      gender=pdata$gender,
                      vital_status=pdata$vital_status,
                      show_annotation_name = TRUE,
                      annotation_name_gp = gpar(fontsize = 7))
oncoPrint(mat,top_annotation = ha,
          alter_fun = alter_fun, col = col, 
          column_title = column_title, heatmap_legend_param = heatmap_legend_param)

#还可以添加柱状图/散点图（这里以柱状图为例子）
ha<-HeatmapAnnotation(stage=pdata$stage,
                      gender=pdata$gender,
                      vital_status=pdata$vital_status,
                      age= anno_barplot(pdata$age_at_initial_pathologic_diagnosis,
                                        border=FALSE, gp = gpar(fill="blue")),
                      show_annotation_name = TRUE,
                      annotation_name_gp = gpar(fontsize = 7))
oncoPrint(mat,top_annotation = ha,
          alter_fun = alter_fun,col=col,
          column_title = column_title, heatmap_legend_param = heatmap_legend_param)

#添加突变注释（可以只展示你想展示的类别）
oncoPrint(mat,
          alter_fun = alter_fun,col=col,
          column_title = column_title, heatmap_legend_param = heatmap_legend_param,
          right_annotation = rowAnnotation(
            row_barplot = anno_oncoprint_barplot(c("Nonsense_Mutation", "Missense_Mutation",
                                                   "Splice_Site"), 
                                                 border = TRUE, height = unit(4, "cm"), 
                                                 axis_param = list(side = "bottom", labels_rot = 90))))

#还可以添加热图、密度图等（这里以热图为例）
#热图数据为随机产生
ht_list = oncoPrint(mat,
                    alter_fun = alter_fun, col = col, row_barplot = anno_oncoprint_barplot(c("Nonsense_Mutation", "Missense_Mutation",
                                                                                             "Splice_Site"), 
                                                                                           border = TRUE, height = unit(2, "cm"), 
                                                                                           axis_param = list(side = "bottom", labels_rot = 90))),
                    column_title = column_title, heatmap_legend_param = heatmap_legend_param) +
  Heatmap(matrix(rnorm(nrow(mat)*5), ncol = 5), name = "expr", width = unit(2, "cm"))
draw(ht_list)

当然，ComplexHeatmap包的功能远远不止以上这些，标签的位置、颜色、字体的大小、位置等全都是可以根据相对应的参数自己灵活改变的。

3.总结

以上，我们简单介绍了画突变全景图的两个包，现在让我们来对比以下这两个包的优缺点吧。

	maftools	ComplexHeatmap
需要的文件格式	maf格式的文件	突变矩阵
画图主要函数	oncoplot（）	oncoPrint（）
代码	相对简单	相对复杂（需要一定的R基础）
灵活度	相对固定	相当灵活，想画哪里写哪里
出图	相对单调（因为是对maf文件画图）	自定义的程度相对更高，可以加入更多元素

同时，小编在学习过程了，也总结了点小窍门，希望能帮助大家少踩雷:

虽然maftools处理的是maf文件，但是不要慌，txt、csv、xlx等格式的文件和maf格式文件是可以相互转换的。（以下这段代码是用cBioProtal网站txt的突变文件转变为maf后，用maftools画图，完全可行的）

不管是使用maftools包还是ComplexHeatmap包，突变信息的Tumor_Sample_Barcode是完整的TCGA代码16位，但是临床信息的Tumor_Sample_Barcode的只有前12位，所以在匹配临床信息时都要记得使用substr(data$Tumor_Sample_Barcode,1,12)先处理条码的不同哦。
虽然ComplexHeatmap包对代码要求相对高一些，不要慌，在此安利大家一个超详细的讲解，网址：https://www.jianshu.com/p/3448ca155440

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-01-28，如有侵权请联系 cloudcommunity@tencent.com 删除

barcode