|
前面我們講過一個R函數(shù)搞定風(fēng)險評估散點圖，熱圖，其中LASSO模型的輸入就是單因素cox分析得到的顯著與生存相關(guān)的基因。今天我們就來探討一下如何使用R來做單因素和多因素cox回歸分析。 我們用R的survival包自帶的一套肺癌的數(shù)據(jù)來舉例 #安裝下面兩個R包install.packages(c("survival", "survminer"))#加載這兩個R包library("survival")library("survminer")#加載肺癌這套數(shù)據(jù)data("lung")#顯示前6行head(lung)這里每一行是一個樣本，從第三列開始每一列是一個特征 1.單因素cox回歸分析 對單個特征進(jìn)行cox回歸分析，看它是否與樣本的生存顯著相關(guān) #單因素cox回歸分析，這里看性別sex這個特征res.cox <- coxph(Surv(time, status) ~ sex, data = lung)res.cox可以看到這里算出來的p值是0.00149，是顯著的 我們在來看一下summary summary(res.cox)這里的exp(coef)就是HR（hazard ratio，風(fēng)險率），lower .95和upper .95為95%的置信區(qū)間 2.批量單因素cox回歸分析 一般我們的關(guān)注的特征都比較多，用上面的代碼一個一個來做單因素cox回歸分析效率太低了，下面我們來看看如何批量做單因素cox回歸分析。 #假設(shè)我們要對如下5個特征做單因素cox回歸分析covariates <- c("age", "sex", "ph.karno", "ph.ecog", "wt.loss")#分別對每一個變量，構(gòu)建生存分析的公式univ_formulas <- sapply(covariates, function(x) as.formula(paste('Surv(time, status)~', x)))#循環(huán)對每一個特征做cox回歸分析univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})#提取HR，95%置信區(qū)間和p值univ_results <- lapply(univ_models, function(x){ x <- summary(x) #獲取p值 p.value<-signif(x$wald["pvalue"], digits=2) #獲取HR HR <-signif(x$coef[2], digits=2); #獲取95%置信區(qū)間 HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2) HR.confint.upper <- signif(x$conf.int[,"upper .95"],2) HR <- paste0(HR, " (", HR.confint.lower, "-", HR.confint.upper, ")") res<-c(p.value,HR) names(res)<-c("p.value","HR (95% CI for HR)") return(res) })#轉(zhuǎn)換成數(shù)據(jù)框，并轉(zhuǎn)置res <- t(as.data.frame(univ_results, check.names = FALSE))as.data.frame(res)write.table(file="univariate_cox_result.txt",as.data.frame(res),quote=F,sep="\t")得到的結(jié)果如下，你會發(fā)現(xiàn)對于sex這個特征來說，結(jié)果跟前面單獨做得到的結(jié)果是一樣的。 3.多因素cox回歸分析 前面是單獨看每一個特征是否跟生存相關(guān)，而多因素cox回歸是同時檢測多個特征是否與生存相關(guān)。一般先通過單因素cox回歸分析找出與生存顯著相關(guān)的特征，然后基于這些特征再去做多因素cox回歸分析，或者做LASSO分析。 res.cox <- coxph(Surv(time, status) ~ age + sex + ph.karno + ph.ecog, data = lung)x <- summary(res.cox)pvalue=signif(as.matrix(x$coefficients)[,5],2)HR=signif(as.matrix(x$coefficients)[,2],2)low=signif(x$conf.int[,3],2)high=signif(x$conf.int[,4],2)multi_res=data.frame(p.value=pvalue, HR=paste(HR," (",low,"-",high,")",sep=""), stringsAsFactors = F)multi_reswrite.table(file="multivariate_cox_result.txt",multi_res,quote=F,sep="\t")得到的結(jié)果如下 參考資料： |
|
|
