小云又为大家来进行知识分享啦！今天给大家带来的是利用单因素Cox回归分析，来筛选有预后意义的特征基因，接下来跟着小云开始今天的学习。

1. 如何进行预后基因筛选？

小伙伴们是不是跟小云存在同样的疑问？不急！让小云来为大家解答这个疑惑！一般利用survival包对带有生存时间和生存状态的基因表达矩阵文件进行单因素Cox回归分析，通过计算的Pvalue值来筛选有预后意义的特征基因，然后用森林图来可视化该结果；这就是基本的分析过程，其实很简单，只要数据格式合适，就可以很快很轻松完成分析，今天小云通过循环的方式批量对每个基因进行Cox回归分析，代码实用性非常不错奥！那就和小云一起开始今天的实操吧！

#安装需要的R包

insatall.packages("survival")

#加载需要的R包

library(survival)

2. 准备数据

cox.txt

#带有生存时间和生存状态的基因表达矩阵文件，行名为样本名，第一列表示生存时间，第二列表示生存状态，其他列表示基因所对应的表达量数据。

3. cox单因素回归分析

#设置一个p值，对cox分析结果进行筛选

pFilter=0.05

#建一个空白数据框，后面for循环输出用

outResult=data.frame()

#建一个向量，后面for循环输出用，因为后面还要用到surstat及surtime，所以先放在向量里

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

#从第3列开始循环，因为1列2列不是gene，是surstat和surtime

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

 #开始逐一循环cox分析

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

  #summary命令对tdcox总结，方便后面提取数据

tdcoxSummary = summary(tdcox)

#提取p值，这个主要是后面提取有意义的gene用

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

  这里我们不需要提取所有基因的数据，只需要有意义的gene的结果，所以设置如果pvalue<0.05才提取数据

if(pvalue<pFilter){

     sigGenes=c(sigGenes,i)

    # 合并行，实际上是对循环结果的合并，前面设置的空白数据框outResult这里用，循环必须有个开始

outResult=rbind(outResult,

             cbind(id=i,#合并列，是每个基因的统计数据,#提取单个基因的HR                           

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

#提取单个基因的HR的95%CI低值

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

#提取单个基因的HR的95%CI高值

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

#提取单个基因的p值

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

                              )

  }

}

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

#以有统计学意义的gene取表达数据的子集

UniCoxSurSigGeneExp=td[,sigGenes]

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

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

#保存结果文件

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

4. 绘制随机森林图

#提取绘图相关数据

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

hr <- sprintf("%.3f",tducs$"HR")#从tducs数据中提取HR值，%.3f指数据保留小数点后3位

hrLow  <- sprintf("%.3f",tducs$"L95CI")#提取HR值95%置信区间的低值

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

Hazard.ratio <- paste0(hr,"(",hrLow,"-",hrHigh,")")#将HR及其置信区间数据组合成一个因子

pValue <- ifelse(tducs$pvalue<0.001, "<0.001", sprintf("%.3f", tducs$pvalue))#提取p值，且当p<0.001时显示为<0.001，而不是确切数值

pdf(file="UniCoxSurForestPlot.pdf", width = 16,height = 3)#pdf绘图命令开始打印过程，设置图片宽度6及高度3（自行根据数据设置）

#可以不执行pdf绘图命令，这样执行命令时图像会一一展示出来，执行pdf命令后图像仅在最后输出至pdf文件中

n <- nrow(tducs)#提取出tducs的行数，也就是有多少个基因

nRow <- n+1 #设置一个比基因数多1的数值

ylim <- c(1,nRow) #暂时设置一个y轴长度值，长度比基因数多1个，

layout(matrix(c(1,2),nc=2),width=c(2.5,2)) # 设置页面排版，即1排两图，左图1，右图2。

#开始画图的左边部分，即图1，森林图的文字及数值标注部分

xlim = c(0,2.5)#设置一个x轴长度

par(mar=c(4,2.5,2,1))#设置图形1的边界，即;下4，左2.5,上2，右1

plot(0,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")#开始绘制图1，此时为空白图，axex=F表示不显示边框，可以设置为T看看

text.cex=0.8 #设置图中文本的大小

text(0,n:1,gene,adj=0,cex=text.cex)#添加文本，即在图上添加一列基因名，从坐标0开始正向添加，n即前面定义的基因数，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从坐标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值，从坐标2.5反向添加

#开始画图的右边部分，即图2，森林图的图形

par(mar=c(4,1,2,1),mgp=c(2,0.5,0))#设置边界，mpg设置横标题及坐标轴的位置

xlim = c(0,max(as.numeric(hrLow),as.numeric(hrHigh))) #设置图2的横坐标

plot(0,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Hazard ratio")#开始画图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)#在横坐标1的地方画一条竖线

boxcolor = ifelse(as.numeric(hr) > 1, 'red', 'green')#设置颜色

points(as.numeric(hr), n:1, pch = 1, col = boxcolor, cex=1.3)#HR处显示为圈，颜色为上面设的颜色。

axis(1)#显示横坐标

#画图结束，输出pdf

dev.off()

5. 结果文件

1. UniCoxSurvival.txt

该结果文件为Cox回归分析计算结果，第一列表示基因名，第二列表示HR值，第三列表示基因的HR的95%CI低值，第四列表示基因的HR的95%CI高值，第五列为Pvalue值。

2. UniCoxSurSigGeneExp.txt

该结果文件为筛选的有预后意义的基因表达矩阵文件，第一列表示样本名，第二列表示生存状态，第三列表示生存时间，其他列为基因对应的表达量数据。

3. UniCoxSurForestPlot.pdf

该结果图片为利用cox回归分析结果来绘制森林图。

今天小云的分享就到这里，其他相关其他分析内容欢迎尝试本公司新开发的云平台生物信息分析小工具，零代码完成分析，云平台网址：http://www.biocloudservice.com/home.html，包括cox回归筛选预后标记基因构建模型(http://www.biocloudservice.com/128/128.php),逐步回归法筛选基因构建Cox风险模型(http://www.biocloudservice.com/124/124.php)等小工具，欢迎小伙伴们来尝试奥！