niche composition in space

explore_niche_composition.Rmd

@@ -0,0 +1,190 @@
+---
+title: "What is inside the niches?"
+author: "Ricardo Ramirez"
+date: "5/4/2021"
+output:
+  html_document: default
+  pdf_document: default
+editor_options:
+  chunk_output_type: console
+---
+
+## Introduction
+
+By performing an integration of spots of all slides, I aimed to define classes per spot that represent a niche (or cell community). My intuition is that spots that are similar share compositions or functions throughout samples. This comes from my assumption that the heart tissue is a mosaic and the idea is to identify the color and size of the tiles.
+
+
+```{r setup, include=FALSE}
+library(tidyverse)
+library(Seurat)
+library(ComplexHeatmap)
+library(circlize)
+col_fun = colorRamp2(c(0, 1), c("black", "yellow"))
+```
+
+## Analyzing compositions from a single cell atlas
+
+Load annotations of visium samples
+
+```{r}
+sample_dict <- read.table("./markers/NEW_PatIDs_visium_overview_allsamples.tsv",
+                                   sep = "\t", header = T) %>%
+  mutate_all(as.character) %>%
+  dplyr::select(Visium, New.Ids) %>%
+  dplyr::rename(orig.ident = Visium,
+                patient_sample = New.Ids) %>%
+  dplyr::mutate(patient = map_chr(strsplit(patient_sample, "_"), 
+                                  ~.x[[1]]),
+                area = map_chr(strsplit(patient_sample, "_"), 
+                               ~.x[[2]]))
+```
+
+Load all of the data and visualize the proportions
+
+```{r}
+# Get atlas cell-info ------------------------------------------------------------------------------------
+atlas_meta <- readRDS("./processed_visium/integration/ps_integrated_slides.rds")[[1]][["annotations"]]
+
+# Cell-type composition -----------------------------------
+cell_info <- atlas_meta %>%
+  group_by(orig.ident, opt_clust_integrated) %>%
+  summarize(ncells = length(opt_clust_integrated)) %>%
+  dplyr::mutate(all_sample_cells = sum(ncells)) %>%
+  dplyr::ungroup() %>%
+  dplyr::mutate(cell_prop = ncells/all_sample_cells) %>%
+  left_join(sample_dict) %>%
+  rename(sample = orig.ident,
+         niche_ct_meanprops_id = opt_clust_integrated) %>%
+  mutate(niche_ct_meanprops_id = paste0("niche_", niche_ct_meanprops_id))
+
+head(cell_info)
+```
+
+```{r}
+cell_info %>%
+  ggplot(aes(x = ncells,
+           y = patient_sample,
+           fill = niche_ct_meanprops_id)) + 
+  geom_bar(position="fill", 
+           stat="identity") +
+  theme(legend.position = "bottom") +
+  ylab("") + xlab("n_spots")  
+```
+
+
+```{r, fig.height=6, fig.width=8}
+cell_info_mat <- cell_info %>%
+  dplyr::select(niche_ct_meanprops_id, 
+                patient_sample,
+                cell_prop) %>%
+  pivot_wider(names_from = niche_ct_meanprops_id, values_from  = cell_prop) %>%
+  column_to_rownames("patient_sample") %>%
+  as.matrix() 
+
+cell_info_mat %>%
+  Heatmap(col = col_fun)
+
+cell_info_mat[is.na(cell_info_mat)] <- 0
+
+cell_info_mat %>%
+  scale() %>%
+  Heatmap()
+```
+
+
+## Compositions of cell types within each niche
+
+```{r}
+visium_folder = "./processed_visium/objects/"
+
+visium_files <- list.files(visium_folder, full.names = F)
+visium_samples <- gsub("[.]rds", "", visium_files)
+
+visium_df <- tibble(visium_file = paste0(visium_folder, 
+                                         visium_files),
+                    sample = visium_samples) %>%
+  mutate(niche_ct_meanprops = map(visium_file, function(f_path) {
+    
+    visium_slide <- readRDS(f_path)
+    Idents(visium_slide) <- "opt_clust_integrated"
+    
+    niche_cors <- map(set_names(levels(Idents(visium_slide))),
+                      function(niche) {
+                        
+                        niche_cells <- WhichCells(visium_slide, idents = c(niche))
+                        
+                        prop_mat <- GetAssayData(visium_slide, 
+                                                 assay = "c2l_major_props") %>%
+                          as.matrix()
+                        
+                        prop_mat <- prop_mat[, niche_cells, drop = F] %>%
+                          rowMeans()
+                      })
+    
+  })) %>%
+  unnest_longer("niche_ct_meanprops")
+
+visium_df <- visium_df %>%
+  dplyr::select(sample, niche_ct_meanprops, niche_ct_meanprops_id) %>%
+  mutate(niche_ct_meanprops = map(niche_ct_meanprops, function(x) {
+    
+    df_props <- x %>%
+      as.data.frame()
+    
+    colnames(df_props) = "mean_prop"
+    
+    df_props %>% rownames_to_column("cell_type")
+    
+  })) %>% unnest() %>%
+  arrange(niche_ct_meanprops_id, cell_type, sample) %>%
+  left_join(cell_info, by = c("sample", "niche_ct_meanprops_id")) %>%
+  mutate(niche_id = paste0(patient_sample, ".", niche_ct_meanprops_id))
+```
+
+### Plotting
+
+```{r, fig.height=16, fig.width=5}
+props_mat <- visium_df %>%
+  dplyr::select(cell_type, niche_id, mean_prop) %>%
+  pivot_wider(values_from = mean_prop,
+              names_from = cell_type) %>%
+  column_to_rownames("niche_id") %>%
+  as.matrix()
+
+#pdf("./processed_visium/integration/niche_compositions/mean_props.pdf",
+#    width = 8, height = 16)
+
+draw(Heatmap(scale(props_mat), cluster_rows = F,
+             row_split = factor(rownames(props_mat) %>% strsplit("[.]") %>% map_chr(., ~.x[2]))))
+
+#dev.off()
+```
+
+
+### What if we filter niches from slides where there's not enough representation?
+
+```{r, fig.height=16, fig.width=5}
+ncells_niche <- 30
+
+props_mat <- visium_df %>%
+  dplyr::filter(ncells > ncells_niche) %>%
+  dplyr::select(cell_type, niche_id, mean_prop) %>%
+  pivot_wider(values_from = mean_prop,
+              names_from = cell_type) %>%
+  column_to_rownames("niche_id") %>%
+  as.matrix()
+
+draw(Heatmap(scale(props_mat), cluster_rows = F,
+             row_split = factor(rownames(props_mat) %>% strsplit("[.]") %>% map_chr(., ~.x[2]))))
+```
+
+```{r, fig.height=16, fig.width=5}
+draw(Heatmap((props_mat), cluster_rows = F, col = col_fun,
+             row_split = factor(rownames(props_mat) %>% strsplit("[.]") %>% map_chr(., ~.x[2]))))
+```
+
+### What if we don't force thee order
+
+```{r, fig.height=16, fig.width=5}
+draw(Heatmap((props_mat), cluster_rows = T, col = col_fun))
+```