批量单因素Cox回归函数封装

欢迎关注交流

Step0：加载包加载数据

load("data/sur_tcga_icgc_os.Rdata")
names(expr_icgc_os)
library(survival)
library(tidyverse)

Step1:小试（单基因）

cox_model <- coxph(Surv(os, Event) ~ ATIC, data = expr_icgc_os)
HR <- round(summary(cox_model)$coef[, "exp(coef)"], 2)
lower <- round(summary(cox_model)$conf.int[,"lower .95"], 2)
upper <- round(summary(cox_model)$conf.int[,"upper .95"], 2)
pvalue <- format(round(summary(cox_model)$coefficients[, "Pr(>|z|)"], digits = 10), scientific=TRUE)

Step2:小试（数据集）

2.1 对每个基因进行单因素 Cox 回归分析，并存储结果

result <- lapply(names(expr_icgc_os)[1:158], function(gene) {
  cox_model <- coxph(Surv(os, Event) ~ expr_icgc_os[[gene]], data = expr_icgc_os)
  HR <- round(summary(cox_model)$coef[, "exp(coef)"], 2)
  lower <- round(summary(cox_model)$conf.int[,"lower .95"], 2)
  upper <- round(summary(cox_model)$conf.int[,"upper .95"], 2)
  p_value <- format(round(summary(cox_model)$coefficients[, "Pr(>|z|)"], digits = 10), scientific=TRUE)
  c(HR = HR, ci_lower = lower, ci_upper = upper, p_value = p_value)
})
names(result) <- names(expr_icgc_os)[1:158]

2.2 整理结果至数据框中

cox_results <- data.frame(Gene = names(result), 
                          HR = sapply(result, function(x) x["HR"]), 
                          CI_lower = sapply(result, function(x) x["ci_lower"]), 
                          CI_upper = sapply(result, function(x) x["ci_upper"]), 
                          P_value = sapply(result, function(x) x["p_value"])) %>%
  mutate(Adjusted_P_value = p.adjust(P_value, method = "fdr"))
view(cox_results)

Step3：定义函数计算

calculate_cox_values <- function(data, time_var, event_var, gene_names) {
  result <- lapply(gene_names, function(gene) {
    cox_model <- coxph(Surv(data[[time_var]], data[[event_var]]) ~ data[[gene]], data = data)
    HR <- round(summary(cox_model)$coef[, "exp(coef)"], 2)
    lower <- round(summary(cox_model)$conf.int[,"lower .95"], 2)
    upper <- round(summary(cox_model)$conf.int[,"upper .95"], 2)
    p_value <- format(round(summary(cox_model)$coefficients[, "Pr(>|z|)"], digits = 10), scientific=TRUE)
    c(HR = HR, ci_lower = lower, ci_upper = upper, p_value = p_value)
  })
  
  names(result) <- gene_names
  
  cox_results <- data.frame(Gene = names(result), 
                            HR = sapply(result, function(x) x["HR"]), 
                            CI_lower = sapply(result, function(x) x["ci_lower"]), 
                            CI_upper = sapply(result, function(x) x["ci_upper"]), 
                            P_value = sapply(result, function(x) x["p_value"])) %>%
    mutate(Adjusted_P_value = p.adjust(P_value, method = "fdr"))
  
  return(cox_results)
}

# 使用示例
result_table <- calculate_cox_values(expr_icgc_os, "os", "Event", names(expr_icgc_os)[1:158])
view(result_table)

Step4：优化代码

uni_cox_values <- function(data, time, event) {
  library(survival)
  library(dplyr)
  
  univariate_cox_var <- setdiff(names(data), c(time, event))
  
  result <- lapply(univariate_cox_var, function(univariate_cox_var) {
    cox_model <- coxph(Surv(data[[time]], data[[event]]) ~ data[[univariate_cox_var]], data = data)
    HR <- format(round(summary(cox_model)$coef[, "exp(coef)"], 2), nsmall = 2)
    lower <- format(round(summary(cox_model)$conf.int[,"lower .95"], 2), nsmall = 2)
    upper <- format(round(summary(cox_model)$conf.int[,"upper .95"], 2), nsmall = 2)
    p_value <- format(round(summary(cox_model)$coefficients[, "Pr(>|z|)"], digits = 6), scientific=TRUE)
    c(HR = HR, ci_lower = lower, ci_upper = upper, p_value = p_value)
  })
  
  names(result) <- univariate_cox_var
  
  cox_results_table <- data.frame(Variable = names(result), 
                                  HR = sapply(result, function(x) x["HR"]), 
                                  ci_lower = sapply(result, function(x) x["ci_lower"]), 
                                  ci_upper = sapply(result, function(x) x["ci_upper"]), 
                                  p_value = sapply(result, function(x) x["p_value"])) %>%
    mutate(adjusted_p_value = p.adjust(p_value, method = "fdr"))
    #mutate(Adjusted_P_value = round(p.adjust(P_value, method = "fdr"), 3))
    
    ##使用gsub函数去掉行名中的.HR
    #rownames(cox_icgc) <- gsub("\\.HR", "", rownames(cox_icgc))
  
  return(cox_results_table)
}
## 使用示例
## icgc
cox_icgc <- uni_cox_values(expr_icgc_os, "os", "Event")

## tcga
cox_tcga <- uni_cox_values(expr_tcga_os, "os", "Event")

## 其他数据............