diff --git a/ven_hma_samples/Figure 0/introduce_data.Rmd b/ven_hma_samples/Figure 0/introduce_data.Rmd new file mode 100644 index 0000000..8538c6f --- /dev/null +++ b/ven_hma_samples/Figure 0/introduce_data.Rmd @@ -0,0 +1,279 @@ + + +```{r include=FALSE} +source("../../synapseUtil.R") +source("../scripts/loading_data.R") +source("../scripts/plot_helper.R") +source("../scripts/compute_lsc17.R") +source("../scripts/predict_helper.R") +source("../scripts/s2n_model.R") + +samples_to_exclude = c("PTRC_Exp28_Prot_33", "PTRC_Exp28_Prot_15", "PTRC_Exp28_Prot_03") +m_exp28 = m_exp28[, setdiff(sampleNames(m_exp28), samples_to_exclude)] + +features_50perc = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) < 19] +complete_features = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) == 0] + +plot_folder = "./many_plots/" +if (!dir.exists(plot_folder)){ + dir.create(plot_folder) +} + +``` + +## Load msigdb + +```{r} + +library(msigdbr) +t2g_hallmark <- msigdbr(species = "Homo sapiens", category = "H") %>% + dplyr::select(gs_name, gene_symbol) + +t2g_reactome <- msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:REACTOME") %>% + dplyr::select(gs_name, gene_symbol) %>% + mutate(gs_name = gsub("_", " ", gs_name)) + +t2g_kegg <- msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:KEGG_LEGACY") %>% + dplyr::select(gs_name, gene_symbol) %>% + mutate(gs_name = gsub("_", " ", gs_name)) + +t2g_gobp <- msigdbr(species = "Homo sapiens", category = "C5", subcategory = "BP") %>% + dplyr::select(gs_name, gene_symbol) %>% + mutate(gs_name = gsub("_", " ", gs_name)) + + +``` + + +SSGSEA + + +```{r} +library(GSVA) + +pathway_df = t2g_kegg %>% unique() %>% + group_by(gs_name) %>% + summarise(genes = list(gene_symbol)) +pathway_list = pathway_df$genes +names(pathway_list) = pathway_df$gs_name + +m_exp28_complete = m_exp28[rowSums(is.na(exprs(m_exp28))) <= 0, ] +my_params <- ssgseaParam(expr = exprs(m_exp28_complete), + geneSets = pathway_list) +results_kegg <- gsva(my_params) + + +pathway_df = t2g_hallmark %>% unique() %>% + group_by(gs_name) %>% + summarise(genes = list(gene_symbol)) +pathway_list = pathway_df$genes +names(pathway_list) = pathway_df$gs_name + +m_exp28_complete = m_exp28[rowSums(is.na(exprs(m_exp28))) <= 0, ] +my_params <- ssgseaParam(expr = exprs(m_exp28_complete), + geneSets = pathway_list) +results_hallmark <- gsva(my_params) + + +pathway_df = t2g_reactome %>% unique() %>% + group_by(gs_name) %>% + summarise(genes = list(gene_symbol)) +pathway_list = pathway_df$genes +names(pathway_list) = pathway_df$gs_name + +m_exp28_complete = m_exp28[rowSums(is.na(exprs(m_exp28))) <= 0, ] +my_params <- ssgseaParam(expr = exprs(m_exp28_complete), + geneSets = pathway_list) +results_reactome <- gsva(my_params) + + +pathway_df = t2g_gobp %>% unique() %>% + group_by(gs_name) %>% + summarise(genes = list(gene_symbol)) +pathway_list = pathway_df$genes +names(pathway_list) = pathway_df$gs_name + +m_exp28_complete = m_exp28[rowSums(is.na(exprs(m_exp28))) <= 0, ] +my_params <- ssgseaParam(expr = exprs(m_exp28_complete), + geneSets = pathway_list) +results_gobp <- gsva(my_params) + + +m_exp28_ssgsea = MSnSet(exprs = rbind(results_kegg, results_reactome, results_hallmark, results_gobp), + pData = pData(m_exp28)) +xx = plot_pca(m_exp28, "subcohort", output_type = "full") +yy = plot_pca(m_exp28_ssgsea, "subcohort", output_type = "full") +zz = plot_pca(m_4pat[, m_4pat$Treatment %in% c("DV", "UT")], "Treatment", output_type = "full") + +``` + + +PCA of SSGSEA + +```{r} +treatment_colors = c("DV" = "#00798c", "UT" = "forestgreen") +subcohort_colors_ = subcohort_colors +subcohort_colors_[[1]] = "#bd7b13" +df.u = xx$data.u +df.u$subcohort = pData(m_exp28)[, "subcohort"] +df.v = xx$data.v + +p <- ggplot(data = df.u, mapping = aes(x = df.u[, 1], y = df.u[, 2], + color = subcohort)) + geom_point(size = 2) + + geom_hline(yintercept = 0, lty = "longdash", color = "darkgrey") + + geom_vline(xintercept = 0, lty = "longdash", color = "darkgrey") + + labs(x = xx$axis_labs[1], y = xx$axis_labs[2]) + + theme_bw() + theme(aspect.ratio = 1, text = element_text(size = 15)) + + stat_ellipse(mapping = aes(fill = subcohort, color = NULL), + geom = "polygon", type = "norm", + level = 0.5, alpha = 0.1, show.legend = FALSE) + + ggtitle("PCA of (5309 features)") + scale_fill_manual(values = subcohort_colors) + + scale_color_manual(values = subcohort_colors) +p +ggsave(plot = p, filename = paste0(plot_folder, "pca_m_exp28.pdf")) + +df.u = yy$data.u +df.u$subcohort = pData(m_exp28)[, "subcohort"] +df.v = yy$data.v + +p <- ggplot(data = df.u, mapping = aes(x = df.u[, 1], y = df.u[, 2], + color = subcohort)) + geom_point(size = 2) + + geom_hline(yintercept = 0, lty = "longdash", color = "darkgrey") + + geom_vline(xintercept = 0, lty = "longdash", color = "darkgrey") + + labs(x = yy$axis_labs[1], y = yy$axis_labs[2]) + + theme_bw() + theme(aspect.ratio = 1, text = element_text(size = 15)) + + stat_ellipse(mapping = aes(fill = subcohort, color = NULL), + geom = "polygon", type = "norm", + level = 0.5, alpha = 0.07, show.legend = FALSE) + + ggtitle("PCA of (8844 pathways)") + scale_fill_manual(values = subcohort_colors) + + scale_color_manual(values = subcohort_colors) +p +ggsave(plot = p, filename = paste0(plot_folder, "pca_m_exp28_sigs.pdf")) + + +df.u = zz$data.u +df.u$Treatment = pData(m_4pat[, m_4pat$Treatment %in% c("DV", "UT")])[, "Treatment"] +df.v = zz$data.v + +p <- ggplot(data = df.u, mapping = aes(x = df.u[, 1], y = df.u[, 2], + color = Treatment)) + geom_point(size = 2) + + geom_hline(yintercept = 0, lty = "longdash", color = "darkgrey") + + geom_vline(xintercept = 0, lty = "longdash", color = "darkgrey") + + labs(x = zz$axis_labs[1], y = zz$axis_labs[2]) + + theme_bw() + theme(aspect.ratio = 1, text = element_text(size = 15)) + + stat_ellipse(mapping = aes(fill = Treatment, color = NULL), + geom = "polygon", type = "norm", + level = 0.5, alpha = 0.1, show.legend = FALSE) + + ggtitle("PCA of 3982 features") + scale_fill_manual(values = treatment_colors) + + scale_color_manual(values = treatment_colors) +p +ggsave(plot = p, filename = paste0(plot_folder, "pca_m_4pat.pdf")) + + +``` + + +Diffexp of SSGSEA + +```{r} +diffexp_results = diffexp_helper(m_exp28, "subcohort", c("Refractory-Response_no_relapse", "Refractory-Relapse")) +diffexp_results_sigs = diffexp_helper(m_exp28_ssgsea, "subcohort", c("Refractory-Response_no_relapse", "Refractory-Relapse")) + +``` + + +```{r} +logFC_cutoff = 3 +plot_df = diffexp_results %>% + mutate(pvalue_ = -log10(t_test_pval), + alpha = case_when(t_test_pval <= 0.001 ~ 1, + logFC > logFC_cutoff | logFC < -logFC_cutoff ~ 1, + TRUE ~ 0.01), + color = case_when(t_test_pval < 0.001 & logFC > 0 ~ "logFC > 0", + t_test_pval < 0.001 & logFC < 0 ~ "logFC < 0", + logFC > logFC_cutoff ~ "logFC > 0", + logFC < -logFC_cutoff ~ "logFC < 0", + TRUE ~ ""), + size = case_when(t_test_pval < 0.001 ~ "A", + abs(logFC) > logFC_cutoff ~ "A", + TRUE ~ "B")) %>% select(-feature) %>% unique() +plot_df$feature = rownames(plot_df) +p = ggplot(plot_df, aes(x = logFC, y = pvalue_, alpha = alpha, color = color, size = size)) + + geom_point() + facet_wrap(~contrast) + scale_color_manual(values = c("grey20", "steelblue4", "firebrick4")) + + geom_hline(aes(yintercept = -log10(0.001)), linetype = "dashed", color = "forestgreen", linewidth = 0.7) + + guides(alpha = "none", size = "none", color = "none") + ylab("pvalue") + xlim(-5.4, 5.4) + + scale_y_continuous(breaks = c(0,2,3,4,6,8), labels = c("1", "0.01", "0.001", "1e-4", "1e-6", "1e-8")) + theme_bw() + + scale_size_manual(values = c("A" = 7, "B" = 3)) + + theme(text = element_text(size = 27), strip.background = element_rect(fill = NA, color = "black")) + + ggrepel::geom_label_repel(data = plot_df %>% filter(t_test_pval < 0.001 | abs(logFC) > logFC_cutoff), + aes(label = feature), min.segment.length = 1, size = 5) +p +plot_path = paste0(plot_folder, "volcano_plots_m_exp28.png") +ggsave(plot = p, filename = plot_path, width = 15, height = 7) +plot_path = paste0(plot_folder, "volcano_plots_m_exp28.pdf") +ggsave(plot = p, filename = plot_path, width = 15, height = 7) + + +logFC_cutoff = 0.2 +plot_df = diffexp_results_sigs %>% + mutate(pvalue_ = -log10(t_test_pval), + alpha = case_when(t_test_pval <= 0.001 ~ 1, + logFC > logFC_cutoff | logFC < -logFC_cutoff ~ 1, + TRUE ~ 0.01), + color = case_when(t_test_pval < 0.001 & logFC > 0 ~ "logFC > 0", + t_test_pval < 0.001 & logFC < 0 ~ "logFC < 0", + logFC > logFC_cutoff ~ "logFC > 0", + logFC < -logFC_cutoff ~ "logFC < 0", + TRUE ~ ""), + size = case_when(t_test_pval < 0.001 ~ "A", + abs(logFC) > logFC_cutoff ~ "A", + TRUE ~ "B")) %>% select(-feature) %>% unique() +plot_df <- plot_df %>% + mutate(feature = rownames(.) %>% as.character()) %>% + mutate(feature = stringr::str_wrap(feature, width = 30)) +p = ggplot(plot_df, aes(x = logFC, y = pvalue_, alpha = alpha, color = color, size = size)) + + geom_point() + facet_wrap(~contrast) + scale_color_manual(values = c("grey20", "steelblue4", "firebrick4")) + + geom_hline(aes(yintercept = -log10(0.001)), linetype = "dashed", color = "forestgreen", linewidth = 0.7) + + guides(alpha = "none", size = "none", color = "none") + ylab("pvalue") + xlim(-0.35, 0.35) + + scale_y_continuous(breaks = c(0,2,3,4,6,8), labels = c("1", "0.01", "0.001", "1e-4", "1e-6", "1e-8")) + theme_bw() + + scale_size_manual(values = c("A" = 7, "B" = 3)) + + theme(text = element_text(size = 27), strip.background = element_rect(fill = NA, color = "black")) + + ggrepel::geom_label_repel(data = plot_df %>% filter(t_test_pval < 0.001 | abs(logFC) > logFC_cutoff), + aes(label = feature), point.padding = unit(0.5, "lines"), nudge_y = 0.7, + min.segment.length = 1, size = 5) +p +plot_path = paste0(plot_folder, "volcano_plots_m_exp28_sigs.png") +ggsave(plot = p, filename = plot_path, width = 15, height = 7) +plot_path = paste0(plot_folder, "volcano_plots_m_exp28_sigs.pdf") +ggsave(plot = p, filename = plot_path, width = 15, height = 7) + + +``` + + + +```{r} + +plot_df = diffexp_results_sigs %>% + mutate(pvalue_ = -log10(t_test_pval), + alpha = case_when(t_test_pval <= 0.05 ~ 1, + TRUE ~ 0.01), + color = case_when(t_test_pval < 0.05 & logFC > 0 ~ "pvalue < 0.05, logFC > 0", + t_test_pval < 0.05 & logFC < 0 ~ "pvalue < 0.05, logFC < 0", + TRUE ~ ""), + size = case_when(t_test_pval < 0.05 ~ "A", + TRUE ~ "B")) +p = ggplot(plot_df, aes(x = logFC, y = pvalue_, alpha = alpha, color = color, size = size)) + + geom_point() + facet_wrap(~contrast) + scale_color_manual(values = c("grey20", "steelblue4", "firebrick4")) + + geom_hline(aes(yintercept = -log10(0.05)), linetype = "dashed", color = "forestgreen", linewidth = 0.7) + + guides(alpha = "none", size = "none") + ylab("pvalue") + theme(text = element_text(size = 15)) + xlim(-0.27, 0.27) + + scale_y_continuous(breaks = c(0,2,4,6,8), labels = c("1", "0.01", "1e-4", "1e-6", "1e-8")) + theme_bw() + + scale_size_manual(values = c("A" = 3, "B" = 1)) + theme(text = element_text(size = 18)) +p +# plot_path = paste0(folder_path, "volcano_plots.png") +# ggsave(plot = p, filename = plot_path, width = 10, height = 7) + +``` + + + diff --git a/ven_hma_samples/Figure 1/cell_type_sigs.Rmd b/ven_hma_samples/Figure 1/cell_type_sigs.Rmd index 601a14f..37bfcb4 100644 --- a/ven_hma_samples/Figure 1/cell_type_sigs.Rmd +++ b/ven_hma_samples/Figure 1/cell_type_sigs.Rmd @@ -5,13 +5,12 @@ Loading data ```{r} source("../../synapseUtil.R") +source("../scripts/loading_data.R") source("../scripts/plot_helper.R") source("../scripts/compute_lsc17.R") source("../scripts/predict_helper.R") source("../scripts/s2n_model.R") -source("../scripts/loading_data.R") - samples_to_exclude = c("PTRC_Exp28_Prot_33", "PTRC_Exp28_Prot_15", "PTRC_Exp28_Prot_03") m_exp28 = m_exp28[, setdiff(sampleNames(m_exp28), samples_to_exclude)] diff --git a/ven_hma_samples/Figure 2/exvivo_signature.Rmd b/ven_hma_samples/Figure 2/exvivo_signature.Rmd index 20a3930..1f4156f 100644 --- a/ven_hma_samples/Figure 2/exvivo_signature.Rmd +++ b/ven_hma_samples/Figure 2/exvivo_signature.Rmd @@ -1,14 +1,14 @@ -```{r} +```{r include=FALSE} source("../../synapseUtil.R") +source("../scripts/loading_data.R") source("../scripts/plot_helper.R") source("../scripts/compute_lsc17.R") source("../scripts/predict_helper.R") source("../scripts/s2n_model.R") -source("../scripts/loading_data.R") samples_to_exclude = c("PTRC_Exp28_Prot_33", "PTRC_Exp28_Prot_15", "PTRC_Exp28_Prot_03") m_exp28 = m_exp28[, setdiff(sampleNames(m_exp28), samples_to_exclude)] @@ -17,20 +17,204 @@ features_50perc = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) < 19] complete_features = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) == 0] plot_folder = "./many_plots/" -if (!file.exists(plot_folder)) - dir.create(plot_folder) +if (!dir.exists(plot_folder)){ + dir.create(plot_folder) +} + +complete_DV_UT_features <- complete.cases(exprs(m_4pat[, m_4pat$Treatment %in% c("UT", "DV")])) +complete_DV_UT_features_ = intersect(featureNames(m_4pat)[complete_DV_UT_features], complete_features) # diffexp_results = diffexp_helper(m_exp28, "subcohort", c("Refractory-Response_no_relapse", "Refractory-Relapse")) ``` +PCA of 4 patient ex vivo data + +```{r} +pat_colors = subcohort_colors +pat_colors = c("DV Resistant" = "#00798c", "Untreated" = "#e8991b") +m_ = m_4pat[, m_4pat$Treatment %in% c("UT", "DV")] +yy = plot_pca(m_, phenotype = "Treatment", output_type = "full") +df.u = yy$data.u +df.u$Treatment = pData(m_)[, "Treatment"] +df.u = df.u %>% + mutate(Group = case_when(Treatment == "DV" ~ "DV Resistant", + TRUE ~ "Untreated")) +df.v = yy$data.v + +p <- ggplot(data = df.u, mapping = aes(x = df.u[, 1], y = df.u[, 2], + color = Group)) + geom_point(size = 2) + + geom_hline(yintercept = 0, lty = "longdash", color = "darkgrey") + + geom_vline(xintercept = 0, lty = "longdash", color = "darkgrey") + + labs(x = yy$axis_labs[1], y = yy$axis_labs[2]) + + theme_bw() + theme(aspect.ratio = 1, text = element_text(size = 15)) + + stat_ellipse(mapping = aes(fill = Group, color = NULL), + geom = "polygon", type = "norm", + level = 0.5, alpha = 0.07, show.legend = FALSE) + + ggtitle("PCA of 3932 features") + scale_fill_manual(values = pat_colors) + + scale_color_manual(values = pat_colors) +p +ggsave(plot = p, filename = paste0(plot_folder, "pca_m_4pat.pdf")) + +``` + + +Volcano plot showing the markers within the exvivo data + +```{r} +diffexp_results_4pat = diffexp_helper(m_4pat[complete_DV_UT_features_, ], "Treatment", c("DV-UT")) + +logFC_cutoff = 4 +diffexp_results = diffexp_results_4pat +plot_df = diffexp_results %>% + mutate(pvalue_ = -log10(t_test_pval), + alpha = case_when(t_test_pval <= 0.001 ~ 1, + logFC > logFC_cutoff | logFC < -logFC_cutoff ~ 1, + TRUE ~ 0.01), + color = case_when(t_test_pval < 0.001 & logFC > 0 ~ "logFC > 0", + t_test_pval < 0.001 & logFC < 0 ~ "logFC < 0", + logFC > logFC_cutoff ~ "logFC > 0", + logFC < -logFC_cutoff ~ "logFC < 0", + TRUE ~ ""), + size = case_when(t_test_pval < 0.001 ~ "A", + abs(logFC) > logFC_cutoff ~ "A", + TRUE ~ "B")) %>% select(-feature) %>% unique() +plot_df$feature = rownames(plot_df) +p = ggplot(plot_df, aes(x = logFC, y = pvalue_, alpha = alpha, color = color, size = size)) + + geom_point() + scale_color_manual(values = c("grey20", "steelblue4", "firebrick4")) + + geom_hline(aes(yintercept = -log10(0.001)), linetype = "dashed", color = "forestgreen", linewidth = 0.7) + + guides(alpha = "none", size = "none", color = "none") + ylab("pvalue") + xlim(-3.3, 3.3) + + scale_y_continuous(breaks = c(0,2,3,4,6,8), labels = c("1", "0.01", "0.001", "1e-4", "1e-6", "1e-8")) + theme_bw() + + scale_size_manual(values = c("A" = 7, "B" = 3)) + ggtitle("DV Resistant - Untreated") + + theme(text = element_text(size = 27), strip.background = element_rect(fill = NA, color = "black")) + + ggrepel::geom_label_repel(data = plot_df %>% filter(t_test_pval < 0.001 | abs(logFC) > logFC_cutoff), + aes(label = feature), min.segment.length = 1, size = 5) +p +plot_path = paste0(plot_folder, "volcano_plots_m_4pat.png") +ggsave(plot = p, filename = plot_path, width = 7.5, height = 7) +plot_path = paste0(plot_folder, "volcano_plots_m_4pat.pdf") +ggsave(plot = p, filename = plot_path, width = 7.5, height = 7) + + +``` + + + +Showing the signature in the exvivo data (training data) + +```{r} +# model <- logistic_model(m_4pat[complete_DV_UT_features_, ], "Treatment", skip_loo = TRUE, +# pred_group = "DV", control_group = "UT", +# ## Tune-able parameters +# N_markers = 12, pval_cutoff = 0.001, +# method = "s2n_max_v2", var.equal = TRUE) + +all_features_df_ <- model$full_model$features +xx = predict_helper(m_4pat, model$full_model, case_group = "DV", control_group = "UT") %>% + left_join(., pData(m_4pat), by = "sample") %>% + mutate(Group = case_when(Treatment == "DV" ~ "DV Resistant", + TRUE ~ "Untreated")) + +## Boxplot with no annotations +q = ggplot(xx %>% filter(Treatment %in% c("UT", "DV")), + aes(x = Group, y = logistic_signal2, fill = Group)) + geom_boxplot(width = 0.6) + + ylab("Dec-Ven resistance signal") + scale_fill_manual(values = pat_colors) + theme_bw() + + theme(text = element_text(size = 25)) + ggtitle("DV Resistance signature") + + guides(fill = "none") + +plot_path = paste0(plot_folder, "training_data_signal_boxplot.png") +ggsave(plot = q, filename = plot_path, width = 8, height = 7) + + +``` + + +Using all the markers we have (21 in total) + +```{r} +model <- logistic_model(m_4pat[complete_DV_UT_features_, ], "Treatment", skip_loo = TRUE, + pred_group = "DV", control_group = "UT", + ## Tune-able parameters + N_markers = 12, pval_cutoff = 0.001, + method = "s2n_max_v2", var.equal = TRUE) + +all_features_df_ <- model$full_model$features +xx = predict_helper(m_exp28, model$full_model, case_group = "DV", control_group = "UT") %>% + left_join(., p_data, by = c("sample" = "sample_name")) %>% + group_by(subcohort) %>% + mutate(corr = cor(logistic_signal2, survival_days)) %>% ungroup() %>% + mutate(subcohort_corr = paste(subcohort, formatC(corr, digits = 2))) %>% + mutate(refractory = case_when(`sub-cohort.PNNL` == "Refractory" ~ 1, + TRUE ~ 0)) +xx_ = xx %>% filter(subcohort != "Paired_relapse_sample") +pred = ROCR::prediction(xx_$logistic_signal2, xx_$refractory) +modelresult1 = list(prob = xx_$logistic_signal2, features = NA, top = NA, + auc = ROCR::performance(pred, "auc")@y.values[[1]], pred = pred) + +xx_ = xx %>% filter(subcohort %in% c("Refractory", "Response_no_relapse")) +pred = ROCR::prediction(xx_$logistic_signal2, xx_$refractory) +modelresult2 = list(prob = xx_$logistic_signal2, features = NA, top = NA, + auc = ROCR::performance(pred, "auc")@y.values[[1]], pred = pred) + + +pooled_t_test = stat_test_(xx %>% mutate(x = subcohort, y = logistic_signal2)) +t_test_pval = sub(" \n.*", "", pooled_t_test) %>% sub(".* = ", "", .) %>% as.numeric() +## Boxplot with no annotations +q = ggplot(xx %>% filter(subcohort != "Paired_relapse_sample"), + aes(x = subcohort, y = logistic_signal2, fill = group)) + geom_boxplot(width = 0.6) + + ylab("Dec-Ven signal") + scale_fill_manual(values = subcohort_colors) + theme_bw() + + ggtitle("DV resistance signature") + + theme(text = element_text(size = 27)) + + +## Boxplot with t_test annotations +q1 = t_test_ann(q, c("Refractory", "Response_no_relapse")) +q2 = t_test_ann(q1, c("Relapse", "Response_no_relapse")) +q2 = t_test_ann(q2, c("Refractory", "Relapse"), y_start = 1.07) +# plots[[paste(N_markers, pval_cutoff, method, sep = "_")]] = list(q2) +q2 +ggsave(plot = q2, filename = paste0(plot_folder, "m_exp28_exvivo_signature_boxplot.png"), width = 10, height = 7) +ggsave(plot = q2, filename = paste0(plot_folder, "m_exp28_exvivo_signature_boxplot.pdf"), width = 10, height = 7) + + +model_features = all_features_df_$feature +yy = plot_pca(m_exp28[model_features, ], phenotype = "subcohort", output_type = "full") + +df.u = yy$data.u +df.u$subcohort = pData(m_exp28)[, "subcohort"] +df.v = yy$data.v + +p <- ggplot(data = df.u, mapping = aes(x = df.u[, 1], y = df.u[, 2], + color = subcohort)) + geom_point(size = 2) + + geom_hline(yintercept = 0, lty = "longdash", color = "darkgrey") + + geom_vline(xintercept = 0, lty = "longdash", color = "darkgrey") + + labs(x = yy$axis_labs[1], y = yy$axis_labs[2]) + + theme_bw() + theme(aspect.ratio = 1, text = element_text(size = 15)) + + stat_ellipse(mapping = aes(fill = subcohort, color = NULL), + geom = "polygon", type = "norm", + level = 0.5, alpha = 0.07, show.legend = FALSE) + + ggtitle("PCA of 21 exvivo markers") + scale_fill_manual(values = subcohort_colors) + + scale_color_manual(values = subcohort_colors) +p +ggsave(plot = p, filename = paste0(plot_folder, "pca_m_exp28_exvivo_markers.pdf")) + + +p <- plotAUC_gg(modelresult1) + + theme(text = element_text(size = 25), + panel.background = element_rect(fill = "white", colour = "black"), + panel.grid.major = element_line(colour = "grey90"), + panel.grid.minor = element_line(colour = "grey95")) +ggsave(plot = p, filename = paste0(plot_folder, "exvivo_signature_refractory_auc.pdf"), width = 7, height = 6.7) + +``` + + +Using different numbers of markers ```{r} ############################## DV - UT signature ############################## # N_markers = 10 -complete_DV_UT_features <- complete.cases(exprs(m_4pat[, m_4pat$Treatment %in% c("UT", "DV")])) -complete_DV_UT_features_ = intersect(featureNames(m_4pat)[complete_DV_UT_features], complete_features) var.equal = TRUE plots = list() results_ttest = data.frame() @@ -215,7 +399,3 @@ p - - - - diff --git a/ven_hma_samples/Figure 3/ven_hma_sig.Rmd b/ven_hma_samples/Figure 3/ven_hma_sig.Rmd new file mode 100644 index 0000000..2f68581 --- /dev/null +++ b/ven_hma_samples/Figure 3/ven_hma_sig.Rmd @@ -0,0 +1,303 @@ + + + + + + +```{r include=FALSE} +source("../../synapseUtil.R") +source("../scripts/loading_data.R") +source("../scripts/plot_helper.R") +source("../scripts/compute_lsc17.R") +source("../scripts/predict_helper.R") +source("../scripts/s2n_model.R") + +samples_to_exclude = c("PTRC_Exp28_Prot_33", "PTRC_Exp28_Prot_15", "PTRC_Exp28_Prot_03") +m_exp28 = m_exp28[, setdiff(sampleNames(m_exp28), samples_to_exclude)] + +features_50perc = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) < 19] +complete_features = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) == 0] + +plot_folder = "./many_plots/" +if (!dir.exists(plot_folder)){ + dir.create(plot_folder) +} + +# diffexp_results = diffexp_helper(m_exp28, "subcohort", c("Refractory-Response_no_relapse", "Refractory-Relapse")) + +``` + + +```{r} +pData(m_exp28) = pData(m_exp28) %>% + mutate(refractory = case_when(subcohort == "Refractory" ~ "Refractory", + TRUE ~ "Not_Refractory")) +m_exp28_ = m_exp28[, m_exp28$subcohort != "Paired_relapse_sample"] +m_4pat_complete_features = featureNames(m_4pat)[(rowSums(is.na(exprs(m_4pat))) == 0)] + +make_title <- function(pval_cutoff, feature_df){ + title_top = paste0("pval_cutoff = ", pval_cutoff, "\n") + M_ = feature_df %>% pull(feature) + n_ = intersect(M_, m_4pat_complete_features) + full_title = paste0(title_top, "using", length(n_), "/", length(M_), "markers") + return(full_title) +} +``` + + +```{r} +diffexp = diffexp_helper(m_exp28_[complete_features, ], contrast_var = "refractory", + contrasts = "Refractory-Not_Refractory") +diffexp2 = diffexp_helper(m_exp28_[complete_features, ], contrast_var = "subcohort", + contrasts = "Refractory-Response_no_relapse") +``` + + + + +```{r} + +for (pval_cutoff in seq(11, 30)/1000){ + ### Refractory - Relapse signature applied to 4 patient cell lines + test = logistic_model(m_exp28_[complete_features, ], + group_var = "subcohort", pred_group = "Refractory", N_markers = 50, mode = "simple", + control_group = "Relapse", combine_markers = TRUE, skip_loo = TRUE, pval_cutoff = pval_cutoff) + xx = predict_helper_(m_4pat, case_markers = test$full_model$features %>% + filter(group_ == "Refractory") %>% pull(feature), + control_markers = test$full_model$features %>% + filter(group_ != "Refractory") %>% pull(feature)) %>% + filter(grepl("UT", sample) | grepl("DV", sample)) %>% + mutate(group = case_when(grepl("DV", sample) ~ "DV", + TRUE ~ "UT")) + treatment_colors = c("DV" = "#00798c", "UT" = "forestgreen") + q = ggplot(xx, aes(x = group, y = signal, fill = group)) + geom_boxplot(width = 0.75) + + ylab("Refractory-Relapse signal") + scale_fill_manual(values = treatment_colors) + + theme(text = element_text(size = 15)) + + ggtitle(make_title(pval_cutoff, test$full_model$features)) + q1 = t_test_ann(q, c("DV", "UT")) + ggsave(plot = q1, filename = paste0("many_plots/refractory_relapse_signature_pval_", pval_cutoff, ".png")) + + + + + ### Refractory - Response_no_relapse signature applied to 4 patient cell lines + test = logistic_model(m_exp28_[complete_features, ], + group_var = "subcohort", pred_group = "Refractory", N_markers = 50, mode = "simple", + control_group = "Response_no_relapse", combine_markers = TRUE, + skip_loo = TRUE, pval_cutoff = pval_cutoff) + xx = predict_helper_(m_4pat, case_markers = test$full_model$features %>% + filter(group_ == "Refractory") %>% pull(feature), + control_markers = test$full_model$features %>% + filter(group_ != "Refractory") %>% pull(feature)) %>% + filter(grepl("UT", sample) | grepl("DV", sample)) %>% + mutate(group = case_when(grepl("DV", sample) ~ "DV", + TRUE ~ "UT")) + treatment_colors = c("DV" = "#00798c", "UT" = "forestgreen") + q = ggplot(xx, aes(x = group, y = signal, fill = group)) + geom_boxplot(width = 0.75) + + ylab("Refractory-Response_no_relapse signal") + scale_fill_manual(values = treatment_colors) + + theme(text = element_text(size = 15)) + + ggtitle(make_title(pval_cutoff, test$full_model$features)) + q1 = t_test_ann(q, c("DV", "UT")) + ggsave(plot = q1, filename = paste0("many_plots/refractory_response_no_relapse_signature_pval_", pval_cutoff, ".png")) + + + + + ### Refractory - response signature applied to 4 patient cell lines + test = logistic_model(m_exp28_[complete_features, ], + group_var = "refractory", pred_group = "Refractory", N_markers = 50, mode = "simple", + control_group = "Not_Refractory", combine_markers = TRUE, + skip_loo = TRUE, pval_cutoff = pval_cutoff) + xx = predict_helper_(m_4pat, case_markers = test$full_model$features %>% + filter(group_ == "Refractory") %>% pull(feature), + control_markers = test$full_model$features %>% + filter(group_ != "Refractory") %>% pull(feature)) %>% + filter(grepl("UT", sample) | grepl("DV", sample)) %>% + mutate(group = case_when(grepl("DV", sample) ~ "DV", + TRUE ~ "UT")) + treatment_colors = c("DV" = "#00798c", "UT" = "forestgreen") + q = ggplot(xx, aes(x = group, y = signal, fill = group)) + geom_boxplot(width = 0.75) + + ylab("Refractory-Relapse signal") + scale_fill_manual(values = treatment_colors) + + theme(text = element_text(size = 15)) + + ggtitle(make_title(pval_cutoff, test$full_model$features)) + q1 = t_test_ann(q, c("DV", "UT")) + ggsave(plot = q1, filename = paste0("many_plots/refractory_response_signature_pval_", pval_cutoff, ".png")) +} + + +``` + + + +```{r} +sig_pval_df = data.frame(pval_cutoff = seq(1, 16)/1000, DV_UT_pvalue = c(0.014, 0.0035, 0.0035, 0.28, 0.59, 1, 0.62, 0.62, + 0.55, 0.55, 0.25, 0.28, 0.33, 0.12, 0.11, 0.096), + N_markers = c("2", "5", "5", "6", "7", "8", "11", "11","12", "12", "14", "15", "16", "18", "19", "21")) +ggplot(sig_pval_df, aes(x = pval_cutoff, y = DV_UT_pvalue, label = N_markers)) + geom_bar(stat = 'identity') + geom_label() + + ggtitle("Ven+HMA signature performance on 4 patient data") + +``` + + + + +```{r} +pval_cutoff = 0.025 +test = logistic_model(m_exp28_[complete_features, ], + group_var = "subcohort", pred_group = "Refractory", N_markers = 50, mode = "simple", + control_group = "Response_no_relapse", combine_markers = TRUE, + skip_loo = TRUE, pval_cutoff = pval_cutoff) +xx = predict_helper_(m_4pat, case_markers = test$full_model$features %>% + filter(group_ == "Refractory") %>% pull(feature), + control_markers = test$full_model$features %>% + filter(group_ != "Refractory") %>% pull(feature)) %>% + filter(grepl("UT", sample) | grepl("DV", sample)) %>% + mutate(group = case_when(grepl("DV", sample) ~ "DV", + TRUE ~ "UT")) +treatment_colors = c("DV" = "#00798c", "UT" = "forestgreen") +q = ggplot(xx, aes(x = group, y = signal, fill = group)) + geom_boxplot(width = 0.67) + + ylab("Signal") + scale_fill_manual(values = treatment_colors) + + theme_bw() + theme(text = element_text(size = 23), plot.title = element_text(hjust = 0.5)) + + ggtitle("Ven+HMA resistance signature") +q1 = t_test_ann(q, c("DV", "UT")) +ggsave(plot = q1, filename = paste0("many_plots/refractory_response_no_relapse_signature_pval_", pval_cutoff, ".png")) +ggsave(plot = q1, filename = paste0("many_plots/refractory_response_no_relapse_signature_pval_", pval_cutoff, ".pdf")) + +pred = ROCR::prediction(xx$signal, xx$group == "DV") +modelresult1 = list(prob = xx$signal, features = NA, top = NA, + auc = ROCR::performance(pred, "auc")@y.values[[1]], pred = pred) + +p <- plotAUC_gg(modelresult1) + + theme(text = element_text(size = 25), + panel.background = element_rect(fill = "white", colour = "black"), + panel.grid.major = element_line(colour = "grey90"), + panel.grid.minor = element_line(colour = "grey95")) +ggsave(plot = p, filename = "many_plots/refractory_signature_DV_auc.pdf", width = 7, height = 6.7) + + +# for (pval_cutoff in c(0.013, 0.014, 0.015, 0.016)){ +# ### Refractory - response signature applied to 4 patient cell lines +# test = logistic_model(m_exp28_[complete_features, ], +# group_var = "refractory", pred_group = "Refractory", N_markers = 20, mode = "simple", +# control_group = "Not_Refractory", combine_markers = TRUE, +# skip_loo = TRUE, pval_cutoff = pval_cutoff) +# xx = predict_helper_(m_4pat, case_markers = test$full_model$features %>% +# filter(group_ == "Refractory") %>% pull(feature), +# control_markers = test$full_model$features %>% +# filter(group_ != "Refractory") %>% pull(feature)) %>% +# filter(grepl("UT", sample) | grepl("DV", sample)) %>% +# mutate(group = case_when(grepl("DV", sample) ~ "DV", +# TRUE ~ "UT")) +# treatment_colors = c("DV" = "#00798c", "UT" = "forestgreen") +# q = ggplot(xx, aes(x = group, y = signal, fill = group)) + geom_boxplot(width = 0.75) + +# ylab("Refractory-Relapse signal") + scale_fill_manual(values = treatment_colors) + +# theme(text = element_text(size = 15)) + +# ggtitle(make_title(pval_cutoff, test$full_model$features)) +# q1 = t_test_ann(q, c("DV", "UT")) +# ggsave(plot = q1, filename = paste0("many_plots/refractory_response_signature_pval_", pval_cutoff, ".png")) +# } + + + +``` + + + + + +```{r} +diffexp_results_venhma = diffexp_helper(m_exp28_[complete_features, ], "subcohort", c("Refractory-Response_no_relapse")) + +logFC_cutoff = 2 +plot_df = diffexp_results_venhma %>% filter(feature %in% m_4pat_complete_features) %>% + mutate(pvalue_ = -log10(t_test_pval), + alpha = case_when(t_test_pval <= pval_cutoff ~ 1, + logFC > logFC_cutoff | logFC < -logFC_cutoff ~ 1, + TRUE ~ 0.01), + color = case_when(t_test_pval < pval_cutoff & logFC > 0 ~ "logFC > 0", + t_test_pval < pval_cutoff & logFC < 0 ~ "logFC < 0", + logFC > logFC_cutoff ~ "logFC > 0", + logFC < -logFC_cutoff ~ "logFC < 0", + TRUE ~ ""), + size = case_when(t_test_pval < pval_cutoff ~ "A", + abs(logFC) > logFC_cutoff ~ "A", + TRUE ~ "B")) %>% select(-feature) %>% unique() +plot_df$feature = rownames(plot_df) +p = ggplot(plot_df, aes(x = logFC, y = pvalue_, alpha = alpha, color = color, size = size)) + + geom_point() + scale_color_manual(values = c("grey20", "steelblue4", "firebrick4")) + + geom_hline(aes(yintercept = -log10(pval_cutoff)), linetype = "dashed", color = "forestgreen", linewidth = 0.7) + + guides(alpha = "none", size = "none", color = "none") + ylab("pvalue") + xlim(-3.3, 3.3) + + scale_y_continuous(breaks = c(0,2,3,4,6,8), labels = c("1", "0.01", "0.001", "1e-4", "1e-6", "1e-8")) + theme_bw() + + scale_size_manual(values = c("A" = 7, "B" = 3)) + ggtitle("Refractory - Response no relapse") + + theme(text = element_text(size = 27), strip.background = element_rect(fill = NA, color = "black")) + + ggrepel::geom_label_repel(data = plot_df %>% filter(t_test_pval < pval_cutoff | abs(logFC) > logFC_cutoff), + aes(label = feature), min.segment.length = 1, size = 5) +p +plot_path = "many_plots/volcano_plots_refractory_response_no_relapse.png" +ggsave(plot = p, filename = plot_path, width = 7.5, height = 7) +plot_path = "many_plots/volcano_plots_refractory_response_no_relapse.pdf" +ggsave(plot = p, filename = plot_path, width = 7.5, height = 7) + + +``` + + + + +```{r} +all_features_df_ <- test$full_model$features +xx = predict_helper(m_exp28, test$full_model, case_group = "Refractory", control_group = "Response_no_relapse") %>% + left_join(., p_data, by = c("sample" = "sample_name")) %>% + group_by(subcohort) %>% + mutate(corr = cor(logistic_signal2, survival_days)) %>% ungroup() %>% + mutate(subcohort_corr = paste(subcohort, formatC(corr, digits = 2))) %>% + mutate(refractory = case_when(`sub-cohort.PNNL` == "Refractory" ~ 1, + TRUE ~ 0)) + +q = ggplot(xx %>% filter(subcohort != "Paired_relapse_sample"), + aes(x = subcohort, y = logistic_signal2, fill = group)) + geom_boxplot(width = 0.6) + + ylab("Signal") + scale_fill_manual(values = subcohort_colors) + theme_bw() + + ggtitle("Refractory - Response no relapse") + + theme(text = element_text(size = 27)) + +ggsave(plot = q, filename = "many_plots/m_exp28_refractory_response_no_relapse_signature_boxplot.pdf", width = 10, height = 7) +``` + + + +```{r} +# seed_ = 654654651 +# set.seed(seed_) +# data2 = data.frame(id = sampleNames(m_exp28), subcohort = m_exp28$subcohort) %>% +# filter(subcohort %in% c("Refractory", "Response_no_relapse")) +# n_folds = 3 +# x = groupdata2::fold(data2, k = n_folds, cat_col = "subcohort") +# +# for (i in range(n_folds)){ +# excluded_samples = x %>% filter(.folds == i) +# train_samples = x %>% filter(.folds != i) +# m_train = m_exp28[complete_features, train_samples$id] +# m_test = m_exp28[complete_features, excluded_samples$id] +# +# model = logistic_model(m_train, group_var = "subcohort", pred_group = "Refractory", N_markers = 10, mode = "s2n_max_v2", +# control_group = "Response_no_relapse", combine_markers = TRUE, skip_loo = TRUE, pval_cutoff = 0.005) +# +# xx = predict_helper(m_test, model$full_model, case_group = "Refractory", control_group = "Response_no_relapse") %>% +# left_join(., p_data, by = c("sample" = "sample_name")) +# } + +``` + + + +```{r} +# feat1 = all_features_df_ # using the exvivo model +# feat2 = all_features_df_ # using the ven+hma model +common = intersect(feat1$feature, feat2$feature) + +feat1 = left_join(feat1, diffexp_results_venhma %>% select(logFC, feature, t_test_pval), by = "feature") +feat2 = left_join(feat2, diffexp_results_4pat %>% select(logFC, feature, t_test_pval), by = "feature") + +``` + + diff --git a/ven_hma_samples/Figure 4/hazard_markers.Rmd b/ven_hma_samples/Figure 4/hazard_markers.Rmd new file mode 100644 index 0000000..d175954 --- /dev/null +++ b/ven_hma_samples/Figure 4/hazard_markers.Rmd @@ -0,0 +1,102 @@ + + + + + + +```{r include=FALSE} +source("../../synapseUtil.R") +source("../scripts/loading_data.R") +source("../scripts/plot_helper.R") +source("../scripts/compute_lsc17.R") +source("../scripts/predict_helper.R") +source("../scripts/s2n_model.R") + +samples_to_exclude = c("PTRC_Exp28_Prot_33", "PTRC_Exp28_Prot_15", "PTRC_Exp28_Prot_03") +m_exp28 = m_exp28[, setdiff(sampleNames(m_exp28), samples_to_exclude)] + +features_50perc = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) < 19] +complete_features = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) == 0] + +plot_folder = "./many_plots/" +if (!dir.exists(plot_folder)){ + dir.create(plot_folder) +} + +# diffexp_results = diffexp_helper(m_exp28, "subcohort", c("Refractory-Response_no_relapse", "Refractory-Relapse")) + +``` + + +```{r} +library("survival") +library("survminer") + +pData(m_exp28) = pData(m_exp28) %>% + mutate(refractory = case_when(subcohort == "Refractory" ~ "Refractory", + TRUE ~ "Not_Refractory"), + status = case_when(vitalStatus == "Dead" ~ 2, # 2=dead, 1=censored + TRUE ~ 1), + time = survival_days) +m_exp28_ = m_exp28[, m_exp28$subcohort != "Paired_relapse_sample"] + +hazards_helper = function(msnset, feature){ + hazards_df = data.frame(time = msnset$time, status = msnset$status, + value = exprs(msnset)[feature, ]) + res_ = coxph(Surv(time, status) ~ value, data = hazards_df) + out = data.frame(feature = feature, wald_pval = summary(res_)$waldtest["pvalue"], + LR_pval = summary(res_)$logtest["pvalue"]) + return(out) +} + +hazards_ = function(msnset){ + res_ = data.frame() + for (feature in featureNames(msnset)){ + res_ = rbind(res_, hazards_helper(msnset, feature)) + } + return(res_) +} + +hazards_res = hazards_(m_exp28_[features_50perc, ]) +m___ = m_exp28_[hazards_res %>% filter(wald_pval < 0.05) %>% pull(feature), ] +plot_pca(m___, "subcohort", output_type = "") +``` + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/ven_hma_samples/Figure 5/drug_auc_validation.Rmd b/ven_hma_samples/Figure 5/drug_auc_validation.Rmd new file mode 100644 index 0000000..eaff768 --- /dev/null +++ b/ven_hma_samples/Figure 5/drug_auc_validation.Rmd @@ -0,0 +1,194 @@ + + +```{r} +source("../../synapseUtil.R") +source("../scripts/loading_data.R") +source("../scripts/plot_helper.R") +source("../scripts/compute_lsc17.R") +source("../scripts/predict_helper.R") +source("../scripts/s2n_model.R") + +# m_exp28 = readRDS("msnset_exp28_global.RDS") +samples_to_exclude = c("PTRC_Exp28_Prot_33", "PTRC_Exp28_Prot_15", "PTRC_Exp28_Prot_03") +m_exp28 = m_exp28[, setdiff(sampleNames(m_exp28), samples_to_exclude)] + +complete_features = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) == 0] +m_210$new = paste(m_210$InitialAMLDiagnosis, m_210$PostChemotherapy) +beataml_clinical = openxlsx::read.xlsx(syn$get("syn74822940")$path) +``` + + + +```{r} +complete_features = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) == 0] +complete_DV_UT_features <- complete.cases(exprs(m_4pat[, m_4pat$Treatment %in% c("UT", "DV")])) +complete_DV_UT_features_ = intersect(featureNames(m_4pat)[complete_DV_UT_features], complete_features) +complete_V_UT_features <- complete.cases(exprs(m_4pat[, m_4pat$Treatment %in% c("UT", "V")])) +complete_V_UT_features_ = intersect(featureNames(m_4pat)[complete_V_UT_features], complete_features) +model <- logistic_model(m_4pat[complete_DV_UT_features_, ], "Treatment", skip_loo = TRUE, + pred_group = "DV", control_group = "UT", + ## Tune-able parameters + N_markers = 12, pval_cutoff = 0.001, + method = "s2n_max_v2", var.equal = TRUE) + +refractory_markers = model$full_model$features %>% filter(group_ == "DV") %>% pull(feature) +longresponse_markers = model$full_model$features %>% filter(group_ != "DV") %>% pull(feature) + +xx = predict_helper_(m_md_anderson, refractory_markers, longresponse_markers) %>% + left_join(., pData(m_md_anderson), by = c("sample" = "sample_name")) + +ggplot(xx, aes(x = subcohort, y = signal, fill = subcohort)) + geom_boxplot() +m_md_anderson$subcohort_ = gsub(" ", "_", m_md_anderson$subcohort) +yy = diffexp_helper(m_md_anderson, "subcohort_", "Long_response-Refractory") + +``` + + + +```{r} +# model <- logistic_model(m_4pat[complete_DV_UT_features_, ], "Treatment", skip_loo = TRUE, +# pred_group = "DV", control_group = "UT", +# ## Tune-able parameters +# N_markers = 12, pval_cutoff = 0.001, +# method = "simple", var.equal = TRUE) + +refractory_markers = model$full_model$features %>% filter(group_ == "DV") %>% pull(feature) +longresponse_markers = model$full_model$features %>% filter(group_ != "DV") %>% pull(feature) + +beataml_210_predictions = predict_helper_(m_210, refractory_markers, longresponse_markers) %>% + left_join(., drug_data_df, by = c("sample" = "sample_id")) %>% + left_join(., pData(m_210), by = c("sample" = "Barcode.ID")) + +plot_df = beataml_210_predictions %>% + filter(inhibitor == "Azacytidine - Venetoclax") +ll = plot_df %>% filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", + "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +a = cor.test(ll$auc, ll$signal) +a +line_slope = coef(lm(ll$signal ~ ll$auc))[[2]] +line_intercept = coef(lm(ll$signal ~ ll$auc))[[1]] +plot_title = paste0("Ex-vivo signature vs Ven+Aza drug AUC\n", "correlation = ", formatC(a$estimate, 2), + ", p_value = ", formatC(a$p.value, 2)) +ggplot(ll, aes(x = auc, y = signal, color = FLT3.ITD, label = sample)) + geom_point(size = 2.5) + + ggtitle(plot_title) + ylab("signal") + xlab("Ven+Aza drug AUC") + + geom_abline(intercept = line_intercept, slope = line_slope, size = 1.5, alpha = 0.2) + + theme_bw() + theme(text = element_text(size = 15)) + +plot_df = beataml_210_predictions %>% + filter(inhibitor == "Venetoclax") +ll = plot_df %>% filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", + "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +a = cor.test(ll$auc, ll$signal) +a +line_slope = coef(lm(ll$signal ~ ll$auc))[[2]] +line_intercept = coef(lm(ll$signal ~ ll$auc))[[1]] +plot_title = paste0("Ex-vivo signature vs Ven drug AUC\n", "correlation = ", formatC(a$estimate, 2), + ", p_value = ", formatC(a$p.value, 2)) +ggplot(ll, aes(x = auc, y = signal, color = FLT3.ITD, label = sample)) + geom_point(size = 2.5) + + ggtitle(plot_title) + ylab("signal") + xlab("Ven drug AUC") + + geom_abline(intercept = line_intercept, slope = line_slope, size = 1.5, alpha = 0.2) + + theme_bw() + theme(text = element_text(size = 15)) + + +``` + + + + + +```{r} +# diffexp2 = diffexp_helper(m_exp28[complete_features, ], contrast_var = "subcohort", +# contrasts = "Refractory-Response_no_relapse") + +refractory_markers = diffexp2 %>% filter(logFC > 0, t_test_pval < 0.025) %>% pull(feature) +longresponse_markers = diffexp2 %>% filter(logFC < 0, t_test_pval < 0.025) %>% pull(feature) + +beataml_210_predictions = predict_helper_(m_210, refractory_markers, longresponse_markers) %>% + left_join(., drug_data_df, by = c("sample" = "sample_id")) %>% + left_join(., pData(m_210), by = c("sample" = "Barcode.ID")) + +plot_df = beataml_210_predictions %>% + filter(inhibitor == "Azacytidine - Venetoclax") +ggplot(plot_df, aes(x = auc, fill = new)) + geom_density(alpha = 0.2) +ll = plot_df %>% filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", + "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +a = cor.test(ll$auc, ll$signal) +a +line_slope = coef(lm(ll$signal ~ ll$auc))[[2]] +line_intercept = coef(lm(ll$signal ~ ll$auc))[[1]] +plot_title = paste0("In-vivo signature vs Ven+Aza drug AUC\n", "correlation = ", formatC(a$estimate, 2), + ", p_value = ", formatC(a$p.value, 2)) +ggplot(ll, aes(x = auc, y = signal, color = FLT3.ITD, label = sample)) + geom_point(size = 2.5) + + ggtitle(plot_title) + ylab("signal") + xlab("Ven+Aza drug AUC") + + geom_abline(intercept = line_intercept, slope = line_slope, size = 1.5, alpha = 0.2) + + theme_bw() + theme(text = element_text(size = 15)) + + +plot_df = beataml_210_predictions %>% + filter(inhibitor == "Venetoclax") +ggplot(plot_df, aes(x = auc, fill = new)) + geom_density(alpha = 0.2) +ll = plot_df %>% filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", + "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +a = cor.test(ll$auc, ll$signal) +a +line_slope = coef(lm(ll$signal ~ ll$auc))[[2]] +line_intercept = coef(lm(ll$signal ~ ll$auc))[[1]] +plot_title = paste0("In-vivo signature vs Ven drug AUC\n", "correlation = ", formatC(a$estimate, 2), + ", p_value = ", formatC(a$p.value, 2)) +ggplot(ll, aes(x = auc, y = signal, color = FLT3.ITD, label = sample)) + geom_point(size = 2.5) + + ggtitle(plot_title) + ylab("signal") + xlab("Ven drug AUC") + + geom_abline(intercept = line_intercept, slope = line_slope, size = 1.5, alpha = 0.2) + + theme_bw() + theme(text = element_text(size = 15)) + + +``` + + + + +```{r} +# model <- logistic_model(m_4pat[complete_V_UT_features_, ], "Treatment", skip_loo = TRUE, +# pred_group = "V", control_group = "UT", +# ## Tune-able parameters +# N_markers = 12, pval_cutoff = 0.001, +# method = "simple", var.equal = TRUE) + +refractory_markers = model$full_model$features %>% filter(group_ == "V") %>% pull(feature) +longresponse_markers = model$full_model$features %>% filter(group_ != "V") %>% pull(feature) + +beataml_210_predictions = predict_helper_(m_210, refractory_markers, longresponse_markers) %>% + left_join(., drug_data_df, by = c("sample" = "sample_id")) %>% + left_join(., pData(m_210), by = c("sample" = "Barcode.ID")) + +plot_df = beataml_210_predictions %>% + filter(inhibitor == "Azacytidine - Venetoclax") +ll = plot_df %>% filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", + "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +a = cor.test(ll$auc, ll$signal) +a +line_slope = coef(lm(ll$signal ~ ll$auc))[[2]] +line_intercept = coef(lm(ll$signal ~ ll$auc))[[1]] +plot_title = paste0("Ex-vivo signature (ven vs Untreated) vs Ven+Aza drug AUC\n", "correlation = ", formatC(a$estimate, 2), + ", p_value = ", formatC(a$p.value, 2)) +ggplot(ll, aes(x = auc, y = signal, color = FLT3.ITD, label = sample)) + geom_point(size = 2.5) + + ggtitle(plot_title) + ylab("signal") + xlab("Ven+Aza drug AUC") + + geom_abline(intercept = line_intercept, slope = line_slope, size = 1.5, alpha = 0.2) + + theme_bw() + theme(text = element_text(size = 15)) + +plot_df = beataml_210_predictions %>% + filter(inhibitor == "Venetoclax") +ll = plot_df %>% filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", + "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +a = cor.test(ll$auc, ll$signal) +a +line_slope = coef(lm(ll$signal ~ ll$auc))[[2]] +line_intercept = coef(lm(ll$signal ~ ll$auc))[[1]] +plot_title = paste0("Ex-vivo signature (ven vs Untreated) vs Ven drug AUC\n", "correlation = ", formatC(a$estimate, 2), + ", p_value = ", formatC(a$p.value, 2)) +ggplot(ll, aes(x = auc, y = signal, color = FLT3.ITD, label = sample)) + geom_point(size = 2.5) + + ggtitle(plot_title) + ylab("signal") + xlab("Ven drug AUC") + + geom_abline(intercept = line_intercept, slope = line_slope, size = 1.5, alpha = 0.2) + + theme_bw() + theme(text = element_text(size = 15)) + +``` + diff --git a/ven_hma_samples/exvivo_signature_panels/exvivo_signature.Rmd b/ven_hma_samples/exvivo_signature_panels/exvivo_signature.Rmd new file mode 100644 index 0000000..3f0f773 --- /dev/null +++ b/ven_hma_samples/exvivo_signature_panels/exvivo_signature.Rmd @@ -0,0 +1,358 @@ + + + +```{r} +source("../../synapseUtil.R") +source("../scripts/loading_data.R") +source("../scripts/plot_helper.R") +source("../scripts/compute_lsc17.R") +source("../scripts/predict_helper.R") +source("../scripts/s2n_model.R") + +samples_to_exclude = c("PTRC_Exp28_Prot_33", "PTRC_Exp28_Prot_15", "PTRC_Exp28_Prot_03") +m_exp28 = m_exp28[, setdiff(sampleNames(m_exp28), samples_to_exclude)] +m_exp28_ = m_exp28[, m_exp28$subcohort != "Paired_relapse_sample"] + +features_50perc = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) < 19] +complete_features = featureNames(m_exp28)[rowSums(is.na(exprs(m_exp28))) == 0] + +plot_folder = "./" + +## Need to apply signature in in-vivo dataset, so we filter by complete features in U, VT AND the in-vivo dataset +complete_DV_UT_features <- complete.cases(exprs(m_4pat[, m_4pat$Treatment %in% c("UT", "DV")])) +complete_DV_UT_features_ = intersect(featureNames(m_4pat)[complete_DV_UT_features], complete_features) + +``` + + + + +```{r} +# stats_df = data.frame() +# for (chosen_N in 1:14){ +# print(chosen_N) +# loocv_model_DV = logistic_model(m_4pat[complete_DV_UT_features_, ], "Treatment", skip_loo = FALSE, +# pred_group = "DV", control_group = "UT", +# ## Tune-able parameters +# N_markers = chosen_N, pval_cutoff = 1, +# method = "simple", var.equal = TRUE) +# DV_markers_up = loocv_model_DV$full_model$features %>% filter(group_ == "DV") %>% pull(feature) +# DV_markers_down = loocv_model_DV$full_model$features %>% filter(group_ == "UT") %>% pull(feature) +# +# ### Exvivo DV signature stats +# invivo_prediction = predict_helper_(m_exp28, DV_markers_up, DV_markers_down) %>% +# left_join(., pData(m_exp28), by = c("sample" = "sample_name")) %>% +# mutate(refractory = case_when(`sub-cohortPNNL` == "Refractory" ~ 1, +# TRUE ~ 0)) +# xx_ = invivo_prediction %>% filter(subcohort != "Paired_relapse_sample") +# +# beataml_210_predictions = predict_helper_(m_210, DV_markers_up, DV_markers_down) %>% +# left_join(., drug_data_df, by = c("sample" = "sample_id")) %>% +# left_join(., pData(m_210), by = c("sample" = "Barcode.ID")) +# +# aza_ven_df = beataml_210_predictions %>% +# filter(inhibitor == "Azacytidine - Venetoclax") %>% +# filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", +# "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +# ven_df = beataml_210_predictions %>% +# filter(inhibitor == "Venetoclax") %>% +# filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", +# "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +# +# stats_df = rbind(stats_df, +# data.frame(N = chosen_N, signature = "Exvivo DV", loocv_auc = loocv_model_DV$auc, +# beataml_ven_corr = cor(ven_df$signal, ven_df$auc), +# beataml_azaven_corr = cor(aza_ven_df$signal, aza_ven_df$auc), +# invivo_auc = ROCR::performance(ROCR::prediction(xx_$signal, xx_$refractory), "auc")@y.values[[1]])) +# +# +# +# loocv_model_V = logistic_model(m_4pat[complete_DV_UT_features_, ], "Treatment", skip_loo = FALSE, +# pred_group = "V", control_group = "UT", +# ## Tune-able parameters +# N_markers = chosen_N, pval_cutoff = 1, +# method = "simple", var.equal = TRUE) +# V_markers_up = loocv_model_V$full_model$features %>% filter(group_ == "V") %>% pull(feature) +# V_markers_down = loocv_model_V$full_model$features %>% filter(group_ == "UT") %>% pull(feature) +# +# +# ### Exvivo Ven signature stats +# invivo_prediction = predict_helper_(m_exp28, V_markers_up, V_markers_down) %>% +# left_join(., pData(m_exp28), by = c("sample" = "sample_name")) %>% +# mutate(refractory = case_when(`sub-cohortPNNL` == "Refractory" ~ 1, +# TRUE ~ 0)) +# xx_ = invivo_prediction %>% filter(subcohort != "Paired_relapse_sample") +# +# beataml_210_predictions = predict_helper_(m_210, V_markers_up, V_markers_down) %>% +# left_join(., drug_data_df, by = c("sample" = "sample_id")) %>% +# left_join(., pData(m_210), by = c("sample" = "Barcode.ID")) +# +# aza_ven_df = beataml_210_predictions %>% +# filter(inhibitor == "Azacytidine - Venetoclax") %>% +# filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", +# "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +# ven_df = beataml_210_predictions %>% +# filter(inhibitor == "Venetoclax") %>% +# filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", +# "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +# +# stats_df = rbind(stats_df, +# data.frame(N = chosen_N, signature = "Exvivo V", loocv_auc = loocv_model_V$auc, +# beataml_ven_corr = cor(ven_df$signal, ven_df$auc), +# beataml_azaven_corr = cor(aza_ven_df$signal, aza_ven_df$auc), +# invivo_auc = ROCR::performance(ROCR::prediction(xx_$signal, xx_$refractory), "auc")@y.values[[1]])) +# } + +# write.table(stats_df, "exvivo_signature_stats.txt", sep = "\t", quote = F, col.names = T, row.names = F) + +# stats_df = read.table("exvivo_signature_stats.txt", sep = "\t", header = T) +stats_df = read.table(syn$get("syn76730225")$path, sep = "\t", header = T) + +``` + + + + +```{r} +chosen_signature = "DV" +plot_df = stats_df %>% + filter(signature == paste("Exvivo", chosen_signature)) + +ggplot(plot_df, aes(x = N, y = beataml_ven_corr)) + + geom_line(linewidth = 2, color = "orange") + + geom_point(size = 3, color = "steelblue") + + theme_bw() + theme(text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + + ggtitle(paste("Correlation of Exvivo", chosen_signature, "signature w/\nVen drug AUC in BeatAML")) + + ylab("Correlation") + xlab("N markers") +ggsave(paste0(plot_folder, "exvivo_", chosen_signature, "_signature_beataml_ven_drug_auc_stats.png")) +ggsave(paste0(plot_folder, "exvivo_", chosen_signature, "_signature_beataml_ven_drug_auc_stats.pdf")) + + +ggplot(plot_df, aes(x = N, y = loocv_auc)) + + geom_bar(stat = 'identity', width = 0.67, color = "black", fill = "steelblue") + + geom_text(aes(label = round(loocv_auc, 2)), vjust = 1.2, color = 'white') + + theme_bw() + theme(text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + + ggtitle(paste("LOOCV AUC - Exvivo", chosen_signature, "signature")) + + ylab("AUC") + xlab("N markers") +ggsave(paste0(plot_folder, "exvivo_", chosen_signature, "_signature_loocv_auc_stats.png")) +ggsave(paste0(plot_folder, "exvivo_", chosen_signature, "_signature_loocv_auc_stats.pdf")) + + + +``` + + + +### Testing performace of exvivo signature in invivo data + + +```{r} + +loocv_model_DV = logistic_model(m_4pat[complete_DV_UT_features_, ], "Treatment", skip_loo = FALSE, + pred_group = "DV", control_group = "UT", + ## Tune-able parameters + N_markers = 9, pval_cutoff = 1, + method = "simple", var.equal = TRUE) +DV_markers_up = loocv_model_DV$full_model$features %>% filter(group_ == "DV") %>% pull(feature) +DV_markers_down = loocv_model_DV$full_model$features %>% filter(group_ == "UT") %>% pull(feature) + +xx = predict_helper_(m_exp28, DV_markers_up, DV_markers_down) %>% + left_join(., p_data, by = c("sample" = "sample_name")) %>% + mutate(refractory = case_when(subcohort == "Refractory" ~ 1, + TRUE ~ 0), + subcohort = case_when(subcohort == "Response_no_relapse" ~ "Response\nno relapse", + subcohort == "Relapse" ~ "Relapsed", + TRUE ~ subcohort)) +xx = xx %>% filter(subcohort != "Paired_relapse_sample") +pred = ROCR::prediction(xx$signal, xx$refractory) +modelresult1 = list(prob = xx$signal, features = NA, top = NA, + auc = ROCR::performance(pred, "auc")@y.values[[1]], pred = pred) + +pooled_t_test = stat_test_(xx %>% mutate(x = subcohort, y = signal)) +t_test_pval = sub(" \n.*", "", pooled_t_test) %>% sub(".* = ", "", .) %>% as.numeric() +## Boxplot with no annotations +q = ggplot(xx %>% filter(subcohort != "Paired_relapse_sample"), + aes(x = subcohort, y = signal, fill = subcohort)) + geom_boxplot(width = 0.6) + + ylab("Dec-Ven signal") + scale_fill_manual(values = subcohort_colors) + theme_bw() + + ggtitle("Exvivo Dec-Ven resistance signature") + + theme(text = element_text(size = 22), plot.title = element_text(hjust = 0.5)) + + +## Boxplot with t_test annotations +q1 = t_test_ann(q, c("Refractory", "Response\nno relapse")) +q2 = t_test_ann(q1, c("Relapsed", "Response\nno relapse")) +q2 = t_test_ann(q2, c("Refractory", "Relapsed"), y_start = 1.035) +q2 +ggsave(plot = q2, filename = paste0(plot_folder, "boxplot_exvivo_DV_signature_applied_invivo.png"), + width = 8, height = 8) +ggsave(plot = q2, filename = paste0(plot_folder, "boxplot_exvivo_DV_signature_applied_invivo.pdf"), + width = 8, height = 8) + + +p = plotAUC_gg(modelresult1) + theme_bw() + ggtitle("Detecting refractory patients \nExvivo DV signature") + + theme(text = element_text(size = 18), plot.title = element_text(hjust = 0.5)) +ggsave(plot = p, filename = paste0(plot_folder, "aucplot_exvivo_DV_signature_applied_invivo.png"), + width = 8, height = 8) +ggsave(plot = p, filename = paste0(plot_folder, "aucplot_exvivo_DV_signature_applied_invivo.pdf"), + width = 8, height = 8) + + +``` + + +```{r} +loocv_model_V = logistic_model(m_4pat[complete_DV_UT_features_, ], "Treatment", skip_loo = FALSE, + pred_group = "V", control_group = "UT", + ## Tune-able parameters + N_markers = 6, pval_cutoff = 1, + method = "simple", var.equal = TRUE) +V_markers_up = loocv_model_V$full_model$features %>% filter(group_ == "V") %>% pull(feature) +V_markers_down = loocv_model_V$full_model$features %>% filter(group_ == "UT") %>% pull(feature) + +xx = predict_helper_(m_exp28, V_markers_up, V_markers_down) %>% + left_join(., p_data, by = c("sample" = "sample_name")) %>% + mutate(refractory = case_when(subcohort == "Refractory" ~ 1, + TRUE ~ 0), + subcohort = case_when(subcohort == "Response_no_relapse" ~ "Response\nno relapse", + subcohort == "Relapse" ~ "Relapsed", + TRUE ~ subcohort)) +xx = xx %>% filter(subcohort != "Paired_relapse_sample") +pred = ROCR::prediction(xx$signal, xx$refractory) +modelresult1 = list(prob = xx$signal, features = NA, top = NA, + auc = ROCR::performance(pred, "auc")@y.values[[1]], pred = pred) + +pooled_t_test = stat_test_(xx %>% mutate(x = subcohort, y = signal)) +t_test_pval = sub(" \n.*", "", pooled_t_test) %>% sub(".* = ", "", .) %>% as.numeric() +## Boxplot with no annotations +q = ggplot(xx %>% filter(subcohort != "Paired_relapse_sample"), + aes(x = subcohort, y = signal, fill = subcohort)) + geom_boxplot(width = 0.6) + + ylab("Ven signal") + scale_fill_manual(values = subcohort_colors) + theme_bw() + + ggtitle("Exvivo Ven resistance signature") + + theme(text = element_text(size = 22), plot.title = element_text(hjust = 0.5)) + + +## Boxplot with t_test annotations +q1 = t_test_ann(q, c("Refractory", "Response\nno relapse")) +q2 = t_test_ann(q1, c("Relapsed", "Response\nno relapse")) +q2 = t_test_ann(q2, c("Refractory", "Relapsed"), y_start = -0.4) +q2 +ggsave(plot = q2, filename = paste0(plot_folder, "boxplot_exvivo_V_signature_applied_invivo.png"), + width = 8, height = 8) +ggsave(plot = q2, filename = paste0(plot_folder, "boxplot_exvivo_V_signature_applied_invivo.pdf"), + width = 8, height = 8) + + +p = plotAUC_gg(modelresult1) + theme_bw() + ggtitle("Detecting refractory patients \nExvivo Ven signature") + + theme(text = element_text(size = 18), plot.title = element_text(hjust = 0.5)) +ggsave(plot = p, filename = paste0(plot_folder, "aucplot_exvivo_V_signature_applied_invivo.png"), + width = 8, height = 8) +ggsave(plot = p, filename = paste0(plot_folder, "aucplot_exvivo_V_signature_applied_invivo.pdf"), + width = 8, height = 8) + + + +``` + + + +### Showing correlation with drug AUC in BeatAML + +```{r} +## Exvivo DV signature +beataml_210_predictions = predict_helper_(m_210, DV_markers_up, DV_markers_down) %>% + left_join(., drug_data_df, by = c("sample" = "sample_id")) %>% + left_join(., pData(m_210), by = c("sample" = "Barcode.ID")) + +plot_df = beataml_210_predictions %>% + filter(inhibitor == "Azacytidine - Venetoclax") +ll = plot_df %>% filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", + "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +a = cor.test(ll$auc, ll$signal) +a +line_slope = coef(lm(ll$signal ~ ll$auc))[[2]] +line_intercept = coef(lm(ll$signal ~ ll$auc))[[1]] +plot_title = paste0("Ex-vivo Dec+Ven signature vs Ven+Aza drug AUC\n", "correlation = ", formatC(a$estimate, 2), + ", p_value = ", formatC(a$p.value, 2)) +p = ggplot(ll, aes(x = auc, y = signal, label = sample)) + geom_point(size = 2.5) + + ggtitle(plot_title) + ylab("signal") + xlab("Ven+Aza drug AUC") + + geom_abline(intercept = line_intercept, slope = line_slope, size = 1.5, alpha = 0.2) + + theme_bw() + theme(text = element_text(size = 15)) +p +ggsave(plot = p, filename = paste0(plot_folder, "beataml_VEN_AZA_drugauc_validation_exvivo_DV_signature.pdf")) +ggsave(plot = p, filename = paste0(plot_folder, "beataml_VEN_AZA_drugauc_validation_exvivo_DV_signature.png")) + + +plot_df = beataml_210_predictions %>% + filter(inhibitor == "Venetoclax") +ll = plot_df %>% filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", + "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +a = cor.test(ll$auc, ll$signal) +a +line_slope = coef(lm(ll$signal ~ ll$auc))[[2]] +line_intercept = coef(lm(ll$signal ~ ll$auc))[[1]] +plot_title = paste0("Ex-vivo Dec+Ven signature vs Ven drug AUC\n", "correlation = ", formatC(a$estimate, 2), + ", p_value = ", formatC(a$p.value, 2)) +p = ggplot(ll, aes(x = auc, y = signal, label = sample)) + geom_point(size = 2.5) + + ggtitle(plot_title) + ylab("signal") + xlab("Ven drug AUC") + + geom_abline(intercept = line_intercept, slope = line_slope, size = 1.5, alpha = 0.2) + + theme_bw() + theme(text = element_text(size = 15)) +p +ggsave(plot = p, filename = paste0(plot_folder, "beataml_VEN_drugauc_validation_exvivo_DV_signature.pdf")) +ggsave(plot = p, filename = paste0(plot_folder, "beataml_VEN_drugauc_validation_exvivo_DV_signature.png")) + + +``` + + + +```{r} + +## Exvivo V signature +beataml_210_predictions = predict_helper_(m_210, V_markers_up, V_markers_down) %>% + left_join(., drug_data_df, by = c("sample" = "sample_id")) %>% + left_join(., pData(m_210), by = c("sample" = "Barcode.ID")) + +plot_df = beataml_210_predictions %>% + filter(inhibitor == "Azacytidine - Venetoclax") +ll = plot_df %>% filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", + "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +a = cor.test(ll$auc, ll$signal) +a +line_slope = coef(lm(ll$signal ~ ll$auc))[[2]] +line_intercept = coef(lm(ll$signal ~ ll$auc))[[1]] +plot_title = paste0("Ex-vivo Ven signature vs Ven+Aza drug AUC\n", "correlation = ", formatC(a$estimate, 2), + ", p_value = ", formatC(a$p.value, 2)) +p = ggplot(ll, aes(x = auc, y = signal, label = sample)) + geom_point(size = 2.5) + + ggtitle(plot_title) + ylab("signal") + xlab("Ven+Aza drug AUC") + + geom_abline(intercept = line_intercept, slope = line_slope, size = 1.5, alpha = 0.2) + + theme_bw() + theme(text = element_text(size = 15)) +p +ggsave(plot = p, filename = paste0(plot_folder, "beataml_VEN_AZA_drugauc_validation_exvivo_V_signature.pdf")) +ggsave(plot = p, filename = paste0(plot_folder, "beataml_VEN_AZA_drugauc_validation_exvivo_V_signature.png")) + + +plot_df = beataml_210_predictions %>% + filter(inhibitor == "Venetoclax") +ll = plot_df %>% filter(Specimen.access.group.concatenated %in% c("Initial Acute Leukemia Diagnosis", + "Initial Acute Leukemia Diagnosis; Post-Chemotherapy")) +a = cor.test(ll$auc, ll$signal) +a +line_slope = coef(lm(ll$signal ~ ll$auc))[[2]] +line_intercept = coef(lm(ll$signal ~ ll$auc))[[1]] +plot_title = paste0("Ex-vivo Ven signature vs Ven drug AUC\n", "correlation = ", formatC(a$estimate, 2), + ", p_value = ", formatC(a$p.value, 2)) +p = ggplot(ll, aes(x = auc, y = signal, label = sample)) + geom_point(size = 2.5) + + ggtitle(plot_title) + ylab("signal") + xlab("Ven drug AUC") + + geom_abline(intercept = line_intercept, slope = line_slope, size = 1.5, alpha = 0.2) + + theme_bw() + theme(text = element_text(size = 15)) +p +ggsave(plot = p, filename = paste0(plot_folder, "beataml_VEN_drugauc_validation_exvivo_V_signature.pdf")) +ggsave(plot = p, filename = paste0(plot_folder, "beataml_VEN_drugauc_validation_exvivo_V_signature.png")) + + +``` + + + + + diff --git a/ven_hma_samples/scripts/compute_lsc17.R b/ven_hma_samples/scripts/compute_lsc17.R new file mode 100644 index 0000000..4b618f6 --- /dev/null +++ b/ven_hma_samples/scripts/compute_lsc17.R @@ -0,0 +1,42 @@ +##### NOTE ########## + +## Obtained from https://github.com/biodev/beataml2_manuscript/blob/main/R/preprocess.R +## This is the function to compute lsc17 scores used in the Dan + Jeff publication https://pmc.ncbi.nlm.nih.gov/articles/PMC9378589/pdf/nihms-1822675.pdf +## paper name is: Integrative Analysis of Drug Response and Clinical Outcome in Acute Myeloid Leukemia +library(data.table) + +compute.lsc17 <- function(exprs){ + + #from https://www.nature.com/articles/nature20598#Tab8 + #(DNMT3B×0.0874) + (ZBTB46×−0.0347) + (NYNRIN×0.00865) + (ARHGAP22×−0.0138) + (LAPTM4B×0.00582) + (MMRN1×0.0258) + (DPYSL3×0.0284) + + #(KIAA0125×0.0196) + (CDK6×−0.0704) + (CPXM1×−0.0258) + (SOCS2×0.0271) + (SMIM24×−0.0226) + (EMP1×0.0146) + (NGFRAP1×0.0465) + + #(CD34×0.0338) + (AKR1C3×−0.0402) + (GPR56×0.0501). + #As above- and below-median scores in the training cohort were associated with adverse and favourable cytogenetic risk, respectively, a median threshold was used to discretize scores into high and low groups. + + lsc17.coefs <- c(DNMT3B=0.0874, ZBTB46=-0.0347, NYNRIN=0.00865, ARHGAP22=-0.0138, LAPTM4B=0.00582, MMRN1=0.0258, DPYSL3=0.0284, + KIAA0125=0.0196, CDK6=-0.0704, CPXM1=-0.0258, SOCS2=0.0271, SMIM24=-0.0226, EMP1=0.0146, NGFRAP1=0.0465, CD34=0.0338, + AKR1C3=-0.0402,GPR56=0.0501) + + common.genes <- intersect(names(lsc17.coefs), colnames(exprs))#16 + + if (length(common.genes) == 16){ + stopifnot(all(setdiff(names(lsc17.coefs), colnames(exprs)) == "SMIM24")) + + #smim24 has an ensembl id of ENSG00000095932 so maps to C19orf77 + + names(lsc17.coefs)[names(lsc17.coefs) == "SMIM24"] <- "C19orf77" + + common.genes <- intersect(names(lsc17.coefs), colnames(exprs))#17 + + stopifnot(length(common.genes) == 17) + } + + # lsc17.score <- exprs[,names(lsc17.coefs)] %*% lsc17.coefs + lsc17.score <- exprs[, common.genes] %*% lsc17.coefs[common.genes] + + #should compute median per cohort + lsc.dt <- data.table(ptid=rownames(lsc17.score), LSC17=lsc17.score[,1]) + + lsc.dt + +} diff --git a/ven_hma_samples/scripts/helper_scripts.R b/ven_hma_samples/scripts/helper_scripts.R new file mode 100644 index 0000000..cd2af9d --- /dev/null +++ b/ven_hma_samples/scripts/helper_scripts.R @@ -0,0 +1,386 @@ +library(MSnSet.utils) +library(dplyr) +library(KSEAapp) +library(ggplot2) +library(scales) +library(ggpubr) + + +## Modification of MSnSetUtils function. +plot_pca <- function(eset, phenotype = NULL, shape = NULL, label = NULL, z_score = TRUE, + princomp_center = TRUE, show_ellipse = TRUE, components = 1:2, biplot = FALSE, + biplot_labels = NULL, standardize = TRUE, save_dfs = NULL, + num_features = 6L, show_NA = TRUE, label_size = 3, output_type = NA, + legend_title = phenotype, + arrow_args = list(), label_args = list(), ...) { + + # Handling coloring by phenotype. Do this first, in case + # rows are removed when show_NA = FALSE + if (!is.null(phenotype)) { + colorBy <- pData(eset)[, phenotype] + # If not showing missing values, remove those samples + if (!show_NA) { + idx <- !is.na(colorBy) + eset <- eset[, idx] + colorBy <- colorBy[idx] + } + } else { + show_ellipse <- FALSE + colorBy <- NULL + } + if (!is.null(shape)){ + shapeBy <- pData(eset)[, shape] + if (!show_NA) { + idx <- !is.na(shapeBy) + eset <- eset[, idx] + shapeBy <- shapeBy[idx] + } + } else { + shapeBy <- NULL + } + + # Check that components are valid + if (length(components) != 2) { + stop(sprintf("components must be a vector of length 2, not %d.", + length(components))) + } + if (!all(components %in% 1:ncol(eset))) { + stop(sprintf("The values of components must be between 1 and %d.", + ncol(eset))) + } + + complete_rows <- complete.cases(exprs(eset)) + + # Check that there are enough complete rows for PCA + if (sum(complete_rows) < 2) { + stop("There are fewer than 2 rows with non-missing data.") + } + + message(sprintf("Subsetting to %d complete rows for PCA.", + sum(complete_rows))) + + # Subset to complete rows + eset <- eset[complete_rows, ] + + # If z_score, convert to Z-Scores by sample (row when transposed) + if (z_score) { + z <- t(scale(exprs(eset), center = TRUE, scale = TRUE)) + } else { + z <- t(exprs(eset)) + } + + ## PCA + # By default, center = TRUE, scale. = FALSE + pca_res <- prcomp(z, center = princomp_center) + + u <- pca_res$x # Scores + v <- pca_res$rotation # Eigenvectors + + if (standardize) { + n <- nrow(u) + lam <- pca_res$sdev * sqrt(n) + + # Scale u down and v up. Product is still the same + u <- t(t(u) / lam) + v <- t(t(v) * lam) + } + + # Determine ratio between scale of v and u + u_range <- apply(u[, components], 2, function(x) abs(range(x))) + v_range <- apply(v[, components], 2, function(x) abs(range(x))) + + ratio <- max(v_range / u_range) # ratio for scaling v and secondary axes + v <- v / ratio # scale v + + if (!is.null(save_dfs)){ + u <- u %>% as.data.frame() + # If colorBy is not NULL, add that column to df + if (!is.null(colorBy)) { + u$colorBy <- colorBy + } + if (!is.null(shapeBy)) { + u$shapeBy <- shapeBy + } + + assign(save_dfs, list("sample_decomposition" = u, "feature_decomposition" = v %>% as.data.frame() %>% + mutate(feature = rownames(.))), envir = globalenv()) + } + + # Data frames for plotting + df.u <- as.data.frame(u[, components]) + df.v <- as.data.frame(v[, components]) + + # Percent of variance explained by each PC + d <- pca_res$sdev # Standard deviations + var_expl <- round(100 * d ^ 2 / sum(d ^ 2), digits = 2)[components] + axis_labs <- sprintf("PC%d (%g%%)", #"%sPC%d (%g%%)", + # ifelse(obs.scale == 0, "Standardized ", ""), + components, + var_expl) + + # If colorBy is not NULL, add that column to df + if (!is.null(colorBy)) { + df.u$colorBy <- colorBy + } + if (!is.null(shapeBy)) { + df.u$shapeBy <- shapeBy + } + + ## Visualization + # Base plot + p <- ggplot(data = df.u, mapping = aes(x = df.u[, 1], y = df.u[, 2], color = colorBy, shape = shapeBy)) + + geom_hline(yintercept = 0, lty = "longdash", color = "darkgrey") + + geom_vline(xintercept = 0, lty = "longdash", color = "darkgrey") + + labs(x = axis_labs[1], y = axis_labs[2]) + + theme_bw() + + theme(aspect.ratio = 1) + + # 50% confidence ellipse layer first so they are + # beneath the layer of points or labels. + if (show_ellipse & !is.numeric(colorBy)) { + p <- p + + stat_ellipse(mapping = aes(fill = colorBy, color = NULL), + geom = "polygon", type = "norm", + level = 0.5, alpha = 0.1, show.legend = TRUE) + } + + # If label is NULL, add points. Otherwise, add labels + if (is.null(label)) { + p <- p + + geom_point(...) + } else { + labels <- pData(eset)[, label] + p <- p + geom_point(...) + + ggrepel::geom_label_repel(mapping = aes(label = labels), + size = label_size, ...) + } + + # Set titles for color and fill legend + p <- p + + guides(color = guide_legend(title = legend_title), + fill = guide_legend(title = legend_title)) + + # If colorBy is numeric, use a colorbar + if (is.numeric(colorBy)) { + p <- p + + guides(color = guide_colorbar(title = legend_title)) + } + + ## Biplot + if (biplot) { + # Get the indices of the top influential features + # from each principal component. num_features determines how + # many to select from each component. + top_features <- lapply(1:2, function(i) { + order(abs(df.v)[, i], decreasing = TRUE)[1:num_features] + }) + top_features <- unique(unlist(top_features)) + + # Subset loadings to top features and rename columns + df.v <- df.v[top_features, ] + colnames(df.v) <- c("xend", "yend") + df.v$x <- df.v$y <- 0 + + # If biplot_labels is not provided, default to row names + if (is.null(biplot_labels)) { + df.v$labels <- rownames(df.v) + } else { + df.v$labels <- fData(eset)[top_features, biplot_labels] + } + + scale_args <- list(expand = expansion(mult = rep(0.1, 2)), + sec.axis = sec_axis(~ . * ratio)) + + # Arguments for geom_segment + arrow_args <- list(mapping = aes(x = x, y = y, xend = xend, yend = yend), + arrow = arrow(length = unit(0.5, "line")), + data = df.v, color = "red3") %>% + # Allow user-supplied args to overwrite defaults + modifyList(val = arrow_args, keep.null = TRUE) + + # Arguments for geom_label_repel + label_args <- list(mapping = aes(x = xend, y = yend, label = labels), + data = df.v, + color = arrow_args[["color"]], + max.overlaps = Inf, + min.segment.length = 0, + fill = alpha("white", 0.5)) %>% + # Allow user-supplied args to overwrite defaults + modifyList(val = label_args, keep.null = TRUE) + + # Add segments with arrows and text labels + p <- p + + # Add extra padding around plot area and secondary axes for v units + do.call(scale_x_continuous, scale_args) + + do.call(scale_y_continuous, scale_args) + + do.call(geom_segment, arrow_args) + + do.call(geom_label_repel, label_args) + + theme(axis.text.y.right = element_text(color = arrow_args[["color"]]), + axis.text.x.top = element_text(color = arrow_args[["color"]]), + axis.ticks.y.right = element_line(color = arrow_args[["color"]]), + axis.ticks.x.top = element_line(color = arrow_args[["color"]])) + } + + if (output_type == "full"){ + return(list("plot" = p, "data.u" = as.data.frame(u[, ]), "data.v" = as.data.frame(v[, ]), + "axis_labs" = axis_labs)) + } else{ + return(p) + } +} + +diffexp_helper <- function(m, contrast_var, contrasts){ + pData(m)$bgd_ <- pData(m)[[contrast_var]] + pData(m)$Sample <- sampleNames(m) + + all_results <- data.frame() + + for (contrast in contrasts){ + contrast_groups = strsplit(contrast, "-")[[1]] + contrast = paste0("bgd_", contrast_groups[[1]], "-bgd_", contrast_groups[[2]]) + limma_res <- limma_contrasts(m, model.str = "~0 + bgd_", + coef.str = "bgd_", contrasts = contrast) %>% as.data.frame() + rownames(limma_res) <- limma_res$feature + counter = 1 + + m_contrast <- m[, m$bgd_ %in% contrast_groups] + p_values_t_test <- vector(mode="character", length = nrow(limma_res)) + p_values_welch_test <- vector(mode="character", length = nrow(limma_res)) + p_values_wilcox_test <- vector(mode="character", length = nrow(limma_res)) + for (feature in limma_res$feature){ + data_df <- data.frame(value = exprs(m_contrast)[feature, ], + Sample = colnames(exprs(m_contrast))) %>% + filter(!is.na(value)) %>% + merge(pData(m_contrast) %>% select(Sample, bgd_), by = "Sample") + + p_values_t_test[[counter]] <- tryCatch({t.test(value ~ bgd_, data = data_df, + alternative = "two.sided", var.equal = TRUE)[[3]]}, + error = function(e) {NA}) + p_values_welch_test[[counter]] <- tryCatch({t.test(value ~ bgd_, data = data_df, + alternative = "two.sided", var.equal = FALSE)[[3]]}, + error = function(e) {NA}) + p_values_wilcox_test[[counter]] <- tryCatch({wilcox.test(value ~ bgd_, data = data_df, + alternative = "two.sided")[[3]]}, + error = function(e) {NA}) + counter = counter + 1 + } + + limma_res <- limma_res %>% + mutate(t_test_pval = as.numeric(p_values_t_test), + t_test_adj = p.adjust(t_test_pval, method = "BH"), + welch_pval = as.numeric(p_values_welch_test), + welch_adj = p.adjust(welch_pval, method = "BH"), + wilcox_pval = as.numeric(p_values_wilcox_test), + wilcox_adj = p.adjust(wilcox_pval, method = "BH")) + all_results <- rbind(all_results, limma_res) + } + + return(all_results) +} + + +plot_features <- function (m, features, feature_name_col = NULL, color_by = NULL, mode = "box", + order_by = color_by){ + p_data <- pData(m) %>% select(-matches("sample name")) + p_data[["sample name"]] <- rownames(p_data) + p_data <- p_data[order(p_data[[order_by]]), ] + p_data[['sample name']] <- factor(p_data[['sample name']], levels = p_data[['sample name']]) + if (is.null(feature_name_col)) + feature_names <- featureNames(m) + else feature_names <- fData(m)[[feature_name_col]] + + idx <- c() + # all(grepl("^.*-[STY]{1}[0-9]+.*$", feature_names)) finds if all features in msnset are phosphosites + # !any(grepl("^.*-[STY]{1}[0-9]+.*$", features)) finds if NONE of the chosen features provided are phosphosites + # when both true, we plot the phosphosites from feature_names belonging to the chosen features (proteins). + if (all(grepl("^.*-[STY]{1}[0-9]+.*$", feature_names)) & !any(grepl("^.*-[STY]{1}[0-9]+.*$", features))){ + print("detected phospho data, collecting phosphosites of the supplied features") + for (feature in features){ + idx <- c(idx, which(startsWith(feature_names, feature))) + } + } else{ + idx <- which(featureNames(m) %in% features) + } + + x <- exprs(m[idx, ]) %>% as.data.frame() %>% tibble::rownames_to_column("feature") %>% + tidyr::pivot_longer(-feature, names_to = "sample name", values_to = "abundance") %>% + inner_join(p_data, by = "sample name") %>% + mutate(`sample name` = factor(`sample name`, levels = p_data[['sample name']])) + if (length(features) > 1){ + text_element <- element_blank() + } else { + text_element <- element_text(angle = 45, + hjust = 1) + } + if (mode == "point"){ + p <- x %>% ggplot() + aes(x = `sample name`, y = abundance) + + geom_point(size = 3) + theme_bw() + theme(axis.text.x = text_element) + + facet_wrap(~ feature, scales = "free") + if (!is.null(color_by)) + p <- p + aes_string(color = color_by) + p + } else if (mode == "box") { + p <- x %>% ggplot() + aes(x = color_by, y = abundance) + + geom_boxplot(notch = TRUE) + theme_bw() + theme(axis.text.x = text_element) + + facet_wrap(~ feature, scales = "free") + if (!is.null(color_by)) + p <- p + aes_string(color = color_by) + p + } + +} + + +volcano_function <- function(diffexp_res, chosen_terms, t2g, label_genes = c("Tkfc")){ + sig_genes <- diffexp_res %>% filter(t_test_adj < 0.05) %>% pull(feature) %>% unique() + t2g_annotation <- t2g %>% + filter(gene_symbol %in% sig_genes, + gs_name %in% chosen_terms) %>% + group_by(gs_name) %>% + mutate(total = n()) %>% + ungroup() %>% + group_by(gene_symbol) %>% + top_n(n = 1, wt = -total) %>% + slice(1) %>% + select(feature = gene_symbol, pathway = gs_name) + + ## Remove later + t2g_annotation$pathway <- gsub("_", " ", t2g_annotation$pathway) + t2g_annotation$pathway <- gsub("HALLMARK ", "", t2g_annotation$pathway) + ### + + gg_color_hue <- function(n) { + hues = seq(15, 375, length = n + 1) + hcl(h = hues, l = 65, c = 100)[1:n] + } + + diffexp_res_annotated <- left_join(diffexp_res, t2g_annotation, by = "feature") %>% + mutate(pathway = case_when(!is.na(pathway) ~ pathway, + TRUE ~ 'NA')) %>% + mutate(pathway = factor(pathway, levels = c(t2g_annotation$pathway %>% unique(), 'NA')), + size = case_when(pathway != 'NA' ~ 4, + TRUE ~ 1.3)) %>% + mutate(feat_lab = case_when(pathway != 'NA' ~ feature, + TRUE ~ '')) %>% + arrange(t_test_adj) + diffexp_res_annotated$rank <- 1:nrow(diffexp_res_annotated) + + diffexp_res_annotated$alpha[diffexp_res_annotated$feat_lab == ''] <- 0.85 + diffexp_res_annotated$alpha[diffexp_res_annotated$feat_lab != ''] <- 1 + diffexp_res_annotated$feat_lab[diffexp_res_annotated$feat_lab == ''] <- NA + diffexp_res_annotated$feat_lab[!(diffexp_res_annotated$feature %in% label_genes)] <- NA + + + p <- ggplot(diffexp_res_annotated, aes(x = logFC, y = -log10(t_test_adj), + color = pathway, alpha = alpha, size = size)) + geom_point() #+ + # scale_color_manual(values = c(gg_color_hue(length(unique(t2g_annotation$pathway))), "grey30"), + # name = "Pathway", limits = unique(t2g_annotation$pathway)) #+ + ## Add again later + # geom_label_repel(aes(label = feat_lab), size = 10, force = 5) + + # scale_size_identity() + scale_alpha(guide = 'none') + + # theme(# axis.text = element_text(size = 20), legend.text = element_text(size = 20), + # text = element_text(size = 17)) + + return(list(p, diffexp_res_annotated)) +} + diff --git a/ven_hma_samples/scripts/loading_data.R b/ven_hma_samples/scripts/loading_data.R new file mode 100644 index 0000000..d1e9f0c --- /dev/null +++ b/ven_hma_samples/scripts/loading_data.R @@ -0,0 +1,114 @@ +library(dplyr) +library(MSnSet.utils) +library(ggplot2) +library(ggpubr) +library(patchwork) + +syn <- synapseLogin() +subcohort_colors = c("#e8991b", "#00798c", "#d1495b", "grey27") +subcohort_colors2 = c("#e8991b", "#00798c", "#d1495b", "grey27") +names(subcohort_colors) = c("Response\nno relapse", "Refractory", "Relapsed", "Paired_relapse_sample") +vital_colors = c('orange3', 'plum4', 'darkgrey') +names(vital_colors) = c("Alive", "Dead", "LTFU") + +## Load exp28 +mat_path = syn$get("syn71717828")$path +xx = read.csv2(mat_path, sep = "\t") +# xx_ = xx %>% mutate(N_missing = rowSums(xx[, c(7:47)] == "")) %>% +# group_by(Genes) %>% mutate(min_ = min(N_missing), max_ = max(N_missing)) %>% +# ungroup() %>% filter(N_missing == min_) +# yy = read.csv2("../Aggregation_report.pg_matrix.tsv", sep = "\t") + +mat_df = xx[, c(3, 7:47)] %>% as.data.frame() ## Take Gene name and samples +colnames(mat_df) = sub("^.*WorkDir1\\.", "", colnames(mat_df)) +colnames(mat_df) = sub("_FAIMS_.*$", "", colnames(mat_df)) +colnames(mat_df) = sub("Astral_", "", colnames(mat_df)) +colnames(mat_df) = sub("_03Nov25_Monty_ES906_3258_50SPD.mzML", "", colnames(mat_df)) +colnames(mat_df) +## Rolling up, since there are duplicate gene names. Taking log2 then mean to rollup. +mat_df = mat_df %>% + tidyr::pivot_longer(-Genes, names_to = "sample_name", values_to = "value") %>% + mutate(value = as.numeric(value), + value = log2(value)) %>% + group_by(Genes, sample_name) %>% + mutate(value_rollup = mean(value)) %>% ungroup() %>% + select(Genes, sample_name, value_rollup) %>% unique() %>% + tidyr::pivot_wider(names_from = "sample_name", values_from = "value_rollup") +mat_df = mat_df %>% filter(Genes != "") +mat = mat_df[, -1] %>% as.matrix() +rownames(mat) = mat_df$Genes + + +meta_path = syn$get("syn69058992")$path +meta = openxlsx::read.xlsx(meta_path) +clinical_summary = openxlsx::read.xlsx(meta_path, sheet = "clinical_summary") %>% + select(labId, ageAtDiagnosis, priorMalignancyType, priorMalignancyChemo, + elnRisk2022, karyotype, FISH, immunophenotype, geneMutationSummary, + dxAtSpecimenAcquisition, specificDxAtAcquisition, specimenGroups, specimenType, + cumulativeTreatmentTypes, cumulativeTreatmentRegimenCount, cumulativeTreatmentRegimens, + percentBlastsBM, percentBlastsPB) +meta = left_join(meta, clinical_summary, by = "labId") +meta$survival_days = meta[[13]] +p_data = meta %>% mutate(sample_name = datasetNameGlobal, + subcohort = gsub(" - ", "_", `sub-cohortPNNL`), + `sub-cohort.PNNL` = subcohort, + subcohort = gsub(" ", "_", subcohort), + group = subcohort, + group = factor(group, levels = c("Response_no_relapse", "Refractory", "Relapse", "Paired_relapse_sample")), + group_surv = case_when(subcohort == "Response_no_relapse" & survival_days >= 500 ~ "response_long", + subcohort == "Response_no_relapse" & survival_days < 500 ~ "response_short", + subcohort == "Refractory" ~ "Refractory", + subcohort == "Relapse" & survival_days < 500 ~ "relapse_short", + subcohort == "Relapse" & survival_days >= 500 ~ "relapse_long", + TRUE ~ subcohort)) +rownames(p_data) = p_data$sample_name +p_data = p_data[colnames(mat), ] + + +## This is the Ven + HMA dataset (41 samples) +m_exp28 = MSnSet(exprs = mat, pData = p_data) +# hist(apply(exprs(m_exp28), 2, mean, na.rm = T)) +# hist(apply(exprs(m_exp28), 1, mean, na.rm = T)) +exprs(m_exp28) = sweep(exprs(m_exp28), 2, apply(exprs(m_exp28), 2, mean, na.rm = T), FUN = '-') + + + + +## Load corrected exp25 dataset, subset to 210 features? + +## This is the ex vivo dataset +m_4pat <- readRDS(syn$get("syn64605135")$path) +m_210 <- readRDS(syn$get("syn72644178")$path) +m_210 <- m_210[complete.cases(exprs(m_210)), ] +m_4pat <- m_4pat[featureNames(m_4pat) %in% featureNames(m_210), ] +drug_data <- read.csv(syn$get("syn51674470")$path) + +# zz = read.csv("../210_cohort_drug_auc_syn51674470.csv") +# table(zz$inhibitor) + +## Inhibitors of interest +drug_data_df <- rbind(drug_data %>% filter(inhibitor == "Venetoclax"), ## 127 samples) + drug_data %>% filter(inhibitor == "Azacytidine - Venetoclax"), ## 20 samples + drug_data %>% filter(inhibitor == "Bortezomib - Venetoclax"), ## 19 samples + drug_data %>% filter(inhibitor == "Dasatinib - Venetoclax"), ## 92 samples + drug_data %>% filter(inhibitor == "Doramapimod - Venetoclax"), ## 110 samples + drug_data %>% filter(inhibitor == "GW-2580 - Venetoclax"), ## 20 samples + drug_data %>% filter(inhibitor == "Idelalisib - Venetoclax"), ## 107 samples + drug_data %>% filter(inhibitor == "Olaparib - Venetoclax"), ## 18 samples + drug_data %>% filter(inhibitor == "Quizartinib - Venetoclax"), ## 80 samples + drug_data %>% filter(inhibitor == "Ruxolitinib - Venetoclax"), ## 111 samples + drug_data %>% filter(inhibitor == "Sorafenib - Venetoclax"), ## 92 samples + drug_data %>% filter(inhibitor == "Trametinib - Venetoclax"), ## 84 samples + drug_data %>% filter(inhibitor == "Venetoclax - Artemisinin"), ## 76 samples + drug_data %>% filter(inhibitor == "Venetoclax - Ibrutinib"), ## 119 samples + drug_data %>% filter(inhibitor == "Venetoclax - JQ1"), ## 85 samples + drug_data %>% filter(inhibitor == "Venetoclax - Palbociclib"), ## 85 samples + drug_data %>% filter(inhibitor == "Venetoclax - Panobinostat"), ## 82 samples + drug_data %>% filter(inhibitor == "Venetoclax - PH797804") ## 11 samples +) + + + + + + diff --git a/ven_hma_samples/scripts/plot_helper.R b/ven_hma_samples/scripts/plot_helper.R new file mode 100644 index 0000000..e7d06a8 --- /dev/null +++ b/ven_hma_samples/scripts/plot_helper.R @@ -0,0 +1,92 @@ +### Function to annotate t_test pvalues in a boxplot (y = value, x = category) comparing two values in the x-axis. +t_test_ann <- function(q, comparison, pval_cutoff = 1.1, y_start = NULL, test = "t_test"){ + x_var = rlang::quo_get_expr(q$mapping$x) + y_var = rlang::quo_get_expr(q$mapping$y) + test_df = q$data %>% dplyr::select(x_var, y_var) %>% + mutate(y = !!y_var, x = !!x_var) + if (!(is.factor(test_df$x))){ + test_df <- test_df %>% + mutate(x = factor(x)) + } + boxplot_df = layer_data(q, 1) + factor_levels = sapply(comparison, function(x){which(levels(test_df$x) == x)}) + factor_levels_inbetween = min(factor_levels):max(factor_levels) + x_start = min(factor_levels) - 0.05 + x_end = max(factor_levels) + 0.05 + seg_y = 0.04 * (max(boxplot_df$ymax) - min(boxplot_df$ymin)) + if (is.null(y_start)){ + y_start = max(boxplot_df[factor_levels_inbetween, "upper"]) + seg_y/2 + } + y_end = y_start + seg_y + + test_df = test_df %>% filter(!!x_var %in% comparison) + t_test = t.test(y ~ x, data = test_df, + alternative = "two.sided", var.equal = TRUE) + w_test = wilcox.test(y ~ x, data = test_df, + alternative = "two.sided", var.equal = TRUE) + if (test == "t_test"){ + stat_test = t_test + } else if (test == "wilcox_test"){ + stat_test = w_test + } + if (stat_test[[3]] < pval_cutoff){ + q = q + annotate("segment", x = x_start, xend = x_end, y = y_end, yend = y_end, linewidth = 0.8) + + annotate("segment", x = x_start, xend = x_start, y = y_start, yend = y_end, linewidth = 0.8) + + annotate("segment", x = x_end, xend = x_end, y = y_start, yend = y_end, linewidth = 0.8) + + annotate("text", x = (x_start + x_end)/2, y = y_end + seg_y, label = formatC(stat_test[[3]], digits = 2), + size = 4.3) + + } + return(q) +} + +# q1 = t_test_ann(q, c("Refractory", "Response_no_relapse")) +# q2 = t_test_ann(q1, c("Relapse", "Response_no_relapse")) +# q2 + +## Helper function to compute t_test between the combined group 'Repsonse_no_relapse + Relapse' vs 'Refractory' +stat_test_ <- function(plot_df){ + test_df = plot_df %>% + mutate(x_ = case_when(x == "Response\nno relapse" ~ "combined_response_group", + x == "Relapsed" ~ "combined_response_group", + TRUE ~ x)) %>% + filter(x_ != "Paired_relapse_sample") + t_test = t.test(y ~ x_, data = test_df, + alternative = "two.sided", var.equal = TRUE) + w_test = wilcox.test(y ~ x_, data = test_df, + alternative = "two.sided", var.equal = TRUE) + message_ = paste("t_test_pval =", formatC(t_test[[3]], digits = 2), "\n") %>% + paste("wilcox_test_pval =", formatC(w_test[[3]], digits = 2)) + print(message_) +} + +## Plot CD14 and CD34 +## Example of how to use the t_test_ann function +# plot_df <- exprs(m_exp28[c("CD14", "CD34"), ]) %>% as.data.frame() %>% tibble::rownames_to_column("feature") %>% +# tidyr::pivot_longer(-feature, names_to = "sample_name", values_to = "abundance") %>% +# inner_join(p_data, by = "sample_name") +# +# q = ggplot(plot_df %>% filter(feature == "CD14"), +# aes(x = subcohort, y = abundance, fill = group)) + geom_boxplot(width = 0.75) + +# ylab("Ven signal") + scale_fill_manual(values = subcohort_colors) + +# ggtitle("CD14") + +# theme(text = element_text(size = 15)) +# q +# +# q = ggplot(plot_df %>% filter(feature == "CD34"), +# aes(x = subcohort, y = abundance, fill = group)) + geom_boxplot(width = 0.75) + +# ylab("Ven signal") + scale_fill_manual(values = subcohort_colors) + +# ggtitle("CD34") + +# theme(text = element_text(size = 15)) +# q +# +# stat_test_(plot_df %>% mutate(x = subcohort, y = abundance)) +# ## Example of how to use the t_test_ann function +# q1 = t_test_ann(q, c("Refractory", "Response_no_relapse")) +# q2 = t_test_ann(q1, c("Relapse", "Response_no_relapse"), y_start = 3.5) +# q2 = t_test_ann(q2, c("Refractory", "Relapse"), y_start = 4) +# q2 +# q1 = t_test_ann(q, c("Refractory", "Response_no_relapse"), test = "wilcox_test") +# q2 = t_test_ann(q1, c("Relapse", "Response_no_relapse"), test = "wilcox_test") +# q2 = t_test_ann(q2, c("Refractory", "Relapse"), test = "wilcox_test", y_start = -0.12) +# q2 \ No newline at end of file diff --git a/ven_hma_samples/scripts/predict_helper.R b/ven_hma_samples/scripts/predict_helper.R new file mode 100644 index 0000000..9a26ad7 --- /dev/null +++ b/ven_hma_samples/scripts/predict_helper.R @@ -0,0 +1,62 @@ + + + +## Function to compute signal value using the markers contained in 'full_model'. This is an element in the output produced by +## the 'logistic_model' function, in the script 's2n_model' which is loaded at the top. +## NOTE: The markers in the 'full_model' objects are intersected with the set of COMPLETE (NO NA) features from the new data. +predict_helper <- function(new_msnset, full_model, case_group, control_group){ + case_marker = full_model$features %>% filter(group_ == case_group) %>% pull(feature) + control_marker = full_model$features %>% filter(group_ == control_group) %>% pull(feature) + + ## This ensures we don't have any NA when computing the signal in the new samples. + case_marker = intersect(case_marker, complete_features) %>% intersect(., featureNames(new_msnset)) + control_marker = intersect(control_marker, complete_features) %>% intersect(., featureNames(new_msnset)) + + case_mat = t(exprs(new_msnset)[case_marker, , drop = FALSE]) %>% as.data.frame() + control_mat = t(exprs(new_msnset)[control_marker, , drop = FALSE]) %>% as.data.frame() + if (ncol(control_mat) > 0){ + control_mat$control_avg = rowMeans(control_mat, na.rm = T) + control_mat <- control_mat %>% select(control_avg) + } + if (ncol(case_mat) > 0){ + case_mat$case_avg = rowMeans(case_mat, na.rm = T) + case_mat <- case_mat %>% select(case_avg) + } + + logistic_df <- data.frame(sample = sampleNames(new_msnset)) %>% + cbind(case_mat, control_mat) + logistic_pred = predict.glm(full_model$model, logistic_df, type = "response") + coefs = full_model$model$coefficients + logistic_df <- logistic_df %>% + mutate(logistic_signal = coefs[['(Intercept)']] + coefs[['case_avg']]*case_avg + coefs[['control_avg']]*control_avg, + logistic_signal2 = case_avg - control_avg, + prob = logistic_pred, + prob2 = 1/(1+exp(-logistic_signal))) + return(logistic_df) +} + + + +## Function to compute signal value using two sets of markers. The signal will be mean(case_markers) - mean(control_markers) +## NOTE: The markers are intersected with the set of COMPLETE (NO NA) features from the new data. +predict_helper_ <- function(new_msnset, case_markers, control_markers){ + complete_features = featureNames(new_msnset)[(rowSums(is.na(exprs(new_msnset))) == 0)] + case_markers = intersect(case_markers, complete_features) + control_markers = intersect(control_markers, complete_features) + case_mat = t(exprs(new_msnset)[case_markers, , drop = FALSE]) %>% as.data.frame() + control_mat = t(exprs(new_msnset)[control_markers, , drop = FALSE]) %>% as.data.frame() + if (ncol(control_mat) > 0){ + control_mat$control_avg = rowMeans(control_mat, na.rm = T) + control_mat <- control_mat %>% select(control_avg) + } + if (ncol(case_mat) > 0){ + case_mat$case_avg = rowMeans(case_mat, na.rm = T) + case_mat <- case_mat %>% select(case_avg) + } + logistic_df <- data.frame(sample = sampleNames(new_msnset)) %>% + cbind(case_mat, control_mat) %>% + mutate(signal = case_avg - control_avg) + return(logistic_df) +} + + diff --git a/ven_hma_samples/scripts/s2n_model.R b/ven_hma_samples/scripts/s2n_model.R index e955060..0ae94ee 100644 --- a/ven_hma_samples/scripts/s2n_model.R +++ b/ven_hma_samples/scripts/s2n_model.R @@ -274,47 +274,271 @@ eln_model <- function(msnset, response, pred.cls, alpha){ } +## Modification of MSnSetUtils function. +plot_pca <- function(eset, phenotype = NULL, shape = NULL, label = NULL, z_score = TRUE, + princomp_center = TRUE, show_ellipse = TRUE, components = 1:2, biplot = FALSE, + biplot_labels = NULL, standardize = TRUE, save_dfs = NULL, + num_features = 6L, show_NA = TRUE, label_size = 3, output_type = NA, + legend_title = phenotype, + arrow_args = list(), label_args = list(), ...) { + + # Handling coloring by phenotype. Do this first, in case + # rows are removed when show_NA = FALSE + if (!is.null(phenotype)) { + colorBy <- pData(eset)[, phenotype] + # If not showing missing values, remove those samples + if (!show_NA) { + idx <- !is.na(colorBy) + eset <- eset[, idx] + colorBy <- colorBy[idx] + } + } else { + show_ellipse <- FALSE + colorBy <- NULL + } + if (!is.null(shape)){ + shapeBy <- pData(eset)[, shape] + if (!show_NA) { + idx <- !is.na(shapeBy) + eset <- eset[, idx] + shapeBy <- shapeBy[idx] + } + } else { + shapeBy <- NULL + } + + # Check that components are valid + if (length(components) != 2) { + stop(sprintf("components must be a vector of length 2, not %d.", + length(components))) + } + if (!all(components %in% 1:ncol(eset))) { + stop(sprintf("The values of components must be between 1 and %d.", + ncol(eset))) + } + + complete_rows <- complete.cases(exprs(eset)) + + # Check that there are enough complete rows for PCA + if (sum(complete_rows) < 2) { + stop("There are fewer than 2 rows with non-missing data.") + } + + message(sprintf("Subsetting to %d complete rows for PCA.", + sum(complete_rows))) + + # Subset to complete rows + eset <- eset[complete_rows, ] + + # If z_score, convert to Z-Scores by sample (row when transposed) + if (z_score) { + z <- t(scale(exprs(eset), center = TRUE, scale = TRUE)) + } else { + z <- t(exprs(eset)) + } + + ## PCA + # By default, center = TRUE, scale. = FALSE + pca_res <- prcomp(z, center = princomp_center) + + u <- pca_res$x # Scores + v <- pca_res$rotation # Eigenvectors + + if (standardize) { + n <- nrow(u) + lam <- pca_res$sdev * sqrt(n) + + # Scale u down and v up. Product is still the same + u <- t(t(u) / lam) + v <- t(t(v) * lam) + } + + # Determine ratio between scale of v and u + u_range <- apply(u[, components], 2, function(x) abs(range(x))) + v_range <- apply(v[, components], 2, function(x) abs(range(x))) + + ratio <- max(v_range / u_range) # ratio for scaling v and secondary axes + v <- v / ratio # scale v + + if (!is.null(save_dfs)){ + u <- u %>% as.data.frame() + # If colorBy is not NULL, add that column to df + if (!is.null(colorBy)) { + u$colorBy <- colorBy + } + if (!is.null(shapeBy)) { + u$shapeBy <- shapeBy + } + + assign(save_dfs, list("sample_decomposition" = u, "feature_decomposition" = v %>% as.data.frame() %>% + mutate(feature = rownames(.))), envir = globalenv()) + } + + # Data frames for plotting + df.u <- as.data.frame(u[, components]) + df.v <- as.data.frame(v[, components]) + + # Percent of variance explained by each PC + d <- pca_res$sdev # Standard deviations + var_expl <- round(100 * d ^ 2 / sum(d ^ 2), digits = 2)[components] + axis_labs <- sprintf("PC%d (%g%%)", #"%sPC%d (%g%%)", + # ifelse(obs.scale == 0, "Standardized ", ""), + components, + var_expl) + + # If colorBy is not NULL, add that column to df + if (!is.null(colorBy)) { + df.u$colorBy <- colorBy + } + if (!is.null(shapeBy)) { + df.u$shapeBy <- shapeBy + } + + ## Visualization + # Base plot + p <- ggplot(data = df.u, mapping = aes(x = df.u[, 1], y = df.u[, 2], color = colorBy, shape = shapeBy)) + + geom_hline(yintercept = 0, lty = "longdash", color = "darkgrey") + + geom_vline(xintercept = 0, lty = "longdash", color = "darkgrey") + + labs(x = axis_labs[1], y = axis_labs[2]) + + theme_bw() + + theme(aspect.ratio = 1) + + # 50% confidence ellipse layer first so they are + # beneath the layer of points or labels. + if (show_ellipse & !is.numeric(colorBy)) { + p <- p + + stat_ellipse(mapping = aes(fill = colorBy, color = NULL), + geom = "polygon", type = "norm", + level = 0.5, alpha = 0.1, show.legend = TRUE) + } + + # If label is NULL, add points. Otherwise, add labels + if (is.null(label)) { + p <- p + + geom_point(...) + } else { + labels <- pData(eset)[, label] + p <- p + geom_point(...) + + ggrepel::geom_label_repel(mapping = aes(label = labels), + size = label_size, ...) + } + + # Set titles for color and fill legend + p <- p + + guides(color = guide_legend(title = legend_title), + fill = guide_legend(title = legend_title)) + + # If colorBy is numeric, use a colorbar + if (is.numeric(colorBy)) { + p <- p + + guides(color = guide_colorbar(title = legend_title)) + } + + ## Biplot + if (biplot) { + # Get the indices of the top influential features + # from each principal component. num_features determines how + # many to select from each component. + top_features <- lapply(1:2, function(i) { + order(abs(df.v)[, i], decreasing = TRUE)[1:num_features] + }) + top_features <- unique(unlist(top_features)) + + # Subset loadings to top features and rename columns + df.v <- df.v[top_features, ] + colnames(df.v) <- c("xend", "yend") + df.v$x <- df.v$y <- 0 + + # If biplot_labels is not provided, default to row names + if (is.null(biplot_labels)) { + df.v$labels <- rownames(df.v) + } else { + df.v$labels <- fData(eset)[top_features, biplot_labels] + } + + scale_args <- list(expand = expansion(mult = rep(0.1, 2)), + sec.axis = sec_axis(~ . * ratio)) + + # Arguments for geom_segment + arrow_args <- list(mapping = aes(x = x, y = y, xend = xend, yend = yend), + arrow = arrow(length = unit(0.5, "line")), + data = df.v, color = "red3") %>% + # Allow user-supplied args to overwrite defaults + modifyList(val = arrow_args, keep.null = TRUE) + + # Arguments for geom_label_repel + label_args <- list(mapping = aes(x = xend, y = yend, label = labels), + data = df.v, + color = arrow_args[["color"]], + max.overlaps = Inf, + min.segment.length = 0, + fill = alpha("white", 0.5)) %>% + # Allow user-supplied args to overwrite defaults + modifyList(val = label_args, keep.null = TRUE) + + # Add segments with arrows and text labels + p <- p + + # Add extra padding around plot area and secondary axes for v units + do.call(scale_x_continuous, scale_args) + + do.call(scale_y_continuous, scale_args) + + do.call(geom_segment, arrow_args) + + do.call(geom_label_repel, label_args) + + theme(axis.text.y.right = element_text(color = arrow_args[["color"]]), + axis.text.x.top = element_text(color = arrow_args[["color"]]), + axis.ticks.y.right = element_line(color = arrow_args[["color"]]), + axis.ticks.x.top = element_line(color = arrow_args[["color"]])) + } + + if (output_type == "full"){ + return(list("plot" = p, "data.u" = as.data.frame(u[, ]), "data.v" = as.data.frame(v[, ]), + "axis_labs" = axis_labs)) + } else{ + return(p) + } +} + + ###SG: added this from ../../drug_treated_samples diffexp_helper <- function(m, contrast_var, contrasts){ - pData(m)$bgd_ <- pData(m)[[contrast_var]] - pData(m)$Sample <- sampleNames(m) - - all_results <- data.frame() - - for (contrast in contrasts){ - contrast_groups = strsplit(contrast, "-")[[1]] - contrast = paste0("bgd_", contrast_groups[[1]], "-bgd_", contrast_groups[[2]]) - limma_res <- limma_contrasts(m, model.str = "~0 + bgd_", - coef.str = "bgd_", contrasts = contrast) %>% as.data.frame() - rownames(limma_res) <- limma_res$feature - counter = 1 - - m_contrast <- m[, m$bgd_ %in% contrast_groups] - p_values_t_test <- vector(mode="character", length = nrow(limma_res)) - p_values_welch_test <- vector(mode="character", length = nrow(limma_res)) - for (feature in limma_res$feature){ - data_df <- data.frame(value = exprs(m_contrast)[feature, ], - Sample = colnames(exprs(m_contrast))) %>% - filter(!is.na(value)) %>% - merge(pData(m_contrast) %>% select(Sample, bgd_), by = "Sample") + pData(m)$bgd_ <- pData(m)[[contrast_var]] + pData(m)$Sample <- sampleNames(m) + + all_results <- data.frame() + + for (contrast in contrasts){ + contrast_groups = strsplit(contrast, "-")[[1]] + contrast = paste0("bgd_", contrast_groups[[1]], "-bgd_", contrast_groups[[2]]) + limma_res <- limma_contrasts(m, model.str = "~0 + bgd_", + coef.str = "bgd_", contrasts = contrast) %>% as.data.frame() + rownames(limma_res) <- limma_res$feature + counter = 1 + + m_contrast <- m[, m$bgd_ %in% contrast_groups] + p_values_t_test <- vector(mode="character", length = nrow(limma_res)) + p_values_welch_test <- vector(mode="character", length = nrow(limma_res)) + for (feature in limma_res$feature){ + data_df <- data.frame(value = exprs(m_contrast)[feature, ], + Sample = colnames(exprs(m_contrast))) %>% + filter(!is.na(value)) %>% + merge(pData(m_contrast) %>% select(Sample, bgd_), by = "Sample") + + p_values_t_test[[counter]] <- tryCatch({t.test(value ~ bgd_, data = data_df, + alternative = "two.sided", var.equal = TRUE)[[3]]}, + error = function(e) {NA}) + p_values_welch_test[[counter]] <- tryCatch({t.test(value ~ bgd_, data = data_df, + alternative = "two.sided", var.equal = FALSE)[[3]]}, + error = function(e) {NA}) + counter = counter + 1 + } - p_values_t_test[[counter]] <- tryCatch({t.test(value ~ bgd_, data = data_df, - alternative = "two.sided", var.equal = TRUE)[[3]]}, - error = function(e) {NA}) - p_values_welch_test[[counter]] <- tryCatch({t.test(value ~ bgd_, data = data_df, - alternative = "two.sided", var.equal = FALSE)[[3]]}, - error = function(e) {NA}) - counter = counter + 1 - } - - limma_res <- limma_res %>% - mutate(t_test_pval = as.numeric(p_values_t_test), - t_test_adj = p.adjust(t_test_pval, method = "BH"), - welch_pval = as.numeric(p_values_welch_test), - welch_adj = p.adjust(welch_pval, method = "BH")) - all_results <- rbind(all_results, limma_res) - } - - return(all_results) + limma_res <- limma_res %>% + mutate(t_test_pval = as.numeric(p_values_t_test), + t_test_adj = p.adjust(t_test_pval, method = "BH"), + welch_pval = as.numeric(p_values_welch_test), + welch_adj = p.adjust(welch_pval, method = "BH")) + all_results <- rbind(all_results, limma_res) + } + + return(all_results) }