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

生信补给站

发布于 2024-04-24 15:52:01

3550

发布于 2024-04-24 15:52:01

当使用各种机器学习方法（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_replace
riskScore_cli2$pathologic_T2 <- str_replace(riskScore_cli2$pathologic_T2, "T1[a-d-c-d]?", "T1")
#str_detect
riskScore_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")
)