歡迎關(guān)注”生信修煉手冊(cè)”!

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

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

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

每一行表示一個(gè)變量，用errorbar展示該變量對(duì)應(yīng)的風(fēng)險(xiǎn)值的大小和置信區(qū)間，并將風(fēng)險(xiǎn)值和p值標(biāo)記在圖上。

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

> library(survminer) > require("survival") > model <- coxph( Surv(time, status) ~ sex + rx + adhere, + data = colon ) > model Call: coxph(formula = Surv(time, status) ~ sex + rx + adhere, data = colon)coef exp(coef) se(coef) z p sex -0.04615 0.95490 0.06609 -0.698 0.484994 rxLev -0.02724 0.97313 0.07690 -0.354 0.723211 rxLev+5FU -0.43723 0.64582 0.08395 -5.208 1.91e-07 adhere 0.29355 1.34118 0.08696 3.376 0.000736Likelihood ratio test=46.51 on 4 df, p=1.925e-09 n= 1858, number of events= 920 > ggforest(model)

效果圖如下

這種方式確實(shí)出圖簡(jiǎn)單，一步到位，但是提供的參數(shù)卻很少，靈活性很小，基本上沒(méi)法修改圖中的元素，另外一種方式，就是使用forest這個(gè)R包，這個(gè)R包靈活性很大，通過(guò)調(diào)參可以實(shí)現(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" + )

效果圖如下

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

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

> 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) + ) + )

效果圖如下

基于上述知識(shí)儲(chǔ)備和函數(shù)的幫助文檔，我們就可以實(shí)現(xiàn)和文章中圖片一致的效果圖了，只需要仔細(xì)鉆研函數(shù)的幫助文檔即可。

·end·

—如果喜歡，快分享給你的朋友們吧—

原創(chuàng)不易，歡迎收藏，點(diǎn)贊，轉(zhuǎn)發(fā)！生信知識(shí)浩瀚如海，在生信學(xué)習(xí)的道路上，讓我們一起并肩作戰(zhàn)！

本公眾號(hào)深耕耘生信領(lǐng)域多年，具有豐富的數(shù)據(jù)分析經(jīng)驗(yàn)，致力于提供真正有價(jià)值的數(shù)據(jù)分析服務(wù)，擅長(zhǎng)個(gè)性化分析，歡迎有需要的老師和同學(xué)前來(lái)咨詢(xún)。

? 更多精彩

• KEGG數(shù)據(jù)庫(kù)，除了pathway你還知道哪些

• 全網(wǎng)最完整的circos中文教程

• DNA甲基化數(shù)據(jù)分析專(zhuān)題

• 突變檢測(cè)數(shù)據(jù)分析專(zhuān)題

• mRNA數(shù)據(jù)分析專(zhuān)題

• lncRNA數(shù)據(jù)分析專(zhuān)題

• circRNA數(shù)據(jù)分析專(zhuān)題

• miRNA數(shù)據(jù)分析專(zhuān)題

• 單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析專(zhuān)題

• chip_seq數(shù)據(jù)分析專(zhuān)題

• Hi-C數(shù)據(jù)分析專(zhuān)題

• HLA數(shù)據(jù)分析專(zhuān)題

• TCGA腫瘤數(shù)據(jù)分析專(zhuān)題

• 基因組組裝數(shù)據(jù)分析專(zhuān)題

• CNV數(shù)據(jù)分析專(zhuān)題

• GWAS數(shù)據(jù)分析專(zhuān)題

• 機(jī)器學(xué)習(xí)專(zhuān)題

• 2018年推文合集

• 2019年推文合集

• 2020推文合集

? 寫(xiě)在最后

轉(zhuǎn)發(fā)本文至朋友圈，后臺(tái)私信截圖即可加入生信交流群，和小伙伴一起學(xué)習(xí)交流。

掃描下方二維碼，關(guān)注我們，解鎖更多精彩內(nèi)容！

一個(gè)只分享干貨的

生信公眾號(hào)

總結(jié)

以上是生活随笔為你收集整理的绘制cox生存分析结果的森林图的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。