Merge pull request #58 from d3b-center/rjcorb/56-breakpoint-survival-…

…interaction

Update braf-fusion survival analyses with interaction models

analyses/braf-fusions/03-survival-by-braf-breakpoints.Rmd

@@ -47,12 +47,13 @@ braf_hist <- braf_hist %>%
     TRUE ~ predicted_ancestry
   )) %>%
   dplyr::mutate(breakpoint_group = fct_relevel(breakpoint_group,
-                                               c("16:9", "15:9", "16:11", "18:10", "rare")))
+                                               c("16:11", "15:9", "16:9", "18:10", "rare")))
 ```
 
 Run Kaplan-Meier analysis of event-free survival by breakpoint type (common vs. rare)
 
 ```{r}
+
 kap_efs_type <- survival_analysis(
   metadata  = braf_hist,
   ind_var = "breakpoint_type",
@@ -70,7 +71,7 @@ km_plot_type <- plotKM(model = kap_efs_type,
 km_plot_type
 
 ggsave(file.path(plot_dir, "pbta_ancestry_km_efs_breakpoint_type.pdf"),
-       width = 10, height = 5)
+       width = 7.5, height = 5)
 ```
 
 Run Kaplan-Meier analysis of event-free survival by breakpoint group
@@ -93,7 +94,7 @@ km_plot_group <- plotKM(model = kap_efs_group,
 km_plot_group
 
 ggsave(file.path(plot_dir, "pbta_ancestry_km_efs_breakpoint_group.pdf"),
-       width = 9, height = 5)
+       width = 7.5, height = 5)
 ```
 
 Run Kaplan-Meier analysis of event-free survival by predicted ancestry among 16:9 breakpoint pts
@@ -116,7 +117,7 @@ km_plot_group <- plotKM(model = kap_efs_group,
 km_plot_group
 
 ggsave(file.path(plot_dir, "pbta_ancestry_km_efs_16_9_predicted_ancestry.pdf"),
-       width = 9, height = 5)
+       width = 7.5, height = 5)
 ```
 
 Run Kaplan-Meier analysis of event-free survival by predicted ancestry among 15:9 breakpoint pts (combine asian ancestries here)
@@ -139,23 +140,42 @@ km_plot_group <- plotKM(model = kap_efs_group,
 km_plot_group
 
 ggsave(file.path(plot_dir, "pbta_ancestry_km_efs_15_9_predicted_ancestry.pdf"),
-       width = 9, height = 5)
+       width = 7.5, height = 5)
 
 ```
 
 Run and plot cox proportional hazards model of EFS including tumor resection and breakpoint type as covariates
 
+Interaction model:
+
+```{r}
+braf_model_efs <- fit_save_model(braf_hist[braf_hist$extent_of_tumor_resection != "Not Reported",],
+                                    terms = "extent_of_tumor_resection*breakpoint_type",
+               file.path(results_dir, "coxph_braf_efs_int_resection_breakpoint_type.RDS"),
+               "multivariate",
+               years_col = "EFS_years", status_col = "EFS_status"
+              )
+
+plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_int_resection_breakpoint_type.RDS")))
+
+ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_int_resection_breakpoint_type.pdf"),
+       width = 8, height = 3)
+```
+
+
+Additive model: 
+
 ```{r}
 braf_model_efs <- fit_save_model(braf_hist[braf_hist$extent_of_tumor_resection != "Not Reported",],
                                     terms = "extent_of_tumor_resection+breakpoint_type",
-               file.path(results_dir, "coxph_braf_efs_resection_breakpoint_type.RDS"),
+               file.path(results_dir, "coxph_braf_efs_add_resection_breakpoint_type.RDS"),
                "multivariate",
                years_col = "EFS_years", status_col = "EFS_status"
               )
 
-plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_resection_breakpoint_type.RDS")))
+plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_add_resection_breakpoint_type.RDS")))
 
-ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_resection_breakpoint_type.pdf"),
+ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_add_resection_breakpoint_type.pdf"),
        width = 8, height = 3)
 ```
 
@@ -170,31 +190,49 @@ braf_hist <- braf_hist %>%
 
 braf_model_efs <- fit_save_model(braf_hist[braf_hist$extent_of_tumor_resection != "Not Reported",],
                                     terms = "extent_of_tumor_resection+predicted_ancestry+breakpoint_type",
-               file.path(results_dir, "coxph_braf_efs_resection_ancestry_breakpoint_type.RDS"),
+               file.path(results_dir, "coxph_braf_efs_add_resection_ancestry_breakpoint_type.RDS"),
                "multivariate",
                years_col = "EFS_years", status_col = "EFS_status"
               )
 
-plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_resection_ancestry_breakpoint_type.RDS")))
+plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_add_resection_ancestry_breakpoint_type.RDS")))
 
-ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_resection_ancestry_breakpoint_type.pdf"),
+ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_add_resection_ancestry_breakpoint_type.pdf"),
        width = 8, height = 4)
 
 ```
 
 Run and plot cox proportional hazards model of EFS including tumor resection and breakpoint group as covariates
 
+Interaction model:
+
+```{r}
+braf_model_efs <- fit_save_model(braf_hist[braf_hist$extent_of_tumor_resection != "Not Reported",],
+                                    terms = "extent_of_tumor_resection*breakpoint_group",
+               file.path(results_dir, "coxph_braf_efs_int_resection_breakpoint_group.RDS"),
+               "multivariate",
+               years_col = "EFS_years", status_col = "EFS_status"
+              )
+
+plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_int_resection_breakpoint_group.RDS")))
+
+ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_int_resection_breakpoint_group.pdf"),
+       width = 8, height = 3)
+
+```
+Additive model:
+
 ```{r}
 braf_model_efs <- fit_save_model(braf_hist[braf_hist$extent_of_tumor_resection != "Not Reported",],
                                     terms = "extent_of_tumor_resection+breakpoint_group",
-               file.path(results_dir, "coxph_braf_efs_resection_breakpoint_group.RDS"),
+               file.path(results_dir, "coxph_braf_efs_add_resection_breakpoint_group.RDS"),
                "multivariate",
                years_col = "EFS_years", status_col = "EFS_status"
               )
 
-plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_resection_breakpoint_group.RDS")))
+plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_add_resection_breakpoint_group.RDS")))
 
-ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_resection_breakpoint_group.pdf"),
+ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_add_resection_breakpoint_group.pdf"),
        width = 8, height = 3)
 
 ```
@@ -204,14 +242,14 @@ Re-run with predicted ancestry as covariate
 ```{r}
 braf_model_efs <- fit_save_model(braf_hist[braf_hist$extent_of_tumor_resection != "Not Reported",],
                                     terms = "extent_of_tumor_resection+predicted_ancestry+breakpoint_group",
-               file.path(results_dir, "coxph_braf_efs_resection_ancestry_breakpoint_group.RDS"),
+               file.path(results_dir, "coxph_braf_efs_add_resection_ancestry_breakpoint_group.RDS"),
                "multivariate",
                years_col = "EFS_years", status_col = "EFS_status"
               )
 
-plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_resection_ancestry_breakpoint_group.RDS")))
+plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_add_resection_ancestry_breakpoint_group.RDS")))
 
-ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_resection_ancestry_breakpoint_group.pdf"),
+ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_add_resection_ancestry_breakpoint_group.pdf"),
        width = 8, height = 4)
 
 ```
@@ -221,14 +259,14 @@ Run and plot cox proportional hazards model of EFS in 16:09 breakpoint patients,
 ```{r}
 braf_model_efs <- fit_save_model(braf_hist[braf_hist$extent_of_tumor_resection != "Not Reported" & braf_hist$breakpoint_group == "16:9",],
                                     terms = "extent_of_tumor_resection+predicted_ancestry",
-               file.path(results_dir, "coxph_braf_efs_16_9_resection_ancestry_breakpoint.RDS"),
+               file.path(results_dir, "coxph_braf_efs_add_16_9_resection_ancestry_breakpoint.RDS"),
                "multivariate",
                years_col = "EFS_years", status_col = "EFS_status"
               )
 
-plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_16_9_resection_ancestry_breakpoint.RDS")))
+plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_add_16_9_resection_ancestry_breakpoint.RDS")))
 
-ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_breakpoint16_9_resection_ancestry.pdf"),
+ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_add_breakpoint16_9_resection_ancestry.pdf"),
        width = 8, height = 3)
 
 ```
@@ -238,14 +276,14 @@ Run and plot cox proportional hazards model of EFS in 15:09 breakpoint patients,
 ```{r}
 braf_model_efs <- fit_save_model(braf_hist[braf_hist$extent_of_tumor_resection != "Not Reported" & braf_hist$breakpoint_group == "15:9",],
                                     terms = "extent_of_tumor_resection+AS_ancestry",
-               file.path(results_dir, "coxph_braf_efs_15_9_resection_ancestry_breakpoint.RDS"),
+               file.path(results_dir, "coxph_braf_efs_add_15_9_resection_ancestry_breakpoint.RDS"),
                "multivariate",
                years_col = "EFS_years", status_col = "EFS_status"
               )
 
-plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_15_9_resection_ancestry_breakpoint.RDS")))
+plotForest(readRDS(file.path(results_dir, "coxph_braf_efs_add_15_9_resection_ancestry_breakpoint.RDS")))
 
-ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_breakpoint15_9_resection_ancestry.pdf"),
+ggsave(file.path(plot_dir, "pbta_ancestry_forest_efs_add_breakpoint15_9_resection_ancestry.pdf"),
        width = 8, height = 3)
 
 ```

analyses/survival/util/survival_models.R

@@ -96,22 +96,22 @@ survival_analysis <- function(metadata,
   # Pull out this data
   status <- metadata %>%
     dplyr::pull({{status_col}})
-
+  
   # Code status variable as 0 (no event) or 1 (event)
   metadata <- metadata %>%
     dplyr::mutate(
       !!status_col := case_when(
         status %in% c("LIVING", "NO EVENT") ~ FALSE,
         status %in% c("DECEASED", "EVENT") ~ TRUE
       ))
-      
-      
-    #   !!status_col := case_when(
-    #     status_col == "OS_status" ~ as.numeric(
-    #       factor(status, levels = c("LIVING", "DECEASED"))),
-    #     status_col == "EFS_status" ~ as.numeric(
-    #       factor(status, levels = c("NO EVENT", "EVENT")))
-    # ))
+
+
+  #   !!status_col := case_when(
+  #     status_col == "OS_status" ~ as.numeric(
+  #       factor(status, levels = c("LIVING", "DECEASED"))),
+  #     status_col == "EFS_status" ~ as.numeric(
+  #       factor(status, levels = c("NO EVENT", "EVENT")))
+  # ))
 
   ############################ Set up the ind data #############################
   # If other ind_data has been supplied, attempt to join it to metadata
@@ -258,38 +258,38 @@ fit_save_model <- function(df,
       metadata_sample_col = "Kids_First_Biospecimen_ID",
       os_days_col = years_col # we want to use years, not days
     )
-      } else if (model_type == "multivariate") {
-        
-        if (years_col == "EFS_years"){
-          
-          # Recode OS_Status (for univariate, `survival_analysis()` does the recoding)
-          df <- df %>%
-            mutate(EFS_status = ifelse(EFS_status == "NO EVENT", 0, 1))
-          
-          # Fit model
-          fitted_multi <- survival::coxph(
-            formula(
-              paste0("survival::Surv(time = EFS_years, event = EFS_status) ~ ", terms)
-            ),
-            data = df
-          )
-          
-        } else if (years_col == "OS_years"){
-          
-          # Recode OS_Status (for univariate, `survival_analysis()` does the recoding)
-          df <- df %>%
-            mutate(OS_status = ifelse(OS_status == "LIVING", 0, 1))
-          
-          # Fit model
-          fitted_multi <- survival::coxph(
-            formula(
-              paste0("survival::Surv(time = OS_years, event = OS_status) ~ ", terms)
-            ),
-            data = df
-          )
-          
-        }
-
+  } else if (model_type == "multivariate") {
+
+    if (years_col == "EFS_years"){
+
+      # Recode OS_Status (for univariate, `survival_analysis()` does the recoding)
+      df <- df %>%
+        mutate(EFS_status = ifelse(EFS_status == "NO EVENT", 0, 1))
+
+      # Fit model
+      fitted_multi <- survival::coxph(
+        formula(
+          paste0("survival::Surv(time = EFS_years, event = EFS_status) ~ ", terms)
+        ),
+        data = df
+      )
+
+    } else if (years_col == "OS_years"){
+
+      # Recode OS_Status (for univariate, `survival_analysis()` does the recoding)
+      df <- df %>%
+        mutate(OS_status = ifelse(OS_status == "LIVING", 0, 1))
+
+      # Fit model
+      fitted_multi <- survival::coxph(
+        formula(
+          paste0("survival::Surv(time = OS_years, event = OS_status) ~ ", terms)
+        ),
+        data = df
+      )
+
+    }
+    
     # Set up list object to match parts of `survival_analysis()` output we need
     fit <- list(
       model = fitted_multi,
@@ -326,7 +326,7 @@ plotKM <- function(model,
       event_type <- "OS"
 
       diff_obj <- survdiff(survival::Surv(OS_days, OS_status) ~ term,  
-                              model$original_data)
+                           model$original_data)
       diff_pvalue <- 1 - pchisq(diff_obj$chisq, length(diff_obj$n) - 1)
       diff_pvalue_formatted <- format(
         signif(diff_pvalue, 2),
@@ -350,7 +350,7 @@ plotKM <- function(model,
       event_type <- "EFS"
 
       diff_obj <- survdiff(survival::Surv(EFS_days, EFS_status) ~ term,  
-                               model$original_data)
+                           model$original_data)
       diff_pvalue <- 1 - pchisq(diff_obj$chisq, length(diff_obj$n) - 1)
       diff_pvalue_formatted <- format(
         signif(diff_pvalue, 2),
@@ -373,10 +373,10 @@ plotKM <- function(model,
     colors <- colorblindr::palette_OkabeIto[1:(length(levels)+1)]
     colors <- colors[-4]
     lines <- c(rep("solid", length(levels)), 
-              rep("dashed", length(levels)))
+               rep("dashed", length(levels)))
     labels <- glue::glue("{event_type}:{levels}")
 
-      
+
     km_plot <- survminer::ggsurvplot(fit = model$model, 
                                      data = model$original_data,
                                      palette = colors,
@@ -395,7 +395,7 @@ plotKM <- function(model,
 
     km_plot_graph <- km_plot$plot + 
       ggplot2::annotate("text", 
-                        200, 0.15, 
+                        3500, 0.9, 
                         label = pvalue_label) +
       theme(legend.text = element_text(size = 16, color = "black")) +
       cowplot::background_grid()
@@ -452,8 +452,8 @@ plotKM <- function(model,
       scientific = FALSE)
 
     os_pvalue_label <- ifelse(diff_os_pvalue_formatted < 0.001, 
-                           "OS P < 0.001",
-                           paste0("OS P = ", diff_os_pvalue_formatted))
+                              "OS P < 0.001",
+                              paste0("OS P = ", diff_os_pvalue_formatted))
 
     diff_efs_obj <- survdiff(survival::Surv(EFS_days, EFS_status) ~ variable_efs,  
                              data_efs)
@@ -463,8 +463,8 @@ plotKM <- function(model,
       scientific = FALSE)
 
     efs_pvalue_label <- ifelse(diff_efs_pvalue_formatted < 0.001, 
-                              "EFS P < 0.001",
-                              paste0("EFS P = ", diff_efs_pvalue_formatted))
+                               "EFS P < 0.001",
+                               paste0("EFS P = ", diff_efs_pvalue_formatted))
 
     km_plot <- survminer::ggsurvplot(fit = fit, 
                                      data = data_efs,
@@ -552,10 +552,14 @@ plotForest <- function(model) {
                     levels = term_order,
                     labels = term_labels)
     ) %>%
-    filter(estimate > 1e-4 & estimate < 1500)
+    filter(estimate > 1e-4 & estimate < 1500) %>%
+    arrange(term) %>%
+    dplyr::mutate(term_display = str_replace(str_replace(term, paste(names(model$xlevels), collapse = "|"), ""), paste(names(model$xlevels), collapse = "|"), "")) %>%
+    dplyr::mutate(term_display = fct_relevel(term_display, unique(term_display)))
+
 
   forest_plot <- ggplot(survival_df) +
-    aes(x = estimate, y = term, fill = significant
+    aes(x = estimate, y = factor(term_display), fill = significant
     ) + 
     # add CI first so line doesn't cover open point
     geom_errorbarh(
@@ -639,5 +643,4 @@ plotForest <- function(model) {
   print(forest_panels)
 }
 
-
-
+