Introduction

Cervical cancer is the second most common malignant tumor of the female reproductive system worldwide. Due to its large population, China accounted for 11.9% of global cervical cancer deaths and 12.3% of global cervical cancer DALYs in 2017.1

Definitive CRT including with external beam radiation therapy (EBRT) and brachytherapy (BT) boost is the current standard of care for cervical cancer patients with FIGO IIB-IVA disease.2

However, the systemic failure and pelvic failure were 17% and 6%, respectively, in retroEMBRACE3 in contrast to conventional BT studies, in which the systemic failure and pelvic failure were 9–12% and 16–17%, respectively.4,5 Checkpoint inhibitors have recently been approved for advanced and/or recurrent disease in combination with chemotherapy, and they are being studied on the front line in combination with chemoradiation.6

Literature reported that the expression profiles of ferroptosis-related genes (FRGs) are closely related to the tumor microenvironment (TME) and the prognostic survival of cervical squamous cell carcinoma patients.7 However, the mechanism is also unclear. Mechanistic modulation of ferroptosis in immune-combination radiotherapy in cervical cancer patients provides new ideas for individualized patient treatment.

In this study, we analyzed samples from radiotherapy-sensitive and insensitive CC patients and compared them with relevant data from the GEO database. We derived relevant prognostic genes through differential expression analysis, PPI network construction, one-way regression analysis, LASSO regression analysis, constructed risk models and column-line diagrams, and performed immune infiltration analysis, transcription factor prediction and ceRNA regulatory network analysis, single-cell analysis, aiming to provide new ideas for guiding the prognosis and treatment of CC patients, network analysis, and single-cell analysis, aiming to provide new ideas to guide the prognosis and treatment of CC patients.

Materials and Methods Data and Sample Source

A total of 23 samples were collected for sequencing, including 12 radiation-sensitive samples from cervical cancer (CC) patients and 11 radiation-resistant samples from CC patients. IRB approved Identifier from the Ethics Committee of Shandong First Medical University Affiliated Cancer Hospital (Shandong Academy of Medical Sciences). Ethical/Copyright Corrections: Our study complies with the Declaration of Helsinki.

Our research team prospectively conducted this research at Shandong First Medical University Affiliated Cancer Hospital, China. The inclusion criteria included: (1) Patients 18 years or older; (2) patients with histologically confirmed cervical cancer received radiation therapy. (3) Received ≥1 cycle of systemic chemotherapy. (4) Eastern Cooperative Oncology Group performance status of 0–1. (5) Adequate hematologic, hepatic, and kidney function profile. (6) Urine or serum pregnancy test within 3 days prior to the first dose in female subjects of childbearing potential. (7) Subjects were willing and able to comply with the schedule of visits, treatment protocols, laboratory tests, and compliance with other requirements of the study.

Exclusion criteria included: (1) Concurrent enrollment in another clinical study, unless it is an observational, non-interventional clinical study or a follow-up period of an interventional study. (2) Known active tuberculosis (TB): Subjects suspected of having active TB need to be excluded by clinical examination. (3) Known active syphilis infection. (4) Known history of mental illness, drug dependence, alcoholism or drug addiction. (5) Any pre-existing or current medical condition, treatment, or laboratory test abnormality that may confound the results of the study, interfere with the subject’s ability to participate in the study in its entirety, or participation in the study may not be in the subject’s best interest. (6) Localized or systemic disease not due to malignancy or disease or condition secondary to a tumor that may result in a higher medical risk and/or uncertainty in the evaluation of survival, such as tumor-like leukemic reaction (white blood cell count >20 × 109/L), malignant manifestations (eg, known weight loss of more than 10% in the 3 months prior to Screening). (7) Be pregnant or breastfeeding, or plan to breastfeed during the study. (8) other conditions that the investigator considers inappropriate for enrollment. (9) All methods were performed in accordance with the relevant guidelines and regulations.

The Cancer Genome Atlas (TCGA) database offered the transcriptome data, clinical details, and survival information for cervical squamous cell carcinoma and adenocarcinoma (TCGA-CESC). A total of 306 CC samples (cervical tissue) were included in the TCGA-CESC, which was served as a validation set. The GSE9750 and GSE44001 datasets were acquired from the Gene Expression Omnibus (GEO) database. The GSE44001 dataset contained 300 CC samples (cervical tissue) including mRNA expression profile and survival information and was designated as a training set, with GPL14951 sequencing platform employed. The GSE9750 dataset including 33 CC samples (cervical tissue, GPL96) and 24 normal samples was used for variance analysis. We obtained 416 ferroptosis related genes (FRGs) from the literature.7 A total of 1,793 immune related genes (IRGs) were acquired from the immunology database and analysis portal (ImmPort) database.

Identification of Differentially Expressed Genes (DEGs) in CC Radiotherapy and Functional Annotation Analysis

The radiation-DEGs (DEGs1) were detected based on RNA sequencing data of radiation-sensitive samples by employing the DESeq2 package (version 1.41.10) (|Log2FC| > 1 and p < 0.05).8 The Limma package (version 3.56.2) was utilized to obtain CC-DEGs (DEGs2) within GSE9750 (|Log2FC| > 1 and p < 0.05).9–11 The DEGs related to CC radiotherapy were identified by the intersection between DEGs1 and DEGs2 utilizing the R package ggVennDiagram (version 1.2.3).12

Utilizing clusterProfiler (version 4.8.2), Gene Ontology (GO), as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed to examine the functional annotation of the DEGs related to CC radiotherapy (FDR < 0.05).13

Identification of FRGs-IRGs-DEGs (FIGs) and Protein–Protein Interactions (PPI) Network

FRG-DEGs were identified by intersection between DEGs related to CC radiotherapy and the FRGs utilizing the R package VennDiagram (version 1.7.1).14 Similarly, the identification of IRG-DEGs was achieved by the intersection between DEGs related to CC radiotherapy and the IRGs. The correlation analysis was ultimately conducted between FRG-DEGs and IRG-DEGs. Genes exhibiting a correlation coefficient exceeding 0.4 were designated as FIGs. The clusterProfiler package (version 4.8.2) was utilized to conduct GO and KEGG pathway enrichment analysis of FIGs (FDR < 0.05).13 The interactions among proteins corresponding FIGs were predicted (with a confidence score threshold of 0.4) by utilizing the STRING database, and the PPI network was visualized using Cytoscape (version 3.8.2).15

Univariate Cox Analysis and Least Absolute Shrinkage and Selection Operator (LASSO) Analysis

Univariate Cox analysis was conducted on candidate genes with survival data in GSE44001 dataset (p < 0.05), and candidate genes related to prognosis were screened. Second, utilizing the glmnet package (version 4.1–8), LASSO analysis was conducted on candidate genes related to prognosis to identify prognostic genes.16 The box plot was employed to illustrate the distinction expression of prognostic genes with normal and CC samples in GSE9750 dataset.

Construction and Validation of Risk Model

In order to develop prognostic risk models, the expressions of prognostic genes were assessed, and regression coefficients were obtained from LASSO analysis. Subsequently, a risk score was computed utilizing the subsequent formula.

The risk coefficient of each gene is denoted by the coef, while the expr indicates the level of gene expression. Meanwhile, the risk scores for 300 CC samples were calculated within GSE44001 dataset. Furthermore, a Kaplan–Meier (K-M) survival analysis was completed to observe the difference in prognosis between the different risk groups.17,18 The time ROC (version 0.4) package was employed to generate receiver operating characteristic curves (ROC) for disease-free survival (DFS) times of 1, 3, and 5 years.

Creation and Evaluation of Nomogram

Univariate and multivariate Cox regression analyses were performed for prognostic genes in the GSE44001 dataset to identify independent prognostic factors. The rms software (version 6.7–0) was employed to establish a nomogram utilizing independent prognostic factors, enabling the prediction of survival rates.19 The accuracy and reliability of the nomogram were assessed by calibration curves and decision curve analysis (DCA).

Gene Set Variation Analysis (GSVA) and Gene Set Enrichment Analysis (GSEA)

To further explore the relevant pathways connecting groups of high and low risk in GSE44001, the GSVA analysis was performed using the R package GSVA (version 1.49.4).20

In the RNA sequencing data set and the GSE44001 dataset, the correlation coefficients of expression amount between all genes and each prognostic gene were calculated as the sequencing standard.13

Function Similarity of Prognostic Genes and GeneMANIA Network

To assess the function similarity of the prognostic genes, the semantic analysis method offered in GOSemSim (version 2.27.2) was employed.21 The GeneMANIA database was utilized to predict genes and biological functions associated with prognostic genes.

Immune Microenvironment Analysis in Different Risk Groups

The relative fraction of 22 different types of immune cells in the samples of GSE44001 was evaluated employing the CIBERSORT method.

The ESTIMATE algorithm was utilized to calculate the immune score, stromal score and ESTIMATE score of CC patients in GSE44001.

Immunotherapy Analysis

The literature provided 41 immunochemokines, 57 immune cell marker genes, as well as both immunosuppressive and immunostimulatory factors.22–24 The tumor immune dysfunction and exclusion (TIDE) scores were computed between two risk groups during the GSE44001 dataset.

Construction of Regulatory Networks, as Well as Drug Prediction Analysis

The NCBI, JASPAR, and UCSC databases were employed to predict TFs and its potential binding sites of prognostic genes. The miRDB and miRWalk were used to predict miRNAs of prognostic genes. The target drugs of prognostic genes were predicted in the Drug-Gene Interaction Database and DrugBank database to conduct drug-gene network.

Single-Cell RNA-Seq Analysis

First, quality control was performed on the single-cell dataset GSE208653 using the Seurat package (version 5.0.1) to filter out cells (min.cells = 3, min.features = 200).25–27 Secondly, FindVariableFeatures function vst method was used to extract genes with large inter-cell variation coefficients.28 Finally, to evaluate the relationship between prognostic genes and cell types, the expression of prognostic genes in the different cell types annotated was further analyzed.

PCR Analysis

Tissue RNA was extracted, and 1ul of RNA was taken for RNA concentration detection using NanoPhotometer N50. Reverse transcription analysis of mRNA was performed using SweScript First Strand cDNA synthesis kit.

Immunohistochemistry

Pathological indicators were detected by pathologists through immunohistochemistry and all slides were assessed in agreement by two observers, who were unaware of the patients’ clinical data.

Statistical Analysis

The R software (version 4.2.1) was utilized to process and analyze the data. The p value <0.05 was considered statistically significant.

Results Identification of DEGs of Radio‑Sensitive Samples in the CC Cohort

A total of 1,782 DEGs1 were identified from RNA sequencing data when comparing radiation-sensitive samples and radiation-resistant samples. Among these DEGs1, 925 were found to be upregulated, while 857 were downregulated (Figure 1A and B). There were 3,905 DEGs2 obtained from the GSE9750 dataset, with 2,045 upregulated DEGs2 and 1,860 downregulated DEGs2 (Figure 1C and D). The 329 DEGs related to CC radiotherapy were identified by the intersection between DEGs1 and DEGs2 (Figure 1E). A total of 188 entries were enriched by GO, including 8 CCs (Cellular Component), 18 MFs (Molecular Function) and 162 BPs (Biological Process). The first ten entries of MFs and BPs and all entries of CCs were selected for GO enrichment. GO enrichment analysis indicated that the functions of DEGs related to CC radiotherapy were substantially enriched in heme binding, integrin binding, and endoderm formation, etc. The KEGG pathways were breast cancer, MAPK signaling pathway PI3K-AKT and tryptophan metabolism in cancer (Figure 1F and G).

Figure 1 (A–G) Identification of DEGs of radio‑sensitive samples in the CC cohort. (A) DEGs1 Volcano Map. (B) DEGs1 heat map. (C) DEGs2 Volcano Map. (D) DEGs2 density heatmap. (E) Venn diagram of intersection between DEGs1 and DEGs2. (F) DEGs GO enrichment map. (G) DEGs KEGG enrichment map.

Differences FRGs-IRGs Identification Associated with Ferroptosis and Tumor Immune Pathways

A total of 15 FRG-DEGs and 52 IRG-DEGs were identified by intersecting DEGs related to CC radiotherapy with the FRGs and the IRGs, respectively (Figure 2A). The correlation analysis of 15 FRG-DEGs and 52 IRG-DEGs identified 62 FIGs with a correlation coefficient greater than 0.4 (Figure 2B).

Figure 2 (A–D) Differences FRGs-IRGs identification. (A) Venn plot of the intersection of candidate gene 1 (left) and candidate gene 2 (right). (B) Heat map of correlation between candidate gene 1 and candidate gene 2. Boxes with *indicate correlations greater than 0.4. (C) FRGs IRGs GO enrichment map. (D) FRGs IRGs KEGG enrichment map.

The functional enrichment analysis of 62 FIGs indicated that the GO functions were significantly enriched in negative regulation of cell adhesion, ERK1 and ERK2 cascade, receptor ligand activity, growth factor activity, etc. (Figure 2C). KEGG pathways were enriched in Coronavirus disease – COVID-19, MAPK signaling pathway, endocrine resistance, toll-like receptor signaling pathway, etc. (Figure 2D).

CALCRL, UCHL1, GNRH1, ACVRL1, and MUC1 Were Identified as Prognostic Genes

A total of 53 candidate genes exhibited interactions with each other via PPI network for subsequent analysis (Figure 3A). After analyzing 53 candidate genes, six candidate prognosis genes were identified in the GSE44001 dataset utilizing univariate Cox regression analysis (p < 0.05) (Figure 3B). Then, LSAAO analysis was utilized to identify five prognostic genes (CALCRL, UCHL1, GNRH1, ACVRL1, and MUC1) from six candidate prognosis genes and construct a risk model (Figure 3C and D). The box plots illustrated that the expression trend of five prognostic genes was consistent in the GSE9750 dataset and RNA sequencing dataset. The CALCRL and ACVRL1 were up-regulated in CC samples and radiation-sensitive samples, whereas GNRH1 UCHL1, and MUC1 were down-regulated (Figure 3E). The GNRH1, MUC1 and UCHL1 were up-regulated in CC samples compared to control group, whereas CALCRL and ACVRL1 were down-regulated compared to control group (Figure 3F).

Figure 3 (A–F) CALCRL, UCHL1, GNRH1, ACVRL1, and MUC1 were identified as prognostic genes. (A) FRGs IRGs PPI network. (B) Single factor Cox forest map. (C and D) LASSO regression analysis for screening characteristic genes. (E and F) Expression of key genes in CC samples. *p<0.05; **p<0.01; ***p<0.001; ns p>0.05.

The Risk Model Demonstrated Favorable Effectiveness in Predicting Outcomes

The risk scores of the patient samples in the dataset GSE44001 were calculated, and the surv_cutpoint function was used to calculate the optimal threshold for the risk scores (0.9262071), which categorized the patients into high-risk and low-risk groups. The scatter plot demonstrated a significant increase in fatality count as the risk scores escalated within the samples. The K-M (Figure 4A and C) curves demonstrated that patients with CC in the high-risk group exhibited a worse survival duration. The risk model (Figure 4B and D) demonstrated favorable effectiveness in predicting outcomes at 1, 3, and 5 years, as evidenced by ROC (Figure 4E and F) curves with AUC values surpassing 0.6. The effectiveness of the risk model was evaluated by the analyses conducted on GSE44001 and validated using TCGA-CESC datasets.

Figure 4 (A–D) The risk model demonstrated favorable effectiveness in predicting outcomes. (A) K-M curve (left), (B) survival curve (middle) and (E) ROC (right) of training set GSE44001. (C) K-M curve (left), (D) survival curve (middle) and (F) ROC (right) of the validation set TCGA-CESC dataset.

Nomogram Exhibited an Accurate Predictive Capability

Univariate and multivariate Cox regression analysis demonstrated that CALCRL, GNRH1, and MUC1 were independent prognostic factors (Figure 5A and B). A nomogram was constructed based on independent prognostic factors (Figure 5C). The survival predictions of 1, 3 and 5 years were closely aligned with the theoretical straight line, which indicated a good predictive performance of the nomogram (Figure 5D). In addition, the nomogram demonstrated greater advantages compared to individual genes, as evident from the three- and five-year DCA curves (Figure 5E).

Figure 5 (A–E) Nomogram exhibited an accurate predictive capability. (A and B) Univariate and multivariate cox regression analysis. (C) Based on independent prognostic factors, a column chart was constructed in the dataset GSE44001 to predict the 1, 3, and 5-year survival rates of CC patients. (D) Nomogram model calibration curve. (E) Nomogram model DCA curve for CC patients in years 1, 3, and 5.

Prognostic Genes Were Enriched in Immune Related Pathways

Based on the sorting by p values, the top 20 pathways of GSVA analysis were demonstrated, such as nagy staga components human, reactome translation of replicase and assembly of the replication transcription complex, thillainadesan znf217 targets up, etc. (Figure 6A).

Figure 6 Continued.

Figure 6 (A–G) Prognostic genes were enriched in immune related pathways. (A) GSVA results are hot. (B and C) Results of GNRH1 enrichment in CC radiation samples; (D and E) Results of GNRH1 enrichment in GSE44001. (F) Box plot of functional similarity of biomarkers. (G) GeneMANIA database result network diagram.

In RNA sequencing data set, the GSEA enrichment analysis results of prognostic genes revealed significant enrichments in immune related functions and pathways, including spliceosome, cell cycle, lysosome, extracellular matrix structural constituent, collagen fibril organization, immunoglobulin complex, etc. (Figure 6B and C, Supplementary Figures 1–8).

In the GSE44001 dataset, prognostic genes were significantly enriched in ferroptosis, Fanconi anemia pathway, gap junction, focal adhesion, prion disease, T cell receptor signaling pathway, etc. (Figure 6D and E, Supplementary Figures 9–16).

GNRH1 Might Be a Relatively Important Gene

The results of functional similarity analysis showed that CALCRL, UCHL1, ACVRL1 and MUC1 had high average functional similarity, while GNRH1 had significantly lower functional similarity than the other four genes, which might be a relatively important gene (Figure 6F). There were 20 genes (RAMP1, GNRH2, PAMP2, etc.) related to the function of five prognostic genes, among which ACVRL1 was related to the regulation of pathway-restricted SMAD protein phosphorylation, pathway-restricted SMAD protein phosphorylation, etc. (Figure 6G).

The Tumor Microenvironment Exhibited Notable Variances in Two Risk Groups

Heat maps showed differences in the immune microenvironment between the high and low-risk groups (Figure 7A). In two risk groups, there were four immune cell categories (activated natural killer (NK) cells, activated memory CD4+ T cells, naive CD4+ T cells, T follicular helper cells) and there existed notable disparities (Figure 7B). The Spearman correlation analysis revealed a significant negative association between activated NK cells with activated memory CD4+ T cells (R=−0.43) and naive CD4+ T cells (R=−0.18), as well as T follicular helper cells with activated memory CD4+ T cells (R=−0.18) (Figure 7C). Additionally, there was a significant negative association observed between ACVRL1 with naive CD4+ T cells (R=−0.21) and a significant positive association with activated memory CD4+ T cells (R=0.29) (Figure 7D). Furthermore, two risk groups exhibited notable variances in ESTIMATE score, indicating a significant difference in the tumor microenvironment between the two groups (Supplementary Figure 17). The risk score was negatively correlated with immune score, stromal score and ESTIMATE score (Supplementary Figure 18).

Figure 7 (A–H) The tumor microenvironment exhibited notable variances in two risk groups. (A) Heat map of immune cell infiltration. (B) Box plot of differences in immune cells between high and low-risk groups in GSE44001. (C) Differential immune infiltration cell correlation heatmap. (D) Biomarkers and differential immune cell correlation heatmaps. (E–H) Heat map of the correlation between biomarkers and different immune regulatory factors. * denotes the difference in immune cells in the high and low risk groups, the more * the more significant the difference: *p<0.05; **p<0.01; ***p<0.001; ns p>0.05.

There Was a Significant Negative Correlation Between MUC1 and GNRH1 with Most of the Immunomodulators, Immune Cell Marker Genes, and Immunochemokines

The correlation analysis of prognostic genes and immunomodulators revealed a significant negative correlation between MUC1 and GNRH1 with most of the immunomodulators, while CALCRL and ACVRL1 exhibited a significant positive correlation with most of the immunomodulators (Figure 7E). At the same time, the risk scores were found to exhibit a significant positive correlation with the majority of immunosuppressive factors while displaying a negative correlation with most immunostimulants (Figure 7F). The correlation between prognostic genes and immune cell marker genes showed the same trend. There were significant negative correlations between MUC1 and GNRH1 with most immune cell marker genes, as well as significant positive correlations between CALCRL and ACVRL1 with most immune cell marker genes. The risk scores were significantly positively correlated with CD68, CEACAM8, HLA−DRA, HLA−DPA1 and HLA−DPB1, as well as negatively correlated with CDR8, FOXP3, NRP1 and TGFB1 (Figure 7G). In the same way, MUC1 and GNRH1 showed significant negative correlations with most immunochemokines, while CALCRL and ACVRL1 exhibited significant positive correlations with most immunochemokines. Risk scores were positively correlated with most immunochemokines (CCL11, CCL15, CXCL17, etc.) and negatively correlated with CCL14, CCL19, CCL5 and CXCL16 (Figure 7H).

ACVRL1 and CALCRL Could Be Crucial Target Genes in ceRNA Regulatory Network

There was a total of 21 TFs associated with prognostic genes. The TFs predicted by ACVRL1 were Prdm5, SPIB, Arx, and PATZ1. CALCRL predicted Gfi1B and ZNF345A as TFs. GNRH1 predicted DMRTC2, ZNF460, ZNF384, PRARD, RXRG, IKZF1, HNF4A, and HNF4G as TFs. MUC1 included ZNF331, Tcf21, Foxq11, INSM1, and ZBED2 as predicted TFs. UCHL1 predicted ZNF24 and ZNF460 as TFs (Supplementary Figures 19–23). The prediction of ACVRL1 identified 18 potential binding sites for TFs, while the prediction of MUC1 revealed five potential binding sites and the prediction of GNRH1 indicated 24 potential binding sites, as well as the prediction of UCHL1 identified four potential binding sites. These results were then visualized using Cytoscape software to construct a network of TF-binding sites (Figure 8A–D). There were 52 miRNAs associated with ACVRL1, 20 miRNAs associated with CALCRL, as well as three miRNAs associated with MUC1 and UCHL1, respectively (Figure 8E, Supplementary Figures 24–27). Moreover, the miRNA corresponding to ACVRL1 predicted 668 lncRNAs, and the miRNA corresponding to CALCRL predicted 144 lncRNAs, indicating that they could be crucial target genes for CC prognosis (Figure 8F).

Figure 8 Continued.

Figure 8 (A–G) ceRNA regulatory network. (A) ACVRL1 transcription factor binding site network diagram. (B) MUC1 transcription factor binding site network diagram. (C) GNRH1 transcription factor binding site network diagram. (D) UCHL1 transcription factor binding site network diagram. (E) MiRNA biomarker network diagram. (F) LncRNA miRNA mRNA network diagram. (G) Biomarker drug network diagram.

In the DGI and DrugBank databases, a total of 35 drugs were predicted by five prognostic genes, and the drugs predicted by the CACRL were found in both databases: Olcegepant, Telcagepant, Erenumab, Ubrogepant and Rimegepant (Figure 8G).

The Expression of ACVRL1 Was Lower in Fibroblasts in CC Samples

After filtering the GSE208653 dataset, there were 5,031 cells remaining. Then, 2,000 highly variable genes were identified with the top 10 genes were demonstrated, and 10 PCs were selected for subsequent analysis (Supplementary Figures 28 and 29). Subsequently, all the cells were grouped into 10 clusters and annotated (Figure 9A). They were natural killer (NK) T cells, epithelial cells, neutrophils, myeloid cells, smooth muscle cells, plasma cells, fibroblasts, mast cells, B cells, and endothelial cells. The expression of these marker genes in each cell was displayed using a dot matrix (Figure 9B). The prognostic gene MUC1 was expressed in epithelial cells, and ACVRL1 was more highly expressed in fibroblasts in normal samples (Figure 9C, Supplementary Figure 30).

Figure 9 (A–C) The tumor microenvironment. (A) UMAP clustering analysis. (B) Expression of marker genes for different cell types. (C) The expression of biomarkers in different cell types.

PCR and IHC Validation Results

Our group performed PCR of tissue samples to verify the above results. The results of PCR showed that the expression of UCHL1, GNRH1 and MUC1 was up-regulated in tumor tissues compared to normal cervical tissues, while the expression of ACVRL1 was down-regulated. The expression of CALCRL was not statistically significant. The expression of ACVRL1 was down-regulated, while the expression of CALCRL was down-regulated, but did not reach statistical significance (Figure 10).

Figure 10 PCR detection of the expression of the above 5 indicators in CC tissue and adjacent normal tissue. *p<0.05; **p<0.01; ***p<0.001; ns p>0.05.

At the same time, IHC for the above indicators were performed (Figure 11). The results of the IHC and the PCR results generally followed the same trend. The expression of UCHL1, GNRH1 and MUC1 was higher in tumor than normal tissue. The expression of ACVRL1 and CALCRL was similar in both tumor and normal tissue.

Figure 11 IHC of the expression of the above 5 indicators in CC tissue and adjacent normal tissue.

Discussion

The biological behavior of cervical cancer is complex, and the individualization of radiotherapy varies widely. Our group detected tumor and normal tissue samples before and after radiotherapy in cervical cancer patients with different sensitivities to radiotherapy by sequencing, compared the differences and analyzed them in comparison with the relevant data from databases, such as GEO, and obtained five prognosis-related genes in cervical cancer: CALCRL, UCHL1, GNRH1, ACVRL1, and MUC1. The flowchart of our study is shown in Figure 12.

Figure 12 Trail profile. The blue color represents the database data and analysis, the red color represents the data analysis of cervical cancer patients’ tissues, and the yellow line indicates the interactive analysis of the two and the validation of the results.

CALCRL was reported predominantly expressed in thyroid carcinomas, parathyroid adenomas, small-cell lung cancers, large-cell neuroendocrine carcinomas of the lung, pancreatic neuroendocrine neoplasms, renal clear-cell carcinomas, pheochromocytomas, lymphomas, and melanomas. The strong expression of CALCRL may represent a useful target for future research.29–31

UCHL1, as a marker of neurodegenerative diseases specifically expressed in the brain and testis, is also associated with the occurrence of tumors.32 In lung cancer cells, UCHL1 promoted PD-L1 expression and inhibition of UCHL1 might suppress immune escape of NSCLC through downregulation of PD-L1 expression in NSCLC cells.33 Therefore, based on the previous studies on UCHL1 as well as the findings of our group, it is suggested by the results of the present study for the first time that UCHL1 may act as an important mechanistic factor in radiotherapy combined with immunotherapy in patients with cervical cancer.

It is well established that GnRH1 and its receptor are expressed in cancer cells derived from reproductive tissues and administration of GnRH1 analogs inhibits their proliferation.34,35 Our study shows for the first time that the differential expression of GNRH1, a prognostic gene, in CC cervical cancer patients is regulated by ferroptosis and immunity. The low functional similarity between GNRH1 and the remaining four genes, CALCRL, UCHL1, ACVRL1 and MUC1, suggests that GNRH1 may be a focus for future studies.

MUC1 is significantly elevated in tumor tissue expression in patients with reproductive system tumors and is an important predictor of diagnosis and prognosis in patients with cervical cancer.36–38 Enhancement of MUC1-specific immune response by co-administration of liposomal DDA/MPLA and lipoglycopeptide has led to new breakthroughs in immunovaccine research for epithelial cancers.39

ACVRL1 is closely related to the TGF-β pathway40 and has been reported in studies related to chemoresistance in colorectal cancer,41 but no clear prognostic or other factors have been reported in patients with cervical cancer. Our study is the first to report that ACVRL1 shows a positive correlation with immune cell marker genes and immunochemokines in cervical cancer patients. However, in our research, the PCR results indicated that ACVRL1 expression is a little lower in CC than in normal tissue. The function and mechanism of ACVRL1 still need to be further evaluated.

The results also showed a significant positive correlation between activated memory CD4 T cells and T cell follicular helper cells and activated NK cells, as well as a significant positive correlation between activated NK cells and naïve CD4 T cells.

The results of this study are different from previous genomic results on ferroptosis prognostic prediction in cervical cancer.7 The main reason for this is that our patient group was selected for factors that are sensitive and resistant to radiotherapy.

In this study, we obtained five ferroptosis and immune-related prognostic genes in cervical cancer, and constructed the ferroptosis and immune-related prognostic risk model for the first time in cervical cancer, which can help to further improve the prognostic evaluation system of cervical cancer patients, and the results of immunoassay and drug prediction analysis can provide references to the treatment of cervical cancer, and the five genes in the model are also expected to be the potential targets for cervical cancer treatment, which can provide new ideas for the guidance of prognosis and treatment of cervical cancer patients. The five genes in the model are also expected to be potential targets for cervical cancer treatment, which can provide new ideas to guide the prognosis and treatment of cervical cancer patients.

Data Sharing Statement

The datasets generated and/or analyzed during the current study are not publicly available now but are available from the corresponding author on reasonable request.

Informed Consent

All participants gave written informed consent.

Approval of the research protocol by an Institutional Reviewer Board:

The study protocol was approved by the Research Ethics Committee of Shandong First Medical University Affiliated Cancer Hospital (Shandong Academy of Medical Sciences) (Approval No. 2022003063).

Acknowledgments

Thanks for Dr. Xin Wang and all nurses’ groups for the help and take care of all patients.

Funding

The project was supported by the China Postdoctoral Science Foundation (Grant No. 2022M721446), Shandong Provincial Natural Science Foundation (Grant No. ZR2022MH288), Shandong Postdoctoral Science Foundation (Grant No. SDCX-ZG-202303064), Wu Jieping Medical Foundation (320.6750.2023-17-26), Xinrui Cancer Foundation (FB2023031673). Shandong Medical Association Clinical Research Foundation (Grant No. YXH2022ZX02151), Bethune Translational Research Foundation of Radiation Oncology (Grant No. flzh202125), Healthy Science Development Plan of Shandong, China (Grant No. 202009030673), Shandong First Medical University Research Project on Education and Teaching (XM2022163), Medical and Health Science and Technology Development Project of Shandong Province (202315021029).

Disclosure

All authors declare no conflict of interest.

References

1. Guo M, Xu J, Du J. Trends in cervical cancer mortality in China from 1989 to 2018: an age-period-cohort study and Joinpoint analysis. BMC Public Health. 2021;21(1):1329. doi:10.1186/s12889-021-11401-8

2. Musunuru HB, Pifer PM, Mohindra P, Albuquerque K, Beriwal S. Advances in management of locally advanced cervical cancer. Indian J Med Res. 2021;154(2):248–261. doi:10.4103/ijmr.IJMR_1047_20

3. Tan LT, Pötter R, Sturdza A, et al. Change in patterns of failure after image-guided brachytherapy for cervical cancer: analysis from the RetroEMBRACE Study. Int J Radiat Oncol Biol Phys. 2019;104(4):895–902. doi:10.1016/j.ijrobp.2019.03.038

4. Chemoradiotherapy for Cervical Cancer Meta-Analysis Collaboration. Reducing uncertainties about the effects of chemoradiotherapy for cervical cancer: a systematic review and meta-analysis of individual patient data from 18 randomized trials. J Clin Oncol. 2008;26(35):5802–5812. doi:10.1200/JCO.2008.16.4368

5. Pötter R, Tanderup K, Schmid MP, et al. MRI-guided adaptive brachytherapy in locally advanced cervical cancer (EMBRACE-I): a multicentre prospective cohort study. Lancet Oncol. 2021;22(4):538–547. doi:10.1016/S1470-2045(20)30753-1

6. Podwika SE, Duska LR. Top advances of the year: cervical cancer. Cancer. 2023;129(5):657–663. doi:10.1002/cncr.34617

7. Han S, Wang S, Lv X, Li D, Feng Y. Ferroptosis-related genes in cervical cancer as biomarkers for predicting the prognosis of gynecological tumors. Front Mol Biosci. 2023;10:1188027. doi:10.3389/fmolb.2023.1188027

8. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. doi:10.1186/s13059-014-0550-8

9. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47. doi:10.1093/nar/gkv007

10. Gustavsson EK, Zhang D, Reynolds RH, Garcia-Ruiz S, Ryten M. ggtranscript: an R package for the visualization and interpretation of transcript isoforms using ggplot2. Bioinformatics. 2022;38(15):3844–3846. doi:10.1093/bioinformatics/btac409

11. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32(18):2847–2849. doi:10.1093/bioinformatics/btw313

12. 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. doi:10.3389/fgene.2021.706907

13. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012;16(5):284–287. doi:10.1089/omi.2011.0118

14. Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinf. 2011;12:35. doi:10.1186/1471-2105-12-35

15. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–2504. doi:10.1101/gr.1239303

16. Engebretsen S, Bohlin J. Statistical predictions with glmnet. Clin Clin Epigenet. 2019;11(1):123. doi:10.1186/s13148-019-0730-1

17. Lei J, Qu T, Cha L, et al. Clinicopathological characteristics of pheochromocytoma/paraganglioma and screening of prognostic markers. J Surg Oncol. 2023;128(4):510–518. doi:10.1002/jso.27358

18. Liu TT, Li R, Huo C, et al. Identification of CDK2-related immune forecast model and ceRNA in lung adenocarcinoma, a pan-cancer analysis. Front Cell Dev Biol. 2021;9:682002. doi:10.3389/fcell.2021.682002

19. Xu J, Yang T, Wu F, Chen T, Wang A, Hou S. A nomogram for predicting prognosis of patients with cervical cerclage. Heliyon. 2023;9(11):e21147. doi:10.1016/j.heliyon.2023.e21147

20. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. 2013;14:7. doi:10.1186/1471-2105-14-7

21. Lee CY, Bu LX, DeBenedetti A, Williams BJ, Rennie PS, Jia WW. Transcriptional and translational dual-regulated oncolytic herpes simplex virus type 1 for targeting prostate tumors. Mol Ther. 2010;18(5):929–935. doi:10.1038/mt.2010.26

22. Zhang W, Shang X, Liu N, et al. ANK2 as a novel predictive biomarker for immune checkpoint inhibitors and its correlation with antitumor immunity in lung adenocarcinoma. BMC Pulm Med. 2022;22(1):483. doi:10.1186/s12890-022-02279-2

23. Pan JH, Zhou H, Cooper L, et al. LAYN is a prognostic biomarker and correlated with immune infiltrates in gastric and colon cancers. Front Immunol. 2019;10:6. doi:10.3389/fimmu.2019.00006

24. Xu M, Kong Y, Chen N, et al. Identification of immune-related gene signature and prediction of CeRNA network in active ulcerative colitis. Front Immunol. 2022;13:855645. doi:10.3389/fimmu.2022.855645

25. Hao Y, Hao S, Andersen-Nissen E, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184(13):3573–3587.e3529. doi:10.1016/j.cell.2021.04.048

26. Lun AT, Bach K, Marioni JC. Pooling across cells to normalize single-cell RNA sequencing data with many zero counts. Genome Biol. 2016;17:75. doi:10.1186/s13059-016-0947-7

27. Wolock SL, Lopez R, Klein AM. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Cell Syst. 2019;8(4):281–291.e289. doi:10.1016/j.cels.2018.11.005

28. Guo C, Qu X, Tang X, et al. Spatiotemporally deciphering the mysterious mechanism of persistent HPV-induced malignant transition and immune remodelling from HPV-infected normal cervix, precancer to cervical cancer: integrating single-cell RNA-sequencing and spatial transcriptome. Clin Transl Med. 2023;13(3):e1219. doi:10.1002/ctm2.1219

29. Tang S, Li S, Shi X, et al. CALCRL induces resistance to daunorubicin in acute myeloid leukemia cells through upregulation of XRCC5/TYK2/JAK1 pathway. Anticancer Drugs. 2024;35(2):163–176. doi:10.1097/CAD.0000000000001547

30. Huang Z, Zhang H, Xing C, et al. Identification and validation of CALCRL-associated prognostic genes in acute myeloid leukemia. Gene. 2022;809:146009. doi:10.1016/j.gene.2021.146009

31. Xing C, Yin H, Yao ZY, Xing XL. Prognostic signatures based on ferroptosis- and immune-related genes for cervical squamous cell carcinoma and endocervical adenocarcinoma. Front Oncol. 2021;11:774558. doi:10.3389/fonc.2021.774558

32. Jia Q, Wang H, Xiao X, et al. UCHL1 acts as a prognostic factor and promotes cancer stemness in cervical squamous cell carcinoma. Pathol Res Pract. 2023;247:154574. doi:10.1016/j.prp.2023.154574

33. Mao R, Tan X, Xiao Y, et al. Ubiquitin C-terminal hydrolase L1 promotes expression of programmed cell death-ligand 1 in non-small-cell lung cancer cells. Cancer Sci. 2020;111(9):3174–3183. doi:10.1111/cas.14529

34. Gründker C, Emons G. Role of gonadotropin-releasing hormone (GnRH) in ovarian cancer. Reprod Biol Endocrinol. 2003;1:65. doi:10.1186/1477-7827-1-65

35. Limonta P, Moretti RM, Montagnani Marelli M, Motta M. The biology of gonadotropin hormone-releasing hormone: role in the control of tumor growth and progression in humans. Front Neuroendocrinol. 2003;24(4):279–295. doi:10.1016/j.yfrne.2003.10.003

36. Hebbar V, Damera G, Sachdev GP. Differential expression of MUC genes in endometrial and cervical tissues and tumors. BMC Cancer. 2005;5:124. doi:10.1186/1471-2407-5-124

37. Zhang J, Wang L, Jiang J, Qiao Z. Elevation of microRNA-512-5p inhibits MUC1 to reduce radioresistance in cervical cancer. Cell Cycle. 2020;19(6):652–665. doi:10.1080/15384101.2019.1711314

38. Yin L, Zhou Y, Hong S, Ding F, Cai H. Strategies for synthesizing and enhancing the immune response of cancer vaccines based on MUC1 glycopeptide antigens. Chembiochem. 2023;24(10):e202200805. doi:10.1002/cbic.202200805

39. Du JJ, Zhou SH, Cheng ZR, et al. MUC1 specific immune responses enhanced by Coadministration of liposomal DDA/MPLA and lipoglycopeptide. Front Chem. 2022;10:814880. doi:10.3389/fchem.2022.814880

40. Goumans MJ, Ten Dijke P. TGF-β signaling in control of cardiovascular function. Cold Spring Harb Perspect Biol. 2018;10(2):a022210. doi:10.1101/cshperspect.a022210

41. Lu X, Liu R, Liao Y, et al. ACVRL1 drives resistance to multitarget tyrosine kinase inhibitors in colorectal cancer by promoting USP15-mediated GPX2 stabilization. BMC Med. 2023;21(1):366. doi:10.1186/s12916-023-03066-4