在生信分析的过程中，单因素、多因素COX回归分析常用来筛选与患者预后有关的因素，相信很多小伙伴对此并不陌生。那么，如何在R中进行单因素、多因素COX回归分析？COX回归分析与KM生存分析的区别在哪？做了单因素COX回归分析之后为什么还要做多因素COX回归分析呢？用于筛选变量的p值要设置为多少呢？这篇文章都会给你答案，快和小碗看下去吧！
话不多说，我们直接上代码：
首先，我们要养成良好的习惯，先设置一下环境，不仅能让代码看起来更加清晰，而且能为后续运行省去很多不必要的麻烦：
#清除工作空间中的所有对象
rm(list=ls())
#清理内存
gc()
#设置工作路径
setwd("..")
#加载包
library(survival)
#加载示例数据
load("rdata/示例数据.rdata")
我们用str()函数看一下示例数据df的格式：
可以看到，示例数据df是一个数据框，有50个观测（行），102个变量（列）。
使用view()函数查看df：
可以看到，df的50个观察就是50个样本，从sample1到sample50；102个变量中的前两个分别是每个样本的生存状态和生存时间，从第三个到最后一个就是我们想要进行单因素、多因素COX回归分析的100个基因了。
接下来，开始正式进行分析：
###单因素COX回归分析----
#设置p值的阈值
pfilter <- 0.05
#新建空白数据框
uniresult <- data.frame()
#使用for循环对输入数据中的100个基因依次进行单因素COX分析
#单因素COX回归分析中p值＜0.05的基因，其分析结果输入到之前新建的空白数据框uniresult中
for(i in colnames(df[,3:ncol(df)])){
unicox <- coxph(Surv(time = os_time, event = os_status) ~ df[,i], data = df)
unisum<- summary(unicox)
pvalue <- round(unisum$coefficients[,5],3)
if(pvalue<pfilter){
uniresult <- rbind(uniresult,
cbind(gene=i,
HR=unisum$coefficients[,2],
L95CI=unisum$conf.int[,3],
H95CI=unisum$conf.int[,4],
pvalue=unisum$coefficients[,5]
))
}
}
#保存单因素COX回归分析结果
write.csv(uniresult,file = "result/单因素COX分析结果.csv",row.names = F)
在上述代码中，我们使用到了for循环对示例数据中的100个基因依次进行单因素COX分析，并将其中p<0.05的基因的基因名、HR值、HR的95%置信区间的低值与高值、p值提取出来保存到uniresult这个对象当中，我们来看一下uniresult：
也是一个数据框，有12行，说明通过单因素COX回归分析，我们从100个基因中筛选出了12个与病人生存有关的基因。
也可以使用view()函数查看一下uniresult，更加直观：
接着，我们就对这12个基因再进行多因素COX回归分析。我们要提取这12个基因在各个样本中的表达情况，并与各样本的生存状态和生存时间组成数据框，作为多因素COX回归分析的输入数据unigene。
#提取各样本的生存状态、生存时间、
#以及单因素COX回归分析中p值＜0.05的基因在各样本中的表达情况，
#作为多因素COX回归分析的输入数据
unigene <- subset(df,select = c("os_status","os_time",uniresult$gene))
我们来看一下unigene：
可以看到unigene也是个数据框，有50个样本，14列（1、2列分别是生存状态、生存时间，3-14列是单因素COX筛选出的12个基因）。
使用view()函数查看一下：
和单因素COX回归分析的输入数据——示例数据df结构一致。
然后，就可以开始进行多因素COX回归分析了；
###多因素COX回归分析----
multicox <- coxph(Surv(time = os_time,event = os_status) ~ ., data = unigene)
multisum <- summary(multicox)
#提取所有基因的多因素COX回归分析结果至multiresult对象中
gene <- colnames(unigene)[3:ncol(unigene)]
HR <- multisum$coefficients[,2]
L95CI <- multisum$conf.int[,3]
H95CI <- multisum$conf.int[,4]
pvalue <- multisum$coefficients[,5]
multiresult <- data.frame(gene=gene,
HR=HR,
L95CI=L95CI,
H95CI=H95CI,
pvalue=pvalue)
multiresult <- multiresult[multiresult$pvalue<pfilter,]
#保存多因素COX回归分析结果
write.csv(multiresult,file = "result/多因素COX分析结果.csv",row.names = F)
多因素COX回归分析是同时分析多个元素（这里是12个基因）与样本的生存的关系，所以代码相对于单因素COX回归分析会更加简单，不用for循环，一行代码就搞定了，结果保存在multiresult这个对象当中，我们来看一下：
可以看到，multiresult是一个12行、5列的数据框。注意，multiresult中是全部进行分析的12个基因的基因名、HR值、HR的95%置信区间的低值与高值、p值，而不仅仅是p<0.05的基因。
用view()再看一下：
可以发现，很不凑巧，这12个基因经过多因素COX回归分析后p值均大于0.05。
以上呢，就是进行单因素、多因素COX回归分析的全部代码了。对于单因素、多因素COX回归分析的结果，我们可以绘制森林图来进行可视化，由于篇幅限制，森林图的代码我们后续再分享，请各位观众老爷点个关注，才不会迷路哦~
讲完代码，我们顺带来解释一下几个常见的问题：
①有没有小伙伴和小碗有过同样的疑问——KM生存分析、COX回归分析都是分析与病人预后有关的因素的方法，那它们两者之间有什么区别呢？
我们来回忆一下KM生存分析的一般流程是根据某种因素比如性别、肿瘤组织学分级、某个基因表达量的高低将样本分成几组，然后比较各个组别之间的生存的差别，计算p值。那就要求我们要分析的变量无论本身是分类变量也好，还是转化成分类变量也好，总之得是分类变量，这样才能对样本进行分组嘛，而且无论你是根据性别分为两组也好，还是根据组织学分级分为3-4组也好，KM生存分析一次只能分析一个变量对生存的影响，也就是说，KM生存分析本质也是一个单变量生存分析。
这样，把KM简单拉一遍之后，COX回归分析的优势就很明显了：首先它能同时分析多个因素与生存的关系，也就是多因素COX回归分析；其次，它能分析连续型变量对生存的影响，我们示例数据中分析的就是基因在样本中的表达情况（连续型变量）对样本生存的影响。
②KM生存分析是单变量分析，单因素COX分析也是单变量分析，那它们的分析结果一样吗？
小碗使用如下代码将示例数据中的100个基因一次进行了KM生存分析，也以p<0.05作为阈值筛选出了与病人生存有关的基因：
这部分代码不做讲解，感兴趣的同学可以自行学习。
我们来看一下KM生存分析、单因素COX回归分析筛选出的基因分别是km_gene、cox_gene：
可以看到KM生存分析共筛选出8个，而单因素COX回归分析筛选出12个，且它们并不完全相同。这是因为，虽然都是分析某个因素，但他们内部的算法是不同的，所以结果有出入也很正常。我们根据文献选择更适合自己研究的分析即可。
写到这里，留个问题给大家：基因表达量不是连续变量吗？怎么才能用KM生存分析呢？答案前面讲过哦。
③为什么大多数生信文章中都是先进行单因素COX回归分析再进行多因素COX回归分析？
我们进行单因素COX回归分析的时候一次只纳入一个变量，就不能排除有这么一种情况存在：变量a影响生存其实是通过影响变量b而导致的，也就是说变量a其实对生存没有直接的影响。那么，就需要将单因素COX 分析筛选出来的显著变量继续纳入多因素COX比例风险模型，从而排除掉类似变量a的这种混杂因素。单因素COX回归分析像是对变量的一个初步筛选，多因素COX回归分析筛选出来的才是独立预后因素。
好了，今天的分享到这里就结束啦，如果对你有帮助的话就给小碗一个免费的赞吧，我们下期再见！
