在生信分析的過(guò)程中，單因素、多因素COX回歸分析常用來(lái)篩選與患者預(yù)后有關(guān)的因素，相信很多小伙伴對(duì)此并不陌生。那么，如何在R中進(jìn)行單因素、多因素COX回歸分析？COX回歸分析與KM生存分析的區(qū)別在哪？做了單因素COX回歸分析之后為什么還要做多因素COX回歸分析呢？用于篩選變量的p值要設(shè)置為多少呢？這篇文章都會(huì)給你答案，快和小碗看下去吧！
話不多說(shuō)，我們直接上代碼：
首先，我們要養(yǎng)成良好的習(xí)慣，先設(shè)置一下環(huán)境，不僅能讓代碼看起來(lái)更加清晰，而且能為后續(xù)運(yùn)行省去很多不必要的麻煩：
#清除工作空間中的所有對(duì)象
rm(list=ls()) ?
#清理內(nèi)存
gc() ?
#設(shè)置工作路徑
setwd("..")
#加載包
library(survival)
#加載示例數(shù)據(jù)
load("rdata/示例數(shù)據(jù).rdata") ?
我們用str()函數(shù)看一下示例數(shù)據(jù)df的格式：
可以看到，示例數(shù)據(jù)df是一個(gè)數(shù)據(jù)框，有50個(gè)觀測(cè)（行），102個(gè)變量（列）。
使用view()函數(shù)查看df：
可以看到，df的50個(gè)觀察就是50個(gè)樣本，從sample1到sample50；102個(gè)變量中的前兩個(gè)分別是每個(gè)樣本的生存狀態(tài)和生存時(shí)間，從第三個(gè)到最后一個(gè)就是我們想要進(jìn)行單因素、多因素COX回歸分析的100個(gè)基因了。
接下來(lái)，開始正式進(jìn)行分析：
###單因素COX回歸分析----
#設(shè)置p值的閾值
pfilter <- 0.05 ??
#新建空白數(shù)據(jù)框
uniresult <- data.frame() ?
#使用for循環(huán)對(duì)輸入數(shù)據(jù)中的100個(gè)基因依次進(jìn)行單因素COX分析
#單因素COX回歸分析中p值＜0.05的基因，其分析結(jié)果輸入到之前新建的空白數(shù)據(jù)框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回歸分析結(jié)果
write.csv(uniresult,file = "result/單因素COX分析結(jié)果.csv",row.names = F)
在上述代碼中，我們使用到了for循環(huán)對(duì)示例數(shù)據(jù)中的100個(gè)基因依次進(jìn)行單因素COX分析，并將其中p<0.05的基因的基因名、HR值、HR的95%置信區(qū)間的低值與高值、p值提取出來(lái)保存到uniresult這個(gè)對(duì)象當(dāng)中，我們來(lái)看一下uniresult：
也是一個(gè)數(shù)據(jù)框，有12行，說(shuō)明通過(guò)單因素COX回歸分析，我們從100個(gè)基因中篩選出了12個(gè)與病人生存有關(guān)的基因。
也可以使用view()函數(shù)查看一下uniresult，更加直觀：
接著，我們就對(duì)這12個(gè)基因再進(jìn)行多因素COX回歸分析。我們要提取這12個(gè)基因在各個(gè)樣本中的表達(dá)情況，并與各樣本的生存狀態(tài)和生存時(shí)間組成數(shù)據(jù)框，作為多因素COX回歸分析的輸入數(shù)據(jù)unigene。
#提取各樣本的生存狀態(tài)、生存時(shí)間、
#以及單因素COX回歸分析中p值＜0.05的基因在各樣本中的表達(dá)情況，
#作為多因素COX回歸分析的輸入數(shù)據(jù)
unigene <- subset(df,select = c("os_status","os_time",uniresult$gene))
我們來(lái)看一下unigene：
可以看到unigene也是個(gè)數(shù)據(jù)框，有50個(gè)樣本，14列（1、2列分別是生存狀態(tài)、生存時(shí)間，3-14列是單因素COX篩選出的12個(gè)基因）。
使用view()函數(shù)查看一下：
和單因素COX回歸分析的輸入數(shù)據(jù)——示例數(shù)據(jù)df結(jié)構(gòu)一致。
然后，就可以開始進(jìn)行多因素COX回歸分析了；
###多因素COX回歸分析----
multicox <- coxph(Surv(time = os_time,event = os_status) ~ ., data = unigene)
multisum <- summary(multicox)
#提取所有基因的多因素COX回歸分析結(jié)果至multiresult對(duì)象中
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回歸分析結(jié)果
write.csv(multiresult,file = "result/多因素COX分析結(jié)果.csv",row.names = F)
多因素COX回歸分析是同時(shí)分析多個(gè)元素（這里是12個(gè)基因）與樣本的生存的關(guān)系，所以代碼相對(duì)于單因素COX回歸分析會(huì)更加簡(jiǎn)單，不用for循環(huán)，一行代碼就搞定了，結(jié)果保存在multiresult這個(gè)對(duì)象當(dāng)中，我們來(lái)看一下：
可以看到，multiresult是一個(gè)12行、5列的數(shù)據(jù)框。注意，multiresult中是全部進(jìn)行分析的12個(gè)基因的基因名、HR值、HR的95%置信區(qū)間的低值與高值、p值，而不僅僅是p<0.05的基因。
用view()再看一下：
可以發(fā)現(xiàn)，很不湊巧，這12個(gè)基因經(jīng)過(guò)多因素COX回歸分析后p值均大于0.05。
以上呢，就是進(jìn)行單因素、多因素COX回歸分析的全部代碼了。對(duì)于單因素、多因素COX回歸分析的結(jié)果，我們可以繪制森林圖來(lái)進(jìn)行可視化，由于篇幅限制，森林圖的代碼我們后續(xù)再分享，請(qǐng)各位觀眾老爺點(diǎn)個(gè)關(guān)注，才不會(huì)迷路哦~
講完代碼，我們順帶來(lái)解釋一下幾個(gè)常見(jiàn)的問(wèn)題：
①有沒(méi)有小伙伴和小碗有過(guò)同樣的疑問(wèn)——KM生存分析、COX回歸分析都是分析與病人預(yù)后有關(guān)的因素的方法，那它們兩者之間有什么區(qū)別呢？
我們來(lái)回憶一下KM生存分析的一般流程是根據(jù)某種因素比如性別、腫瘤組織學(xué)分級(jí)、某個(gè)基因表達(dá)量的高低將樣本分成幾組，然后比較各個(gè)組別之間的生存的差別，計(jì)算p值。那就要求我們要分析的變量無(wú)論本身是分類變量也好，還是轉(zhuǎn)化成分類變量也好，總之得是分類變量，這樣才能對(duì)樣本進(jìn)行分組嘛，而且無(wú)論你是根據(jù)性別分為兩組也好，還是根據(jù)組織學(xué)分級(jí)分為3-4組也好，KM生存分析一次只能分析一個(gè)變量對(duì)生存的影響，也就是說(shuō)，KM生存分析本質(zhì)也是一個(gè)單變量生存分析。
這樣，把KM簡(jiǎn)單拉一遍之后，COX回歸分析的優(yōu)勢(shì)就很明顯了：首先它能同時(shí)分析多個(gè)因素與生存的關(guān)系，也就是多因素COX回歸分析；其次，它能分析連續(xù)型變量對(duì)生存的影響，我們示例數(shù)據(jù)中分析的就是基因在樣本中的表達(dá)情況（連續(xù)型變量）對(duì)樣本生存的影響。
②KM生存分析是單變量分析，單因素COX分析也是單變量分析，那它們的分析結(jié)果一樣嗎？
小碗使用如下代碼將示例數(shù)據(jù)中的100個(gè)基因一次進(jìn)行了KM生存分析，也以p<0.05作為閾值篩選出了與病人生存有關(guān)的基因：
這部分代碼不做講解，感興趣的同學(xué)可以自行學(xué)習(xí)。
我們來(lái)看一下KM生存分析、單因素COX回歸分析篩選出的基因分別是km_gene、cox_gene：
可以看到KM生存分析共篩選出8個(gè)，而單因素COX回歸分析篩選出12個(gè)，且它們并不完全相同。這是因?yàn)?，雖然都是分析某個(gè)因素，但他們內(nèi)部的算法是不同的，所以結(jié)果有出入也很正常。我們根據(jù)文獻(xiàn)選擇更適合自己研究的分析即可。
寫到這里，留個(gè)問(wèn)題給大家：基因表達(dá)量不是連續(xù)變量嗎？怎么才能用KM生存分析呢？答案前面講過(guò)哦。
③為什么大多數(shù)生信文章中都是先進(jìn)行單因素COX回歸分析再進(jìn)行多因素COX回歸分析？
我們進(jìn)行單因素COX回歸分析的時(shí)候一次只納入一個(gè)變量，就不能排除有這么一種情況存在：變量a影響生存其實(shí)是通過(guò)影響變量b而導(dǎo)致的，也就是說(shuō)變量a其實(shí)對(duì)生存沒(méi)有直接的影響。那么，就需要將單因素COX 分析篩選出來(lái)的顯著變量繼續(xù)納入多因素COX比例風(fēng)險(xiǎn)模型，從而排除掉類似變量a的這種混雜因素。單因素COX回歸分析像是對(duì)變量的一個(gè)初步篩選，多因素COX回歸分析篩選出來(lái)的才是獨(dú)立預(yù)后因素。
好了，今天的分享到這里就結(jié)束啦，如果對(duì)你有幫助的話就給小碗一個(gè)免費(fèi)的贊吧，我們下期再見(jiàn)！
