生信分析（1）：单变量+多变量COX分析

从TCGA上下载数据库和临床数据之后，往往需要进行COX分析，一般的分析思路是先进行单变量，在进行多变量的分析。然而，当关注的基因比较多是，手动输入就会比较麻烦。接下来介绍一种利用循环的方法，快速的对多个变量进行分析。

首先是导入数据，包括基因表达counts数据和临床数据sur，及autophage基因集（来自HADb : Human Autophagy Database，参考文章《A risk signature with four autophagy-related genes for predicting survival of glioblastoma multiforme》），可根据需要替换为其他需要分析的基因列表，以及要用到的包：

setwd("D:/A1/Rdata/Autophage/胰腺癌")
library("survival")
library("survminer")
library("clusterProfiler")
library("DO.db")

#输入数据，data为数据矩阵，sur为生存数据，Autophage为自噬相关基因表
data <- read.csv("data.csv",header = T)
sur <- read.csv("sur.csv",header = T)
autophage <- read.csv("Autophage.csv",header = T)

#去掉重复值
index <- duplicated(sur$Sample) 
sur <- sur[!index,]

data格式如下：

sur数据格式如下：

对data进行标准化：

#log2标准化
row <- row.names(data)
data <- as.matrix(data)
data <- apply(data,2, as.numeric)
qx <- as.numeric(quantile(data, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0)
if (LogC) {
  data[which(data <= 0)] <- NA
  data <- log2(data+1)
}
data <- as.data.frame(data)
row.names(data) <- row

因为data中的基因格式是EMSEMBL，所以将autophage中的基因SYMBOL转换为EMSEMBL：

autogn <- autophage$Symbol
autogn <- bitr(autogn, fromType = "SYMBOL",toType = c("ENSEMBL"),OrgDb = "org.Hs.eg.db")

随后进行单变量cox分析：

基本思路是，将autophage中的基因EMSEMBL一个个取出，与data中的基因对应，并将该基因按中位数分为高表达组（1）和低表达组（0），与生存数据存入同一个数据框中（sur_r）,随后进行单因素cox分析，并将结果存入unicox中。

a <- colnames(data)
a <- gsub("\\.","-",a) #替换字符串
i <- length(a)
while (i > 0) {a[i] <- substr(a[i],1,12);i <- i-1}
r <- length(row.names(autogn))
k <- length(data)
data$name <- row.names(data)
unicox <- data.frame(geneID = NA,HR = NA,downCI = NA,upCI = NA,P.val = NA)
s = 1

#自动将data中的自噬基因匹配并进行单因素cox分析
while(r > 0){
  gene <- autogn[r,2]
  if (gene %in% row.names(data)){
    b <- data[data$name==gene,1:k]
    b <- as.numeric(b)
    ab <- data.frame(Sample = a,group = b)
    sur_r <- merge(sur,ab,by = "Sample")
    avg <- mean(b)
    sur_r$group[sur_r$group < avg] <- 0
    sur_r$group[sur_r$group >= avg] <- 1
    n <- sur_r$days_to_death
    n <- as.numeric(n)
    sur_r$days_to_death <- n
    n <- sur_r$status
    n <- as.numeric(n)
    sur_r$status <- n
    n <- sur_r$group
    n <- as.numeric(n)
    sur_r$group <- n
    res.cox <- coxph(Surv(days_to_death, status) ~ group, data = sur_r)
    res <- summary(res.cox)
    unicox[s,1] <- gene
    unicox[s,2] <- res$conf.int[1]
    unicox[s,3] <- res$conf.int[3]
    unicox[s,4] <- res$conf.int[4]
    unicox[s,5] <- res$waldtest[3]
    s = s + 1
  }
  ;r <- r-1
}

unicox格式如下：

存储了基因名，HR，95%CI和P值，下面需要选出P值<0.05的基因进入多因素cox分析，最开始纳入了16个变量，经过筛选后仅剩下6个：

#筛选有统计学意义的变量
unicox_sub <- unicox[unicox$P.val < 0.05,]
write.csv(unicox,"unicox_new.csv",row.names = F)
write.csv(unicox_sub,"unicox_select_new.csv",row.names = F)

#多元cox回归分析
p_max <- 1
muticox_list <- unicox_sub$geneID
i <- length(muticox_list)
l <- 19
c <- colnames(sur_r)
while (i > 0) {
  gene <- muticox_list[i]
  b <- data[data$name==gene,1:k]
  b <- as.numeric(b)
  avg <- mean(b)
  b[b < avg] <- 0
  b[b >= avg] <- 1
  ab <- data.frame(Sample = a,group = b)
  sur_me <- merge(sur,ab,by = "Sample")
  sur_r[,l] <- sur_me$group
  c[l] <- gene
  l <- l + 1
  ;i <- i - 1
}
colnames(sur_r) <- c
l <- length(muticox_list)
muticox_list[l+1] <- "days_to_death"
muticox_list[l+2] <- "status"
sur_muti <- subset(sur_r,select = muticox_list)
res.muticox <- coxph(Surv(days_to_death, status) ~ ., data = sur_muti)
sum_muticox <- summary(res.muticox)
muticox_result <- cbind(sum_muticox[["conf.int"]],sum_muticox[["coefficients"]])
muticox_result <- as.data.frame(muticox_result)
colnames(muticox_result)[9] <- "p.val"
p_max <- max(muticox_result$p.val)
while (p_max > 0.15){
  muticox_result <- muticox_result[order(muticox_result$p.val,decreasing = T),]
  muticox_result <- muticox_result[-1,]
  muticox_list <- row.names(muticox_result)
  i <- length(muticox_list)
  l <- 19
  c <- colnames(sur_r)
  while (i > 0) {
    gene <- muticox_list[i]
    b <- data[data$name==gene,1:k]
    b <- as.numeric(b)
    avg <- median(b)
    b[b < avg] <- 0
    b[b >= avg] <- 1
    ab <- data.frame(Sample = a,group = b)
    sur_me <- merge(sur,ab,by = "Sample")
    sur_r[,l] <- sur_me$group
    c[l] <- muticox_list[i]
    l <- l + 1
    ;i <- i - 1
  }
  colnames(sur_r) <- c
  l <- length(muticox_list)
  muticox_list[l+1] <- "days_to_death"
  muticox_list[l+2] <- "status"
  sur_muti <- subset(sur_r,select = muticox_list)
  res.muticox <- coxph(Surv(days_to_death, status) ~ ., data = sur_muti)
  sum_muticox <- summary(res.muticox)
  muticox_result <- cbind(sum_muticox[["conf.int"]],sum_muticox[["coefficients"]])
  muticox_result <- as.data.frame(muticox_result)
  colnames(muticox_result)[9] <- "p.val"
  p_max <- max(muticox_result$p.val)
}

gene_df <- row.names(muticox_result)
gene_df <- bitr(gene_df, fromType = "ENSEMBL",toType = c("SYMBOL"),OrgDb = "org.Hs.eg.db")
row.names(muticox_result) <- gene_df$SYMBOL
write.csv(muticox_result,"muticox_result_new.csv")

代码比较长，基本思路是先进行一次多因素cox分析，得到模型中P值最大的变量，随后将其存在p_max中，并将改变了剔除，随后利用循环不断剔除p_max最大的变量，直到p_max<0.15，也就是多变量分析的向后选择法，剔除标准为0.15，随后得到结果：

最终剩余6个变量，对这些变量绘制森林图：

muticox_list[1:l] <- gene_df$SYMBOL
colnames(sur_muti) <- muticox_list
res.muticox <- coxph(Surv(days_to_death, status) ~ ., data = sur_muti)

#绘制森林图
tiff(file="muticox_result_new.tiff",
     width = 20,          
     height = 20,            
     units ="cm",
     compression="lzw",
     bg="white",
     res=500)
ggforest(res.muticox)
dev.off()

 以上便是全过程了，巧用循环语句，可以减少生信分析中很多重复的工作。

最后还有一点疑问，如何根据上述的结果构建预测癌症预后的风险模型呢？欢迎大家参与讨论！