课前准备---空间转录组数据分析之分子niche

#! usr/R
####zhaoyunfei
####20240716

suppressMessages({
library(Seurat)
library(compositions)
library(tidyverse)
library(clustree)
library(uwot)
library(cluster)
library(RColorBrewer)
})

adata = readRDS('Muscle.spatial.rds')

cluster_info <- as.data.frame(adata$seurat_clusters)

colnames(cluster_info) = 'mol_niche'

cluster_info$row_id = rownames(cluster_info)

integrated_compositions <- t(adata@assays$predictions@data)

integrated_compositions = integrated_compositions[,-which(colnames(integrated_compositions) %in% c("max",'MTJ','RegMyon'))]

niche_summary_pat <- integrated_compositions %>% as.data.frame() %>% rownames_to_column("row_id") %>% pivot_longer(-row_id,values_to = "ct_prop", names_to = "cell_type") %>% left_join(cluster_info) %>% mutate(orig.ident = strsplit(row_id, "[..]") %>% map_chr(., ~ .x[1])) %>% group_by(orig.ident, mol_niche, cell_type) %>% summarize(median_ct_prop = median(ct_prop))

niche_summary <- niche_summary_pat %>% ungroup() %>% group_by(mol_niche, cell_type) %>% summarise(patient_median_ct_prop = median(median_ct_prop))

niche_summary$mol_niche = paste('niche',niche_summary$mol_niche,sep = '_')

niche_summary_pat %>% ungroup() %>% group_by(mol_niche, cell_type) %>% summarise(patient_median_ct_prop = median(median_ct_prop))

niche_summary_mat <- niche_summary %>% pivot_wider(values_from = patient_median_ct_prop, names_from =  cell_type, values_fill = 0) %>% column_to_rownames("mol_niche") %>% as.matrix()

niche_order <- hclust(dist(niche_summary_mat))

niche_order <- niche_order$labels[niche_order$order]

ct_order <- hclust(dist(t(niche_summary_mat)))

ct_order <- ct_order$labels[ct_order$order]

# Find characteristic cell types of each niche
# We have per patient the proportion of each cell-type in each niche

run_wilcox_up <- function(prop_data) {
  
    prop_data_group <- prop_data[["mol_niche"]] %>% unique() %>% set_names()
  
    map(prop_data_group, function(g) {
    
    test_data <- prop_data %>% mutate(test_group = ifelse(mol_niche == g,"target", "rest")) %>% mutate(test_group = factor(test_group,levels = c("target", "rest")))
    
    wilcox.test(median_ct_prop ~ test_group, data = test_data, alternative = "greater") %>% broom::tidy()
  }) %>% enframe("mol_niche") %>% unnest()
  
}

wilcoxon_res <- niche_summary_pat %>%
  ungroup() %>%
  group_by(cell_type) %>%
  nest() %>%
  mutate(wres = map(data, run_wilcox_up)) %>%
  dplyr::select(wres) %>%
  unnest() %>%
  ungroup() %>%
  mutate(p_corr = p.adjust(p.value))

wilcoxon_res <- wilcoxon_res %>% mutate(significant = ifelse(p_corr <= 0.15, "*", ""))

wilcoxon_res$mol_niche = paste('niche',wilcoxon_res$mol_niche,sep = '_')

####head(niche_summary_pat)

mean_ct_prop_plt <- niche_summary %>%
  left_join(wilcoxon_res, by = c("mol_niche", "cell_type")) %>%
  mutate(cell_type = factor(cell_type, levels = ct_order),
         mol_niche = factor(mol_niche, levels = niche_order)) %>%
  ungroup() %>%
  group_by(cell_type) %>%
  mutate(scaled_pat_median = (patient_median_ct_prop - mean(patient_median_ct_prop))/sd(patient_median_ct_prop)) %>%
  ungroup() %>%
  ggplot(aes(x = cell_type, y = mol_niche, fill = scaled_pat_median)) +
  geom_tile(color="black",size=0.5) +
  geom_text(aes(label = significant)) +
  theme(axis.text.x = element_text(angle = 180, hjust = 1, vjust = 0.5, size = 12),
        legend.position = "bottom",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.text.y = element_text(size=12)) +
  scale_fill_gradient2(high = "red",mid = 'white', low = "blue") + theme_classic()


# Finally describe the proportions of those niches in all the data
cluster_counts <- cluster_info %>%
  dplyr::select_at(c("row_id", "mol_niche")) %>%
  group_by(mol_niche) %>%
  summarise(nspots = length(mol_niche)) %>%
  mutate(prop_spots = nspots/sum(nspots))

write_csv(cluster_counts, file = "./mol_niche_prop_summary.csv")

barplts <- cluster_counts %>%
  mutate(ct_niche = factor(mol_niche, levels = niche_order)) %>%
  ggplot(aes(y = mol_niche, x = prop_spots)) +
  geom_bar(stat = "identity") +
  theme_classic() + ylab("") +
  theme(axis.text.y = element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.text.x = element_text(size=12))

niche_summary_plt <- cowplot::plot_grid(mean_ct_prop_plt, barplts, align = "hv", axis = "tb")

pdf('characteristic_mol_niches.pdf',width = 15,height = 5.5)

plot(niche_summary_plt)

dev.off()

adata = adata[,cluster_info$row_id]

Idents(adata) = cluster_info$ct_niche

sample_cols = unique(c(brewer.pal(8,"Dark2"),brewer.pal(9,"Set1"),brewer.pal(8,"Set2"),brewer.pal(12,"Set3")))
defined_cols = c('#e6194b', '#4363d8', '#f58231', '#911eb4', '#46f0f0', '#f032e6', '#bcf60c', '#fabebe', '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000', '#aaffc3', '#808000', '#ffd8b1', '#000075', '#808080', sample_cols)

p = SpatialDimPlot(adata, label = TRUE, label.size = 3)

for (i in 1:length(p)){p[[i]] = p[[i]] + guides(fill=guide_legend(override.aes = list(size=5))) + scale_fill_manual(values = defined_cols) }

pdf('niche.spatial.pdf',width = 15)

print(p)

dev.off()