Gene expression profiles were meticulously selected based on predefined criteria from two distinct sources: 300 samples of CESC tissues were obtained from the Gene Expression Omnibus (GEO: GSE44001), and an additional 304 profiles, inclusive of three samples from normal adjacent tissues, were sourced from The Cancer Genome Atlas (TCGA-CESC). The inclusion criteria for these datasets were as follows: The cohort size needs to be large enough, including patients of different ages and stages, regardless of race, to be included in this study. All data underwent standardized normalization and batch effect removal processes to ensure consistency and reliability in subsequent analyses.

A comprehensive dataset of 514 ARGs was acquired from the GeneCards database (https://www.genecards.org/) and the Harmonizome portal, referenced in publications (13, 14). Following the initial data acquisition, a comprehensive analysis was performed on the TCGA-CESC cohort using the ‘limma’ package (Version: 3.58.1), a component of the Bioconductor software suite (Version: 3.18). This analytical process resulted in the identification of 170 differentially expressed genes (DEGs), achieved through a comparative evaluation of the expression profiles of the 514 ARGs between CESC tumor tissues and adjacent normal tissues.

To elucidate distinct patterns in anoikis regulation, we employed an advanced consensus clustering algorithm, utilizing the k-means method to categorize anoikis-related gene expression profiles. Subsequently, sophisticated techniques for dimensionality reduction were used, specifically t-distributed Stochastic Neighbor Embedding (t-SNE) and Uniform Manifold Approximation and Projection (UMAP). These methods facilitated the validation of the clustering’s robustness and reliability, utilizing the comprehensive capabilities of the R package ‘ggplot2’ for intricate data visualization and analysis.

Gene Set Variation Analysis (GSVA) was chosen for its ability to provide a comprehensive assessment of gene set enrichment variations within CESC tissues, offering insights into the underlying molecular mechanisms. The gene set file ‘c2.cp.kegg.v7.4.symbols.gmt’ required for this analysis was obtained from the Molecular Signatures Database (MSigDB). The enrichment analysis was then meticulously executed using the ‘GSVA’ R package (version 1.50.0). Subsequently, the enrichment analysis was meticulously executed using the ‘GSVA’ R package (version 1.50.0), ensuring a rigorous and precise evaluation of the gene sets (15).

The patient cohort was systematically partitioned into two distinct groups: a training set for the initial construction of the risk model and a validation set dedicated to its subsequent verification. Initial identification of genes associated with survival outcomes commenced with univariate Cox regression analysis. Subsequently, a refined analysis employing a minor absolute shrinkage and selection operator (LASSO) regression technique was performed, utilizing the ‘glmnet’ package (Version: 4.1–8) in R for implementation (16). Determination of the optimal penalty regularization parameter λ was achieved through an exhaustive 10-fold cross-validation methodology. In the ensuing stage, a sophisticated multivariate Cox regression approach was employed to delineate essential genes and ascertain their respective coefficients meticulously. A selection of nine pivotal anoikis-associated gene signatures was made to formulate the risk signatures, guided by the most favorable lambda values and their corresponding coefficients. The computation of each patient’s ARG signature risk score was RiskScore = ∑ 1 9 C o e f i ∗ E x p i .

In this model, ‘Coefi’ denotes the risk coefficient for each gene, while ‘Expi’ represents the gene’s expression level. To ascertain the predictive efficacy of this model, we utilized Kaplan-Meier (KM) survival analysis and time-dependent Receiver Operating Characteristic (ROC) curves. The dual approach facilitated a thorough evaluation of the model’s aptitude in forecasting patient outcomes. Then, a set of nine anoikis-related DEGs exhibiting a significant correlation with overall survival (OS) were discerned through meticulous multivariate Cox regression and LASSO analyses within both the GSE44001 and TCGA-CESC cohorts.

In this research, the CIBERSORT algorithm and single-sample Gene Set Enrichment Analysis (ssGSEA) R scripts were employed for the quantitative assessment of the relative abundances of immune cell subsets infiltrating the tumor microenvironment (17). Specifically, CIBERSORT was utilized to ascertain the composition of various immune cell types within the groups with lower and higher prognostic risk. This method guaranteed that the aggregate score attributed to all inferred immune cell types in each specimen was normalized to one. Additionally, Spearman’s rank correlation methodology was utilized to elucidate the relationships between the computed risk scores and the degrees of immune cell infiltration.

An intricately crafted predictive nomogram integrating clinicopathological features with risk scores has been developed. A calibration plot is employed for internal validation purposes to enhance the model’s accuracy. Additionally, the nomogram’s prognostic capabilities were verified through the application of a time-cumulative index, providing a robust measure of its predictive performance. Furthermore, a decision curve analysis (DCA) was conducted meticulously to evaluate the net clinical benefits of the nomogram (18).

The Tumor Immune Single-Cell Hub (TISCH; http://tisch.comp-genomics.org), a comprehensive online repository specializing in scRNA-sequencing data pertinent to the TME, was utilized as a pivotal resource in this study (19). Subsequently, Hallmark analysis was employed to assess apoptosis level changes within various cell clusters. This extensive database facilitated a systematic exploration of the heterogeneity within the TME across diverse datasets and cellular phenotypes.

To meticulously corroborate the differential expression of the nine selected genes, including BCL2, BAX, IGF1, PLAU, EDA2R, ABL1, MIR200A, FASN, and NTRK3, in control versus cervical cancer tissues, immunohistochemical staining data from the Human Protein Atlas (HPA, http://www.proteinatlas.org) were systematically analyzed (20). This approach enabled a comprehensive visualization of the distinct expression patterns of these genes in tumorous tissue compared to normal tissue. The immunohistochemical analysis provided clear evidence of aberrant gene expression in the cancerous specimens, thereby substantiating the potential role of these genes in cervical carcinogenesis.

The Ect1/e6e7 and HeLa cell lines used in this study were purchased from the Shanghai Cell Bank (Shanghai, China). The cells were cultured in a complete medium (DMEM + 10% fetal bovine serum + 1% penicillin-streptomycin) and maintained in a 37°C incubator with 5% CO₂ and saturated humidity. The medium was refreshed every 2–3 days. Cells were passaged when they reached over 90% confluence.

Cells were first rinsed with PBS after discarding the culture medium to prepare for total RNA extraction. RNA was isolated using the Trizol reagent (AG21102, Precision Biotechnology). The cDNA synthesis was performed with the Evo M-MLV RT Reverse Transcription Kit II (AG11711, Accurate Biotechnology). RT-qPCR was carried out using the SYBR Green Pro Taq HS premixed qPCR kit (AG11701, Accurate Biotechnology) on a LightCycler^® 96 system (Roche Ltd, Switzerland). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) served as the internal control. Gene expression fold changes were calculated using the 2^-ΔΔCT method. Primer sequences are listed in Table 1 .

Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/) has RNA sequencing data of tumor tissues and normal tissues from TCGA and Genotype-Tissue Expression Project (GTEx), through which the platform not only allows access to gene expression in different tumors but also allows single-gene survival analysis (21). TCGA tumors vs TCGA normal + GTEx normal were compared using one-way ANOVA with a p-value cutoff of 0.05.

Comprehensive statistical analyses were conducted utilizing R software (Version: 4.3.2). Specific statistical tests, including t-tests and ANOVA, were applied to meticulously analyze the data. The threshold for statistical significance was established at a P value of less than 0.05, and the false discovery rate (FDR) was rigorously controlled with a cut-off set at q<0.05. The experimental data are presented as mean ± standard deviation.

Initially, 563 ARGs were identified from the GeneCards and Harmonizome portals and narrowed down to 514 ARGs through Venn diagram analysis ( Figure 1A ). Comparison with normal adjacent tissues revealed 170 DEGs in CESC samples. Data from the TCGA-CESC cohort and GSE44001 were integrated, forming a consolidated ‘CESC-GSE44001’ cohort with 17,073 genes. Univariate Cox regression analysis identified 33 out of these 170 DEGs as significantly associated with patient survival (P < 0.05) ( Figure 1B ). Apart from BCL2, IGF1, GLI2, ITGA8, NTRK3, and ONECUT1, which indicated varied prognoses, the remaining 27 genes were predominantly correlated with poor prognosis. A network diagram revealed tight interconnections among these 33 genes ( Figure 1C ). Considering the frequent chromosomal aberrations in CESC, copy number variation data was extracted from the TCGA database ( Figure 1D ). Chromosomal changes in ARGs were thoroughly analyzed to locate each gene accurately. Chromosome one exhibited the highest occurrence of ARGs, including CLIC4, SLC2A1, KIF14, NRAS, IRF6, and PARP1. A significant genomic ‘GAIN’ was observed in COL4A2 on chromosome 13, coupled with a notable ‘LOSS’ of the BRCA2 gene on the same chromosome ( Figure 1E ).

To elucidate the functional implications of ARGs in CESC, we utilized the Consensus Cluster Plus R package (Version: 1.66.0) to perform consensus clustering of 33 prognosis-related DEGs (P <0.05). Optimal stratification into two distinct subtypes was achieved at k = 2 ( Figure 2A ). Subsequent overall survival (OS) analysis (P < 0.001) revealed a statistically notable prognostic disparity between these subtypes ( Figure 2B ). The validity of this bifurcation was corroborated through t-SNE and UMAP analyses, confirming the distinct segregation at k = 2 ( Figures 2C, D ). Heat maps delineating ARGs expression across these subtypes, alongside their clinicopathological attributes, indicated ITGA8 as a potential prognostic enhancer ( Figure 2E ). Furthermore, KEGG pathway enrichment analysis conducted for both subgroups identified significant enrichment of the histidine metabolism pathway in subgroup B, which exhibited a favorable prognosis ( Figure 2F ). Additionally, Gene Set Enrichment Analysis (GSEA) underscored reduced cytokine-cytokine receptor interactions, enrichment of focal adhesion, and JAK-STAT signaling pathways in subgroup B, pivotal pathways implicated in tumor cell migration and colonization of new anchorage sites ( Figure 2F ).

Box plots show the expression patterns of ARGs in the two subgroups ( Figure 3A ). IGF1, GLI2, MIR200A, ITGA8, and NTRK3 levels were lower in subgroup A than in subgroup B. Conversely, the other genes were highly expressed in subgroup A. Given their correlation with overall survival, these DEGs may encode pivotal molecules that significantly influence the prognosis of patients with CESC. Furthermore, they present as potential targets for developing targeted therapeutic strategies. There were also significant differences in the extent of immune cell infiltration, and the proportions of activated CD4⁺ T cells, dendritic cells, and other immune cells in subgroup A were significantly higher than those in subgroup B ( Figure 3B ).

In our investigation into the clinical relevance of ARGs, 33 ARGs (p < 0.05) were subjected to Lasso-penalized Cox regression analysis ( Figures 4A, B ). The optimal lambda (λ) value for Lasso regression, corresponding to the minimum cross-validation error, was determined to be 17. To refine our predictive model, multivariate Cox regression was employed, leading to the selection of 9 pivotal ARGs from the initial 17. These were instrumental in constructing the prognostic model. The resultant risk score, based on this nine-ARG signature, was designated as the “ARG score” with corresponding correlation coefficients detailed in Supplementary Table S1 . The Prognostic Index (PI) was calculated as follows: PI = (0.929 × BAX expression) + (0.177 × PLAU expression) + (0.729 × EDA2R expression) + (0.782 × ABL1 expression) + (0.197 × MIR200A expression) + (0.317 × FASN expression) − (0.616 × BCL2 expression) − (0.362 × IGF1 expression) − (0.978 × NTRK3 expression).

Kaplan-Meier analysis indicated a poorer prognosis for patients in the high-risk group, a trend corroborated in the TCGA-CESC validation cohort ( Figures 4C, D ). Time-dependent ROC curves for the model, assessing 1-, 3-, and 5-year overall survival, demonstrated its robust predictive accuracy ( Figures 4E, F ). A significant disparity in risk scores was observed between the previously identified subtypes ( Figures 4G, H ). An alluvial diagram was constructed to illustrate transitions among ARG clusters, changes in ARG scores, and survival status.

We investigated the differences in the TME of CESC patients classified into high- and low-risk categories. We conducted a quantitative analysis of immune cell proportions using CIBERSORT R Scripts, ordering CESC samples from lowest to highest risk scores to reflect the variability in immune cell types ( Figure 5A ). Notably, a rising trend in activated mast cells correlated with increasing risk scores (R = 0.25) ( Figure 5B ). Additionally, eight other immune cell types showed correlations with risk scores ( Supplementary Figure S1 ). In high-risk patients, activated mast cells constituted a more significant portion of the TME ( Figure 5C ), suggesting their influence on adverse prognoses. This association of various immune cells with CESC risk categories enhances our understanding of TME composition in these tumors ( Figure 5D ). The nine-gene ARG score model demonstrated distinct expression patterns between risk groups, corresponding to variations in immune cell infiltration ( Figures 5E, F ). Differential stromal and immune scores, derived from expression profile evaluations, were observed between the high- and low-risk groups ( Figure 5G ). Moreover, using the ‘oncoPredict’ R package (Version: 0.2), we identified differential drug susceptibility profiles across these groups, highlighting potential therapeutic avenues ( Supplementary Table S2 ; Supplementary Figures S2, S3 ).

We integrated the ARG score with the clinical and pathological characteristics of CESC patients to develop a comprehensive nomogram ( Figure 6A ). The accuracy and prognostic efficacy of this nomogram were confirmed by the calibration plot ( Figure 6B ). The cumulative risk curve ( Figure 6C ) showed a progressive increase in overall survival risk for CESC patients with higher nomogram scores, underscoring its prognostic value. Decision Curve Analysis demonstrated the ability of the nomogram to predict short- and long-term survival outcomes in CESC patients ( Figure 6D ). The forest plot identified the T-stage, N-stage, and risk score as crucial risk factors within the nomogram ( Figure 6E ). Survival analyses comparing T1N0 vs. T1N1 and T2N0 vs. T2N1 indicated significantly lower survival rates for patients with N1 during T1 compared to N0 (P<0.001) ( Figure 6F ). These findings collectively validate our ARG-based nomogram as a robust tool for clinical prognosis in CESC patients, significantly contributing to personalized patient care strategies.

Mapping the expression of anoikis-associated genes to different cell types in the tumor microenvironment elucidated the interactions between these genes in tumor, stromal, and immune cells, crucial for understanding anoikis resistance and cancer progression. We analyzed the expression profiles of eight ARGs in the TME using the CESC_GSE168652 single-cell dataset from TISCH. This dataset identified 21 cell clusters across seven distinct cell types ( Figure 7A ). ABL1, BCL2, and IGF1 exhibited high expression in fibroblasts and smooth muscle cells (SMCs). BAX and PLAU were highly expressed in malignant cells, while FASN showed broad expression across multiple cell types. EDA2R had low expression in various cell types, whereas NTRK3 was highly expressed in endometrial stromal cells ( Figures 7B, C ). Hallmark analysis revealed lower apoptosis levels in tumor cells, marked by BAX and PLAU, indicating their unique roles ( Figure 7A ). Since MIR200A was not expressed in the cell clusters, it was excluded from the analysis.

Immunohistochemical analysis utilizing the HPA database revealed a notable disparity in gene expression within cervical cancer tissues. Specifically, the genes ABL1, BAX, FASN, and PLAU exhibited abnormally high expression levels in cervical cancer specimens, in stark contrast to BCL2, IGF1, and NTRK3, which were markedly under-expressed ( Figure 8A ). This differential expression pattern was further corroborated by survival analysis. The survival curves, based on the expression of each gene, provided additional evidence supporting the observed expression trends in relation to patient survival outcomes ( Figure 8B ). Further in vitro experiments indicated that, compared to Ect1/e6e7 cells, Hela cells exhibited upregulated mRNA levels of BAX, FASN, and PLAU, while BCL2, IGF1, and NTRK3 mRNA levels were downregulated. ABL1 and EDA2R did not show significant differences ( Figure 9A ). To validate these findings further, we analyzed GEPIA data, which showed similar results. Compared with normal tissue, BAX, FASN, and PLAU were significantly upregulated in CESC, but ABL1 and EDA2R did not achieve the expected results ( Figure 9B ). Meanwhile, BCL2, IGF1, and NTRK3 expression were significantly reduced in the tumors.