R语言：多因素cox

> library(survival)

> library(survminer)
> rt=read.table("input.txt",header=T,sep="\t",check.names=F,row.names=1)
> head(rt)

 futime fustat   SCN11A    CCM2L    CRHBP P2RX5-TAX1BP3    RNF166 RP11-452H21.4 MIR600HG   ZNF589 C9orf139   CARMIL2   BCL11A
TCGA-F2-6880 0.79178082      0 3.953795 7.548042 3.458474      4.500245  9.529264      1.929494 7.687410 9.097719 3.744485  5.742071 8.924075
TCGA-RL-AAAS 0.02465753      0 4.143991 8.702790 5.379176      5.227610 10.235375      1.632377 7.493609 9.036441 5.450246  8.710139 9.358276
TCGA-YY-A8LH 5.52328767      0 3.467007 6.520029 1.892517      4.049564  9.067794      1.895853 6.670473 9.020213 5.087604  6.708852 7.414977
TCGA-2J-AABO 1.20547945      0 3.616952 7.048115 3.424015      5.162443  8.898039      2.120082 6.885002 8.541811 3.952987  7.552598 6.146567
TCGA-FB-AAQ0 1.29589041      1 3.020632 5.271652 1.776686      2.530078  9.238666      2.193365 6.858516 9.650095 4.924908  5.956430 5.437092
TCGA-XN-A8T5 1.97260274      0 6.163678 8.479497 4.831885      5.917886 10.613918      3.978881 8.982204 9.404336 5.834611 10.138566 9.425766
                 BTNL9  TMSB4XP1 RP11-147L13.15     LIPE    B3GNT3    CYFIP2     TOB1   TSPAN1     STYK1    FAM53B     IFFO1 HIST1H2BD
TCGA-F2-6880 10.492229 2.5968970       6.576671 9.642247  9.359292 10.960196 12.32776  8.73397  3.111328 10.221345  9.482621  8.555508
TCGA-RL-AAAS  7.542891 3.0526159       7.306673 8.598290 11.331874 11.324258 11.24907 12.21470  7.558937 11.246579 10.281999  7.574798
TCGA-YY-A8LH  9.794470 4.0576097       5.431500 9.428263 13.300859 10.011144 13.54537 14.72629 10.137633 10.705148  8.410821  9.299678
TCGA-2J-AABO  8.373017 0.2320884       6.630489 6.942621 12.079437  9.940927 12.12967 13.57408  8.572091 10.574610  9.844777 11.569551
TCGA-FB-AAQ0  8.717819 1.1665470       4.926700 8.639384 12.675229  8.487557 12.71227 14.19719 10.053513  9.499449  9.181525  9.402634
TCGA-XN-A8T5  8.521171 2.4907681       8.000762 8.493964 11.402173 12.048411 11.71136 13.03989  8.271611 11.477607 10.192197  9.737990

> rt[,3:ncol(rt)]=log2(rt[,3:ncol(rt)]+1)
> rt[,"futime"]=rt[,"futime"]/365

#cox模型建立

> cox <- coxph(Surv(futime, fustat) ~ ., data = rt)
> cox=step(cox,direction = "both")
> riskScore=predict(cox,type="risk",newdata=rt)
> risk=as.vector(ifelse(riskScore>median(riskScore),"high","low"))
>write.table(cbind(id=rownames(cbind(rt[,1:2],riskScore,risk)),cbind(rt[,1:2],riskScore,risk)),file="risk.txt",sep="\t",quote=F,row.names=F)
> cox=read.table("risk.txt",header=T,sep="\t",check.names=F,row.names=1)
head(cox)

 futime fustat riskScore risk
TCGA-F2-6880 0.79178082      0 0.1407025  low
TCGA-RL-AAAS 0.02465753      0 0.7242240  low
TCGA-YY-A8LH 5.52328767      0 1.0683394  low
TCGA-2J-AABO 1.20547945      0 0.8469659  low
TCGA-FB-AAQ0 1.29589041      1 1.7719856 high
TCGA-XN-A8T5 1.97260274      0 0.2416835  low

> fit <- survfit(Surv(futime, fustat) ~ risk, data = cox)
> ggsurvplot(fit, data = cox, risk.table = TRUE)

交流学习