library(Seurat)
library(stringr)
library(ggplot2)
library(openxlsx)
library(ggpubr)

#import metadata from  Excel file
Seurat_metadata <- "SingleCell_seurat_metadata.xlsx"
t <- read.xlsx(xlsxFile = Seurat_metadata, sheet = "metadata_Teff_memcells")

# Extract relevant columns for each score analysed (e.g, Activation, Exhaustion, SASP)
t_activation_filt <- t[,c("RNA_Group", "Cluster", "Activation Score1")]
t_sasp_filt <- t[,c("RNA_Group", "Cluster", "SASP Score1")]
t_exhaustion_filt <- t[,c("RNA_Group", "Cluster", "Exhaustion_Score1")]


# Create the activation score plot + ANOVA and T test statistic (exact pvalue)
activation_score_plot<- ggplot(data = t_filt, mapping = aes(x = RNA_Group, y = Activation_Score1, fill = RNA_Group)) + #color = RNA_Group
  theme_bw(base_size = 16) +
  theme(legend.position = "none") +
  xlab("") +
  labs(title="", y="Activation Score", fill="Condition")+
  geom_boxplot(alpha = .5, outlier.size = 0) +
  scale_fill_manual(values = c("FER" = "lightgrey", "OAT" = "khaki", "NOA" = "dodgerblue"))+
  facet_grid(.~Cluster) + #switch = "x"
  stat_compare_means(method = "anova", label.y = 0.45) +
  stat_compare_means(comparisons = list(c("NOA", "OAT"), c("FER", "OAT"),c("FER", "NOA")), exact=T, method = "t.test",tip.length = 0.01, bracket.size = 0.5,
                     p.adjust.method = "bonferroni", label = "p.format", label.y = c(0.30, 0.35, 0.40))  

# Save the plot as a PDF
pdf(file = "activation_score_exact_value.pdf", width = 16, height = 6)
activation_score_plot
dev.off()


# Create the exhaustion score plot + ANOVA and T test statistic (exact pvalue)
exhaustion_score_plot <-ggplot(data = t_exhaustion_filt, mapping = aes(x = RNA_Group, y = Exhaustion_Score1, fill = RNA_Group)) + #color = RNA_Group
  theme_bw(base_size = 16) +
  theme(legend.position = "none") +
  xlab("") +
  labs(title="", y="Exhaustion Score", fill="Condition")+
  geom_boxplot(alpha = .5, outlier.size = 0) +
  scale_fill_manual(values = c("FER" = "lightgrey", "OAT" = "khaki", "NOA" = "dodgerblue"))+
  facet_grid(.~Cluster) + #switch = "x"
  stat_compare_means(method = "anova", label.y = 0.85) +
  stat_compare_means(comparisons = list(c("NOA", "OAT"), c("FER", "OAT"),c("FER", "NOA")), exact=T, method = "t.test",tip.length = 0.01, bracket.size = 0.5,
                     p.adjust.method = "bonferroni", label = "p.format", label.y = c(0.50, 0.60, 0.70)) 

# Save the plot as a PDF
pdf(file = "exhaustion_score_exact_value.pdf", width = 16, height = 6)
exhaustion_score_plot
dev.off()

# Create the sasp score plot + ANOVA and T test statistic (exact pvalue)
sasp_activation_score_plot<- ggplot(data = t_sasp_filt, mapping = aes(x = RNA_Group, y = SASP_Score1, fill = RNA_Group)) + #color = RNA_Group
  theme_bw(base_size = 16) +
  theme(legend.position = "none") +
  xlab("") +
  labs(title="", y="SASP Score", fill="Condition")+
  geom_boxplot(alpha = .5, outlier.size = 0) +
  scale_fill_manual(values = c("FER" = "lightgrey", "OAT" = "khaki", "NOA" = "dodgerblue"))+
  facet_grid(.~Cluster) + #switch = "x"
  stat_compare_means(method = "anova", label.y = 0.75) +
  stat_compare_means(comparisons = list(c("NOA", "OAT"), c("FER", "OAT"),c("FER", "NOA")), exact=T, method = "t.test",tip.length = 0.01, bracket.size = 0.5,
                     p.adjust.method = "bonferroni", label = "p.format", label.y = c(0.45, 0.55, 0.65)) 

# Save the plot as a PDF
pdf(file = "sasp_score_exact_value.pdf", width = 16, height = 6)
sasp_activation_score_plot
dev.off()