RNAseq | ComplexHeatmap绘制临床数据热图（所见即所得）

生信补给站 生信补给站 2024-04-19 09:00

当使用各种机器学习方法（CoxBoost，Lasso，SuperPC，randomForestSRC，Elastic Net等）完成预后模型后，除了在组学层面（ IPS评分，药物反应预测，WGCNA等）进行一系列的分析外，还可以将定义的风险得分和临床指标进行比较。如Expression of hub genes of endothelial cells in glioblastoma-A prognostic model for GBM patients integrating single-cell RNA sequencing and bulk RNA sequencing中下图所示

图片

最初我完成该图的方法是用含有基因表达的热图，然后截图或者PS成只有临床指标。这里介绍使用ComplexHeatmap直接完成该图。

一 载入R包，数据


使用前面系列推文的TCGA-SKCM的临床数据和随访数据，以及经过lasso模型计算的风险评分结果 。

后台回复“临床”既可获取Expr_phe_cli_riskscore.RData的示例数据。

library(tidyverse)library(ComplexHeatmap)library(ggsci) #颜色library(circlize) #连续颜色
#载入数据load("Expr_phe_cli_riskscore.RData")
riskScore_cli2 <- riskScore_cli %>%   inner_join(phe) %>%   column_to_rownames("sample") %>% select(- "_PATIENT")head(riskScore_cli2)

图片

图片

只要数据框中含有想展示的表型数据即可，一般会有风险得分，生存信息以及重要的临床指标，当然也可以其他重点关注的指标：（1）重点基因突变与否，例如KRAS突变 （2）某个CNV有无（3）TMB ，MSI，IDH等等你想展示的指标。如果添加基因表达量的话那就是正常的热图即可。

2，临床数据处理

在TCGA下载的临床数据需要进行一些处理，可以在excel中完成，当然也可以使用R完成。

包括但不限于以下（1）连续数值按照某个阈值转为分类 （2）向量和因子的转化 （3）将数据中的T1a ，T1b，T1 统一为T1期 类似的整理。

（1）和（2）比较简单，如下

#连续数值，按需转为分类riskScore_cli2$Age =ifelse(riskScore_cli2$age > 60,">60","<=60")
#字符串转为因子riskScore_cli2$OS <- as.factor(riskScore_cli2$OS)riskScore_cli2$tumor_stage <- as.factor(riskScore_cli2$tumor_stage)

（3）可以使用多种方式完成数据整理

A ：T分期使用直接指定的方法

注意%in% c("T1a","T1b","T1")的向量中要列出所有想转化的，假设有T1c的话 也需要加上。

table(riskScore_cli2$pathologic_T)# T0  T1 T1a T1b  T2 T2a T2b  T3 T3a T3b  T4 T4a T4b Tis  TX # 23  10  21  10  30  32  15  14  39  37  13  26 109   7  44
riskScore_cli2$pathologic_T[riskScore_cli2$pathologic_T %in% c("T1a","T1b","T1")] <- "T1"riskScore_cli2$pathologic_T[riskScore_cli2$pathologic_T %in% c("T2a","T2b","T2")] <- "T2"riskScore_cli2$pathologic_T[riskScore_cli2$pathologic_T %in% c("T3a","T3b","T3")] <- "T3"riskScore_cli2$pathologic_T[riskScore_cli2$pathologic_T %in% c("T4a","T4b","T4")] <- "T4"table(riskScore_cli2$pathologic_T)

B：N分期，使用gsub替换的方式

table(riskScore_cli2$pathologic_N)# N0  N1 N1a N1b  N2 N2a N2b N2c  N3  NX #226  17  19  37   6  13  21   9  56  34 
riskScore_cli2$pathologic_N <- gsub("N1[abc]?", "N1", riskScore_cli2$pathologic_N)riskScore_cli2$pathologic_N <- gsub("^N2.*", "N2", riskScore_cli2$pathologic_N)riskScore_cli2$pathologic_N <- gsub("^N3.*", "N3", riskScore_cli2$pathologic_N)

C：M分期，使用grepl的方法

table(riskScore_cli2$pathologic_M)# M0  M1 M1a M1b M1c #407   5   5   5   9 
riskScore_cli2$pathologic_M <- ifelse(grepl("^M1", riskScore_cli2$pathologic_M),                                       "M1", riskScore_cli2$pathologic_M)riskScore_cli2$pathologic_M <- ifelse(grepl("^M0", riskScore_cli2$pathologic_M),                                       "M0", riskScore_cli2$pathologic_M)

D：还可以使用str_replace 或者 str_detect等方法进行转化，这里示例展示一下，不运行不影响推文的后续操作 。

riskScore_cli2$pathologic_T2 <- riskScore_cli2$pathologic_T# str_replaceriskScore_cli2$pathologic_T2 <- str_replace(riskScore_cli2$pathologic_T2, "T1[a-d-c-d]?", "T1")#str_detectriskScore_cli2$pathologic_T2 <- ifelse(str_detect(riskScore_cli2$pathologic_T2, "^T1"), "T1", riskScore_cli2$pathologic_T2)

以上就完成了本次分析需要的数据处理部分。

二 临床指标热图可视化


1，直接绘制

使用ComplexHeatmap绘制临床数据注释图 ，重点在于构建一个和临床数据相同列的0矩阵 。

# 提取想展示的临床数据riskScore_cli2 <- riskScore_cli2 %>%   select(riskScore:tumor_stage,Age) %>%   select(- "age")# 构建列注释块ha=HeatmapAnnotation(df=riskScore_cli2)# 构建zero矩阵zero_row_mat=matrix(nrow=0, ncol=nrow(riskScore_cli2))#绘制热图Hm <- Heatmap(zero_row_mat, top_annotation=ha)Hm#调整legend的位置和大小draw(Hm, merge_legend = TRUE,      heatmap_legend_side = "bottom",      annotation_legend_side = "bottom",     width = unit(16, "cm"), height = unit(1, "cm")     )

图片

2，图形优化调整

上面可以顺利的完成图形可视化，相较文献还可以在（1）表型内容排序，比如优先Score高低排序，然后Stage排序（2）表型注释的顺序，比如先展示Score，然后OS，stage等 和（3）每种表型进行自定义的颜色设置 上进行优化和调整。

（1）表型内部排序  ，使用arrange 进行排序，可以依次选择多个指标

riskScore_cli3 <- riskScore_cli2 %>%  arrange(riskScore2,OS,tumor_stage,gender,OS.time,Age)

（2）和（3）一起在HeatmapAnnotation注释中解决，如果为省事未展示T M N分期 ，可以自行添加。

library(circlize)#连续性变量的颜色设置col_fun_time <- colorRamp2(  c(0, 3000, 11000),  #根据值的范围设置  c("#DC0000FF", "grey", "#1f78b4"))#ha <- HeatmapAnnotation(  Score = riskScore_cli3$riskScore2,  Stage = riskScore_cli3$tumor_stage ,  OS.Status = riskScore_cli3$OS,  OS.Time = riskScore_cli3$OS.time,  Gender = riskScore_cli3$gender ,  Age = riskScore_cli3$Age,  col = list(     Score = c("High" = "#BC3C29FF", "Low" = "#0072B5FF"),    OS.Status = c("0" = "#E18727FF", "1" = "#20854EFF"), #分类    OS.Time = col_fun_time , #连续    Gender = c("female" = "#AB3282", "male" = "#3A6963"),    Age = c("<=60" = "#712820", ">60" = "#E4C755"),    Stage = c("0" = "#E64B35FF", "1" = "#4DBBD5FF",              "2" = "#00A087FF", "3" = "#3C5488FF",              "4" = "#DC0000FF", "NA" = "#8491B4FF")      ))

可视化展示

Hm <- Heatmap(zero_row_mat, top_annotation=ha)draw(Hm, merge_legend = TRUE,      heatmap_legend_side = "bottom",      annotation_legend_side = "bottom",     width = unit(16, "cm"), height = unit(1, "cm"))

图片

以上就完成了风险得分和临床指标的热图，拿去发文章吧。


RNAseq12
生信文章常见图表63
RNAseq · 目录
上一篇RNAseq-ML | SuperPC 算法构建预后模型 并预测

微信扫一扫
关注该公众号

人划线