当使用各种机器学习方法（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 44riskScore_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 34riskScore_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 9riskScore_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"))
以上就完成了风险得分和临床指标的热图，拿去发文章吧。
◆ ◆ ◆ ◆ ◆
精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）
RNAseq纯生信挖掘思路分享？不，主要是送你代码！（建议收藏）
