Robust and efficient annotation of cell states through gene signature scoring

ANS method

We propose an improvement to the widely used gene signature scoring method originally introduced by Tirosh et al. (2016), which had been implemented in Seurat (Satija et al. 2015) and, with minor modifications, in SCANPY (Wolf et al. 2018). In brief, the Tirosh method computes gene signature scores for each cell by subtracting a control expression from the signature genes’ expressions. The major drawback of this approach comes from the control gene selection, which uses fixed binning of genes according to their average gene expression across cells (default: 25 bins) (Supplemental Fig. S1). As verified by our benchmarking, such selection may bias the estimates of expression gain, specifically for highly expressed genes, which are frequently used in gene signatures. Our proposed method, ANS, solves this issue by creating control genes for each signature gene separately, by selecting genes in an average-expression neighborhood closest to the signature genes’ average expression (Methods).

Design of the benchmark

In our work, we compared eight gene signature scoring methods (Table 1), including ANS, the original Tirosh method (as implemented in Seurat [Satija et al. 2015] and SCANPY [Wolf et al. 2018]), its variants using all genes (Seurat_AG) or least variable genes (Seurat_LVG) as controls, JASMINE (Jasmine_LH and Jasmine_OR settings) (Noureen et al. 2022), and UCell (Andreatta and Carmona 2021). Like Seurat and SCANPY, Seurat_AG and Seurat_LVG first compute the average expression of the genes, sort them in increasing order, and divide them into 25 equally sized expression bins. Although Seurat_AG selects all genes from an expression bin as control genes, Seurat_LVG chooses genes with the smallest dispersion. For a unified, Python-based benchmark, all methods developed in R (Seurat, JASMINE, UCell) were reimplemented in Python (Fig. 1A). We confirmed the consistency between the R and Python versions of ANS (Methods) (Fig. 1B; Supplemental Table S1), ensuring reproducible and reliable benchmarking.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

Performance and robustness analysis of gene signature scoring methods. (A) Benchmark questions for the gene signature scoring methods and availability of the methods in Python (blue) and R (red). Reimplementations of the original R packages in Python are indicated with a white asterisk. (B) The scores obtained from R and Python implementations compared for two cancer data sets (CRC, ESCC). Each subplot shows the Pearson correlation coefficient (R) between implementations. Discrepancies in Seurat scores between R (“AddModuleScore” of Seurat) and Python are attributed to randomization in the scoring method. (C) Control gene selection bias of the scoring methods for the top 8% of highly expressed genes in B memory kappa cells within the PBMC data set. The x-axis shows genes sorted by their average expression levels; the y-axis shows the mean score across all cells for a single-gene signature. Vertical dashed lines indicate the expression bin boundary and the top 50 highly expressed genes. The bias of a scoring method is indicated by how far the mean score of a gene deviates from zero. (D) The influence of data set composition and batch effect on scoring CRC cells using a 100-gene signature associated with malignant cells. Each dot represents the mean score for all cells within a sample, grouped by cell type (malignant in blue or nonmalignant in orange), sequencing chemistry type (SC3Pv2 in gray or SC3Pv3 in beige), and scoring mode (scoring all the samples together or individually). The black horizontal bar represents the mean value of all dots within each group. P-value annotation: (ns) P > 0.05, (*) 0.01 < P ≤ 0.05, (**) 1 × 10⁻³ < p ≤ 0.01, (***) 1 × 10⁻⁴ < p ≤ 1 × 10⁻³, (****) p ≤ 1 × 10⁻⁴.

We benchmarked the methods using varied scRNA-seq data sets, including peripheral blood mononuclear cells (PBMCs) (Hao et al. 2021), colorectal carcinoma (CRC) (Pelka et al. 2021) esophageal squamous cell carcinoma (ESCC) (Zhang et al. 2021a), breast carcinoma (BRCA) (Wu et al. 2021), high-grade serous ovarian cancer (HGSOC) (Vázquez-García et al. 2022), cutaneous squamous cell carcinoma (cSCC) (Ji et al. 2020), and lung adenocarcinoma (LUAD) (Kim et al. 2020), as well as a neuronal differentiation data set (Jerber et al. 2021). Further details on the benchmark design are provided in the Supplemental Methods.

Control gene selection bias in gene signature scoring

To evaluate the impact of control gene selection, we compared the approaches used in ANS, Seurat, Seurat_AG, and Seurat_LVG across diverse cell subtypes within the PBMC data set. Specifically, we evaluated the behavior of the scores for the top 8% genes (about 900 genes) with the highest average expression values corresponding to the last two bins of the Tirosh-based methods (Methods) (Fig. 1C; Supplemental Fig. S2). On homogeneous data sets, scores of single-gene signatures are expected to be distributed around zero; therefore, a deviation of the score from zero indicates higher control-gene selection bias. Although we observed minimal score biases in the last but one expression bin, the last bin showed large biases in the score values for all published methods owing to the significant variance of the genes’ average expressions (Fig. 1C). Whereas Seurat, Seurat_AG, and SCANPY first underestimated two-thirds of the gene scores and overestimated one-third, Seurat_LVG overestimated almost all scores in the last expression bin (Fig. 1C). Importantly, all methods showed notable score biases for the top 50 expressed genes; these genes get typically excluded from gene signatures scored by ANS owing to the impossibility of constructing a valid control gene set. Overall, ANS induced the smallest bias, confirming the proposed strategy for control gene selection as the optimal one.

The influence of cell-type proportions, batch effects, signature length, and inclusion of irrelevant genes on gene signature scoring

The mean control values and, thus, the choice of control genes depend on the data set composition of the tissue or sample analyzed. Consequently, when utilizing methods that employ control genes, we anticipated greater score variability when scoring was conducted on a per-sample basis compared with assessing all cells collectively.

To verify this hypothesis, we calculated scores of 100-gene signatures associated with the malignant cell phenotype of CRC and ESCC data sets (Methods). We compared the distributions of the mean scores per sample and cell type under two scenarios: when scoring was conducted on a per-sample basis and when all samples are together (Fig. 1D; Supplemental Fig. S3). We observed a reduction of score variance for all scoring methods and both CRC and ESCC data sets, except for Seurat_LVG and UCell in the ESCC data set, when scoring was performed on cells of all the samples together (Supplemental Fig. S4), leading to more comparable score ranges over all samples. Additionally, we assessed how scoring on all samples together contributed to diminishing the batch effects for the Tirosh-based scoring methods. For this, we evaluated the differences in score distributions across two batches representing different sequencing chemistry, called SC3Pv2 and SC3Pv3 (Fig. 1D; Supplemental Fig. S5; Supplemental Table S2). Batch effects significantly affected score values of all Tirosh-based scoring methods when the scoring procedure was performed on a per-sample basis. In contrast, we observed decreased distribution shifts between chemistry-related batches when scoring on all samples. Therefore, our results suggest that scoring should be performed on all cells in the data set simultaneously for all methods that use control gene sets, including ANS, as this strategy mitigates both biological and technical batch effects while ensuring comparable score ranges across samples.

Next, we assessed scoring methods’ robustness to changes in the signature length using the CRC and ESCC data sets. Methods were compared based on the minimum number of signature genes needed to perfectly classify malignant and nonmalignant cells (Methods) (Supplemental Figs. S3, S6). Most methods, including ANS, achieved perfect classification with only 12–15 marker genes, selected as top differentially expressed genes between malignant and nonmalignant cells (Supplemental Figs. S3, S6). Jasmine_OR required about 11 times more genes for an accurate cell-type annotation in CRC. The small length of these effective signatures might be owing to the relative ease of distinguishing malignant from nonmalignant cells, that is, that the signatures contain genes with high specificity for the two populations.

Further, to test the methods’ robustness to the presence of irrelevant genes in a signature, we progressively replaced genes in a “pure” 100-gene malignant signature with random genes (Methods). All methods, including ANS and except Seurat_LVG, maintained high performance (area under the ROC curve [AUCROC] > 0.9) with up to 85% noise, requiring only about 15 informative genes in both data sets tested (Supplemental Figs. S3, S6; Supplemental Tables S3, S4). In this case, ANS reached on average 0.88 accuracy, showing the first- and third-best performance for CRC and ESCC, respectively, whereas Seurat_LVG consistently underperformed, needing at least twice as many correct genes for comparable results in both data sets.

We also evaluated computational time requirements of the eight methods as a function of both the number of cells and the number of genes in the signature (Supplemental Fig. S7). All methods exhibited linear scalability with data set size. For short signatures (up to 10 genes), Seurat, SCANPY, and ANS were the fastest (<10 sec for 150,000 cells), whereas Seurat and SCANPY remained the quickest for longer signatures (1000 genes). In contrast, UCell and JASMINE required approximately seven times longer than Seurat and SCANPY, whereas ANS was about three times slower in this case. To evaluate the practical impact of the increased computation time of ANS, we assessed the length distribution of 35,134 gene signatures from the MSigDB human collection (Supplemental Fig. S8; Liberzon et al. 2015). The majority (64%) of analyzed signatures had lengths within the range of 100 genes, indicating that in practice ANS offers efficient computation, within seconds to a few minutes depending on the data set size, while providing increased performance.

Assessing score-range comparability for score-based cell-type and cell-state annotation

Reliable score-based annotation requires scores with high information content and comparable ranges. To assess the information quantity and score comparability produced by ANS and other methods, we analyzed the accuracy of cell-type and cell-state annotation across four cancer data sets (BRCA, HGSOC, LUAD, and sCC) and four PBMC cell subsets, using specific nonoverlapping gene signatures (Methods) (Supplemental Table S5). We accessed the fairness of produced scores by evaluating the accuracy of an unsupervised cell annotation, based on the argmax assignment of cells to states from signature scores, against the ground-truth cell labels. Additionally, we evaluated the information quantity of scores produced by each scoring method using a cross-validated supervised logistic regression predicting cell state based on signature scores (Methods).

We observed that score-based cell annotation was accurate for distinguishing cell types (B cells, monocytes, and NK cells) in the PBMC data set using scores produced by most methods (Fig. 2A). All methods achieved high balanced accuracies ranging from 0.921 (Jasmine_OR) to 0.999 (SCANPY and UCell) and F1-scores from 0.958 (Seurat) to 0.999 (UCell), demonstrating robust performance in distinguishing these distinct cell types, with ANS presenting the second-best balanced accuracy and F1-score (0.997 and 0.996, respectively) (Supplemental Table S6). Indeed, most scoring methods produced comparable score ranges, with the highest scores observed for cells of matched cell type; JASMINE was an exception, showing overlapping ranges for signature scores across B cells (Fig. 2A).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Comparative analysis of scoring methods for cell-type and cell-state annotation. Only nonoverlapping cell-type- or cell-state-specific signatures were used. (A) Score distributions for cell-type-specific signatures (B cells, monocytes, and NK cells) separated by true cell-type annotations, calculated for each scoring method. (B) Score distributions for B cell subtype signatures separated by true cell subtypes, calculated by each scoring method. (C) Row-normalized confusion matrix of B cell subtype annotation based on the highest scores. (D) Relationship between hard-labeling performance and score information quantity in cancer and PBMC data sets. Scatterplots show balanced accuracy (x-axis) against score information quantity (y-axis) for various scoring method–data set combinations. Balanced accuracy quantifies hard-labeling performance; a score information quantity indicates the scores’ effectiveness in subtype classification. The diagonal line indicates perfect metric alignment, with vertical distances from this line representing scale imbalance. (E) Quantitative analysis of scale imbalance across scoring methods and tissue types (cancer and PBMC). The mean and standard deviation of scale imbalance for each method are shown. Scale imbalance is the absolute difference between score information quantity and balanced accuracy in direct label assignment. The method with the lowest mean scale imbalance, indicating optimal consistency between information content and labeling accuracy, is highlighted in bold. The intensity of the red background corresponds to the increase in scale imbalance. (F) Cell-state and cell-type annotation performance overview for all eight data sets and scoring methods.

However, when we applied the score-based annotation approach to assign cell states—naive, intermediate, and memory B cells—we observed limitations in using scores produced by all methods except ANS. Indeed, ANS was the only approach to produce informative and comparable score ranges for the evaluated cell states (Fig. 2B). Tirosh-based methods (Seurat, Seurat_AG, Seurat_LVG, and SCANPY) consistently scored highest for the B memory cell signature, regardless of the true cell state, resulting in poor labeling performance (Fig. 2C; Supplemental Fig. S9). UCell showed higher B memory cell scores than B intermediate ones when applied to B intermediate cells. Conversely, JASMINE scored highest for the B intermediate cell signature in B intermediate cells but failed to produce distinguishable score ranges for other cell states.

To comprehensively evaluate observed patterns, we expanded our analysis with four data sets of malignant cells with ground-truth assignments of cells to malignant cell states (BRCA, HGSOC, LUAD, and sCC) and two additional PBMC cell subsets comprising CD4 (CD4⁺) and CD8A (CD8⁺) cell-state annotations. Across these six data sets, ANS maintained high performance, providing unbiased gene signature scores for cell annotation, achieving the highest balanced accuracy in five cases (BRCA: 0.560; HGSOC: 0.638; cSCC: 0.843; CD4⁺ T cells: 0.645; CD8⁺ T cells: 0.823) and performing competitively in LUAD (0.897 vs. SCANPY's 0.914) (Fig. 2F; Supplemental Figs. S9, S10; Supplemental Table S6).

Further, to assess score information content, we compared the unsupervised score-based cell classification accuracy with a supervised cross-validated linear classifier approach. We define the latter's accuracy as the information quantity of the gene signature scores, representing the feasibility of accurate cell annotation (Methods). We define the difference between the information quantity and the unsupervised accuracy as the scale imbalance, which indicates score ranges that vary substantially across cell types (Supplemental Methods). Across the eight data sets (seven for cell-state and one for cell-type annotations), ANS consistently ranked among the best-performing methods, showing the closest alignment of balanced accuracy of the unsupervised score-based annotation with information quantity (Fig. 2D). Notably, ANS demonstrated substantially better performance for sCC, CD8⁺, and B cell subtype classification compared with other scoring methods. Moreover, across both cancer and PBMC data sets, ANS exhibited minimal scale imbalance, showing the smallest standard deviation, suggesting comparable score ranges (Fig. 2E; Supplemental Table S6). Although these results are based on nonoverlapping gene signatures, similar outcomes were observed for signatures that shared certain genes (Supplemental Figs. S11, S12; Supplemental Table S7).

Assessing score-based annotation of rare and transitioning cell populations

To test the performance of scoring methods in a biological context involving state transitions and rare cell populations, we assessed the feasibility of score-based cell labeling using a neuronal differentiation data set from (Jerber et al. 2021). The data set is composed of 215 human induced pluripotent stem cell (iPSC) lines undergoing differentiation toward a midbrain neural fate (Fig. 3; Supplemental Figs. S13–S15). Specifically, single cells were profiled at three maturation stages (days 11, 30, and 52), capturing a continuum of differentiation states from progenitor-like to mature neuronal states, including floor plate progenitors (FPP), proliferating floor plate progenitors (P_FPP), neuroblasts (NB), dopaminergic neurons (DA), serotonergic-like neurons (Serts), proliferating serotonergic-like neurons (P_Serts), astrocyte-like cells (Astro), ependymal-like 1 (Epen1), ependymal-like 2 (Epen2), and three uncharacterized neuronal cell subtypes (U_Neur1–3).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Benchmarking of score-based cell labeling using neuronal differentiation single-cell data set. (A) Overall balanced accuracy of score-based labels assigned based on the 12 cell signatures. (B) Cell-type-specific balanced accuracy across all benchmarked scoring methods. (C) Cell-type proportion distribution across the analyzed data set, including floor plate progenitors (FPP), proliferating floor plate progenitors (P_FPP), neuroblasts (NB), dopaminergic neurons (DA), serotonergic-like neurons (Serts), proliferating serotonergic-like neurons (P_Serts), astrocyte-like (Astro), ependymal-like 1 (Epen1), ependymal-like 2 (Epen2), and three unknown neuron groups 1–3 (U_Neur1, U_Neur2, U_Neur3). (D) PHATE embeddings generated based on the ANS scores of the 12 cell types. The plots are colored by ground-truth cell-type annotations, sequencing time points (days), and pseudotime estimated from the embeddings.

We observed substantial variation in the score ranges across tested cell-state signatures, with most methods, except for ANS and JASMINE, showing consistently high scores for FPP and NB signatures regardless of the underlying cell type (Supplemental Fig. S13). ANS achieved overall the highest balanced accuracy of score-based annotation across all 12 cell types, with other Tirosh-based approaches showing substantially worse performance (Fig. 3A). Of note, although the data set contains several rare cell populations, including NB, P_Serts, Epen1, and U_Neur1–2, each representing only 0.2% to 1.3% of all cells, ANS-based annotation was the best for three out of five of those cell types and second best for the remaining two (Fig. 3B,C). Moreover, ANS-based annotations showed the highest pairwise balanced accuracy for closely related populations (Epen1 and Epen2, FPP and P_FPP, DA and U_Neur1, and Sert and U_Neur3) (Fig. 3B; Supplemental Fig. S14), highlighting its ability to discriminate transcriptionally similar cells.

We further explored ANS-derived scores to verify if they accurately capture continuous cellular transitions. Specifically, we used potential of heat diffusion for affinity-based transition embedding (PHATE) (Moon et al. 2019), a dimensionality reduction method that preserves both local and global structure in the data, to project analyzed cells into a low-dimensional manifold. The resulting ANS score–based embeddings accurately reconstructed the expected differentiation trajectories and were consistent with both the ground-truth cell-type labels and sequencing time points (Fig. 3D). Importantly, ANS score–based embeddings presented clearer lineages compared with embeddings derived from raw gene expression data (Supplemental Fig. S15). Although these results further illustrate the practical utility of ANS for capturing meaningful biological signals, a systematic comparison with other dimensionality reduction approaches warrants further investigation. Nevertheless, ANS scores effectively denoise the signal and provide robust features for trajectory inference of cells in continuous states.

Assessing the effect of removing highly expressed signature genes for score-based cell annotation

As described in the previous section, all Tirosh-based approaches suffer from bias when control genes are selected for highly expressed signature genes. To mitigate this issue, ANS filters out signature genes that fall within the top c/2 expressed genes, where c is the size of the control gene set. To directly evaluate the benefit of our control gene selection strategy, we performed additional benchmarking experiments on five PBMC subsets, including both broad cell types (B cells, NK cells, and monocytes) and finer cell states (B cell, CD4⁺ T cell, and CD8⁺ T cell subtypes), under five different signature gene filtering strategies. Specifically, we applied (1) no filtering; (2) removal of the top 50 and (3) top 100 expressed genes, corresponding to ANS-based filtering with control sets equal to 100 or 200 genes, respectively; and (4) removal of genes in the top expression bin when genes were split into 25 bins or (5) into 50 bins. We compared the balanced accuracy of score-based cell assignments using ANS and the four Tirosh-based approaches across all filtering settings.

Although all methods achieved almost perfect accuracy for broad cell-type identification (i.e., B cells, NK cells, and monocytes) across the filtering strategies, we observed considerable differences for cell-state annotation, with ANS consistently achieving the highest balanced accuracy (Supplemental Fig. S16). Tirosh-based methods reached comparable performance only after removing signature genes within the top 100 expressed genes, but their accuracy declined again when filtering by top expression bin, highlighting their sensitivity to the selected filtering threshold.

Notably, our analysis showed that a large proportion of signature genes often falls within the top expression bins (Supplemental Fig. S17; Supplemental Tables S8–S12), suggesting that bin-based filtering may eliminate too many biologically informative genes and thereby reduce the robustness of score-based cell annotation. In contrast, ANS-based annotation reached optimal accuracy when excluding only the 50 most highly expressed genes and maintained high accuracy even without gene filtering, directly demonstrating the benefits of the adaptive control gene selection. Together, these results confirm that our proposed control gene selection strategy enhances score robustness for cell annotation beyond what can be achieved through simple filtering of highly expressed signature genes.

Devising cancer cell–specific EMT signature: a case study for the ANS application

To showcase the practical utility of ANS, we used it to infer a pancancer malignant cell-specific EMT signature. The difficulty in constructing such a gene signature comes from the high similarity of transcriptional states of malignant EMT-exhibiting cells with malignant non-EMT cells, on the one hand, and with nonmalignant mesenchymal cells, such as cancer-associated fibroblasts (CAFs), on the other hand.

Indeed, when we explored the behavior of previously published EMT signatures using the ANS scoring method on single cells from the ESCC (Zhang et al. 2021a), lung adenocarcinoma by Xing et al. (2021) (LUAD_Xing), and CRC (Pelka et al. 2021) and BRCA (Wu et al. 2021) tumors, we observed that the published pancancer EMT signatures, such as the Hallmark EMT (Subramanian et al. 2005; Liberzon et al. 2015), the pEMT signature (Barkley et al. 2022), and five other EMT signatures (Gröger et al. 2012; Tan et al. 2014; Mak et al. 2016; Foroutan et al. 2017; Hollern et al. 2018), resulted in extremely high scores in CAFs in addition to the mesenchymal-like (MES-like) malignant cells (Fig. 4A; Supplemental Fig. S18), whereas the signatures of Foroutan et al. (2017), Gröger et al. (2012), Mak et al. (2016), and Tan et al. (2014) resulted in high scores for pericytes in CRC and ESCC, and the former also scored high in endothelial cells for all data sets (Supplemental Fig. S18).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Application of ANS to devise EMT signature specific to malignant cells. (A) Score distributions of the Hallmark EMT, the pEMT, and our proposed ESCC-specific and LUAD_Xing-specific cancer EMT gene signatures in different cell types of ESCC, LUAD_Xing, CRC, and BRCA. (B) Venn diagram for the overlap in gene lists for the Hallmark EMT, pEMT signatures, and the signature we designed. (C) Association of ESCC-specific and LUAD_Xing-specific cancer EMT signature scores and histological subtypes in TCGA (only cancer types with at least one histotype are included). Significant associations (Q value < 0.05) are represented by dots above the dotted line (yellow). The FDR-adjusted P-values from the two-sided Kruskal–Wallis test were −log₁₀-transformed and sorted in decreasing value, with a higher y-value indicating higher significance. (D) Multivariate Cox survival analysis with disease stage, patient age, and cancer cell–specific EMT signature scores calculated in bulk TCGA RNA-seq data. The model was stratified by cancer type. Variable-specific P-value significance, (***) P-value < 0.001, (**) P-value < 0.01, (*) P-value < 0.05. (E) Association between the cancer cell–specific EMT scores in bulk TCGA tumors and patient treatment response. Association was calculated with Spearman's rank correlation between the scores and treatment response ranked one (complete response) to four (clinical progressive disease; Methods).

To devise a cancer cell–specific EMT signature, we scored established EMT signatures using ANS in the ESCC and LUAD_Xing data sets. Briefly, we selected genes that were significantly upregulated in cancer cells with the highest EMT scores (classified as cancer EMT cells) compared with CAFs in the two data sets (Methods). Overall, the signature included many known EMT-related genes, as well as genes previously not included in canonical EMT signatures (Fig. 4B; Supplemental Table S13). By construction, the devised signature scored highly for the MES-like malignant cells in the training data sets (ESCA and LUAD_Xing), and clearly distinguished cells in the validation data sets (CRC and BRCA) (Fig. 4A; Supplemental Table S14). Specifically, the accuracy of our gene signature discriminating MES-like malignant cells versus all remaining cells measured with the area under the precision-recall curve (AUCPRC) on CRC and BRCA data sets was 0.701 and 0.511, respectively; it was slightly outperformed on CRC only by the gene signature proposed by Hollern et al. (2018) (AUCPRC 0.707) (Supplemental Table S14). However, this signature performed suboptimally in our ESCC and LUAD_Xing training data sets (AUCPRC 0.028 vs. 0.592 of our gene signature in ESCC and 0.244 vs. 0.513 in LUAD_Xing).

Calculation of cancer cell–specific EMT phenotypes using bulk RNA-seq data offers a promising alternative to assess their link with clinical characteristics in cases when corresponding single-cell data are not available. To verify the applicability of our devised signature in bulk tumors, we used RNA-seq from The Cancer Genome Atlas (TCGA) data sets and scored both our cancer-cell EMT and hallmark EMT signatures.

Scores of the ANS-derived signature showed higher correlation with tumor purity than the hallmark EMT ones in 29/33 cancer types (Supplemental Fig. S19), indicating higher sensitivity of our signature to capture EMT signals within the malignant compartment than with the nonmalignant one. Next, we compared the score distributions across cancer types (Supplemental Fig. S20). Hematological malignancies, including acute myeloid leukemia (LAML) and diffuse large B cell lymphoma (DLBC), showed the lowest scores, consistent with their nonepithelial origin and the absence of EMT in liquid cancers. In contrast, epithelial cancers, such as head and neck squamous cell carcinoma (HNSC), pancreatic adenocarcinoma (PAAD), and cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC) showed the highest scores, aligning with the well-established involvement of EMT in these tumor types (Lee et al. 2008; Puram et al. 2017; Yang et al. 2023). The Hallmark EMT signature showed much higher scores in sarcomas (SARC) than our cancer cell–specific EMT signature, an observation consistent with the mesenchymal rather than epithelial origin of this tumor type. The elevated Hallmark EMT scores in SARC therefore likely reflect the inclusion of stromal and mesenchymal gene expression in bulk tumor samples, further supporting high specificity of our cancer cell–specific EMT signature. Finally, we observed significantly higher cancer-cell EMT scores in primary melanoma tumors (P-value < 0.0001) (Supplemental Fig. S21), an observation consistent with the hypothesis that EMT is necessary to gain motility for invasion, whereas mesenchymal-to-epithelial transition (MET) is necessary for new site establishment (Tsai et al. 2012). In the same data set, hallmark EMT scores showed no significant differences, further highlighting higher specificity of our proposed signature (Supplemental Fig. S21).

Next, we used TCGA bulk RNA-seq data to investigate clinical significance of the proposed cancer cell–specific EMT signature. We found our signature scores significantly associated with the histological subtype in 10/24 cancer data sets for which at least two subtypes were available (Fig. 4C; Supplemental Table S15). The scores were significantly associated with worse patient outcome in multivariate Cox regression analysis stratified by cancer type, with patient age and tumor stage as covariates (P-value < 0.05) (Fig. 4D). The association remained significant also in a model including hallmark EMT scores, providing evidence that our devised signature captures additional clinically relevant signal beyond that of the hallmark EMT signature (Supplemental Fig. S22). Finally, we assessed the association between patient treatment response and the calculated cancer cell–specific EMT scores to test the hypothesis that EMT-like cancer cells within tumors lead to higher treatment resistance (Shibue and Weinberg 2017; Brabletz et al. 2018). We found positive association between the scores and worse treatment response of patients treated with chemotherapy in 15 of 25 cancer types, as well as in two of eight cancer types for which patients were treated with immunotherapy (Fig. 4E). Together, these results highlight the biological significance of the proposed signature and further demonstrate potential of ANS in revealing biologically and clinically relevant information.