TCGA+biomarker——Cox回归森林图

森林图简介

森林图（forest plot），从定义上讲，它一般是在平面直角坐标系中，以一条垂直于X轴的无效线（通常坐标X=1或0）为中心，用若干条平行于X轴的线段，来表示每个研究的效应量大小及其95%可信区间，并用一个棱形来表示多个研究合并的效应量及可信区间，它是Meta分析中最常用的结果综合表达形式，现在也广泛应用在biomarker此类研究中。

森林图的科研用途

提到森林图，很多人的第一反应就是Meta分析。实际上，除了Meta分析，森林图还有诸多用处。森林图可以直观的反映出效应量（例如RR、OR、HR或者WMD）大小及其95% CI，这些效应量指标通常都是通过采用多因素回归分析所得，因此我们同样可以把森林图借鉴过来，用于展示多因素回归分析的结果。总结来说，森林图的科研用途主要用于Meta和临床实验。

临床实验普通分析（重点）——常规森林图

下图就是常规Cox回归结果的森林图展示，主要体现了变量、病人数量、P值和HR值，同时在每个变量内部设置reference（便于直观比较，变量为字符串-因子化属性）。比如: 年纪older变量位于无效线（即中间的那条竖线）右侧，说明年老有助于死亡。森林图在常规情况下事件结局是"生/死"这种两分类，但有时候事件结局是"有效/无效"、"治疗/未治疗"等等其他二分类情况，评估事件是好事还是坏事。比如生存（生：0；死：1，数值小的默认为事件结局），就是好事件，当位于无效线左侧的变量，说明这些变量有助于事件发生，是保护因素；位于无效线右侧的变量，说明这些变量不利于事件发生，是危险因素；当与无效线相交时，说明这些变量与事件发生之间关系不强！在整体数据上，用来评估这些变量因素对事件结局的影响！

image

PS：若是不想在多分类变量中设置refence，变量字符串形式即可，不要因子化。

image

临床实验分组分析（重点）——分组森林图

很多时候，文章展现的是亚组森林图。如下图：生存为事件结局，添加阿司匹林和阿哌沙班两个用药组。当位于无效线左侧的变量，可认为阿司匹林用药组的结局事件发生率大于阿哌沙班用药组；位于无效线右侧的变量，说明阿司匹林用药组的结局事件发生率大于阿哌沙班用药组。从整体数据来看，阿司匹林效果比阿哌沙班好！

image

森林图的结果解释

1. 普通森林图解

image

图列解释：

第1列为多个分类变量；

第2列为各个变量下样本数量：以stage为例，stage_I数量为100，stage_II为100，stage_III为100；

第3列多因素cox回归的统计学P值，

第4列为HR及其置信区间的图形表达，图中中间的竖线为无效线，即OR=1，表示所研究因素和结局无统计学关联；每条横线为该研究的95%可信区间，横线中央的小方块(或其他图形)为OR值的点估计。若某个研究95%可信区间的横线与无效竖线有交叉，可认为所研究因素和结局无统计学关联；若该横线落在无效竖线的左侧，可认为所研究因素有利于结局的发生；若该横线落在无效线的右侧，可认为所研究因素不利于结局的发生。

第5列为HR及其置信区间的数值化表示。

图行解释：

以stage为例，其中stage_I为参照，stage_II和stage_III都位于无效线的右侧，说明不利于生存。与stage_I相比，stage_II患者是stage_I患者死亡风险的2.26倍，他的上下可信区间是1.78_{2.58，并且经过统计是P=0.01有意义；同时stage_III患者是stage_I患者死亡风险的1.86倍，他的上下可信区间是1.58}1.98，并且经过统计是P=0.02也有意义。

2. 分组森林图解(评估多个分类变量对S-1+Doc和S-1两种治疗方案的作用，结局事件：PFS)

image

图列解释：

第1列为多个分类变量；

第2列和第3列表示S-1+Doc和S-1这两组的基本情况：如“87/454”，前者为发生某事件的例数、后者为各组样本量；

第5列为HR及其置信区间的数值化表示；

第6列为交互检测统计学P值，按照不同的变量分层分析，交互作用的P值不显著，表明不同分层结果一致，文章的结果可靠。

图行解释：

首行是总人群Total的分析结果，HR=0.632，表明联合治疗组与S-1单药相比，其复发风险降低37%；此时横线线段不与无效线相交，提示联合治疗组与单药治疗组总体人群中的差异有统计学意义。

下面都是分层分析的结果，在不同的性别、年龄、分期等多个亚组人群中RFS(无复发生存)的变化情况。以性别为例：男性（Male）中的HR为0.605，表明在男性胃癌患者接受联合治疗与接受单独S-1治疗相比，其复发风险降低了39.5%，此时横线线段不与无效线相交，也提示联合治疗组与单药治疗组在男性这组人群中的差异有统计学意义。女性人群中的HR=0.704，表明在女性胃癌患者接受联合治疗与接受单独S-1治疗相比，其复发的风险降低了29.6%，此时横线线段与无效线相交，提示联合治疗组与单药治疗组在女性这组人群中的差异无统计学意义。最后出现的交互作用P值为0.6058，大于0.05，说明男性和女性中药物两组方案的药物治疗与疾病复发关系的差异并不明显，也就是说明男性HR（0.605）和女性的HR（0.705）差异不显著，不能认为药物治疗和复发的关系在不同性别中不一样。

如何绘制森林图？

第一种情况: 单因素cox回归结果

针对单个变量，做cox回归绘制森林图，ggforest——直接用单因素cox回归返回的结果做森林图，ggforest(res.cox, data = lung) ，以cox回归结果和数据源两个参数即可，这种方式森林图主要只展现变量、样本数量、置信区间及线图（样式比较局限）。

image

针对单个变量（多分类），如果结果只想显示整体，那么将该多分类变量数值化表示；如果想协变量都展示，那么将该多分类变量字符化，并转化为因子，首位因子被默认为参照！

image

针对多个变量****（****常见形式****）的单因素回归分析做森林图，或者想在森林图中添加其他信息时需要将多个变量的结果进行汇总，整理成表格（或其他格式），然后采用forestplot函数做森林图，整理数据结构如下：变量名、样本数量、占比等信息源自临床信息汇总整理；置信区间和P值源自cox回归结果。
image
针对特殊情况——亚组森林图，同样也需要将所有信息整理成表格形式，然后用forestplot绘图。变量名、样本数量、占比、某因素条件分组（如下图的FOLFOX和FOLFOX plus SIRT两种条件，其下数据可以是样本数量、事件发生率等）源自临床信息汇总整理；置信区间源自cox回归结果；显著性视情况，有文章不放也有文章放（常见两种P value和P interaction）。分组森林图的做法，就是在固定的某种因素条件下，探索各种变量的分组条件下的单因素cox回归，如年龄：不同年龄亚组；如性别：男女不同亚组。简单点理解：其实就是在各种不同变量的亚群数据集中，探索某个固定因素对结局事件的影响。

image

第二种情况：多因素cox回归森林图

针对多因素cox回归，满足普通需求，用ggforest即可（只能展示变量、样本数量、置信区间及线图这几个信息）

image

对于多分类变量，如果结果只想显示整体，那么将该多分类变量数值化表示；如果想协变量都展示，那么将该多分类变量字符化，并转化为因子，首位因子被默认为参照！（如sex）

image

具体绘制方法：

1. ggforest函数绘制森林图（以1个单因素cox为例）

##------------------------sex整体展示-------------------
rm(list = ls())
library("survival")
library("survminer")
#载入并查看数据集
data("lung")
head(lung)
str(lung)#该数据将所有变量都转换为数值型，包括性别(1,2表示)，分期(1,2,3,4表示)等。若是字符型的话，结果会有所不同！
res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
res.cox
summary(res.cox)
#绘制森林图(ggforest),绘制单个变量的单因素cox森林图
ggforest(res.cox, data = lung,  #数据集
         main = 'Hazard ratio',  #标题
         cpositions = c(0.05, 0.15, 0.35),  #前三列距离
         fontsize = 1, #字体大小
         noDigits = 3 #保留HR值以及95%CI的小数位数
)

#-------------------------sex分类分开展示--------------------
rm(list = ls())
library("survival")
library("survminer")
#载入并查看数据集
data("lung")
head(lung)
str(lung)#该数据都为数值型，如性别(1,2表示).要分类展示cox回归需要将分类变量因子化
#将分类变量从数值型改为因子
lung <- within(lung, {
  sex <- factor(sex, labels = c('female', 'male'))})  # female为对照
res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
res.cox
summary(res.cox)
#绘制森林图(ggforest),多分类变量分开展示
ggforest(res.cox, data = lung,  #数据集
         main = 'Hazard ratio',  #标题
         cpositions = c(0.05, 0.15, 0.35),  #前三列距离
         fontsize = 1, #字体大小
         noDigits = 3 #保留HR值以及95%CI的小数位数
) 
##是否分开展示看结果，整体展示和分开展示哪个好说明好解释用哪个！

2. ggforest函数绘制森林图（以多因素cox为例）

##----------------------适用于数值型变量,即每个变量cox结果只有一行-------------------
rm(list = ls())
library("survival")
library("survminer")
#载入并查看数据集
data("lung")
head(lung)
str(lung)#该数据将所有变量都转换为数值型
#cox 回归分析
res.cox <- coxph(Surv(time, status) ~ age+sex+ph.ecog+ph.karno+pat.karno, data = lung)
res.cox
summary(res.cox)
#绘制森林图(ggforest)
ggforest(res.cox, data = lung,  #数据集
         main = 'Hazard ratio',  #标题
         cpositions = c(0.05, 0.15, 0.35),  #前三列距离
         fontsize = 1, #字体大小
         noDigits = 3 #保留HR值以及95%CI的小数位数
)

#------------------适用于多分类变量，每个变量下的子变量都对应一行数据----------------------

rm(list = ls())
library("survival")
library("survminer")
#载入并查看数据集
data("lung")
head(lung)
str(lung)#该数据都为数值型，如性别(1,2表示). 要分类展示cox回归需要将分类变量因子化
#将分类变量sex从数值型改为因子
lung <- within(lung, {
  sex <- factor(sex, labels = c('female', 'male'))  # female为对照
})
#cox回归分析
res.cox <- coxph(Surv(time, status) ~ age+sex+ph.ecog+ph.karno+pat.karno, data = lung)
res.cox
summary(res.cox)
#绘制森林图(ggforest),多分类分开展示
ggforest(res.cox, data = lung,  #数据集
         main = 'Hazard ratio',  #标题
         cpositions = c(0.05, 0.15, 0.35),  #前三列距离
         fontsize = 1, #字体大小
         refLabel = '1', #相对变量的数值标签，也可改为1
         noDigits = 3 #保留HR值以及95%CI的小数位数
)

##是否分开展示看结果，整体展示和分开展示哪个好说明好解释用哪个！！

3. Forestplot函数绘制森林图（多个单因素cox为例）

Forestplot R包绘制森林图更灵活，可根据自身需求绘制。它主要的难点在于这个R包依赖于已经整理好的结果，因此需要提前整理好相关数据。绘制森林图的数据包括两部分：
1. 图中的文字部分：变量、样本数及占比、分组条件下各自事件发生数量或比例等，这些需要自行整理！
2. 置信区间线：数据源自cox回归结果中的HR值、上限值、下限值以及P值
大概数据要求如下：

image

具体绘制代码可参考：R-forestplot包| HR结果绘制森林图


rm(list = ls())
#载入R函数包
library(survival)
library(survminer)
library(forestplot)
#需要将多个自变量的单因素cox回归结果进行整理
#森林图的数据包括两部分：
##1.图中的文字部分：变量、样本数及占比、分组条件下各自事件发生数量或比例等，这些需要自行整理！
##2. 置信区间线：数据源自cox回归结果中的HR值、上限值、下限值以及P值
##以整理好的uni_data为例
data <- read.csv("uni_data.csv", stringsAsFactors=FALSE)
#查看数据
head(data) #数据中包含HR值、上限值、下限值以及P值，以及分组信息
#对数据进行部分修改，方便行名和列名字输出：
## 构建tabletext，更改列名称，展示更多信息：
np <- ifelse(!is.na(data$Count), paste(data$Count," (",data$Percent,")",sep=""), NA)
tabletext <- cbind(c("Subgroup","\n",data$Variable),
                   c("No. of Patients (%)","\n",np),
                   c("4-Yr Cum. Event Rate\n PCI","\n",data$PCI.Group),
                   c("4-Yr Cum. Event Rate\n Medical Therapy","\n",data$Medical.Therapy.Group),
                   c("P Value","\n",data$P.Value))

##绘制森林图：
forestplot(labeltext=tabletext, graph.pos=5,
           mean=c(NA,NA,data$Point.Estimate),
           lower=c(NA,NA,data$Low), upper=c(NA,NA,data$High),
           title="Hazard Ratio",##定义x轴
           xlab="    <---PCI Better---   ---Medical Therapy Better--->",
##根据亚组的位置，设置线型，宽度造成“区块感”：
           hrzl_lines=list("3" = gpar(lwd=1, col="#99999922"),
"7" = gpar(lwd=60, lineend="butt", columns=c(2:6), col="#99999922"),
"15" = gpar(lwd=60, lineend="butt", columns=c(2:6), col="#99999922"),
"23" = gpar(lwd=60, lineend="butt", columns=c(2:6), col="#99999922"),
"31" = gpar(lwd=60, lineend="butt", columns=c(2:6), col="#99999922")),
#fpTxtGp函数中的cex参数设置各个组件的大小：
           txt_gp=fpTxtGp(label=gpar(cex=1.25),ticks=gpar(cex=1.1),xlab=gpar(cex = 1.2),title=gpar(cex = 1.2)),
##fpColors函数设置颜色：
           col=fpColors(box="#1c61b6", lines="#1c61b6", zero = "gray50",summary = "red"),
#箱线图中基准线的位置：
           zero=1,
           cex=0.9, lineheight = "auto",
           colgap=unit(8,"mm"), #列间距：
#误差线的宽度，箱子大小：
           lwd.ci=2, boxsize=0.5,
#箱线图两端添加小竖线及高度：
           ci.vertices=TRUE, ci.vertices.height = 0.4)

往期回顾
TCGA+biomarker——常见结果展示
 TCGA+biomarker——Sample基线表
 TCGA+biomarker——单因素Cox回归

更多内容可关注公共号“YJY技能修炼”~~~