rm(list = ls())
setwd("/raid1/shiq/cellModule/normal/08_st/06_xenium")

library(dplyr)
library(ggplot2)
library(patchwork)
library(ggpubr)
library(ComplexHeatmap)
library(RColorBrewer)
library(circlize)
library(ggsci)
library(viridis)
library(extrafont)
# font_import()
# fonts()
color <- readRDS("/raid1/shiq/resource/color/color_list.rds")

## prepare
m_ls <- list()
for (i in 1:6) {
    cm02 <- read.csv(paste0("./Xenium_V1_hColon_Non_diseased_Add_on_FFPE/Selection_", i, "_cells_stats_CM02.csv"), skip = 2, header = T)
    cm03 <- read.csv(paste0("./Xenium_V1_hColon_Non_diseased_Add_on_FFPE/Selection_", i, "_cells_stats_CM03.csv"), skip = 2, header = T)

    colnames(cm02) <- c("cellID", "cluster", "cm02", "area")
    colnames(cm03) <- c("cellID", "cluster", "cm03", "area")

    m <- merge(cm02, cm03)
    m <- m %>%
        group_by(cluster) %>%
        summarise_if(is.numeric, .funs = sum) %>%
        as.data.frame()

    m$cm02ave <- m$cm02 / m$area
    m$cm03ave <- m$cm03 / m$area
    sampleID <- paste0("Selection_", i)
    m$sampleID <- sampleID
    m_ls[[sampleID]] <- m
}

df <- Reduce(f = rbind, x = m_ls)

color$xenium <- c(
    "Epi" = rgb(29 / 255, 72 / 255, 122 / 255),
    "LP" = rgb(86 / 255, 86 / 255, 86 / 255)
)

## plot
ggplot(df, aes(x = cluster, y = cm02ave * 1e2, color = cluster)) +
    geom_boxplot(outlier.shape = NA, linewidth = 0.3, alpha = 0) +
    geom_jitter(size = 1, width = 0.3, shape = 16) +
    # geom_line(aes(group = sampleID)) +
    scale_color_manual(values = color$xenium) +
    labs(x = "", y = "Transcript density\n(1e-2 per um2)", title = "CM02") +
    # coord_cartesian(ylim = c(0, 9)) +
    stat_compare_means(
        # Use only p.format as label. Remove method name.
        aes(label = paste0("P = ", after_stat(p.format))),
        method = "wilcox.test", size = 2, label.y.npc = 1
    ) +
    theme_classic() +
    theme(
        line = element_line(linewidth = 0.3),
        text = element_text(size = 6, family = "Arial", color = "black"),
        axis.text = element_text(size = 6, family = "Arial", color = "black"),
        # axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(size = 7, hjust = 0.5),
        legend.position = "none",
        strip.text.x = element_text(size = 6),
        strip.background = element_blank(),
        aspect.ratio = 1.5 / 1
    ) -> p1

ggplot(df, aes(x = cluster, y = cm03ave * 1e2, color = cluster)) +
    geom_boxplot(outlier.shape = NA, linewidth = 0.3, alpha = 0) +
    geom_jitter(size = 1, width = 0.3, shape = 16) +
    # geom_line(aes(group = sampleID)) +
    scale_color_manual(values = color$xenium) +
    labs(x = "", y = "Transcript density\n(1e-2 per um2)", title = "CM03") +
    # coord_cartesian(ylim = c(0, 7.5)) +
    stat_compare_means(
        # Use only p.format as label. Remove method name.
        aes(label = paste0("P = ", after_stat(p.format))),
        method = "wilcox.test", size = 2, label.y.npc = 1
    ) +
    theme_classic() +
    theme(
        line = element_line(linewidth = 0.3),
        text = element_text(size = 6, family = "Arial", color = "black"),
        axis.text = element_text(size = 6, family = "Arial", color = "black"),
        # axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(size = 7, hjust = 0.5),
        legend.position = "none",
        strip.text.x = element_text(size = 6),
        strip.background = element_blank(),
        aspect.ratio = 1.5 / 1
    ) -> p2

plt <- p1 + p2 + plot_layout(nrow = 1)
ggsave(filename = "./Figure3e.pdf", plot = plt, width = 3, height = 2)