Skip to content
This repository has been archived by the owner on Oct 2, 2024. It is now read-only.

Update braf-fusion survival analyses with interaction models #58

Merged
merged 2 commits into from
Apr 15, 2024
Merged

Update braf-fusion survival analyses with interaction models #58

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6f9e3b6
add interaction models, change breakpoint ref
rjcorb Feb 13, 2024
910d3b8
add interaction models and plots
rjcorb Feb 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
84 changes: 61 additions & 23 deletions analyses/braf-fusions/03-survival-by-braf-breakpoints.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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)
```

Expand All @@ -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)

```
Expand All @@ -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)

```
Expand All @@ -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)

```
Expand All @@ -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)

```
Expand Down
170 changes: 102 additions & 68 deletions analyses/braf-fusions/03-survival-by-braf-breakpoints.html
View file
Open in desktop

Large diffs are not rendered by default.

Binary file modified analyses/braf-fusions/plots/breakpoint_group_ancestry_ct_enr_heatmap.pdf
View file
Open in desktop
Binary file not shown.
Binary file modified analyses/braf-fusions/plots/breakpoint_group_methyl_subtype_ct_enr_heatmap.pdf
View file
Open in desktop
Binary file not shown.
Binary file modified analyses/braf-fusions/plots/breakpoint_group_region_ct_enr_heatmap.pdf
View file
Open in desktop
Binary file not shown.
Binary file modified analyses/braf-fusions/plots/breakpoint_group_resection_ct_enr_heatmap.pdf
View file
Open in desktop
Binary file not shown.
Binary file renamed ...efs_breakpoint15_9_resection_ancestry.pdf → ...add_breakpoint15_9_resection_ancestry.pdf
View file
Open in desktop
Binary file not shown.
Binary file renamed ...efs_breakpoint16_9_resection_ancestry.pdf → ...add_breakpoint16_9_resection_ancestry.pdf
View file
Open in desktop
Binary file not shown.
Binary file renamed ...s_resection_ancestry_breakpoint_group.pdf → ...d_resection_ancestry_breakpoint_group.pdf
View file
Open in desktop
Binary file not shown.
Binary file renamed ...fs_resection_ancestry_breakpoint_type.pdf → ...dd_resection_ancestry_breakpoint_type.pdf
View file
Open in desktop
Binary file not shown.
Binary file renamed ...forest_efs_resection_breakpoint_group.pdf → ...st_efs_add_resection_breakpoint_group.pdf
View file
Open in desktop
Binary file not shown.
Binary file renamed ..._forest_efs_resection_breakpoint_type.pdf → ...est_efs_add_resection_breakpoint_type.pdf
View file
Open in desktop
Binary file not shown.
Binary file added analyses/braf-fusions/plots/pbta_ancestry_forest_efs_int_resection_breakpoint_group.pdf
View file
Open in desktop
Binary file not shown.
Binary file added analyses/braf-fusions/plots/pbta_ancestry_forest_efs_int_resection_breakpoint_type.pdf
View file
Open in desktop
Binary file not shown.
Binary file modified analyses/braf-fusions/plots/pbta_ancestry_km_efs_15_9_predicted_ancestry.pdf
View file
Open in desktop
Binary file not shown.
Binary file modified analyses/braf-fusions/plots/pbta_ancestry_km_efs_16_9_predicted_ancestry.pdf
View file
Open in desktop
Binary file not shown.
Binary file modified analyses/braf-fusions/plots/pbta_ancestry_km_efs_breakpoint_group.pdf
View file
Open in desktop
Binary file not shown.
Binary file modified analyses/braf-fusions/plots/pbta_ancestry_km_efs_breakpoint_type.pdf
View file
Open in desktop
Binary file not shown.
Binary file modified analyses/braf-fusions/plots/rare_breakpoints_by_ancestry.pdf
View file
Open in desktop
Binary file not shown.
Binary file renamed ...fs_15_9_resection_ancestry_breakpoint.RDS → ...dd_15_9_resection_ancestry_breakpoint.RDS
View file
Open in desktop
Binary file not shown.
Binary file renamed ...fs_16_9_resection_ancestry_breakpoint.RDS → ...dd_16_9_resection_ancestry_breakpoint.RDS
View file
Open in desktop
Binary file not shown.
Binary file renamed ...s_resection_ancestry_breakpoint_group.RDS → ...d_resection_ancestry_breakpoint_group.RDS
View file
Open in desktop
Binary file not shown.
Binary file renamed ...fs_resection_ancestry_breakpoint_type.RDS → ...dd_resection_ancestry_breakpoint_type.RDS
View file
Open in desktop
Binary file not shown.
Binary file renamed ...h_braf_efs_resection_breakpoint_group.RDS → ...af_efs_add_resection_breakpoint_group.RDS
View file
Open in desktop
Binary file not shown.
Binary file renamed ...ph_braf_efs_resection_breakpoint_type.RDS → ...raf_efs_add_resection_breakpoint_type.RDS
View file
Open in desktop
Binary file not shown.
Binary file added analyses/braf-fusions/results/coxph_braf_efs_int_resection_breakpoint_group.RDS
View file
Open in desktop
Binary file not shown.
Binary file added analyses/braf-fusions/results/coxph_braf_efs_int_resection_breakpoint_type.RDS
View file
Open in desktop
Binary file not shown.
111 changes: 57 additions & 54 deletions analyses/survival/util/survival_models.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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),
Expand All @@ -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),
Expand All @@ -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,
Expand All @@ -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()
Expand Down Expand Up @@ -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)
Expand All @@ -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,
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -639,5 +643,4 @@ plotForest <- function(model) {
print(forest_panels)
}




Loading