如何确定连续变量的最佳切割点
这里有很详细介绍https://www.rdocumentation.org/packages/survminer/versions/0.4.6/topics/surv_cutpoint
setwd("C:\\Users\\1\\Desktop\\ska3\\12.survival")
gene="SKA3"
library(survival)
library("survminer")
rt=read.table("surviva.txt",header=T,sep="\t",check.names=F)
rt$futime=rt$futime/365
res.cut <- surv_cutpoint(rt, time = "futime", event = "fustat",
variables = c("SKA3"),minprop = 0.25)
summary(res.cut)
plot(res.cut, "SKA3", palette = "npg")
res.cat <- surv_categorize(res.cut)
head(res.cat)
library("survival")
fit <- survfit(Surv(futime, fustat) ~SKA3, data = res.cat)
ggsurvplot(fit, data = res.cat, risk.table = TRUE,,pval = TRUE, conf.int = TRUE)
哭 p值没意义o(╥﹏╥)o
精致图文版
setwd("/Users/zhongyue/Desktop/生物信息学/MAL2/12.survival") #¹¤×÷Ŀ¼£¨ÐèÐ޸ģ©
gene="MAL2"
library(survival)
library("survminer")
rt=read.table("survival.txt",header=T,sep="\t",check.names=F)
rt$futime=rt$futime/365
res.cut <- surv_cutpoint(rt, time = "futime", event = "fustat",
variables = c("MAL2"),minprop = 0.25)
summary(res.cut)
plot(res.cut, "MAL2", palette = "npg")
res.cat <- surv_categorize(res.cut)
head(res.cat)
library("survival")
pdf(file="survival.pdf",
width=6,
height=6)
fit <- survfit(Surv(futime, fustat) ~MAL2, data = res.cat)
ggsurv <- ggsurvplot(fit,pval = TRUE, conf.int = TRUE,ncensor.plot = TRUE,ncensor.plot.height = 0.35,
risk.table = TRUE, risk.table.col = "strata",risk.table.y.text.col = TRUE,xlab = "Time in years",legend.labs=c("MAL2-High","MAL2-Low"),ggtheme = theme_bw())#theme_light()两个主题，risk.table.y.text = FALSE，y轴去掉名字
ggsurv$plot <- ggsurv$plot + labs(title = "Over survival group by MAL2 in old")+theme(plot.title =element_text(size=12,hjust=0.5))#标题居中
ggsurv$table<- ggsurv$table+theme(plot.title =element_text(size=12,hjust=0.5))
ggsurv$ncensor.plot <- ggsurv$ncensor.plot+theme(plot.title =element_text(size=12,hjust=0.5))
print(ggsurv)
dev.off()
中位值划分：
#install.packages("survival")
#install.packages("survminer")
library(survival)
library(survminer)
setwd("C:\\Users\\1\\Desktop\\ska3\\12.survival")
rt=read.table("surviva.txt",header=T,sep="\t",check.names=F,row.names=1)
gene=colnames(rt)[3]
a=ifelse(rt[,gene]<=median(rt[,gene]),"low","high")
diff=survdiff(Surv(futime, fustat) ~a,data = rt)
pValue=1-pchisq(diff$chisq,df=1)
fit=survfit(Surv(futime, fustat) ~ a, data = rt)
if(pValue<0.001){
pValue="<0.001"
}else{
pValue=paste0("=",round(pValue,3))
}
surPlot=ggsurvplot(fit,
data=rt,
conf.int=TRUE,
pval=paste0("p",pValue),
pval.size=6,
risk.table=T,
legend.labs=c("high","low"),
legend.title=paste0(gene," level"),
xlab="Time(years)",
break.time.by = 1,
risk.table.title="",
palette=c("red", "blue"),
risk.table.height=.25)
pdf(file=paste(gene,".survival.pdf",sep=""), width = 6.5, height = 5.5,onefile = FALSE)
print(surPlot)
dev.off()
summary(fit)
