EMTscore Analysis and Plots
2026-08-09
1 Introduction
1.1 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 (Groves et al. 2024) 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.
1.2 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.
1.4 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.
2 Installation
EMTscore can be installed from Bioconductor using BiocManager:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("EMTscore")3 Example walkthrough (bulk RNA-seq)
3.1 Load cell/sample annotation
Here we provide a cell annotation file Groves et al. 2023. 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.
data(cell_annotation_file)
head(cell_annotation_file)
#> name source celltype_annotation
#> 1 c.COLO668 CCLE A2
#> 2 c.CORL24 CCLE A
#> 3 c.CORL279 CCLE A
#> 4 c.CORL311 CCLE P
#> 5 c.CORL47 CCLE A
#> 6 c.CORL88 CCLE A23.2 Load gene expression data
Here we use a bulk gene expression matrix from Groves et al. 2023. 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.)
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)
#> m.DMS153 m.NCIH60 m.NCIH69 m.NCIH82 m.NCIH128
#> A1BG 2.95519898 3.239727011 3.599770170 3.055231198 4.579506618
#> A1CF 0.07721193 0.167166221 0.117398456 0.536178723 0.088160029
#> A2M 2.26750213 0.122302456 0.085891273 1.586479424 0.064499801
#> A2ML1 0.01101060 -0.007241361 -0.005085505 -0.004471288 -0.003818945
dim(geneExp)
#> [1] 15950 1203.3 Load gene signature data
Here, we use a built-in gene set Panchy_et_al_E_signature and Panchy_et_al_M_signature Panchy et al. 2022.
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 Tan et al. 2014, MSigDB_Hallmark Liberzon et al. 2015 and GO Gene Ontology Consortium 2021.
data("Panchy_et_al_E_signature", package = "EMTscore")
data("Panchy_et_al_M_signature", package = "EMTscore")
head(Panchy_et_al_M_signature)
#> GeneName
#> 1 AKAP12
#> 2 AKAP2
#> 3 AKT3
#> 4 ANGPTL2
#> 5 ANK2
#> 6 AP1S23.4 Prepare gene sets and save as GMT
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"))3.5 Compute scores using multiple methods
Gene set scores can be computed with one or more of the following methods:
nnPCA: nnPCA is a nonnegative PCA based method (Panchy et al. 2021) . It is preferred because of its high efficiency and the capacity of generating muliple scores for a gene set.
AUCell: AUCell is a method based on area under the curve for gene ranks (Aibar et al. 2017) .
ssGSEA: ssGSEA is a widely used method based on Kolmogorov-Smirnov test (Barbie et al. 2009) .
GSVA: GSVA is a method that estimates variation of gene set activity over a sample population in an unsupervised manner (Hänzelmann et al. 2013) .
JASMINE: JASMINE is a recently developed method based on gene ranks (Noreen et al. 2022) .
SCSE: SCSE is another normalized sum-based method (Pont et al. 2019) .
- Calculate EMT scores with single EMT gene set
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")
#> Genes in the gene sets NOT available in the dataset:
#> geneSet1: 5 (2% of 200)
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")- Calculate EMT scores with multiple EMT gene sets
gmt_file <- system.file("extdata", "EM_signature.gmt", package = "EMTscore")
nnPCA_Result_multiple <- Execute_nnPCA_parallel(geneExp, gmt_file, dimension = 1, cores = 1)
#> 121
AUCell_Result_multiple <- Execute_AUCell_parallel(geneExp, gmt_file, cores = 1)
#> 121
#> Processing pathway 1/2: Panchy_et_al_E_signature
#> Genes in the gene sets NOT available in the dataset:
#> geneSet1: 9 (4% of 232)
#> Processing pathway 2/2: Panchy_et_al_M_signature
#> Genes in the gene sets NOT available in the dataset:
#> geneSet1: 17 (9% of 193)
ssGSEA_Result_multiple <- Execute_ssGSEA_parallel(geneExp, gmt_file, cores = 1)
#> 121
GSVA_Result_multiple <- Execute_GSVA_parallel(geneExp, gmt_file, cores = 1)
#> 121
JASMINE_Result_multiple <- Execute_JASMINE_parallel(geneExp, gmt_file, cores = 1)
#> 121
SCSE_Result_multiple <- Execute_SCSE_parallel(geneExp, gmt_file, cores = 1)
#> 1213.6 Prepare Data for Plotting
Here we use nnPCA scores for our plot examples.
data_for_plot <- data_prepare(cell_annotation_file, nnPCA_Result_multiple, merge_colname = "name")3.7 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.
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")
)
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> ℹ The deprecated feature was likely used in the EMTscore package.
#> Please report the issue at <https://github.com/wenmm/EMTscore/issues>.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
#> ℹ Please use the `linewidth` argument instead.
#> ℹ The deprecated feature was likely used in the EMTscore package.
#> Please report the issue at <https://github.com/wenmm/EMTscore/issues>.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
plot1
#> `height` was translated to `width`.
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.
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)# 自动添加方法名称作为主标题
})
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
names(plot_list) <- names(method_list)
combined_plot <- Arrange_plots(plot_list, ncol_per_row = 3)
#> `height` was translated to `width`.
#> `height` was translated to `width`.
#> `height` was translated to `width`.
#> `height` was translated to `width`.
#> `height` was translated to `width`.
#> `height` was translated to `width`.
combined_plot
3.8 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.
# 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")
)
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
plot2
#> `height` was translated to `width`.
3.9 Combined Plot
We can combine both types of scatter plot to give an overall picture.
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"
)
#> `height` was translated to `width`.
#> `height` was translated to `width`.
#> `height` was translated to `width`.
print(combined_plot)
3.10 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.
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
#> `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
3.11 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.
# Example heatmap
plot_heatmap_function(t(geneExp), Panchy_et_al_M_signature)
4 Example Walkthrough: Single-Cell RNA-seq
Next, we demonstrate the analysis workflow using a separate single-cell RNA-seq dataset Cook et al. 2020.
4.1 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).
eh <- ExperimentHub::ExperimentHub()
query(eh, "EMTscoreData")
#> ExperimentHub with 12 records
#> # snapshotDate(): 2026-07-30
#> # $dataprovider: David P. Cook
#> # $species: Homo sapiens
#> # $rdataclass: SingleCellExperiment
#> # additional mcols(): taxonomyid, genome, description,
#> # coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
#> # rdatapath, sourceurl, sourcetype
#> # retrieve records with, e.g., 'object[["EH10282"]]'
#>
#> title
#> EH10282 | MCF7_TNF.rda
#> EH10283 | MCF7_EGF.rda
#> EH10284 | MCF7_TGFB1.rda
#> EH10285 | OVCA420_TNF.rda
#> EH10286 | OVCA420_EGF.rda
#> ... ...
#> EH10289 | DU145_EGF.rda
#> EH10290 | DU145_TGFB1.rda
#> EH10291 | A549_TNF.rda
#> EH10292 | A549_EGF.rda
#> EH10293 | A549_TGFB1.rda
A549_TNF <- eh[["EH10291"]]
#> see ?EMTscoreData and browseVignettes('EMTscoreData') for documentation
#> loading from cache
A549_EGF <- eh[["EH10292"]]
#> see ?EMTscoreData and browseVignettes('EMTscoreData') for documentation
#> loading from cache
A549_TGFB1 <- eh[["EH10293"]]
#> see ?EMTscoreData and browseVignettes('EMTscoreData') for documentation
#> loading from cacheEach 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.
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
)4.2 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.
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")
#> Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
#> ℹ Please use tidy evaluation idioms with `aes()`.
#> ℹ See also `vignette("ggplot2-in-packages")` for more information.
#> ℹ The deprecated feature was likely used in the EMTscore package.
#> Please report the issue at <https://github.com/wenmm/EMTscore/issues>.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
p_nnPCA
4.3 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.
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)
#> 121
#> 121
#> 121
#> 121
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)
4.4 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.
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)4.4.1 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.
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))
#> Warning in geom_richtext(aes(label = Pathway, y = 0), hjust = 0, fill = NA, :
#> Ignoring unknown parameters: `background_colour`
4.4.2 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.
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))
#> Warning in geom_richtext(aes(label = Pathway, y = 0), hjust = 1, fill = NA, :
#> Ignoring unknown parameters: `background_colour`
4.5 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:
4.5.1 Stemness score
Calculated using the stemsig.tsv gene signature file, reflecting cellular pluripotency or stem-like characteristics.
signature_file <- system.file("extdata", "stemsig.tsv", package = "EMTscore")
scores <- compute_Signature_score(geneExp, signature_file , score_name = "stemness_score")
head(scores)
#> stemness_score
#> m.DMS153 0.5344621
#> m.NCIH60 0.1719915
#> m.NCIH69 0.3539699
#> m.NCIH82 0.4988987
#> m.NCIH128 0.4462295
#> m.NCIH146 0.34359804.5.2 Senescence score
Calculated using the cellular_senescence_sig.tsv gene signature file, reflecting the cellular aging state.
signature_file <- system.file("extdata", "cellular_senescence_sig.tsv", package = "EMTscore")
scores <- compute_Signature_score(geneExp, signature_file, score_name = "senescence_score")
head(scores)
#> senescence_score
#> m.DMS153 0.3504790
#> m.NCIH60 0.1974465
#> m.NCIH69 0.1508204
#> m.NCIH82 0.2459747
#> m.NCIH128 0.2579179
#> m.NCIH146 0.3323265These 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.
4.6 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 (Jing et al. 2025) .
and the cellular senescence signature (cellular_senescence_sig.tsv) was derived from (Malta T. 2018) .
.
4.6.1 Compute stemness and senescence scores for each cell using predefined gene signatures
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"))
#> Converting SCE to Seurat: A549_TGFB1
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper structure
#> Updating matrix keys for DimReduc 'PCA'
#> Updating matrix keys for DimReduc 'UMAP'
#> Updating matrix keys for DimReduc 'UMAP_PSEUDO'
#> Ensuring keys are in the proper structure
#> Ensuring feature names don't have underscores or pipes
#> Updating slots in RNA
#> Updating slots in PCA
#> Updating slots in UMAP
#> Setting UMAP DimReduc to global
#> Updating slots in UMAP_PSEUDO
#> Setting UMAP_PSEUDO DimReduc to global
#> Validating object structure for Assay 'RNA'
#> Validating object structure for DimReduc 'PCA'
#> Validating object structure for DimReduc 'UMAP'
#> Validating object structure for DimReduc 'UMAP_PSEUDO'
#> Object representation is consistent with the most current Seurat version
#> Converting SCE to Seurat: A549_EGF
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper structure
#> Updating matrix keys for DimReduc 'PCA'
#> Updating matrix keys for DimReduc 'UMAP'
#> Updating matrix keys for DimReduc 'UMAP_PSEUDO'
#> Ensuring keys are in the proper structure
#> Ensuring feature names don't have underscores or pipes
#> Updating slots in RNA
#> Updating slots in PCA
#> Updating slots in UMAP
#> Setting UMAP DimReduc to global
#> Updating slots in UMAP_PSEUDO
#> Setting UMAP_PSEUDO DimReduc to global
#> Validating object structure for Assay 'RNA'
#> Validating object structure for DimReduc 'PCA'
#> Validating object structure for DimReduc 'UMAP'
#> Validating object structure for DimReduc 'UMAP_PSEUDO'
#> Object representation is consistent with the most current Seurat version
#> Converting SCE to Seurat: A549_TNF
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper structure
#> Updating matrix keys for DimReduc 'PCA'
#> Updating matrix keys for DimReduc 'UMAP'
#> Updating matrix keys for DimReduc 'UMAP_PSEUDO'
#> Ensuring keys are in the proper structure
#> Ensuring feature names don't have underscores or pipes
#> Updating slots in RNA
#> Updating slots in PCA
#> Updating slots in UMAP
#> Setting UMAP DimReduc to global
#> Updating slots in UMAP_PSEUDO
#> Setting UMAP_PSEUDO DimReduc to global
#> Validating object structure for Assay 'RNA'
#> Validating object structure for DimReduc 'PCA'
#> Validating object structure for DimReduc 'UMAP'
#> Validating object structure for DimReduc 'UMAP_PSEUDO'
#> Object representation is consistent with the most current Seurat version
head(result$A549_TGFB1[[]])
#> orig.ident nCount_RNA nFeature_RNA percent.mito
#> Mix4b_GCTTCCACAGTAAGAT Mix4b 12304 2973 0.04425497
#> Mix4a_CGAGCACTCGTCTGCT Mix4a 15755 3305 0.03817617
#> Mix2_CGGCTAGTCAGTTTGG Mix2 9789 2560 0.06295276
#> Mix4b_CCACGGATCGAGCCCA Mix4b 12267 2835 0.04853420
#> Mix3b_TCTTCGGCAATAAGCA Mix3b 17831 3868 0.04494508
#> Mix3b_CATTCGCAGCCAGGAT Mix3b 19465 3817 0.05425243
#> Sample CellLine Treatment Time Doublet
#> Mix4b_GCTTCCACAGTAAGAT A549_TGFB1_1d A549 TGFB1 1d Singlet
#> Mix4a_CGAGCACTCGTCTGCT A549_TGFB1_3d A549 TGFB1 3d Singlet
#> Mix2_CGGCTAGTCAGTTTGG A549_TGFB1_1d A549 TGFB1 1d Singlet
#> Mix4b_CCACGGATCGAGCCCA A549_TGFB1_3d A549 TGFB1 3d Singlet
#> Mix3b_TCTTCGGCAATAAGCA A549_TGFB1_8h A549 TGFB1 8h Singlet
#> Mix3b_CATTCGCAGCCAGGAT A549_TGFB1_7d A549 TGFB1 7d Singlet
#> S.Score G2M.Score Phase Mix SCT_snn_res.0.8
#> Mix4b_GCTTCCACAGTAAGAT -0.03511477 -1.0326206 G1 Mix4b 9
#> Mix4a_CGAGCACTCGTCTGCT -0.31431997 -1.9355842 G1 Mix4a 9
#> Mix2_CGGCTAGTCAGTTTGG 0.12863995 -0.9834930 S Mix2 9
#> Mix4b_CCACGGATCGAGCCCA 0.01644399 -1.2497361 S Mix4b 9
#> Mix3b_TCTTCGGCAATAAGCA 0.40030832 -0.7211847 S Mix3b 3
#> Mix3b_CATTCGCAGCCAGGAT 0.25042823 -2.3383846 S Mix3b 9
#> SCT_snn_res.0.1 RNA_snn_res.0.1 seurat_clusters
#> Mix4b_GCTTCCACAGTAAGAT 0 0 2
#> Mix4a_CGAGCACTCGTCTGCT 0 0 0
#> Mix2_CGGCTAGTCAGTTTGG 0 0 2
#> Mix4b_CCACGGATCGAGCCCA 0 0 0
#> Mix3b_TCTTCGGCAATAAGCA 0 0 1
#> Mix3b_CATTCGCAGCCAGGAT 0 0 0
#> RNA_snn_res.0.05 Cluster RNA_snn_res.0.5 Pseudotime
#> Mix4b_GCTTCCACAGTAAGAT 0 0 2 0.5356527
#> Mix4a_CGAGCACTCGTCTGCT 0 0 0 0.6486413
#> Mix2_CGGCTAGTCAGTTTGG 0 0 2 0.4916917
#> Mix4b_CCACGGATCGAGCCCA 0 0 0 0.5998631
#> Mix3b_TCTTCGGCAATAAGCA 0 0 1 0.2758947
#> Mix3b_CATTCGCAGCCAGGAT 0 0 0 0.6584807
#> ident Stemness_Score Senescence_Score
#> Mix4b_GCTTCCACAGTAAGAT 2 0.5546042 0.3371816
#> Mix4a_CGAGCACTCGTCTGCT 0 0.4168301 0.6174615
#> Mix2_CGGCTAGTCAGTTTGG 2 0.2696344 0.4621716
#> Mix4b_CCACGGATCGAGCCCA 0 0.4572547 0.3467303
#> Mix3b_TCTTCGGCAATAAGCA 1 0.6890547 0.3302015
#> Mix3b_CATTCGCAGCCAGGAT 0 0.4079503 0.42601924.6.2 Compute EMT Scores
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)
#> 121
#> Converting SCE to Seurat: A549_TGFB1
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper structure
#> Updating matrix keys for DimReduc 'PCA'
#> Updating matrix keys for DimReduc 'UMAP'
#> Updating matrix keys for DimReduc 'UMAP_PSEUDO'
#> Ensuring keys are in the proper structure
#> Ensuring feature names don't have underscores or pipes
#> Updating slots in RNA
#> Updating slots in PCA
#> Updating slots in UMAP
#> Setting UMAP DimReduc to global
#> Updating slots in UMAP_PSEUDO
#> Setting UMAP_PSEUDO DimReduc to global
#> Validating object structure for Assay 'RNA'
#> Validating object structure for DimReduc 'PCA'
#> Validating object structure for DimReduc 'UMAP'
#> Validating object structure for DimReduc 'UMAP_PSEUDO'
#> Object representation is consistent with the most current Seurat version
#> 121
#> Converting SCE to Seurat: A549_EGF
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper structure
#> Updating matrix keys for DimReduc 'PCA'
#> Updating matrix keys for DimReduc 'UMAP'
#> Updating matrix keys for DimReduc 'UMAP_PSEUDO'
#> Ensuring keys are in the proper structure
#> Ensuring feature names don't have underscores or pipes
#> Updating slots in RNA
#> Updating slots in PCA
#> Updating slots in UMAP
#> Setting UMAP DimReduc to global
#> Updating slots in UMAP_PSEUDO
#> Setting UMAP_PSEUDO DimReduc to global
#> Validating object structure for Assay 'RNA'
#> Validating object structure for DimReduc 'PCA'
#> Validating object structure for DimReduc 'UMAP'
#> Validating object structure for DimReduc 'UMAP_PSEUDO'
#> Object representation is consistent with the most current Seurat version
#> 121
#> Converting SCE to Seurat: A549_TNF
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper structure
#> Updating matrix keys for DimReduc 'PCA'
#> Updating matrix keys for DimReduc 'UMAP'
#> Updating matrix keys for DimReduc 'UMAP_PSEUDO'
#> Ensuring keys are in the proper structure
#> Ensuring feature names don't have underscores or pipes
#> Updating slots in RNA
#> Updating slots in PCA
#> Updating slots in UMAP
#> Setting UMAP DimReduc to global
#> Updating slots in UMAP_PSEUDO
#> Setting UMAP_PSEUDO DimReduc to global
#> Validating object structure for Assay 'RNA'
#> Validating object structure for DimReduc 'PCA'
#> Validating object structure for DimReduc 'UMAP'
#> Validating object structure for DimReduc 'UMAP_PSEUDO'
#> Object representation is consistent with the most current Seurat version
#> 121
head(EMscore_result$A549_TGFB1[[]])
#> orig.ident nCount_RNA nFeature_RNA percent.mito
#> Mix4b_GCTTCCACAGTAAGAT Mix4b 12304 2973 0.04425497
#> Mix4a_CGAGCACTCGTCTGCT Mix4a 15755 3305 0.03817617
#> Mix2_CGGCTAGTCAGTTTGG Mix2 9789 2560 0.06295276
#> Mix4b_CCACGGATCGAGCCCA Mix4b 12267 2835 0.04853420
#> Mix3b_TCTTCGGCAATAAGCA Mix3b 17831 3868 0.04494508
#> Mix3b_CATTCGCAGCCAGGAT Mix3b 19465 3817 0.05425243
#> Sample CellLine Treatment Time Doublet
#> Mix4b_GCTTCCACAGTAAGAT A549_TGFB1_1d A549 TGFB1 1d Singlet
#> Mix4a_CGAGCACTCGTCTGCT A549_TGFB1_3d A549 TGFB1 3d Singlet
#> Mix2_CGGCTAGTCAGTTTGG A549_TGFB1_1d A549 TGFB1 1d Singlet
#> Mix4b_CCACGGATCGAGCCCA A549_TGFB1_3d A549 TGFB1 3d Singlet
#> Mix3b_TCTTCGGCAATAAGCA A549_TGFB1_8h A549 TGFB1 8h Singlet
#> Mix3b_CATTCGCAGCCAGGAT A549_TGFB1_7d A549 TGFB1 7d Singlet
#> S.Score G2M.Score Phase Mix SCT_snn_res.0.8
#> Mix4b_GCTTCCACAGTAAGAT -0.03511477 -1.0326206 G1 Mix4b 9
#> Mix4a_CGAGCACTCGTCTGCT -0.31431997 -1.9355842 G1 Mix4a 9
#> Mix2_CGGCTAGTCAGTTTGG 0.12863995 -0.9834930 S Mix2 9
#> Mix4b_CCACGGATCGAGCCCA 0.01644399 -1.2497361 S Mix4b 9
#> Mix3b_TCTTCGGCAATAAGCA 0.40030832 -0.7211847 S Mix3b 3
#> Mix3b_CATTCGCAGCCAGGAT 0.25042823 -2.3383846 S Mix3b 9
#> SCT_snn_res.0.1 RNA_snn_res.0.1 seurat_clusters
#> Mix4b_GCTTCCACAGTAAGAT 0 0 2
#> Mix4a_CGAGCACTCGTCTGCT 0 0 0
#> Mix2_CGGCTAGTCAGTTTGG 0 0 2
#> Mix4b_CCACGGATCGAGCCCA 0 0 0
#> Mix3b_TCTTCGGCAATAAGCA 0 0 1
#> Mix3b_CATTCGCAGCCAGGAT 0 0 0
#> RNA_snn_res.0.05 Cluster RNA_snn_res.0.5 Pseudotime
#> Mix4b_GCTTCCACAGTAAGAT 0 0 2 0.5356527
#> Mix4a_CGAGCACTCGTCTGCT 0 0 0 0.6486413
#> Mix2_CGGCTAGTCAGTTTGG 0 0 2 0.4916917
#> Mix4b_CCACGGATCGAGCCCA 0 0 0 0.5998631
#> Mix3b_TCTTCGGCAATAAGCA 0 0 1 0.2758947
#> Mix3b_CATTCGCAGCCAGGAT 0 0 0 0.6584807
#> ident Escore Mscore
#> Mix4b_GCTTCCACAGTAAGAT 2 -1.252093 0.74292141
#> Mix4a_CGAGCACTCGTCTGCT 0 -2.090583 1.41903314
#> Mix2_CGGCTAGTCAGTTTGG 2 -1.026221 -0.47984540
#> Mix4b_CCACGGATCGAGCCCA 0 -1.088944 2.38074330
#> Mix3b_TCTTCGGCAATAAGCA 1 2.162446 -1.03455609
#> Mix3b_CATTCGCAGCCAGGAT 0 -1.625282 -0.05151672
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)
#> Converting SCE to Seurat: A549_TGFB1
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper structure
#> Updating matrix keys for DimReduc 'PCA'
#> Updating matrix keys for DimReduc 'UMAP'
#> Updating matrix keys for DimReduc 'UMAP_PSEUDO'
#> Ensuring keys are in the proper structure
#> Ensuring feature names don't have underscores or pipes
#> Updating slots in RNA
#> Updating slots in PCA
#> Updating slots in UMAP
#> Setting UMAP DimReduc to global
#> Updating slots in UMAP_PSEUDO
#> Setting UMAP_PSEUDO DimReduc to global
#> Validating object structure for Assay 'RNA'
#> Validating object structure for DimReduc 'PCA'
#> Validating object structure for DimReduc 'UMAP'
#> Validating object structure for DimReduc 'UMAP_PSEUDO'
#> Object representation is consistent with the most current Seurat version
#> Converting SCE to Seurat: A549_EGF
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper structure
#> Updating matrix keys for DimReduc 'PCA'
#> Updating matrix keys for DimReduc 'UMAP'
#> Updating matrix keys for DimReduc 'UMAP_PSEUDO'
#> Ensuring keys are in the proper structure
#> Ensuring feature names don't have underscores or pipes
#> Updating slots in RNA
#> Updating slots in PCA
#> Updating slots in UMAP
#> Setting UMAP DimReduc to global
#> Updating slots in UMAP_PSEUDO
#> Setting UMAP_PSEUDO DimReduc to global
#> Validating object structure for Assay 'RNA'
#> Validating object structure for DimReduc 'PCA'
#> Validating object structure for DimReduc 'UMAP'
#> Validating object structure for DimReduc 'UMAP_PSEUDO'
#> Object representation is consistent with the most current Seurat version
#> Converting SCE to Seurat: A549_TNF
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper structure
#> Updating matrix keys for DimReduc 'PCA'
#> Updating matrix keys for DimReduc 'UMAP'
#> Updating matrix keys for DimReduc 'UMAP_PSEUDO'
#> Ensuring keys are in the proper structure
#> Ensuring feature names don't have underscores or pipes
#> Updating slots in RNA
#> Updating slots in PCA
#> Updating slots in UMAP
#> Setting UMAP DimReduc to global
#> Updating slots in UMAP_PSEUDO
#> Setting UMAP_PSEUDO DimReduc to global
#> Validating object structure for Assay 'RNA'
#> Validating object structure for DimReduc 'PCA'
#> Validating object structure for DimReduc 'UMAP'
#> Validating object structure for DimReduc 'UMAP_PSEUDO'
#> Object representation is consistent with the most current Seurat version
head(Mscore_result$A549_TGFB1[[]])
#> orig.ident nCount_RNA nFeature_RNA percent.mito
#> Mix4b_GCTTCCACAGTAAGAT Mix4b 12304 2973 0.04425497
#> Mix4a_CGAGCACTCGTCTGCT Mix4a 15755 3305 0.03817617
#> Mix2_CGGCTAGTCAGTTTGG Mix2 9789 2560 0.06295276
#> Mix4b_CCACGGATCGAGCCCA Mix4b 12267 2835 0.04853420
#> Mix3b_TCTTCGGCAATAAGCA Mix3b 17831 3868 0.04494508
#> Mix3b_CATTCGCAGCCAGGAT Mix3b 19465 3817 0.05425243
#> Sample CellLine Treatment Time Doublet
#> Mix4b_GCTTCCACAGTAAGAT A549_TGFB1_1d A549 TGFB1 1d Singlet
#> Mix4a_CGAGCACTCGTCTGCT A549_TGFB1_3d A549 TGFB1 3d Singlet
#> Mix2_CGGCTAGTCAGTTTGG A549_TGFB1_1d A549 TGFB1 1d Singlet
#> Mix4b_CCACGGATCGAGCCCA A549_TGFB1_3d A549 TGFB1 3d Singlet
#> Mix3b_TCTTCGGCAATAAGCA A549_TGFB1_8h A549 TGFB1 8h Singlet
#> Mix3b_CATTCGCAGCCAGGAT A549_TGFB1_7d A549 TGFB1 7d Singlet
#> S.Score G2M.Score Phase Mix SCT_snn_res.0.8
#> Mix4b_GCTTCCACAGTAAGAT -0.03511477 -1.0326206 G1 Mix4b 9
#> Mix4a_CGAGCACTCGTCTGCT -0.31431997 -1.9355842 G1 Mix4a 9
#> Mix2_CGGCTAGTCAGTTTGG 0.12863995 -0.9834930 S Mix2 9
#> Mix4b_CCACGGATCGAGCCCA 0.01644399 -1.2497361 S Mix4b 9
#> Mix3b_TCTTCGGCAATAAGCA 0.40030832 -0.7211847 S Mix3b 3
#> Mix3b_CATTCGCAGCCAGGAT 0.25042823 -2.3383846 S Mix3b 9
#> SCT_snn_res.0.1 RNA_snn_res.0.1 seurat_clusters
#> Mix4b_GCTTCCACAGTAAGAT 0 0 2
#> Mix4a_CGAGCACTCGTCTGCT 0 0 0
#> Mix2_CGGCTAGTCAGTTTGG 0 0 2
#> Mix4b_CCACGGATCGAGCCCA 0 0 0
#> Mix3b_TCTTCGGCAATAAGCA 0 0 1
#> Mix3b_CATTCGCAGCCAGGAT 0 0 0
#> RNA_snn_res.0.05 Cluster RNA_snn_res.0.5 Pseudotime
#> Mix4b_GCTTCCACAGTAAGAT 0 0 2 0.5356527
#> Mix4a_CGAGCACTCGTCTGCT 0 0 0 0.6486413
#> Mix2_CGGCTAGTCAGTTTGG 0 0 2 0.4916917
#> Mix4b_CCACGGATCGAGCCCA 0 0 0 0.5998631
#> Mix3b_TCTTCGGCAATAAGCA 0 0 1 0.2758947
#> Mix3b_CATTCGCAGCCAGGAT 0 0 0 0.6584807
#> ident Mscore_PC1 Mscore_PC2
#> Mix4b_GCTTCCACAGTAAGAT 2 1.3830702 1.23193708
#> Mix4a_CGAGCACTCGTCTGCT 0 2.6292531 -0.02719862
#> Mix2_CGGCTAGTCAGTTTGG 2 0.3742408 0.62699570
#> Mix4b_CCACGGATCGAGCCCA 0 1.9583972 0.56242490
#> Mix3b_TCTTCGGCAATAAGCA 1 -1.9880096 -0.14310837
#> Mix3b_CATTCGCAGCCAGGAT 0 1.3065442 0.30206067
plot_EMT_from_objects(Mscore_result, col_name = "Pseudotime", emt_score_col = "Mscore_PC1")
#> `geom_smooth()` using formula = 'y ~ x'
plot_EMT_from_objects(Mscore_result, col_name = "Pseudotime", emt_score_col = "Mscore_PC2")
#> `geom_smooth()` using formula = 'y ~ x'
4.6.3 E vs M scatter plot
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")
)
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
plot1
#> `height` was translated to `width`.
4.6.4 Mscore PC1 vs Mscore PC2 scatter plot
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")
)
#> Scale for colour is already present.
#> Adding another scale for colour, which will replace the existing scale.
plot2
#> `height` was translated to `width`.
### Combined Plot
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"
)
#> `height` was translated to `width`.
#> `height` was translated to `width`.
#> `height` was translated to `width`.
print(combined_plot)
4.6.5 Examine Correlation Between Stemness and Senescence
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")
)
#> `geom_smooth()` using formula = 'y ~ x'
4.6.6 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# 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)
}
}
#> `geom_smooth()` using formula = 'y ~ x'
#> `geom_smooth()` using formula = 'y ~ x'
#> `geom_smooth()` using formula = 'y ~ x'
#> `geom_smooth()` using formula = 'y ~ x'
5 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.
options(old_opts)sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 grid stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] Seurat_5.5.1 SeuratObject_5.4.0
#> [3] sp_2.2-3 SingleCellExperiment_1.35.2
#> [5] SummarizedExperiment_1.43.0 Biobase_2.73.2
#> [7] GenomicRanges_1.65.1 Seqinfo_1.3.0
#> [9] IRanges_2.47.2 S4Vectors_0.51.6
#> [11] MatrixGenerics_1.25.0 matrixStats_1.5.0
#> [13] ExperimentHub_3.3.1 AnnotationHub_4.3.2
#> [15] BiocFileCache_3.3.0 dbplyr_2.6.0
#> [17] BiocGenerics_0.59.11 generics_0.1.4
#> [19] EMTscoreData_1.1.0 RColorBrewer_1.1-3
#> [21] ggalluvial_0.12.6 mclust_6.1.3
#> [23] GSA_1.03.3 AUCell_1.35.0
#> [25] EMTscore_0.99.10 ComplexHeatmap_2.29.0
#> [27] ggthemes_5.2.0 paletteer_1.7.0
#> [29] circlize_0.4.18 pheatmap_1.0.13
#> [31] ggpubr_1.0.0 gridExtra_2.3.1
#> [33] dplyr_1.2.1 ggtext_0.1.2
#> [35] ggplot2_4.0.3 nsprcomp_0.5.1-2
#> [37] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] GSVA_2.7.12 spatstat.sparse_3.2-0
#> [3] httr_1.4.8 doParallel_1.0.17
#> [5] tools_4.6.1 sctransform_0.4.3
#> [7] backports_1.5.1 R6_2.6.1
#> [9] HDF5Array_1.41.1 mgcv_1.9-4
#> [11] uwot_0.2.4 rhdf5filters_1.25.4
#> [13] GetoptLong_1.1.1 litedown_0.10
#> [15] withr_3.0.3 progressr_1.0.0
#> [17] cli_3.6.6 Cairo_1.7-0
#> [19] spatstat.explore_3.8-2 fastDummies_1.7.6
#> [21] isoband_0.3.0 labeling_0.4.3
#> [23] sass_0.4.10 S7_0.2.2
#> [25] spatstat.data_3.1-9 ggridges_0.5.7
#> [27] pbapply_1.7-4 commonmark_2.0.0
#> [29] R.utils_2.13.0 dichromat_2.0-1
#> [31] parallelly_1.48.0 RSQLite_3.53.3
#> [33] shape_1.4.6.1 spatstat.random_3.5-1
#> [35] ica_1.0-3 car_3.1-5
#> [37] Matrix_1.7-6 abind_1.4-8
#> [39] R.methodsS3_1.8.2 lifecycle_1.0.5
#> [41] yaml_2.3.12 carData_3.0-6
#> [43] rhdf5_2.57.9 SparseArray_1.13.2
#> [45] Rtsne_0.17 blob_1.3.0
#> [47] promises_1.5.0 crayon_1.5.3
#> [49] miniUI_0.1.2 lattice_0.22-9
#> [51] beachmat_2.29.0 cowplot_1.2.0
#> [53] annotate_1.91.0 KEGGREST_1.53.6
#> [55] magick_2.9.1 pillar_1.11.1
#> [57] knitr_1.51 rjson_0.2.23
#> [59] future.apply_1.20.2 codetools_0.2-20
#> [61] glue_1.8.1 spatstat.univar_3.2-0
#> [63] data.table_1.18.4 memuse_4.2-3
#> [65] vctrs_0.7.3 png_0.1-9
#> [67] spam_2.11-4 gtable_0.3.6
#> [69] rematch2_2.1.2 cachem_1.1.0
#> [71] xfun_0.60 S4Arrays_1.13.0
#> [73] mime_0.13 survival_3.8-9
#> [75] iterators_1.0.14 tinytex_0.60
#> [77] fitdistrplus_1.2-6 ROCR_1.0-12
#> [79] nlme_3.1-170 bit64_4.8.2
#> [81] filelock_1.0.3 RcppAnnoy_0.0.23
#> [83] bslib_0.12.0 irlba_2.3.7
#> [85] KernSmooth_2.23-26 otel_0.2.0
#> [87] colorspace_2.1-3 DBI_1.3.0
#> [89] tidyselect_1.2.1 curl_7.1.0
#> [91] bit_4.6.0 compiler_4.6.1
#> [93] httr2_1.3.0 graph_1.91.0
#> [95] h5mread_1.5.0 xml2_1.6.0
#> [97] DelayedArray_0.39.4 plotly_4.12.1
#> [99] bookdown_0.47 scales_1.4.0
#> [101] lmtest_0.9-40 rappdirs_0.3.4
#> [103] goftest_1.2-3 stringr_1.6.0
#> [105] SpatialExperiment_1.23.0 digest_0.6.39
#> [107] spatstat.utils_3.2-4 rmarkdown_2.31
#> [109] XVector_0.53.0 htmltools_0.5.9
#> [111] pkgconfig_2.0.3 sparseMatrixStats_1.25.0
#> [113] fastmap_1.2.0 rlang_1.3.0
#> [115] GlobalOptions_0.1.4 htmlwidgets_1.6.4
#> [117] shiny_1.14.0 DelayedMatrixStats_1.35.0
#> [119] farver_2.1.2 jquerylib_0.1.4
#> [121] zoo_1.9-0 jsonlite_2.0.0
#> [123] BiocParallel_1.47.0 R.oo_1.27.1
#> [125] BiocSingular_1.29.0 magrittr_2.0.5
#> [127] Formula_1.2-6 dotCall64_1.2
#> [129] patchwork_1.3.2 Rhdf5lib_2.1.0
#> [131] Rcpp_1.1.2 reticulate_1.46.0
#> [133] stringi_1.8.9 MASS_7.3-66
#> [135] plyr_1.8.9 parallel_4.6.1
#> [137] listenv_1.0.0 ggrepel_0.9.8
#> [139] deldir_2.0-4 Biostrings_2.81.6
#> [141] splines_4.6.1 tensor_1.5.1
#> [143] gridtext_0.1.6 igraph_2.3.3
#> [145] spatstat.geom_3.8-2 markdown_2.0
#> [147] ggsignif_0.6.4 RcppHNSW_0.7.0
#> [149] reshape2_1.4.5 ScaledMatrix_1.21.0
#> [151] BiocVersion_3.24.0 XML_3.99-0.23
#> [153] evaluate_1.0.5 BiocManager_1.30.27
#> [155] foreach_1.5.2 httpuv_1.6.17
#> [157] polyclip_1.10-7 RANN_2.6.2
#> [159] tidyr_1.3.2 purrr_1.2.2
#> [161] future_1.75.0 clue_0.3-68
#> [163] scattermore_1.2 BiocBaseUtils_1.15.1
#> [165] rsvd_1.0.5 broom_1.0.13
#> [167] xtable_1.8-8 RSpectra_0.16-2
#> [169] rstatix_1.1.0 later_1.4.8
#> [171] viridisLite_0.4.3 tibble_3.3.1
#> [173] memoise_2.0.1 AnnotationDbi_1.75.2
#> [175] cluster_2.1.8.3 globals_0.19.1
#> [177] GSEABase_1.75.0