Identifying and validating of prognostic genes associated with myeloid cell differentiation in cervical cancer: development of a risk model based on single-cell RNA sequencing combined with bulk RNA sequencing data

Introduction

Cervical cancer (CESC) is a leading gynecological malignancy, accounting for approximately 660,000 new cases and 350,000 deaths annually worldwide (1,2). Despite advances in human papillomavirus (HPV) vaccination and screening, it remains a significant global public health challenge, with a disproportionately high burden in developing countries—accounting for 84% of new cases and 88% of related deaths (3,4). The 5-year survival rate in regions like sub-Saharan Africa can be as low as 35%, far below the global average (5). This dire situation underscores the critical need for improved prognostic tools and more effective therapeutic strategies.

The progression of CESC is a complex process governed by both cancer cells and the tumor microenvironment (TME). Within the TME, myeloid cells are key innate immune players, whose functional polarization profoundly influences cancer progression (6-8). Recruited and educated by tumor-derived factors (e.g., CSF-1, CCL2), myeloid precursors undergo epigenetic and metabolic reprogramming (9). These cells predominantly differentiate into immunosuppressive M2-like tumor-associated macrophages (TAMs) (10), which secrete factors like transforming growth factor-beta (TGF-β) and interleukin-10 (IL-10) to inhibit cytotoxic T cells, while also promoting angiogenesis and stromal support via vascular endothelial growth factor (VEGF), thereby fostering a pro-tumorigenic microenvironment (11,12). Aberrant myelopoiesis within the TME drives the accumulation of myeloid-derived suppressor cells (MDSCs) and polarization of macrophages toward an M2 phenotype, collectively establishing a potent immunosuppressive niche (13). These myeloid populations suppress cytotoxic T cell activity (14), facilitate immune evasion (15), and promote metabolic reprogramming and epithelial-mesenchymal transition (EMT) (16), ultimately contributing to tumor progression and poor clinical outcomes (15,17). Consequently, myeloid cell differentiation (MCD) is central to immune escape in CESC and represents a promising therapeutic frontier. However, current prognostic models for CESC largely rely on bulk transcriptomic data, which obscure the substantial cellular heterogeneity within the myeloid compartment. While the overall role of TAMs has been well documented, the dynamic gene expression programs governing MCD and their context-specific regulatory networks in CESC remain incompletely understood. The lack of high-resolution, differentiation-associated signatures consequently limits the precision of individualized risk stratification.

Therefore, we aimed to identify and validate myeloid cell differentiation-associated genes (MCD-AGs) with prognostic value in CESC. We integrated transcriptomic data from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO), and applied bioinformatics pipelines—including differential expression, survival, and machine learning analyses—to construct a reliable risk prediction model. Furthermore, we leveraged single-cell RNA sequencing (scRNA-seq) to delineate the expression dynamics of these genes during myeloid differentiation. Our findings were experimentally validated using clinical samples via reverse transcription quantitative polymerase chain reaction (RT-qPCR) and immunohistochemistry (IHC). Our study thus provides a novel prognostic tool and offers mechanistic insights into CESC progression from the perspective of MCD, potentially informing future immunotherapeutic strategies. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0174/rc).

Methods

Data acquisition

The overall research workflow is summarized in Figure 1. We obtained RNA-seq, clinical, mutational, and survival data for 275 CESC patients and three matched normal samples from TCGA-CESC dataset (18). For scRNA-seq analysis, we acquired the GSE168652 dataset (GPL24676) from the GEO, comprising one tumor and one matched normal sample (19). The E-MTAB-11948 dataset was retrieved from ArrayExpress, which included three tumor and two normal samples after quality control (one sample was excluded due to an exceptionally low number of detected genes).

Figure 1 Research flowchart. DEG, differentially expressed gene; GO, Gene Ontology; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; MCD-RG, myeloid cell differentiation-related gene; MSigDB, Molecular Signatures Database; PPI, protein-protein interaction; RT-qPCR, reverse transcription quantitative polymerase chain reaction.

MCD-related genes (MCD-RGs) were obtained from Molecular Signatures Database (MSigDB; v2023.2, C5:GO:BP), selecting “GOBP_MYELOID_CELL_DIFFERENTIATION” (GO:0030099) (20). This GO-based set includes 423 genes governing myeloid lineage commitment and differentiation (e.g., monocytes, macrophages, granulocytes, dendritic cells) (table available at https://cdn.amegroups.cn/static/public/tcr-2026-1-0174-1.xlsx), ensuring a biologically robust basis for model construction.

Differential expression analysis

In the TCGA-CESC dataset, differentially expressed genes (DEGs) were identified by comparing CESC tissues to adjacent normal tissues using the DESeq2 package (v1.42.0) (21). The thresholds were set at an adjusted P value <0.05 and |log2 fold change (FC)| >2. A volcano plot was generated using the ggplot2 package (v3.4.4) (https://link.springer.com/book/10.1007/978-0-387-98141-3), and a heatmap was created with the ComplexHeatmap package (v2.16.0) (22).

Identification, functional analysis, and protein-protein interactions (PPI) of candidate genes

Candidate genes were identified as the intersection of the DEGs and MCD-RGs using the VennDiagram package (v1.2.3) (23). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were then conducted on these candidate genes using the clusterProfiler package (v4.10.0) (24), with a significance threshold of adjusted P value <0.05. PPI analysis was performed for the candidate genes using the STRING database (http://www.string-db.org/), setting a minimum confidence score threshold of >0.4. Finally, a PPI network was constructed and visualized with Cytoscape (v3.10.2) (25).

Spearman correlation analysis was conducted using TCGA-CESC transcriptome data to assess associations between four prognostic genes (TNF, PTPN6, FASN, TFRC) and myeloid cell lineages. Gene expression was correlated with markers of monocytes (CD14, FCGR3A), total macrophages (ITGAM, CSF1R), M1 (CD80, CCR7), M2 (CD163, MSR1), neutrophils (CEACAM8, MPO), and MDSCs (CD33). Significance was defined as |r| >0.3 with false discovery FDR-adjusted P<0.05.

Establishment and confirmation of risk model

For prognostic gene screening, survival data from patients with CESC in the TCGA-CESC dataset were analyzed. Genes associated with overall survival were identified using univariate Cox regression analysis via the coxph function of the survival package (v3.5-7), with the criteria of a statistically significant association (P<0.01). The proportional hazards (PH) assumption was verified using the cox.zph function, requiring P>0.05. Prognostic genes were further refined through least absolute shrinkage and selection operator (LASSO) Cox regression analysis using the glmnet package (v4.1-4) (26), with the family parameter set to “cox” and ten-fold cross-validation to determine the optimal lambda value.

To evaluate the prognostic risk model, the 275 CESC samples from TCGA were randomly divided into training and validation sets at a 7:3 ratio. This resulted in 193 samples (70%) in the training cohort for model development and 82 samples (30%) in the validation cohort for model testing. The risk score for each patient was calculated using the formula: riskscore=∑i=1kβi∗expi, where β represents the linear coefficients derived from LASSO regression and exp denotes the expression level of each prognostic gene.

Patients in both training and validation sets were stratified into high-risk group (HRG) and low-risk group (LRG) based on the optimal risk score threshold determined by the surv_cutpoint function from the survminer package. The distribution of risk scores and survival status was visualized using the ggplot2 package (v3.4.4), while expression patterns of prognostic genes were displayed using the ComplexHeatmap package (v2.16.0). Kaplan-Meier survival analysis was performed with the survminer package (v0.4.9) to compare overall survival between HRG and LRG, with statistical significance defined as P<0.05. Finally, the time-dependent receiver operating characteristic (ROC) curve was generated using the timeROC package (v0.4) (27) to assess the model’s predictive accuracy, with an area under the curve (AUC) >0.6 considered indicative of acceptable performance.

Building and assessing of a nomogram

Within the training cohort, we further integrated the risk score with clinical covariates to identify independent prognostic predictors. First, univariate Cox regression analysis (using the survival package, v3.5-7) was performed, retaining variables for further consideration with a liberal significance threshold of P<0.2. The PH assumption for each variable was confirmed using the cox.zph function (P>0.05). Subsequently, all selected variables were included in a multivariate Cox regression analysis to identify independent prognostic factors, applying the same significance threshold (P<0.2).

A nomogram was then constructed using the rms package (v6.7-0) to predict 1-, 2-, and 3-year overall survival probabilities, incorporating these independent factors. The predictive accuracy and clinical utility of the nomogram were rigorously evaluated. This included generating calibration curves (rms package, v6.7-0) to assess agreement between predicted and observed outcomes, plotting time-dependent ROC curves (timeROC package, v0.4) to quantify discriminative ability, and performing decision curve analysis (DCA) (ggDCA package, v1.1) to evaluate net clinical benefit across different decision thresholds.

Ultimately, the nomogram and the optimized risk score threshold provide a concrete tool for clinical stratification, translating the “risk” into specific “stratification criteria”. This facilitates personalized clinical decision-making for individual patients.

Clinical features analysis and gene set enrichment analysis (GSEA)

To investigate the underlying biological differences, we performed a differential expression analysis (HRG vs. LRG) on the TCGA-CESC dataset using DESeq2 (v1.42.0). The resulting genes, ranked by their log2FC, were then subjected to GSEA using the clusterProfiler package (v4.10.0). GSEA was executed against the “c2.cp.kegg.v7.5.1.symbols.gmt” (KEGG pathways) and “c5.go.v7.5.1.symbols.gmt” (GO) gene sets from MSigDB (v7.5.1), with gene sets considered significantly enriched at P<0.05.

Immunoinfiltration analysis

The immune landscape of the training cohort was analyzed using the single-sample GSEA (ssGSEA) algorithm, implemented via the GSVA R package (v1.49.4) (28), to quantify the relative abundances of 28 immune cell types (29). The Wilcoxon rank-sum test was then applied to compare the ssGSEA enrichment scores of these immune cells between the HRG and LRG, with a statistical significance threshold of P<0.05. Furthermore, Spearman correlation analysis was performed using the psych package (v2.2.5) to investigate the relationships among the differentially expressed immune cells. This analysis was based on the survival data from the training cohort, and correlations were considered significant at an absolute value of the correlation coefficient (|r|) >0.3 and P<0.05.

Tumor mutation and drug sensitivity analyses

To investigate somatic mutations between the HRG and LRG, mutation data from the CESC training set were downloaded, with Mutect2 selected as the variant caller. The data in Mutation Annotation Format (MAF) was summarized, analyzed, and visualized using the Maftools package (v2.22.0) (30). Waterfall plots were generated to display the mutational landscape and the spectrum of variant genes.

The sensitivity of patients with CESC to common chemotherapeutic agents was assessed using the Genomics of Drug Sensitivity in Cancer (GDSC) database, which includes 138 drugs. The half-maximum inhibitory concentration (IC50) values for these drugs were predicted for the HRG and LRG in the training set using the pRRophetic package (v0.5) (31). The Wilcoxon rank-sum test was then applied to compare the estimated IC50 values between the two risk groups for each drug, with a significance threshold set at P<0.05.

Drug forecasts and molecular docking analyses

To identify drugs associated with the prognostic genes, the Drug-Gene Interaction Database (DGIdb, https://dgidb.org) was queried to predict compounds targeting these genes, with a filter applied to select only “Approved” drugs. The resulting drug-biomarker relationships were visualized using Cytoscape (v3.10.2). For each biomarker, the drug with the highest interaction score was selected as a key candidate for subsequent molecular docking.

To assess the binding potential between the protein targets of the prognostic genes and the key drugs, the three-dimensional structures of the key drugs were first retrieved in Structure Data File (SDF) format from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The crystal structure of the target protein [Protein Data Bank (PDB) ID: 3hhd] was obtained from the PDB (https://www.rcsb.org/). Molecular docking was then performed using the CB-Dock2 web server (https://cadd.labshare.cn/cb-dock2/index.php) to predict the binding poses and affinities between the protein and the drugs. A Vina score of ≤−5 kcal/mol was considered indicative of strong binding affinity.

scRNA-seq analysis

To obtain high-quality cells for analysis, scRNA-seq data from the E-MTAB-11948 and GSE168652 datasets were processed using the Seurat package (v4.4.0) (32). A Seurat object was created with the parameters min.cells =3 and min.features =200. Potential doublets were identified and removed using the scDblFinder package (v1.16.0) (33). Cells were further filtered based on the following quality control criteria: percent of mitochondrial reads (percent.mt) <20%, number of detected genes (100≤ nFeature_RNA ≤6,000), and total RNA counts (200≤ nCount_RNA ≤20,000).

After quality control, highly variable genes (HVGs) were identified using the FindVariableFeatures function. The two datasets were then integrated using the IntegrateData function to merge the samples. The integrated data was scaled and subjected to principal component analysis (PCA) using the RunPCA function. To correct for batch effects between the E-MTAB-11948 and GSE168652 datasets, we applied the Harmony algorithm via the RunHarmony function (theta =2, lambda =1). Significant principal components (PCs) for downstream analysis were selected based on the Jackstraw procedure, which uses permutation testing to assess the significance of each PC. Using these significant PCs, unsupervised cell clustering was performed with the FindNeighbors and FindClusters functions. A resolution parameter of 0.2 was selected for clustering by optimizing the silhouette coefficient. Finally, the cell clusters were visualized in two dimensions using t-distributed stochastic neighbor embedding (t-SNE) with the RunTSNE function.

Cell types were annotated by comparing the expression of established marker genes from the literature (34) across the clusters. The expression patterns of these marker genes were visualized in a dot plot, and the relative abundance of each cell type in CESC versus paracancerous tissues was displayed in a bar plot.

To evaluate the expression of the prognostic genes within the single-cell atlas, we used the Wilcoxon rank-sum test to identify genes that were differentially expressed between CESC and paired normal samples at a significance level of P<0.05. Cell clusters that consistently showed significant upregulation of multiple prognostic genes were prioritized as key cell types for further investigation.

Cellular communication and pseudotime trajectory analyses

To investigate cell-cell communication between the key cell types and other annotated populations, we applied the CellChat package (v1.6.1) (35) to the CESC and paracancerous samples. Bubble plots were used to visualize the strength of ligand-receptor interactions involving the key cell types, with interactions considered significant at a probability of P<0.05.

To explore the heterogeneity within each key cell type, we performed sub-clustering on the cells belonging to each type (resolution =0.2). Marker genes for each resulting subcluster were identified using the FindAllMarkers function with the following parameters: min.pct =0.5, logfc.threshold =0.25, and test.use =“auc” (36). Based on the expression of these marker genes, the key cells were categorized into distinct subtypes. The developmental trajectories of these subtypes were then reconstructed using pseudotime analysis with the Monocle2 package (v2.30.0) (37). Finally, the plot_cells function was used to visualize the expression dynamics of the prognostic genes along the inferred pseudotime trajectory.

Expression level analysis of prognostic genes

The Wilcoxon test was conducted to compare prognostic gene expression between CESC and matched normal samples from TCGA, with significance set at P<0.05.

RT-qPCR

Five CESC tissue samples and five matched adjacent normal cervical tissue samples (as controls) were collected from The First Affiliated Hospital of Guizhou University of Traditional Chinese Medicine. All samples were immediately snap-frozen and stored at −80 ℃ until RNA extraction. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of The First Affiliated Hospital of Guizhou University of Traditional Chinese Medicine (No. KL2025-004). All participants, including those who provided both CESC and paired adjacent normal tissue samples, provided written informed consent after being fully informed of the study protocol. The manuscript contains no identifiable personal information. Total RNA was extracted from the tissue samples using TRIzol reagent (R401-01, Vazyme Biotech, Nanjing, China) following the manufacturer’s instructions. cDNA was synthesized from the extracted RNA using the Hifair^® III 1st Strand cDNA Synthesis Supermix for qPCR (11141ES60, Yeasen Biotechnology, Shanghai, China). Quantitative PCR was then performed using the 2× Universal Blue SYBR Green qPCR Master Mix (G3326-05, Servicebio Technology, Wuhan, China). The expression of prognostic genes was validated, with GAPDH serving as the endogenous control. The primer sequences used are provided in Table S1. Relative gene expression levels were calculated using the 2^−ΔΔCt method, and statistical differences between groups were assessed using a two-tailed Student’s t-test, with a significance threshold of P<0.05.

Statistical analysis

Data analysis was performed using R software (v4.2.2). Continuous data are presented as mean ± standard deviation (SD). Differences between groups were analyzed using the Wilcoxon rank-sum test or Student’s t-test, as appropriate, with statistical significance defined as P<0.05.

Results

Identification, functional enrichment, and PPI network of candidate genes in CESC

A total of 2,780 DEGs were identified between CESC and paracancerous specimens from the TCGA-CESC dataset. Among these, 1,681 were upregulated and 1,099 were downregulated in CESC samples. The top 10 upregulated and downregulated DEGs, ranked by adjusted P value, were visualized in a volcano plot (Figure 2A), while a heatmap illustrated the overall gene expression patterns across all samples (Figure 2B). The intersection of these DEGs with 423 MCD-RGs yielded 73 candidate genes for subsequent analysis (Figure 2C).

Figure 2 Identification of candidate genes, functional enrichment analysis, and PPI network in CESC. (A) Volcano plot of DEGs. Each dot represents a gene, with blue dots indicating downregulated genes and red dots indicating upregulated genes. The two dashed lines represent the thresholds used for screening DEGs. (B) Heatmap displaying the distribution of expression levels across samples, with colors transitioning from blue to red to indicate increasing expression levels. (C) Venn diagram of candidate genes. Orange represents the DEGs between cervical cancer samples and paracancerous tissues, while light purple indicates genes associated with myeloid cell differentiation. (D) GO enrichment results for candidate genes (top 10). (E) KEGG enrichment results for candidate genes (top 10). In both (D) and (E), the x-axis represents the number of genes enriched in each pathway, and the y-axis represents the specific pathways. The color indicates the level of statistical significance. (F) Construction of the PPI network for candidate genes. Blue boxes represent the candidate genes, and gray lines illustrate the interactions among these genes. BP, biological process; CC, cellular component; CESC, cervical cancer; DEGs, differentially expressed genes; FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MCD-RG, myeloid cell differentiation-related gene; MF, molecular function; PPI, protein-protein interaction.

Functional enrichment analysis of the 73 candidate genes revealed significant involvement in 1,164 GO terms and 47 KEGG pathways (adjusted P<0.05, Figure 2D,2E; table available at https://cdn.amegroups.cn/static/public/tcr-2026-1-0174-2.xlsx, Tables S2). The GO terms comprised 1,136 biological processes (BPs), 25 molecular functions (MFs), and 3 cellular components (CCs). Notable enriched terms included MCD and its regulatory processes in the BP category; the RNA polymerase II transcription regulator complex and the plasma membrane external side in the CC category; and cytokine receptor binding and kinase regulator activity in the MF category. Significant KEGG pathways included the JAK-STAT signaling pathway and inflammatory bowel disease.

A PPI network was constructed, comprising 73 nodes and 237 interactions, with 15 proteins remaining as isolated nodes (Figure 2F). Key interactions within the network included those between FASN and FLNA, as well as between PTPN6 and LILRB4.

Development and verification of risk models in CESC

The univariate Cox regression analysis identified four genes that were significantly associated with the prognosis of CESC (P<0.01, Figure 3A). The PH assumption for these genes was verified using the cox.zph function, with all genes satisfying the assumption (P>0.05, Figure 3B). These four candidate genes—TNF, PTPN6, FASN, and TFRC—were subsequently refined as the final prognostic signature through LASSO Cox regression, selected at the optimal penalty parameter [log(lambda.min) =−3.998] that minimized the cross-validation error (Figure 3C,3D). To further link this signature to myeloid cell dynamics (Figure S1), correlation analysis showed that only PTPN6 was significantly associated with myeloid markers (|r|>0.3, P_adj<0.05), strongest with CD33 (r=0.48, P_adj=5.74×10⁻¹⁶), and also with ITGAM (0.44), CSF1R (0.42), CCR7 (0.42), CD80 (0.36), CD163 (0.30), MSR1 (0.37), CD14 (0.40), and FCGR3A (0.40), whereas TNF, FASN, and TFRC showed no significant correlations (|r|≤0.3), supporting a specific link between PTPN6 and myeloid subtypes, particularly MDSCs and TAMs.

Figure 3 Screening diagram of CESC risk model. (A) Forest plot of univariate Cox analysis. A HR greater than 1 indicates an increased risk, while an HR less than 1 suggests a decreased risk. (B) PH test. The solid line shows the smooth fit of the residuals with respect to time. (C) LASSO regression. The x-axis represents log(lambda), and the y-axis represents the cross-validation error. Each line corresponds to a gene, with the dashed line on the left indicating the point where the cross-validation error reaches its minimum value. (D) Screening of prognostic genes through LASSO regression. The x-axis represents log(lambda), and the y-axis represents the coefficient of each gene. The dashed line marks the optimal lambda value [log(lambda.min) =−3.998], where four variables with non-zero coefficients are selected. CESC, cervical cancer; CI, confidence interval; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; PH, proportional hazards.

A risk score for each patient was then calculated using a linear combination of the expression levels of these four genes, weighted by their respective regression coefficients derived from the LASSO analysis (Table S3). The formula was as follows: Risk score = (0.366 × TNF) + (−0.612 × PTPN6) + (0.245 × FASN) + (0.267 × TFRC). Using dataset-specific optimal cutoffs determined by the surv_cutpoint function, patients in the entire TCGA-CESC cohort (cutoff =1.063), the training set (cutoff =1.131), and the validation set (cutoff =0.969) were stratified into HRG and LRG.

The distribution of risk scores and survival status demonstrated a higher proportion of deceased patients in the HRG across all datasets, indicating a strong association between the risk score and poor survival outcomes (Figure 4A). Consistent with this, expression heatmaps revealed a distinct pattern characterized by lower PTPN6 and higher TNF, FASN, and TFRC expression in the HRG (Figure 4B). Kaplan-Meier survival analysis confirmed that patients in the HRG had significantly shorter overall survival in the full TCGA cohort (P<0.0001), the training set (P<0.0001), and the validation set (P=0.047, Figure 4C). Finally, time-dependent ROC analysis showed that the risk model possessed good predictive accuracy for 1-, 2-, and 3-year overall survival, with all AUC values exceeding 0.6 (Figure 4D), underscoring the model's robust prognostic potential.

Figure 4 Construction and validation plot of the CESC risk model. (A) Risk score distribution and survival status in the training and validation sets. Red dots indicate high-risk patients and blue dots indicate low-risk patients; survival status is shown below (red, deceased; blue, alive). (B) Heatmaps showing expression patterns of the four prognostic genes in the training and validation sets (red, high expression; blue, low expression). (C) Kaplan-Meier survival curves comparing overall survival between high- and low-risk groups in the training and validation sets. Shaded areas represent 95% confidence intervals. (D) Time-dependent ROC curves evaluating the predictive performance of the model at 1, 2, and 3 years in the training and validation sets. AUC, area under the curve; CESC, cervical cancer; CI, confidence interval; FP, false positive; ROC, receiver operating characteristic; TP, true positive.

Independent prognostic and nomogram construction of CESC

Univariate Cox regression analysis was performed to screen for prognostic factors, using a liberal P value threshold of <0.2 to minimize the omission of potentially relevant variables. Both the risk score and age met this criterion and were therefore included in the subsequent multivariate analysis (Figure 5A). The PH assumption was confirmed for both variables (P>0.05, Figure 5B). Multivariate Cox regression analysis, applying the same screening threshold, identified both the risk score and age as independent prognostic factors for CESC (Figure 5C).

Figure 5 Independent prognostic factors of CESC and construction of nomogram. (A) Forest plot of univariate analysis. (B) Results of the univariate PH test. (C) Forest plot of multivariate analysis. (D) Nomogram. Independent prognostic factors are presented on the left side of the figure. The scale on the right side represents the range of available values for each factor. The length of the line segment indicates the magnitude of each factor’s contribution to cervical cancer. “Point” refers to the individual score corresponding to each independent prognostic factor under different value settings. “Total score” is the cumulative score obtained by summing individual scores of all factors. “Pr” denotes the probability of mortality. (E) Calibration curve of the nomogram. (F) ROC curve of the nomogram. (G) DCA curve of the nomogram. The abscissa represents the threshold probability, and the ordinate represents the net benefit, calculated by subtracting the disadvantages from the advantages. Different slanted curves correspond to different clinical diagnostic models. The “None” curve represents no treatment (Pi < Pt), resulting in a net benefit of zero. The “All” curve indicates all samples are treated, with a net benefit represented by a negatively sloped line. The “Model” curve represents the nomogram model. AUC, area under the curve; CESC, cervical cancer; CI, confidence interval; DCA, decision curve analysis; FP, false positive; OS, overall survial; PH, proportional hazards; ROC, receiver operating characteristic; TP, true positive.

These independent factors were incorporated to construct a nomogram for predicting 1-, 2-, and 3-year overall survival probabilities in patients with CESC (Figure 5D). The calibration curves for these time points showed excellent agreement between the nomogram’s predictions and the observed survival outcomes, indicating high calibration accuracy (Figure 5E). The time-dependent ROC analysis demonstrated that the nomogram had good predictive discrimination, with AUC values of 0.758, 0.877, and 0.693 for 1, 2, and 3 years, respectively (Figure 5F). Furthermore, DCA revealed that the use of this nomogram provided a superior net clinical benefit across a wide range of risk thresholds compared to alternative strategies, underscoring its potential clinical utility (Figure 5G). Collectively, these results indicate that the nomogram is a robust tool for prognostic prediction in CESC.

Clinical features, functional enrichment and immune microenvironment of CESC

The clinical heatmap illustrates the distribution of factors such as age, stage, follow-up duration, and survival status between the HRG and LRG in CESC (Figure 6A). GSEA was performed to compare the HRG and LRG, revealing significant differences in biological pathway activity. A total of 55 KEGG pathways were significantly enriched (FDR <0.05), with notable examples including the TGF-βsignaling pathway, primary immunodeficiency, and retinol metabolism (Figure 6B, table available at https://cdn.amegroups.cn/static/public/tcr-2026-1-0174-3.xlsx). Similarly, 1,216 GO terms were significantly enriched, with key terms related to synapse organization, antigen receptor-mediated signaling, and leukocyte-mediated cytotoxicity (Figure 6C, table available at https://cdn.amegroups.cn/static/public/tcr-2026-1-0174-4.xlsx). These findings associate the TGF-β signaling pathway and synaptic organization in shaping an unfavorable immune landscape, where they may facilitate the formation of an MCD-derived immunosuppressive barrier in high-risk CESC patients.

Figure 6 Clinical features, functional enrichment, and immune microenvironment map of CESC. (A) Correlation between risk scores and various clinical characteristics of patients with cervical cancer. (B) Top 5 KEGG pathway enrichment results from GSEA. (C) Top 5 GO pathway enrichment results from GSEA. The upper part of (B) and (C) shows the ES curve. The peak values correspond to the ES of the respective pathway, while the middle section indicates the gene positions within the pathway along the ranked list, with each vertical line representing a gene. (D) Heatmap of immune cell content, calculated using ssGSEA. Each grid represents the relative content of cells in the sample, with red indicating higher content. (E) Differential immune cell analysis between high- and low-risk groups. (F) Correlation between differential immune cells. Circle size represents the magnitude of correlation coefficients. ns, not significant; *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. CESC, cervical cancer; ES, enrichment score; GO, Gene Ontology; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MDSC, myeloid-derived suppressor cell; N, node; ssGSEA, single-sample GSEA; T, tumor.

Analysis of immune cell infiltration revealed significant differences in the abundance of 11 immune cell types between the HRG and LRG (P<0.05), including activated B cells, eosinophils, and activated CD8+ T cells (Figure 6D). The relative scores for these cell types are visualized in a separate panel (Figure 6E). Correlation analysis further showed strong positive relationships (correlation coefficient r>0.3, P<0.05) among several of the differentially infiltrated immune cells. Notably, a very strong correlation was observed between immature B cells and activated B cells (r=0.89, P<0.001), as well as between MDSCs and immature B cells (r=0.91, P<0.05) (Figure 6F, Tables S4,S5). This correlative pattern suggests a potential coordinated role for immature B cells, activated B cells, and MDSCs in the context of CESC prognosis.

Tumor mutations and drug sensitivity of CESC

Somatic mutation analysis of the TCGA-CESC cohort revealed that the top three most frequently mutated genes were consistent between the HRG and LRG, namely TTN, PIK3CA, and KMT2C, with comparable mutation frequencies (Figure 7A,7B). In both groups, missense mutations were the most predominant variant classification. These results indicate a broadly similar mutational landscape between the two risk groups.

Figure 7 Tumor mutation and drug sensitivity analysis of CESC. (A) Somatic mutation spectrum of the high-risk group. (B) Somatic mutation map of the low-risk group. The number of samples (as a percentage) is on the right vertical axis. (C) Correlation analysis between high- and low-risk groups and IC₅₀. (D) Box plot of the top 10 differences in drug IC₅₀ values between high- and low-risk groups. The plot shows the IC₅₀ data for various compounds in both risk groups. ****, P<0.0001. CESC, cervical cancer; FDR, false discovery rate; IC₅₀, half-maximum inhibitory concentration; TMB, tumor mutational burden.

We next sought to identify potential therapeutic agents by comparing the predicted drug sensitivity (IC50 values) between the HRG and LRG. The analysis revealed that 85 drugs exhibited significantly different IC50 distributions between the two groups (FDR <0.05). Specifically, patients in the HRG were predicted to be more sensitive to several agents, including A.770041, BIRB.0796, and BMS.708163 (Figure 7C,7D; Table S6). These findings highlight a distinct drug sensitivity profile for the HRG and suggest that these pharmacological agents represent candidate compounds for targeted therapy in this patient subset.

Drugs and molecular docking of 4 prognostic genes

Querying the DGIdb identified 59 drugs that are known or predicted to interact with the protein products of TNF, FASN, and PTPN6, including compounds such as risperidone and sorafenib (Figure 8A, Table S7). To evaluate the binding potential between these prognostic genes and their predicted drugs, molecular docking was performed. The analysis revealed strong binding affinities between several pairs: FASN exhibited the strongest binding affinity for pretomanid, with a binding free energy of -8.9 kcal/mol. The interaction involved hydrogen bonds with amino acid residues LYS-193, THR-196, and GLU-115, with bond lengths of 3.0, 3.0, and 2.4 Å, respectively (Figure 8B, Table 1). PTPN6 demonstrated strong binding to tofacitinib, with a binding free energy of −8.2 kcal/mol. Hydrogen bonds were formed with THR-501 and GLN-504, with bond lengths of 3.3 and 3.0 Å, respectively (Figure 8C, Table 1). TNF showed robust binding to meropenem anhydrous, with a binding free energy of −7.1 kcal/mol. The interaction featured hydrogen bonds with the GLY-121 residue, with bond lengths of 2.9, 2.1, and 3.2 Å (Figure 8D, Table 1). Collectively, these molecular docking results, particularly the strong binding affinity observed between FASN and pretomanid, suggest that these identified drug-target pairs warrant further investigation as potential therapeutic strategies for CESC.

Figure 8 The drug-gene interaction network and molecular docking diagram of prognostic genes. (A) Interaction network of genes and drugs. The orange squares represent genes, while the blue squares represent drugs. (B) Molecular docking of pretomanid with FASN. (C) Molecular docking of tofacitinib with PTPN6. (D) Molecular docking of meropenem anhydrous with TNF. (B-D) Proteins are depicted using the cartoon model, while drugs are shown using the ball-and-stick model. The hydrogen bonds formed between the proteins and the drugs are colored yellow.

Identification of key cells

The initial scRNA-seq dataset comprised profiles of 84,457 cells and 32,248 genes (Figure S2A). Following rigorous quality control, a high-quality dataset of 67,964 cells and 32,248 genes was retained for downstream analysis (Figure S2B). The data were normalized, and the top 2,000 HVGs were selected to capture biological heterogeneity (Figure S2C). Dimensionality reduction was performed using the top 30 PCs, which were determined to be statistically significant. Unsupervised clustering based on these PCs revealed 19 distinct cell clusters (Figure S2D-S2G).

These clusters were annotated into five major cell types using canonical marker genes: epithelial cells, lymphocytes, endothelial cells, macrophages, and fibroblasts (Figure 9A). A dot plot visually confirms the expression patterns of these marker genes across the annotated cell types (Figure 9B). A bar plot further illustrates the cellular composition, showing that epithelial cells and lymphocytes constituted the largest proportions within the CESC samples (Figure 9C).

Figure 9 Identification of crucial cells. (A) Clustering diagram after cell annotation. (B) Dot plot of marker gene expression in each cell cluster after annotation. The horizontal axis represents marker genes, the vertical axis represents annotated cells, and each dot indicates the gene expression level within the cell. The darker the dot color, the higher the expression level, and the size of the dot reflects the proportion of the gene’s expression in the cell. (C) Proportion of cell clusters in normal and disease groups. (D) Expression differences of prognostic genes in different cell types. ns, not significant; **, P<0.01; ***, P<0.001. t-SNE, t-distributed stochastic neighbor embedding.

To identify the key cell type mediating the prognostic signature, we compared the expression of the four prognostic genes (TNF, PTPN6, FASN, TFRC) between CESC and normal samples within each cell type. This analysis revealed that macrophages exhibited significant differential expression of all four prognostic genes (Wilcoxon rank-sum test, P<0.05). Specifically, TNF, PTPN6, and TFRC were significantly upregulated in CESC-associated macrophages, whereas FASN was downregulated (Figure 9D). Notably, this downregulation of FASN in macrophages presents a contrast to its overall upregulation at the tissue level (see section “Expression level assessment and validation of 4 prognostic genes in CESC”), highlighting the cell-type-specific expression patterns within the TME. Based on this pronounced dysregulation, macrophages were therefore identified as the key cell type for subsequent cellular communication and subcluster analysis.

Cellular communication and pseudotime trajectory in CESC

Cell communication analysis revealed distinct interaction networks between CESC and control samples. In CESC, macrophages were the most interactive cell type, engaging in frequent and strong communications with endothelial cells, epithelial cells, and fibroblasts (Figure 10A). In contrast, within control samples, fibroblasts served as the central communicators, exhibiting prominent interactions with other fibroblasts, endothelial cells, and epithelial cells (Figure 10B). At the molecular level, the key outgoing signal from macrophages in CESC samples was the CXCL8-ACKR1 ligand-receptor pair targeting endothelial cells (Figure 10C). In control samples, fibroblasts primarily signaled to endothelial cells via the CCL2-ACKR1 pair (Figure 10D). These results underscore a shift in the cellular communication network in CESC, with macrophages emerging as a central signaling hub.

Figure 10 Analysis of cellular communication and macrophage pseudotemporal trajectories in the microenvironment of CESC. (A) Intercellular communication network in normal samples. (B) Intercellular communication network in disease samples. Nodes and colors represent different cell types. The size of the circular nodes reflects the number of cells in a given cell type. The thickness of the lines indicates the intensity of cell communication, with the line color corresponding to the cell type color. (C) Interaction pathways of each cell type in the normal group. (D) Interaction pathways of each cell type in the disease group. (E) Pseudotemporal analysis. The darker the color, the earlier the stage; the lighter the color, the later the stage. (F) Dimension reduction and clustering analysis. (G) Pseudotemporal trajectory analysis of different subgroups. Different colors represent different subgroups. (H) Expression of prognostic genes at different differentiation stages of key cells. Each dot represents a cell, with different colors denoting different stages. CESC, cervical cancer; t-SNE, t-distributed stochastic neighbor embedding.

Pseudotime trajectory analysis of macrophages inferred a continuous differentiation process (Figure 10E). Unsupervised clustering identified seven distinct macrophage subtypes (cluster 0–6) along this trajectory (Figure 10F). Pseudotime values were used to quantify differentiation progression, with darker-to-lighter color gradients indicating early-to-late states; clusters 3, 4, and 6 were annotated as early differentiation states, while clusters 0, 1, and 2 represented more advanced states (Figure 10G). Analysis of prognostic gene expression dynamics along pseudotime revealed distinct patterns: FASN expression remained relatively stable, TFRC expression progressively decreased, PTPN6 expression exhibited a biphasic pattern (initial increase followed by a decrease and stabilization), and TNF expression gradually declined and then plateaued (Figure 10H).

Collectively, the cell communication and pseudotime analysis results suggest that the prognostic genes TFRC, PTPN6, and TNF are associated with macrophage differentiation states and may contribute to the pathogenic role of macrophages in CESC.

Expression level assessment and validation of 4 prognostic genes in CESC

The Wilcoxon test confirmed the significant upregulation of the four prognostic genes in CESC samples from the TCGA-CESC dataset (P<0.05, Figure 11A). This differential expression pattern was further validated at the transcriptional level using RT-qPCR on independent clinical samples. Consistent with the bioinformatics analysis, the expression levels of TNF, FASN, TFRC, and PTPN6 were all significantly elevated in CESC tissues compared to matched normal controls (P<0.05, Figure 11B).

Figure 11 Evaluation and validation of the expression levels of four prognostic genes in cervical cancer. (A) The abscissa represents gene names, and the ordinate represents gene expression levels. The red color denotes the cervical cancer group, while the blue color represents the control group. (B) RT-qPCR detection of prognostic gene expression levels. The abscissa represents the groups, and the ordinate represents the relative gene expression level (gene/GAPDH). *, P<0.05; **, P<0.01; ****, P<0.0001. CESC, cervical cancer; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; RT-qPCR, reverse transcription quantitative polymerase chain reaction.

The consistent upregulation of these genes, coupled with their established association with patient prognosis, strongly supports their potential as viable therapeutic targets for CESC.

Discussion

MCD in CESC progression

CESC continues to be a major cause of cancer-related mortality among women worldwide. In this study, we developed a prognostic model underscoring the critical role of dysregulated MCD in CESC progression. MCD fosters an immunosuppressive TME by promoting angiogenesis and immune evasion (38). Our findings link the aberrant expression of MCD-RGs to high-risk disease phenotypes. This aligns with established mechanisms where TAMs drive immunosuppression via TNF secretion (39), and the TGF-β pathway, activated by factors like CDR1as, induces EMT (40). Furthermore, TGF-β polarizes MDSCs, which can collaborate with immature B cells to sustain immune evasion (41,42). Collectively, our results and these literature findings underscore the potential of targeting MCD to overcome immunotherapy resistance in CESC.

HPV infection and metabolic alterations in TME

While persistent high-risk HPV infection is the primary etiological factor for CESC (43), our multi-omics data reveal that additional metabolic and microenvironmental perturbations contribute to disease pathogenesis. We identified abnormalities in retinoic acid metabolism, which is known to alter myeloid differentiation via RARα/RXR heterodimers and inhibit anti-tumor M2 macrophage polarization (44,45). Additionally, the enrichment of synaptic organization pathways suggests the involvement of neuroimmune interactions, wherein nociceptive nerves modulate the TME via molecules such as NGF-TrkA and CGRP (46). These findings expand the understanding of non-canonical drivers of CESC progression beyond HPV alone.

Functional characterization of prognostic genes and model

Building upon the prognostic significance of MCD, we experimentally validated the upregulation of four key MCD-related genes in CESC. The pro-inflammatory cytokine TNF enhances CD8⁺ T cell activation but paradoxically sustains a chronic inflammatory TME conducive to tumor growth (47). PTPN6, encoding SHP-1, is a key regulator of JAK/STAT signaling, and its overexpression correlates with poor survival and macrophage differentiation (48-50). In our study, PTPN6 emerged as a key gene in the prognostic model, strongly correlating with myeloid markers, especially CD33 (MDSCs) and ITGAM/CSF1R (macrophages), suggesting a role in immunosuppressive myeloid recruitment/differentiation in CESC. Its downregulation in the HRG (Figure 4B), consistent with its inhibitory function, may reflect altered myeloid composition/function and impact survival, supporting its inclusion in a myeloid differentiation-based model. FASN promotes lymphangiogenesis and metastasis via lipid raft-dependent pathways and recruits M2 macrophages to aid immune evasion (51-54). Finally, TFRC mediates ferroptosis through HIF-1 and YAP/ACSL4 axes, highlighting its role as a therapeutic vulnerability (55,56). The upregulation of these genes, validated by RT-qPCR in our clinical samples, aligns consistently with TCGA data, reinforcing their collective role in CESC pathogenesis. It is noteworthy that FASN exhibited a cell-type-specific pattern of downregulation in TAMs despite overall tissue upregulation, suggesting a complex, context-dependent metabolic reprogramming in the CESC TME.

Collectively, these findings suggest that MCD-related genes may not only regulate tumor-intrinsic processes but also reshape the broader signaling landscape of the TME.

The significant enrichment of synapse organization-related pathways in the HRG suggests a potential involvement of neuro-immune interactions in CESC progression. Emerging evidence indicates that synaptic proteins and neurotransmitters within the TME can modulate the recruitment and polarization of myeloid cells, potentially steering them toward an immunosuppressive M2-like phenotype (57,58). The enrichment of these pathways in the HRG may reflect a more complex regulatory network where ‘neural-like’ signaling facilitates immune evasion and tumor stroma remodeling, ultimately contributing to the poor prognosis observed in high-risk patients (59,60).

Single-cell insights into macrophage dynamics

To further elucidate the expression patterns and functional context of these prognostic genes, we performed single-cell analysis, which revealed profound macrophage heterogeneity in CESC. Our cell communication analysis indicated that in the CESC TME, macrophages predominantly signal to endothelial cells via the CXCL8-ACKR1 axis, thereby fostering angiogenesis (61). In contrast, fibroblasts in normal tissue utilize CCL2-ACKR1 for homeostatic interactions (62). Pseudotime trajectory analysis uncovered dynamic expression of our prognostic genes during macrophage differentiation. Specifically, early TNF upregulation aligns with pro-inflammatory features, followed by a gradual decline toward later states, which may reflect a shift from anti-tumor activity to more immunosuppressive phenotypes (63). The progressive decrease of TFRC along pseudotime, with lower expression in later stages, may imply a reduced demand for iron metabolism as macrophages reach more differentiated states (64). PTPN6 exhibits an initial elevation along pseudotime, followed by stabilization at intermediate-to-late stages, suggesting a potential role in modulating the transition from inflammatory to immune-regulatory states (65,66). In contrast, the relatively stable expression of FASN across pseudotime indicates that lipid metabolic processes are maintained throughout macrophage differentiation, potentially supporting essential cellular functions (67). These findings posit that targeting specific macrophage subsets or their key communication pathways (e.g., CXCL8-ACKR1) could effectively disrupt pro-tumor niches.

Conclusions

In conclusion, our study delineates a prognostic signature derived from MCD patterns by integrating single-cell and bulk RNA sequencing data. Our results demonstrate that the established risk model serves as a reliable bioinformatic tool for predicting the survival of CESC patients and is closely associated with the tumor immune microenvironment.

The MCD-related prognostic model offers significant potential for clinical application by facilitating individualized risk stratification. In practice, this signature could assist clinicians in identifying high-risk patients who may benefit from intensified surveillance or specific immunotherapeutic interventions targeting the TGF-β pathway. Future prospective studies are warranted to validate the clinical utility of this model in diverse patient cohorts.

Acknowledgments

We would like to acknowledge that this study utilized data and resources from public databases including The Cancer Genome Atlas (TCGA), the Gene Expression Omnibus (GEO), the Functional Genomics Data (ArrayExpress), and the Molecular Signatures Database (MSigDB). We would also like to thank The First Affiliated Hospital of Guizhou University of Traditional Chinese Medicine for its support. Finally, we would like to extend our sincere gratitude to all participants and researchers who contributed to this study.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0174/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0174/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0174/prf

Funding: This study was supported by the Wang Qi's Inheritance Studio of Renowned Traditional Chinese Medicine in Guizhou Province (No. [2024] 38 of the Guizhou Traditional Chinese Medicine Administration).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2026-1-0174/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of The First Affiliated Hospital of Guizhou University of Traditional Chinese Medicine (No. KL2025-004), and all participants provided written informed consent.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Filho AM, Laversanne M, Ferlay J, et al. The GLOBOCAN 2022 cancer estimates: Data sources, methods, and a snapshot of the cancer burden worldwide. Int J Cancer 2025;156:1336-46. [Crossref] [PubMed]
3. Xu M, Cao C, Wu P, et al. Advances in cervical cancer: current insights and future directions. Cancer Commun (Lond) 2025;45:77-109. [Crossref] [PubMed]
4. Abrar SS, Hairon SM, Yaacob NM, et al. Global Research Trends in Cervical Cancer Survival in the Last Two Decades: A Bibliometric Analysis. Malays J Med Sci 2024;31:103-11. [Crossref] [PubMed]
5. Emagneneh T, Mulugeta C, Ejigu B, et al. Survival status of women with cervical cancer in Sub-Saharan Africa: a systematic review and meta-analysis, 2024. Front Oncol 2024;14:1491840. [Crossref] [PubMed]
6. Greene JT, Brian BF 4th, Senevirathne SE, et al. Regulation of myeloid-cell activation. Curr Opin Immunol 2021;73:34-42. [Crossref] [PubMed]
7. van Vlerken-Ysla L, Tyurina YY, Kagan VE, et al. Functional states of myeloid cells in cancer. Cancer Cell 2023;41:490-504. [Crossref] [PubMed]
8. Ding J, Hu H, Zhu Y, et al. Inhibiting CMTM4 reverses the immunosuppressive function of myeloid-derived suppressor cells and augments immunotherapy response in cervical cancer. J Immunother Cancer 2025;13:e011776. [Crossref] [PubMed]
9. Aizaz M, Khan A, Khan F, et al. The cross-talk between macrophages and tumor cells as a target for cancer treatment. Front Oncol 2023;13:1259034. [Crossref] [PubMed]
10. Chen Z, Zhao B. The role of tumor-associated macrophages in HPV induced cervical cancer. Front Immunol 2025;16:1586806. [Crossref] [PubMed]
11. Wang D, Han X, Liu HL. The role and research progress of tumor-associated macrophages in cervical cancer. Am J Cancer Res 2024;14:5999-6011. [Crossref] [PubMed]
12. Che Y, Yang Y, Suo J, et al. Induction of systemic immune responses and reversion of immunosuppression in the tumor microenvironment by a therapeutic vaccine for cervical cancer. Cancer Immunol Immunother 2020;69:2651-64. [Crossref] [PubMed]
13. Liu X, Kang X, Kang H, et al. The immunosuppressive role of MDSCs in HCC: mechanisms and therapeutic opportunities. Cell Commun Signal 2025;23:155. [Crossref] [PubMed]
14. He S, Zheng L, Qi C. Myeloid-derived suppressor cells (MDSCs) in the tumor microenvironment and their targeting in cancer therapy. Mol Cancer 2025;24:5. [Crossref] [PubMed]
15. Yi M, Li T, Niu M, et al. Exploiting innate immunity for cancer immunotherapy. Mol Cancer 2023;22:187. [Crossref] [PubMed]
16. Wang Y, Jia A, Bi Y, et al. Metabolic Regulation of Myeloid-Derived Suppressor Cell Function in Cancer. Cells 2020;9:1011. [Crossref] [PubMed]
17. Li Y, Gao X, Huang Y, et al. Tumor microenvironment promotes lymphatic metastasis of cervical cancer: its mechanisms and clinical implications. Front Oncol 2023;13:1114042. [Crossref] [PubMed]
18. Shi H, Zhong F, Yi X, et al. Application of an Autophagy-Related Gene Prognostic Risk Model Based on TCGA Database in Cervical Cancer. Front Genet 2020;11:616998. [Crossref] [PubMed]
19. Li J, Wan C, Li X, et al. Characterization of tumor microenvironment and tumor immunology based on the double-stranded RNA-binding protein related genes in cervical cancer. J Transl Med 2023;21:647. [Crossref] [PubMed]
20. Wu D, Liu Y, Liu J, et al. Myeloid cell differentiation-related gene signature for predicting clinical outcome, immune microenvironment, and treatment response in lung adenocarcinoma. Sci Rep 2024;14:17460. [Crossref] [PubMed]
21. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014;15:550. [Crossref] [PubMed]
22. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 2016;32:2847-9. [Crossref] [PubMed]
23. Gao CH, Yu G, Cai P. ggVennDiagram: An Intuitive, Easy-to-Use, and Highly Customizable R Package to Generate Venn Diagram. Front Genet 2021;12:706907. [Crossref] [PubMed]
24. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
25. Szklarczyk D, Morris JH, Cook H, et al. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res 2017;45:D362-8. [Crossref] [PubMed]
26. Zhang X, Zhao M, Chu T, et al. Comprehensive bioinformatics analysis of EXOSC family genes in head and neck squamous cell carcinoma. Sci Rep 2025;15:30361. [Crossref] [PubMed]
27. Zhou X, Qian Y, Ling C, et al. An integrated framework for prognosis prediction and drug response modeling in colorectal liver metastasis drug discovery. J Transl Med 2024;22:321. [Crossref] [PubMed]
28. Bull C, Byrne RM, Fisher NC, et al. Dual gene set enrichment analysis (dualGSEA); an R function that enables more robust biological discovery and pre-clinical model alignment from transcriptomics data. Sci Rep 2024;14:30202. [Crossref] [PubMed]
29. Xia MD, Yu RR, Chen DM. Identification of Hub Biomarkers and Immune-Related Pathways Participating in the Progression of Antineutrophil Cytoplasmic Antibody-Associated Glomerulonephritis. Front Immunol 2021;12:809325. [Crossref] [PubMed]
30. Chikaishi Y, Matsuoka H, Sugihara E, et al. Mutation Analysis of TMB-High Colorectal Cancer: Insights Into Molecular Pathways and Clinical Implications. Cancer Sci 2025;116:1082-93. [Crossref] [PubMed]
31. Zhang S, Yan J, Pan W, et al. Assessment of cancer-associated fibroblast signature genes in ovarian cancer patients: impact on immunity, drug resistance, and prognosis. Mol Cell Probes 2025;83:102038. [Crossref] [PubMed]
32. Stuart T, Butler A, Hoffman P, et al. Comprehensive integration of single-cell data. Cell 2019;177:1888-1902.e21. [Crossref] [PubMed]
33. Wolock SL, Lopez R, Klein AM. Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic Data. Cell Syst 2019;8:281-291.e9. [Crossref] [PubMed]
34. Li C, Guo L, Li S, et al. Single-cell transcriptomics reveals the landscape of intra-tumoral heterogeneity and transcriptional activities of ECs in CC. Mol Ther Nucleic Acids 2021;24:682-94. [Crossref] [PubMed]
35. Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun 2021;12:1088. [Crossref] [PubMed]
36. Wang H, Li H, Zou H, et al. Comprehensive single-cell sequencing reveals the tumor microenvironment and tumor-specific characteristics in trachea squamous cell carcinoma. Front Oncol 2025;15:1575647. [Crossref] [PubMed]
37. Veo B, Wang D, DeSisto J, et al. Single-cell multi-omics identifies metabolism-linked epigenetic reprogramming as a driver of therapy-resistant medulloblastoma. Nat Commun 2025;16:10470. [Crossref] [PubMed]
38. Goswami S, Anandhan S, Raychaudhuri D, et al. Myeloid cell-targeted therapies for solid tumours. Nat Rev Immunol 2023;23:106-20. [Crossref] [PubMed]
39. Verona F, Di Bella S, Schirano R, et al. Cancer stem cells and tumor-associated macrophages as mates in tumor progression: mechanisms of crosstalk and advanced bioinformatic tools to dissect their phenotypes and interaction. Front Immunol 2025;16:1529847. [Crossref] [PubMed]
40. Zhong G, Zhao Q, Chen Z, et al. TGF-β signaling promotes cervical cancer metastasis via CDR1as. Mol Cancer 2023;22:66. [Crossref] [PubMed]
41. Safavi P, Moghadam KB, Haghighi Z, et al. Interplay between LncRNA/miRNA and TGF-β Signaling in the Tumorigenesis of Gynecological Cancer. Curr Pharm Des 2024;30:352-61. [Crossref] [PubMed]
42. Jianyi D, Haili G, Bo Y, et al. Myeloid-derived suppressor cells cross-talk with B10 cells by BAFF/BAFF-R pathway to promote immunosuppression in cervical cancer. Cancer Immunol Immunother 2023;72:73-85. [Crossref] [PubMed]
43. Wang R, Zhang Y, Zhu Y, et al. Clinical Characteristics and Prognosis of Cervical Cancer Patients with Human Immunodeficiency Virus Infection: A Retrospective Study. Gynecol Obstet Invest 2022;87:324-32. [Crossref] [PubMed]
44. Liu J, Yang Z, Liu F, et al. Retinoic acid receptor gamma (RARγ) drives M2-like macrophage polarisation via CFI to promote thyroid cancer progression. Int Immunopharmacol 2026;173:116324. [Crossref] [PubMed]
45. Henning P, Westerlund A, Horkeby K, et al. Vitamin A enhanced periosteal osteoclastogenesis is associated with increased number of tissue-derived macrophages/osteoclast progenitors. J Biol Chem 2024;300:107308. [Crossref] [PubMed]
46. Zhang Y, Lin C, Liu Z, et al. Cancer cells co-opt nociceptive nerves to thrive in nutrient-poor environments and upon nutrient-starvation therapies. Cell Metab 2022;34:1999-2017.e10. [Crossref] [PubMed]
47. Zhou R, Xie Y, Wang Z, et al. Single-cell transcriptomic analysis reveals CD8 + T cell heterogeneity and identifies a prognostic signature in cervical cancer. BMC Cancer 2025;25:498. [Crossref] [PubMed]
48. Cui P, Lian J, Liu Y, et al. Pan-cancer analysis of the prognostic and immunological roles of SHP-1/ptpn6. Sci Rep 2024;14:23083. [Crossref] [PubMed]
49. Eswaran S, Adiga D, Khan G N, et al. Comprehensive analysis of the exocytosis pathway genes in cervical cancer. Am J Med Sci 2022;363:526-37. [Crossref] [PubMed]
50. Hu J, Wang S, Zhang X, et al. A genetic variant in the TAPBP gene enhances cervical cancer susceptibility by increasing m(6)A modification. Arch Toxicol 2024;98:3425-38. [Crossref] [PubMed]
51. Falchook G, Infante J, Arkenau HT, et al. First-in-human study of the safety, pharmacokinetics, and pharmacodynamics of first-in-class fatty acid synthase inhibitor TVB-2640 alone and with a taxane in advanced tumors. EClinicalMedicine 2021;34:100797. [Crossref] [PubMed]
52. Ping P, Li J, Lei H, et al. Fatty acid metabolism: A new therapeutic target for cervical cancer. Front Oncol 2023;13:1111778. [Crossref] [PubMed]
53. Du Q, Liu P, Zhang C, et al. FASN promotes lymph node metastasis in cervical cancer via cholesterol reprogramming and lymphangiogenesis. Cell Death Dis 2022;13:488. [Crossref] [PubMed]
54. Vanauberg D, Schulz C, Lefebvre T. Involvement of the pro-oncogenic enzyme fatty acid synthase in the hallmarks of cancer: a promising target in anti-cancer therapies. Oncogenesis 2023;12:16. [Crossref] [PubMed]
55. Xu X, Liu T, Wu J, et al. Transferrin receptor-involved HIF-1 signaling pathway in cervical cancer. Cancer Gene Ther 2019;26:356-65. [Crossref] [PubMed]
56. Zhao MY, Liu P, Sun C, et al. Propofol Augments Paclitaxel-Induced Cervical Cancer Cell Ferroptosis In Vitro. Front Pharmacol 2022;13:816432. [Crossref] [PubMed]
57. Zhang L, Jing W, Li Y, et al. Neurotransmitters: Key regulators of the tumor immune microenvironment. Semin Immunol 2026;81:102015. [Crossref] [PubMed]
58. Yang J, Wu Y, Lv X, et al. Neurotransmitters: an emerging target for therapeutic resistance to tumor immune checkpoint inhibitors. Mol Cancer 2025;24:216. [Crossref] [PubMed]
59. Yang J, Wu Y, Lv X, et al. Neurotransmitters: an emerging target for therapeutic resistance to tumor immune checkpoint inhibitors. Mol Cancer 2025;24:216. [Crossref] [PubMed]
60. Liang Y, Li H, Gan Y, et al. Shedding Light on the Role of Neurotransmitters in the Microenvironment of Pancreatic Cancer. Front Cell Dev Biol 2021;9:688953. [Crossref] [PubMed]
61. Du S, Qian J, Tan S, et al. Tumor cell-derived exosomes deliver TIE2 protein to macrophages to promote angiogenesis in cervical cancer. Cancer Lett 2022;529:168-79. [Crossref] [PubMed]
62. Guo F, Kong W, Li D, et al. M2-type tumor-associated macrophages upregulated PD-L1 expression in cervical cancer via the PI3K/AKT pathway. Eur J Med Res 2024;29:357. [Crossref] [PubMed]
63. Saeed AF. Tumor-Associated Macrophages: Polarization, Immunoregulation, and Immunotherapy. Cells. 2025;14:741. [Crossref] [PubMed]
64. Cong B, Zhang Q, Cao X. The function and regulation of TET2 in innate immunity and inflammation. Protein Cell 2021;12:165-73. [Crossref] [PubMed]
65. Hao F, Wang C, Sholy C, et al. Strategy for Leukemia Treatment Targeting SHP-1,2 and SHIP. Front Cell Dev Biol 2021;9:730400. [Crossref] [PubMed]
66. Zhou Z, Yao J, Wu D, et al. Type 2 cytokine signaling in macrophages protects from cellular senescence and organismal aging. Immunity 2024;57:513-527.e6. [Crossref] [PubMed]
67. Lim SA, Wei J, Nguyen TM, et al. Lipid signalling enforces functional specialization of T(reg) cells in tumours. Nature 2021;591:306-11. [Crossref] [PubMed]