R语言画KM生存曲线如何标中位生存时间 km生存曲线p值

转载

网线小游侠 2024-08-24 22:01:43

需要的文件

表型文件（包括生存状态等），表达值（以某个基因表达量作为分组的话）

如果是手动下载的文件需要读取成可用的格式，用fread函数

phe=fread("phe.txt",sep="\t",header = F, data.table = F)

exp=fread("GSE12417-GPL96_series_matrix.txt",
          skip=72,data.table = F)

##下载的文件需要先查看一下，phe.txt是手动从matrix.txt文件里复制粘贴创建的（characteristic行）

第一环节：数据预处理:形成一个带有生存状态、生存时间、分组基因表达量三列数据的文件

第一步：从phe文件里提取生存状态和生存时间

需要用到字符分割、合并

s=phe$characteristics_ch1
s1=strsplit(s,split = ";",fixed = T)
s1[[1]]

查看

R语言画KM生存曲线如何标中位生存时间 km生存曲线p值_数据

##把分割后的字符串的第一个元素提取出来，合并成为一个新的向量
os.time = sapply(s1,function(x){x[3]})
status= sapply(s1,function(x){x[4]})
os.time[1:5]

R语言画KM生存曲线如何标中位生存时间 km生存曲线p值_图例_02

os.time1=strsplit(os.time,split = " ",fixed = T)
os.time2= sapply(os.time1,function(x){x[4]})
os.time2=as.numeric(os.time2)

status[1:5]
status1=strsplit(status,split = " ",fixed = T)
status2=sapply(status1,function(x){x[4]})
status2=as.numeric(status2)

s2=data.frame(os.time2,status2)
rownames(s2)=rownames(phe)###行名变成样本名，好与表达矩阵merge

第二步：提取表达矩阵中感兴趣的基因作为分组，合成生存分析所需要的文件，以EZH2为例

exp3.e=exp3["EZH2",]
?t
exp3.e <-data.frame(t(exp3.e))
s3=merge(exp3.e,s2,by.x = 0,by.y = 0)

此时的s3预览

R语言画KM生存曲线如何标中位生存时间 km生存曲线p值_r语言_03

第二环节：分组：划分基因表达为高或低分组

ifelse函数

##############################################################################
###四分位生存曲线
quantile(s3$EZH2)
s3$zz=ifelse(s3$EZH2>9.53,'high', 
                 ifelse( s3$EZH2>9.23 & s3$EZH2<9.53,'h.stable', 
                         ifelse( s3$EZH2>8.88& s3$EZH2<9.23,'l.stable','down') ))
############################################################################
##中位数分组
s2.exp$EZH2 = ifelse(s2.exp$EZH2>median(s2.exp$EZH2),'high','low')
s2.exp$EZH2 = ifelse(s2.exp$EZH2>9.235157,'H','low')
s3$group=ifelse(s3$EZH2>median(s3$EZH2),"H","low")
##s3$EHZ2=blabla意思是新增一列

（s2.exp和s3是一个文件）

第三环节：生存分析、画图

library(survival)
install.packages("survminer")
library(survminer)

#######################s3和s2.exp是一样的，s3是我自己搞的，后者是粘贴的代码###
Surv(os.time2,status2)
?survfit
x1=survfit(Surv(os.time2, status2)~group, data=s3)
###在s3这个文件里，对group这一列进行生存分析
summary(x1)
class(x1)
#Surv()函数创建生存数据对象（主要输入生存时间和状态逻辑值）

#survfit()函数对生存数据对象拟合生存函数，创建KM(Kaplan-Meier)生存曲线

#survdiff()用于不同组的统计检验

summary(fit)$table 查看生存分析的数据

然后画图：

ggsurvplot(x1, data = s3, pval = T,ylab="Survival probability",conf.int=T,
           xlab="Time(Days)",linetype=c("dashed","solid"),palette="lancet",risk.table=T,pval.coord = c(90, 0.9))


ggsurvplot(x1,   
           main = "Survival curve", # 添加标题
           font.main = c(16, "bold", "darkblue"), # 设置标题字体大小、格式和颜色
           font.x = c(14, "bold.italic", "red"), # 设置x轴字体大小、格式和颜色
           font.y = c(14, "bold.italic", "darkred"), # 设置y轴字体大小、格式和颜色
           font.tickslab = c(12, "plain", "darkgreen")) # 设置坐标轴刻度字体大小、格式和颜色


ggsurvplot(x1, 
           surv.median.line = "hv", # 添加中位数生存时间线
           
           # Change legends: title & labels
           legend.title = "EZH2", # 设置图例标题
           legend.labs = c("high", "low"), # 指定图例分组标签
           
           # Add p-value and tervals
           pval = TRUE, # 设置添加P值
           pval.method = FALSE, #设置添加P值计算方法
           conf.int = TRUE, # 设置添加置信区间
           
           # Add risk table
           risk.table = TRUE, # 设置添加风险因子表
           tables.height = 0.2, # 设置风险表的高度
           tables.theme = theme_gray(), # 设置风险表的主题
           
           # Color palettes. Use custom color: c("#E7B800", "#2E9FDF"),
           # or brewer color (e.g.: "Dark2"), or ggsci color (e.g.: "jco")
           palette = c("#34ebcc", "#4f34eb"), # 设置颜色画板
           ggtheme = theme_classic() # Change ggplot2 theme
)

还有几种不同的自定义设置查看当时R脚本，不一一粘贴了

绘制累积风险曲线和总患者生存曲线，还有分面生存曲线（ggsurvplot_facet()函数），

批量输出一批基因作为分组的生存曲线等。