R语言多基因最佳分组生存曲线

生信喵

背景

批量画生存分析Best separation的图，所有基因保存到一个文件，或每个基因1个文件。

只保留p value<0.05的基因

在图上标注HR、p value、95% CI

设置分组cutoff，例如组内sample数量不小于30%

找best separation用的是survminer的函数，之前我们做过画出最佳分组的生存曲线

使用场景

在展示基因表达水平（连续变量）对生存期的影响时找到最佳分组，多个基因，一次出图

后面可以接 双因子生存分析，对比两个基因的高低组合对预后的影响。

输入数据

第一列是生存时间、第二列是追踪情况、第三列以后是基因表达水平

svdata <- read.csv("easy_input.csv",header = T,row.names = 1)
dim(svdata) #一共18个基因

svdata前几行如下：

开始画图

##加载必需的包
library(survival)
library(survminer)

##对数据集的基因进行bestSeparation统计
res.cut <- surv_cutpoint(svdata, time = "futime", 
                         event = "fustat", 
                         variables = names(svdata)[3:ncol(svdata)], 
                         minprop = 0.3) #默认组内sample不能低于30%

##按照bestSeparation分高低表达
res.cat <- surv_categorize(res.cut)

##统计作图
my.surv <- Surv(res.cat$futime, res.cat$fustat)
pl<-list()
for (i in colnames(res.cat)[3:ncol(svdata)]) {
  group <- res.cat[,i] 
  survival_dat <- data.frame(group = group)
  fit <- survfit(my.surv ~ group)
  
  ##计算HR以及95%CI
  ##修改分组参照
  group <- factor(group, levels = c("low", "high"))
  data.survdiff <- survdiff(my.surv ~ group)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
  up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
  low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
  
  #只画出p value<=0.05的基因，如果不想筛选，就删掉下面这行
  if (p.val>0.05) next
  
  HR <- paste("Hazard Ratio = ", round(HR,2), sep = "")
  CI <- paste("95% CI: ", paste(round(low95,2), round(up95,2), sep = " - "), sep = "")
  
  #按照基因表达量从低到高排序，便于取出分界表达量
  svsort <- svdata[order(svdata[,i]),]
  
  pl[[i]]<-ggsurvplot(fit, data = survival_dat ,
           #ggtheme = theme_bw(), #想要网格就运行这行
           conf.int = F, #不画置信区间，想画置信区间就把F改成T
           #conf.int.style = "step",#置信区间的类型，还可改为ribbon
           censor = F, #不显示观察值所在的位置
           palette = c("#D95F02","#1B9E77"), #线的颜色对应高、低
         
           legend.title = i,#基因名写在图例题目的位置
           font.legend = 11,#图例的字体大小
           #font.title = 12,font.x = 10,font.y = 10,#设置其他字体大小

           #在图例上标出高低分界点的表达量，和组内sample数量
           legend.labs=c(paste0(">",round(svsort[fit$n[2],i],2),"(",fit$n[1],")"),
                           paste0("<",round(svsort[fit$n[2],i],2),"(",fit$n[2],")")),
         
           #在左下角标出pvalue、HR、95% CI
           #太小的p value标为p < 0.001
           pval = paste(pval = ifelse(p.val < 0.001, "p < 0.001", 
                                             paste("p = ",round(p.val,3), sep = "")),
                        HR, CI, sep = "\n"))
  
  #如果想要一个图保存为一个pdf文件，就把下面这行前面的“#”删掉
  ggsave(paste0(i,".pdf"),width = 4,height = 4)
}

length(pl)

只保留p<=0.05的10个基因，筛掉了p>0.05的8个基因

批量出图

用survminer包自带的函数组图

res <- arrange_ggsurvplots(pl, 
                           print = T,
                           ncol = 3, nrow = 4)#每页纸画几列几行

#保存到pdf文件
ggsave("bestSurvPlot.pdf",res,width = 12,height = 16)

生信喵

文中所用数据可以关注公众号“生信喵实验柴”

发送关键词“20250415”获取