【原】RNAseq|批量单因素生存分析 + 绘制森林图

生存分析作为转录组文章中的VIP，太常见了，那么如何批量得到所有候选基因的单因素结果以及可视化结果呢？

本文将分别使用循环方式 和ezcox进行批量单基因生存分析，以及使用ggplot2 和forestplot绘制单因素生存分析森林图。

一 载入R包，数据

仍然使用之前处理过的TCGA的SKCM数据，此外需要读入生存数据和临床数据

library(tidyverse)library(openxlsx)library("survival")library("survminer")load("RNAseq.SKCM.RData")
#选取部分基因作示例data.mat <- t(expr[order(apply(expr, 1, mad), decreasing = T)[1:1000],])#读取生存数据surv <- read.table("TCGA-SKCM.survival.tsv",sep = "\t" , header = T,                   stringsAsFactors = FALSE ,                   check.names = FALSE)phe <- read.xlsx("cli.xlsx",sheet = 1)

真实项目中，data.mat部分可以是某些通路，功能的基因（免疫，铜死亡，铁死亡，细胞凋亡等），差异分析的基因，WGCNA找到的hub基因，总之可以是各种方法找到的候选基因。

二 批量单因素分生存分析

1，使用循环的方式进行分析
首先处理表达数据，注意基因名字的处理，tidyverse包非常值的狠狠学

module_exp <- as.data.frame(data.mat) %>%  rownames_to_column("sample") %>%   inner_join(surv) %>%  #添加生存数据  select(sample,OS,OS.time,`_PATIENT`,everything()) %>% #将生存列放到前面  select_all(~str_replace_all(., "-", "_"))  #基因ID不规范会报错,下划线替换-dim(module_exp)
#指定待分析的基因module_expr.cox <- module_expcovariates <- names(module_expr.cox[,5:ncol(module_expr.cox)])#构建单因素模型univ_formulas <- sapply(covariates,function(x) as.formula(paste('Surv(OS.time, OS)~', x)))univ_models <- lapply( univ_formulas, function(x){coxph(x, data = module_expr.cox)})# 提取结果univ_results <- lapply(univ_models,                       function(x){                          x <- summary(x)                         p.value<-signif(x$wald["pvalue"], digits=3)                         wald.test<-signif(x$wald["test"], digits=3)                         beta<-signif(x$coef[1], digits=3);#coeficient beta                         HR <-signif(x$coef[2], digits=3);#exp(beta)                         HR.confint.lower <- signif(x$conf.int[,"lower .95"], 3)                         HR.confint.upper <- signif(x$conf.int[,"upper .95"], 3)                         HR <- paste0(HR, " (",                                       HR.confint.lower, "-", HR.confint.upper, ")")                         res<-c(beta, HR, wald.test, p.value)                         names(res)<-c("beta", "HR (95% CI for HR)", "wald.test", "p.value")                         return(res)                       })res <- t(as.data.frame(univ_results, check.names = FALSE))

注：univ_results部分中的内容可以自定义，包括添加更多的信息，调整小数位数等；2，使用 ezcox 一行输出
ezcox是ShixiangWang大佬 开发的R包，一行代码输出所需结果 ，参考使用 ezcox 进行批量 Cox 模型处理 - 知乎 (zhihu.com)

对比两种方式，结果是一致的。

三 绘制森林图

对于单因素的结果，经常出现的可视化方式就是绘制森林图 。可以使用经典的forestplot-R包绘制（封装），或者使用ggplot2绘制（自由设置）。1 ，forestplot包绘制
‍forestplot绘制的关键就在于构建tabletext信息。因为基因太多，随机选择30个基因进行绘制

set.seed(1234)res_2_test <- res_2[sample(nrow(res_2), 30), ]
#构建tabletextsample <- as.data.frame(res_2_test)tabletext1<-as.character(sample[,1])tabletext2<-round(as.numeric(sample[,'HR']),5)tabletext3<-paste(round(as.numeric(sample[,'HR']),3),round(as.numeric(sample[,'lower_95']),2),sep="(")tabletext4<-paste(tabletext3,round(as.numeric(sample[,'upper_95']),2),sep="-")tabletext5<-paste0(tabletext4,sep=")")tabletext<-cbind(tabletext1,tabletext2,tabletext5)
forestplot(labeltext=tabletext, #文本信息             mean = round(sample[,'HR'],3),##HR值lower = round(sample[,"lower_95"],2),##95%置信区间upper = round(sample[,"upper_95"],2),#95%置信区间           boxsize = 0.8,##大小           graph.pos=4,#图在表中的列位置           graphwidth = unit(0.4,"npc"),#图在表中的宽度比例           fn.ci_norm="fpDrawDiamondCI",#box类型选择钻石,可以更改fpDrawNormalCI；fpDrawCircleCI等col=fpColors(box="steelblue", lines="black", zero = "black"),#颜色设置           lwd.ci=2,ci.vertices.height = 0.1,ci.vertices=TRUE,#置信区间用线宽、高、型           zero=1,#zero线横坐标           lwd.zero=2,#zero线宽           grid=T,           lwd.xaxis=2,#X轴线宽           title="Hazard Ratio",           xlab="",#X轴标题           clip=c(-Inf,3),#边界           colgap = unit(0.5,"cm")   )

更多参数设置可以使用 ??forestplot 查看 或者 R-forestplot包| HR结果绘制森林图

2， ggplot2 方式绘制自由度较高，需要对ggplot2有基本的了解，ggplot2|详解八大基本绘图要素

ggplot(res_2_test, aes(HR, Variable))+  ##定义X轴和Y轴，以类型分类  geom_point(size=2.5)+  #点的大小  geom_errorbarh(aes(xmax =upper_95, xmin = lower_95), height = 0.4)+ ##95%置信区间，误差线  geom_vline(aes(xintercept = 1))+ #以1为分界线  xlab('HR(95%CI)') + ylab(' ')+  #定义横纵坐标  theme_bw(base_size = 12)

假如部分基因有特定信息，比如是重点的免疫基因，有的是凋亡基因，也可以通过颜色和符号标识出来 （此为示例，并非真实免疫和凋亡基因）。

res_2_test$type <- c(rep("immune",5),rep("",10) ,rep("apoptosis",7) ,rep("",8) )#排序 ，标识type信息ggplot(res_2_test, aes(HR, reorder(Variable,HR),col=type,shape=type))+    geom_point(size=2.5)+    geom_errorbarh(aes(xmax =upper_95, xmin = lower_95), height = 0.4)+   geom_vline(aes(xintercept = 1))+   xlab('HR(95%CI)') + ylab(' ')+    theme_bw(base_size = 12)+   scale_color_manual(values = c("gray", "steelblue", "red")) #设置颜色

四 单因素预后显著基因

根据二中得到的所有基因的单因素生存分析结果，可以根据阈值（p < 0.05）筛选 预后显著的基因集，

KM_sig <- res_2$res %>%   filter(p.value <= 0.05)

提取显著的基因进行后续的分析，一般是lasso 或者其他机器学习的方法进行筛选，然后构建多因素COX模型（风险模型），最终获取每个样本的risk score ，这样再使用更多的数据集验证其稳定性以及其他方面的分析。

◆ ◆ ◆  ◆ ◆

更多精心内容详见：精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）