独立预后分析：通过使用生存分析中的Cox比例风险模型（Cox proportional hazards model），评估基因表达与生存时间之间的关联。该分析可以帮助研究人员确定哪些基因对生存时间具有显著的影响。
森林图绘制：森林图是一种常用的图形表示方式，用于展示多个预后因素的效应量（如风险比或危险比）和显著性水平（如p值）。通过绘制森林图，可以直观地比较不同基因或预后因素对生存时间的影响。
单因素独立预后分析：代码使用coxph函数计算每个基因（或预后因素）与生存时间的关联，得到效应量（HR）和显著性水平（p值）。将结果保存到单因素预后分析结果输出文件（uniOutFile），然后调用bioForest函数绘制对应的单因素森林图。 多因素独立预后分析：根据单因素预后分析结果，选择p值小于1的基因，将它们作为预后因素进行多因素独立预后分析。这样可以考虑多个因素对生存时间的综合影响。多因素预后分析结果保存到多因素预后分析结果输出文件（multiOutFile），并绘制多因素森林图（调用bioForest函数）。
#install.packages(‘survival’)
library(survival) #引用包
setwd(“C:\\biowolf\\Gene\\16.indep”) #设置工作目录
定义森林图函数（bioForest）：用于绘制森林图，展示独立预后分析结果的效应量（Hazard Ratio）和显著性水平（p-value）。函数接受三个参数：coxFile（独立预后分析结果文件），forestFile（森林图输出文件），forestCol（森林图颜色）。函数读取独立预后分析结果文件，提取相关信息，并绘制森林图。
bioForest=function(coxFile=null, forestFile=null, forestCol=null){
#读取输入文件
rt <- read.table(coxFile, header=T, sep=”\t”, check.names=F, row.names=1)
gene <- rownames(rt)
hr <- sprintf(“%.3f”,rt$”HR”)
hrLow <- sprintf(“%.3f”,rt$”HR.95L”)
hrHigh <- sprintf(“%.3f”,rt$”HR.95H”)
Hazard.ratio <- paste0(hr,”(“,hrLow,”-“,hrHigh,”)”)
pVal <- ifelse(rt$pvalue<0.001, “<0.001”, sprintf(“%.3f”, rt$pvalue))
#输出图形
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))
#绘制森林图左边的临床信息
xlim = c(0,3)
par(mar=c(4,2.5,2,1))
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)
#绘制右边的森林图
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()
}
############定义森林图函数############
定义独立预后分析函数（indep）：该函数用于进行基因表达数据的独立预后分析。函数接受六个参数：expFile（表达数据文件），cliFile（临床数据文件），uniOutFile（单因素预后分析结果输出文件），multiOutFile（多因素预后分析结果输出文件），uniForest（单因素森林图输出文件），multiForest（多因素森林图输出文件）。
#定义独立预后分析函数
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) #读取表达文件
cli=read.table(cliFile, header=T, sep=”\t”, check.names=F, row.names=1) #读取临床文件
#数据合并，读取表达数据文件和临床数据文件，然后将它们合并为一个数据框。
sameSample=intersect(row.names(cli),row.names(exp))
exp=exp[sameSample,]
cli=cli[sameSample,]
rt=cbind(exp, cli)
#进行单因素独立预后分析，计算每个基因与生存风险的相关性，并将结果保存到单因素预后分析结果输出文件。
uniTab=data.frame()
for(i in colnames(rt[,3:ncol(rt)])){
cox <- coxph(Surv(futime, fustat) ~ rt[,i], data = rt)
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函数绘制单因素森林图，将预后分析结果可视化。
bioForest(coxFile=uniOutFile, forestFile=uniForest, forestCol=”green”)
#多因素独立预后分析，函数进行多因素独立预后分析，计算多个因素对生存风险的综合影响，并将结果保存到多因素预后分析结果输出文件。
uniTab=uniTab[as.numeric(uniTab[,”pvalue”])<1,]
rt1=rt[,c(“futime”, “fustat”, as.vector(uniTab[,”id”]))]
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”)
}
#调用bioForest函数绘制多因素森林图，展示多因素预后分析结果的效应量和显著性水平
#调用函数，进行独立预后分析。调用indep函数，传入相应的参数进行独立预后分析。根据输入的表达数据文件和临床数据文件，代码将进行单因素和多因素独立预后分析，并将结果以文件和森林图的形式输出。
indep(expFile=”expTime.txt”,
cliFile=”clinical.txt”,
uniOutFile=”uniCox.txt”,
multiOutFile=”multiCox.txt”,
uniForest=”uniForest.pdf”,
multiForest=”multiForest.pdf”)
来给大家展示展示结果图吧，使用firefox命令。
独立预后分析是一种常用的生存分析方法，用于评估基因表达等因素对生存时间的影响。通过该分析，可以识别与生存风险相关的基因，并了解基因对生存时间的影响程度。森林图则提供了一种可视化展示独立预后分析结果的方式，方便对结果进行观察和比较。
下期将为你带来更多R语言的骚操作技巧，以下推荐的是一个多功能的生信平台。
云生信平台链接：http://www.biocloudservice.com/home.html。
云生信平台链接：http://www.biocloudservice.com/home.html。