【原】繪制cox生存分析結(jié)果的森林圖

生信修煉手冊 2022-05-30 發(fā)布于廣東

展開全文

在之前meta分析的文章中我們介紹了森林圖的畫法，典型的森林圖如下所示

每一行表示一個study，用errorbar展示log odds ratio值的分布，并將p值和m值標記在圖中。森林圖主要用于多個study的分析結(jié)果的匯總展示。

在構(gòu)建預后模型時，通常會先對所有基因進行單變量cox回歸，然后篩選其中顯著的基因進行多變量cox回歸來建模，對于cox回歸的結(jié)果，每個基因也都會有一hazard ratio和對應(yīng)的p值，也可以用森林圖的形式來展現(xiàn)，比如NAD+的文獻中就采用了這樣的一張森林圖

每一行表示一個變量，用errorbar展示該變量對應(yīng)的風險值的大小和置信區(qū)間，并將風險值和p值標記在圖上。

根據(jù)cox生存分析的結(jié)果繪制森林圖有多種方式，使用survminer包的ggforest函數(shù)，是最簡便的一種，代碼如下

> library(survminer)> require("survival")> model <- coxph( Surv(time, status) ~ sex + rx + adhere,+                 data = colon )> modelCall:coxph(formula = Surv(time, status) ~ sex + rx + adhere, data = colon)

              coef exp(coef) se(coef)      z        psex       -0.04615   0.95490  0.06609 -0.698 0.484994rxLev     -0.02724   0.97313  0.07690 -0.354 0.723211rxLev+5FU -0.43723   0.64582  0.08395 -5.208 1.91e-07adhere     0.29355   1.34118  0.08696  3.376 0.000736

Likelihood ratio test=46.51  on 4 df, p=1.925e-09n= 1858, number of events= 920> ggforest(model)

效果圖如下

這種方式確實出圖簡單，一步到位，但是提供的參數(shù)卻很少，靈活性很小，基本上沒法修改圖中的元素，另外一種方式，就是使用forest這個R包，這個R包靈活性很大，通過調(diào)參可以實現(xiàn)很多自定義效果，基本用法如下

> row_names <- list(list("test = 1", expression(test >= 2)))> test_data <- data.frame(+   coef = c(1.59, 1.24),+   low = c(1.4, 0.78),+   high = c(1.8, 1.55)+ )> forestplot(+   labeltext = row_names,+   mean = test_data$coef,+   lower = test_data$low,+   upper = test_data$high,+   zero = 1,+   cex  = 2,+   lineheight = "auto",+   xlab = "Lab axis txt"+ )

效果圖如下

雖然輸出很簡陋大，但是從基本用法可以看出，我們可以自定義變量名稱，指定風險值的大小，這樣我們只需要從cox回歸的結(jié)果中提取我們需要繪圖的元素進行繪制即可。

基本用法之外中添加的變量是單列注釋，如果要實現(xiàn)文獻中圖片的多列注釋效果，可以參考下面這個例子

> test_data <- data.frame(+   coef1 = c(1, 1.59, 1.3, 1.24),+   coef2 = c(1, 1.7, 1.4, 1.04),+   low1 = c(1, 1.3, 1.1, 0.99),+   low2 = c(1, 1.6, 1.2, 0.7),+   high1 = c(1, 1.94, 1.6, 1.55),+   high2 = c(1, 1.8, 1.55, 1.33)+ )>> col_no <- grep("coef", colnames(test_data))> row_names <- list(+   list("Category 1", "Category 2", "Category 3", expression(Category >= 4)),+   list(+     "ref",+     substitute(+       expression(bar(x) == val),+       list(val = round(rowMeans(test_data[2, col_no]), 2))+     ),+     substitute(+       expression(bar(x) == val),+       list(val = round(rowMeans(test_data[3, col_no]), 2))+     ),+     substitute(+       expression(bar(x) == val),+       list(val = round(rowMeans(test_data[4, col_no]), 2))+     )+   )+ )>> coef <- with(test_data, cbind(coef1, coef2))> low <- with(test_data, cbind(low1, low2))> high <- with(test_data, cbind(high1, high2))> forestplot(row_names, coef, low, high,+   title = "Cool study",+   zero = c(0.98, 1.02),+   grid = structure(c(2^-.5, 2^.5),+     gp = gpar(col = "steelblue", lty = 2)+   ),+   boxsize = 0.25,+   col = fpColors(+     box = c("royalblue", "gold"),+     line = c("darkblue", "orange"),+     summary = c("darkblue", "red")+   ),+   xlab = "The estimates",+   new_page = TRUE,+   legend = c("Treatment", "Placebo"),+   legend_args = fpLegend(+     pos = list("topright"),+     title = "Group",+     r = unit(.1, "snpc"),+     gp = gpar(col = "#CCCCCC", lwd = 1.5)+   )+ )