純代碼實(shí)操,帶你做基因與臨床的生存分析與校準(zhǔn)
我們經(jīng)常對基因表達(dá)數(shù)據(jù)與臨床數(shù)據(jù)進(jìn)行生存分析,今天小云給大家?guī)聿灰粯拥姆治鲱A(yù)測,他就是Nomogram。我們先了解一下概念理論再進(jìn)行操作。
Nomogram,也常稱為諾莫圖或者列線圖,在醫(yī)學(xué)領(lǐng)域的期刊出現(xiàn)頻率越來愈多,常用于評估腫瘤學(xué)和醫(yī)學(xué)的預(yù)后情況,可將Logistic回歸或Cox回歸的結(jié)果進(jìn)行可視化呈現(xiàn)。
校準(zhǔn)曲線(Calibration Curve)量度著預(yù)測準(zhǔn)確性。有助于確定在不同時(shí)間點(diǎn)上模型的預(yù)測能力如何。這樣的分析對于基因與臨床因素在腫瘤研究中的應(yīng)用具有重要意義。
言歸正傳,實(shí)操部分來了。下面是示例樣本的表達(dá)量信息和臨床數(shù)據(jù)文件。


#install.packages("survival")
#install.packages("regplot")
#install.packages("rms")
#引用包
library(survival)
library(regplot)
library(rms)
?
expFile="expTime.txt" ??????#表達(dá)數(shù)據(jù)文件
cliFile="clinical.txt" ?????#臨床數(shù)據(jù)文件
setwd() ????#設(shè)置工作目錄
#定義了幾個變量,包括"expFile"和"cliFile",分別為表達(dá)數(shù)據(jù)文件和臨床數(shù)據(jù)文件的文件名。
#另外,通過"setwd"函數(shù)設(shè)置了工作目錄
?
exp=read.table(expFile, header=T, sep="\t", check.names=F, row.names=1)
#讀取名為"expFile"的表達(dá)數(shù)據(jù)文件,該文件是一個以制表符分隔的文本文件,包含基因表達(dá)矩陣。
#"header=T"表示第一行是列名,"sep="\t""表示使用制表符作為分隔符,"check.names=F"表示不檢查列名的合法性,
#"row.names=1"表示使用第一列作為行名。讀取后的數(shù)據(jù)被保存在變量"exp"中,是一個包含基因表達(dá)矩陣的數(shù)據(jù)框。
?
cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)
cli=cli[apply(cli,1,function(x)any(is.na(match('unknow',x)))),,drop=F]
cli$Age=as.numeric(cli$Age)
#讀取名為"cliFile"的臨床數(shù)據(jù)文件,該文件也是一個以制表符分隔的文本文件,包含與每個樣本相關(guān)的臨床信息。
#"header=T"表示第一行是列名,"sep="\t""表示使用制表符作為分隔符,"check.names=F"表示不檢查列名的合法性,
#"row.names=1"表示使用第一列作為行名。
#然后,代碼使用apply函數(shù)過濾出含有"unknow"的行,再將"Age"列轉(zhuǎn)換為數(shù)值類型。
?
samSample=intersect(row.names(exp), row.names(cli))
exp1=exp[samSample,,drop=F]
cli=cli[samSample,,drop=F]
rt=cbind(exp1, cli)
#合并表達(dá)數(shù)據(jù)"exp"和臨床數(shù)據(jù)"cli"。首先,找出兩個數(shù)據(jù)框中共同具有的樣本(行名),并將這些樣本作為"exp"和"cli"的子集。
#然后使用cbind函數(shù)將"exp1"和"cli"按列合并,得到名為"rt"的新數(shù)據(jù)框,其中包含了基因表達(dá)數(shù)據(jù)和臨床信息。
?
res.cox=coxph(Surv(futime, fustat) ~ . , data = rt)
nom1=regplot(res.cox,
??????????????plots = c("density", "boxes"),
??????????????clickable=F,
??????????????title="",
??????????????points=TRUE,
??????????????droplines=TRUE,
??????????????observation=rt[2,],
??????????????rank="sd",
??????????????failtime = c(1,3,5),
??????????????prfail = F)
#進(jìn)行生存分析。首先,使用coxph函數(shù)進(jìn)行Cox比例風(fēng)險(xiǎn)模型擬合,其中"Surv(futime, fustat)"是生存時(shí)間和事件(死亡或事件)的生存對象。
#擬合結(jié)果保存在"res.cox"變量中。
#然后,使用"regplot"函數(shù)繪制列線圖(Nomogram圖),顯示基因表達(dá)和臨床因素與生存風(fēng)險(xiǎn)的關(guān)系。
?
nomoRisk=predict(res.cox, data=rt, type="risk")
rt=cbind(exp1, Nomogram=nomoRisk)
outTab=rbind(ID=colnames(rt), rt)
write.table(outTab, file="nomoRisk.txt", sep="\t", col.names=F, quote=F)

#輸出列線圖的風(fēng)險(xiǎn)打分文件,計(jì)算并輸出列線圖(Nomogram圖)中的風(fēng)險(xiǎn)得分。
#使用predict函數(shù)基于Cox模型計(jì)算每個樣本的風(fēng)險(xiǎn)得分,并將結(jié)果保存在"nomoRisk"變量中。
#接著,將風(fēng)險(xiǎn)得分與表達(dá)數(shù)據(jù)合并為新的數(shù)據(jù)框"rt",再將數(shù)據(jù)框的列名和數(shù)據(jù)寫入名為"nomoRisk.txt"的文本文件中。
#校準(zhǔn)曲線
pdf(file="calibration.pdf", width=5, height=5)
#繪制了校準(zhǔn)曲線(Calibration Curve),用于評估Nomogram模型的預(yù)測性能
#1年校準(zhǔn)曲線
f <- cph(Surv(futime, fustat) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=1)
cal <- calibrate(f, cmethod="KM", method="boot", u=1, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1),
?xlab="Nomogram-predicted OS (%)", ylab="Observed OS (%)", lwd=1.5, col="green", sub=F)
#3年校準(zhǔn)曲線
f <- cph(Surv(futime, fustat) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=3)
cal <- calibrate(f, cmethod="KM", method="boot", u=3, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="", lwd=1.5, col="blue", sub=F, add=T)
#5年校準(zhǔn)曲線
f <- cph(Surv(futime, fustat) ~ Nomogram, x=T, y=T, surv=T, data=rt, time.inc=5)
cal <- calibrate(f, cmethod="KM", method="boot", u=5, m=(nrow(rt)/3), B=1000)
plot(cal, xlim=c(0,1), ylim=c(0,1), xlab="", ylab="", ?lwd=1.5, col="red", sub=F, add=T)
legend('bottomright', c('1-year', '3-year', '5-year'),
???col=c("green","blue","red"), lwd=1.5, bty = 'n')
dev.off()
#使用cph函數(shù)建立Cox模型,將Nomogram作為解釋變量,并設(shè)定不同的時(shí)間間隔。
#然后,使用calibrate函數(shù)生成校準(zhǔn)曲線。代碼將校準(zhǔn)曲線繪制在一個PDF文件中,
#分別繪制1年、3年和5年的校準(zhǔn)曲線,并在圖例中顯示不同時(shí)間間隔的顏色。最后,通過dev.off()函數(shù)關(guān)閉PDF繪圖設(shè)備。


綜合分析基因表達(dá)數(shù)據(jù)和臨床數(shù)據(jù),然后利用生存分析和校準(zhǔn)曲線評估基因表達(dá)和臨床信息對腫瘤患者預(yù)后的預(yù)測能力,預(yù)測患者預(yù)后以及發(fā)現(xiàn)潛在的生物標(biāo)志物。

