510_T_cell_plotting_umap_comp_heatmap.Rmd

---
title: "MSK SPECTRUM freeze major subset deep dives"
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_output_type: console
---


# Figure 5 upper

```{r chunk_510_010}

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

coi <- "T.super"
cell_sort <- "CD45+"
cell_type_major <- "T.cell"
louvain_resolution <- 0.3
louvain_cluster <- "RNA_snn_res.0.3"
pcut <- 0.05

# p1 <- ggplot(data.frame(x = 1:5, y = 1:5), aes(x, y)) + geom_point()
# vp_outer <- viewport(x = 0.5, y = 0.5, width = 0.5, height = 0.5, angle = 90)
# vp_inner <- viewport(x = 0.5, y = 0.5, width = 1, height = 1, angle = 0)
# grid.newpage()
# pushViewport(vp_outer)
# grid.rect()
# print(p1, vp = vp_inner)
# popViewport(1)
# g1 <- grid.grab()
# ggdraw() + draw_grob(g1)


```

```{r chunk_510_020}

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

# scrna_meta_tbl
# clrs
# markers_v7
# markers_v7_super
# cell_type_super_lookup

```

```{r chunk_510_030}

myfeatures <- c("UMAP_1", "UMAP_2", "umapharmony_1", "umapharmony_2", "sample", "doublet", "nCount_RNA", "nFeature_RNA", "percent.mt", "doublet_score", "cell_type", "cluster_label_sub")

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

=======
gois <- c("CD4", "CD8A", "CD8B", "FOXP3", "GNLY", "MKI67")

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

>>>>>>> 3b877a7c59812e1ad259a7239a3b4b8e355e7ff4
marker_tbl <- read_tsv("/work/shah/users/uhlitzf/data/SPECTRUM/freeze/v7/supplementary_tables/T.super_marker_table_annotated_full.tsv")

marker_tbl_top <- marker_tbl %>% 
  filter(avg_logFC > 0.5,
         p_val_adj < 0.01,
         pct.1 > 0.2,
         pct.2 < 0.8,
         !is.na(cluster_label_sub),
         !str_detect(gene, "^RPS|^RPL")) %>% 
  group_by(cluster_label_sub) %>% 
  slice(1:50)

marker_sheet <- read_tsv("/work/shah/users/uhlitzf/data/SPECTRUM/freeze/v7/supplementary_tables/T.super_marker_sheet_full.tsv")

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

plot_data_sub <- as_tibble(FetchData(seu_obj_sub, c(myfeatures))) %>% 
  left_join(scrna_meta_tbl, by = "sample") %>% 
  mutate(cell_type_super = cell_type_super_lookup[cell_type]) %>% 
  mutate(sort_short = str_remove_all(sort_parameters, "singlet, live, ")) %>% 
  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_sub = ordered(cluster_label_sub, levels = names(clrs$cluster_label_sub[[coi]]))) %>% 
  mutate(umapscaled_1 = scales::rescale(umapharmony_1),
         umapscaled_2 = scales::rescale(umapharmony_2))
  
plot_data_sub <- filter(plot_data_sub, sort_short_x == cell_sort, !is.na(tumor_supersite))

plot_data_genes <- as_tibble(FetchData(seu_obj_sub, c("cell_id", "umapharmony_1", "umapharmony_2", "sample", "cell_type", gois))) %>% 
  gather(gene, value, -umapharmony_1, -umapharmony_2, -sample, -cell_type, -cell_id) %>% 
  left_join(scrna_meta_tbl, by = "sample") %>% 
  mutate(cell_type_super = cell_type_super_lookup[cell_type]) %>% 
  mutate(sort_short = str_remove_all(sort_parameters, "singlet, live, ")) %>% 
  mutate(sort_short_x = ifelse(sort_short == "U" & cell_type_super == "Immune", 
                               "CD45+", ifelse(sort_short == "U" & cell_type_super == "Stromal", "CD45-", sort_short)))
  
plot_data_genes <- filter(plot_data_genes, sort_short_x == cell_sort, !is.na(tumor_supersite))

```


## umaps 

```{r chunk_510_50, fig.width=7, fig.height=4}

alpha_lvl <- ifelse(nrow(plot_data_sub) < 20000, 0.2, 0.1)
pt_size <- ifelse(nrow(plot_data_sub) < 20000, 0.2, 0.05)

common_layers_disc <- list(  
  ggrastr::geom_point_rast(size = pt_size, alpha = alpha_lvl, 
                           raster.dpi = 50, raster.width = 10, raster.height = 10),
  NoAxes(),
  guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))),
  labs(color = ""),
  theme(aspect.ratio = 1)
)

common_layers_cont <- list(  
  ggrastr::geom_point_rast(size = pt_size, alpha = alpha_lvl, 
                           raster.dpi = 50, raster.width = 10, raster.height = 10),
  NoAxes(),
  scale_color_gradientn(colors = viridis(9)),
  guides(color = guide_colorbar()),
  theme(aspect.ratio = 1)
)

umap_coord_anno <- function(size) {
  ggplot(tibble(group = c("UMAP1", "UMAP2"),
                x = c(0, 0), xend = c(1, 0),
                y = c(0, 0), yend = c(0, 1),
                lx = c(0.5, -0.15), ly = c(-0.15, 0.5),
                angle = c(0, 90))) +
    geom_segment(aes(x, y, xend = xend, yend = yend, group = group),
                 arrow = arrow(angle = 20, type = "closed", length = unit(0.1, "npc")),
                 size = size/4, lineend = "round") +
    geom_text(aes(lx, ly, label = group, angle = angle), size = size) +
    theme_void() +
    coord_fixed(xlim = c(-0.3, 1), ylim = c(-0.3, 1))
}

add_umap_coord <- function(gg_obj, x = 0, y = 0, width = 0.3, height = 0.3, size = 4.5) {
  p <- ggdraw() + 
    draw_plot(gg_obj, x = 0, y = 0, width = 1, height = 1) +
    draw_plot(umap_coord_anno(size = size), x = x, y = y, width = width, height = height)
  return(p)
}

```

### cluster umap

```{r chunk_510_60, fig.width=6, fig.height=5}

plot_data_label <- plot_data_sub %>% 
  group_by(cluster_label_sub) %>% 
  summarise(umapharmony_1 = median(umapharmony_1),
            umapharmony_2 = median(umapharmony_2),
            umapscaled_1 = median(umapscaled_1),
            umapscaled_2 = median(umapscaled_2)) %>% 
  mutate(cluster_number = as.numeric(cluster_label_sub))

umap_cell_type <- ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = cluster_label_sub)) + 
  common_layers_disc +
  #facet_wrap(~cluster_label_sub) + 
  geom_point(aes(umapharmony_1, umapharmony_2), color = "white",
             data = plot_data_label, alpha = 0.5, size = 6) +
  geom_text(aes(umapharmony_1, umapharmony_2, label = cluster_number), color = "black",
             data = plot_data_label) +
  scale_color_manual(values = clrs$cluster_label_sub[[coi]]) +
  guides(color = guide_legend(override.aes = list(size = 2, alpha = 1), nrow = 12)) +
  annotate(geom = "text", 
           x = (max(plot_data_genes$umapharmony_1) + min(plot_data_genes$umapharmony_1))/2,
           y = min(plot_data_genes$umapharmony_2), color = "black",
           label = paste0("n = ", nrow(distinct(plot_data_genes, cell_id))))

umap_site <- ggplot(plot_data_sub, aes(umapharmony_1, umapharmony_2, color = tumor_supersite)) + 
  common_layers_disc +
  scale_color_manual(values = clrs$tumor_supersite)

umap_cell_type_void <- umap_cell_type + guides(color = F)
umap_site_void <- umap_site + guides(color = F)

umap_cell_type_legend <- cowplot::get_legend(umap_cell_type)
umap_site_legend <- cowplot::get_legend(umap_site)

umap_site <- ggplot(filter(plot_data_sub, tumor_supersite %in% c("Ascites", "Adnexa")),
                    aes(umapharmony_1, umapharmony_2, color = tumor_megasite)) + 
  ggrastr::geom_point_rast(aes(umapharmony_1, umapharmony_2), color = "grey80",
                           size = 0.01, alpha = 0.1,
                           raster.dpi = 50, raster.width = 10, raster.height = 10,
                           data = select(plot_data_sub, -tumor_supersite)) + 
  common_layers_disc +
  # facet_wrap(~tumor_supersite) + 
  scale_color_manual(values = clrs$tumor_megasite)

add_umap_coord(umap_cell_type_void)
ggsave_pdf("figures/510_T_cell_plotting_umap_comp_heatmap/004_umap_cluster_label_sub.pdf", width = 6, height = 5)
ggsave_png("figures/510_T_cell_plotting_umap_comp_heatmap/004_umap_cluster_label_sub.png", width = 6, height = 5)
```


### gene umaps

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

umap_gois <- ggplot(plot_data_genes) +
  ggrastr::geom_point_rast(aes(umapharmony_1, umapharmony_2),
                           size = pt_size, alpha = alpha_lvl, color = "grey80",
                           raster.dpi = 50, raster.width = 10, raster.height = 10,
                           data = filter(plot_data_genes, value == 0))  + 
  ggrastr::geom_point_rast(aes(umapharmony_1, umapharmony_2, color = value),
                           size = pt_size, alpha = alpha_lvl, 
                           raster.dpi = 50, raster.width = 10, raster.height = 10,
                           data = filter(plot_data_genes, value > 0)) + 
  NoAxes() + 
  scale_color_gradientn(colors = viridis(9)) +
  guides(color = guide_colorbar()) +
  theme(aspect.ratio = 1) +
  facet_wrap(~gene) +
  labs(color = "Expression")

umap_gois
ggsave_pdf("figures/510_T_cell_plotting_umap_comp_heatmap/004_umap_gois.pdf", width = 8, height = 4)
ggsave_png("figures/510_T_cell_plotting_umap_comp_heatmap/004_umap_gois.png", width = 8, height = 4)


```

### kernel umap

#### site

```{r chunk_510_62, fig.width=5, fig.height=5}

source("src/umap_kernel.R")

kernel_tbl <- kde2d_contrast(
  data_tbl = plot_data_sub,
  x = umapharmony_1, y = umapharmony_2, kernel_group = tumor_megasite,
  kernel_subsets_a = c("Adnexa"), kernel_subsets_b = c("Other", "Ascites")
) %>%
  filter(value.x > 0.001 & value.y > 0.001)

rdbu <- rev(RColorBrewer::brewer.pal(7, "RdBu"))

kernel_umap <- ggplot() +
  ggrastr::geom_point_rast(aes(umapscaled_1, umapscaled_2), color = "grey80", size = 0.01, alpha = 0.02,
                           data = plot_data_sub,
                           raster.dpi = 50, raster.width = 10, raster.height = 10) +
  ggrastr::geom_point_rast(aes(x, y, color = value_quenched), data = kernel_tbl, size = 0.01, 
                           raster.dpi = 50, raster.width = 10, raster.height = 10) +
  geom_point(aes(umapscaled_1, umapscaled_2), color = "white",
             data = plot_data_label, alpha = 0.5, size = 6) +
  geom_text(aes(umapscaled_1, umapscaled_2, label = cluster_number), color = "black",
            data = plot_data_label) +
  scale_color_gradientn(colours = rdbu,
                        values = scales::rescale(c(min(kernel_tbl$value_quenched), 0,
                                                   max(kernel_tbl$value_quenched)))) +
  guides(color = guide_colorbar(label.position = "right",
                                title.position = "top",
                                title.hjust = 0,
                                title.vjust = 1,
                                direction = "vertical")) +
  theme(aspect.ratio = 1,
        legend.key.height = unit(0.05, "npc"),
        legend.key.width = unit(0.03, "npc"),
        legend.position = c(-0.1, 0.95),
        legend.justification = c("left", "top"),
        plot.title = element_text(hjust = 0.5, vjust = 0.5, face = "plain", size = 18)) + 
  labs(color = "Enrichment\nin adnexa") +
  remove_xaxis +
  remove_yaxis

add_umap_coord(kernel_umap)
ggsave_pdf("figures/510_T_cell_plotting_umap_comp_heatmap/004_umap_kernel.pdf", width = 6, height = 5)
ggsave_png("figures/510_T_cell_plotting_umap_comp_heatmap/004_umap_kernel.png", width = 6, height = 5)

```

#### mutsig

```{r chunk_510_63, fig.width=5, fig.height=5}

kernel_umap_wrapper <- function(kernel_subsets_a = c("FBI"),
                                kernel_subsets_b = c("HRD-Dup", "HRD-Del", "TD"),
                                quench = 0.015,
                                color_title = paste0("Enrichm.\nin ", kernel_subsets_a)) {
  
  kernel_tbl <- kde2d_contrast(
    data_tbl = plot_data_sub %>% filter(tumor_supersite != "Ascites"),
    x = umapharmony_1, y = umapharmony_2, kernel_group = consensus_signature,
    kernel_subsets_a = kernel_subsets_a, kernel_subsets_b = kernel_subsets_b,
    quench = quench
  ) %>%
    filter(value.x > 0.001 & value.y > 0.001)
  
  kernel_umap <- ggplot() +
    ggrastr::geom_point_rast(aes(umapscaled_1, umapscaled_2), 
                             color = "grey80", size = 0.01, alpha = 0.02, data = plot_data_sub,
                             raster.dpi = 50, raster.width = 10, raster.height = 10) +
    ggrastr::geom_point_rast(aes(x, y, color = value_scaled), data = kernel_tbl, size = 0.01, 
                             raster.dpi = 50, raster.width = 10, raster.height = 10) +
    geom_point(aes(umapscaled_1, umapscaled_2), color = "white",
               data = plot_data_label, alpha = 0.1, size = 6) +
    geom_text(aes(umapscaled_1, umapscaled_2, label = cluster_number), 
              color = "black",
              data = plot_data_label) +
    scale_color_gradientn(colours = rdbu,
                          values = scales::rescale(c(min(kernel_tbl$value_scaled), 0,
                                                     max(kernel_tbl$value_scaled)))) +
    guides(color = guide_colorbar(label.position = "right",
                                  title.position = "top",
                                  title.hjust = 0,
                                  title.vjust = 1,
                                  direction = "vertical")) +
    theme(aspect.ratio = 1,
          legend.key.height = unit(0.05, "npc"),
          legend.key.width = unit(0.03, "npc"),
          legend.position = c(-0.1, 0.95),
          legend.justification = c("left", "top"),
          plot.title = element_text(hjust = 0.5, vjust = 0.5, face = "plain", size = 18)) + 
    labs(color = color_title,
         title = paste0(kernel_subsets_a, collapse = ", ")) +
    remove_xaxis +
    remove_yaxis
  
  return(kernel_umap)
}

k1 <- kernel_umap_wrapper("HRD-Dup", c("FBI", "HRD-Del", "TD"), 0.015, "Enrichm.")
k2 <- kernel_umap_wrapper("HRD-Del", c("FBI", "HRD-Dup", "TD"), 0.015) + remove_guides
k3 <- kernel_umap_wrapper("FBI", c("HRD-Dup", "HRD-Del", "TD"), 0.015) + remove_guides
k4 <- kernel_umap_wrapper("TD", c("FBI", "HRD-Dup", "HRD-Del"), 0.015) + remove_guides

kernel_grid <- ggdraw() +
  draw_plot(plot_grid(k1 + remove_guides, k2, k3, k4, ncol = 2, align = "hv"),
            width = 0.8, height = 1) +
  draw_grob(get_legend(k1), x = 0.88, y = -0.25)

kernel_umap_wrapper(c("HRD-Dup", "HRD-Del"), c("FBI", "TD"), 0.01, "Enrichm. in\nHRD vs HRP") + labs(title = "")

ggsave_pdf("figures/510_T_cell_plotting_umap_comp_heatmap/004_umap_kernel_hrd_hrp.pdf", width = 6, height = 5)
ggsave_png("figures/510_T_cell_plotting_umap_comp_heatmap/004_umap_kernel_hrd_hrp.png", width = 6, height = 5)

```


## module heatmap

```{r chunk_510_64, fig.width=15, fig.height=5}

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

plot_data_modules <- as_tibble(FetchData(seu_obj_sub, c(myfeatures, module_names))) %>% 
  gather(pathway, 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 = str_remove_all(sort_parameters, "singlet, live, ")) %>% 
  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_sub = ordered(cluster_label_sub, levels = names(clrs$cluster_label_sub[[coi]]))) %>% 
  mutate(key_group = ifelse(str_detect(pathway, "pathway"), "pathway", "module")) %>% 
  mutate(pathway = str_remove_all(pathway, ".pathway|.module")) %>% 
  mutate(cluster_label_super = case_when(
    str_detect(cluster_label_sub, "^CD4") ~ "CD4",
    str_detect(cluster_label_sub, "^CD8") ~ "CD8",
    str_detect(cluster_label_sub, "^NK") ~ "NK",
    str_detect(cluster_label_sub, "^Cycling") ~ "Cycling",
    T ~ "Other"
    )) %>% 
  mutate(cluster_label_super = ordered(cluster_label_super, levels = c("CD4", "CD8", "Other", "NK", "Cycling")))
  
plot_data_modules <- filter(plot_data_modules, sort_short_x == cell_sort, !is.na(tumor_supersite))

plot_data_modules_summary_sample <- plot_data_modules %>% 
  group_by(pathway, key_group, sample, patient_id_short, cluster_label_sub, cluster_label_super) %>% 
  summarise(mean_score = mean(value, na.rm = T)) %>% 
  ungroup

plot_data_modules_summary_cluster <- plot_data_modules %>% 
  group_by(pathway, key_group, cluster_label_sub, cluster_label_super) %>% 
  summarise(mean_score = mean(value, na.rm = T)) %>% 
  ungroup %>% 
  mutate(cluster_number = as.numeric(cluster_label_sub)) %>% 
  mutate(cluster_label_sub = fct_rev(cluster_label_sub))

```

```{r chunk_510_65, fig.width=12, fig.height=4.5}

rdbu <- rev(RColorBrewer::brewer.pal(7, "RdBu"))

common_heat_layers <- list(
  # scale_fill_gradient2(low = scales::muted("blue"), high = scales::muted("red"), 
  #                      na.value = "grey10", 
  #                      # breaks = c(-cut_value, 0, cut_value), 
  #                      # labels = c(paste0("≤-", cut_value), "0", paste0("≥", cut_value)),
  #                      # limits = c(-cut_value, cut_value)
  #                      ),
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        axis.ticks = element_blank(),
        axis.line = element_blank(),
        axis.title = element_blank(),
        plot.margin = margin(0, 0, 0, 0),
        strip.text = element_blank())
)

pw_plot1 <- plot_data_modules_summary_cluster %>% 
  filter(key_group == "pathway", pathway %in% c("JAK.STAT", "TNFa", "TGFb")) %>% 
  ggplot() +
  geom_tile(aes(pathway, cluster_label_sub, fill = mean_score)) +
  common_heat_layers +
  facet_grid(cluster_label_super~., scales = "free", space = "free") + 
  labs(x = "Cluster", y = "Pathway", fill = "Mean\nPROGENy\nscore") +
  scale_fill_gradientn(colours = rdbu,
                        values = scales::rescale(c(min(filter(plot_data_modules_summary_cluster, pathway %in% c("JAK.STAT", "TNFa", "TGFb"))$mean_score), 0, max(filter(plot_data_modules_summary_cluster, pathway %in% c("JAK.STAT", "TNFa", "TGFb"))$mean_score)))) + 
  scale_x_discrete(position = "top") +
  theme(axis.text.y = element_blank(),
        axis.text.x.top = element_text(angle = 90, hjust = 0, vjust = 0.5),
        legend.position = "bottom") +
  guides(fill = guide_colorbar(title.position = "top"))

pw_plot2 <- plot_data_modules_summary_cluster %>% 
  filter(key_group == "module") %>% 
  mutate(pathway = str_remove_all(pathway, "CD8.")) %>% 
  mutate(pathway = ordered(pathway, levels = c("Naive", "Predysfunctional", "Cytotoxic", "Dysfunctional"))) %>% 
  ggplot() +
  geom_tile(aes(pathway, cluster_label_sub, fill = mean_score)) +
  common_heat_layers +
  facet_grid(cluster_label_super~., scales = "free", space = "free") + 
  labs(x = "Cluster", y = "Module", fill = "Mean\nmodule\nscore") +
  scale_fill_gradientn(colours = rdbu,
                        values = scales::rescale(c(min(filter(plot_data_modules_summary_cluster, key_group == "module")$mean_score), 0, max(filter(plot_data_modules_summary_cluster, key_group == "module")$mean_score)))) + 
  scale_x_discrete(position = "top") +
  theme(axis.text.y = element_blank(),
        axis.text.x.top = element_text(angle = 90, hjust = 0, vjust = 0.5),
        legend.position = "bottom") +
  guides(fill = guide_colorbar(title.position = "top"))

pw_plot1_anno <- plot_data_modules_summary_cluster %>% 
  distinct(cluster_label_sub, cluster_number, cluster_label_super) %>% 
  mutate(facet_helper = "") %>% 
  ggplot() +
  geom_tile(aes(facet_helper, cluster_label_sub, fill = cluster_label_sub)) +
  geom_point(aes(facet_helper, cluster_label_sub), color = "white", size = 6, alpha = 0.5) +
  geom_text(aes(facet_helper, cluster_label_sub, label = cluster_number)) +
  common_heat_layers + 
  scale_fill_manual(values = clrs$cluster_label_sub[[coi]]) +
  facet_grid(cluster_label_super~., scales = "free", space = "free") +
  theme(axis.text.x = element_blank(),
        strip.text = element_blank()) +
  guides(fill = F)

# heat_module_grid <- plot_grid(pw_plot1, pw_plot2, pw_plot1_anno, rel_heights = c(0.45, 0.15, 0.4), ncol = 1, align = "v", axis = "lrtb")
# heat_module_grid <- plot_grid(pw_plot1, pw_plot2, pw_plot1_anno, rel_heights = c(0.2, 0.24, 0.56), ncol = 1, align = "v", axis = "lrtb")
heat_module_grid <- plot_grid(pw_plot1_anno, 
                              ggdraw(),
                              pw_plot2 + remove_guides, 
                              ggdraw(),
                              pw_plot1 + remove_guides, 
                              rel_widths = c(0.57, 0.01, 0.22, 0.03, 0.17), 
                              nrow = 1, align = "h", axis = "tb")

heat_module_grid

ggsave_pdf("figures/510_T_cell_plotting_umap_comp_heatmap/004_module_heatmap.pdf", heat_module_grid, width = 4.5, height = 12)
ggsave_png("figures/510_T_cell_plotting_umap_comp_heatmap/004_module_heatmap.png", heat_module_grid, width = 4.5, height = 12)

```

[figures/510_T_cell_plotting_umap_comp_heatmap/004_module_heatmap.pdf](figures/510_T_cell_plotting_umap_comp_heatmap/004_module_heatmap.pdf)

[figures/510_T_cell_plotting_umap_comp_heatmap/004_module_heatmap.png](figures/510_T_cell_plotting_umap_comp_heatmap/004_module_heatmap.png)

## marker heatmap 

```{r chunk_510_70, fig.width=15, fig.height=5}

plot_data_markers <- as_tibble(FetchData(seu_obj_sub, c(myfeatures, unique(marker_tbl_top$gene)))) %>% 
  gather(gene, value, -c(1:(length(myfeatures)+1))) %>% 
  left_join(scrna_meta_tbl, by = "sample") %>% 
  mutate(cluster_label_sub = ordered(cluster_label_sub, levels = my_subtypes)) %>% 
  group_by(cluster_label_sub, gene) %>% 
  summarise(value = mean(value, na.rm = T)) %>% 
  group_by(gene) %>% 
  mutate(value = scales::rescale(value)) %>% 
  left_join(select(marker_tbl_top, cluster_label_sub_x = cluster_label_sub, gene), by = "gene") %>% 
  mutate(cluster_label_sub_x = ordered(cluster_label_sub_x, levels = rev(names(clrs$cluster_label_sub[[coi]])))) %>% 
  na.omit()

```

```{r chunk_510_75, fig.width=15, fig.height=5}

highlight_genes <- marker_tbl_top %>% 
  group_by(cluster_label_sub) %>% 
  slice(1:2) %>% 
  mutate(cluster_label_sub_x = ordered(cluster_label_sub, levels = rev(names(clrs$cluster_label_sub[[coi]])))) %>% 
  ungroup() %>% 
  select(cluster_label_sub_x, gene) %>% 
  na.omit %>% 
  mutate(highlight = T)

plot_data_markers_mat <- plot_data_markers %>% 
  spread(cluster_label_sub, value) %>% 
  left_join(highlight_genes, by = c("gene", "cluster_label_sub_x")) %>% 
  arrange(desc(cluster_label_sub_x), gene) %>% 
  select(cluster_label_sub_x, gene, highlight, everything())

ha_row <- rowAnnotation(
  `Cell type` = plot_data_markers_mat$cluster_label_sub_x,
  col = list(`Cell type` = clrs$cluster_label_sub[[coi]]),
  show_legend = F,
  annotation_name_side = "bottom"
)

ha_col <- columnAnnotation(
  `Cell type` = anno_block(gp = gpar(fill = clrs$cluster_label_sub[[coi]], col = NA),
                           labels = as.character(1:(ncol(plot_data_markers_mat)-3)),
                           labels_rot = 90),
  show_legend = F,
  annotation_name_side = "left",
  annotation_name_rot = 0
)

gene_idx <- plot_data_markers_mat$highlight == TRUE

ha_genes <- rowAnnotation(
  link = anno_mark(
    at = which(gene_idx),
    labels = plot_data_markers_mat$gene[which(gene_idx)],
    labels_gp = gpar(fontsize = 10), padding = unit(1, "mm"),
    side = "left",
    labels_rot = 180
  )
)

marker_heatmap <- Heatmap(
  as.matrix(plot_data_markers_mat[,-c(1:3)]), 
  heatmap_legend_param = list(
    title = "Scaled expression", 
    title_position = "leftcenter-rot"
  ),
  row_order = 1:length(plot_data_markers_mat$cluster_label_sub_x),
  row_split = plot_data_markers_mat$cluster_label_sub_x, 
  column_split = 1:length(colnames(plot_data_markers_mat)[-c(1:3)]),
  column_order = colnames(plot_data_markers_mat)[-c(1:3)], 
  column_names_side = "bottom",
  right_annotation = ha_row,
  left_annotation = ha_genes,
  bottom_annotation = ha_col,
  cluster_rows = F, 
  row_title = NULL,
  column_title = NULL,
  col = viridis(9)
)

# marker_heatmap

```

## composition

### per site

```{r chunk_510_080, fig.width=3, fig.height=9}

source("src/comp_plot.R")

comp_tbl_sample_sort <- plot_data_sub %>%
  filter(therapy == "pre-Rx") %>%
  group_by(sample, tumor_subsite, tumor_supersite, tumor_megasite, patient_id_short,
           therapy, sort_short_x, consensus_signature, cluster_label_sub) %>%
  tally %>%
  group_by(sample, tumor_subsite, tumor_supersite, tumor_megasite, patient_id_short,
           therapy, sort_short_x, consensus_signature) %>%
  mutate(nrel = n/sum(n)*100,
         log10n = log10(n)) %>%
  mutate(sample_id = sample) %>%
  mutate(tumor_supersite = ordered(tumor_supersite, levels = rev(names(clrs$tumor_supersite)))) %>% 
  ungroup %>% 
  mutate(cell_type_naive = ifelse(str_detect(cluster_label_sub, "mem"), "T naive/mem", "other")) %>% 
  mutate(cell_type_dysfunc = ifelse(str_detect(cluster_label_sub, "dysfunc"), "T dysfunctional", "other"))# %>%
  # # Pad cell type combinations with zero counts in the data
  # complete(cell_type_naive,
  #          nesting(sample, tumor_subsite, tumor_supersite, tumor_megasite,
  #                  patient_id_short, therapy, sort_short_x, consensus_signature),
  #          fill = list(n = 0)) %>%
  # complete(cell_type_dysfunc,
  #          nesting(sample, tumor_subsite, tumor_supersite, tumor_megasite,
  #                  patient_id_short, therapy, sort_short_x, consensus_signature),
  #          fill = list(n = 0))

plist1 <- default_comp_grid_list(
  comp_tbl_sample_sort, cell_type_naive, "T naive/mem",
  cluster_label_sub, super_type_sub = "T.super", vec_plot = T, nmax = 5000, highlight = T)

plist2 <- default_comp_grid_list(
  comp_tbl_sample_sort, cell_type_dysfunc, "T dysfunctional",
  cluster_label_sub, super_type_sub = "T.super", vec_plot = T, nmax = 5000, highlight = T)

comp_grid1 <- plot_grid(plotlist = plist1[-3], ncol = 1, align = "v", rel_heights = c(0.15, 0.15, 0.2, 0.5))

comp_grid2 <- plot_grid(plotlist = plist2[-3], ncol = 1, align = "v", rel_heights = c(0.15, 0.15, 0.2, 0.5))

comp_grid1
comp_grid2
ggsave_pdf("figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_memory.pdf", comp_grid1, width = 3, height = 9)
ggsave_pdf("figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_dysfunc.pdf", comp_grid2, width = 3, height = 9)
ggsave_png("figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_memory.png", comp_grid1, width = 3, height = 9)
ggsave_png("figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_dysfunc.png", comp_grid2, width = 3, height = 9)

```

[figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_memory.pdf](figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_memory.pdf)

[figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_memory.png](figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_memory.png)

[figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_dysfunc.pdf](figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_dysfunc.pdf)

[figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_dysfunc.png](figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_comp_dysfunc.png)

```{r chunk_510_110, fig.width=20, fig.height=8}

heatmap_grob <- grid.grabExpr(draw(marker_heatmap), width = 8, height = 16)
vp_outer <- viewport(x = 0.5, y = 0.5, width = 0.5, height = 2, angle = 270)
vp_inner <- viewport(x = 0.5, y = 0.5, width = 1, height = 1, angle = 0)
grid.newpage()
pushViewport(vp_outer)
grid.draw(heatmap_grob)
popViewport(1)
heatmap_grob_rot <- grid.grab()

umap_heat_grid <- ggdraw() +
  draw_plot(add_umap_coord(umap_cell_type_void, size = 3.5, x = 0),
            x = 0, y = 0.5, width = 0.2, height = 0.5) +
  # draw_plot(add_umap_coord(umap_site_void, size = 3.5, x = 0),
  #           x = 0, y = 0, width = 0.2, height = 0.5) +
  draw_grob(heatmap_grob_rot,
            x = 0.2, y = 0, width = 0.8, height = 1)
  # draw_plot_label(c("A", "B", "C", "D"), x = c(0, 0.63, 0, 0), y = c(0.995, 0.995, 0.72, 0.45))

umap_heat_grid

ggsave_png("figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_umap_heat.png", umap_heat_grid, width = 20, height = 8)
ggsave_pdf("figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_umap_heat.pdf", umap_heat_grid, width = 20, height = 8)

```

## mini umap inserts

```{r chunk_510_120, fig.width=6, fig.height=1.5}

plot_data_mini_umap <- plot_data_modules %>% 
  distinct(umapharmony_1, umapharmony_2, cluster_label_sub,
           cluster_label_super) %>% 
  mutate(dysfunc_naivemem = case_when(
    str_detect(cluster_label_sub, "dysfunc") ~ "Dysfunctional",
    str_detect(cluster_label_sub, "naive|mem") ~ "Naive/memory",
    T ~ "Other",
  )) %>% 
  mutate(dysfunc_naivemem = ordered(dysfunc_naivemem, 
                                    levels = c("Naive/memory", "Dysfunctional", "Other")))

umap_mini_cluster_label_sub <- ggplot() + 
  ggrastr::geom_point_rast(
    aes(umapharmony_1, umapharmony_2), 
    data = select(plot_data_mini_umap, -cluster_label_super),
    alpha = 0.02, size = 0.01, color = "grey80",
    raster.dpi = 15, raster.width = 40, raster.height = 40
  ) +
  ggrastr::geom_point_rast(
    aes(umapharmony_1, umapharmony_2, color = cluster_label_sub), 
    data = select(plot_data_mini_umap, -cluster_label_super),
    alpha = 0.02, size = 0.01, 
    raster.dpi = 15, raster.width = 40, raster.height = 40
  ) +
  ggrastr::geom_point_rast(
    aes(umapharmony_1, umapharmony_2, color = cluster_label_sub), 
    data = filter(plot_data_mini_umap, cluster_label_super != "Other"),
    alpha = 0.5, size = 0.01,
    raster.dpi = 15, raster.width = 40, raster.height = 40
  ) +
  facet_wrap(~cluster_label_super, ncol = 1) + 
  theme_void() +
  scale_color_manual(values = clrs$cluster_label_sub$T.super) +
  remove_guides + 
  theme(aspect.ratio = 1, 
        panel.spacing = unit(0.1, "npc"), 
        strip.text = element_text(size = 14))

umap_mini_cluster_label_sub
ggsave_pdf("figures/510_T_cell_plotting_umap_comp_heatmap/005_umap_mini_cluster_label_sub.pdf", umap_mini_cluster_label_sub, width = 1.5, height = 6)
ggsave_png("figures/510_T_cell_plotting_umap_comp_heatmap/005_umap_mini_cluster_label_sub.png", umap_mini_cluster_label_sub, width = 1.5, height = 6)

```

```{r chunk_510_130, fig.width=3, fig.height=1}

umap_mini_dysfunc_naivemem <- ggplot() + 
  ggrastr::geom_point_rast(
    aes(umapharmony_1, umapharmony_2), 
    data = select(plot_data_mini_umap, -dysfunc_naivemem),
    alpha = 0.02, size = 0.01, color = "grey80",
    raster.dpi = 15, raster.width = 40, raster.height = 40
  ) +
  ggrastr::geom_point_rast(
    aes(umapharmony_1, umapharmony_2, color = cluster_label_sub), 
    data = select(plot_data_mini_umap, -dysfunc_naivemem),
    alpha = 0.02, size = 0.01, 
    raster.dpi = 15, raster.width = 40, raster.height = 40
  ) +
  ggrastr::geom_point_rast(
    aes(umapharmony_1, umapharmony_2, color = cluster_label_sub), 
    data = filter(plot_data_mini_umap, dysfunc_naivemem != "Other"),
    size = 0.01,
    raster.dpi = 15, raster.width = 40, raster.height = 40
  ) +
  facet_grid(~dysfunc_naivemem) + 
  theme_void() +
  scale_color_manual(values = clrs$cluster_label_sub$T.super) +
  remove_guides + 
  theme(aspect.ratio = 1, 
        panel.spacing.x = unit(0.2, "npc"),
        strip.text = element_text(size = 12))

umap_mini_dysfunc_naivemem

ggsave_pdf("figures/510_T_cell_plotting_umap_comp_heatmap/005_umap_mini_dysfunc_naivemem.pdf", umap_mini_dysfunc_naivemem, width = 3, height = 1)
ggsave_png("figures/510_T_cell_plotting_umap_comp_heatmap/005_umap_mini_dysfunc_naivemem.png", umap_mini_dysfunc_naivemem, width = 3, height = 1)


```

## full heatmap grid with mini umaps

```{r chunk_510_140, fig.width=6, fig.height=1.5}

ggdraw() +
  draw_plot(umap_mini_cluster_label_sub, x = -0.05, y = 0.15, width = 3/7.5, height = 0.68) +
  draw_plot(heat_module_grid, x = 0.3, y = 2/14, width = 4.5/7.5, height = 12/14) + 
  draw_grob(get_legend(pw_plot2), x = 0.62, y = -0.42) +
  draw_grob(get_legend(pw_plot1), x = 0.82, y = -0.42)

ggsave_pdf("figures/510_T_cell_plotting_umap_comp_heatmap/004_module_heatmap_mini_umap.pdf", width = 7.5, height = 14)
ggsave_png("figures/510_T_cell_plotting_umap_comp_heatmap/004_module_heatmap_mini_umap.png", width = 7.5, height = 14)

```


[figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_umap_heat.png](figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_umap_heat.png)

[figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_umap_heat.pdf](figures/510_T_cell_plotting_umap_comp_heatmap/004_T.cell_umap_heat.pdf)

## session info

```{r chunk_510_999}

devtools::session_info()

```