add UMAP script

analyses/molecular-subtyping-MB/06-mb-shh-umap.Rmd

@@ -0,0 +1,351 @@
+---
+title: 'Create MB SHH methylation UMAP'
+output: 
+  html_document:
+  toc: TRUE
+toc_float: TRUE
+author: Ryan Corbett
+date: "2024"
+---
+
+Load libraries and set directory paths
+```{r}
+suppressPackageStartupMessages({
+  library(tidyverse)
+  library(umap)
+  library(ggplot2)
+  library(devtools)
+  library(gdata)
+})
+
+root_dir <- rprojroot::find_root(rprojroot::has_dir(".git"))
+
+data_dir <- file.path(root_dir, "data")
+analysis_dir <- file.path(root_dir, "analyses", "molecular-subtyping-MB")
+results_dir <- file.path(analysis_dir, "results")
+input_dir <- file.path(analysis_dir, "input")
+plots_dir <- file.path(analysis_dir, "plot")
+```
+
+Set file paths
+```{r}
+hist_file <- file.path(data_dir, "histologies.tsv")
+methyl_file <- file.path(data_dir, "v14", "methyl-beta-values.rds")
+mb_shh_file <- file.path(results_dir, "mb_shh_subtypes_w_molecular_data.tsv")
+```
+
+Wrangle data
+```{r get methyl ids}
+hist <- read_tsv(hist_file)
+
+mb_shh_subtypes <- read_tsv(mb_shh_file)
+```
+
+Filter hist for mb shh methyl samples, and append to subtype df
+```{r}
+hist_mb_methyl <- hist %>%
+  dplyr::filter(pathology_diagnosis == "Medulloblastoma",
+                experimental_strategy == "Methylation") %>%
+  dplyr::rename(Kids_First_Biospecimen_ID_methyl = Kids_First_Biospecimen_ID)
+
+mb_shh_subtypes <- read_tsv(mb_shh_file) %>%
+  left_join(hist_mb_methyl %>%
+              dplyr::select(match_id, Kids_First_Biospecimen_ID_methyl)) %>%
+  dplyr::filter(!is.na(Kids_First_Biospecimen_ID_methyl)) %>%
+  distinct(match_id, Kids_First_Biospecimen_ID_methyl, .keep_all = TRUE) %>%
+  # redefine un-subtyped samples as "unk"
+  dplyr::mutate(molecular_subtype = case_when(
+    molecular_subtype == "MB, SHH" ~ "MB, SHH unk",
+    TRUE ~ molecular_subtype
+  )) %>%
+  dplyr::mutate(molecular_subtype = fct_relevel(molecular_subtype,
+                                               c("MB, SHH alpha", "MB, SHH beta",
+                                                 "MB, SHH gamma", "MB, SHH delta",
+                                                 "MB, SHH unk")))
+```
+
+Get number of samples by MB SHH subtype
+```{r}
+table(mb_shh_subtypes$SHH_subtype)
+```
+
+Load methylation data and filter for ids in `mb_shh_subtypes`
+```{r load methyl}
+methyl <- readRDS(methyl_file)
+
+mb_methyl <- methyl[,colnames(methyl) %in% c("Probe_ID", hist_mb_methyl$Kids_First_Biospecimen_ID_methyl)]
+
+mb_methyl <- mb_methyl %>%
+  distinct(Probe_ID, .keep_all = TRUE) %>%
+  column_to_rownames("Probe_ID")
+```
+
+Identify 10k most variable probes among MB samples
+```{r}
+mb_methyl_var <- apply(mb_methyl, 1, var, na.rm = TRUE)
+
+mb_var_probes <- names(sort(mb_methyl_var, decreasing = TRUE)[1:20000])
+```
+
+```{r}
+set.seed(2024)
+
+mb_umap_results <- umap::umap(t(mb_methyl[mb_var_probes, ]))
+mb_umap_plot_df <- data.frame(mb_umap_results$layout) %>%
+  tibble::rownames_to_column("Kids_First_Biospecimen_ID_methyl") %>%
+  left_join(hist_mb_methyl)
+```
+
+Plot UMAP with molecular subtype and age range 
+```{r}
+mb_umap_plot_df %>%
+  ggplot(aes(x = X1, 
+             y = X2,
+             color = molecular_subtype)) +
+  geom_point(alpha = 0.7) +
+  labs(color = "molecular subtype") +
+  theme_bw() +
+  xlab("UMAP1") +
+  ylab("UMAP2") 
+
+ggsave(file.path(plots_dir, "umap_mb.pdf"),
+       width = 5.5, height = 3.5)
+```
+
+Identify 10k most variable probes among MB Group 3/4 samples
+```{r}
+g34_samples <- hist_mb_methyl %>%
+  dplyr::filter(molecular_subtype %in% c("MB, Group3", "MB, Group4")) %>%
+  pull(Kids_First_Biospecimen_ID_methyl)
+
+mb_g34_methyl_var <- apply(mb_methyl[,colnames(mb_methyl) %in% g34_samples], 1, var, na.rm = TRUE)
+
+mb_g34_var_probes <- names(sort(mb_g34_methyl_var, decreasing = TRUE)[1:20000])
+```
+
+```{r}
+set.seed(2024)
+
+mb_g34_umap_results <- umap::umap(t(mb_methyl[mb_g34_var_probes, colnames(mb_methyl) %in% g34_samples]))
+mb_g34_umap_plot_df <- data.frame(mb_g34_umap_results$layout) %>%
+  tibble::rownames_to_column("Kids_First_Biospecimen_ID_methyl") %>%
+  left_join(hist_mb_methyl)
+```
+
+Plot UMAP with molecular subtype and age range 
+```{r}
+mb_g34_umap_plot_df %>%
+  dplyr::filter(grepl("MB_G34", dkfz_v12_methylation_subclass))  %>%
+  ggplot(aes(x = X1, 
+             y = X2,
+             color = dkfz_v12_methylation_subclass,
+             shape = molecular_subtype)) +
+  geom_point(alpha = 0.7) +
+  labs(color = "methylation subtype",
+       shape = "molecular subtype") +
+  theme_bw() +
+  xlab("UMAP1") +
+  ylab("UMAP2") 
+
+ggsave(file.path(plots_dir, "umap_mb_group34.pdf"),
+       width = 5.5, height = 3.5)
+```
+
+
+Identify 10k most variable probes among MB SHH samples
+```{r}
+mb_shh_methyl_var <- apply(mb_methyl[,colnames(mb_methyl) %in% mb_shh_subtypes$Kids_First_Biospecimen_ID_methyl], 1, var, na.rm = TRUE)
+
+mb_shh_var_probes <- names(sort(mb_shh_methyl_var, decreasing = TRUE)[1:20000])
+```
+
+Generate UMAP df 
+```{r}
+set.seed(2024)
+
+mb_shh_umap_results <- umap::umap(t(mb_methyl[mb_shh_var_probes, colnames(mb_methyl) %in% mb_shh_subtypes$Kids_First_Biospecimen_ID_methyl]))
+mb_shh_umap_plot_df <- data.frame(mb_shh_umap_results$layout) %>%
+  tibble::rownames_to_column("Kids_First_Biospecimen_ID_methyl") %>%
+  inner_join(mb_shh_subtypes)
+
+mb_shh_umap_plot_df <- mb_shh_umap_plot_df %>%
+  dplyr::mutate(age_range = case_when(
+    age_at_diagnosis_years <= 5 ~ "0-5",
+    age_at_diagnosis_years <= 10 ~ "5-10",
+    age_at_diagnosis_years <= 15 ~ "10-15",
+    TRUE ~ ">15"
+  )) %>%
+  dplyr::mutate(age_range = fct_relevel(age_range,
+                                        c("0-5", "5-10",
+                                          "10-15", ">15"))) %>%
+  dplyr::mutate(consensus_CN_MYCN = case_when(
+    is.na(consensus_CN_MYCN) ~ "neutral",
+    TRUE ~ consensus_CN_MYCN
+  )) %>%
+    dplyr::mutate(consensus_CN_GLI2 = case_when(
+    is.na(consensus_CN_GLI2) ~ "neutral",
+    TRUE ~ consensus_CN_GLI2
+  )) %>%
+    dplyr::mutate(consensus_CN_CCND2 = case_when(
+    is.na(consensus_CN_CCND2) ~ "neutral",
+    TRUE ~ consensus_CN_CCND2
+  )) %>%
+    dplyr::mutate(consensus_CN_PTEN = case_when(
+    is.na(consensus_CN_PTEN) ~ "neutral",
+    TRUE ~ consensus_CN_PTEN
+  )) %>%
+  dplyr::mutate(classification_source = case_when(
+    classification_source == "Genomic/Expression" ~ "Molecular",
+    is.na(classification_source) ~ "Unavailable",
+    TRUE ~ classification_source
+  )) %>%
+  write_tsv(file.path(results_dir, "mb_shh_subtypes_w_molecular_umap_data.tsv"))
+```
+
+Plot UMAP with molecular subtype, classification source, and age range 
+```{r}
+mb_shh_umap_plot_df %>%
+  ggplot(aes(x = X1, 
+             y = X2,
+             color = molecular_subtype,
+             size = age_range,
+             shape = classification_source)) +
+  geom_point(alpha = 0.7) +
+  labs(color = "molecular subtype",
+       size = "age range (years)",
+       shape = "classifcation source") +
+  theme_bw() +
+  xlab("UMAP1") +
+  ylab("UMAP2") +
+  # colors to match subtypes in Garcia-Lopez 2020 review
+  scale_color_manual(values = c("aquamarine3", "goldenrod2", 
+                                "royalblue1", "plum4",
+                                "gray")) +
+  guides(color = guide_legend(order = 1),
+         shape  = guide_legend(order = 2),
+         size  = guide_legend(order = 3))
+
+ggsave(file.path(plots_dir, "umap_mb_shh.pdf"),
+       width = 6.5, height = 4.5)
+```
+
+Plot UMAP with methylation subtype and methylation score
+
+```{r}
+mb_shh_umap_plot_df %>%
+  dplyr::filter(grepl("MB", dkfz_v12_methylation_subclass_collapsed)) %>%
+
+  ggplot(aes(x = X1, 
+             y = X2,
+             color = dkfz_v12_methylation_subclass_collapsed,
+             alpha = dkfz_v12_methylation_subclass_score_mean)) +
+  geom_point(size = 3) +
+  labs(color = "methylation subtype",
+       alpha = "methylation subtype score") +
+  theme_bw() +
+  xlab("UMAP1") +
+  ylab("UMAP2") +
+  # colors to match subtypes in Garcia-Lopez 2020 review
+  scale_color_manual(values = c("goldenrod2", "royalblue1",
+                                "aquamarine3", "plum4"))
+
+ggsave(file.path(plots_dir, "umap_mb_shh_methylation_subtype.pdf"),
+       width = 5.5, height = 3.5)
+```
+
+Plot UMAP with CN status for MYCN, GLI2, CCND2, and PTEN
+
+```{r}
+umap_plot_cn_df <- mb_shh_umap_plot_df %>%
+  dplyr::select(molecular_subtype, X1, X2,
+                consensus_CN_MYCN,
+                consensus_CN_GLI2,
+                consensus_CN_CCND2,
+                consensus_CN_PTEN) %>%
+  dplyr::rename(MCYN = consensus_CN_MYCN,
+                GLI2 = consensus_CN_GLI2,
+                CCND2 = consensus_CN_CCND2,
+                PTEN = consensus_CN_PTEN) %>%
+  gather(key = "gene_name", value = "CN_status",
+         -molecular_subtype, -X1, -X2)
+
+umap_plot_cn_df %>%
+  ggplot(aes(x = X1, 
+             y = X2,
+             color = molecular_subtype,
+             shape = CN_status)) +
+  geom_point(alpha = 0.7, size = 3) +
+  labs(color = "methylation subtype",
+       shape = "CN status") +
+  facet_wrap(~gene_name, nrow = 2) +
+  theme_bw() +
+  xlab("UMAP1") +
+  ylab("UMAP2") +
+  # colors to match subtypes in Garcia-Lopez 2020 review
+  scale_color_manual(values = c("aquamarine3", "goldenrod2", 
+                                "royalblue1", "plum4",
+                                "gray"))
+
+ggsave(file.path(plots_dir, "umap_mb_shh_cn_status.pdf"),
+       width = 8, height = 5.5)
+```
+
+Plot UMAP with TP53 alteration status
+
+```{r}
+mb_shh_umap_plot_df %>%
+  ggplot(aes(x = X1, 
+             y = X2,
+             color = molecular_subtype,
+             shape = tp53_status)) +
+  geom_point(alpha = 0.7, size = 3) +
+  labs(color = "methylation subtype",
+       shape = "TP53 status") +
+  theme_bw() +
+  xlab("UMAP1") +
+  ylab("UMAP2") +
+  # colors to match subtypes in Garcia-Lopez 2020 review
+  scale_color_manual(values = c("aquamarine3", "goldenrod2", 
+                                "royalblue1", "plum4",
+                                "gray"))
+
+ggsave(file.path(plots_dir, "umap_mb_shh_tp53_status.pdf"),
+       width = 5.5, height = 3.5)
+```
+
+Identify 10k most variable probes among MB samples
+```{r}
+wnt_samples <- hist_mb_methyl %>%
+  dplyr::filter(molecular_subtype %in% c("MB, WNT")) %>%
+  pull(Kids_First_Biospecimen_ID_methyl)
+
+mb_wnt_methyl_var <- apply(mb_methyl[,colnames(mb_methyl) %in% wnt_samples], 1, var, na.rm = TRUE)
+
+mb_wnt_var_probes <- names(sort(mb_wnt_methyl_var, decreasing = TRUE)[1:20000])
+```
+
+```{r}
+set.seed(2024)
+
+mb_wnt_umap_results <- umap::umap(t(mb_methyl[mb_wnt_var_probes, colnames(mb_methyl) %in% wnt_samples]))
+mb_wnt_umap_plot_df <- data.frame(mb_wnt_umap_results$layout) %>%
+  tibble::rownames_to_column("Kids_First_Biospecimen_ID_methyl") %>%
+  left_join(hist_mb_methyl)
+
+```
+
+Plot UMAP with molecular subtype and age range 
+```{r}
+mb_wnt_umap_plot_df %>%
+  ggplot(aes(x = X1, 
+             y = X2,
+             color = dkfz_v12_methylation_subclass)) +
+  geom_point(alpha = 0.7, size = 3) +
+  labs(color = "methylation subtype") +
+  theme_bw() +
+  xlab("UMAP1") +
+  ylab("UMAP2") 
+
+ggsave(file.path(plots_dir, "umap_mb_wnt.pdf"),
+       width = 5.5, height = 3.5)
+```