如何用GEO数据集进行批量基因的COX回归分析

生信菜鸟团

发布于 2021-03-23 15:01:39

5.7K40

代码可运行

文章被收录于专栏：生信菜鸟团生信菜鸟团

运行总次数：0

代码可运行

在进行数据挖掘过程中，我们往往会有对于所筛选出来的目标基因判断他们与预后之间的关系，这是我们就需要进行COX回归分析。下面以GEO数据库GSE62254这部分胃癌数据为例，分析其基本过程。

STEP1：获取目标数据GSE62254的基因表达矩阵expr及预后信息survival_file

基因表达矩阵的获取这里有两种方式一种如下图所示直接通过网页进行下载，

第二种则是通过R 包实现，代码如下：

library(GEOquery) #需提前安装

gset <- getGEO('GSE62254',destdir=".",
               AnnotGPL = F,     ## 注释文件
               getGPL = F)
save(gset,file="gset.rda")
gset <- gset[[1]] #降级处理
exprSet <- exprs(gset) #获取表达矩阵
dim(exprSet)

这样我们就得到了样本和基因组成的表达矩阵：

进而可以根据自己的需求只保留自己的目标基因。

预后信息的获取则比较灵活，在数据库网页可能存在下载链接也有可能像本例一样存在于数据库所属文章的附属文件里

对于预后信息我们只需关注与生存死亡以及生存时间相关的两列OS及OS.time,所以我们需要整理预后信息对样本信息及其对应的OS及OS.time进行保留，并且读入我们的工作环境。

library(readxl)
survival_file <- read_excel("D:/生信/bio paper/3rd 数据挖掘 gastric cancer/survival_file.xlsx")
View(survival_file)

继而通过merge函数，通过GSM_ID将目标基因表达矩阵以及预后信息进行融合，得到可以进行回归分析的目标矩阵data

survival_file <-survival_file[row.names(survival_file)%in%colnames(expr),] #对两者进行匹配
expr <- expr[,colnames(expr)%in%row.names(survival_file)]
dim(survival_file)
dim(expr)
expr <- as.data.frame(expr)
data <-merge.data.frame(expr,survival_file,by.x="sample",by.y ="GEO_ID")

STEP2 COX 回归分析及森林图绘制

通过一个for循环对所有目标基因进行回归分析，并且以dataframe的形式对结果进行输出：

for(i in colnames(data[,4:ncol(data)])){
  cox<- coxph(Surv(OS.time, OS) ~ get(i), data = data)
 coxSummary = summary(cox)
 result=rbind(result,
               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|)"]))
}
result[,2:5] <-apply(result[,2:5],2,as.numeric)

通过P值以及HR对有预后意义的基因进行筛选

table(result$pvalue<0.05)

森林图绘制：#读取输入文件

if(T){
  rt<- result3
  rt2<- cbind(rt[,c(2,3,4)],rt[,5])
 names(rt2) <- c('HR','HR.95L','HR.95H','pvalue')
  rt<- rt2
 gene <- gsub('[``]','',result3$id)
  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))

#输出图形

 if(T){

    n<- nrow(rt)
   nRow <- n+1
   ylim <- c(1,nRow)
   layout(matrix(c(1,2),nc=2),width=c(3,2))
   
    #绘制森林图左边的基因信息
   xlim = c(0,2.6)
   par(mar=c(4,2.5,2,0))
   plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")
   text.cex=0.6
    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.15,n+1,'pvalue',cex=text.cex,font=2,adj=1)
   text(2.5,n:1,Hazard.ratio,adj=1,cex=text.cex);text(2.45,n+1,'Hazardratio',cex=text.cex,font=2,adj=1,)

#绘制森林图

 par(mar=c(4,0,2,1),mgp=c(2,0.5,0))
   xlim = c(0,5)
   plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Hazardratio")
   arrows(as.numeric(hrLow),n:1,as.numeric(hrHigh),n:1,angle=0,code=3,length=0.05,col="black",lwd=2.5)
   abline(v=1,col="gray",lty=2,lwd=1.5)
   boxcolor = '#104E8B'
   points(as.numeric(hr), n:1, pch = 15, cex=0.6,col = boxcolor)
   axis(1)
  }}