2024-06-16

小师妹带你处理批量生存分析

小师妹生信果 2023-09-02 19:00:34

转自公众号：R语言学徒
http://mp.weixin.qq.com/s?__biz=Mzg5MDk3Mzg4OA==&mid=2247485303&idx=1&sn=e54a06d5236144bfd15e5eac724cb37b

点击蓝字关注我们

生存分析是一种用于研究患者生存时间和相关因素的统计方法。在肿瘤研究中，生存分析常用于评估患者的生存率、预后和治疗效果。生存率估计：通过生存曲线，可以估计患者在不同时间点上存活的概率。生存曲线通常Kaplan-Meier方法绘制，根据患者在观察期间的事件（例如死亡）和失访情况来计算。预后评估：生存分析可以评估不同因素对患者预后的影响。通过将患者分为不同的亚组，可以比较不同亚组之间的生存差异。

常见的方法包括Log-rank检验和Cox比例风险模型，用于确定与生存相关的预测因素。预测生存：基于已知的临床和分子特征，生存分析可以构建生存预测模型。这些模型可以根据患者的特征预测其生存概率，并用于指导治疗决策。在进行生存分析时，通常需要以下数据：患者的生存时间：表示从诊断或治疗开始至死亡或最后随访的时间长度。事件数据：表示患者是否发生了预定义的事件，例如死亡。如果事件没有发生，或者患者在最后随访时仍然存活，就被视为失访。预测因素：包括临床特征（如年龄、性别、病理分级等）和分子特征（如基因表达、蛋白质标记物等）。

生存分析的核心是生存函数和生存曲线。生存函数描述了在给定时间点上个体存活下来的概率，生存曲线则是将生存函数可视化后得到的曲线图。通过生存曲线，我们可以观察不同组别或条件下个体的生存情况是否有显著差异。

#首先，加载所需的包library(survival)library(survminer)library(data.table)library(tidyverse)#使用fread()函数从文件"trainRiskOut.txt"中读取数据，并将其转换为数据框对象。接着，将第一列作为行名，并移除第一列。rt=fread(paste0("./trainRiskOut.txt"),header=T,check.names=F)

rt=as.data.frame(rt)rownames(rt)=rt[,1]rt=rt[,-1]

rt$futime=rt$futime/365#将存活时间（futime）从天数转换为以年为单位，以便进行生存分析。

#对于数据框中的第3到第6列（基因表达数据），进行以下操作：将存活时间（futime）转换为以年为单位。对于每个基因，根据表达值的中位数，将患者分为"high"（高风险）和"low"（低风险）两组。使用survdiff()函数计算两组之间的生存差异，并计算相应的p值。使用survfit()函数拟合生存曲线。使用ggsurvplot()函数生成生存曲线图，包括风险组别的标注、p值的显示和风表。最后，将每个基因的生存曲线图保存为PDF文件。#循环遍历基因列进行生存分析和绘图：for (i in 3:6){gene=colnames(rt[i])group=ifelse(rt[,i]>median(rt[,i]),"high","low")#gene = colnames(rt[i])：获取当前基因的名称。#group = ifelse(rt[, i] > median(rt[, i]), "high", "low")：根据当前基因的表达值，将患者分为"high"（高风险）和"low"（低风险）两组。diff=survdiff(Surv(futime, fustat) ~group,data = rt)#计算两组之间的生存差异，返回相关的统计结果。pValue=1-pchisq(diff$chisq,df=1)if(pValue<0.001){pValue="p<0.001"}else{pValue=paste0("p=",sprintf("%.03f",pValue))}fit <- survfit(Surv(futime, fustat) ~ group, data = rt)#拟合生存曲线模型。surPlot=ggsurvplot(fit,data=rt,title=paste0("gene:",gene=colnames(rt[i]))   ,conf.int=F,pval=pValue,pval.size=6,legend.title="Risk",legend.labs=c("High risk", "Low risk"),xlab="Time(years)",ylab="Overall survival",break.time.by = 3,palette=c("red", "blue"),risk.table=TRUE,risk.table.title="",risk.table.col = "strata",risk.table.height=.25)#使用ggsurvplot()函数生成生存曲线图，其中包括生存曲线、p值、风险表等信息。pdf(file=paste0("Surv_",gene,".pdf"), width = 6.5, height =5.5, onefile = FALSE)print(surPlot)dev.off()  #关闭PDF设备}}