本文首发于“生信补给站”公众号 https://mp.weixin.qq.com/s/3kLClCRiBOLbDbawydVjUg
生存分析作为转录组文章中的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包非常值的狠狠学
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)
注：univ_results部分中的内容可以自定义，包括添加更多的信息，调整小数位数等；2，使用 ezcox 一行输出
ezcox是ShixiangWang大佬 开发的R包，一行代码输出所需结果 ，参考使用 ezcox 进行批量 Cox 模型处理 - 知乎 (zhihu.com)
对比两种方式，结果是一致的。
三 绘制森林图
对于单因素的结果，经常出现的可视化方式就是绘制森林图 。可以使用经典的forestplot-R包绘制（封装），或者使用ggplot2绘制（自由设置）。1 ，forestplot包绘制
forestplot绘制的关键就在于构建tabletext信息。因为基因太多，随机选择30个基因进行绘制
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_exp covariates <- 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))
更多参数设置可以使用 ??forestplot 查看 或者 R-forestplot包| HR结果绘制森林图
2， ggplot2 方式绘制自由度较高，需要对ggplot2有基本的了解，ggplot2|详解八大基本绘图要素
set.seed(1234) res_2_test <- res_2[sample(nrow(res_2), 30), ] #构建tabletext sample <- 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") )
假如部分基因有特定信息，比如是重点的免疫基因，有的是凋亡基因，也可以通过颜色和符号标识出来 （此为示例，并非真实免疫和凋亡基因）。
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)
四 单因素预后显著基因
根据二中得到的所有基因的单因素生存分析结果，可以根据阈值（p < 0.05）筛选 预后显著的基因集，
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")) #设置颜色
提取显著的基因进行后续的分析，一般是lasso 或者其他机器学习的方法进行筛选，然后构建多因素COX模型（风险模型），最终获取每个样本的risk score ，这样再使用更多的数据集验证其稳定性以及其他方面的分析。
