Methylation data from 312 cervical cancer samples, generated using the Illumina Infinium HumanMethylation450 BeadChip array, were obtained from the UCSC Xena database (https://xenabrowser.net/datapages/). The $\beta$-value quantified the methylation status of individual probes, ranging from 0 (unmethylated) to 1 (fully methylated). Probes with more than 70% missing values were first excluded from the analysis. The k-nearest neighbors (KNN) algorithm was then applied to impute the remaining missing data. To ensure data quality, unstable genomic regions were removed, including CpG sites on sex chromosomes and those overlapping with single-nucleotide polymorphisms (SNPs). After quality control and normalization, CpGs in promoter regions—defined as 2 kb upstream to 0.5 kb downstream of the transcription start site—were selected. Methylation changes in these regions have been reported to regulate gene expression.

RNA-Seq transcriptome datasets from 309 cervical cancer samples were also obtained from the UCSC Xena database and preprocessed. Finally, 287 samples with complete gene expression profiles and clinical data were selected for further analysis. This cohort was randomly divided into training and testing sets at a 60:40 ratio.

For external validation, the GSE44001 dataset was obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Preprocessing included KNN imputation for missing values, removal of low-expression and low-variance genes, log2 transformation of expression values, and Z-score standardization.

In this study, Cox proportional hazards regression models were used to identify classification characteristics associated with cervical cancer outcomes. Initially, univariate Cox models were applied, incorporating the methylation level of each CpG site together with age, clinical stage, ethnicity, Tumor-Node-Metastasis (TNM) stage, and histological type, to identify risk factors for survival outcomes. CpG sites with statistical significance in the univariate Cox analysis were subsequently entered into a multivariate Cox model, with clinical stage and pathologic M and N stages as covariates to control for potential confounding. All these covariates were statistically significant in the univariate model. Finally, CpGs that remained statistically significant in the multivariate model were retained as subsequent categorical features.

Model construction was performed using the ‘coxph’ function from the survival package (version 3.8.3) in R (version 4.4.2; R Foundation for Statistical Computing, Vienna, Austria). For each CpG site, the multivariate Cox model was defined as follows:

h(t,x)_i = h0(t)exp($\beta$_methymethy_i + $\beta$_stagestage + $\beta$_MM + $\beta$_NN)

Here, methy_i is the methylation level at the CpG site in the sample; M, N, and stage denote the M stage, N stage, and clinical stage, respectively. The terms $\beta$_methy, $\beta$_M, $\beta$_N, and $\beta$_stage represent the regression coefficients for each corresponding variable.

Methylation-based subtypes of cervical cancer were identified using consensus clustering analysis with the ConsensusClusterPlus package (version 1.70.0) in R. The algorithm generated sub-datasets from the original methylation data matrix through random sampling. Each random sample contained 80% of the cases and all feature variables. Each sub-dataset was clustered using the Partitioning Around Medoids algorithm, with the number of clusters (k) set from 2 to 10. The procedure was repeated 100 times. Pairwise consistency between samples was calculated as the proportion of times two samples clustered together across all iterations. These values were compiled into a consensus matrix to quantify clustering stability.

For each k, a delta area plot was generated to assess the incremental increase in the area under the cumulative distribution function (CDF) curve between adjacent k values. When the delta value declined and reached a stable plateau, increasing the number of clusters was considered to have minimal impact on clustering performance.

The consensus matrix heatmap displayed clustering consistency between samples through color intensity, with higher intensity indicating greater stability. For interpretation, DNA methylation clusters were annotated with histological type, TNM stage, clinical stage, race, and age. Heatmaps were created using the pheatmap package (version 1.0.13) in R. The optimal number of cervical cancer molecular subtypes was determined by integrating information from the delta area plot and consensus matrix heatmap. This approach ensured robust and reproducible clustering while minimizing the impact of random variation from single clustering runs.

Survival analysis was performed using the survival package in R. Overall survival among different clusters was visualized with Kaplan-Meier curves. The log-rank test was applied to assess the statistical significance of survival differences between clusters. A p-value 0.05 was considered statistically significant. Lymphocyte infiltration was estimated by calculating the mean expression of cytotoxic T lymphocyte (CTL) marker genes. A list of 547 CD8+ T cell-related genes was obtained from the CIBERSORT database (http://CIBERSORT.stanford.edu) and used for this calculation.

In this study, we employed the Quantitative Differentially Methylated Regions (QDMR) tool (version v1.0; Group of Computational Epigenetic Research, College of Bioinformatics Science and Technology, Harbin Medical University, Harbin, China) identify distinct DNA methylation sites in cervical cancer subgroups. The QDMR quantitative approach was applied to genome-wide methylation profiles to detect differentially methylated regions (DMRs) by adjusting Shannon entropy. This method enabled the identification of specific hypermethylated or hypomethylated sites in cervical cancer clusters [14].

To explore the biological roles of the genes corresponding to the differential CpGs in different methylation clusters, we first mapped the CpGs to their associated genes using the R annotated package “IlluminaHumanMethylationEPICanno.ilm10b4.hg19” (version 0.6.0). When a site corresponded to multiple genes, we retained only the first gene name. After that, the corresponding sets of unique genes in each cluster were extracted separately. Subsequently, we performed functional enrichment analyses, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG), on the genes in each cluster using the “clusterProfiler” package (version 4.14.6). Enrichment results was considered significant when the Benjamini-Hochberg adjusted q-value was below 0.05. For visualization, we retained only pathways with more than one enriched gene and presented the results as bar plots using the ggplot2 package (version 3.5.2).

A Bayesian model was constructed using the training set data and applied to classify samples in the test dataset. External validation was conducted with the GSE44001 dataset from GEO database [15]. The clustering approach from the training phase was applied to assign category labels to all samples in the external validation cohort. The pROC package (version 1.19.0.1) in R was used to perform the receiver operating characteristic (ROC) curve analysis.

To examine variations in the immune microenvironment among methylation-based clusters, multiple algorithms were applied for a comprehensive assessment. First, immune cell infiltration levels were evaluated using the CIBERSORT method, which generated a matrix of 22 immune cell subsets, including B cells, T cells, and macrophages. The relative proportions of these 22 immune cell types were estimated across different clusters. Second, the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm was applied to quantify immune enrichment scores for 28 immune cell types, including subpopulations with immunosuppressive or cytotoxic functions. Third, the ESTIMATE package (version 1.0.13) in R was used to calculate Immune Score, Stromal Score, and Tumor Purity across clusters. Finally, differences in the expression of immune checkpoint genes (Cytotoxic T-Lymphocyte-Associated Protein 4, CTLA-4; Programmed Death 1, PD-1; Programmed Death Ligand 1, PD-L1; Programmed Death Ligand 2, PD-L2) among the methylation clusters were analyzed. To compare variations in immune function indices across the four groups, the Kruskal-Wallis test was performed for all statistical analyses.

A total of 2013 immune-related genes were obtained from the publicly available ImmPort portal (https://www.immport.org/home). Afterwards, the intersections were taken with genes corresponding to differentially methylated CpGs, and 23 differentially methylated immune genes were finally screened out. The average methylation level of each gene was calculated for each sample. The “ggplot2” package in R was then used to generate violin plots for comparing methylation differences among gene clusters.

Three commonly used chemotherapeutic agents—Cisplatin, Paclitaxel, and Docetaxel—were selected for drug sensitivity analysis. The oncoPredict package (version 1.2) in R was used to estimate the half-maximal inhibitory concentration (IC50) for each agent across patient groups [16]. The IC50 represents the concentration at which a drug achieves 50% of its maximal inhibitory effect and serve as a standard measure of drug efficacy. Independent t-tests were then performed to assess differences in drug sensitivity among the molecular clusters.

We downloaded 142 cervical cancer pathology images from The Cancer Genome Atlas (TCGA) through the GDC data portal (https://portal.gdc.cancer.gov/) [17]. The images were preprocessed using Python (version 3.9.7; Python Software Foundation, Wilmington, DE, USA) in PyCharm Professional Edition (version 2024.3.3; JetBrains s.r.o., Prague, Czech Republic), and converted to PNG format. Pathology feature extraction was then performed using CellProfiler (version 4.2.8; Broad Institute of MIT and Harvard, Cambridge, MA, USA) [18]. Automatic image segmentation was carried out with the “IdentifyPrimaryObjects” and “IdentifySecondaryObjects” modules to delineate nuclei and cytoplasmic regions. Six analytical modules—IdentifyPrimaryObjects, IdentifySecondaryObjects, MeasureObjectSizeShape, MeasureTexture, MeasureObjectIntensity and MeasureObjectIntensityDistribution—were applied to extract quantitative cellular features, including size, shape, texture, and intensity. The software generated corresponding measurements for each identified cell. In addition, we compared the distribution of tumor-infiltrating lymphocytes (TILs) among the four methylation clusters [19]. TILs were classified into four types based on morphological characteristics: “Brisk Band-like”, “Brisk Diffuse”, “Non-Brisk Focal”, and “Non-Brisk Multifocal”.

Molecular subtypes related to prognosis in cervical cancer were identified through unsupervised clustering of DNA methylation data from the UCSC Xena database. Prior to clustering, preprocessing steps included removing missing values, imputing data with the KNN method, excluding CpG sites on sex chromosomes or overlapping with single nucleotide polymorphisms (SNPs), performing quality control, normalizing data, and extracting CpG sites located in promoter regions. The dataset was then randomly divided into training and test cohorts in a 60:40 ratio.

Univariate Cox proportional hazards regression was performed for each CpG site in the training cohort of 168 tumor samples, incorporating both CpG methylation levels and relevant clinical variables. In total, 13,945 CpG loci showed statistically significant associations (p 0.05) with patient survival. Clinical stage (p = 0.00083) and the M (p = 0.01459), and N (p = 0.00931) components of TNM stage were also significantly associated with survival. The significant CpGs identified above were then included in a multivariate Cox proportional hazards regression analysis. Significant clinical variables from the univariate analysis were used as covariates to further identify independent prognosis-associated methylation sites. Ultimately, 1842 CpG sites with statistically significant results were selected as the final categorization features (Supplementary Table 1).

Next, consensus clustering was performed using the $\beta$-values of 1842 independent prognostically relevant CpG sites to identify molecular subtypes of cervical cancer DNA methylation. As shown in Fig. 1A, the delta area plateaued after five categories. Therefore, k = 5 was selected as the optimal number of clusters for downstream analysis. In addition, the consensus matrix provided a visual assessment of cluster robustness and aided in determining the optimal number of clusters. A binary color scale was applied, with white indicating a value of 0 and dark blue indicating a value of 1, assuming that entries from the same cluster were arranged adjacently within the matrix. In this configuration, an ideal consensus matrix was represented as a heatmap with prominent blue blocks along the diagonal, indicating perfect sample clustering against a neutral white background. The consensus matrix generated from consensus clustering is shown in Fig. 1B. Fig. 1B, corresponding to k = 5, displayed a clearly delineated five-cluster structure. The heatmap associated with the resulting dendrogram was generated using the pheatmap function and annotated with DNA methylation-based classification, histological type, TNM stage, clinical stage, race, and age (Fig. 1C). Due to the extremely small sample size in Cluster 5 (containing only a single case), this group was regarded as an outlier and excluded from further analyses.

Fig. 1.

Selection of the optimal cluster number and consensus clustering results for DNA methylation. (A) Delta area plot showing the relative change in the area under the cumulative distribution function (CDF) curve for each cluster number k. (B) Consensus matrix heatmap at k = 5. (C) Heatmap of DNA methylation clusters annotated with histological type, Tumor-Node-Metastasis (TNM) stage, clinical stage, race, and age.

To further explore prognostic heterogeneity, we conducted a comparative survival analysis among the four remaining clusters. Kaplan-Meier curves revealed significant differences in overall survival among the clusters (Fig. 2A). Pairwise log-rank tests indicated that Cluster 1 had significantly better survival than Cluster 2 (p = 0.027) and Cluster 3 (p = 0.045). Cluster 4 also exhibited significantly better survival than Cluster 2 (p = 0.027) and Cluster 3 (p = 0.047). The remaining pairwise comparisons were not statistically significant.

Fig. 2.

Survival and clinical/immune characteristics of DNA methylation-defined clusters. (A) Kaplan-Meier survival curves for the four clusters in the training cohort, with statistical significance assessed using the log-rank test. (B) Distribution of lymphocyte infiltration scores across clusters. (C–E) Distribution of clinical stage, metastasis status, and histologic type within clusters. (F–H) Alternative views of panels (C–E) to facilitate cluster-wise comparisons.

Differences in clinical stage, metastases status, and histologic type among cervical cancer patients are associated with substantial variations in incidence, survival, and treatment response. Identical DNA methylation subgroups included patients at different clinical stages (Fig. 2C), suggesting that methylation patterns may be independent of disease stage. Similar patterns have been reported in colon adenocarcinoma, where methylation-based molecular subtypes did not fully align with clinical stages, supporting the view that methylation profiling can offer stage-independent insights into tumor heterogeneity [20]. The same was true for metastatic status, histologic type (Fig. 2D,E). Thus, DNA methylation profiling may serve as a reliable biomarker for cervical cancer, independent of clinical stage. Patients within the same clinical stage exhibited distinct DNA methylation profiles (Fig. 2F). For example, the four methylation-based clusters showed heterogeneous methylation profiles despite sharing the same stage. Similar patterns were also observed for metastasis status and histological subtypes (Fig. 2G,H). These findings suggest that DNA methylation status offers a more refined molecular classification of cervical cancer and provides deeper insight into its underlying heterogeneity.

Lymphocyte infiltration into the tumor microenvironment is generally regarded as a favorable prognostic indicator in cancer patients [21]. CTL expression profiles are often used as indirect markers of lymphocytic infiltration, with higher expression levels linked to improved survival [21]. In this study, lymphocyte infiltration scores were calculated for each sample in the training dataset. Cluster 1 showed the highest immune infiltration among all subgroups (Fig. 2B). This finding suggests a relatively better prognosis for this cluster. In the analysis of metastasis among the four clusters, Cluster 1 had the fewest metastatic samples (Fig. 2D). This result was consistent with its higher lymphocyte infiltration level.

To characterize hypermethylated and hypomethylated CpG loci within specific DNA methylation subgroups of cervical cancer, we employed QDMR software, a quantitative analysis tool. A total of 1842 CpG sites from the 4 subgroups were analyzed to identify subgroup-specific CpGs. For each subgroup, the average DNA methylation level across the 1842 CpG sites was calculated for each sample, generating a dataset with dimensions of 1842 $\times$ 4, which served as input for QDMR analysis. To identify CpG sites specifically associated with each subgroup, we reduced the standard deviation (SD) threshold to 0.01. This adjustment yielded 501 subgroup-specific hypermethylated or hypomethylated CpG sites, corresponding to 443 unique genes. These CpG sites may serve as potential DNA methylation biomarkers for distinguishing molecular subtypes of cervical cancer. The number of subgroup-specific CpG sites varied across subgroups, ranging from 9 to 357 (Fig. 3A), with Cluster 4 showing the highest number of hypermethylated sites (Table 1). The methylation patterns of these CpG sites across the four clusters are presented in Fig. 3B.

Fig. 3.

Specific hyper- and hypomethylated CpG sites in DNA methylation clusters. (A) Selected CpG sites for each prognostic DNA methylation subtype. (B) Heatmap showing the methylation patterns of selected CpG sites across the four identified DNA methylation clusters.

Notably, several well-established biomarkers in cervical cancer research were among the differentially methylated genes identified by QDMR. For instance, promoter hypermethylation of paired box 1 (PAX1) and SRY-box transcription factor 1 (SOX1) has been reported to have high sensitivity and specificity for detecting CIN3⁺ lesions and cervical cancer [22]. The presence of these widely recognized methylation markers in our analysis supports the biological relevance of our subtype classification and highlights its potential utility in cervical cancer diagnosis and prognosis.

We performed enrichment analysis on genes associated with the specific CpG sites identified in each DNA methylation cluster. The analysis revealed that genes unique to Cluster 1 were predominantly involved in toxin metabolism and RNA catabolism (Supplementary Fig. 1a). In Cluster 2, genes were primarily associated with cellular movement and transport (Supplementary Fig. 1b). Cluster 3 was characterized by genes involved in diverse metabolic processes (Supplementary Fig. 1c). Genes in Cluster 4 were largely associated with muscle-related processes and with the development and regulation of the nervous system (Supplementary Fig. 1d). These results indicate that distinct DNA methylation subgroups correspond to different biological functions.

To assess the discriminatory power of the identified CpG sites for each DNA methylation cluster, we constructed a Bayesian model using the training set and evaluated its performance via 10-fold cross-validation, with 501 specific CpGs as features. As shown in Fig. 4A, the area under the ROC curve (AUC) reached 0.824, demonstrating good discriminative ability. The classification performance of the Bayesian model is summarized in the confusion matrix in Table 2, where columns denote predicted cluster labels and rows represent actual cluster assignments. Bold values indicate the number of instances in which predicted clusters matched the actual clusters.

Fig. 4.

Prognostic model performance and predictive outcomes. (A) Receiver operating characteristic (ROC) curve showing the model’s discriminative ability (area under the ROC curve (AUC) = 0.824). (B) Kaplan-Meier survival curves for the four clusters predicted by the model in the test cohort. (C) Lymphocyte infiltration scores across the four clusters in the test cohort. (D–F) Clinical stage, metastasis, and histological type distribution in each DNA methylation cluster. (G–I) Alternative views of panels (D–F) to facilitate cluster-wise comparisons.

We then applied this predictive model to classify samples in the test cohort. Each sample was assigned to a cluster corresponding to the subgroups defined in the training set. Survival analysis showed that the predicted clusters had distinct and statistically significant differences in overall survival (Fig. 4B). These findings indicate that the specific CpG sites identified in this study could serve as potential biomarkers for cervical cancer. Clinical characteristics in the test set were assessed using the same approach as in the training set. Lymphocytic infiltration was also evaluated using the same methodology, revealing highly similar patterns between the training and test cohorts (Fig. 4C). As in the training set, Cluster 1 was the predominant subgroup among both early-stage (I/II) and advanced-stage (III/IV) patients. Among patients with metastasis, Cluster 3 was the most frequent subgroup. In contrast, Cluster 1 consistently predominated among patients without metastasis (Fig. 4D–I). For histological types, all four clusters included a high proportion of squamous cell carcinoma. These findings further demonstrate the reliable predictive performance of the model and the robustness of its identified features.

The inherent limitations of cervical cancer survival data from TCGA may reduce the generalizability of our findings. Therefore, an externally validated dataset comprising 300 cervical cancer samples from GEO (accession number GSE44001) was introduced. Differentially methylated sites identified by QDMR software were first mapped to genes. The mapped genes were then intersected with genes in the GEO dataset, and the overlapping genes were used for further analysis. The external validation cohort was classified using the clustering strategy applied to the training dataset. After cluster assignment, survival analysis was performed for the four identified groups to evaluate potential prognostic differences. Survival analysis revealed significant differences among the four clusters (p = 0.035, Supplementary Fig. 2).

These results indicate that genes mapped from methylation sites can effectively distinguish molecular subtypes with different prognoses at the transcriptome level. This finding demonstrates that the typing system developed in this study has cross-histological stability. It also suggests a potential functional association between DNA methylation patterns and gene expression, which may collectively influence tumor prognosis.

The tumor microenvironment (TME) influences immune responses and affects outcomes of cancer therapy. Next, we examined differences in immune cell infiltration among four clusters. Immune cell abundance in the clusters was assesses using the CYBERSORT algorithm, revealing relatively high proportions of B cells, macrophages, and dendritic cells in each cluster (Fig. 5A). These findings suggested that B cells, macrophages, and dendritic cells could serve as potential therapeutic targets for cervical cancer.

Fig. 5.

Immune characteristics of cervical cancer subtypes. (A) A stacked bar chart displayed the relative proportions of tumor-infiltrating immune cell populations within each cluster. (B) Infiltration levels of 28 immune cell subtypes across the four clusters were quantified through ssGSEA. (C) Immune score, stromal score, estimate score, and tumor purity of the 4 clusters were assessed by the ESTIMATE algorithm. (D) Differential expression patterns of immune checkpoint molecules Programmed Death 1 (PD-1), Programmed Death Ligand 1 (PD-L1), Programmed Death Ligand 2 (PD-L2), and Cytotoxic T-Lymphocyte-Associated Protein 4 (CTLA-4) were compared among the four clusters. Statistical significance was indicated by asterisks: *p 0.05, **p 0.01, ***p 0.001, ****p 0.0001.

The infiltration profiles of 28 immune cell types were evaluated across clusters using the ssGSEA algorithm (Fig. 5B). Cluster 1 exhibited increased infiltration of central memory CD8+ T cells and effector memory CD4+ T cells. This pattern indicated potentially enhanced responsiveness to immune checkpoint inhibitor therapy in these patients. Cluster 2 showed higher levels of central memory CD4+ T cells, eosinophils, and natural killer (NK) cells, along with reduced infiltration of effector memory CD8+ T cells and dendritic cells. These results suggested a distinct immune landscape. Cluster 3 was enriched in $\gamma$$\delta$ T cells and natural killer T (NKT) cells. Cluster 4 showed increased infiltration of CD56+ NK cells. However, Cluster 4 showed relatively low levels of helper T cells, macrophages, and myeloid-derived suppressor cells (MDSCs). This pattern could impair antigen presentation and suppress immune activation, potentially attenuating the antitumor immune response. Furthermore, the ESTIMATE algorithm was used to calculate immune scores, stromal scores, and tumor purity for each cluster (Fig. 5C). Cluster 1 had the highest stromal score. Cluster 3 showed the greatest immune score, and Cluster 4 displayed the highest tumor purity. Collectively, these findings suggested that patients in Clusters 1 and 3 were more likely to benefit from immune checkpoint blockade therapies.

In recent years, immune checkpoint inhibition has become a standard strategy in many cancer therapies. We examined the expression levels of several immune checkpoint genes, including PDCD1 (PD-1), CD274 (PD-L1), PDCD1LG2 (PD-L2), and CD152 (CTLA-4). Cluster 1 showed markedly higher expression of immune checkpoint genes, including PD-L1 and CTLA-4. In contrast, Cluster 4 exhibited the lowest expression of PD-1, PD-L1, PD-L2, and CTLA-4 (Fig. 5D). The pronounced upregulation of these immune checkpoint genes in Cluster 1 indicates that patients in this subgroup could potentially benefit from immune checkpoint blockade therapy.

Differentially methylated sites identified by QDMR software were first mapped to genes. The resulting gene set was then intersected with 2013 immune-associated genes obtained from the ImmPort database, yielding 23 differentially methylated immune-related genes. The DNA methylation profiles of these differential genes were compared across the four clusters. The analysis revealed that the DNA methylation levels of most immune genes differed among the clusters (Fig. 6). Notably, most differential immune genes in Cluster 4 exhibited high DNA methylation levels. This pattern could suppress the expression of these genes, potentially impairing immune cell function and activity, and may contribute to tumor immune escape. In contrast, most specific immune genes in Cluster 3 displayed relatively low DNA methylation levels, suggesting increased gene expression activity. This may support the maintenance of normal immune function. These findings provide a more detailed understanding of how epigenetic alterations influence immune dynamics in the tumor microenvironment and may guide future research on tumor immunity.

Fig. 6.

Comparison of immune gene–specific DNA methylation profiles among the four identified clusters. Statistical significance was indicated by asterisks and “ns”: ***p 0.001, ****p 0.0001, ns = not significant.

Chemotherapy remains a cornerstone in the treatment of cervical cancer. To assess potential differences in drug response among molecular subgroups, we used the ‘oncoPredict’ R package to estimate the sensitivity of the four DNA methylation clusters to three commonly used chemotherapeutic agents. As shown in Fig. 7A, Cluster 4 exhibited relatively high predicted drug sensitivity values for three agents: cisplatin (p = 0.0132), paclitaxel (p = 0.0144), and docetaxel (p = 0.0135). In general, higher predicted sensitivity values indicate better responsiveness to chemotherapeutic agents, suggesting that patients in Cluster 4 may respond more favorably to these treatments. In terms of IC50 levels (Fig. 7B), docetaxel showed the lowest IC50 values in all clusters. A lower IC50 value indicates that a lower drug concentration is required to achieve 50% inhibition of tumor cell growth, reflecting better chemosensitivity. Accordingly, docetaxel showed the greatest predicted chemosensitivity among the three drugs.

Fig. 7.

Drug sensitivity across clusters. (A) Box plot illustrating differences in predicted drug sensitivity to three chemotherapeutic agents (cisplatin, docetaxel, and paclitaxel) among clusters. (B) Half-maximal inhibitory concentration (IC50) values of three chemotherapeutic drugs in four clusters are shown. **** denotes statistically significant differences (p 0.0001).

Given that variations in the immune microenvironment across clusters could be reflected in pathology images, we used “CellProfiler” software to extract features from preprocessed cervical cancer pathology section images. A total of 385 image-based metrics were obtained for further analysis. The Kruskal-Wallis test identified significant differences in several image features among the four clusters. Specifically, in the four features of nucleus area, length of the transverse axis of the nucleus, maximum intensity of hematoxylin staining in the nucleus, and the proportion of staining material in the cytoplasm, Cluster 4 exhibited higher numerical features compared with the other clusters (Fig. 8A). These findings indicated that cells in Cluster 4 exhibited heterogeneity and increased proliferative activity, suggesting that the tumor tissue cells corresponding to this cluster had a higher degree of malignancy, and their biological behaviors might be more invasive.

Fig. 8.

Differences in pathological features between subtypes. (A) Comparisons among the four clusters in nucleus area, length of the transverse axis of the nucleus, maximum intensity of hematoxylin staining in the nucleus, and staining material in the cytoplasm. (B) Variations in the types of tumor-infiltrating lymphocytes (TILs) among the four clusters.

In addition, the distribution of TILs in pathology images was used to assess lymphocyte localization and the strength of immune response. Fig. 8B shows variations in the composition of TIL types among clusters. In all four clusters, TILs were predominantly of the “Brisk Diffuse” type, consistent with previous descriptions in cervical cancer [19], which supports the reliability of our results. Taken together, these observations indicate that the molecular differences among subpopulations were also reflected in the spatial organization and morphological characteristics of cells. This provides morphological evidence for understanding the heterogeneity of the immune microenvironment in cervical cancer and its potential role in tumor progression.