Add analysis for medulloblastoma subtypes, immune cell distributions.

Dockerfile

@@ -217,6 +217,7 @@ RUN R -e 'BiocManager::install(c( \
   "ggplot2", \
   "ggpubr", \
   "ggsignif", \
+  "ggstatsplot", \
   "glmnet", \
   "glmnetUtils", \
   "gplots", \

analyses/immune-deconv/immune-deconv-mb-subtypes.R

@@ -0,0 +1,27 @@
+# Load Libraries
+suppressPackageStartupMessages({
+  library(optparse)
+  library(tidyverse)
+})
+source("util/heatmap_by_subtype.R")
+
+# parse parameters
+option_list <- list(
+  make_option(c("--xcell_output_path"), type = "character",
+              help = "path to xCell output"),
+  make_option(c("--output_heatmap_path"), type = "character",)
+  )
+
+opt <- parse_args(OptionParser(option_list = option_list))
+xcell_output_path <- opt$xcell_output_path
+output_heatmap_path <- opt$output_heatmap_path
+
+# Read in the output of 01-immune-deconv.R 
+xcell_output <- readRDS(xcell_output_path)
+
+# Select Medulloblastoma samples 
+mb <- xcell_output %>% filter(cancer_group == 'Medulloblastoma')
+
+# Make heatmap
+heatmap_by_subtype(mb, annot_colors = NA, output_file = output_heatmap_path)
+

analyses/immune-deconv/mb-subtype-exploration.Rmd

@@ -0,0 +1,100 @@
+---
+title: "mb-subtype-exploration"
+author: "Nicholas Rodriguez"
+date: "2024-12-15"
+output:
+  html_document: default
+  pdf_document: default
+editor_options:
+  chunk_output_type: inline
+---
+
+## Set up
+
+### Libraries and functions
+
+```{r}
+library(tidyverse)
+library(readr)
+library(ggstatsplot) 
+source("util/heatmap_by_subtype.R")
+```
+
+### Read in immune deconv output
+```{r}
+resultsDir <- "results"
+# Read in quantiseq output
+qsData <- read_rds(file.path(resultsDir, "quantiseq_output.rds"))
+xcellData <- read_rds(file.path(resultsDir, "xcell_output.rds"))
+```
+
+### Dataset
+
+#### Sample size
+```{r}
+# Include only medulloblastoma cancer group
+xcellData <- xcellData %>% filter(cancer_group=="Medulloblastoma")
+qsData <- qsData %>% filter(cancer_group=="Medulloblastoma")
+qsData %>% select(Kids_First_Biospecimen_ID) %>% unique() %>% count() 
+```
+
+#### Subtype counts of medulloblastoma (MB)
+
+```{r}
+qsData %>% count(molecular_subtype)
+```
+
+
+#### Cell types identified by quantiseq
+
+```{r}
+qsData %>% select(cell_type) %>% unique()
+```
+#### Cell types identified by xCell 
+```{r}
+xcellData %>% select(cell_type) %>% unique()
+```
+
+### Boxplot of quantiseq fractions
+```{r, fig.height=8}
+qsDataFilteredNoUCCell <- qsData %>% filter(cell_type != "uncharacterized cell")
+qsDataFilteredNoUCCell %>% ggplot(aes(x = cell_type, y = fraction)) + 
+  geom_boxplot() +
+  theme(axis.text.x = element_text(angle = 45, hjust = 1), ) + 
+  scale_y_continuous(limits = c(0, 1))
+```
+
+### Quantiseq Heatmap 
+```{r}
+heatmap_by_subtype(qsDataFilteredNoUCCell, annot_colors = NA, output_file = NA)
+```
+### xCell Heatmap
+
+```{r, fig.width=10}
+xcellFilteredNoUCCell <- xcellData %>% filter(cell_type != "uncharacterized cell")
+heatmap_by_subtype(xcellFilteredNoUCCell, annot_colors = NA, output_file = NA)
+```
+
+### Anova 
+```{r}
+summary(aov(fraction ~ molecular_subtype,
+  data = qsDataFilteredNoUCCell,
+))
+```
+
+
+```{r, fig.height=6}
+
+ggbetweenstats(
+  data = xcellFilteredNoUCCell,
+  x = molecular_subtype,
+  y = fraction,
+  type = "parametric", # for ANOVA
+  plot.type = "violin",
+  pairwise.comparisons = TRUE,
+  pairwise.display = "significant",
+  centrality.plotting = TRUE,
+  bf.message = TRUE
+)
+```
+

analyses/immune-deconv/run-immune-deconv-mb-subtypes.sh

@@ -0,0 +1,22 @@
+#!/bin/bash
+# Module author: Komal S. Rathi, updated Kelsey Keith
+# 2022-07
+
+# This script generates a heatmap for medulloblastoma sutypes. 
+
+set -e
+set -o pipefail
+
+# This script should always run as if it were being called from
+# the directory it lives in.
+script_directory="$(perl -e 'use File::Basename;
+  use Cwd "abs_path";
+  print dirname(abs_path(@ARGV[0]));' -- "$0")"
+cd "$script_directory" || exit
+
+### 
+# generate deconvolution output
+echo "Heatmap MB Subtypes"
+Rscript --vanilla immune-deconv-mb-subtypes.R \
+--xcell_output_path 'results/xcell_output.rds' \
+--output_heatmap_path 'results/heatmap_mb_subtypes.png'

analyses/immune-deconv/util/heatmap_by_subtype.R

@@ -0,0 +1,44 @@
+# Create a heatmap from `molecular_subtype` of immunedeconv output.
+# Modified from https://github.com/rokitalab/OpenPedCan-Project/blob/dev/analyses/immune-deconv/util/heatmap_by_group.R
+# Assumes deconv data is filtered for a single `cancer_group`. 
+
+# load libraries
+suppressPackageStartupMessages({
+  library(tidyverse)
+  library(pheatmap)
+})
+
+# function to create heatmap of average immune scores per cell type per cancer and gtex group 
+heatmap_by_subtype <- function(deconv_output, annot_colors, output_file) {
+
+  # create a generalized group column
+  deconv_output <- deconv_output %>%
+    mutate(group = ifelse(cohort == "GTEx", gtex_group, cancer_group))
+
+  # create labels: count of samples per group
+  deconv_output <- deconv_output %>%
+    group_by(group, cell_type) %>%
+    unique() %>%
+    mutate(label = n()) %>%
+    mutate(label = molecular_subtype)
+
+  # calculate mean scores per cell type per histology
+  deconv_output <- deconv_output %>%
+    filter(!cell_type %in% c("microenvironment score", "stroma score", "immune score")) %>%
+    group_by(cell_type, label) %>%
+    dplyr::summarise(mean = mean(fraction))
+
+  # convert into matrix of cell type vs histology
+  deconv_output <- deconv_output %>%
+    spread(key = label, value = mean) %>%
+    column_to_rownames('cell_type')
+
+  # heatmap of average immune scores per cell type per cancer and gtex group
+  deconv_output %>%
+    t() %>%
+    pheatmap(scale = "column", angle_col = 45,
+             main = "Average Immune scores",
+             annotation_row = NA, annotation_colors = annot_colors,
+             annotation_legend = T, cellwidth = 14, cellheight = 10,
+             filename = output_file)
+}