在生信分析的過程中，

單因素、多因素COX回歸

分析常用來篩選與患者預(yù)后有關(guān)的因素，相信很多小伙伴對此并不陌生。 那么，如何在R中進(jìn)行單因素、多因素COX回歸分析？之前分享過繪制KM曲線，諾莫圖展示COX結(jié)果，如何使用survminer包繪制Cox生存分析結(jié)果的

森林圖

。 這篇文章都會(huì)給你答案，小果接下來就給大家進(jìn)行超詳細(xì)的代碼實(shí)操展示。

下面是示例樣本的生存數(shù)據(jù)信息和臨床數(shù)據(jù)文件

#install.packages('survival') library(survival)??????#引用包

############定義森林圖函數(shù)############

#定義了一個(gè)名為bioForest的函數(shù)， #該函數(shù)有三個(gè)參數(shù)：coxFile、forestFile和forestCol，這些參數(shù)可以在函數(shù)被調(diào)用時(shí)傳遞實(shí)際值。 bioForest=function(coxFile=null, forestFile=null, forestCol=null){ #讀取輸入文件

rt <- read.table(coxFile, header=T, sep="\t", check.names=F, row.names=1)

#用read.table函數(shù)從coxFile中讀取數(shù)據(jù)，數(shù)據(jù)以Tab分隔，包含表頭（header=T）， #第一列作為行名（row.names=1），并將數(shù)據(jù)保存到名為rt的數(shù)據(jù)框中。 gene <- rownames(rt) #將rt中的行名（基因名稱）保存到gene變量中。 hr <- sprintf("%.3f",rt$"HR") #將rt數(shù)據(jù)框中"HR"列的值保留三位小數(shù)并保存到hr變量中。 hrLow?<- sprintf("%.3f",rt$"HR.95L") hrHigh <- sprintf("%.3f",rt$"HR.95H") Hazard.ratio <- paste0(hr,"(",hrLow,"-",hrHigh,")") #使用paste0函數(shù)將"HR"、"HR.95L"和"HR.95H"列的值合并為一個(gè)字符串， #形成風(fēng)險(xiǎn)比和置信區(qū)間的表示，保存到Hazard.ratio變量中。 pVal <- ifelse(rt$pvalue<0.001, "<0.001", sprintf("%.3f", rt$pvalue)) #用ifelse函數(shù)，如果rt數(shù)據(jù)框中的"pvalue"列小于0.001， #則將pVal變量設(shè)置為"<0.001"，否則將其設(shè)置為保留三位小數(shù)的字符串表示。 #輸出圖形pdf函數(shù)創(chuàng)建一個(gè)新的PDF繪圖設(shè)備，指定輸出文件名為forestFile pdf(file=forestFile, width=6.5, height=4.5) n <- nrow(rt) nRow <- n+1 ylim <- c(1,nRow) layout(matrix(c(1,2),nc=2),width=c(3,2.5)) #計(jì)算繪制圖形所需的總行數(shù)，包括基因和每個(gè)基因?qū)?yīng)的數(shù)據(jù)行。 設(shè)置y軸的范圍，從1到nRow。 #使用layout函數(shù)設(shè)置繪圖的布局，將繪制分為兩列，第一列寬度為3英寸，第二列寬度為2.5英寸。 #繪制森林圖左邊的臨床信息 xlim = c(0,3) par(mar=c(4,2.5,2,1)) #用par函數(shù)設(shè)置圖形的繪圖參數(shù)。mar參數(shù)用于設(shè)置邊緣空白的大小，包含4個(gè)值，依次為下、左、上、右的邊緣空白大小。 plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="") text.cex=0.8 text(0,n:1,gene,adj=0,cex=text.cex) text(1.5-0.5*0.2,n:1,pVal,adj=1,cex=text.cex);text(1.5-0.5*0.2,n+1,'pvalue',cex=text.cex,font=2,adj=1) text(3.1,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3.1,n+1,'Hazard ratio',cex=text.cex,font=2,adj=1) #圖中的坐標(biāo)(1.5-0.50.2, n:1)處添加p-value的值，并通過adj=1參數(shù)將文本右對齊，cex=text.cex設(shè)置字體縮放。 #然后在坐標(biāo)(1.5-0.50.2, n+1)處添加'pvalue'文本，并設(shè)置字體加粗。 #繪制右邊的森林圖 par(mar=c(4,1,2,1),mgp=c(2,0.5,0)) xlim = c(0,max(as.numeric(hrLow),as.numeric(hrHigh))) plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Hazard ratio") arrows(as.numeric(hrLow),n:1,as.numeric(hrHigh),n:1,angle=90,code=3,length=0.05,col="darkblue",lwd=3) abline(v=1, col="black", lty=2, lwd=2) boxcolor = ifelse(as.numeric(hr) > 1, forestCol, forestCol) points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=2) axis(1) dev.off() }

############定義森林圖函數(shù)############

? #定義獨(dú)立預(yù)后分析函數(shù) #定義了一個(gè)名為indep的函數(shù)，該函數(shù)有六個(gè)參數(shù)：expFile、cliFile、uniOutFile、multiOutFile、uniForest和multiForest， #這些參數(shù)可以在函數(shù)被調(diào)用時(shí)傳遞實(shí)際值。 indep=function(expFile=null,cliFile=null,uniOutFile=null,multiOutFile=null,uniForest=null,multiForest=null){

exp=read.table(expFile, header=T, sep="\t", check.names=F, row.names=1)?????#讀取表達(dá)文件

cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)?????#讀取臨床文件

#用read.table函數(shù)從expFile中讀取數(shù)據(jù)，數(shù)據(jù)以Tab分隔，包含表頭（header=T）， #第一列作為行名（row.names=1），并將數(shù)據(jù)保存到名為exp的數(shù)據(jù)框中。 #用read.table函數(shù)從cliFile中讀取數(shù)據(jù)，數(shù)據(jù)以Tab分隔，包含表頭（header=T）， #第一列作為行名（row.names=1），并將數(shù)據(jù)保存到名為cli的數(shù)據(jù)框中。 #數(shù)據(jù)合并 sameSample=intersect(row.names(cli),row.names(exp)) #獲取兩個(gè)數(shù)據(jù)框cli和exp行名（樣本名稱）的交集，并保存到sameSample變量中。 exp=exp[sameSample,] cli=cli[sameSample,] rt=cbind(exp, cli)

#單因素獨(dú)立預(yù)后分析

uniTab=data.frame() for(i in colnames(rt[,3:ncol(rt)])){ ?cox <- coxph(Surv(futime, fustat) ~ rt[,i], data = rt) ?#使用coxph函數(shù)進(jìn)行Cox比例風(fēng)險(xiǎn)模型分析，其中Surv(futime, fustat)表示生存時(shí)間和事件發(fā)生情況， ?#rt[,i]表示當(dāng)前循環(huán)的列作為預(yù)測因子，data=rt表示使用rt數(shù)據(jù)框。 ?coxSummary = summary(cox) ?uniTab=rbind(uniTab, ??????????????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|)"]) ??????????????) }

write.table(uniTab,file=uniOutFile,sep="\t",row.names=F,quote=F)

bioForest(coxFile=uniOutFile, forestFile=uniForest, forestCol="green")

#多因素獨(dú)立預(yù)后分析

uniTab=uniTab[as.numeric(uniTab[,"pvalue"])<1,] rt1=rt[,c("futime", "fustat", as.vector(uniTab[,"id"]))] #從rt數(shù)據(jù)框中選擇"futime"、"fustat"和uniTab數(shù)據(jù)框中的"pvalue"列對應(yīng)的預(yù)測因子列， #并將它們合并為一個(gè)新的數(shù)據(jù)框rt1。 multiCox=coxph(Surv(futime, fustat) ~ ., data = rt1) multiCoxSum=summary(multiCox) multiTab=data.frame() multiTab=cbind( ?????????????HR=multiCoxSum$conf.int[,"exp(coef)"], ?????????????HR.95L=multiCoxSum$conf.int[,"lower .95"], ?????????????HR.95H=multiCoxSum$conf.int[,"upper .95"], ?????????????pvalue=multiCoxSum$coefficients[,"Pr(>|z|)"]) multiTab=cbind(id=row.names(multiTab),multiTab)

write.table(multiTab,file=multiOutFile,sep="\t",row.names=F,quote=F)

bioForest(coxFile=multiOutFile, forestFile=multiForest, forestCol="red")

}

#調(diào)用indep函數(shù)進(jìn)行獨(dú)立預(yù)后分析，傳遞相關(guān)的輸入文件和輸出文件名作為參數(shù)。 indep(expFile="expTime.txt", ??????cliFile="clinical.txt", ??????uniOutFile="uniCox.txt", ??????multiOutFile="multiCox.txt", ??????uniForest="uniForest.pdf", ??????multiForest="multiForest.pdf") ?

下期將為你帶來更多R語言的騷操作技巧