5f_bulk_hnsc_single_cell_validation.Rmd

---
title: "Figure 6F: Enrichment of TCGA-HNSC deconvolution results in single cell data"
output:
  html_document:
    df_print: paged
    self_contained: yes
---


```{r}
source("../R/figure_utils.R")
source("../R/setup.R")
library(Seurat)
```

```{r}
save_dir <- "../out/dualsimplex_save_tcga_hnscc_4kg_9ct_no_hk_neighbours_alexey/"
custom_cols <- c(
  "#1f78b4", "#33a02c", "#e31a1c", "#ff7f00", "#6a3d9a",
  "#b15928", "#a6cee3", "#b2df8a", "#fb9a99", "#fdbf6f"
)
```

```{r}
lo2 <- DualSimplexSolver$from_state(save_dir, save_there = F)
```

Markers for enrichment in single cell (100 per cell type)
```{r fig.width = 14, fig.height = 6}
lo2$plot_projected("markers", use_dims = NULL, with_legend = T, with_history = F)
```
```{r}
plt <- lo2$plot_projected("zero_distance", use_dims = NULL, with_legend = T, with_history = F)
ggsave("../out/5c_zero__disance_with_legend.svg", plt, width = 6.5, height = 5, device = svglite)

```

```{r fig.width=6.5, fig.height=5}
seurat_obj <- readRDS("../data/large/GSE181919_integrated.rds")
so_malignant <- readRDS("../data/large/GSE181919_malignant_integrated.rds")
so_epithelial <- readRDS("../data/large/GSE181919_epithelial_integrated.rds")

plt <- Seurat::DimPlot(
  seurat_obj,
  group.by = "cell.type",
  reduction = "tsne",
  cols = custom_cols,
  raster = TRUE
) + Seurat::NoAxes() + ggtitle("GSE181919 dataset")
show(plt)
Seurat::DimPlot(so_malignant, group.by = "cell.type", reduction = "umap")
Seurat::DimPlot(so_epithelial, group.by = "cell.type", reduction = "umap")

ggsave("../out/6f_gse181919_dataset.svg", plt, width = 6.5, height = 5, device = svg)
```


```{r fig.width=20, fig.height=15}
seurat_obj <- add_list_markers(seurat_obj, lo2$get_marker_genes())
so_malignant <- add_list_markers(so_malignant, lo2$get_marker_genes())
so_epithelial <- add_list_markers(so_epithelial, lo2$get_marker_genes())

plot_marker_enrichment(seurat_obj, lo2$get_ct_names(), limits = NULL)
plot_marker_enrichment(so_malignant, lo2$get_ct_names(), limits = NULL)
plot_marker_enrichment(so_epithelial, lo2$get_ct_names(), limits = NULL)
```

```{r fig.width=10, fig.height=15}
cancer_cts <- c("cell_type_4", "cell_type_5", "cell_type_6", "cell_type_8", "cell_type_9")

plot_list <- plot_marker_enrichment(
  so_malignant,
  cancer_cts,
  limits = NULL,
  ncol = 2,
  ggadd = function(x, y) { x + Seurat::NoAxes() + Seurat::NoLegend() },
  raster = T,
  pt.size = 3, wrap = F
)

res_plot <- ggarrange(plotlist = plot_list, ncol=2, nrow=3, common.legend = T, legend = "none")

ggsave("../out/5f_malignant_subset.svg", res_plot, width = 10, height = 15, device = svg)
res_plot
```
```{r fig.width=10, fig.height=15}
plot_list <- plot_marker_enrichment(
  so_malignant,
  cancer_cts,
  limits = NULL,
  ncol = 2,
  ggadd = function(x, y) { x + Seurat::NoAxes() },
  raster = T,
  pt.size = 3, wrap = F
)

res_plot <- ggarrange(plotlist = plot_list, ncol=2, nrow=3, common.legend = T, legend = "right")

ggsave("../out/5f_malignant_subset_with_legend.svg", res_plot, width = 10, height = 15, device = svglite)
res_plot

```

```{r fig.width=4.7, fig.height=4}
plt <- Seurat::DimPlot(so_malignant, group.by = "hpv", raster = T,   pt.size = 3) + Seurat::NoAxes() + ggtitle("")
ggsave("../out/6f_malignant_subset_hpv.svg", plt, width = 4.7, height = 4, device = svg)
plt
```


```{r}
table(seurat_obj$cell.type)
```

```{r, fig.width=7, fig.height=4}
enrichment_map <- as.data.frame(seurat_obj@meta.data %>% group_by(cell.type) %>% summarize(
  across(paste0(lo2$get_ct_names(), 1), list(function(x) {
      tc <- 1000
    mean(tail(sort(x), tc))
  })), .groups = "drop"
))

rownames(enrichment_map) <- enrichment_map$cell.type
enrichment_map$cell.type <- NULL
colnames(enrichment_map) <- substr(colnames(enrichment_map), 1, nchar(colnames(enrichment_map)) - 3)

col_order <- c(
  "cell_type_1", "cell_type_2", "cell_type_3", "cell_type_7", cancer_cts
)

row_order <- c(
  "Macrophages", "Myocytes", "Epithelial.cells", "Fibroblasts", "Malignant.cells"
)
enrichment_map <- enrichment_map[row_order, col_order]

#enrichment_map_norm <- as.data.frame(t(apply(enrichment_map, 1, function(x) (x - min(x)) / (max(x) - min(x)))))
enrichment_map_norm <- as.data.frame(apply(enrichment_map, 2, function(x) (x - min(x)) / (max(x) - min(x))))

plt <- pheatmap::pheatmap(
  enrichment_map_norm,
  display_numbers = round(enrichment_map, 2),
  cluster_rows = F,
  cluster_cols = F,
  treeheight_row = 0,
  treeheight_col = 0,
  legend = F
)
ggsave("../out/6f_heatmap.svg", plt, width = 7, height = 4, device = svglite::svglite)
plt
```