流程更新----空间细胞聚类及配受体共现分析（针对visium、bin模式的Stereo-seq、以及HD）

# Loading required packages
library(ISCHIA)
library(robustbase)
library(data.table)
library(ggplot2)
library(Seurat)
library(dplyr)
library(factoextra)
library(cluster)
library(showtext)
library(gridExtra)
library(pdftools)

# Set random seed for reproducibility
set.seed(123)

# Load data
pdac <- readRDS("/path/to/pdac_mets_rctd.rds")
assay_matrix <- pdac[["rctd_tier1"]]@data
norm_weights <- as.data.frame(t(assay_matrix))

# Elbow Method
k.values <- 1:20
wss_values <- sapply(k.values, function(k) kmeans(norm_weights, k, nstart = 10)$tot.withinss)

pdf("1_elbow_plot.pdf")
plot(k.values, wss_values, type = "b", pch = 19, frame = FALSE, 
     xlab = "Number of clusters K", ylab = "Total within-cluster sum of squares",
     main = "Elbow Method for Optimal K")
dev.off()

# Gap Statistic
gap_stat <- function(k) {
  km.res <- kmeans(norm_weights, k, nstart = 10)
  if (k == 1) return(NA)
  obs_disp <- sum(km.res$withinss)
  reference_disp <- mean(replicate(10, {
    km.null <- kmeans(matrix(rnorm(nrow(norm_weights) * ncol(norm_weights)), 
                             ncol = ncol(norm_weights)), k, nstart = 10)
    sum(km.null$withinss)
  }))
  log(reference_disp) - log(obs_disp)
}

gap_stat_values <- sapply(k.values, gap_stat)

pdf("2_gap_statistic_plot.pdf")
plot(k.values, gap_stat_values, type = "b", pch = 19, frame = FALSE, 
     xlab = "Number of Clusters (K)", ylab = "Gap Statistic",
     main = "Gap Statistic: Determining Optimal K")
dev.off()

# Calinski-Harabasz Index
calinski_harabasz_index <- function(data, labels) {
  num_clusters <- length(unique(labels))
  num_points <- nrow(data)
  centroids <- tapply(data, labels, FUN = colMeans)
  between_disp <- sum(sapply(1:num_clusters, function(i) {
    cluster_points <- data[labels == i, ]
    nrow(cluster_points) * sum((colMeans(cluster_points) - centroids[i, ]) ^ 2)
  }))
  within_disp <- sum(sapply(1:num_clusters, function(i) {
    cluster_points <- data[labels == i, ]
    sum((cluster_points - centroids[i, ]) ^ 2)
  }))
  (between_disp / (num_clusters - 1)) / (within_disp / (num_points - num_clusters))
}

ch_values <- sapply(k.values, function(k) {
  km.res <- kmeans(norm_weights, k, nstart = 10)
  calinski_harabasz_index(norm_weights, km.res$cluster)
})

pdf("3_calinski_harabasz_plot.pdf")
plot(k.values, ch_values, type = "b", pch = 19, frame = FALSE,
     xlab = "Number of Clusters (K)", ylab = "Calinski-Harabasz Index",
     main = "Calinski-Harabasz Index: Determining Optimal K")
dev.off()

# ISCHIA Analysis
pdf("4_composition_cluster_k_plot.pdf")
Composition.cluster.k(norm_weights, 20)
dev.off()

pdac <- Composition.cluster(pdac, norm_weights, 12)
pdac$cc_12 <- pdac$CompositionCluster_CC

# Spatial Dimension Plot
image_names <- c("IU_PDA_T1", "IU_PDA_T2", "IU_PDA_HM2", "IU_PDA_HM2_2", "IU_PDA_NP2", 
                 "IU_PDA_T3", "IU_PDA_HM3", "IU_PDA_T4", "IU_PDA_HM4", "IU_PDA_HM5", 
                 "IU_PDA_T6", "IU_PDA_HM6", "IU_PDA_LNM6", "IU_PDA_LNM7", "IU_PDA_T8", 
                 "IU_PDA_HM8", "IU_PDA_LNM8", "IU_PDA_T9", "IU_PDA_HM9", "IU_PDA_T10", 
                 "IU_PDA_HM10", "IU_PDA_LNM10", "IU_PDA_NP10", "IU_PDA_T11", "IU_PDA_HM11", 
                 "IU_PDA_NP11", "IU_PDA_T12", "IU_PDA_HM12", "IU_PDA_LNM12", "IU_PDA_HM13")

paletteMartin <- c('#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#911eb4', 
                   '#46f0f0', '#f032e6', '#bcf60c', '#fabebe', '#008080', '#e6beff')

all_ccs <- unique(pdac$CompositionCluster_CC)
color_mapping <- setNames(paletteMartin[1:length(all_ccs)], all_ccs)

pdf("5_spatial_plots_K12.pdf", width = 10, height = 7)
for (image_name in image_names) {
  plot <- SpatialDimPlot(pdac, group.by = "CompositionCluster_CC", images = image_name) +
    scale_fill_manual(values = color_mapping) +
    theme_minimal() +
    ggtitle(image_name)
  print(plot)
}
dev.off()

# Enriched Cell Types
save_cc_plot <- function(cc) {
  plot <- Composition_cluster_enrichedCelltypes(pdac, cc, as.matrix(norm_weights))
  pdf_name <- paste0(cc, ".pdf")
  pdf(file = pdf_name)
  print(plot)
  dev.off()
}

ccs <- paste0("CC", 1:12)
for (cc in ccs) {
  save_cc_plot(cc)
}

pdf_files <- paste0("CC", 1:12, ".pdf")
pdf_combine(pdf_files, output = "6_enrichedCelltypes_CC_12.pdf")

# UMAP
pdac.umap <- Composition_cluster_umap(pdac, norm_weights)
pdf("7_umap_pie_chart.pdf")
print(pdac.umap$umap.deconv.gg)
dev.off()

# Add UMAP to Seurat object
emb.umap <- pdac.umap$umap.table
emb.umap$CompositionCluster_CC <- NULL
emb.umap$Slide <- NULL
emb.umap <- as.matrix(emb.umap)
colnames(emb.umap) <- c("UMAP1", "UMAP2")

pdac[['umap.ischia12']] <- CreateDimReducObject(embeddings = emb.umap, key = 'umap.ischia12_', assay = 'rctd_tier1')

pdf("8_seurat_ischia_umap_12.pdf")
DimPlot(pdac, reduction = "umap.ischia12", label = FALSE, group.by="cc_12")
dev.off()

# Bar plots
pdf("9_barplot_SampVsorig_12.pdf", height=12, width=20)
dittoBarPlot(pdac, "orig.ident", group.by = "cc_12")
dev.off()

pdf("10_barplot_origVsSamp_12.pdf", height=10, width=20)
dittoBarPlot(pdac, "cc_12", group.by = "orig.ident")
dev.off()

# Cell type co-occurrence
CC4.celltype.cooccur <- spatial.celltype.cooccurence(spatial.object=pdac, deconv.prob.mat=norm_weights, 
                                                     COI="CC4", prob.th= 0.05, 
                                                     Condition=unique(pdac$orig.ident))
pdf("11_celltype_cooccurrence_CC4.pdf")
plot.celltype.cooccurence(CC4.celltype.cooccur)
dev.off()

### This script performs the wilcoxon rank sum test and hierarchical clustering on the RCTD tier data to identify the significant abundant cell types between the clusters.
  
#### Load necessary packages
 ```{r}
library(Seurat)
library(compositions)
library(tidyverse)
library(clustree)
library(patchwork)
library(uwot)
library(scran)
library(cluster)
library(ggrastr)
library(cowplot)
# library(conflicted) # to be loaded in case of a conflict arises.
 
config <- config::get()
 
# source(here::here("pdac_nac", "visualization", "eda.R"))
 
### Load the seurat object and get the proportions data
so <- readRDS(here::here(config$data_processed, "06-pdac_CC10_msig.rds"))
 
# Join with metadata if needed
metadata <- so@meta.data %>%
  select(orig.ident, patient_id, neoadjuvant_chemo, CompositionCluster_CC) %>%
  rownames_to_column("row_id")
 
 
# Get the proportions data
rctd_tier2 <- t(so@assays$rctd_tier2@data)
 
# Ensure the data is in the right format
rownames(rctd_tier2) <- make.unique(rownames(rctd_tier2))
 
# Log transformation of rctd_tier1
log_comps <- log10(rctd_tier2)
## Perform the summary statistics
We perform the summary statistics for the RCTD tier data. We perform hierarchical clustering and do the wilcoxon rank sum test to identify the differentially abundant cell types between the clusters. 
#### Prepare the data for the summary statistics
# Prepare data for summary statistics
cluster_summary_pat <- rctd_tier2 %>%
  as.data.frame() %>%
  rownames_to_column("row_id") %>%
  left_join(metadata, by = "row_id") %>%  # Join with meta_data using row_id as the key
  pivot_longer(-c(row_id, orig.ident, patient_id, neoadjuvant_chemo, CompositionCluster_CC), values_to = "ct_prop", names_to = "cell_type") %>%
  group_by(orig.ident, patient_id, neoadjuvant_chemo, CompositionCluster_CC, cell_type) %>%
  summarize(median_ct_prop = median(ct_prop, na.rm = TRUE))

# Aggregate data for median ct prop
cluster_summary <- cluster_summary_pat %>%
  ungroup() %>%
  group_by(CompositionCluster_CC, cell_type) %>%
  summarize(patient_median_ct_prop = median(median_ct_prop, na.rm = TRUE))
 
# Prepare matrix for hierarchical clustering
cluster_summary_mat <- cluster_summary %>%
  pivot_wider(values_from = patient_median_ct_prop, names_from = cell_type, values_fill = list(patient_median_ct_prop = 0)) %>%
  column_to_rownames("CompositionCluster_CC") %>%
  as.matrix()
 
# conflicts_prefer(stats::"dist") # resolve conflicts between %*% functions

# Perform hierarchical clustering
cluster_order <- hclust(dist(cluster_summary_mat))$labels[hclust(dist(cluster_summary_mat))$order] # use this if you want to order the clusters based on the hierarchical clustering
ct_order <- hclust(dist(t(cluster_summary_mat)))$labels[hclust(dist(t(cluster_summary_mat)))$order]
 
# Order Clusters in ascending order
# cluster_order1 <- c("CC1", "CC2", "CC3", "CC4", "CC5", "CC6", "CC7", "CC8", "CC9", "CC10")
 
# Wilcoxon test for characteristic cell types
run_wilcox_up <- function(prop_data) {
  prop_data_group <- prop_data[["CompositionCluster_CC"]] %>% unique() %>% set_names()
  map(prop_data_group, function(g) {
    test_data <- prop_data %>%
      mutate(test_group = ifelse(CompositionCluster_CC == 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("CompositionCluster_CC") %>% unnest()
}

 
wilcoxon_res <- cluster_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)) %>%
  mutate(significant = ifelse(p_corr <= 0.15, "*", ""))

 
#### Save the summary statistics

# give the path to save the summary statistics
file_path_cluster_summ <- here::here(config$data_interim,  "summary_of_clusters.txt")
file_path_wilcox_res <- here::here(config$data_interim, ß "wilcoxon_res_cells_clusters.txt")
 
# Save the summary statistics
write.table(cluster_summary_pat, file = file_path_cluster_summ, col.names = TRUE, row.names = FALSE, quote = FALSE, sep = "\t")
write.table(wilcoxon_res, file = file_path_wilcox_res, col.names = TRUE, row.names = FALSE, quote = FALSE, sep = "\t")

 
#### Plot the summary statistics
 
# Plotting mean ct prop and barplots
mean_ct_prop_plt <- cluster_summary %>%
  left_join(wilcoxon_res, by = c("CompositionCluster_CC", "cell_type")) %>%
  mutate(cell_type = factor(cell_type, levels = ct_order), CompositionCluster_CC = factor(CompositionCluster_CC, levels = cluster_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 = CompositionCluster_CC, fill = scaled_pat_median)) +
  geom_tile(color = "black") +
  geom_text(aes(label = significant)) +
  theme(axis.text.x = element_text(angle = 90, 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()
 
cluster_counts <- cluster_info %>%
  dplyr::select_at(c("row_id", "CompositionCluster_CC")) %>%
  group_by(CompositionCluster_CC) %>%
  summarize(nspots = length(CompositionCluster_CC)) %>%
  mutate(prop_spots = nspots / sum(nspots))
 
file_path_cluster_prop_summ <- here::here(config$data_interim, "cluster_prop_summary.csv")
 
write_csv(cluster_counts, file_path_cluster_prop_summ)
 
#barplots for cluster counts
barplts <- cluster_counts %>%
  mutate(CompositionCluster_CC = factor(CompositionCluster_CC, levels = cluster_order)) %>%
  ggplot(aes(y = CompositionCluster_CC, 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))
 
cluster_summary_plt <- cowplot::plot_grid(mean_ct_prop_plt, barplts, align = "hv", axis = "tb") # use if barplots are needed to show the spot counts otherwise directly use mean_ct_prop_plt for the plot

 
#### plot the summary clusters
 
pdf_path_summ_clust <- here::here(config$plots, "wilcox_summary_clusters.pdf")
 
pdf(pdf_path_summ_clust, width = 20, height = 10)
plot(cluster_summary_plt)
dev.off()
 
#box plot for median ct prop
pdf_path_boxplot <- here::here(config$plots,  "wilcox_bboxplot_median_ct_prop.pdf")
 
plt <- cluster_summary_pat %>%
  ggplot(aes(x = CompositionCluster_CC, y = median_ct_prop)) +
  geom_boxplot() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  facet_wrap(. ~ cell_type, ncol = 3, scales = "free_y")
 
pdf(pdf_path_boxplot, width = 20, height = 10)
plot(plt)
dev.off()