2020-01-08单因素COX分析和建立模型

cox_t<-cox_train[,-c(1:9)]

cox_t<-log2(cox_t+0.1)###记得查看结果，是否有NAN，很奇快，没项独立运输没有NAN，但这一句一次行运算，很多NAN。后面风险模型评估就没法建立。

cox_train_univar<-cbind(cox_train[,c(8,6)],cox_t)

cox_train_univar[,2]<-cox_train_univar[,2]/365

colnames(cox_train_univar)[1]<-"status"

colnames(cox_train_univar)[2]<-"life_span"

###单因素cox分析

univar_out<-data.frame(matrix(NA,(ncol(cox_train_univar)-2),5))

rownames(univar_out)<-colnames(cox_train_univar)[-c(1:2)]

colnames(univar_out)<-c("Coeffcient","HR","lower.95","upper.95","P-Value")

for (i in colnames(cox_train_univar)[-c(1:2)]) {

  cox<-coxph(Surv(life_span,status)~cox_train_univar[,i],data = cox_train_univar)

  cox_summary<-summary(cox)

  univar_out[i,1]<-cox_summary$coefficients[,1]

  univar_out[i,2]<-cox_summary$coefficients[,2]

  univar_out[i,3]<-cox_summary$conf.int[,3]

  univar_out[i,4]<-cox_summary$conf.int[,4]

  univar_out[i,5]<-cox_summary$coefficients[,5]

  }

数据格式(行名是所有的肿瘤样本名）

   status    life_span  gene   gene   gene ......

     0             365       ...         ...       .....

   1              235     .....          ....        .....

##通过随机森林筛选关键数个基因

uni001<-univar_out[univar_out$`P-Value`<0.01,]

uni005<-univar_out[univar_out$`P-Value`<0.05,]

bbb=colnames(cox_train_univar)[-c(1,2)]%in%rownames(uni005)###记得删除【-c(1,2)】,否则后面的取值全部错两列

cox_train_multi<-cbind(cox_train_univar[,c(1,2)],cox_train_univar[,bbb])

diseaserf<-rfsrc(Surv(life_span,status)~.,data = cox_train_multi,ntree=500,tree.err = T,importance = T)

outrf<-var.select(object = diseaserf,vdv,method = "vh",nrep = 200)##这一步多跑几次，每次结果不一样，而且30可以换更大值

outrf$topvars

###多因素COX分析 找出是否有统计学意义

sig4genename<-c("ENSG00000275149","ENSG00000245526","ENSG00000230838","ENSG00000237423")#我只找出四个

cox_train_multi_sig<-cbind(cox_train_univar[,c(1,2)],cox_train_univar[,sig5genename])

cox_4<-coxph(Surv(life_span,status)~.,data = cox_train_multi_sig)##建立模型，其实如果有意义的基因不多的话，可以省去随机森林这一步

cox_4<-step(cox_4,direction = "both")

save(cox_4,file = "cox_4.Rdata")

riskscore_train<-predict(cox_4,type = "risk",newdata = cox_train_multi)#通过模型，给每个病例重新计算分值

risktrain<-as.vector(ifelse(riskscore_train>median(riskscore_train),"high","low"))##记着下次测试组时也用median(riskscore_train)

rownames(cox_train_multi_sig)<-cox_train$clinical.submitter_id

risk_traindata<-cbind(id=rownames(cox_train_multi_sig),cox_train_multi_sig[,c(1,2)],riskscore_train,risktrain)