---
title: "EMTscore Analysis and Plots"
author: "Haimei Wen"
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
    number_sections: true
vignette: >
  %\VignetteIndexEntry{EMT Score Vignette}
  %\VignettePackage{EMTscore}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r 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%"
)
```

# Introduction
##  What is EMTscore?
Scores of EMT gene sets have been widely used to quantify the extent of the epithelial-mesenchmal transition.
These scores are particularly useful for analyzing single-cell or bulk transcriptome samples in the context of
EMT continuum and intermediate EMT states. Analyses with EMT scores are challenged by complexity arising from
divergence of EMT progression <a href="#ref1">(Groves et al. 2024)</a> and related gene sets, as well as algorithms underlying the gene set scores. In addition,
it is often difficult to visualize and interpret EMT scores and their relation to other cellular processes.
EMTscore is a package for analyzing and visualizing the expression of
multiple EMT genesets in bulk or single-cell RNA-sequencing data.

##  Key objectives.
EMTscore computes scores for E and M gene sets. It has options for 4 commonly used
E/M gene sets and allow users to use their own gene lists. It also has options for multiple algorithms for computing EMT scores.
Importantly, it leverages nonnegative PCA for identifying multiple mesenchymal (M) scores (based on leading PCs).
It computes correlations with EMT scores with other gene set scores among samples.

##  Related Bioconductor packages.
Several Bioconductor packages provide general gene-set or signature scoring that overlaps in scope with
EMTscore, including `GSVA`, `AUCell`, `singscore`, `UCell`, `escape`, `GSEABase`, and `decoupleR`. These
packages offer excellent, method-agnostic machinery for scoring arbitrary gene sets across bulk or single-cell
data, and EMTscore builds on this ecosystem rather than replacing it (indeed, it uses `AUCell` and `GSVA`
internally as two of its scoring options). EMTscore differs by being focused specifically on the
epithelial–mesenchymal transition: it bundles curated EMT gene sets, lets users compute and directly compare
scores from several algorithms (ssGSEA, GSVA, AUCell, JASMINE, SCSE, and nnPCA) within a single workflow, and
provides EMT-oriented visualizations of epithelial versus mesenchymal states. It also adds an nnPCA-based
method for resolving divergent mesenchymal (M) programs into multiple scores, which complements the
single-score output typical of the more general tools. Users interested in scoring gene sets outside the EMT
context may find the packages above more suitable.

##  Input files.
In our examples, we assume that the input data is processed RNA-seq data stored as r
data files. For bulk data, the file contains the expression matrix.
 For single-cell data, the data file contains a seurat object.
A cell/sample annotation file should also be provided.

# Installation
EMTscore can be installed from Bioconductor using `BiocManager`:

```{r install, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("EMTscore")
```

# Example walkthrough (bulk RNA-seq)

## Load cell/sample annotation
Here we provide a `cell annotation file` <a href="#refGroves2023">Groves et al. 2023</a>. The two most important columns are: one named `name`, which represents the cell name, and another named `celltype_annotation`, which indicates the cell type classification.

```{r emtscore-2}
data(cell_annotation_file)
head(cell_annotation_file)
```

## Load gene expression data
Here we use a bulk gene expression matrix from <a href="#refGroves2023">Groves et al. 2023</a>. The full matrix (15,950 genes x 120 samples) is retrieved from Zenodo and cached locally with `BiocFileCache`, so it is downloaded only once. Row names are gene symbols and column names are sample names. (A small 544-gene subset is also bundled with the package as `data(geneExp)` for the runnable examples on the manual pages.)
```{r 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)
```

## Load gene signature data
Here, we use a built-in gene set `Panchy_et_al_E_signature` and `Panchy_et_al_M_signature` <a href="#ref2">Panchy et al. 2022</a>.

There are several additional options that users can choose from, such as: `Tan_et_al_tumor_E_signature`, `Tan_et_al_tumor_M_signature`, `Tan_et_al_cell_line_E_signature`, `Tan_et_al_cell_line_M_signature` <a href="#ref3">Tan et al. 2014</a>, `MSigDB_Hallmark` <a href="#refMSigDB">Liberzon et al. 2015</a> and `GO` <a href="#refGO">Gene Ontology Consortium 2021</a>.
```{r emtscore-4}
data("Panchy_et_al_E_signature", package = "EMTscore")
data("Panchy_et_al_M_signature", package = "EMTscore")
head(Panchy_et_al_M_signature)
```

## Prepare gene sets and save as GMT

```{r 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"))
```

## Compute scores using multiple methods

Gene set scores can be computed with one or more of the following methods:<br>
`nnPCA`: nnPCA is a nonnegative PCA based method <a href='#ref4'> (Panchy et al. 2021) </a>. It is preferred because of its high efficiency and the capacity of generating muliple scores for a gene set.<br>
`AUCell`: AUCell is a method based on area under the curve for gene ranks <a href='#ref5'> (Aibar et al. 2017) </a>.<br>
`ssGSEA`: ssGSEA is a widely used method based on Kolmogorov-Smirnov test <a href='#ref6'> (Barbie et al. 2009) </a>. <br>
`GSVA`: GSVA is a method that estimates variation of gene set activity over a sample population in an unsupervised manner <a href='#ref_GSVA'> (Hänzelmann et al. 2013) </a>.<br>
`JASMINE`: JASMINE is a recently developed method based on gene ranks  <a href='#ref7'> (Noreen et al. 2022) </a>.<br>
`SCSE`: SCSE is another normalized sum-based method <a href='#ref8'> (Pont et al. 2019) </a>.<br>


1. Calculate EMT scores with single EMT gene set
```{r 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")
```

2. Calculate EMT scores with multiple EMT gene sets
```{r 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)
```

## Prepare Data for Plotting
Here we use nnPCA scores for our plot examples.
```{r emtscore-8}
data_for_plot <- data_prepare(cell_annotation_file, nnPCA_Result_multiple, merge_colname = "name")
```

## E score vs M score Plot
As a fundamental visualization, we plot the M scores against the E scores for all samples. Sample annotations are overlaid to highlight trends across different groups, enabling us to observe patterns of epithelial–mesenchymal variation.
```{r 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
```

To comprehensively compare EMT scoring across different computational methods, we generated E–M scatter plots for all methods, including nnPCA, AUCell, ssGSEA, JASMINE, and SCSE. For each method, samples are projected onto the E (epithelial) and M (mesenchymal) score space, and annotations are overlaid to reveal group-level patterns. The main title of each plot indicates the corresponding method.
```{r 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
```

## dimension 1 M score vs dimension 2 M score Plot
One useful application of nnPCA-based EMT scores is to examine the potential divergence among M genes. To explore this, we extract the data projected onto multiple principal components, allowing us to capture the main axes of variation within the mesenchymal gene set. We then visualize the samples in this reduced-dimensional space, which helps reveal patterns of heterogeneity, potential subpopulations, or outlier behaviors that may be masked in the full gene expression matrix. This approach provides a more nuanced understanding of EMT dynamics and highlights subtle differences in mesenchymal programs across samples.
```{r 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
```

## Combined Plot
We can combine both types of scatter plot to give an overall picture.
```{r 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)
```

## E/M Score Distributions
To explore the distribution of epithelial (E) and mesenchymal (M) scores across different cell types, we plotted separate histograms for each score and cell type. This allows us to compare the overall EMT profiles between cell populations, observe shifts in epithelial or mesenchymal programs, and identify cell types enriched for higher or lower E or M scores.
```{r 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
```

## Heatmap
Heatmaps provide an effective way to visualize the expression patterns of individual genes that contribute to each principal component, highlighting variation and clustering across samples.
```{r emtscore-14}
# Example heatmap
plot_heatmap_function(t(geneExp), Panchy_et_al_M_signature)
```

# Example Walkthrough: Single-Cell RNA-seq
Next, we demonstrate the analysis workflow using a separate single-cell RNA-seq dataset <a href="#refCook2020">Cook et al. 2020</a>.

## Loading the single-cell datasets
The single-cell datasets are distributed through the companion `EMTscoreData` package
and retrieved from `ExperimentHub`. Each dataset is a `SingleCellExperiment` object for
one A549 EMT-induction condition (TGFB1, EGF, or TNF).
```{r emtscore-load-sc}
eh <- ExperimentHub::ExperimentHub()
query(eh, "EMTscoreData")

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

Each condition contains a few thousand cells. To keep the build fast we
down-sample each dataset; increase `n` (or skip this step) to use more cells.
```{r 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
)
```

## Data Importing and Basic Calculations
Users can import single-cell data either as Seurat objects (`RDS` format) or as SingleCellExperiment objects (`RDA` format). Gene sets can be provided using GMT files for input. Different scoring methods can be specified via the `method` argument. Calculated EMT scores can then be plotted against pseudotime to examine dynamic changes in epithelial–mesenchymal states.
```{r 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
```


## Building Gaussian Mixture Models (GMM) in E-M space
We can build simple GMMs for all cells in the E-M space so that we can identify different EMT states.
We can also map the existing annotations/clusters to these GMM states. This is particularly helpful to
identify clusters that are extreme E/M states or intermediate EMT states.
```{r 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)
```


## Identifying Pathways Most Associated with the EMT Process
To understand which signaling pathways are most strongly associated with the epithelial–mesenchymal transition (EMT) process, we can compute pathway activity scores for each sample and correlate them with EMT scores.
```{r 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)
```
### Visualizing the Top 10 Positively Correlated Pathways

To highlight the pathways most strongly associated with the target EMT score, we extracted the top 10 positively correlated pathways from the correlation results. A horizontal bar chart was generated using `ggplot2`, which provides clear visualization of pathway names and allows straightforward comparison of correlation strengths across pathways.

This plot offers an intuitive summary of the biological processes showing the strongest positive association with the EMT score, enabling users to quickly identify pathways that may play key roles in the underlying biological context.
```{r 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))
```


### Visualizing the Top 10 Negatively Correlated Pathways

To identify pathways that show strong inverse associations with the EMT score, we also selected the top 10 pathways exhibiting the most negative correlation values.
```{r 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))
```

## Computing Additional Scores (e.g., Stemness and Senescence Scores)
In addition to EMT scores, we can compute other relevant biological scores to further characterize cellular states, such as stemness and senescence. These scores help to assess differences in cell pluripotency, proliferative potential, and aging status.

Here, we use the `compute_Signature_score` function to calculate scores for each sample based on predefined gene signatures:

### Stemness score
Calculated using the `stemsig.tsv` gene signature file, reflecting cellular pluripotency or stem-like characteristics.
```{r emtscore-20}
signature_file <- system.file("extdata", "stemsig.tsv", package = "EMTscore")
scores <- compute_Signature_score(geneExp, signature_file , score_name = "stemness_score")
head(scores)
```

### Senescence score
Calculated using the `cellular_senescence_sig.tsv` gene signature file, reflecting the cellular aging state.
```{r 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)
```
These scores can be used alongside EMT scores for downstream analyses, such as cell-type characterization, correlation analysis, or visualization, providing a more comprehensive view of cellular states.

## Exploring the Relationship Between Stemness, Senescence, and EMT Scores in Single-Cell Data
In this section, we analyze single-cell RNA-seq datasets to compute stemness, senescence and EMT scores for individual cells and investigate their relationship. We first load multiple single-cell datasets and compute the scores using predefined gene signatures:
The stemness signature (`stemsig.tsv`) was obtained from <a href='#ref_9'> (Jing et al. 2025) </a>.<br> and the cellular senescence signature (`cellular_senescence_sig.tsv`) was derived from <a href='#ref_10'> (Malta T. 2018) </a>.<br>.

### Compute stemness and senescence scores for each cell using predefined gene signatures
```{r 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[[]])
```

### Compute EMT Scores
```{r 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")
```

### E vs M scatter plot

```{r 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
```


### Mscore PC1 vs Mscore PC2 scatter plot
```{r 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
```
### Combined Plot
```{r 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)
```

### Examine Correlation Between Stemness and Senescence
```{r 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")
)
```

### Integration and Visualization of Stemness, Senescence, and EMT Scores
In this section, we integrate stemness, senescence, and EMT scores at the single-cell level to systematically explore their relationships. By merging the computed scores for each cell, we can visualize how stemness and senescence individually correlate with the epithelial (E) and mesenchymal (M) components of EMT.

Specifically, we:

  1.Merge metadata from stemness/senescence computations and EMT score calculations to obtain a single dataset containing all    relevant scores for each cell.

  2.Visualize relationships using scatter plots, where each plot shows one pairing:

      Stemness vs. E-score

      Stemness vs. M-score

      Senescence vs. E-score

      Senescence vs. M-score
```{r 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)
  }
}
```

# Session information

Finally, we restore the user's original global options that were modified in
the setup chunk, so the vignette leaves no lingering changes to the session.
```{r restore-options}
options(old_opts)
```

```{r emtscore-29}
sessionInfo()
```


# References
<h6 id='ref1'> Groves, S. M., Panchy, N., Tyson, D. R., Harris, L. A., Quaranta, V., & Hong, T. (2023). Involvement of Epithelial–Mesenchymal Transition Genes in Small Cell Lung Cancer Phenotypic Plasticity. Cancers, 15(5), 1477. </h6>

<h6 id='ref2'> Panchy, N., Watanabe, K., Takahashi, M., Willems, A., & Hong, T. (2022). Comparative single-cell transcriptomes of dose and time dependent epithelial–mesenchymal spectrums. NAR Genom. Bioinform. 4, lqac072.  </h6>

<h6 id="refGroves2023">
Groves SM, Ildefonso GV, McAtee CO, Ozawa PMM, Ireland AS, Stauffer PE, Wasdin PT, Huang X, Qiao Y, Lim JS, Bader J, Liu Q, Simmons AJ, Lau KS, Iams WT, Hardin DP, Saff EB, Holmes WR, Tyson DR, Lovly CM, Quaranta V.
Archetype tasks link intratumoral heterogeneity to plasticity and cancer hallmarks in small cell lung cancer.
<i>Cell</i>. 2023;186(20):4306–4328.e18.
</h6>


<h6 id='ref3'> Tan, T. Z., Miow, Q. H., Miki, Y., Noda, T., Mori, S., Huang, R. Y., & Thiery, J. P. (2014).
Epithelial–mesenchymal transition spectrum quantification and its efficacy in deciphering survival and drug responses of cancer patients. EMBO Molecular Medicine, 6(10), 1279–1293. </h6>

<h6 id="refMSigDB"> Liberzon A, Birger C, Thorvaldsdóttir H, et al. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Systems. 2015;1(6):417–425. </h6>

<h6 id="refGO"> The Gene Ontology Consortium. The Gene Ontology resource: enriching a gold mine. Nucleic Acids Research. 2021;49(D1):D325–D334. </h6>


<h6 id='ref4'> Panchy, N., Watanabe, K., & Hong, T. (2021). Interpretable, scalable, and transferrable functional projection of large-scale transcriptome data using constrained matrix decomposition. Frontiers in Genetics, 12, 719099. </h6>

<h6 id='ref5'> Aibar, S., González-Blas, C. B., Moerman, T., Huynh-Thu, V. A., Imrichova, H., Hulselmans, G., ... & Aerts, S. (2017). SCENIC: single-cell regulatory network inference and clustering. Nature methods, 14(11), 1083-1086.</h6>

<h6 id='ref_GSVA'> Hänzelmann, S., Castelo, R., & Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics, 14, 7.</h6>

<h6 id='ref6'> Barbie D.A., Tamayo P., Boehm J.S., Kim S.Y., Moody S.E., Dunn I.F., Schinzel A.C., Sandy P., Meylan E., Scholl C.et al.. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. 2009; 462:108–112.  </h6>

<h6 id='ref7'> Noureen N., Ye Z., Chen Y., Wang X., Zheng S.. Signature-scoring methods developed for bulk samples are not adequate for cancer single-cell RNA sequencing data. Elife. 2022; 11:e71994. </h6>

<h6 id='ref8'> Pont F., Tosolini M., Fournié J.J.. Single-Cell Signature Explorer for comprehensive visualization of single cell signatures across scRNA-seq datasets. Nucleic Acids Res. 2019; 47:e133. </h6>

<h6 id="refCook2020"> Cook DP, Vanderhyden BC. Context specificity of the EMT transcriptional response. Nat Commun. 2020 May 1;11(1):2142. </h6>

<h6 id="ref9"> Wang, J., Zhou, X., Yu, P. et al. A transcriptome-based human universal senescence index (hUSI) robustly predicts cellular senescence under various conditions. Nat Aging. 2025 5, 1159–1175. </h6>

<h6 id="ref10"> Malta T, Sokolov A, Gentles A. Machine Learning Identifies Stemness Features Associated with Oncogenic Dedifferentiation. Cell. 2018 173, 338-354.e15. </h6>