Network Perturbation Analysis Identifies Three Molecular Subtypes With Distinct Clinical Outcomes and Therapeutic Vulnerabilities in Endometrial Cancer

We performed network perturbation analysis on transcriptomic data from endometrial cancer patients across multiple cohorts and identified three molecular subtypes with distinct clinical outcomes. The most aggressive subtype, C1, exhibits comprehensive molecular and immunological alterations, including dysregulated hormone signaling, immune exclusion, and specific pathway activation. Multi-dimensional profiling integrating transcriptional, immunological, mutational, and pharmacological analyses revealed subtype-specific therapeutic vulnerabilities, highlighting the potential for targeted intervention in high-risk endometrial cancer patients. Endometrial cancer (EC) represents a major global health burden characterized by substantial molecular and clinical heterogeneity [1]. Although existing molecular classifications, including The Cancer Genome Atlas (TCGA) molecular subtypes, have advanced our understanding of EC diversity, these approaches primarily rely on static genomic features—mutations, copy number alterations, and gene expression [2, 3] —treating molecular components as independent entities while overlooking the interconnected nature of biological systems. Accumulating evidence demonstrates that tumorigenesis involves systematic perturbation of gene regulatory networks rather than isolated molecular alterations [4]. These network-level disruptions encompass not only individual gene expression changes but also fundamental rewiring of regulatory interactions governing cellular phenotypes. Network perturbation analysis, which quantifies the deviation of tumor-associated molecular interactions from homeostatic states in normal tissues, has successfully identified clinically relevant subtypes across multiple cancer types [4-6]. These precedents suggest that network-based approaches could reveal previously unrecognized heterogeneity in endometrial cancer. We conducted a comprehensive network perturbation analysis of endometrial cancer using 783 samples from 4 cohorts and identified three molecular subtypes with distinct survival outcomes and therapeutic vulnerabilities. Analysis of 799 significantly perturbed edges and 694 nodes revealed that the most aggressive subtype (C1) exhibits loss of hormone receptor signaling, an immune-excluded microenvironment, and specific sensitivities to microtubule-targeting agents and proteasome inhibitors, whilst showing resistance to DNA-damaging agents. Integration of transcriptional, immunological, and pharmacological profiles uncovered subtype-specific therapeutic vulnerabilities with immediate translational potential. We processed RNA-seq data from endometrial carcinoma tumor samples (n = 537) and matched normal endometrial samples (n = 35) retrieved from the TCGA database. Validation cohorts included GSE17025, GSE146889, and GSE106191 from the GEO database. Following established methodologies [6], we conducted network perturbation analysis on the TCGA cohort. Specifically, genes were ranked by expression level in each sample, and for each gene, we calculated a perturbation score by comparing its expression rank with the average rank of its interacting partners' interaction networks. These rank-based perturbation scores were then normalized using Z-score transformation relative to the mean and standard deviation observed in normal samples. Normal tissue samples, with their stable interaction patterns, served as the baseline for measuring tumor-specific deviations, allowing us to quantify the extent to which each gene's network context was disrupted in tumor samples. This approach generated sample-level network perturbation profiles that capture coordinated dysregulation patterns across the entire gene regulatory network (Figure 1A). Whereas normal samples showed minimal perturbation, tumor samples exhibited markedly elevated perturbation levels with diverse molecular disruption patterns (Figure 1B,C). Using the intersection of the top 15,000 interactions ranked by smallest p values and largest standard deviations, we identified 799 significantly perturbed edges and 694 nodes that differentiated tumor from normal samples. Functional enrichment analysis revealed that perturbed genes were primarily enriched in therapeutic response pathways, including proteasome function, EGFR tyrosine kinase inhibitor resistance, ubiquitin-mediated proteolysis, neurotrophin signaling, and endocrine resistance (Figure 1D). This convergence on protein homeostasis and therapeutic resistance mechanisms suggests that network disruptions in EC predominantly affect drug response pathways—a finding with immediate treatment implications. Unsupervised clustering of the perturbation network revealed three distinct patient subgroups. Evaluation across cluster numbers (K = 2–9) demonstrated optimal separation at K = 3 (Figure 1E). The robustness of this three-cluster solution was confirmed through multiple metrics: consensus scores exceeded 0.75 for all groups (Figure 1F), cumulative distribution function curves indicated substantial subtype stability (Figure 1G), and proportion of ambiguous clustering analysis yielded minimal ambiguity at k = 3 (Figure 1H). Principal component analysis further validated the clustering structure, revealing clear spatial segregation among the three groups (Figure 1I). Kaplan–Meier analysis demonstrated statistically significant differences in overall survival, with C1 displaying the poorest prognosis (Figure 1J). To validate reproducibility, we identified 300 signature genes characteristic of each cluster. Application of the nearest template prediction algorithm [7] enabled subtype assignment in three additional cohorts, with consistent transcriptional profiles across all four cohorts confirming classification robustness (Figure 1K). Finally, we systematically explored the associations between clinical characteristics and the three subtypes. Importantly, the C1 subtype showed significant enrichment in high-risk clinical features, including advanced clinical stage (Figure S1A), high histological grade (Figure S1B), microsatellite instability (MSI)-low status (Figure S1C), and copy number-high subtype (Figure S1D). Using a predictive model for progesterone receptor (PR) status developed in our previous work [8], we found that C1 patients predominantly exhibited low PR status. Collectively, these associations indicate that our functional phenotype-based classification captures biologically aggressive tumor characteristics that align with established prognostic indicators. The C1 subtype displayed distinctive molecular characteristics that explain its aggressive clinical behavior. Gene set enrichment analysis (GSEA) revealed three molecularly distinct subtypes of endometrial carcinoma with unique biological characteristics: the C1 subtype characterized by cell cycle dysregulation and DNA repair abnormalities (Figure S2A), the C2 subtype featuring enhanced invasiveness and immune pathway activation (Figure S2B), and the C3 subtype displaying metabolic reprogramming with prominent mitochondrial function (Figure S2C). The convergence of pathway enrichment patterns and drug sensitivity profiles reinforces the biological validity of this classification system and provides a robust foundation for precision medicine approaches in endometrial cancer. C1 tumors showed significantly suppressed expression of hormone response markers ESR1, PGR, FOXA2, and IHH, indicating loss of hormone sensitivity and predicting endocrine therapy resistance [8, 9]. Conversely, this subtype exhibited marked upregulation of PDK1 and CTNNB1, revealing metabolic reprogramming toward glycolysis and Wnt/β-catenin pathway activation [10] (Figure 2A). Furthermore, differential expression analysis (Table S1) revealed that C1 exhibited marked upregulation of mitochondrial respiratory chain components, including multiple complex I subunits (NDUFS7/8, NDUFA2/3/11/13, NDUFB7/10), complex III components (UQCRQ, UQCR11, UQCRC1), and complex IV (COX4I1), alongside ATP synthase (ATP5F1D). Concurrent elevation of metabolic regulators (ADIPOQ, SIRT6), antioxidant machinery (GPX4), and ribosomal proteins suggested enhanced oxidative phosphorylation and protein synthesis capacity. This expression profile may be characteristic of metabolically active cancer-associated fibroblasts and is consistent with TGF-β-driven stromal activation [11]. The metabolic reprogramming observed in C1 aligns with extracellular matrix remodeling and physical barriers that restrict T cell infiltration, providing a mechanistic basis for the immune-excluded phenotype. Immune microenvironment analysis revealed striking patterns. Quantitative assessment using single-sample GSEA (ssGSEA) demonstrated that C1 exhibited the lowest infiltration of antitumour effector populations, including activated CD8⁺ T cells (Figure 2B). Paradoxically, despite showing elevated Type I interferon response signatures and increased activated dendritic cell infiltration, C1 displayed the poorest overall immune infiltration—a finding that warrants mechanistic consideration. This subtype showed decreased expression of immunomodulatory molecules, including HLA class I molecules, TLR4, CTLA-4, and PD-L1 (Figure 2C). ESTIMATE algorithm [12] analysis confirmed these findings, demonstrating that C1 possessed the lowest immune and stromal scores. The Tumor Immune Dysfunction and Exclusion algorithm [13] analysis confirmed minimal cytolytic activity, the lowest immune dysfunction, and the highest level of immune exclusion (Figure 2D). This apparent contradiction can be explained by the “double-edged sword” effect of chronic Type I interferon signaling. Whereas acute interferon responses promote beneficial immune activation, chronic interferon stimulation induces immunosuppressive effects through multiple mechanisms: promoting terminal T cell exhaustion via metabolic dysfunction and lipid peroxidation, upregulating immunosuppressive molecules including IL-10 and PD-L1, and decreasing costimulatory molecules (CD80, CD86, CD40) on dendritic cells whilst upregulating inhibitory receptors. This creates a state of “frustrated immune activation,” wherein immune cells are recruited but functionally paralyzed. The markedly low MSI score in C1, consistent with decreased neoantigen burden and diminished HLA class I expression, further reinforces an immunologically “cold” phenotype where chronic interferon-mediated immunosuppression prevents effective antitumour immunity despite initial immune cell recruitment. These findings suggest that C1 patients may require combination immunotherapy targeting multiple immunosuppressive pathways rather than single-agent checkpoint blockade. Transcription factor activation drives cellular phenotype formation. We observed significantly increased activities of E2F4, TFDP1, STAT2, E2F2, SNAI2, MNT, NR2C2, IRF9, and TEAD2 (Figure 2E). This profile indicates simultaneous activation of cell cycle progression (E2F4/E2F2/TFDP1) [14] and epithelial–mesenchymal transition programmes (SNAI2) [15], suggesting coordinated proliferative and invasive capabilities. Mutational landscape and pharmacological profiling provided critical therapeutic insights. C1 tumors harbored the highest TP53 mutation frequency but lowest PTEN and CTNNB1 mutation rates, reflecting unique molecular evolution (Figure S3A). Drug screening using the Genomics of Drug Sensitivity in Cancer database revealed subtype-specific therapeutic vulnerabilities (Figure S3B). C1 demonstrated marked sensitivity to mitotic inhibitors and cell cycle-targeting agents. Notably, sensitivity to paclitaxel (microtubule stabilizer) [16] supports taxane-based chemotherapy as frontline treatment. Vulnerability to BI-2536 (PLK1 inhibitor) [17] and YM155 (survivin inhibitor) [18] reflects dependency on mitotic machinery. Furthermore, sensitivity to bortezomib (proteasome inhibitor) and MS-275 (HDAC inhibitor) suggests reliance on protein homeostasis and epigenetic regulation. Conversely, C1 showed resistance to DNA-damaging agents (mitomycin C) and growth factor pathway inhibitors (erlotinib, bexarotene), indicating intact DNA repair mechanisms and independence from conventional growth signaling—crucial information for avoiding ineffective treatments. These findings provide crucial insights for clinical trial design, suggesting that C1 subtype patients may benefit from therapies targeting mitotic machinery, such as paclitaxel or BI-2536, whilst avoiding DNA-damaging agents or growth factor inhibitors. Stratified treatment strategies could therefore optimize therapeutic efficacy by tailoring treatments based on the specific vulnerabilities and resistances observed in this subtype. In summary, our study establishes a novel molecular classification strategy for endometrial cancer. This classification system shows significant correlations with, yet distinctions from, previous classification schemes. However, despite validating the stability of this classification across multiple public cohorts, further validation in larger real-world clinical cohorts is necessary. Moreover, future studies will require additional experimental evidence to identify novel therapeutic targets and effective treatment strategies based on this classification. This network perturbation analysis establishes a robust molecular classification system for endometrial cancer with immediate clinical applications. The identification of C1 as a high-risk subtype with specific therapeutic vulnerabilities enables rational treatment selection. Molecular profiling at diagnosis could rapidly identify patients requiring intensified treatment with effective agents whilst avoiding resistant therapies. This systems-level approach represents a significant advance toward precision oncology in endometrial cancer, providing a practical framework for personalized treatment strategies. Prospective clinical trials are warranted to validate these findings and test subtype-specific therapeutic approaches. RNA-seq count matrices underwent log2 transformation after adding a pseudocount to TPM-normalized values (log2[TPM + 1]). Microarray probe signals were consolidated at the gene level by averaging expression values when multiple probes mapped to identical gene symbols. Cross-array normalization was performed using the limma package's normalizeBetweenArrays function. Each data set underwent independent processing without inter-data set batch correction. For patient stratification, we implemented consensus clustering using Consensus Cluster Plus [19] (v1.0) with k-means partitioning and Euclidean distance metrics. We explored cluster solutions ranging from K = 2 to K = 9, performing 1000 bootstrap iterations with 80% feature sampling per iteration to ensure statistical robustness. Clustering stability and optimal K selection were evaluated through consensus matrices and cumulative distribution functions visualized using ggplot2. This unsupervised approach enabled identification of molecularly distinct sepsis subtypes based on their characteristic network perturbation signatures rather than individual gene expression patterns. We validated the identified molecular subtypes across external datasets using Nearest Template Prediction (NTP) [7], a robust single-sample classifier that accommodates cross-platform heterogeneity and provides probabilistic confidence estimates for multi-class assignments. For each subtype, we first derived molecular signatures comprising the 300 most discriminative genes identified through differential expression analysis with limma. These subtype-specific gene signatures served as classification templates for patient assignment in independent validation cohorts. The NTP algorithm computed similarity scores between individual patient profiles and each subtype template, enabling probabilistic classification with associated confidence metrics. Default NTP parameters were retained to ensure reproducibility. Classification results and associated metadata were visualized as annotated heatmaps using ComplexHeatmap [20], facilitating assessment of subtype conservation across diverse patient populations and experimental platforms. All computational analyses were performed in R (v4.4.0). For group comparisons, we first evaluated data distributions using Shapiro–Wilk normality tests and Levene's test for variance homogeneity. Parametric tests (Student's t-test for two-group comparisons, one-way ANOVA for multiple groups) were applied to normally distributed data with equal variances. Non-parametric alternatives (Wilcoxon rank-sum test and Kruskal–Wallis test) were employed when parametric assumptions were violated. Chi-square tests evaluated associations between categorical variables. Statistical significance was defined as p < 0.05 for all analyses. Jie Liu: conceptualization, formal analysis, writing – original draft, writing – review and editing. Yuxuan Song: methodology, software, formal analysis. Yiqin Wang: methodology, software, formal analysis. Shaolei Li: formal analysis. Yuanyuan Ma: formal analysis. Meng Zhang: investigation, visualization. Keran Sun: formal analysis. Yibo He: investigation. Hengrui Liu: project administration, conceptualization, formal analysis, writing – original draft. Hao Chi: project administration, conceptualization, formal analysis, writing – original draft. Jingyuan Ning: project administration, conceptualization, formal analysis, writing – original draft, writing – review and editing. We acknowledge Sogen Biotechnology, Zhengzhou, Henan Province, and the SolvingLab team for their technical support and assistance. The authors declare no conflicts of interest. The authors have nothing to report. All data are available in the main Excel or the supporting materials. Transcriptomic and clinical data for all patients were obtained from TCGA and GEO database. Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.