## ----setup, include=FALSE-----------------------------------------------------
# Load required libraries
library(nsprcomp)
library(ggplot2)
library(ggtext)
library(dplyr)
library(gridExtra)
library(ggpubr)
library(pheatmap)
library(grid)
library(circlize)
library(paletteer)
library(ggthemes)
library(ComplexHeatmap)
library(EMTscore)
library(AUCell)
library(GSA)
library(mclust)
library(ggalluvial)
library(RColorBrewer)
# Single-cell example data and object infrastructure
library(EMTscoreData)
library(ExperimentHub)
library(SingleCellExperiment)
library(SummarizedExperiment)
library(Seurat)

# Use a single core so the vignette is reproducible on all builders.
# Save the user's original option values so they can be restored at the end.
old_opts <- options(mc.cores = 1)
# Use the cairo bitmap backend for headless rendering when available
if (capabilities("cairo")) old_opts <- c(old_opts, options(bitmapType = "cairo"))

# Global chunk options
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>",
  dev = "png", dev.args = list(type = "cairo"),
  fig.width = 8, fig.height = 6, dpi = 100, out.width = "100%"
)

## ----install, eval=FALSE------------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("EMTscore")

## ----emtscore-2---------------------------------------------------------------
data(cell_annotation_file)
head(cell_annotation_file)

## ----emtscore-3---------------------------------------------------------------
bfc <- BiocFileCache::BiocFileCache()
geneExp_path <- BiocFileCache::bfcrpath(
  bfc, "https://zenodo.org/records/19487376/files/geneExp.rda"
)
load(geneExp_path) # loads the full `geneExp` matrix (15950 genes x 120 samples)
head(geneExp[, 1:5], 4)
dim(geneExp)

## ----emtscore-4---------------------------------------------------------------
data("Panchy_et_al_E_signature", package = "EMTscore")
data("Panchy_et_al_M_signature", package = "EMTscore")
head(Panchy_et_al_M_signature)

## ----emtscore-5---------------------------------------------------------------
gene_sets <- list(
Panchy_et_al_E_signature = Panchy_et_al_E_signature$GeneName,
Panchy_et_al_M_signature = Panchy_et_al_M_signature$GeneName
)
write_gmt(gene_sets, file.path(tempdir(), "EM_signature.gmt"))

## ----emtscore-6---------------------------------------------------------------
gmt_file <- system.file("extdata", "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION.v2025.1.Hs.gmt", package = "EMTscore")
nnPCA_Result_EMT <- Execute_nnPCA(geneExp, gmt_file,dimension = 1, score_names = "EMT")
AUCell_Result_EMT <- Execute_AUCell(geneExp, gmt_file, score_names = "EMT")
GSVA_Result_EMT <- Execute_GSVA(geneExp, gmt_file, score_names = "EMT")
ssGSEA_Result_EMT <- Execute_ssGSEA(geneExp, gmt_file, score_names = "EMT")
JASMINE_Result_EMT <- Execute_JAS(geneExp, gmt_file, score_names = "EMT")
SCSE_Result_EMT <- Execute_SCSE(geneExp, gmt_file, score_names = "EMT")

## ----emtscore-7---------------------------------------------------------------
gmt_file <- system.file("extdata", "EM_signature.gmt", package = "EMTscore")
nnPCA_Result_multiple <- Execute_nnPCA_parallel(geneExp, gmt_file, dimension = 1, cores = 1)
AUCell_Result_multiple <- Execute_AUCell_parallel(geneExp, gmt_file, cores = 1)
ssGSEA_Result_multiple <- Execute_ssGSEA_parallel(geneExp, gmt_file, cores = 1)
GSVA_Result_multiple <- Execute_GSVA_parallel(geneExp, gmt_file, cores = 1)
JASMINE_Result_multiple <- Execute_JASMINE_parallel(geneExp, gmt_file, cores = 1)
SCSE_Result_multiple <- Execute_SCSE_parallel(geneExp, gmt_file, cores = 1)

## ----emtscore-8---------------------------------------------------------------
data_for_plot <- data_prepare(cell_annotation_file, nnPCA_Result_multiple, merge_colname = "name")

## ----emtscore-9---------------------------------------------------------------
plot1 <- Execute_E_M_plot(
data_for_plot,
E_colname = "Panchy_et_al_E_signature",
M_colname = "Panchy_et_al_M_signature",
celltype_colname = "celltype_annotation",
colors = c("#F87189", "#CE9031", "#A48CF5", "#97A430", "#39A7D0", "#E57D5F",
"#84C7B9", "#E1AF64", "#C26CCF", "#B0BF43", "#57C3E8", "#F29D9E", "#92AAE6")
)
plot1

## ----emtscore-10--------------------------------------------------------------
method_list <- list(
  nnPCA = nnPCA_Result_multiple,
  AUCell = AUCell_Result_multiple,
  ssGSEA = ssGSEA_Result_multiple,
  JASMINE = JASMINE_Result_multiple,
  SCSE = SCSE_Result_multiple
)

plot_list <- lapply(names(method_list), function(method_name) {
  data_for_plot <- data_prepare(
    cell_annotation_file,
    method_list[[method_name]],
    merge_colname = "name"
  )
  p <- Execute_E_M_plot(
    data_for_plot,
    E_colname = "Panchy_et_al_E_signature",
    M_colname = "Panchy_et_al_M_signature",
    celltype_colname = "celltype_annotation",
    colors = c("#F87189", "#CE9031", "#A48CF5", "#97A430", "#39A7D0", "#E57D5F",
               "#84C7B9", "#E1AF64", "#C26CCF", "#B0BF43", "#57C3E8",
               "#F29D9E", "#92AAE6")
  )
  p +  labs(x = "E score", y = "M score") + ggtitle(method_name)# 自动添加方法名称作为主标题
})

names(plot_list) <- names(method_list)
combined_plot <- Arrange_plots(plot_list, ncol_per_row = 3)
combined_plot

## ----emtscore-11--------------------------------------------------------------
# Prepare M gene set and compute nnPCA
data(cell_annotation_file)
gene_set <- list(Panchy_et_al_M_signature = Panchy_et_al_M_signature$GeneName)
m_gmt <- system.file("extdata", "M_signature.gmt", package = "EMTscore")

nnPCA_Mscore <- Execute_nnPCA(geneExp, m_gmt, dimension=2, score_names=c('M1_score','M2_score'))
data_for_plot <- data_prepare(cell_annotation_file, nnPCA_Mscore, merge_colname = "name")

plot2 <- Execute_M_dimension_plot(
data_for_plot = data_for_plot,
M1_colname = "M1_score",
M2_colname = "M2_score",
celltype_colname = "celltype_annotation",
colors = c("#F87189", "#CE9031", "#A48CF5", "#97A430", "#39A7D0", "#E57D5F",
"#84C7B9", "#E1AF64", "#C26CCF", "#B0BF43", "#57C3E8", "#F29D9E", "#92AAE6")
)
plot2

## ----emtscore-12--------------------------------------------------------------
combined_plot <- Arrange_plots(
plots_list = list(plot1, plot2),
ncol_per_row = 2,
subtitles = c("E vs M", "M1 vs M2"),
fig_title = "Panchy_et_al"
)
print(combined_plot)

## ----emtscore-13--------------------------------------------------------------
colors <- c("#F87189", "#CE9031", "#A48CF5", "#97A430", "#39A7D0", "#E57D5F",
"#84C7B9", "#E1AF64", "#C26CCF", "#B0BF43", "#57C3E8", "#F29D9E", "#92AAE6")
p_hist <- data_for_plot %>%
     ggplot( aes(x=M1_score, fill=celltype_annotation)) +
     geom_histogram(alpha=0.6, position = 'identity') +
     scale_fill_manual(values=colors) +
     theme_classic() +
     labs(fill="")
p_hist

## ----emtscore-14--------------------------------------------------------------
# Example heatmap
plot_heatmap_function(t(geneExp), Panchy_et_al_M_signature)

## ----emtscore-load-sc---------------------------------------------------------
eh <- ExperimentHub::ExperimentHub()
query(eh, "EMTscoreData")

A549_TNF   <- eh[["EH10291"]]
A549_EGF   <- eh[["EH10292"]]
A549_TGFB1 <- eh[["EH10293"]]

## ----emtscore-scobjects-------------------------------------------------------
set.seed(1)
subset_cells <- function(sce, n = 1000) {
  n <- min(n, ncol(sce))
  sce[, sample(seq_len(ncol(sce)), n)]
}
A549_TGFB1 <- subset_cells(A549_TGFB1, n = 1000)
A549_EGF   <- subset_cells(A549_EGF,   n = 1000)
A549_TNF   <- subset_cells(A549_TNF,   n = 1000)

# Organize the datasets into a named list. Each element corresponds to one
# experimental condition. The EMTscore functions accept SingleCellExperiment
# objects directly and convert them to Seurat objects internally, so no manual
# conversion is required.
objects <- list(
  A549_TGFB1 = A549_TGFB1,
  A549_EGF   = A549_EGF,
  A549_TNF   = A549_TNF
)

## ----emtscore-15, message=FALSE-----------------------------------------------
gmt_file <- system.file("extdata", "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION.v2025.1.Hs.gmt", package = "EMTscore")

seurat_objs <- add_EMT_score(objects, gmt_file = gmt_file, emt_name = "EMT_score", method = "nnPCA", nnPCA_dim = 1)
p_nnPCA <- plot_EMT_from_objects(seurat_objs, col_name = "Pseudotime", emt_score_col = "EMT_score")
p_nnPCA

## ----emtscore-16, message=FALSE-----------------------------------------------
gmt_file <- system.file("extdata", "EM_signature.gmt", package = "EMTscore")

emt_names <- c("Escore", "Mscore")
result <- add_EMT_score_multiple(objects, gmt_file, emt_names, method = "nnPCA", nnPCA_dim = 1, cores = 1)

plot_all_clusters <- function(result, method = c("Kmeans", "GMM"), emt_names, n_clusters = 3) {
  method <- match.arg(method)

  for (name in names(result)) {
    message("Processing: ", name)
    obj <- result[[name]]

    # Extract E and M signature data
    sig_df <- obj[[emt_names]]
    colnames(sig_df) <- c("Escore", "Mscore")

    # Run clustering
    if (method == "GMM") {
      cl <- predict_cluster_labels(sig_df, method = "GMM", n_clusters = 3, PC_name = c("Escore", "Mscore"))
    } else if (method == "Kmeans") {
      cl <- predict_cluster_labels(sig_df, method = "Kmeans", n_clusters = 3, PC_name = c("Escore", "Mscore"))
    }

    # Scatter plot data
    scatter_df <- data.frame(sig_df, Cluster = as.factor(cl))
    p1 <- ggplot(scatter_df, aes(x = Escore, y = Mscore, color = Cluster)) +
      geom_point(size = 3, alpha = 0.8) +
      scale_color_brewer(palette = "Set1") +
      theme_classic(base_size = 14) +
      labs(title = paste0(name, " - ", method, " clustering"), x = "Escore", y = "Mscore")
    print(p1)

    # Sankey data
    cl_df <- data.frame(cell = rownames(sig_df), Cluster = as.character(unlist(cl)))
    true_df <- data.frame(cell = colnames(obj),
                          TrueLabel = as.character(unlist(obj$Time)))
    df_merge <- inner_join(cl_df, true_df, by = "cell")
    df_count <- df_merge %>%
      group_by(Cluster, TrueLabel) %>%
      summarise(Freq = n(), .groups = "drop")

    # Sankey plot
    df_count$TrueLabel <- factor(
      df_count$TrueLabel,
      levels = c("0d", "8h", "1d", "3d", "7d", "8h_rm", "1d_rm", "3d_rm")
      )
    p2 <- ggplot(df_count, aes(axis1 = Cluster, axis2 = TrueLabel, y = Freq)) +
      geom_alluvium(aes(fill = Cluster), width = 1/12) +
      geom_stratum(width = 1/8, fill = "grey90", color = "black") +
      geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
      scale_x_discrete(limits = c(method, "TrueLabel"), expand = c(.1, .1)) +
      scale_fill_brewer(palette = "Set1") +
      theme_minimal(base_size = 14) +
      labs(title = paste("Sankey Diagram:", method, "Cluster vs True Cell Labels -", name),
           y = "Number of Cells", x = "")
    print(p2)
  }
}

plot_all_clusters(result, method = "GMM", emt_names)

## ----emtscore-17, message=FALSE, warning=FALSE, results='hide'----------------
filtered_file <- system.file("extdata", "filtered.c2.gmt", package = "EMTscore")

nnPCA_Result_multiple <- Execute_nnPCA_parallel(geneExp, filtered_file, dimension = 1, cores = 1)
result <- correlate_sample_scores(score_mat1 = nnPCA_Result_multiple, score_mat2 = nnPCA_Result_EMT, method = "pearson")
head(result)

## ----emtscore-18--------------------------------------------------------------
df_top10 <- head(result, 10)

df_top10$Pathway <- factor(df_top10$Pathway_in_score_mat1,
                           levels = df_top10$Pathway_in_score_mat1)

x <- as.character(df_top10$Pathway)
x <- tolower(gsub("_", " ", x))
df_top10$Pathway <- paste0(toupper(substr(x, 1, 1)), substr(x, 2, nchar(x)))


ggplot(df_top10, aes(x = reorder(Pathway, Correlation), y = Correlation)) +
  geom_col(fill = "#84C7B9", width = 0.8) +
  geom_richtext(
    aes(label = Pathway, y = 0),
    hjust = 0,
    fill = NA,
    label.color = NA,
    colour = "black",
    size = 5,
    label.padding = unit(0.15, "lines"),
    background_colour = "black"
  ) +
  coord_flip() +
  labs(title = "Top 10 Positive Correlated Pathways", x = "Pathways", y = "Correlation") +
  theme_bw(base_size = 12) +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.grid.major.y = element_blank(),
    axis.title.x = element_text(color = "black"),
    axis.title.y = element_text(color = "black"),
    axis.text.x  = element_text(color = "black")
  )+
  theme(axis.ticks.length.x = unit(0.15, 'cm'),
           axis.ticks.length.y = unit(0.15, 'cm'),axis.title = element_text(size = 16),
           axis.text  = element_text(size = 14))

## ----emtscore-19--------------------------------------------------------------
sig_result <- subset(result, P_value < 0.05)
df_top_N10 <- head(sig_result[order(sig_result$Correlation), ], 10)

df_top_N10$Pathway <- factor(df_top_N10$Pathway_in_score_mat1,
                             levels = df_top_N10$Pathway_in_score_mat1)

x <- as.character(df_top_N10$Pathway)
x <- tolower(gsub("_", " ", x))
df_top_N10$Pathway <- paste0(toupper(substr(x, 1, 1)), substr(x, 2, nchar(x)))

ggplot(df_top_N10, aes(x = reorder(Pathway, Correlation), y = Correlation)) +
  geom_col(fill = "#84C7B9", width = 0.8) +
  geom_richtext(
    aes(label = Pathway, y = 0),
    hjust = 1,
    fill = NA,
    label.color = NA,
    colour = "black",
    size = 5,
    label.padding = unit(0.15, "lines"),
    background_colour = "black"
  ) +
  coord_flip() +
  labs(title = "Top 10 Negative Correlated Pathways", x = "Pathways", y = "Correlation") +
  theme_bw(base_size = 12) +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.grid.major.y = element_blank(),
    axis.title.x = element_text(color = "black"),
    axis.title.y = element_text(color = "black"),
    axis.text.x  = element_text(color = "black")
  ) +
  theme(axis.ticks.length.x = unit(0.15, 'cm'),
           axis.ticks.length.y = unit(0.15, 'cm'),axis.title = element_text(size = 16),
           axis.text  = element_text(size = 14))

## ----emtscore-20--------------------------------------------------------------
signature_file <- system.file("extdata", "stemsig.tsv", package = "EMTscore")
scores <- compute_Signature_score(geneExp, signature_file , score_name = "stemness_score")
head(scores)

## ----emtscore-21--------------------------------------------------------------
signature_file <- system.file("extdata", "cellular_senescence_sig.tsv", package = "EMTscore")
scores <- compute_Signature_score(geneExp, signature_file, score_name = "senescence_score")
head(scores)

## ----emtscore-22--------------------------------------------------------------
signature_file1 <- system.file("extdata", "stemsig.tsv", package = "EMTscore")
signature_file2 <- system.file("extdata", "cellular_senescence_sig.tsv", package = "EMTscore")

signature_files <- c(signature_file1, signature_file2)

result <- compute_Signature_score_SingleCell(objects, signature_files, score_name = c("Stemness_Score", "Senescence_Score"))
head(result$A549_TGFB1[[]])

## ----emtscore-23--------------------------------------------------------------
gmt_file <- system.file("extdata", "EM_signature.gmt", package = "EMTscore")

EMscore_result <- add_EMT_score_multiple(objects, gmt_file,emt_names = c("Escore", "Mscore"),
                          method = "nnPCA",
                          nnPCA_dim = 1, cores = 1)

head(EMscore_result$A549_TGFB1[[]])

gmt_file <- system.file("extdata", "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION.v2025.1.Hs.gmt", package = "EMTscore")

Mscore_result <- add_EMT_score(objects, gmt_file,emt_name = c("Mscore_PC1", "Mscore_PC2"),
                          method = "nnPCA",
                          nnPCA_dim = 2)

head(Mscore_result$A549_TGFB1[[]])
plot_EMT_from_objects(Mscore_result, col_name = "Pseudotime", emt_score_col = "Mscore_PC1")
plot_EMT_from_objects(Mscore_result, col_name = "Pseudotime", emt_score_col = "Mscore_PC2")

## ----emtscore-24--------------------------------------------------------------
df <- EMscore_result$A549_TGFB1[[]]
plot1 <- Execute_E_M_plot(
df,
E_colname = "Escore",
M_colname = "Mscore",
celltype_colname = "Time",
colors = c("#F87189", "#CE9031", "#A48CF5", "#97A430", "#39A7D0", "#E57D5F",
"#84C7B9", "#E1AF64", "#C26CCF", "#B0BF43", "#57C3E8", "#F29D9E", "#92AAE6")
)
plot1

## ----emtscore-25--------------------------------------------------------------
df <- Mscore_result$A549_TGFB1[[]]
plot2 <- Execute_M_dimension_plot(
df,
M1_colname = "Mscore_PC1",
M2_colname = "Mscore_PC2",
celltype_colname = "Time",
colors = c("#F87189", "#CE9031", "#A48CF5", "#97A430", "#39A7D0", "#E57D5F",
"#84C7B9", "#E1AF64", "#C26CCF", "#B0BF43", "#57C3E8", "#F29D9E", "#92AAE6")
)
plot2

## ----emtscore-26--------------------------------------------------------------
combined_plot <- Arrange_plots(
plots_list = list(plot1, plot2),
ncol_per_row = 2,
subtitles = c("E vs M", "M1 vs M2"),
fig_title = "Cook_et_al"
)
print(combined_plot)

## ----emtscore-27--------------------------------------------------------------
df <- result$A549_TGFB1[[]]
cor_test <- cor.test(df$Senescence_Score, df$Stemness_Score, method = "pearson")
R_val <- round(cor_test$estimate, 2)
p_val <- format(cor_test$p.value, scientific = TRUE, digits = 2)
n_clust <- length(unique(df$Time))
palette_colors <- colorRampPalette(RColorBrewer::brewer.pal(8, "Set2"))(n_clust)
ggplot(df, aes(x = Senescence_Score, y = Stemness_Score, color = Time)) +
     geom_point(size = 2.5, alpha = 0.8) +
     geom_smooth(method = "lm", se = FALSE, size = 1) +
     scale_color_manual(values = palette_colors) +
     annotate("text",
              x = min(df$Senescence_Score),
              y = max(df$Stemness_Score),
              label = paste0("R = ", R_val, ", p < ", p_val),
              hjust = 0, size = 5) +
    theme_bw() +
  theme(
    panel.grid = element_blank(),
    panel.border = element_rect(
      colour = "black",
      fill = NA,
      size = 1
    )) + theme(
  text = element_text(colour = "black")
) + theme(
  axis.title = element_text(size = 14),
  axis.text  = element_text(size = 14),
  legend.title = element_text(size = 13),
  legend.text = element_text(size = 11),
  plot.title = element_text(size = 12, face = "bold")
)

## ----emtscore-28--------------------------------------------------------------
# Extract metadata
df1 <- result$A549_TGFB1[[]]
df2 <- EMscore_result$A549_TGFB1[[]]

# Add cell names for merging
df1$cell <- rownames(df1)
df2$cell <- rownames(df2)

# Merge metadata
df <- df1 %>%
  select(cell, Stemness_Score, Senescence_Score, Time) %>%
  left_join(df2 %>% select(cell,
                           Escore,
                           Mscore), by = "cell")

# Automatically generate color palette for Time groups
n_groups <- length(unique(df$Time))
palette_colors <- colorRampPalette(brewer.pal(8, "Set2"))(n_groups)

# Define combinations to plot
x_vars <- c("Stemness_Score", "Senescence_Score")
y_vars <- c("Escore", "Mscore")

# Loop through all combinations
for (x_var in x_vars) {
  for (y_var in y_vars) {

    # Pearson correlation
    cor_test <- cor.test(df[[x_var]], df[[y_var]])
    R_val <- round(cor_test$estimate, 2)
    p_val <- format(cor_test$p.value, scientific = TRUE, digits = 2)

    # Plot
    p <- ggplot(df, aes_string(x = x_var, y = y_var, color = "Time")) +
      geom_point(size = 2.5, alpha = 0.8) +
      geom_smooth(method = "lm", se = FALSE, size = 1) +
      scale_color_manual(values = palette_colors) +
      annotate("text",
               x = min(df[[x_var]]),
               y = max(df[[y_var]]),
               label = paste0("R = ", R_val, ", p < ", p_val),
               hjust = 0,
               size = 5) +
      labs(
        x = gsub("_", " ", x_var),
        y = gsub("_", " ", y_var),
        title = paste(gsub("_", " ", y_var), "vs", gsub("_", " ", x_var))
      ) +
      theme_bw() +
  theme(
    panel.grid = element_blank(),
    panel.border = element_rect(
      colour = "black",
      fill = NA,
      size = 1
    )) + theme(
  text = element_text(colour = "black")
) + theme(
  axis.title = element_text(size = 14),
  axis.text  = element_text(size = 14),
  legend.title = element_text(size = 13),
  legend.text = element_text(size = 11),
  plot.title = element_text(size = 12, face = "bold")
)

    print(p)
  }
}

## ----restore-options----------------------------------------------------------
options(old_opts)

## ----emtscore-29--------------------------------------------------------------
sessionInfo()