rm(list=ls())
gc() #清列表清内存
#单因素和多因素COX回归
setwd("C:/AW/COX")
library(survival)
data=read.csv("UCSC.csv",header=T) #载入需要做COX的数据，包括存活状态，时间，和对存活有影响的变量
unimodel=coxph(Surv(OSTIME, OS) ~ AGE, data) #对AGE分组做单因素COX
print(unimodel) # 单因素结果，exp(coef)即HR, 保存exp(coef)和p值
exp(confint(unimodel)) #得到HR的95%置信区间，保存
summary(coxmodel) #详细版本的单因素结果
coxmodel=coxph(Surv(OSTIME, OS) ~ AGE + T + STAGE + M + N , data) #多因素分析
print(coxmodel) #多因素结果
exp(confint(coxmodel)) #多因素置信区间
#做森林图
cox=read.csv("COX.csv",header=T) #森林图需要HR值，和HR的95%置信区间上下限
cox$group = ifelse(cox$lower >1|cox$upper <1,"red3", "black") #对有意义的变量标红
library(ggplot2)
ggplot(data=cox)+
aes(x=hr,y=reorder(X,hr))+ #根据HR值对纵坐标的变量排序
geom_errorbarh(aes(xmax=upper,xmin=lower),color=cox$group,height=0,size=0.8)+ #画出HR值上下限
geom_point(size=1.5,shape=15,color=cox$group)+ #画出HR值
geom_vline(xintercept=1,linetype="dashed",size=0.6,color="darkblue")+ #在1处标蓝色虚线
coord_trans(x="log2")+ #对X轴取log
scale_x_continuous(limits=c(0.2,10),breaks=c(0.5,1,2,3,5,10))+ #规定X轴上下限和需要显示的坐标
labs(x="Hazard ratio",y="")+
theme_bw()+
theme(panel.grid.minor = element_blank()) #去除多余杂线
