1. library(dplyr)
  2. library(tibble)
  3. library(ggplot2)
  4. library(survival)
  5. library(survminer)
  6. library(data.table)
  7. library(stringr)

load phenotype

  1. survival_data <- fread("../../Result/phenotype/common_survival_data.tsv")

clinical parameters

  1. get_clinical <- function(x){
  3.   dat_clin <- x
  4.   # pathologic_T 
  5.   dat_clin$Te1 <- str_extract(dat_clin$TNM, "T\\d[a-z]+")
  6.   dat_clin$Te2 <- str_extract(dat_clin$TNM, "T\\d")
  7.   dat_clin$Te1 <- with(dat_clin, ifelse(is.na(Te1), Te2, Te1))
  8.   # table(dat_clin$Te1)
  9.   # table(dat_clin$Te2)
  10.   dat_clin$Te2[dat_clin$Te2 == "T1"] <- "T1-T2"
  11.   dat_clin$Te2[dat_clin$Te2 == "T2"] <- "T1-T2"
  12.   dat_clin$Te2[dat_clin$Te2 == "T3"] <- "T1-T2"
  13.   dat_clin$Te2[dat_clin$Te2 == "T4"] <- "T3-T4"
  14.   Te2 <- table(dat_clin$Te2)
  16.   # pathologic_N
  17.   dat_clin$N <- str_extract(dat_clin$TNM, "N\\d")
  18.   dat_clin$N <- str_extract(dat_clin$N, "\\d")
  19.   # table(dat_clin$N)
  21.   dat_clin$N[dat_clin$N == "0"] <- "NO"
  22.   dat_clin$N[dat_clin$N == "1"] <- "Yes"
  23.   dat_clin$N[dat_clin$N == "2"] <- "Yes"
  24.   dat_clin$N[dat_clin$N == "3"] <- "Yes"
  25.   N <- table(dat_clin$N)
  27.   # pathologic_M
  28.   dat_clin$M <- str_extract(dat_clin$TNM, "M\\d")
  29.   dat_clin$M <- str_extract(dat_clin$M, "\\d")
  30.   # table(dat_clin$M)
  32.   dat_clin$M[dat_clin$M == "0"] <- "NO"
  33.   dat_clin$M[dat_clin$M == "1"] <- "Yes"
  34.   M <- table(dat_clin$M)
  36.   # Age
  37.   dat_clin$Age <- ifelse(dat_clin$Age > 60, "old", "young")
  38.   Age <- table(dat_clin$Age)
  40.   # Gender
  41.   dat_clin$Gender <- ifelse(dat_clin$Gender == "MALE", "Male", "Female")
  42.   Gender <- table(dat_clin$Gender)
  44.   # stage
  45.   dat_clin$Stage1 <- str_trim(str_extract(dat_clin$Stage, "\\s[H-Z]+"),
  46.                               side = c("both", "left", "right"))
  47.   # table(dat_clin$Stage1)
  48.   dat_clin$Stage1[dat_clin$Stage1 == "I"] <- "Stage I-II"
  49.   dat_clin$Stage1[dat_clin$Stage1 == "II"] <- "Stage I-II"
  50.   dat_clin$Stage1[dat_clin$Stage1 == "III"] <- "Stage III-IV"
  51.   dat_clin$Stage1[dat_clin$Stage1 == "IV"] <- "Stage III-IV"
  52.   Stage <- table(dat_clin$Stage1)
  54.   dat_cal <- data.frame(Clinicopathological=c(
  55.                           rep("Age", length(Age)),
  56.                           rep("Gender", length(Gender)),
  57.                           rep("Stage", length(Stage)),
  58.                           rep("T", length(Te2)),
  59.                           rep("N", length(N)),
  60.                           rep("M", length(M))), 
  61.                         Name = c(names(Age), names(Gender), 
  62.                                  names(Stage), names(Te2),
  63.                                  names(N), names(M)),
  64.                         Number = c(as.numeric(Age), as.numeric(Gender), 
  65.                                  as.numeric(Stage), as.numeric(Te2),
  66.                                  as.numeric(N), as.numeric(M))
  67.                         )
  68.   return(dat_cal)  
  69. }
  71. DT::datatable(get_clinical(survival_data))

curation

  1. get_plot <- function(dat=survival_data, tag="Gender"){
  3.   # dat=survival_data
  4.   # tag="Gender"
  6.   # pvalue
  7.   colnames(dat)[which(colnames(dat) == tag)] <- "group_info"
  8.   dat$group_info <- factor(dat$group_info)
  9.   factors <- unique(as.character(dat$group_info))
  11.   cox <- coxph(Surv(OS.Time, OS) ~ group_info, data = dat)
  12.   tmp <- summary(cox)
  13.   tmp.wald <- data.frame(t(tmp$waldtest)) %>%
  14.         setNames(c("Wald_test", "Wald_df", "Wald_pvlaue"))
  15.   tmp.lg <- data.frame(t(tmp$logtest)) %>%
  16.         setNames(c("lg_rank", "lg_rank_df", "lg_rank_pvlaue"))
  17.   tmp.total <- cbind(tmp.wald, tmp.lg) 
  19.   pvalue <- paste(paste0("Log-Rank P=", signif(tmp.lg$lg_rank_pvlaue, 3)), 
  20.                   paste0("Cox P=", signif(tmp.wald$Wald_pvlaue, 3)), sep = "\n")
  22.   # plot
  23.   fit <- survfit(Surv(OS.Time, OS) ~ group_info, data = dat)
  24.   info <- data.frame(time = fit$time,
  25.                   n.risk = fit$n.risk,
  26.                   n.event = fit$n.event,
  27.                   n.censor = fit$n.censor,
  28.                   surv = fit$surv,
  29.                   upper = fit$upper,
  30.                   lower = fit$lower)
  32.   pl <- ggsurvplot(fit,
  33.              data = dat,  
  34.              surv.median.line = "hv",
  35.              add.all = TRUE,
  36.              palette = "aaas",
  37.              risk.table = TRUE, 
  38.              xlab = "Follow up time(Years)", 
  39.              legend = c(0.8, 0.2), 
  40.              legend.title = "",
  41.              legend.labs = c("all", factors), 
  42.              break.x.by = 2,
  43.              font.legend = c(10, "italic"),
  44.              ggtheme = theme_bw())
  45.   pl$plot <- pl$plot + 
  46.              annotate("text", x=3, y=0.2, label=pvalue)
  48.   res <- list(info=info, pl=pl)
  49.   return(res)
  50. }

plot

  1. # gender 
  2. gender_plot <- get_plot(tag = "Gender")
  3. gender_plot$pl
  5. # stage 
  6. stage_plot <- get_plot(tag = "Stage")
  7. stage_plot$pl

version

  1. sessionInfo()
  1. R version 4.0.2 (2020-06-22)
  2. Platform: x86_64-w64-mingw32/x64 (64-bit)
  3. Running under: Windows 10 x64 (build 19042)
  5. Matrix products: default
  7. locale:
  8. [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
  9. [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
  10. [5] LC_TIME=English_United States.1252    
  11. system code page: 936
  13. attached base packages:
  14. [1] stats     graphics  grDevices utils     datasets  methods   base     
  16. other attached packages:
  17. [1] stringr_1.4.0     impute_1.62.0     data.table_1.13.2 survminer_0.4.8   ggpubr_0.4.0     
  18. [6] survival_3.2-7    ggplot2_3.3.2     tibble_3.0.3      dplyr_1.0.2      
  20. loaded via a namespace (and not attached):
  21.  [1] Rcpp_1.0.5        lattice_0.20-41   tidyr_1.1.2       zoo_1.8-8         digest_0.6.25    
  22.  [6] R6_2.5.0          cellranger_1.1.0  plyr_1.8.6        backports_1.1.10  evaluate_0.14    
  23. [11] pillar_1.4.6      rlang_0.4.7       curl_4.3          readxl_1.3.1      rstudioapi_0.11  
  24. [16] car_3.0-10        Matrix_1.2-18     DT_0.16           rmarkdown_2.5     splines_4.0.2    
  25. [21] foreign_0.8-80    htmlwidgets_1.5.2 munsell_0.5.0     broom_0.7.2       compiler_4.0.2   
  26. [26] xfun_0.19         pkgconfig_2.0.3   htmltools_0.5.0   tidyselect_1.1.0  gridExtra_2.3    
  27. [31] km.ci_0.5-2       rio_0.5.16        reshape_0.8.8     crayon_1.3.4      withr_2.3.0      
  28. [36] grid_4.0.2        jsonlite_1.7.1    xtable_1.8-4      gtable_0.3.0      lifecycle_0.2.0  
  29. [41] magrittr_1.5      KMsurv_0.1-5      scales_1.1.1      zip_2.1.1         stringi_1.5.3    
  30. [46] carData_3.0-4     ggsignif_0.6.0    ellipsis_0.3.1    survMisc_0.5.5    generics_0.1.0   
  31. [51] vctrs_0.3.4       openxlsx_4.2.2    tools_4.0.2       forcats_0.5.0     glue_1.4.2       
  32. [56] purrr_0.3.4       hms_0.5.3         crosstalk_1.1.0.1 rsconnect_0.8.16  abind_1.4-5      
  33. [61] yaml_2.2.1        colorspace_1.4-1  rstatix_0.6.0     knitr_1.30        haven_2.3.1

Reference

  1. 解决生存分析和临床参数相关分析