小云又為大家來(lái)進(jìn)行知識(shí)分享啦！今天給大家?guī)?lái)的是利用單因素Cox回歸分析，來(lái)篩選有預(yù)后意義的特征基因，接下來(lái)跟著小云開(kāi)始今天的學(xué)習(xí)。

1.?如何進(jìn)行預(yù)后基因篩選？

小伙伴們是不是跟小云存在同樣的疑問(wèn)？不急！讓小云來(lái)為大家解答這個(gè)疑惑！一般利用survival包對(duì)帶有生存時(shí)間和生存狀態(tài)的基因表達(dá)矩陣文件進(jìn)行單因素Cox回歸分析，通過(guò)計(jì)算的Pvalue值來(lái)篩選有預(yù)后意義的特征基因，然后用森林圖來(lái)可視化該結(jié)果；這就是基本的分析過(guò)程，其實(shí)很簡(jiǎn)單，只要數(shù)據(jù)格式合適，就可以很快很輕松完成分析，今天小云通過(guò)循環(huán)的方式批量對(duì)每個(gè)基因進(jìn)行Cox回歸分析，代碼實(shí)用性非常不錯(cuò)奧！那就和小云一起開(kāi)始今天的實(shí)操吧！

#安裝需要的R包

insatall.packages("survival")

#加載需要的R包

library(survival)

2.?準(zhǔn)備數(shù)據(jù)

cox.txt

#帶有生存時(shí)間和生存狀態(tài)的基因表達(dá)矩陣文件，行名為樣本名，第一列表示生存時(shí)間，第二列表示生存狀態(tài)，其他列表示基因所對(duì)應(yīng)的表達(dá)量數(shù)據(jù)。

3.?cox單因素回歸分析

#設(shè)置一個(gè)p值，對(duì)cox分析結(jié)果進(jìn)行篩選

pFilter=0.05

#建一個(gè)空白數(shù)據(jù)框，后面for循環(huán)輸出用

outResult=data.frame()

#建一個(gè)向量，后面for循環(huán)輸出用，因?yàn)楹竺孢€要用到surstat及surtime，所以先放在向量里

sigGenes=c("surstat","surtime")

#從第3列開(kāi)始循環(huán)，因?yàn)?列2列不是gene，是surstat和surtime

for(i in colnames(td[,3:ncol(td)])){

?#開(kāi)始逐一循環(huán)cox分析

tdcox <- coxph(Surv(surtime, surstat) ~ td[,i], data = td)

? #summary命令對(duì)tdcox總結(jié)，方便后面提取數(shù)據(jù)

tdcoxSummary = summary(tdcox)

#提取p值，這個(gè)主要是后面提取有意義的gene用

?pvalue=tdcoxSummary$coefficients[,"Pr(>|z|)"]

? 這里我們不需要提取所有基因的數(shù)據(jù)，只需要有意義的gene的結(jié)果，所以設(shè)置如果pvalue<0.05才提取數(shù)據(jù)

if(pvalue<pFilter){

?????sigGenes=c(sigGenes,i)

????#?合并行，實(shí)際上是對(duì)循環(huán)結(jié)果的合并，前面設(shè)置的空白數(shù)據(jù)框outResult這里用，循環(huán)必須有個(gè)開(kāi)始

outResult=rbind(outResult,

???????????? cbind(id=i,#合并列，是每個(gè)基因的統(tǒng)計(jì)數(shù)據(jù),#提取單個(gè)基因的HR ??????????????????????????

HR=tdcoxSummary$conf.int[,"exp(coef)"],

#提取單個(gè)基因的HR的95%CI低值

L95CI=tdcoxSummary$conf.int[,"lower .95"],

#提取單個(gè)基因的HR的95%CI高值

H95CI=tdcoxSummary$conf.int[,"upper .95"],

#提取單個(gè)基因的p值

pvalue=tdcoxSummary$coefficients[,"Pr(>|z|)"])

??????????????????????????????)

??}

}

write.table(outResult,file="UniCoxSurvival.txt",sep="\t",row.names=F,quote=F)

#以有統(tǒng)計(jì)學(xué)意義的gene取表達(dá)數(shù)據(jù)的子集

UniCoxSurSigGeneExp=td[,sigGenes]

#以id也就是樣品名命名行名

UniCoxSurSigGeneExp=cbind(id=row.names(UniCoxSurSigGeneExp),UniCoxSurSigGeneExp)

#保存結(jié)果文件

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

4.?繪制隨機(jī)森林圖

#提取繪圖相關(guān)數(shù)據(jù)

gene <- rownames(tducs)#提取基因名

hr <- sprintf("%.3f",tducs$"HR")#從tducs數(shù)據(jù)中提取HR值，%.3f指數(shù)據(jù)保留小數(shù)點(diǎn)后3位

hrLow ?<- sprintf("%.3f",tducs$"L95CI")#提取HR值95%置信區(qū)間的低值

hrHigh <- sprintf("%.3f",tducs$"H95CI")#提取HR值95%置信區(qū)間的高值

Hazard.ratio <- paste0(hr,"(",hrLow,"-",hrHigh,")")#將HR及其置信區(qū)間數(shù)據(jù)組合成一個(gè)因子

pValue <- ifelse(tducs$pvalue<0.001, "<0.001", sprintf("%.3f", tducs$pvalue))#提取p值，且當(dāng)p<0.001時(shí)顯示為<0.001，而不是確切數(shù)值

pdf(file="UniCoxSurForestPlot.pdf", width = 16,height = 3)#pdf繪圖命令開(kāi)始打印過(guò)程，設(shè)置圖片寬度6及高度3（自行根據(jù)數(shù)據(jù)設(shè)置）

#可以不執(zhí)行pdf繪圖命令，這樣執(zhí)行命令時(shí)圖像會(huì)一一展示出來(lái)，執(zhí)行pdf命令后圖像僅在最后輸出至pdf文件中

n <- nrow(tducs)#提取出tducs的行數(shù)，也就是有多少個(gè)基因

nRow <- n+1 #設(shè)置一個(gè)比基因數(shù)多1的數(shù)值

ylim <- c(1,nRow) #暫時(shí)設(shè)置一個(gè)y軸長(zhǎng)度值，長(zhǎng)度比基因數(shù)多1個(gè)，

layout(matrix(c(1,2),nc=2),width=c(2.5,2)) # 設(shè)置頁(yè)面排版，即1排兩圖，左圖1，右圖2。

#開(kāi)始畫(huà)圖的左邊部分，即圖1，森林圖的文字及數(shù)值標(biāo)注部分

xlim = c(0,2.5)#設(shè)置一個(gè)x軸長(zhǎng)度

par(mar=c(4,2.5,2,1))#設(shè)置圖形1的邊界，即;下4，左2.5,上2，右1

plot(0,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")#開(kāi)始繪制圖1，此時(shí)為空白圖，axex=F表示不顯示邊框，可以設(shè)置為T(mén)看看

text.cex=0.8 #設(shè)置圖中文本的大小

text(0,n:1,gene,adj=0,cex=text.cex)#添加文本，即在圖上添加一列基因名，從坐標(biāo)0開(kāi)始正向添加，n即前面定義的基因數(shù)，gene即基因名

text(1.4,n:1,pValue,adj=1,cex=text.cex);text(1.4,n+1,'pvalue',cex=text.cex,font=2,adj=1)#添加p值，adj=1從坐標(biāo)1.4反向添加

text(2.5,n:1,Hazard.ratio,adj=1,cex=text.cex);text(2.5,n+1,'Hazard ratio',cex=text.cex,font=2,adj=1,)#添加HR值，從坐標(biāo)2.5反向添加

#開(kāi)始畫(huà)圖的右邊部分，即圖2，森林圖的圖形

par(mar=c(4,1,2,1),mgp=c(2,0.5,0))#設(shè)置邊界，mpg設(shè)置橫標(biāo)題及坐標(biāo)軸的位置

xlim = c(0,max(as.numeric(hrLow),as.numeric(hrHigh))) #設(shè)置圖2的橫坐標(biāo)

plot(0,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Hazard ratio")#開(kāi)始畫(huà)圖2，axex=F表示不顯示邊框

arrows(as.numeric(hrLow),n:1,as.numeric(hrHigh),n:1,angle=90,code=3,length=0.05,col="darkblue",lwd=2.5)#HR的橫條

abline(v=1,col="black",lty=2,lwd=2)#在橫坐標(biāo)1的地方畫(huà)一條豎線

boxcolor = ifelse(as.numeric(hr) > 1, 'red', 'green')#設(shè)置顏色

points(as.numeric(hr), n:1, pch = 1, col = boxcolor, cex=1.3)#HR處顯示為圈，顏色為上面設(shè)的顏色。

axis(1)#顯示橫坐標(biāo)

#畫(huà)圖結(jié)束，輸出pdf

dev.off()

5.?結(jié)果文件

1.?UniCoxSurvival.txt

該結(jié)果文件為Cox回歸分析計(jì)算結(jié)果，第一列表示基因名，第二列表示HR值，第三列表示基因的HR的95%CI低值，第四列表示基因的HR的95%CI高值，第五列為Pvalue值。

2.?UniCoxSurSigGeneExp.txt

該結(jié)果文件為篩選的有預(yù)后意義的基因表達(dá)矩陣文件，第一列表示樣本名，第二列表示生存狀態(tài)，第三列表示生存時(shí)間，其他列為基因?qū)?yīng)的表達(dá)量數(shù)據(jù)。

3.?UniCoxSurForestPlot.pdf

該結(jié)果圖片為利用cox回歸分析結(jié)果來(lái)繪制森林圖。

今天小云的分享就到這里，其他相關(guān)其他分析內(nèi)容歡迎嘗試本公司新開(kāi)發(fā)的云平臺(tái)生物信息分析小工具，零代碼完成分析，云平臺(tái)網(wǎng)址：http://www.biocloudservice.com/home.html，包括cox回歸篩選預(yù)后標(biāo)記基因構(gòu)建模型(http://www.biocloudservice.com/128/128.php),逐步回歸法篩選基因構(gòu)建Cox風(fēng)險(xiǎn)模型(http://www.biocloudservice.com/124/124.php)等小工具，歡迎小伙伴們來(lái)嘗試奧！