batch cox-unicox and multicox

载入数据

BRCA_clinical<-read.csv("F:/Bioinfor_project/Breast/AS_research/AS/clinical/TNBC_clinical_final.csv",header=TRUE,sep=",")
rownames(BRCA_clinical)<-BRCA_clinical$patient_id
BRCA_clinical<-BRCA_clinical[,-1]
BRCA_clinical[1:5,1:5]

##              OS_time OS_event  X       ER PR_Status_nature2012
## TCGA_A1_A0SK     967        1 NA Negative             Negative
## TCGA_A2_A04P     548        1 NA Negative             Negative
## TCGA_A2_A0CM     754        1 NA Negative             Negative
## TCGA_A2_A0T2     255        1 NA Negative             Negative
## TCGA_A8_A09X     426        1 NA Negative             Negative

创建生存对象

library(survival)
require(survminer)

## Loading required package: survminer

## Loading required package: ggplot2

## Loading required package: ggpubr

## Loading required package: magrittr

y<-Surv(as.numeric(BRCA_clinical[,"OS_time"]),as.numeric(BRCA_clinical[,"OS_event"]))

## 临床数据因子化并添加标签
BRCA_clinical<-within(BRCA_clinical,{
  T<-factor(T,labels = c("T1/T2","T1/T2","T3/T4","T3/T4"))
  N<-factor(N,labels = c("N0/N1","N0/N1","N2/N3","N2/N3"))
  M<-factor(M,labels = c("M0","M1"))
  age<- factor(ifelse(age>=median(age),">=56","<56"),labels = c(">=56","<56"))
})
table(BRCA_clinical$T)

## 
## T1/T2 T3/T4 
##   101    13

table(BRCA_clinical$age)

## 
## >=56  <56 
##   56   58

单因素cox回归
函数的运行结果并不会自动到环境变量中

BRCA_T_cox.fit_dan <- coxph(y~BRCA_clinical[["N"]]) #Cox比例风险回归模型
BRCA_T_coxresult<-summary(BRCA_T_cox.fit_dan)#cox的结果
BRCA_T_coxresult

## Call:
## coxph(formula = y ~ BRCA_clinical[["N"]])
## 
##   n= 114, number of events= 18 
## 
##                             coef exp(coef) se(coef)     z Pr(>|z|)   
## BRCA_clinical[["N"]]N2/N3 1.7936    6.0108   0.5633 3.184  0.00145 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                           exp(coef) exp(-coef) lower .95 upper .95
## BRCA_clinical[["N"]]N2/N3     6.011     0.1664     1.993     18.13
## 
## Concordance= 0.663  (se = 0.065 )
## Likelihood ratio test= 7.93  on 1 df,   p=0.005
## Wald test            = 10.14  on 1 df,   p=0.001
## Score (logrank) test = 13  on 1 df,   p=3e-04

##BRCA_T_coxresult$coefficients[5]调P值,p=1.629449e-05
##BRCA_T_coxresult$conf.int[1]为HR，3为95%lower，4为95%upper

单因素cox批量计算函数unicox包装

输入感兴趣的变量，批量进行单因素cox分析，保存结果

##############################################################
### 包装批量绘制单因素cox的函数
unicox<-function(vars=c("T","N"),data=BRCA_clinical){
  library(survival)
  require(survminer)
   y<-Surv(as.numeric(data[,"OS_time"]),
           as.numeric(data[,"OS_event"]))
  
  ## for循环
  pvalue<-c()
  HR<-c()
  lower<-c()
  upper<-c()
  varname<-c()
  
  for(i in vars ){
    cox.fit_uni<-coxph(y~data[[i]])
    uni_res <-summary(cox.fit_uni)
    pvalue[[i]]<-round(uni_res$waldtest[[3]],3)#pvalue
    HR[[i]]<-round(uni_res$conf.int[[1]],2)# HR提取
    lower[[i]] <-round(uni_res$conf.int[[3]],2)#lower
    upper[[i]] <-round(uni_res$conf.int[[4]],2)# upper
  }
  ## 保存结果到数据框中
  univar_res<-data.frame(
    varname<-vars,
    HR<-as.numeric(HR),
    CI=paste(as.numeric(lower),"-",as.numeric(upper),sep = ""),
    pvalue<-as.numeric(pvalue)
  )
  colnames(univar_res)<-c("varname","HR","95%CI","pvalue")
  univar_res## 最后的函数运行结果
}
## 批量进行单因素cox回归分析
univar_res<-unicox(vars=c("T","N","M","age"),data=BRCA_clinical)
univar_res

##   varname    HR       95%CI pvalue
## 1       T  2.99   0.98-9.14  0.055
## 2       N  6.01  1.99-18.13  0.001
## 3       M 53.32 4.83-588.11  0.001
## 4     age  1.23   0.49-3.12  0.659

#write.csv(univar_res,file="F:/Bioinfor_project/Breast/AS_research/AS/result/univar_res_clinical.csv")

#####################################################################

TNM多因素cox

BRCA_TNM_cox.fit_dan <- coxph(y~T+N+M+age,data = BRCA_clinical)
BRCA_TNM_coxresult<-summary(BRCA_TNM_cox.fit_dan)#cox的结果
BRCA_TNM_coxresult

## Call:
## coxph(formula = y ~ T + N + M + age, data = BRCA_clinical)
## 
##   n= 114, number of events= 18 
## 
##           coef exp(coef) se(coef)      z Pr(>|z|)   
## TT3/T4  0.8282    2.2891   0.6652  1.245  0.21314   
## NN2/N3  1.6330    5.1191   0.6226  2.623  0.00872 **
## MM1     2.1104    8.2519   1.3762  1.534  0.12515   
## age<56 -0.1859    0.8304   0.5191 -0.358  0.72027   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## TT3/T4    2.2891     0.4368    0.6215     8.431
## NN2/N3    5.1191     0.1953    1.5110    17.343
## MM1       8.2519     0.1212    0.5561   122.459
## age<56    0.8304     1.2043    0.3002     2.297
## 
## Concordance= 0.749  (se = 0.073 )
## Likelihood ratio test= 12.43  on 4 df,   p=0.01
## Wald test            = 17.7  on 4 df,   p=0.001
## Score (logrank) test = 43.75  on 4 df,   p=7e-09

ggforest(BRCA_TNM_cox.fit_dan,data = BRCA_clinical)

## Warning: Removed 4 rows containing missing values (geom_errorbar).

#ggsave(filename = "F:/Bioinfor_project/Breast/AS_research/AS/result/cli_multi_cox.pdf",  width=6,height = 4)
##结果保留为clinical_TNM_cox_result.xls

包装多因素cox回归函数multicox

multicox<-function(vars=c("T","N","M","age"),data=BRCA_clinical,forest=T){
  library(survival)
  require(survminer)
  y<-Surv(as.numeric(data[,"OS_time"]),
           as.numeric(data[,"OS_event"]))
   ## 构建公式
   FM<-as.formula(paste0("y~",paste(vars,collapse = "+")))
   cox.fit_multi <- coxph(FM,data = data)
   munivar_res<-summary(cox.fit_multi)#cox的结果
   pvalue<-round(munivar_res$coefficients[,"Pr(>|z|)"],3)#pvalue
   HR<-round(munivar_res$coefficients[,"exp(coef)"],2)# HR
   lower<-round(munivar_res$conf.int[,3],2)
   upper<-round(munivar_res$conf.int[,4],2)
   ## 保存结果到数据框中
  munivar_res<-data.frame(
    varname<-vars,
    HR<-as.numeric(HR),
    CI=paste(as.numeric(lower),"-",as.numeric(upper),sep = ""),
    pvalue<-as.numeric(pvalue)
  )
  colnames(munivar_res)<-c("varname","HR","95%CI","pvalue")
  ## 保存pdf到工作目录
  if (forest==T){
  ggforest(cox.fit_multi,data = data)
  ggsave(filename = "mult_cox.pdf")}
  munivar_res## 最后的函数运行结果
  
}
###
res<-multicox(vars=c("T","N","M","age"),data=BRCA_clinical,forest = F)
res

##   varname   HR       95%CI pvalue
## 1       T 2.29   0.62-8.43  0.213
## 2       N 5.12  1.51-17.34  0.009
## 3       M 8.25 0.56-122.46  0.125
## 4     age 0.83     0.3-2.3  0.720

#write.csv(res,file = "F:/Bioinfor_project/Breast/AS_research/AS/result/munivar_cox.csv")

batch cox-unicox and multicox

Chard Liu

2019年7月18日

载入数据

单因素cox批量计算函数unicox包装

包装多因素cox回归函数multicox