從TCGA上下載數據庫和臨床數據之后，往往需要進行COX分析，一般的分析思路是先進行單變量，在進行多變量的分析。然而，當關注的基因比較多是，手動輸入就會比較麻煩。接下來介紹一種利用循環的方法，快速的對多個變量進行分析。

首先是導入數據，包括基因表達counts數據和臨床數據sur，autophage是我下的一個自噬基因集，可根據需要替換為其他需要分析的基因列表，以及要用到的包：

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)] <- NAdata <- 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] <- 0sur_r$group[sur_r$group >= avg] <- 1n <- sur_r$days_to_deathn <- as.numeric(n)sur_r$days_to_death <- nn <- sur_r$statusn <- as.numeric(n)sur_r$status <- nn <- sur_r$groupn <- as.numeric(n)sur_r$group <- nres.cox <- coxph(Surv(days_to_death, status) ~ group, data = sur_r)res <- summary(res.cox)unicox[s,1] <- geneunicox[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] <- 0b[b >= avg] <- 1ab <- data.frame(Sample = a,group = b)sur_me <- merge(sur,ab,by = "Sample")sur_r[,l] <- sur_me$groupc[l] <- genel <- 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 <- 19c <- 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] <- 0b[b >= avg] <- 1ab <- data.frame(Sample = a,group = b)sur_me <- merge(sur,ab,by = "Sample")sur_r[,l] <- sur_me$groupc[l] <- muticox_list[i]l <- l + 1;i <- i - 1}colnames(sur_r) <- cl <- 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()

?以上便是全過程了，巧用循環語句，可以減少生信分析中很多重復的工作。

最后還有一點疑問，如何根據上述的結果構建預測癌癥預后的風險模型呢？歡迎大家參與討論！

總結

以上是生活随笔為你收集整理的生信分析（1）：单变量+多变量COX分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。