To better understand the research process, Figure 1 illustrates the workflow diagram of this study. RNA-seq data were acquired from publicly accessible Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases (12, 13). The TCGA dataset contains clinical information from 296 cancer patients, including disease stage and survival outcomes, along with 3 normal tissue controls. Key characteristics of both datasets are systematically summarized in Table 1 . DEGs were rigorously identified using the “limma” package in R, with selection criteria of |log2 FC|>1 and adjusted P<0.05 ( 14).

DEGs were subjected to functional enrichment analysis using the DAVID database (https://david.ncifcrf.gov/) for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations (15). The “clusterProfiler” package in R facilitated statistical evaluation and visualization of enriched biological terms (16, 17). GO enrichment analysis includes three key biological domains: biological processes (BP), cellular components (CC), and molecular functions (MF) (18), while KEGG pathway annotation elucidated predominant metabolic and signaling pathways. Enrichment significance was determined at P<0.05.

Feature selection employed LASSO regression and Support Vector Machine-Recursive Feature Elimination (SVM-RFE) (19, 20). The LASSO algorithm (“glmnet” package) identified CESC-associated genes through penalized regression (α=1), with optimal regularization parameter λ determined by five-fold cross-validation. SVM-RFE (“e1071” package) iteratively eliminated low-contribution features to extract disease-specific signatures. Prognostically significant differentially expressed genes were subsequently screened via univariate Cox regression (P<0.01), followed by multivariate Cox regression to construct the final predictive nomogram (21).

The risk score, derived from the constructed prognostic model, represents the prognostic risk for each CESC patient. The calculation formula for the model’s gene risk score is as follows: gene a coefficient multiplied by gene an expression, plus gene b coefficient multiplied by gene b expression,…, plus gene i coefficient multiplied by gene i expression (where a, b, and i represent specific genes). Training and testing cohorts were stratified into high-/low-risk groups using cohort-specific median thresholds. Survival disparities between risk strata were statistically validated through Kaplan-Meier analysis.

Immune cell infiltration levels of 22 distinct subtypes were quantified through transcriptomic deconvolution using the “CIBERSORT” algorithm implemented in R. Comparative analysis of immune profiles between high- and low-risk cohorts was subsequently performed with the “limma” package.

Gene set enrichment analysis was conducted using the Molecular Signatures Database (MSigDB) collections “c5.go.v7.4.symbols” and “c2.cp.kegg.v7.4.symbols”. Hub genes associated with specific immune cell populations were functionally annotated through the “clusterProfiler” and “enrichplot” packages in R, with statistically significant terms identified at P<0.05 ( 16).

This retrospective analysis included 6 hysterectomy patients stratified into two cohorts: cervical squamous carcinoma (CESC; n=3, mean age 57.02 ± 9.21 years) and benign gynecological conditions (adenomyosis/uterine fibroids/prolapse; n=3, 57.56 ± 9.81 years), all confirmed by histopathological examination. Surgical specimens were snap-frozen in liquid nitrogen within 30 min post-resection for RNA preservation. Exclusion criteria encompassed pre-existing metabolic disorders (diabetes mellitus) or renal dysfunction. The research protocol was approved by the Second Affiliated Hospital of Guangdong Medical University (Approval No. PJKT2024-134).

Gene expression analysis was performed using TRIzol reagent (Invitrogen, USA) to extract 1 μg total RNA per manufacturer’s protocol. RNA was reverse-transcribed into cDNA and amplified via qRT-PCR under standardized cycling conditions: 95°C for 5 min, followed by 40 cycles of 95°C for 15 s and 60°C for 30 s. Primer sequences used for qPCR are detailed in Table 2 . Gene expression levels were quantified using the 2^–ΔΔCT method and normalized with GAPDH.

A prognostic nomogram integrating independent risk factors identified through univariate and multivariate Cox regression analyses was developed using the rms package in R, enabling visualization of 3-, 5-, and 8-year overall survival (OS) predictions for CESC patients. Model validation comprised temporal discrimination assessment via time-dependent receiver operating characteristic (ROC) curves with area under the curve (AUC) quantification, calibration curve analysis for prediction accuracy evaluation.

Clinical utility of the nomogram was systematically evaluated through decision curve analysis (DCA). Using the ROC curve, we determine the optimal threshold for each patient’s risk score via the nomogram. Subsequently, patients in the training and validation queues are classified into high-risk and low-risk categories based on their calculated risk scores. Kaplan-Meier survival curves are employed to analyze OS between these groups in two cohorts to assess survival differences.

Statistical analyses were conducted using R software (version 4.3.1). Continuous and categorical variables were compared between groups using Student’s t-tests and Pearson’s chi-square tests, respectively. Survival outcomes were evaluated through Kaplan-Meier methodology with log-rank testing for group comparisons, supplemented by Cox proportional hazards regression modeling. Prognostic predictors were identified through univariate/multivariate Cox regression. Model discriminative ability was quantified using C-index, while calibration curves assessed prediction-observation agreement. Clinical decision-making utility was evaluated through DCA. Statistical significance in this study was determined using P<0.05.

We retrieved mRNA expression profiles from GEO databases (GSE9750, GSE39001, GSE122697) to identify DEGs in CESC. Following batch calibration and normalization, differential expression analysis of the GEO dataset revealed 666 DEGs, comprising 434 upregulated genes (log2 FC>1) and 232 downregulated genes (log2 FC< -1), as shown in ( Figures 2A, B ). To further explore DEGs in CESC, we performed a detailed differential expression analysis using the TCGA dataset, comparing cancer tissues with adjacent normal tissues. This analysis identified 5,784 significantly differentially expressed genes, including 3,452 upregulated genes (log2 FC>1) and 2,332 downregulated genes (log2 FC<-1), as presented in ( Figures 2C, D ). By integrating the GEO and TCGA datasets, we identified 112 overlapping DEGs, consisting of 20 upregulated genes (log2 FC> 1) and 92 downregulated genes (log2 FC< -1), illustrated in ( Figure 2E ). This study conducted integrative differential gene expression analysis across GEO and TCGA cohorts, characterizing tumor-specific transcriptional dysregulation through systematic comparison of neoplastic and normal tissues.

To investigate the functional characteristics of 112 DEGs in CESC, systematic enrichment analyses were performed through the DAVID database. GO and KEGG analyses were conducted with FDR<0.05. Expressly, BP analysis indicated significant involvement of these DEGs in processes such as DNA replication, chromosome segregation, cell division, and the mitotic cell cycle. CC analysis revealed their primary association with spindle poles, kinetochores, and chromosome regions. Furthermore, MF analysis demonstrated that these DEGs primarily exhibit binding abilities to single-stranded DNA, microtubules, and protein kinases ( Figure 3A ). KEGG pathway analysis ( Figure 3B ) identified core mechanisms including Cell cycle, DNA replication, Oocyte meiosis, Progesterone-mediated oocyte maturation, Motor proteins, p53 signaling pathway, Cellular senescence, Mismatch repair, and Base excision repair. These findings collectively implicate the DEGs in essential oncogenic processes: genomic stability maintenance through DNA replication/repair, mitotic regulation via spindle apparatus components, and potential involvement in cellular senescence mechanisms in CESC pathogenesis.

Disease-associated genes were identified through a multi-step feature selection process. LASSO regression analysis selected 23 characteristic genes (e.g., CXCL8, EFNA1, EZH2) by minimizing prediction error, with error magnitude and gene number relationships visualized ( Figures 4A, B ). SVM-RFE analysis subsequently refined candidate genes through recursive elimination, identifying 51 optimal features at minimal cross-validation error ( Figure 4C ). Intersection analysis of both methods revealed eight core disease-characteristic genes (CXCL8, EFNA1, EZH2, PAQR4, SLC27A6, SPINK5, SYCP2, YEATS2), confirmed by Venn diagram ( Figure 4D ).

Feature selection was conducted through LASSO regression (“glmnet” package) to identify feature genes demonstrating minimal prediction error. The optimal regularization parameter (λ) was determined via cross-validation, selecting hub genes with the lowest error rate as disease-specific biomarkers ( Figures 5A, B ). Furthermore, a risk-scoring model for cervical cancer was constructed and used to stratify 296 CESC patients into low-risk (n=148) and high-risk (n=148) groups based on their median disease risk, as shown in ( Figure 5C ). For further insight, ( Figure 5D ) displays the risk score distribution and corresponding survival status of patients within these risk groups. Finally, univariate Cox regression analysis was conducted to identify four key genes (F10, PPP1R14A, FBLN5, ABCA8), and a gene model was constructed based on these findings. The analysis results are presented clearly in ( Figure 5E ).

To investigate the prognostic significance of 12 key genes (EFNA1, CXCL8, EZH2, PAQR4, SLC27A6, SPINK5, SYCP2, YEATS2, F10, PPP1R14A, FBLN5, ABCA8) in CESC, we conducted a comprehensive survival analysis. Notably, EFNA1 ( Figure 6A ), CXCL8 ( Figure 6B ), and PPP1R14A ( Figure 6C ) demonstrated statistically significant predictive value for overall survival (P<0.05). Conversely, EZH2, PAQR4, SLC27A6, SPINK5, SYCP2, YEATS2, F10, FBLN5, and ABCA8 did not exhibit statistically significant differences in overall survival prognosis ( Figures 6D–L ). These findings suggest that EFNA1, CXCL8, and PPP1R14A may serve as potential prognostic markers for CESC, warranting further investigation for their roles in disease progression and therapeutic targeting.

Time-dependent ROC analysis evaluated the prognostic accuracy of a three-gene signature (EFNA1, CXCL8, PPP1R14A) for cervical squamous cell carcinoma survival outcomes. The EFNA1 ( Figure 7A ) demonstrated moderate discrimination with progressive improvement in AUC from 0.650 (1-year) to 0.715 (10-year), consistently exceeding random prediction. CXCL8 ( Figure 7B ) showed superior and stable performance across all intervals, maintaining AUC between 0.682-0.719 with minimal temporal variation. Notably, PPP1R14A ( Figure 7C ) exhibited temporal performance heterogeneity - achieving peak discriminative capacity at 1-year (AUC=0.726) followed by progressive decline to 0.505 at 10-year follow-up. This temporal analysis revealed two distinct prognostic patterns: CXCL8 maintained consistent predictive reliability throughout the observation period, while PPP1R14A and EFNA1 showed opposing temporal trajectories. The composite model demonstrated the strongest predictive utility for early-stage prognosis (1–3 year AUC range: 0.650-0.726), with diminished accuracy in long-term predictions (5–10 year AUC range: 0.505-0.715). These findings position this multigene signature as a potential tool for short-term survival stratification, though longitudinal prognostic accuracy may require temporal recalibration or incorporation of complementary biomarkers.

Comprehensive immune microenvironment profiling of CESC was performed through analysis of 22 immune cell signatures and infiltration patterns ( Figures 8A, B ). The heatmap visually represents the abundance and proportional differences of various immune cell subtypes in CESC samples, highlighting a higher proportion of CD8+ T cells, plasma cells, Tregs, follicular helper T cells, and macrophages (M1, M0, and M2). To explore the correlation between the three key genes (EFNA1, CXCL8, and PPP1R14A) and immune cell infiltration, we stratified the gene expression levels into high- and low-risk groups for analysis ( Figures 8C-E ). The findings revealed significant differences in dendritic cell activation between high- and low-risk groups for EFNA1 (P<0.05). CXCL8 significantly influenced the proportions of Tregs, NK cells, dendritic cell activation, mast cells, and neutrophils (P<0.05), and changes in PPP1R14A were significantly associated with γδ T cells, macrophages (M0 and M2), dendritic cell activation, and the proportion of neutrophils (P<0.05). These mechanistic insights propose clinically actionable strategies: CXCL8 inhibition may disrupt immunosuppressive networks while PPP1R14A modulation could target macrophage polarization, with EFNA1-mediated dendritic cell activation serving as a potential predictive biomarker for immune checkpoint inhibitor response.

Gene co-expression network analysis revealed distinct interaction patterns among EFNA1, CXCL8, and PPP1R14A. Pearson correlation analysis demonstrated a statistically robust positive correlation between EFNA1 and CXCL8 expression (R=0.21, P<0.001; Figure 9A ), suggesting potential co-regulatory mechanisms or functional synergy. In contrast, non-significant correlations were observed for EFNA1-PPP1R14A (R=0.053, P=0.36; Figure 9B ) and CXCL8-PPP1R14A pairs (R=0.042, P=0.46; Figure 9C ), as evidenced by nullcline-proximal data distributions. This differential correlation profile implies PPP1R14A’s functional autonomy from the EFNA1-CXCL8 axis within the studied biological context.

GSEA systematically decoded functional networks of EFNA1, CXCL8, and PPP1R14A within the CESC immune microenvironment. EFNA1 exhibited dual regulatory capacity: GO terms implicated its involvement in keratinocyte differentiation and T cell receptor complex assembly ( Figure 10A ), while KEGG pathways connected it to MAPK signaling and cytoskeletal remodeling ( Figure 10B ). CXCL8 exhibited dual regulatory roles in cervical carcinogenesis, demonstrating pro-proliferative effects on endothelial/epithelial lineages ( Figure 10C ) concurrent with inflammasome activation through NOD-like receptor signaling pathways ( Figure 10D ). PPP1R14A emerged as a stromal interface regulator, with GO enrichment in fibroblast growth factor signaling and extracellular matrix (ECM) organization ( Figure 10E ), corroborated by KEGG pathways encompassing ECM-receptor interactions and cancer-associated adhesion molecules ( Figure 10F ). Multi-dimensional analysis reveals three synergistic pathological mechanisms in cervical squamous carcinogenesis: EFNA1-mediated immune-structural interactions regulate differentiation processes, CXCL8 coordinates proliferative-inflammatory equilibrium, and PPP1R14A drives stromal-ECM remodeling. These interconnected networks collectively promote tumor progression through differentiation dysregulation, immune evasion, and metastatic niche formation, identifying microenvironment-specific targets for precision immunotherapy development.

To investigate the role of key genes in cervical cancer survival prognosis, we comprehensively analyzed the correlation between the expression levels of EFNA1, CXCL8, and PPP1R14A, and various clinical features, including age, stage, tumor grade (T), lymph node status (N), and metastasis status (M). EFNA1 demonstrated lymph node-specific regulation with elevated expression in NX versus N1 cases (P=0.045, Figures 11A, E ), showing no significant associations with age, stage, grade, or metastasis ( Figures 11B-D, F ). CXCL8 exhibited progressive upregulation in advanced disease stages (T3 vs T2, P=0.032; G3 vs G2, P=0.016; Figures 11G-J ), independent of age or metastatic status ( Figures 11K, L ). PPP1R14A displayed age-dependent suppression (>65 years, P=0.032) and paradoxical elevation in TX-grade tumors (vs T1-2, P=0.023; Figures 11M, N, P ), with no significant correlations to staging or metastasis ( Figures 11O, Q, R ). In short, EFNA1 demonstrates lymph node-specific biomarker potential, CXCL8 correlates with histopathological progression through stage/grade associations, and PPP1R14A exhibits age-related suppression and grade-dependent paradoxical expression. These differential clinical signatures collectively highlight their translational value for precision staging systems and biomarker-driven therapeutic stratification in CESC management.

A prognostic nomogram integrating three molecular markers (EFNA1, CXCL8, PPP1R14A) with clinicopathological parameters (age, TNM stage) was developed to predict 3-/5-/8-year overall survival in cervical squamous cell carcinoma ( Figure 12A ). Tumor grade demonstrated the strongest prognostic weight (C-index: 0.785 training, 0.750 validation), followed by CXCL8 expression and metastasis status. The nomogram showed robust temporal discrimination (3-year AUC: 0.762 training, 0.689 validation; Figures 12B, C ) and precise calibration ( Figures 12D-I ), with DCA confirming clinical utility across 10-45% risk thresholds ( Figures 12J-O ). This multivariable tool enables personalized survival prediction and risk-stratified therapeutic decision-making for CESC patients.

Cervical cancer primarily consists of squamous cell carcinoma (70%), adenocarcinoma (25%), and other rare types (5%) (22). Multimodal analysis of cervical carcinogenesis progression was conducted through histopathological evaluation and molecular profiling. Hematoxylin-eosin staining revealed progressive histoarchitectural alterations ( Figure 13A ): Control tissues maintained normal stratified squamous epithelium, LSIL specimens exhibited basal layer expansion with koilocytic changes, HSIL showed full-thickness dysplasia with nuclear hyperchromasia, and invasive carcinoma displayed complete architectural disarray with stromal invasion. Complementary qRT-PCR analysis demonstrated significant transcriptional upregulation of oncogenic effectors - EFNA1, CXCL8, and PPP1R14A in carcinoma versus control tissues (P<0.05) ( Figure 13B ). This coordinated overexpression pattern correlates with histopathological progression from premalignant lesions to invasive carcinoma, suggesting their synergistic role in cervical carcinogenesis.