2万个基因cox回归分析

Cox比例风险模型（考克斯，1972年）是常用的统计在医学研究调查的患者和一个或多个预测变量的存活时间之间的关联回归模型。
基础知识可见：
R-回归分析

一、R语言中基础公式

fit <- coxph(Surv(OS.time, OS) ~ X, data)

**注释：**coxph()为回归拟合函数；Surv(OS.time,OS)为结果变量，X为解释变量；data为数据集

fit2 <- coxph( Surv(stop, event) ~ size, data = bladder )
fit2
summary(fit2)

fit2.png
summary(fit2).png

fit结果展示

Cox回归结果可以解释为：
统计意义: 标记为“ z”的列给出了Wald统计值。它对应于每个回归系数与其标准误差的比率（z = coef / se（coef））。Wald统计量会评估beta给定变量的系数在统计上显着不同于0。
回归系数: Cox模型结果中要注意的第二个特征是回归系数（coef）的符号。正号表示该变量值较高的受试者的危险（死亡风险）较高，因此预后较差。
危险率: 指数系数（exp（coef）= exp（-0.53）= 0.59），也称为危险比，给出了协变量的影响大小。
危险比的置信区间: 摘要输出还给出了危险比（exp（coef））的上限和下限95％置信区间，下限95％和上限95％。
该模型的全局统计意义: 最后，输出为模型的总体重要性提供了三个备选检验的p值：似然比检验，Wald检验和得分对数秩统计。这三种方法渐近等效。对于足够大的N，它们将给出相似的结果。对于较小的N，它们可能会有所不同。对于小样本量，似然比测试具有更好的行为，因此通常是首选。

二、实战

1、数据整理

rm(list = ls())
library(tidyverse)
library(xlsx)

cli <- read.xlsx2("../2.clinical/cli.xlsx",sheetIndex = 1)
load(file = "../3.DEGs/tpm_exp.Rdata")

clin <- cli[cli$Sample.ID%in%group$Sample.ID,c(1,13,14)]
colnames(clin) <- c("sample.ID","RFS.time","RFS")
clin$RFS <- ifelse(clin$RFS == "1:Recurred/Progressed",1,0)
clin$RFS.time <- as.numeric(clin$RFS.time)
str(clin)

exp <- tpm_exp[,str_sub(colnames(tpm_exp),1,15)%in%clin$sample.ID]
colnames(exp)<- str_sub(colnames(exp),1,15)
clin <- clin[match(colnames(exp),clin$sample.ID),]

save(clin,exp,file = "pre_modle.Rdata")

clin.png
exp.png

数据整理之后，下面的工作程式化输入即可

2、2万个基因的cox回归分析（循环）

rm(list = ls())
library(survival)
library(survminer)
library(tidyverse)

#--------------------------------------------------------
load("pre_modle.Rdata")
exp <- t(scale(t(exp)))
precox <- merge(clin, t(exp), by.x = 1,by.y = 0)
uni_result=data.frame()
for(i in colnames(precox[,4:ncol(precox)])){
  cox <- coxph(Surv(RFS.time, RFS) ~ get(i), data = precox)
  coxSummary = summary(cox)
  uni_result=rbind(uni_result,
                   cbind(id=i,
                         HR=coxSummary$conf.int[,"exp(coef)"],
                         HR.95L=coxSummary$conf.int[,"lower .95"],
                         HR.95H=coxSummary$conf.int[,"upper .95"],
                         pvalue=coxSummary$coefficients[,"Pr(>|z|)"]))
}
uni_result[,2:5] <- apply(uni_result[,2:5],2,as.numeric)
table(uni_result$pvalue<0.05) 
write.table(uni_result ,'TCGA-Unicox-result.xls',
            row.names = F,sep = '\t',quote = F)

数据展示：

unicox结果展示