620_Macrophage_umaps_and_boxplots.Rmd

---
title: "MSK SPECTRUM freeze Macrophages"
author: "Florian Uhlitz"
date: "`r Sys.Date()`"
output:
  bookdown::html_document2:
    includes:
      in_header: header.html
  bookdown::gitbook:
    includes:
      in_header: header.html
editor_options: 
  chunk_620_output_type: console
---


# Figure 6 lower

```{r chunk_620_010}

library(magrittr)
library(tidyverse)
library(Seurat)
library(readxl)
library(cowplot)
library(colorblindr)
library(viridis)
library(ComplexHeatmap)
library(grid)
library(ggpubr)

coi <- "Myeloid.super"
cell_sort <- "CD45+"
cell_type_major <- c("Myeloid.cell", "Dendritic.cell", "Mast.cell")
louvain_resolution <- 0.3
louvain_cluster <- "RNA_snn_res.0.3"
pcut <- 0.05

```

```{r chunk_620_020}

## load global vars: 
source("src/global_vars.R")

# scrna_meta_tbl
# clrs
# markers_v7
# markers_v7_super
# cell_type_super_lookup

```

```{r chunk_620_030}

seu_obj_ml <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/Myeloid.super_processed_filtered_annotated.rds")

seu_obj_mp <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/Macrophages_processed.rds")

seu_obj_dc <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/DCs_processed.rds")

seu_obj_cc <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/Ovarian.cancer.super_processed_filtered_annotated.rds")

seu_obj_tc <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/T.super_processed_filtered_annotated.rds")


gois_mp <- c("S100A8", "S100A9", "CXCL10", "CD274", "C1QC", "FN1", "MARCO", "SELENOP")
gois_dc <- grep("^IRF", rownames(seu_obj_dc), value = T)
gois <- c(gois_mp, gois_dc)
sois <- c("HRD-Dup", "HRD-Del", "FBI", "TD")

seu_obj_dc <- AddModuleScore(seu_obj_dc, list(IRF.genes.module = c("IRF1", "IRF7", "IRF8", "IRF9")), name = "IRF.genes.module")
seu_obj_dc$IRF.genes.module <- seu_obj_dc$IRF.genes.module1
seu_obj_dc$IRF.genes.module1 <- NULL

module_names <- grep("pathway|CD8|module", colnames(seu_obj_dc@meta.data), value = T)

module_names_sub <- c("MARCO", "CXCL10", "CD274", "JAK.STAT.pathway")

module_names_dc <- c("Interferon.regulators.module", "IRF.genes.module", "IRF1", "IRF7", "IRF8", "IRF9")

myfeatures <- c("cell_id", "umapharmony_1", "umapharmony_2", "sample", "cell_type", "cluster_label", "Phase")

my_subtypes <- names(clrs$cluster_label[[coi]])

plot_data <- as_tibble(FetchData(seu_obj_mp, c(myfeatures, module_names, gois))) %>%
  gather(key, value, -c(1:length(myfeatures))) %>% 
  left_join(scrna_meta_tbl, by = "sample") %>%
  mutate(cell_type_super = cell_type_super_lookup[cell_type]) %>%
  mutate(sort_short_x = ifelse(sort_short == "U" & cell_type_super == "Immune",
                               "CD45+", ifelse(sort_short == "U" & cell_type_super == "Stromal", "CD45-", sort_short))) %>% 
  mutate(cluster_label = ordered(cluster_label, levels = my_subtypes),
         goi_expressed = value > 0) %>% 
  group_by(sample, key) %>% 
  mutate(n = length(sample)) %>% 
  mutate(pct.expr = sum(goi_expressed)/n*100) %>% 
  ungroup

plot_data_uniq <- distinct(plot_data, cell_id, .keep_all = T)

plot_data_dc <- as_tibble(FetchData(seu_obj_dc, c(myfeatures, module_names, gois))) %>%
  gather(key, value, -c(1:length(myfeatures))) %>% 
  left_join(scrna_meta_tbl, by = "sample") %>%
  mutate(cell_type_super = cell_type_super_lookup[cell_type]) %>%
  mutate(sort_short_x = ifelse(sort_short == "U" & cell_type_super == "Immune",
                               "CD45+", ifelse(sort_short == "U" & cell_type_super == "Stromal", "CD45-", sort_short))) %>% 
  mutate(cluster_label = ordered(cluster_label, levels = my_subtypes),
         goi_expressed = value > 0) %>% 
  group_by(sample, key) %>% 
  mutate(n = length(sample)) %>% 
  mutate(pct.expr = sum(goi_expressed)/n*100) %>% 
  ungroup

plot_data_uniq_dc <- distinct(plot_data_dc, cell_id, .keep_all = T)


```


## Macrophage umap 

```{r chunk_620_120, fig.width=8, fig.height=4}

dc_plot_cluster <- ggplot(plot_data_uniq) +
  ggrastr::geom_point_rast(aes(umapharmony_1, umapharmony_2, color = cluster_label), 
                           size = 0.01, alpha = 1, raster.dpi = 150) +
  # theme_void() +
  scale_color_manual(values = clrs$cluster_label[[coi]]) +
  guides(color = guide_legend(override.aes = list(size = 3, alpha = 1))) +
  theme(aspect.ratio = 1,
        axis.line = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank()) + 
  labs(color = "Cluster")

plot_data_gois <- plot_data %>% 
  filter(key %in% gois_mp) %>% 
  mutate(key = ordered(key, levels = gois_mp)) %>% 
  mutate(value = ifelse(value > 5, 5, value))

dc_plot_genes <- ggplot() +
  ggrastr::geom_point_rast(aes(umapharmony_1, umapharmony_2), 
                           size = 0.1, alpha = 1, color = "grey80", raster.dpi = 75,
                           data = filter(plot_data_gois, value == 0)) +
  ggrastr::geom_point_rast(aes(umapharmony_1, umapharmony_2, color = value), 
                           size = 0.1, alpha = 1, raster.dpi = 75,
                           data = filter(plot_data_gois, value > 0)) +
  facet_wrap(~key, nrow = 2) + 
  theme(axis.text = element_blank(),
        axis.title = element_blank(),
        axis.line = element_blank(),
        axis.ticks = element_blank()) +  
  scale_color_gradientn(colors = viridis(9), breaks = c(0, 5), 
                        labels = c("low", "high")) +
  theme(aspect.ratio = 1) +
  labs(color = "Expression")

dc_plot_genes
ggsave_pdf("figures/620_Macrophage_umaps_and_boxplots/005_macrophage_gois_umap.pdf", dc_plot_genes, width = 8, height = 4)
ggsave_png("figures/620_Macrophage_umaps_and_boxplots/005_macrophage_gois_umap.png", dc_plot_genes, width = 8, height = 4)

```

## violins

### macrophages

```{r chunk_620_140, fig.width=6.5, fig.height=3.5}

plot_data_vio_mp <- plot_data %>%
  filter(consensus_signature %in% sois) %>%
  filter(tumor_supersite != "Ascites") %>%
  mutate(key = str_replace_all(key, "JAK.STAT.pathway", "JAK/STAT\npathway"))

plot_data_vio_mp_summary <- plot_data_vio_mp %>%
  group_by(sample, key) %>%
  mutate(value = mean(value),
         pct.expr = mean(pct.expr)) %>%
  distinct(sample, key, .keep_all = T)

my_comparisons <- list(
  c("HRD-Dup", "HRD-Del"), 
  c("HRD-Dup", "FBI"), 
  c("HRD-Del", "FBI")
)

my_comparisons2 <- list(
  c("HRD", "HRP")
)


boxplot_wrapper <- function(plot_data_vio, plot_data_vio_summary, group_column, y, 
                            keys, comps, group_var_label, y_label, ylims, ncol) {
  
  group_column <- enquo(group_column)
  y <- enquo(y)
  yrange <- ylims[2]-ylims[1]
  ylabels <- yrange*c(0.95,0.85,0.75)+ylims[1]
  clrs$cluster_label <- clrs$cluster_label$Ovarian.cancer.super
  plot_data_vio_summary <- filter(plot_data_vio_summary, key %in% keys) %>% 
    mutate(key = ordered(key, levels = keys))
  p <- plot_data_vio %>% 
    filter(key %in% keys) %>% 
    mutate(key = ordered(key, levels = keys)) %>% 
    ggplot(aes(!!group_column, !!y)) +
    geom_violin(aes(!!group_column, !!y, fill = !!group_column), color = "white", adjust = 2, alpha = 0.5, width = 1.5) +
    geom_boxplot(aes(!!group_column, !!y, color = !!group_column),
                 width = 0.5, size = 0.75, outlier.shape = NA) +
    geom_boxplot(aes(!!group_column, !!y, fill = !!group_column),
                 color = "white", width = 0.5, outlier.shape = NA, size = 0.5) +
    # geom_jitter(aes(!!group_column, value, group = !!group_column),
    #             color = "black", size = 0.01, width = 0.25, alpha = 0.5,
    #             data = filter(plot_data_vio_summary, key %in% keys)) +
    facet_wrap(~key, ncol = ncol) +
    scale_color_manual(values = clrs[[as_label(group_column)]]) +
    scale_fill_manual(values = clrs[[as_label(group_column)]]) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
          axis.title.x = element_blank()) + 
    coord_cartesian(clip = "off", ylim = ylims) +
    labs(y = y_label, color = group_var_label, 
         fill = group_var_label) +
    remove_guides

  if (any(unlist(comps) != ".all.")) {
    p <- p + 
      stat_compare_means(comparisons = comps, 
                         data = plot_data_vio_summary,
                         label = "p.signif", label.y = ylabels,
                         hide.ns = F)
  }
  
  if (all(unlist(comps) == ".all.")) {
    p <- p + 
      stat_compare_means(ref.group = ".all.", 
                         method.args = list(alternative = "greater"),
                         data = plot_data_vio_summary, 
                         label.y = ylabels[1], hide.ns = T,
                         label = "p.signif")
  }
  
  return(p)
  
}

```

```{r chunk_620_141, fig.width=5, fig.height=4}

violins_genes_mutsig_mp_pdl1 <- boxplot_wrapper(
  plot_data_vio_mp, plot_data_vio_mp_summary, 
  consensus_signature, pct.expr, 
  c("CD274"), 
  my_comparisons, "Mutational\nsignature", "% Pos. Macroph.", c(0,25), 5
)

violins_genes_mutsig_mp_cxcl10 <- boxplot_wrapper(
  plot_data_vio_mp, plot_data_vio_mp_summary, 
  consensus_signature, pct.expr, 
  c("CXCL10"), 
  my_comparisons, "Mutational\nsignature", "% Pos. Macroph.", c(0,100), 5
)

violins_pw_mutsig_mp <- boxplot_wrapper(
  plot_data_vio_mp, plot_data_vio_mp_summary, 
  consensus_signature, value, 
  "JAK/STAT\npathway", 
  my_comparisons, "Mutational\nsignature", "PROGENy score", c(-1, 5), 5
)

violins_grid_mutsig_mp <- plot_grid(
  violins_genes_mutsig_mp_pdl1, violins_genes_mutsig_mp_cxcl10, violins_pw_mutsig_mp,
  nrow = 1, align = "h")

violins_grid_mutsig_mp

ggsave_pdf("figures/620_Macrophage_umaps_and_boxplots/005_macrophage_boxplots.pdf", violins_grid_mutsig_mp, width = 5, height = 3)
ggsave_png("figures/620_Macrophage_umaps_and_boxplots/005_macrophage_boxplots.png", violins_grid_mutsig_mp, width = 5, height = 3)

```

```{r chunk_620_142, fig.width=5, fig.height=4}

violins_genes_mutsig_short_mp_pdl1 <- boxplot_wrapper(
  plot_data_vio_mp, plot_data_vio_mp_summary, 
  consensus_signature_short, pct.expr, 
  c("CD274"), 
  my_comparisons2, "Mutational\nsignature", "% Pos. Macroph.", c(0,25), 5
)

violins_genes_mutsig_short_mp_cxcl10 <- boxplot_wrapper(
  plot_data_vio_mp, plot_data_vio_mp_summary, 
  consensus_signature_short, pct.expr, 
  c("CXCL10"), 
  my_comparisons2, "Mutational\nsignature", "% Pos. Macroph.", c(0,100), 5
)

violins_pw_mutsig_short_mp <- boxplot_wrapper(
  plot_data_vio_mp, plot_data_vio_mp_summary, 
  consensus_signature_short, value, 
  "JAK/STAT\npathway", 
  my_comparisons2, "Mutational\nsignature", "PROGENy score", c(-1, 5), 5
)

violins_grid_mutsig_short_mp <- plot_grid(
  violins_genes_mutsig_short_mp_pdl1, violins_genes_mutsig_short_mp_cxcl10, violins_pw_mutsig_short_mp,
  nrow = 1, align = "h")

violins_grid_mutsig_short_mp

ggsave_pdf("figures/620_Macrophage_umaps_and_boxplots/005_macrophage_boxplots_mutsig_short.pdf", violins_grid_mutsig_short_mp, width = 5, height = 3)
ggsave_png("figures/620_Macrophage_umaps_and_boxplots/005_macrophage_boxplots_mutsig_short.png", violins_grid_mutsig_short_mp, width = 5, height = 3)

```


### DCs

```{r chunk_620_145, fig.width=5, fig.height=3.5}

plot_data_vio_dc <- plot_data_dc %>%
  filter(consensus_signature %in% sois) %>%
  filter(tumor_supersite != "Ascites") %>%
  mutate(key = str_replace_all(key, "Interferon.regulators.module", "IFN\nregulators"))

plot_data_vio_dc_summary <- plot_data_vio_dc %>%
  group_by(sample, key) %>%
  mutate(value = mean(value),
         pct.expr = mean(pct.expr)) %>%
  ungroup %>% 
  distinct(sample, key, .keep_all = T)

violins_genes_mutsig_dc <- boxplot_wrapper(
  plot_data_vio_dc, plot_data_vio_dc_summary, 
  consensus_signature, pct.expr, 
  c("IRF1", "IRF7", "IRF8", "IRF9"), 
  my_comparisons, "Mutational\nsignature", "% Pos. DCs", c(0,100), 5
)

violins_pw_mutsig_dc <- boxplot_wrapper(
  plot_data_vio_dc, plot_data_vio_dc_summary, 
  consensus_signature, value, 
  "IFN\nregulators", 
  my_comparisons, "Mutational\nsignature", "Module score", c(-1, 2), 5
)

violins_genes_mutsig_dc_short <- boxplot_wrapper(
  plot_data_vio_dc, plot_data_vio_dc_summary, 
  consensus_signature_short, pct.expr, 
  c("IRF1", "IRF7", "IRF8", "IRF9"), 
  my_comparisons2, "Mutational\nsignature", "% Pos. DCs", c(0,100), 5
)

violins_pw_mutsig_dc_short <- boxplot_wrapper(
  plot_data_vio_dc, plot_data_vio_dc_summary, 
  consensus_signature_short, value, 
  "IFN\nregulators", 
  my_comparisons2, "Mutational\nsignature", "Module score", c(-1, 2), 5
)



```

```{r chunk_620_146, fig.width=5, fig.height=3.5}

plot_grid(violins_genes_mutsig_dc, violins_pw_mutsig_dc, 
          nrow = 1, rel_widths = c(0.73, 0.27), align = "h")

ggsave_pdf("figures/620_Macrophage_umaps_and_boxplots/005_DC_boxplots.pdf", width = 6, height = 3)
ggsave_png("figures/620_Macrophage_umaps_and_boxplots/005_DC_boxplots.png", width = 6, height = 3)

```

```{r chunk_620_147, fig.width=5, fig.height=3.5}

plot_grid(violins_genes_mutsig_dc_short, violins_pw_mutsig_dc_short, 
          nrow = 1, rel_widths = c(0.73, 0.27))

ggsave_pdf("figures/620_Macrophage_umaps_and_boxplots/005_DC_boxplots_mutsig_short.pdf", width = 6, height = 3)
ggsave_png("figures/620_Macrophage_umaps_and_boxplots/005_DC_boxplots_mutsig_short.png", width = 6, height = 3)

```


```{r chunk_620_150, fig.width=6, fig.height=4}

patient_boxplot_wrapper <- function(xkey) {

  xkey <- enquo(xkey)

  patient_lvls <- plot_data %>%
    mutate(key = !!xkey) %>%
    group_by(patient_id_short) %>%
    summarise(median = median(value)) %>%
    ungroup() %>%
    arrange(median) %>%
    pull(patient_id_short)

  plot_data %>%
    mutate(patient_id_short = ordered(patient_id_short, levels = patient_lvls)) %>%
    ggplot() +
    geom_violin(aes(patient_id_short, !!xkey, fill = consensus_signature),
                color = "white") +
    # geom_boxplot(aes(consensus_signature, value, color = consensus_signature),
    #              width = 0.2, size = 2) +
    geom_boxplot(aes(patient_id_short, !!xkey, fill = consensus_signature),
                 color = "white", width = 0.2, outlier.shape = NA) +
    scale_color_manual(values = clrs$consensus_signature) +
    scale_fill_manual(values = clrs$consensus_signature) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
          axis.title.x = element_blank()) +
    labs(y = "Module score", color = "Mutational\nsignature", title = xkey)
}

# patient_boxplot_list <- lapply(module_names_sub, patient_boxplot_wrapper) %>%
#   setNames(module_names_sub)

```

## correlation

```{r chunk_620_170, fig.width=8, fig.height=3.25}

plot_data_y <- FetchData(seu_obj_cc, c("sample", "cell_type", "JAK.STAT.pathway")) %>%
  bind_rows(FetchData(seu_obj_tc, c("sample", "cell_type", "JAK.STAT.pathway"))) %>% 
  bind_rows(FetchData(seu_obj_mp, c("sample", "cell_type", "JAK.STAT.pathway"))) %>% 
  as_tibble() %>% 
  gather(key_y, value, -sample, -cell_type) %>%
  left_join(scrna_meta_tbl, by = "sample") %>%
  filter((sort_short == "CD45+" & cell_type %in% c("T.cell", "Monocyte")) | (sort_short == "CD45-" & cell_type == "Ovarian.cancer.cell")) %>% 
  group_by(sample, cell_type, key_y, patient_id_short, tumor_supersite, consensus_signature) %>%
  summarise(mean_value_y = mean(value)) %>%
  ungroup %>%
  mutate(sample = str_remove_all(sample, "_CD45N|_CD45P"))

plot_data_cor <- plot_data_vio_dc_summary %>% 
  select(sample, key_x = key, value) %>% 
  left_join(scrna_meta_tbl, by = "sample") %>%
  filter(sort_short == "CD45+") %>%
  group_by(sample, key_x) %>%
  summarise(mean_value_x = mean(value)) %>%
  ungroup %>%
  mutate(sample = str_remove_all(sample, "_CD45N|_CD45P")) %>%
  inner_join(plot_data_y, by = "sample")
 
cor_plot1 <- plot_data_cor %>%
  distinct(sample, cell_type, key_x, key_y, .keep_all = T) %>% 
  filter(key_x == "IFN\nregulators", key_y == "JAK.STAT.pathway") %>%
  mutate(cell_type = case_when(cell_type == "Ovarian.cancer.cell" ~ "Ov cancer cells",
                               cell_type == "T.cell" ~ "T cells",
                               cell_type == "Monocyte" ~ "Macrophages",
                               T ~ "NA")) %>% 
  mutate(cell_type = ordered(cell_type, levels = c("Ov cancer cells", "T cells", "Macrophages"))) %>% 
  ggplot(aes(mean_value_x, mean_value_y, color = consensus_signature)) +
  geom_point(aes(mean_value_x, mean_value_y, color = consensus_signature)) +
  geom_smooth(aes(mean_value_x, mean_value_y), method = "lm", se = T,
              linetype = 2, size = 0.5, color = "black") +
  stat_cor(aes(mean_value_x, mean_value_y),
           method = "spearman", color = "black") +
  scale_color_manual(values = clrs$consensus_signature) +
  facet_grid(~cell_type) + 
  labs(y = "JAK/STAT pathway",
       x = "IFN regulators module score\n(Dendritic cells)",
       color = "Mutational\nsignature")
 
# cor_plot2 <- {cor_plot1 + 
#     labs(y = "JAK-STAT score\n(CD8+ T cells)", 
#          x = "JAK-STAT score\n(cancer cells)",
#          color = "Mutational\nsignature")
# } %+%
#   filter(plot_data_cor, key_tc == "JAK.STAT.pathway", 
#          key_cc == "JAK.STAT.pathway")
# 
# cor_plot_grid <- plot_grid(cor_plot1 + guides(color = F), cor_plot2,
#                            align = "hv", axis = "x", 
#                            rel_widths = c(0.4, 0.6))
# 
# cor_plot_grid
# 

cor_plot1
ggsave_pdf("figures/620_Macrophage_umaps_and_boxplots/005_DC_cancer_cell_correlation.pdf", cor_plot1, width = 8, height = 3.25)
ggsave_png("figures/620_Macrophage_umaps_and_boxplots/005_DC_cancer_cell_correlation.png", cor_plot1, width = 8, height = 3.25)


```


## session info

```{r chunk_620_999}

devtools::session_info()

```