Development and Validation of an Anoikis-Related Gene Signature for Prognostic Prediction in Cervical Cancer

Introduction

Cervical cancer plays a significant role in global morbidity and mortality rates among women, with 604,127 new cases and 341,831 deaths in 2020 globally.¹ Current screening methods, such as cytology, human papillomavirus testing, and colposcopy, have limitations in their specificity and sensitivity, leading to a high rate of false-negative results.² Therefore, novel biomarkers with improved accuracy and reliability are urgently needed to enhance the early detection, personalized treatment, and prevention of misdiagnosis and underdiagnosis of cervical cancer.³ This research area is of great importance as it aims to improve the management and outcomes of cervical cancer patients.

Anoikis is a specialized form of programmed cell death that occurs when cells lose their attachment to the extracellular matrix (ECM) or neighboring cells.⁴ This mechanism plays a critical role in maintaining cellular homeostasis and facilitating proper tissue development. In typical physiological conditions, cells rely on these connections with the ECM and adjacent cells to regulate their proliferation and differentiation. However, when these connections are compromised or disrupted, anoikis is initiated to ensure the preservation of normal physiological states.⁵ Anoikis is primarily regulated through the interplay of two apoptotic pathways, involving the activation of cell surface death receptors and interference with mitochondrial function.⁶ Notably, in certain pathological conditions, such as cancer, cells can develop resistance to anoikis. Extensive research has shown that cancer cells’ resistance to anoikis enables them to detach from their original anchorage-dependent growth, migrate through the lymphatic and circulatory systems to distant sites, and continue to grow and proliferate.^7–9 Consequently, Anoikis resistance is considered a critical initial step in tumor invasion and metastasis. Furthermore, this resistance can reconfigure the tumor microenvironment, facilitating immune evasion and rendering the tumor cells resistant to chemotherapy.^9–11 Nonetheless, the specific mechanisms underlying the role of anoikis-related genes (ARGs) and pathways in cervical cancer remain poorly understood. Therefore, there is a significant need to identify cervical cancer-specific anoikis-related biomarkers, investigate their underlying mechanisms, and explore their implications for prognosis. This endeavor has substantial clinical relevance and with the potential to significantly enhance the early diagnosis and personalized treatment of cervical cancer.

The study aimed to identify and validate key ARGs that are associated with prognosis in cervical cancer based on public datasets from The Cancer Genome Atlas (TCGA) and the Cancer Genome Characterization Initiative (CGCI). We also developed a risk score model to facilitate individualized prediction of survival outcomes. The main objective was to identify novel biomarkers and methodological foundations for personalized screening, targeted therapy, and prognostic assessment for cervical cancer.

Materials and Methods

Data Acquisition and Organization

We obtained the gene expression data (RNA-seq) for TCGA cervical squamous cell carcinoma and endocervical adenocarcinoma (TCGA-CESC) and CGCI HIV+ Tumor Molecular Characterization Project (CGCI-HTMCP-CC) from the GDC database maintained by the National Cancer Institute, with pertinent clinical and survival data (https://portal.gdc.cancer.gov/repository as of October 14, 2023). Table 1 provides comprehensive clinical details for bothcohorts. To harmonize the expression data from these two datasets, we applied the R package ‘sva’ to mitigate batch effects inherent in different cohorts. Furthermore, we retrieved 700 ARGs from the Genecards (https://www.genecards.org/, as of October 2, 2023) and the Harmonizome portal (https://maayanlab.cloud/Harmonizome/, as of October 2, 2023). The whole-exome sequencing (WES) and copy number data were also downloaded from the TCGA database. The expression array datasets (GSE9750 and GSE52903) were downloaded from Gene Expression Omnibus (GEO) database. The immunohistochemistry result was retrieved from Human Protein Altas (HPA) database.

Table 1 Clinical Characteristics of TCGA-CESC and CGCI-HTMCP-CC Datasets

Identification of Anoikis-Related Patterns

In this study, we initially screened 503 candidate ARGs using univariate Cox regression, a standard method for preliminary prognostic biomarker identification.^12,13 Although univariate Cox analysis cannot account for gene interactions, it is an efficient initial screening tool to identify candidate genes significantly associated with survival (p < 0.05), reducing dimensionality for subsequent modeling. The primary purpose was to establish a prognostically relevant gene set as input for LASSO regression analysis and further exploration. This initial screening identified 83 prognostically relevant ARGs from the RNA-Seq data. We used the expression profiles of these ARGs to investigate potential anoikis-related molecular subtypes using the R package ‘ConsensusClusterPlus’. The optimal number of clusters was determined using a consistency matrix. Principal Component Analysis (PCA) was employed to assess subtype stability, and survival differences were evaluated using Kaplan-Meier analysis. An expression heatmap was also generated to illustrate differences between the subtypes.

Construction and Evaluation of the Anoikis-Related Prognostic Signature

To establish and appraise prognostic markers associated with anoikis, we randomly allocated 385 patients into a training set (n=231) and a validation set (n=154), employing the R package ‘caret’. In the training set, the Least Absolute Shrinkage and Selection Operator (LASSO) was employed on the 83 prognostically relevant ARGs to select key predictors, balancing model robustness and overfitting control. The ‘glmnet’ package was used, and the optimal regularization parameter (lambda) was determined via 10-fold cross-validation based on the minimum partial likelihood deviance. This process identified 10 ARGs with non-zero coefficients. The risk score calculation formula is delineated as follows:

The specific formula derived was: risk score = (0.098 * SPP1) + (0.116 * CXCL8) + (−0.239 * TP73) + (−0.280 * GZMB) + (0.310 * ITGA5) + (0.314 * ENDOG) + (−0.356 * CSK) + (−0.377 * MAP2K2) + (−0.397 * PRKACA) + (0.538 * KDM3A). Patients in both the training and validation sets were categorized into high-risk and low-risk groups based on the median risk score calculated in the training set. The prognostic accuracy of the risk score was evaluated using Kaplan-Meier survival curves (Log rank test) and time-dependent receiver-operating characteristic (ROC) curve analysis at 3 and 5 years in both sets.

Comparative Analysis of Prognostic Models and Selection of the Final Model

To rigorously validate the reliability of the model constructed using our selected methodology, we conducted a comparative analysis. This analysis utilized a 10-fold cross-validation strategy, repeated 10 times, on the combined TCGA-CESC and CGCI-HTMCP-CC dataset (n=385 after filtering). Crucially, all models compared were trained and evaluated using the expression data of the same 10 ARGs (SPP1, CXCL8, TP73, GZMB, ITGA5, ENDOG, CSK, MAP2K2, PRKACA, KDM3A) identified through our initial feature selection process described before. This ensured a direct comparison of the algorithms’ performance based on the final 10 gene signature. The following models were evaluated:

Cox Proportional Hazards Model: Trained using the ‘survival::coxph’ function in R. This standard survival analysis method estimates the hazard ratio assuming proportional hazards over time.¹⁴

Random Survival Forest (RSF): An ensemble method based on survival trees, implemented using the ‘randomForestSRC’ package.¹⁵ We used 100 trees (ntree=100) per forest. RSF can capture non-linear effects and interactions without the proportional hazards assumption.

XGBoost Survival: A gradient boosting machine algorithm adapted for survival analysis using the Cox partial likelihood objective function, implemented with the ‘xgboost’ package.¹⁶ We used 100 boosting rounds (nrounds=100).

Gradient Boosting Machine (GBM) Cox: Another gradient boosting approach using the Cox partial likelihood, implemented with the ‘gbm’ package.¹⁷ We used 1000 trees (n.trees=1000) and an interaction depth of 3.

Model performance was assessed in each fold of the cross-validation using time-dependent Area Under the Curve (tdAUC) at 1, 3, and 5 years (calculated using the timeROC package) on both the training and testing sets. The average performance across all repeats and folds was calculated.

Based on this comparative analysis, considering model robustness (similarity between training and test performance), predictive accuracy, and interpretability, the LASSO-Cox approach (used for initial feature selection) leading to the 10-gene Cox model was chosen as the final prognostic model for subsequent analyses and nomogram construction.

Immunohistochemistry

To evaluate the prognostic function of Granzyme B (GZMB), we used the tissue microarray (TMA) obtained from Shanghai Outdo Biotech Co., Ltd. (Shanghai, China). It consisted of 126 cervical cancer tissues and 42 pair-matched normal samples, and the corresponding clinical information is available online. The TMA slide was dried at 60°C for 1 hour, dewaxed in xylene, and dehydrated in a gradient ethanol series. Antigen retrieval was performed by heating the slide in a microwave oven in a container filled with EDTA antigen retrieval buffer (pH 9.0). Subsequently, the slide was immersed in 3% hydrogen peroxide for 30 minutes to inhibit the activity of endogenous peroxides. Next, the TMA tissue was incubated with 10% goat serum (Boster Biological Technology Co., Ltd., Wuhan, China) and blocked at room temperature for 1 hour. Then, the TMA slide was incubated with anti-GZMB (1: 100 dilution, 13588‐1‐AP, Proteintech Group, Wuhan, China) for 2 hours at room temperature. The tissue was incubated at room temperature for 50 minutes with horseradish peroxidase (HRP)-labeled goat anti-rabbit secondary antibody (1:200 dilution, GB23303, ServiceBio, Wuhan, China). Then, the tissue was stained with diaminobenzidine (DAB) for 90 seconds. Finally, the slide was counterstained with hematoxylin, dehydrated, and mounted. To evaluate IHC staining, semi-quantitative scoring criteria were used. The stained slide was scored by 2 observers using a semi-quantitative hybrid (H) scoring method: H-score = (% cells of 0 intensity) + (% cells of 1+ intensity × 1) + (% cells of 2 + intensity × 2) + (% cells of 3 + intensity × 3) + (% cell of 4 + intensity × 4). The staining intensity was classified as follows: 0 (negative), 1 (weak), 2 (moderate), or 3 (strong). The staining area was defined as the distribution of positively stained tumor cells.

Clinical Factor Stratification Analysis and Development and Validation of Nomogram

We conducted an analysis to assess the correlation between the risk score and various clinical characteristics, including patient age, grade, TNM and FIGO stage. For analyses involving clinical variables (age, grade, T stage, FIGO stage), patients with missing data in any of these variables were excluded from the corresponding analysis. Additionally, we performed stratified survival analysis based on tumor stage (T and FIGO) to evaluate differences in survival between high and low-risk groups. To enhance the clinical applicability of the risk score, we integrated it with clinical features to develop a nomogram for predicting 3- and 5-year survival in patients with cervical cancer. The accuracy of the nomogram was assessed through ROC curves, C-index, calibration curves, and DCA curves.

Identification of Differentially Expressed Genes and Functional Enrichment Analysis

Differentially expressed genes (DEGs) between high- and low-risk subgroups were identified using the ‘limma’ R package, with thresholds of adjusted p-value<0.05 and |logFC| > 1. We performed GO, KEGG, and GSEA analyses using the R package ‘clusterProfiler’ to elucidate the functions and pathways associated with DEGs concerning cervical cancer prognosis.

Tumor Immune Microenvironment Assessment in Risk Subgroups

Using the ESTIMATE algorithm, we calculated immune, stromal, and tumor purity scores for each patient. Single-sample GSEA (ssGSEA) was performed on immune gene sets with the ‘GSVA’ R package. The CIBERSORT algorithm was employed to determine the differential abundance of 22 tumor-infiltrating immune cells (TICs) in both high and low-risk subgroups.

Statistical Analysis

All statistical analyses were performed using R software (version 4.3.1). We conducted univariate and multivariate Cox regression analyses using the ‘survminer’ package and assessed the survival curves via the Kaplan-Meier Log rank test. The nomogram was created using the ‘rms’ package. Decision curve analysis was performed using the ‘ggDCA’ package. For data visualization, we used tools such as ‘ggplot2’, ‘ggpubr’, and ‘pheatmap’. Group comparisons were made using t-tests or Mann–Whitney U-tests, and differences within groups were assessed through one-way ANOVA. Statistically significant was defined as a p-value < 0.05 (*p < 0.05; **p < 0.01; ***p < 0.001).

Results

Identification of Anoikis-Related Genes with Prognostic Value

The study’s workflow is illustrated in Supplementary Figure 1. Initially, we obtained data from the TCGA-CESC and CGCI-HTMCP-CC cohorts from the GDC database, encompassing a total of 27,733 genes. To investigate the connection between ARGs and survival, we further analyzed these genes by intersecting them with the acquired ARGs, ultimately isolating 503 ARGs with expression data. Through univariate Cox regression analysis, we pinpointed 83 statistically significant (p <0.05) ARGs that exerted a substantial impact on survival prognosis (Supplementary Table 1). To explore the biological functions represented by the prognostically significant ARGs, we performed GO and KEGG functional enrichment analyses on the 83 genes significantly associated with prognosis in univariate Cox regression analysis. Results showed that these genes were significantly enriched in biological processes and signaling pathways closely related to tumor progression, including apoptosis regulation, extracellular matrix remodeling, immune response, and angiogenesis (Supplementary Tables 2 and 3). These functional enrichment results support the biological rationale of our prognostic model based on ARGs.

Discovery of Anoikis Subtypes in Cervical Cancer

To investigate whether anoikis-related patterns could be employed for molecular subtyping of cervical cancer patients, we conducted unsupervised cluster analysis based on the expression profiles of the 83 prognostically relevant ARGs. Our findings suggested that k = 2 was the optimal number of clusters (Figure 1A), indicating the categorization of cervical cancer patients into two distinct anoikis-related subtypes characterized by ARGs, namely subtype A (n = 194) and subtype B (n = 233). Subsequent PCA analysis revealed clear differentiation in the distribution of patients between these two subtypes (Figure 1B). Additionally, we generated expression heatmaps to depict the ARG expression in these subtypes (Figure 1C). Kaplan-Meier survival curves demonstrated a significantly lower survival rate (p<0.001) for subtype A compared to subtype B (Figure 1D).

Figure 1 Identification and prognostic significance of anoikis-related subtypes. (A) Consensus clustering analysis identified two clusters (k = 2). (B) PCA analysis of distinct anoikis subtypes. (C) Heatmap of the ARGs expression in the two subtypes. (D) Survival analysis of distinct anoikis subtypes.

Prognostic Risk Model Construction and Validation

In order to pinpoint the most robust prognostic predictive features, we executed LASSO and multivariate Cox analyses using anoikis gene expression data from the training set (n=232) for screening prognostically relevant ARGs (Figure 2A and B). Consequently, we identified 10 genes (SPP1, CXCL8, TP73, GZMB, ITGA5, ENDOG, CSK, MAP2K2, PRKACA, KDM3A) for the development of prognostic signature (Figure 2C). The risk score was computed using the formula: risk score = (0.098* SPP1) + (0.116* CXCL8) + (−0.239* TP73) + (−0.280* GZMB) + (0.310* ITGA5) + (0.314* ENDOG) + (−0.356* CSK) + (−0.377* MAP2K2) + (−0.397* PRKACA) + (0.538* KDM3A). We further incorporated WES data from the TCGA database and analyzed mutations and copy number alterations of these genes. Among 304 cervical cancer samples, 25 (8.2%) carried mutations of ARGs, with the highest mutational frequency observed in KDM3A (3.3%) (Supplementary Figure 2A). Additionally, exploring CNV alteration frequency determined a higher incidence of copy number amplification in the CSK, CXCL8, and PRKACA, whereas MAP2K2 and TP73 exhibited a greater frequency of copy number loss (Supplementary Figure 2B).

Figure 2 The establishment of ARGs risk signature. (A and B) LASSO analysis to screen candidate genes. (C) Coefficient of selected 10 ARGs for risk signature. Distribution of risk scores and patient status in (D) training cohort and (E) validation cohort. Survival curves for the (F) training cohort and (G) validation cohort. Time dependent ROC curve of the risk model in (H) training cohort and (I) validation cohort.

Subsequently, we categorized patients into a high-risk group (n=115) and a low-risk group (n=116) based on the median risk score from the training set (Figure 2D and Supplementary Table 4). The results from the survival analysis demonstrated that the high-risk group had significantly shorter survival times compared to the low-risk group (p<0.001). The Kaplan-Meier curves confirmed these results (Figure 2F). ROC curve analysis revealed that the risk signature of anoikis exhibited AUC values of 0.822 and 0.834 at 3 and 5 years, respectively (Figure 2H).

In the validation set (n=154), the same risk score formula was employed to calculate risk scores, and patients were categorized into high-risk (65 cases) and low-risk (89 cases) groups based on the median risk scores from the training set (Figure 2E and Supplementary Table 4). Consistent with the training set results, patients in the high-risk group exhibited significantly worse prognoses than those in the low-risk group (p<0.001) (Figure 2G). ROC curve analysis demonstrated that the risk profile’s performance in the validation set was also excellent, with AUC values of 0.722 and 0.815 at 3 and 5 years, respectively (Figure 2I), establishing it as a valuable prognostic indicator.

Immunohistochemical staining was performed to experimentally verify the expression of GZMB, and representative staining fields were shown in Figure 3A–D. The results showed that GZMB was significantly overexpressed in cervical cancer compared to normal tissue (Figure 3E). High GZMB expression was defined as an H-score greater than the median. We conducted stratified survival analysis of GZMB expression and demonstrated the protective prognostic value of GZMB, as the samples with high GZMB expression had significantly better survival outcomes than the samples with low GZMB expression (Figure 3F–G). Additionally, we preliminarily validated the expression of these genes in cervical cancer using two expression profile chips (GES9750 and GSE52903) and the HPA database. The results showed that these genes were either upregulated or exhibited an upward trend in cervical cancer (Supplementary Figure 3).

Figure 3 Experimental validation of immunohistochemical staining of GZMB in TMA cohort and its prognostic value. (A–D) The intensities of GZMB immunostaining were (A) negative, (B) weakly positive, (C) moderately positive and (D) strong positive in cervical cancer tissues. Left image, magnification ×50; right image, magnification ×400; the red squares indicate the area shown in higher magnification. (E) The H-score of GZMB in cervical cancer and adjacent normal tissues. (F and G) Kaplan–Meier survival curves to compare (F) OS and (G) DFS between the high GZMB expression and low GZMB expression patients.

Comparative Analysis of Prognostic Models

To evaluate the performance of the identified 10-gene signature using different modeling algorithms, we compared its performance against several machine learning algorithms (RSF, XGBoost, GBM) using repeated 10-fold cross-validation, with all models trained on the same 10-gene signature. The average performance metrics are summarized in Supplementary Table 5. While some machine learning models like XGBoost showed high AUC values on the training data (eg, 5-year AUC = 1.00), they exhibited a noticeable drop in performance on the validation data (eg, 5-year Test AUC ≈ 0.66), suggesting potential overfitting within this dataset. The Cox model demonstrated more stable performance across training and validation sets (eg, 5-year Train AUC ≈ 0.824, 5-year Test AUC ≈ 0.803), indicating better generalization. These results support the selection of the 10-gene Cox model for its balance of predictive accuracy, robustness, and interpretability in this context.

Stratified Analysis of Clinical Factors

We conducted a correlation analysis between risk scores and clinical characteristics (Figures 4A–F). Notably, we found that none of the differences in risk scores were statistically significant for patients of different ages, grade, N and M stage. However, patients with advanced T stage (T3-T4) and FIGO stage exhibited higher risk scores, indicating that our risk score model might encompass key information related to the biological behavior of cervical cancer. Subsequently, we performed stratified risk analysis, and the Kaplan-Meier results demonstrated that the risk signature exhibited excellent predictive capability in advanced T stage (T3-T4) and FIGO stages (p < 0.001), underscoring its reliability and clinical applicability as an independent risk factor (Figure 4G–J).

Figure 4 Relationship between ARGs risk scores and clinical characteristics. (A–F) Differences in risk scores of patients with regard to age, FIGO stage, grade, T, N and M stage. (G and H) Survival analysis of patients with OS regrouped according to T stage. (I and J) Survival analysis of patients with OS regrouped according to FIGO stage.

Establishment of Risk Nomogram

Both univariate Cox regression analysis and multivariate Cox regression analysis demonstrated that risk score, T stage and FIGO stage were all significantly linked to overall survival (Figure 5A and B). The findings strongly support the independent prognostic value of the risk signature in cervical cancer patients. Subsequently, we generated a nomogram illustrating the relationship between grade, age, T stage, FIGO stage and risk score in relation to the calculation of survival probability (Figure 5C). The ROC curve analysis indicated that the area under the ROC was 0.819 and 0.834 at 3 and 5 years, respectively, underscoring its strong prognostic validity (Figure 5D). Additionally, the C-index of the nomogram was 0.781. Calibration curves demonstrated significant agreement between the actual and predicted survival probabilities at 3 and 5 years (Figure 5E). Moreover, the results from the decision curve analysis at 3 and 5 years highlighted the superior predictive performance of the nomogram compared to other independent clinical factors (Figure 5F and G).

Figure 5 Establishment and evaluation of a predictive nomogram. (A and B) Results of the univariate Cox regression analysis and multivariate Cox regression analysis regarding OS of the risk signature. (C) The nomogram based on ARGs score and clinical characteristics to predict the prognosis at 3 and 5 years. (D) The ROC curves of the nomogram. (E) The calibration curve for evaluating the accuracy of the nomogram model. The grey dashed diagonal line represents the ideal nomogram. (F and G) The DCA curves of the nomogram at 3 and 5 years.

We also compared the performance of our ARG model with several previously published gene prognostic models for cervical cancer. Du et al developed an ferroptosis-based model with an AUC of 0.78 (5-year).¹⁸ Shi et al’s autophagy-based model demonstrated lower predictive ability with test AUCs of 0.603 (5-year).¹⁹ Wang et al’s immunology-based signature reported an AUC of 0.815.²⁰ Our ARG model (AUC = 0.83) demonstrates competitive predictive performance compared to these benchmarks. However, we acknowledge that this comparison is indirect due to differences in patient cohorts, data types, and validation methodologies across studies. Further validation in independent external cohorts is needed to definitively establish the superiority of our model.

Differentially Expressed Genes and Functional Enrichment Analysis

To unveil differences in biological functions, we compared gene expression between patients in the high-risk and low-risk groups, ultimately identifying 193 differentially expressed genes (DEGs), with detailed information available in Supplementary Table 6. The GO enrichment results demonstrated that these DEGs are highly involved in immune system activation, production of immunoglobulins and immune mediators, antigen recognition and presentation, and other immune response processes (Figure 6A and Supplementary Table 7). The enriched KEGG pathways demonstrated that the differentially expressed genes are strongly correlated with activations of immune cell development and differentiation, autoimmunity, inflammatory cytokine production, pathogen recognition, cell adhesion and migration, antigen presentation, and infections (Figure 6B and Supplementary Table 8). The GSEA analysis based on Hallmark gene sets demonstrated increased immune responses similar to transplant rejection, activations of interferon signaling, enhanced epithelial-mesenchymal transition, and upregulated UV DNA damage responses in the samples (Figure 6C and Supplementary Table 9). These processes have established associations with immune inflammation and cancer progression. The results provide insights into the biological mechanisms represented by the differential gene expression profile.

Figure 6 Biological functional and pathway enrichment analysis of high-risk group and low-risk group based on the risk signature. (A) Top 30 most enriched GO terms according to differential genes between high- and low- risk populations. (B) Top 30 most enriched KEGG pathways of the common differently genes between high- and low-risk populations (C) GSEA enrichment analysis based on hallmark gene set showing the activation states of biological pathways in high- and low-risk populations.

Immunosuppressive Microenvironment in the High-Risk Group

To explore the distinctions in the tumor immune microenvironment between the high-risk and low-risk groups, we utilized the ESTIMATE algorithm to assess stromal scores, immune scores, and tumor purity in cervical cancer patients (Figure 7A). The results revealed significant differences in the tumor microenvironment characteristics between the two groups. Specifically, the low-risk group exhibited significantly higher Immune Scores and Stromal Scores (indicating greater non-tumor cell infiltration) compared to the high-risk group, while tumor purity showed the inverse trend. Given the pronounced differences in immune scores, we further evaluated the immune infiltration patterns of 22 immune cell types in cervical cancer patients using the CIBERSORT algorithm (Figure 7B). The high-risk group showed higher proportions of T cells CD4 memory resting, macrophages M0, activated mast cells and eosinophils, whereas lower proportions of B cells naive, plasma cells, CD8 T cells, activated CD4 memory T cells, T cell follicular helpers, regulatory T cells, macrophage M1, resting dendritic cells, and resting mast cells were noted. Additionally, further immune cell enrichment analysis using the ssGSEA algorithm indicated generally lower immune cell infiltration in the high-risk group, aligning with the lower immune scores observed in this group (Figure 7C).

Figure 7 Analysis of tumor immune microenvironment in cervical cancer patients. (A) Stroma, immune, and ESTIMATE scores based on ESTIMATE algorithm. (B) The violin plots for the infiltration of 22 immune cells according to CIBERSORT analysis. (C) The boxplots for the tumor infiltrating lymphocytes assessments by the ssGSEA algorithm. *p < 0.05, **p < 0.01, and ***p < 0.001.

Discussion

Cervical cancer stands as one of the most prevalent malignant tumors affecting the female reproductive system, with a persistently high global incidence and mortality rate.²¹ Early detection and diagnosis are paramount for improving the prognosis of cervical cancer patients. While widely utilized screening methods such as cytology and HPV testing have achieved some degree of success, their suffer from limited specificity and sensitivity, underscoring the need for more effective biomarkers to enhance cervical cancer diagnosis and patient risk stratification.^22,23 In this study, we identified and validated a novel anoikis-related gene signature, which has significant prognostic value in cervical cancer and effectively stratifies patients into distinct risk groups, potentially aiding in personalized risk assessment. These findings suggest that anoikis-related genes could serve as promising prognostic biomarkers and provide valuable insights into the biological mechanisms underlying cervical cancer progression. Anoikis, a critical form of programmed cell death, regulates cell number and tissue structure by inducing cells to detach from extracellular matrix connections, thus preventing abnormal proliferation and inhibiting detached cells from reattaching to unsuitable matrices. It plays a pivotal role in tumorigenesis, progression, and metastasis.²⁴ In recent years, the scientific community has increasingly recognized the significance of anoikis, particularly its links to properties such as anchorage-dependent growth²⁵ and epithelial-mesenchymal transition,²⁶ both critical in tumor progression and cancer cell dissemination. Furthermore, studies in various malignancies, including osteosarcoma, lung adenocarcinoma, renal clear cell carcinoma, colon adenocarcinoma, hepatocellular carcinoma, and bladder cancer, have demonstrated the pivotal roles of anoikis regulators.^27–32 However, the full extent of the mechanisms and prognostic value of ARGs in cervical cancer remains to be elucidated.

Previous studies have indicated that ARGs possess potential predictive value in the prognosis of various cancers, and their expression patterns are closely associated with tumor aggressiveness, metastasis, and patient survival. For instance, a study by Frisch et al (2016) demonstrated that ARG expression patterns could predict disease-free survival and overall survival in breast cancer patients, with patients expressing specific ARG sets exhibiting poorer responses to chemotherapy.²⁶ Tan et al (2021) developed an ARG-based risk scoring model that effectively predicted the prognosis of lung adenocarcinoma patients and was significantly correlated with tumor mutation burden (TMB) and immune infiltration.³³ Huang et al (2022) discovered that ARG expression profiles could serve as independent predictors of survival and metastasis in pancreatic cancer patients, with EMT-related ARGs being particularly important.³⁴

By integrating expression profiles of cervical cancer tissue and survival data from the TCGA-CESC and CGCI-HTMCP-CC datasets, we initially screened for ARGs sourced from the Genecards website and the Harmonizome portal. Through univariate Cox prognostic analysis, we identified a set of prognostic relevant ARGs, based on which we divided samples into subtype A and subtype B. The significantly poorer prognosis of subtype A suggested that ARGs may play a substantial role in cervical cancer prognosis. Using the identified prognostically relevant ARGs, we developed a robust risk score model featuring 10 key genes to predict the prognosis of cervical cancer patients. In selecting our final prognostic model, we rigorously compared the performance of the LASSO-Cox derived 10-gene signature against several machine learning approaches using repeated cross-validation (Supplementary Table 5). While machine learning models hold promise for capturing complex interactions, our analysis indicated they were prone to overfitting on this dataset, showing high performance on training data but diminished accuracy on unseen test data. The Cox model provided a more robust and generalizable performance, while also offering clear interpretability via hazard ratios for each gene. This balance justified its selection as our primary model. Our subsequent analyses demonstrated that this model exhibits excellent stratification and predictive performance, effectively categorizing cervical cancer patients into high-risk and low-risk groups. Patients in the high-risk group exhibited a significantly poorer prognosis. We found a strong association between high-risk scores and T stage and FIGO stage, and the nomogram proved that our risk model had a better performance than other clinical factors alone in prognosis prediction, underscoring itspotential clinical utility in predicting cervical cancer outcomes. Building on these findings, we further demonstrated the predictive power of the ARG-based model in both training and validation cohorts. Specifically, our analysis revealed that patients in the low-risk group had significantly higher survival rates compared to the high-risk group, further confirming the close association between ARG expression patterns and the malignant progression of cervical cancer. Furthermore, GO and KEGG enrichment analyses indicated that these ARGs are primarily involved in crucial biological processes and signaling pathways such as cell apoptosis, cell adhesion, tumor microenvironment remodeling, and immune evasion. The aberrant activation of these pathways may enhance invasiveness and metastatic potential, enabling cancer cell to evade anoikis and drive cervical cancer progression.

Previous studies have shown that ITGA5 is closely associated with the survival and progression of cervical cancer patients and effectively promotes angiogenesis via the AKT/VEGFA pathway in cervical cancer.³⁵ Researchers have also reported that E6/E7 of HPV16 and 18 regulates the overexpression of CXCL8 in cervical cancer tissues, and that CXCL8 is associated with malignant status and prognosis in cervical cancer patients.^36,37 TP73 was overexpressed in cervical cancer and associated with favorable patient prognosis.³⁸ GZMB is involved in programmed cell death and natural killer cell mediated cytotoxicity, and plays a key role in LINC02474 inhibiting apoptosis in colorectal cancer.³⁹ And as expected, GZMB, overexpressed in cervical cancer tissue, was a significantly protective prognosis marker in our experimental validation.

This study employed univariate Cox regression as an initial screening method for ARGs, a commonly used approach in biomarker studies.⁴⁰ Although univariate analysis has inherent limitations, such as its inability to account for confounding factors and interactions between variables, it offers the advantages of being intuitive and computationally efficient as an initial screening tool. When screening for potential prognostic biomarkers from high-dimensional gene data, the strategy of combining univariate Cox regression with subsequent multivariate analysis is widely adopted.⁴¹ Using this analytical strategy, we successfully identified 83 ARGs significantly associated with cervical cancer prognosis, laying the foundation for subsequent model construction. To construct the prognostic model, we systematically evaluated the applicability of various statistical and machine learning methods. Among them, traditional Cox proportional hazards regression, as a standard method for survival analysis, offers direct clinical interpretability, providing hazard ratios (HRs) and confidence intervals (CIs), while effectively handling censored data.⁴² However, the Cox model faces significant limitations when dealing with high-dimensional data and is strictly dependent on the proportional hazards assumption. To overcome these limitations, we employed the LASSO regression method. LASSO (Least Absolute Shrinkage and Selection Operator) effectively enhances model performance with high-dimensional data and in the presence of collinearity by introducing an L1 regularization term to simultaneously achieve feature selection and coefficient estimation.^43,44 We also assessed the potential of various machine learning approaches for prognostic modeling. Random Forest, as a tree-based ensemble method, can effectively capture non-linear relationships and complex interactions, is not restricted by the proportional hazards assumption, and exhibits strong robustness to outliers.¹⁵ Gradient Boosting Machine (GBM) and its optimized version XGBoost optimize predictive performance by sequentially building weak learners, and are particularly suitable for handling mixed-type data.⁴⁵ However, these machine learning methods generally suffer from “black box” characteristics, and limited clinical interpretability, and often require larger sample sizes to avoid overfitting. After comprehensive evaluation, we selected LASSO-Cox regression as the primary modeling method for this study, mainly based on the following considerations: (1) LASSO-Cox can effectively control the risk of overfitting while maintaining a degree of clinical interpretability; (2) Through its built-in feature selection process, LASSO can identify the most prognostically valuable gene biomarkers, providing clues for subsequent biological mechanism research; (3) Under the current sample size conditions, the machine learning models we implemented (Random Forest, XGBoost, GBM) did not significantly outperform the LASSO-Cox model, possibly due to the relatively limited sample size.^46,47 It is worth noting that when the sample size increases and multi-omics data are incorporated, methods such as deep learning may demonstrate greater advantages.⁴⁸

While our risk scoring model demonstrates promising predictive performance, it is currently a research tool and not yet poised for direct clinical application. Further research is needed to translate these findings into clinical utility. Future studies should focus on prospective validation of our gene signature in independent patient cohorts, ideally in a clinical trial setting. Moreover, exploring the development of a clinically applicable diagnostic assay based on our risk scoring model could be a valuable next step. Such an assay could potentially be used to identify high-risk patients who may benefit from more intensive treatment or closer monitoring, and to guide personalized treatment strategies based on individual risk profiles. In addition, functional validation of the identified anoikis-related genes and further investigation of their role in cervical cancer biology are warranted to fully understand their clinical implications.

Our ARG prognostic model demonstrated superior predictive performance compared to both a clinical-factor-only model and previously published multi-gene prognostic models for cervical cancer.^27,32,41 This suggests that anoikis-related genes may capture unique molecular characteristics associated with cervical cancer prognosis. Notably, the nomogram integrating the ARG model with clinical factors further improved prediction accuracy (AUC=0.89), highlighting the value of integrating molecular biomarkers with traditional clinical factors, which aligns with the principles of predictive model development proposed by Steyerberg et al.⁴⁹

Nevertheless, translating biomarkers into clinical applications still faces multiple challenges, including standardization of techniques, cost-effectiveness analysis, and prospective validation.⁵⁰ In the future, we plan to validate the predictive value of the ARG model through multi-center clinical studies and explore its potential applications in treatment decisions and patient stratification, such as identifying high-risk patients who may benefit from intensified treatment or low-risk patients suitable for de-escalated therapy. Furthermore, combining the ARG model with emerging liquid biopsy and circulating tumor DNA technologies may offer new avenues for non-invasive monitoring and personalized treatment of cervical cancer.⁵¹

However, several limitations must be acknowledged. First, our sample size was relatively small, potentially affecting the statistical robustness of the results. Future studies with larger sample sizes are needed to enhance the reliability of the findings. Second, while our study employed an observational design involving survival analysis to evaluate the prognostic value of the markers, experimental studies are needed to validate the functional roles of these markers. Third, this study represents a bioinformatics exploration, and the specific mechanisms require further elucidation. Fourth, it is important to note that this study was retrospective and based on publicly available databases, which could introduce sample selection bias. Future research should aim to integrate data from multiple sources to improve result reliability. Fifth, considering the multifactorial nature of cervical cancer, this study did not comprehensively incorporate genetic and environmental influences, which presents a notable limitation. Finally, the markers identified in this study are still in the early stages of clinical translation, warranting further validation studies with larger sample sizes and the development of precise detection technologies.

In summary, our study systematically constructed and evaluated a risk scoring model for cervical cancer based on 10 ARGs, demonstrating high accuracy in predicting the prognosis of cervical cancer patients within the analyzed datasets. Moreover, we identified potential novel biomarkers and provided methodological foundations that may facilitate future personalized screening, targeted therapy, and prognostic assessment for cervical cancer. We also revealed differences in immune phenotypes and tumor microenvironments between high- and low-risk cervical cancers, which may provide insights for the development of immunotherapies and personalized treatment strategies. Our study enhances the understanding of the role of ARGs in cervical cancer prognosis and provides a foundation for future research aimed at reducing mortality rates and improving patient outcomes.