A novel programmed-cell-death-related prognostic risk model for cervical cancer based on mitochondrial genes

Introduction

Cervical cancer is a malignant neoplasm that arises from cervix, the lower part of the uterus that connects to the vagina (1). This type of cancer is predominantly associated with persistent infection by high-risk human papillomavirus (HPV) strains, which are responsible for the majority of cervical cancer cases (2). The development of cervical cancer is a multifactorial process influenced by various risk factors, including sexual behavior, immune response, smoking, and coinfections with other sexually transmitted infections (3). The clinical presentation of cervical cancer can be asymptomatic in its early stages, but as the disease progresses, patients may experience symptoms such as abnormal vaginal bleeding, vaginal discharge, and pelvic pain (4). The diagnosis of cervical cancer typically involves a combination of cytological screening (Pap smear), HPV DNA testing, and, if necessary, colposcopy and biopsy (5). Treatment options for cervical cancer encompass a range of surgical, radiation, and chemotherapy approaches, depending on the stage of the disease and the patient’s overall health (6). Despite advances in prevention and treatment, cervical cancer remains a significant global health challenge, particularly in regions with limited access to screening and vaccination programs (7,8).

Programmed cell death (PCD), is a regulated biological process (BP) that eliminates damaged or unnecessary cells in a controlled manner, maintaining tissue homeostasis and contributing to the development of multicellular organisms (9). PCD plays a vital role in cancer growth and metastasis (10). PCD-index has been established as a prognostic indicator for cervical cancer (11). The expression profiles of genes associated with PCD were found to have a significant correlation with the clinical manifestations observed in patients with cervical squamous cell carcinoma (12). Additionally, PCD is regulated by various organelles, among which mitochondria play a central role (13). Mitochondrial (mt) dysfunction, such as reduced oxidative phosphorylation and accumulation of reactive oxygen species, can trigger multiple PCD pathways, including apoptosis, necroptosis, and ferroptosis (13,14). As the “energy factories” of the cell, mitochondria exert complex and critical functions in tumor initiation, progression, metastasis, and therapy resistance (15,16). Therefore, integrating PCD- and mitochondria-related genes to construct prognostic models may provide a theoretical basis for the identification of novel biomarkers in the context of precision oncology. At present, studies on prognostic models based on PCD- and mitochondria-related genes in cervical cancer remain limited. Accordingly, this study aims to develop a robust risk scoring model to more effectively predict the prognosis of patients with cervical cancer. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2806/rc).

Methods

Research object

From The Cancer Genome Atlas (TCGA; https://tcga-data.nci.nih.gov/tcga/) database, mRNA expression profile data and corresponding clinical information for 309 cases of cervical cancer (TCGA-CESC) patients were downloaded, including 306 cervical cancer samples and three adjacent normal samples. Additionally, somatic mutation data were obtained.

From the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) database, the dataset with the accession number GSE173097 was downloaded, containing five cervical cancer samples and six normal samples. The E-MTAB-11948 cervical cancer single-cell dataset, containing three cervical cancer samples and two normal samples, was downloaded from the European Molecular Biology Laboratory’s European Bioinformatics Institute (EMBL-EBI) database (https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-11948).

From the MitoCarta 3.0 database (17), 1,136 human protein-coding genes were obtained, and 1,705 mt-related genes were downloaded from the MitoProteome database (http://mitoproteome.org/). The union of these two datasets resulted in a collection of 2,030 mt-related genes. Then a collection of 1,692 genes associated with PCD was compiled by integrating two datasets (18,19), and these mt and PCD-related gene sets are presented in table available at https://cdn.amegroups.cn/static/public/tcr-2025-1-2806-1.xlsx. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Differential gene expression analysis

Using the R language-based “limma” package (20), differentially expressed genes (DEGs) were screened based on the criteria of absolute value of log-transformed fold change (Log₂FC) greater than 0.6 and adjusted P value less than 0.05.

Functional enrichment analysis

Based on the DEGs between cervical cancer samples and normal samples, we utilized the “clusterProfiler” package (21) to conduct functional enrichment analysis. This analysis targeted Gene Ontology (GO) terms, encompassing BP, molecular function (MF), and cellular component (CC), as well as pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG). The identification of GO terms and KEGG pathways that showed significant enrichment was based on a P value cutoff of below 0.05.

Least absolute shrinkage and selection operator (LASSO) Cox regression analysis

Candidate genes were first subjected to univariate Cox regression analysis, and those significantly associated with cervical cancer prognosis (P<0.05) were retained. To avoid overfitting and address collinearity among candidate variables, we then applied LASSO Cox regression using the “glmnet” R package. The penalty parameter (λ) was optimized via 10-fold cross-validation, and the final eight-gene signature was selected based on the λ.min criterion (minimum cross-validation error). To ensure reproducibility, a fixed random seed (set.seed(123)) was used, which consistently yielded the same gene combination across multiple runs. This signature also demonstrated the optimal AIC among candidate models and maintained stable predictive performance in subsequent bootstrap validation. A RiskScore for each sample was then calculated using the genes selected, based on the formula provided:

Riskscore=∑i=1nCoefi∗Xi[1]

In the analysis, Coefi represents the risk coefficient calculated for each factor in the model, while Xi denotes the expression value of each factor, specifically referring to the mRNA expression levels in this study. Subsequently, the R packages “survival”, “survminer”, and two-sided log-rank tests were employed to determine the values of the RiskScore. Patients were then categorized into low risk group (LRG) and high risk group (HRG) based on the median value of the RiskScore.

Gene set enrichment analysis (GSEA)

The GSEA analysis was conducted between the HRG and LRG groups using the “clusterProfiler” (21) package in R, and pathways with p.adjust <0.05 were filtered for significant enrichment.

Survival analysis

The overall survival (OS) rates of HRG and LRG were estimated based on the Kaplan-Meier method using the “survival” package and “survminer” package in R, and the significance of the differences in survival rates between groups was tested using the log-rank test. The ability of the RiskScore to predict the survival of cervical cancer patients independently of other factors was analyzed using a multivariate Cox regression model. The receiver operating characteristic (ROC) curve was generated using the “timeROC” package (22) in R, and the area under the curve (AUC) was calculated.

Nomogram prognostic prediction model

Nomograms are widely used tools for predicting cancer prognosis. In our research, we leveraged “rms” package (https://CRAN.R-project.org/package=rms) to create a Nomogram incorporating significant prognostic variables identified by a multivariate Cox proportional hazards model. To ensure its predictive accuracy, we produced a calibration plot illustrating the correlation between the nomogram’s estimated probabilities and the actual survival rates. Each patient’s profile was scored by projecting lines vertically from their factor values on the nomogram to the “total points” scale, further predicting 1-, 3-, and 5-year survival rates for cervical cancer patients.

Mutation analysis

The GSCA database at https://guolab.wchscu.cn/GSCA/ was utilized to analyze the copy number variation (CNV) differences of the genes in the RiskScore model. Then we analyzed tumor mutation burden (TMB) using “maftools” package based on MAF files of somatic mutations from TCGA database.

Immune infiltration analysis

To determine the presence and proportions of various immune cell types within the samples based on their gene expression profiles, we utilized the Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) (23). This software employs a mathematical deconvolution method and references a curated list of 547 specific gene markers, known as barcode genes. Leveraging this approach, we were able to estimate the relative infiltration ratios of 22 distinct immune cell populations within the tissue samples.

The “ESTIMATE” package (24) was utilized in our study to predict the levels of stromal and immune cell infiltration within the tumor tissues based on gene expression data, yielding the stromal score and the immune score.

Single cell sequencing analysis

The “Seurat” package (25) was employed for the quality control of single-cell data within the dataset. Cells with low quality were filtered out based on criteria of having fewer than 200 features and a mt percentage exceeding 10%. Subsequently, the “NormalizeData” function with the “LogNormalize” method was applied to standardize the data. The “FindVariableFeatures” function was then utilized to identify genes with high cell-to-cell variability. The top 2,000 most variable genes were selected for further scaling of the data using the “ScaleData” function, achieving normalization. Principal component analysis (PCA) was conducted for linear dimensionality reduction, followed by the application of the uniform manifold approximation and projection (UMAP) method for nonlinear dimensionality reduction. DEGs within each cluster were identified using the “FindAllMarkers” function. Cellular clusters were annotated using “SingleR” (26) and the CellMarker 2.0 database.

Statistical analysis

The Wilcoxon rank-sum test was utilized to compare gene expression differences and immune cell infiltration differences between different groups. Pearson correlation analysis was conducted using the R language to assess the relationships between variables. A P value of less than 0.05 was considered to indicate statistical significance. All statistical analyses were performed using R software (version 4.3.3).

Results

Identification of PCD-related mt genes in cervical cancer

Differential gene expression analysis was performed on tumor versus normal samples utilizing the GSE173097 dataset. A total of 5,614 DEGs were detected in tumor samples compared to normal samples, including 3,142 upregulated genes and 2,472 downregulated genes (Figure 1A,1B, table available at https://cdn.amegroups.cn/static/public/tcr-2025-1-2806-1.xlsx). Subsequently, functional enrichment analysis was conducted on these DEGs, revealing a total of 101 BP terms, including chromosome segregation and leukocyte chemotaxis, 9 CC terms such as condensed chromosome, 17 MF terms including cytokine activity and chemokine activity, and 49 KEGG pathways such as cytokine-cytokine receptor interaction and interleukin-17 (IL-17) signaling pathway that were significantly enriched. The top 5 significantly enriched GO terms of each category and the top 15 KEGG pathways with notable enrichment are displayed in (Figure 1C,1D), with detailed enrichment results available in table available at https://cdn.amegroups.cn/static/public/tcr-2025-1-2806-2.xlsx. An intersection of the DEGs, mt genes, and PCD-related gene sets yielded 119 candidate genes (Figure 1E, table available at https://cdn.amegroups.cn/static/public/tcr-2025-1-2806-1.xlsx).

Figure 1 Identification of candidate PCD-related mt genes in cervical cancer. (A) Differential gene expression between cancer samples and normal samples based on GSE173097 dataset. (B) The heatmap of top 50 DEGs between cancer samples and normal samples based on GSE173097 dataset. (C) Top 5 significantly enriched GO terms (BP, CC, MF). (D) Top 15 significantly enriched KEGG pathways. (E) Intersection of PCD-related genes, mt genes and DEGs based on GSE173097 dataset. BP, biological process; CC, cellular component; DEGs, differentially expressed genes; FC, fold change; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; mt, mitochondrial; PCD, programmed cell death.

Construction of a prognostic model for cervical cancer

Based on the candidate genes, a prognostic model for cervical cancer was established. The sample data from the TCGA cohort were randomly divided into a training set and a validation set in a ratio of 7:3. The expression values of the 119 candidate genes in the training set samples were used as continuous variables for univariate Cox regression analysis, and the hazard ratio (HR) for each gene was calculated. Genes were screened with a threshold of P value <0.05, resulting in nine genes that met the criteria, including IL1B, NAMPT, PRKAB2, SLC2A1, ACAA2, E2F1, G0S2, GZMB, and BNIP3 (Figure 2A). Subsequently, LASSO Cox regression analysis was performed on these nine genes using the training set samples. Based on the lambda values corresponding to different numbers of genes in the LASSO Cox analysis, the optimal number of genes was determined to be eight, corresponding to the smallest lambda value (Figure 2B,2C). A prognostic model was constructed using the selected 8 genes along with the regression coefficients from the LASSO Cox regression analysis, from which a risk score was calculated. The risk score formula was:

Riskscore=0.00269770218544548×IL1B+0.00743162329528295×NAMPT+0.00929530180158692×PRKAB2+0.000353525288866902×SLC2A1+0.0149427336222402×ACAA2−0.00640561519647769×E2F1−0.00811333089202734×GZMB+0.00171520123991157×BNIP3[2]

Figure 2 Establishment of a prognostic model for cervical cancer. (A) Univariate Cox regression analysis of IL1B, NAMPT, PRKAB2, SLC2A1, ACAA2, E2F1, G0S2, GZMB, and BNIP3. (B) The relationship between LASSO penalty and regression coefficient change. (C) Determination of λ in LASSO Cox model by cross validation. (D) KM survival curve of LRG and HRG in training set from TCGA dataset. (E) RiskScore distribution (upper panel), survival status (middle panel), and expression levels in the RiskScore model (bottom panel) based on the training set. (F) KM survival curve of LRG and HRG in validation set from TCGA dataset. (G) RiskScore distribution, survival status, and expression levels in the model based on the validation set. (H) DCA. CI, confidence interval; DCA, decision curve analysis; KM, Kaplan-Meier; HR, hazard ratio; HRG, high risk group; LASSO, least absolute shrinkage and selection operator; LRG, low risk group; TCGA, The Cancer Genome Atlas.

For the training set, the RiskScore for each patient was calculated, and the samples were divided into LRG and HRG based on the median RiskScore. Survival analysis revealed that HRG samples had poorer OS compared to LRG samples (Figure 2D,2E). Moreover, we performed 5-fold cross-validation to reassess model performance. Results showed that in folds 2, 3, and 5, patients classified into the HRG exhibited significantly poorer prognosis compared to those in the LRG, and the mean AUC for the RiskScore was 0.773 across all five folds (Figure S1). The prognostic model developed from the training set was validated using the validation set. The results showed that the OS of patients in the HRG was significantly lower than that of the LRG, with a P value of 0.049 (Figure 2F,2G). Additionally, decision curve analysis demonstrated that this model provides positive clinical net benefits in survival prediction at 1 year (threshold 0.2–0.3), 3 years (threshold 0.5–0.8), and 5 years (threshold 0.7–0.9), and outperforms both the ‘treat all’ and ‘treat none’ reference strategies at all time points (Figure 2H).

Validation and independence of RiskScore model

To evaluate the predictive efficacy of the prognostic model, we assessed its performance across various time points and its integration with established clinical parameters. An ROC curve was plotted for the training dataset, yielding an AUC of 0.85 at 1 year, 0.793 at 3 years, and 0.726 at 5 years,
suggesting that the predictive model could effectively forecast patient prognosis (Figure 3A). The forecasting capability was further validated using the validation set, where the ROC curves indicated AUCs of 0.803 at 1 year, 0.633 at 3 years, and 0.68 at 5 years (Figure 3B).

Figure 3 Validation and independence of RiskScore model. (A,B) ROC curve based on the training set and validation set. (C) Multivariate Cox regression analysis of RiskScore, age, and TNM stages. (D-G) Comparison of RiskScores in different groups of age, T, M, and N stages, respectively. *, P<0.05; **, P<0.01; ***, P<0.001. AIC, Akaike Information Criterion; AUC, area under the curve; CI, confidence interval; HR, hazard ratio; ROC, receiver operating characteristic; TNM, tumor-node-metastasis.

Subsequently, a multivariate Cox regression analysis was performed on the prognostic model, incorporating RiskScores along with age, tumor-node-metastasis (TNM) staging as factors. The results indicated that the RiskScore could serve as an independent factor to predict prognosis of cervical cancer (P<0.001), with the T4 subgroup exhibiting statistical significance (P<0.05), while the P values for age, other subgroups of T staging, and M staging exceeded 0.05 (Figure 3C). The relationship between RiskScores and various age groups, as well as TNM staging, was analyzed, revealing no significant differences in different age groups, M stages, and N stages. However, within the T staging, significant differences were observed in RiskScores between T1 and T3, T2 and T3, and T1 and T4 subgroups, with higher stages corresponding to elevated RiskScores (Figure 3D-3G).

Nomogram model for cervical cancer patient survival prediction

After identifying the RiskScore and T stage as independent prognostic factors for cervical cancer patients, a Nomogram model was constructed to predict the 1-, 3-, and 5-year survival rates based on these two factors (Figure 4A). Calibration curves for the 1-, 3-, and 5-year predictions were plotted. The calibration curves for the respective years closely resembled the ideal curve (a 45-degree line passing through the origin with a slope of 1), indicating that the nomogram model’s predictions for 1-, 3-, and 5-year survival rates were in good agreement with the actual outcomes (Figure 4B-4D).

Figure 4 Nomogram model for cervical cancer patients’ survival prediction. (A) Nomogram for predicting 1-, 3-, and 5-year OS rate in cervical cancer patients. (B-D) Calibration curves of the nomogram model predicting 1-, 3-, and 5-year OS rate. *, P<0.05; ***, P<0.001.OS, overall survival; T, tumor.

Potential pathways influenced by RiskScore

Having established the prognostic significance of the RiskScore, we next sought to further elucidate the functional differences between the risk groups and explore the potential BPs driving the prognostic disparity. Thus, GSEA analysis was performed between LRG and HRG in TCGA training set due to its relatively large sample size. A total of 82 pathways were found to be enriched (P value <0.05). In HRG, the significantly activated pathways included extracellular matrix (ECM)-receptor interaction, focal adhesion, microRNAs in cancer, cytoskeleton in muscle cells, the IL-17 signaling pathway, etc. (Figure 5A). In contrast, the significantly downregulated pathways comprised autoimmune thyroid disease, allograft rejection, antigen processing and presentation, T helper 1 (Th1) and Th2 cell differentiation, the intestinal immune network for immunoglobulin A (IgA) production, etc. (Figure 5B). Detailed enrichment results are presented in table available at https://cdn.amegroups.cn/static/public/tcr-2025-1-2806-3.xlsx.

Figure 5 Potential pathways influenced by RiskScore. (A) Top 5 significantly enriched activated pathways analyzed by GSEA. (B) Top 5 significantly downregulated pathways. ECM, extracellular matrix; GSEA, gene set enrichment analysis; IL, interleukin.

Tumor mutation analysis based on RiskScore

Cancer gene mutations are predominantly observed across all cancer types, with patients presenting with tumors that harbor distinct combinations of gene mutation profiles exhibiting varied survival outcomes (27). This suggests the significant value of mutation analysis in oncology.

We analyzed the CNVs of the eight genes in the RiskScore model using the GSCA database and found that BNIP3 and ACAA2 had a relatively high proportion of heterozygous deletions, while the proportion of heterozygous amplifications in PRKAB2 and E2F1 was relatively high (Figure 6A, table available at https://cdn.amegroups.cn/static/public/tcr-2025-1-2806-4.xlsx). Concurrently, the mutation frequency and types of the eight genes were detected using somatic mutation data, and a low mutation rate for these genes was observed (Figure 6B). Furthermore, the somatic mutation data between LRG and HRG was analyzed, focusing on the top 30 most mutated genes. The analysis revealed that in HRG, the most frequently mutated genes were PIK3CA, TTN, MUC16, KMT2C, SYNE1, FBXW7, and FLG (Figure 6C); whereas in LRG, they were TTN, PIK3CA, DMD, KMT2C, MUC16, CREBBP, and EP300 (Figure 6D). Additionally, the correlation between the RiskScore and TMB was accessed, which displayed a significant negative correlation (R=−0.2, P=0.005) (Figure 6E).

Figure 6 Tumor mutation analysis based on RiskScore. (A) CNV analysis of the eight genes in RiskScore model. (B) Mutation analysis of the eight genes in RiskScore model. (C,D) Mutation analysis of top 30 genes in HRG and LRG. (E) Correlation of TMB and RiskScore. CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; CNV, copy number variation; HRG, high risk group; LRG, low risk group; NA, not available; TMB, tumor mutation burden.

Tumor immune microenvironment (TME) analysis based on RiskScore

Cancer is increasingly recognized as a disease influenced by the TME, where the processes of tumor formation, advancement, and metastasis are intrinsically linked to the TME (28). Using the TCGA training set dataset, the CIBERSORT algorithm was applied to estimate the relative infiltrating proportions of 22 distinct immune cell types (Figure 7A). An examination of the variation in immune cell infiltration between HRG and LRG revealed significant differences in infiltration of 13 immune cell types. Specifically, CD4 memory resting T cells, resting natural killer (NK) cells, M0 macrophages, activated dendritic cells, activated mast cells, and eosinophils exhibited higher infiltration ratios in HRG than LRG. Conversely, resting mast cells, activated NK cells, plasma cells, CD4 memory activated T cells, CD8 T cells, follicular helper T cells, regulatory T cells (Tregs) infiltrating significantly lower in HRG compared to LRG (Figure 7B). Further analysis of the correlation between the RiskScore and these 13 types of immune cells using Pearson correlation demonstrated significant negative associations between the RiskScore and the levels of CD4 memory activated T cells, CD8 T cells, and follicular helper T cells (P<0.05). In contrast, significant positive associations were observed with M0 macrophages, activated mast cells, resting NK cells, and CD4 memory resting T cells (Figure S2).

Figure 7 TME analysis based on RiskScore. (A) Relative infiltrating ratios of 22 types of immune cells. (B) Comparison of infiltrating ratios of 22 immune cell types between HRG and LRG. (C) Comparison of stromal score between HRG and LRG. (D) Comparison of immune score between HRG and LRG. (E) Expression levels of 8 immune checkpoints in HRG and LRG. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. HRG, high risk group; LRG, low risk group; NK, natural killer; TME, tumor immune microenvironment.

Next, the ESTIMATE algorithm was utilized to predict the levels of stroma and immune cell infiltration within the tumor tissues, resulting in the calculation of stromal scores and immune scores. It was observed that while the stromal scores did not show significant differences between HRG and LRG, the immune scores were significantly higher in LRG compared to HRG (Figure 7C,7D).

We further explored the expression levels of eight immune checkpoint genes, including CD274, PD-1, CTLA4, CD80, CD86, LAG3, PDCD1LG2, and TIGIT, between HRG and LRG. The results indicated that expression levels of CD86, CTLA4, LAG3, PD-1, and TIGIT were significantly higher in LRG than HRG (Figure 7E).

Single cell analysis

Single cell analysis was employed to dissect the cellular composition underlying the RiskScore due to its potential to reveal the intricate interactions and behaviors of various cell types within the TME, which may contribute to the observed clinical outcomes.

After stringent quality control of the single cell dataset, 39,807 cells were subjected to analysis. Following standardization, normalization, and dimensionality reduction techniques, 20 distinct cell clusters were identified. Utilizing SingleR and the CellMarker 2.0 database for annotation, this process culminated in the classification of eight cell clusters, as illustrated by the UMAP results (Figure 8A). Subsequent analysis of the eight genes within the RiskScore model revealed that GZMB was specifically highly expressed in T/NK cells, while NAMPT and IL1B were specifically highly expressed in neutrophils, and BNIP3 was specifically highly expressed in epithelial cells (Figure 8B,8C).

Figure 8 Single cell analysis. (A) UMAP of cell subtype identification. (B) Expression distribution of the eight RiskScore model in all cells. (C) The expression levels of the eight RiskScore model genes in each cell subtype. NK, natural killer; UMAP, uniform manifold approximation and projection.

Given the finding that GZMB was specifically overexpressed in tumor tissues and showed a high level of expression in T/NK cells, a sub-cluster analysis was conducted on T and NK cells, yielding six sub-clusters (Figure 9A). Pseudotime analysis of these cells indicated that sub-cluster 4 predominantly emerged in the early stages, and sub-cluster 3 followed a distinct trajectory compared to the other sub-clusters (Figure 9B,9C). Furthermore, the expression dynamics of key genes in the risk model along the pseudotime trajectory demonstrated that GZMB expression initially decreased, then increased, and subsequently decreased again (Figure 9D). Subsequently, focusing the analysis on NK cells, due to the higher expression level of GZMB in these cells, resulted in the identification of three sub-clusters (Figure 9E). The pseudotime analysis revealed that sub-cluster 0 was primarily observed in the early stages (Figure 9F,9G). Examination of the expression changes of key genes from the RiskScore model along the pseudotime trajectory showed that GZMB expression first increased, then decreased, and ultimately increased again (Figure 9H).

Figure 9 Single cell analysis focusing on GZMB gene and T/NK cells. (A) UMAP of cell subtype identification in T/NK cells. (B) Pseudotime analysis based on T/NK cells. (C) Pseudotime analysis of each subtype in T/NK cells. (D) Expression of GZMB in pseudotime based on T/NK cells. (E) UMAP of cell subtype identification in NK cells. (F) Pseudotime analysis based on NK cells. (G) Pseudotime analysis of each subtype in NK cells. (H) Expression of GZMB in pseudotime NK cells. NK, natural killer; UMAP, uniform manifold approximation and projection.

Discussion

In recent years, multiple prognostic models have been developed for cervical cancer based on distinct functional gene sets, including immune checkpoint HLA-G-driven DEGs (29), apoptosis-related genes (30), and lysosome-related genes (31). Hu et al. constructed a pyroptosis-related three-gene signature (CHMP4C, GZMB, TNF) associated with the tumor immune microenvironment, achieving AUC values of 0.733, 0.710, and 0.673 for 1-, 2-, and 3-year survival, respectively (32). Similarly, Yang et al. established a senescence-related six-gene prognostic signature (RPS6KA6, ABI3, PTTG1, E2F1, CBX7, SPOP) that was significantly correlated with OS in cervical cancer. Their model yielded AUCs of 0.764, 0.743, and 0.780 for 1-, 3-, and 5-year survival, and when integrated with clinicopathological factors, the nomogram achieved an AUC of 0.860 for 5-year OS (33). In this study, we constructed a RiskScore model based on eight mt-related genes involved in PCD, namely, IL1B, NAMPT, PRKAB2, SLC2A1, ACAA2, E2F1, GZMB, and BNIP3, to predict prognosis in cervical cancer patients. In the training dataset, our model yielded AUC values of 0.85, 0.793, and 0.726 for 1-, 3-, and 5-year survival, respectively. In contrast to existing prognostic models that primarily focus on specific pathways such as apoptosis, pyroptosis, or senescence, we constructed a model centered on genes associated with both mitochondria and PCD. This biology-driven model, grounded in mt metabolism, offers a novel perspective for deciphering tumor metabolic reprogramming and its intricate crosstalk with the TME.

In cervical cancer tissues, IL-1B expression levels are significantly higher than in normal tissues, and its high expression is closely correlated with poor patient prognosis (34). NAMPT is upregulated in various cancers such as glioma and lung cancer (35), and has been established as a marker of cervical intraepithelial neoplasia and squamous cell carcinoma (36). As a key metabolic enzyme, treatment with the NAMPT-specific inhibitor FK866 in gastric cancer cells significantly suppresses the expression of pro-oncogenic molecules such as NF-κB (37). Moreover, inhibition of the NF-κB/IL-1B signaling axis has been demonstrated to effectively prevent tumor metastasis (38). These findings suggest that the pro-tumor effects of NAMPT and IL-1B may be associated with NF-κB signaling. Importantly, evidence has shown that chronic inflammation plays a critical role in cancer development (39). NF-κB not only serves as a regulator governing inflammation and immune responses, but also exhibits profound bidirectional crosstalk with cellular metabolic reprogramming (40). As the principal transporter mediating basal glucose uptake, SLC2A1 (also known as GLUT1) drives glycolysis to meet the energy demands of cell growth. SLC2A1 is also highly expressed in cervical cancer tissues (41). Meanwhile, studies have shown that NF-κB can upregulate the expression of SLC2A1 (42), suggesting a potential relationship between NF-κB and glycolysis. However, when glucose levels is depleted in the environment, cells shift their energy production to fatty acid oxidation. Glucose and fatty acids are two crucial metabolites, which makes glycolysis and fatty acid oxidation essential metabolic pathways for generating energy. These pathways influence one another, and in the absence of one substrate, cells can still ensure adequate energy production due to their metabolic adaptability (43). During conditions of glucose deprivation, the key glycolytic enzyme PFKP triggers fatty acid oxidation via a non-glycolytic mechanism, specifically, by binding directly to AMPK and bringing it to the mitochondria. This mechanism contributes to sustaining energy balance and redox homeostasis, thereby enhancing cancer cell survival during glucose scarcity (44). Interestingly, PRKAB2 functions as a regulatory subunit of the AMPK complex, contributing to energy sensing and metabolic homeostasis (45). Meanwhile, ACAA2, a key enzyme in the mt fatty acid β-oxidation pathway, catalyzes the final thiolysis reaction to generate acetyl-CoA (46). In triple-negative breast cancer, enhanced ACAA2 activity directly promotes fatty acid metabolism and free fatty acid accumulation, providing the material and energetic basis for tumor growth (47). Furthermore, under hypoxic conditions, BNIP3 drives metabolic reprogramming in uveal melanoma by inducing mitophagy, thereby promoting tumor cell survival and metastasis (48). Together, these genes (e.g. SLC2A1, PRKAB2, ACAA2, BNIP3) enable cancer cells to flexibly utilize glucose and fatty acids even under hypoxic or nutrient-deprived conditions, supporting the energy and biosynthetic demands of rapid proliferation.

Furthermore, E2F1 is a key transcription factor governing cell cycle progression and DNA replication, and its activity is regulated by metabolic signals. Evidence indicates that mt stress (e.g., accumulation of reactive oxygen species) can activate the E2F1 signaling pathway, triggering a cascade that leads to apoptosis (49). GZMB is a core effector molecule through which cytotoxic T lymphocytes and NK cells execute their killing functions (50,51), and its high expression in T/NK cells may reflect an immune-activated state (52). In this study, single-cell transcriptome analysis has demonstrated that GZMB is almost exclusively and specifically expressed in T/NK cells. These findings indicate that the role of GZMB in tumors does not originate from the tumor cells, but rather serves as an effective indicator of the cytotoxic functional state of T/NK cells.

Thus, these eight genes collectively constitute a coherent biological narrative: from IL1B/NAMPT-driven inflammatory priming and basal metabolic remodeling, through SLC2A1/PRKAB2/ACAA2/BNIP3-mediated metabolic flexibility, to the modulation of E2F1 signaling by metabolic cues and the immune effector dysfunction reflected by GZMB, these multi-layered, multi-pathway synergistic interactions together shape a tumor microenvironment characterized by hyperactive metabolism, and enhanced cell growth. This synergistic network, rather than the isolated effects of individual genes, provides a more comprehensive explanation for the poor prognosis observed in HRG patients. Future research should focus on elucidating the cross-talk mechanisms among these pathways and exploring their translational potential as targets for combination therapies.

In this study, we analyzed the association between the RiskScore and the tumor immune microenvironment (as shown in Figure 7). The results revealed that the HRG exhibited significantly higher infiltration levels of M0 macrophages, a phenotype typically indicative of a tumor microenvironment prone to supporting tumor progression (53). In contrast, the LRG showed elevated infiltration of activated CD4⁺ memory T cells and CD8⁺ T cells, suggesting a more robust immune surveillance system that may contribute to more effective clearance of abnormal cells. Furthermore, compared with the HRG, the LRG displayed higher expression levels of several immune checkpoint molecules. Evidence indicates that high expression of PLEK2 or IFI6 in esophageal squamous cell carcinoma is associated with an immunosuppressive tumor microenvironment, which may limit the clinical benefit from immunotherapy alone (54). Similarly, high expression of PRDX4/5/6 has been associated with unfavorable outcomes, whereas low expression of PRDX1/5 correlates with enhanced immune cell infiltration, upregulation of immune-related molecules, and a greater likelihood of response to anti-PD-1 therapy (55). Consistently with these findings, in the present research, patients in the LRG, characterized by enhanced infiltration of CD4⁺ T and CD8⁺ T cells and higher expression of immune checkpoint molecules, may represent a more suitable population for immune checkpoint inhibitor therapy.

Furthermore, we observed that the mutation frequencies of the eight genes (IL1B, NAMPT, PRKAB2, SLC2A1, ACAA2, E2F1, GZMB, and BNIP3) were generally low, which in fact reinforces their functional essentiality. Loss of genomic stability is associated with elevated mutation rates (56). A low mutation rate often indicates high functional conservation of a gene (57). In general, when a gene is extensively conserved among various species, this similarity across species indicates that the gene fulfills essential roles critical for a wide range of life forms, thus being preserved throughout the course of evolution (57). Our results demonstrate that all eight genes exhibit consistently low mutation frequencies, indicating that their prognostic significance arises not from genetic alterations, but from dynamic changes in expression levels. This low mutation burden may constitute a key advantage of our model: the high-risk state does not represent a fixed, irreversible “genotype”, but rather a reversible phenotypic state shaped at the transcriptional level. This implies that the risk status of patients may be modifiable through targeting gene expression regulation, positioning our model as a potential tool for monitoring of therapeutic response.

There are several important limitations in this study. First, its retrospective design based on public datasets means our findings require validation in independent, prospective clinical cohorts. Second, while we linked the model genes to cancer processes through bioinformatic analysis, experimental studies are needed to confirm their biological roles. Meanwhile, in future work, we also plan to validate the expression of key model genes using cervical cancer clinical specimens or cell lines via quantitative real-time polymerase chain reaction (qRT-PCR) or Western blot. Third, our TME and immunotherapy analyses provide only a broad view; deeper investigation using single-cell or spatial methods would better reveal the tumor-immune interactions. Forth, while TMB is one of several established predictors of immunotherapy response (58-60), our study identified a statistically significant but weak negative correlation between the risk score and TMB (R=−0.2, P=0.005). This result suggests only a potential trend and does not provide sufficient evidence to conclude that patients in the LRG would definitively derive greater benefit from immunotherapy. In future work, we plan to directly evaluate the prognostic and diagnostic value of the RiskScore in the context of immunotherapy, rather than inferring it indirectly through TMB, using cervical cancer immunotherapy cohorts with long-term follow-up data. Fifth, HPV is a key driver of cervical carcinogenesis, and incorporating HPV-related genes or expression profiles is indeed important for refining prognostic models. Unfortunately, the TCGA-CESC dataset used in this study does not provide HPV genotyping or viral gene expression data, which precluded us from including this information in the current analysis. This limitation reflects a common challenge in public database-driven research. To address this, we plan to integrate external cervical cancer cohorts that contain HPV status for further validation, and collect clinical samples with HPV testing in our future experimental work to optimize the current model in the future. Similarly, regarding treatment and geographic information, these data were not available in the TCGA-CESC public dataset and thus were not included in the current analysis. Recognizing that these factors may impact prognostic performance, future validation in independent cohorts with comprehensive clinical annotations will be essential to assess the model’s applicability across diverse treatment settings and populations.

Conclusions

To summarize, this study developed a RiskScore model based on PCD-related mt genes, IL1B, NAMPT, PRKAB2, SLC2A1, ACAA2, E2F1, GZMB, and BNIP3. The prognostic model demonstrated accuracy and independency in predicting cervical cancer prognosis. The nomogram model, incorporating the RiskScore and T stage, provided a reliable tool for predicting 1-, 3-, and 5-year survival rates. Furthermore, analysis of TMB revealed potential links between the RiskScore and genomic instability, while investigation of the TME uncovered distinct immune infiltration patterns associated with patient risk. The single-cell analysis deepened our understanding of cellular heterogeneity within TME. Our findings contribute a step to improving cervical cancer prognosis prediction, facilitating personalized treatment strategies, and identifying potential targets for future research and clinical applications.

Acknowledgments

None.

Footnote

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

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

Funding: None.

Conflicts of Interest: Both authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-1-2806/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. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

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. Cervical Cancer Treatment2024. Available online: https://www.ncbi.nlm.nih.gov/books/NBK65985/
2. Okunade KS. Human papillomavirus and cervical cancer. J Obstet Gynaecol 2020;40:602-8. [Crossref] [PubMed]
3. Luvián-Morales J, Gutiérrez-Enríquez SO, Granados-García V, et al. Risk factors for the development of cervical cancer: analysis of the evidence. Front Oncol 2024;14:1378549. [Crossref] [PubMed]
4. Burmeister CA, Khan SF, Schäfer G, et al. Cervical cancer therapies: Current challenges and future perspectives. Tumour Virus Res 2022;13:200238. [Crossref] [PubMed]
5. Koliopoulos G, Nyaga VN, Santesso N, et al. Cytology versus HPV testing for cervical cancer screening in the general population. Cochrane Database Syst Rev 2017;8:CD008587. [Crossref] [PubMed]
6. Palumbo M, Della Corte L, Ronsini C, et al. Surgical Treatment for Early Cervical Cancer in the HPV Era: State of the Art. Healthcare (Basel) 2023;11:2942. [Crossref] [PubMed]
7. Omotoso O, Teibo JO, Atiba FA, et al. Addressing cancer care inequities in sub-Saharan Africa: current challenges and proposed solutions. Int J Equity Health 2023;22:189. [Crossref] [PubMed]
8. Kaur S, Sharma LM, Mishra V, et al. Challenges in Cervical Cancer Prevention: Real-World Scenario in India. South Asian J Cancer 2023;12:9-16. [Crossref] [PubMed]
9. Programmed Cell Death (Apoptosis). Available online: https://www.ncbi.nlm.nih.gov/books/NBK26873/
10. Ma Z, Wang Y, Yu Y, et al. Exploring diverse programmed cell-death patterns to develop a novel gene signature for predicting the prognosis of lung adenocarcinoma patients. J Thorac Dis 2024;16:911-23. [Crossref] [PubMed]
11. Wang W, Chen P, Yuan S. Programmed cell death-index (PCDi) as a prognostic biomarker and predictor of drug sensitivity in cervical cancer: a machine learning-based analysis of mRNA signatures. J Cancer 2024;15:1378-96. [Crossref] [PubMed]
12. Feng S, Wang Z, Zhang H, et al. Identification of prognostic biomarkers for cervical cancer based on programmed cell death-related genes and assessment of their immune profile and response to drug therapy. J Gene Med 2024;26:e3643. [Crossref] [PubMed]
13. Nguyen TT, Wei S, Nguyen TH, et al. Mitochondria-associated programmed cell death as a therapeutic target for age-related disease. Exp Mol Med 2023;55:1595-619. [Crossref] [PubMed]
14. Gao Q, Han X, Wang J, et al. Crosstalk between mitochondrial quality control and novel programmed cell death in pulmonary diseases. Biomed Pharmacother 2025;189:118335. [Crossref] [PubMed]
15. Sharma A, Virmani T, Kumar G, et al. Mitochondrial signaling pathways and their role in cancer drug resistance. Cell Signal 2024;122:111329. [Crossref] [PubMed]
16. van der Merwe M, van Niekerk G, Fourie C, et al. The impact of mitochondria on cancer treatment resistance. Cell Oncol (Dordr) 2021;44:983-95. [Crossref] [PubMed]
17. Rath S, Sharma R, Gupta R, et al. MitoCarta3.0: an updated mitochondrial proteome now with sub-organelle localization and pathway annotations. Nucleic Acids Res 2021;49:D1541-7. [Crossref] [PubMed]
18. Cao K, Zhu J, Lu M, et al. Analysis of multiple programmed cell death-related prognostic genes and functional validations of necroptosis-associated genes in oesophageal squamous cell carcinoma. EBioMedicine 2024;99:104920. [Crossref] [PubMed]
19. Qin H, Abulaiti A, Maimaiti A, et al. Integrated machine learning survival framework develops a prognostic model based on inter-crosstalk definition of mitochondrial function and cell death patterns in a large multicenter cohort for lower-grade glioma. J Transl Med 2023;21:588. [Crossref] [PubMed]
20. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
21. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
22. Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med 2013;32:5381-97. [Crossref] [PubMed]
23. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
24. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
25. Hao Y, Stuart T, Kowalski MH, et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol 2024;42:293-304. [Crossref] [PubMed]
26. Aran D, Looney AP, Liu L, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol 2019;20:163-72. [Crossref] [PubMed]
27. Sinkala M. Mutational landscape of cancer-driver genes across human cancers. Sci Rep 2023;13:12742. [Crossref] [PubMed]
28. Wang Q, Shao X, Zhang Y, et al. Role of tumor microenvironment in cancer progression and therapeutic strategy. Cancer Med 2023;12:11149-65. [Crossref] [PubMed]
29. Xu HH, Wang HL, Xing TJ, et al. A Novel Prognostic Risk Model for Cervical Cancer Based on Immune Checkpoint HLA-G-Driven Differentially Expressed Genes. Front Immunol 2022;13:851622. [Crossref] [PubMed]
30. Zhang L, Zheng S, Chen P. Prognostic model for cervical cancer based on apoptosis-related genes. Comput Methods Biomech Biomed Engin 2025; Epub ahead of print. [Crossref]
31. Ji H, Zheng J, Liu L, et al. LAMP3 signature affects cervical cancer progression through autophagy. BMC Cancer 2025;25:1206. [Crossref] [PubMed]
32. Hu H, Yang M, Dong W, et al. A Pyroptosis-Related Gene Panel for Predicting the Prognosis and Immune Microenvironment of Cervical Cancer. Front Oncol 2022;12:873725. [Crossref] [PubMed]
33. Yang W, An L, Li Y, et al. A cellular senescence-related genes model allows for prognosis and treatment stratification of cervical cancer: a bioinformatics analysis and external verification. Aging (Albany NY) 2023;15:9408-25. [Crossref] [PubMed]
34. Wang L, Zhao W, Hong J, et al. Association between IL1B gene and cervical cancer susceptibility in Chinese Uygur Population: A Case-Control study. Mol Genet Genomic Med 2019;7:e779. [Crossref] [PubMed]
35. Lin TC. Updated Functional Roles of NAMPT in Carcinogenesis and Therapeutic Niches. Cancers (Basel) 2022;14:2059. [Crossref] [PubMed]
36. Vora M, Alattia LA, Ansari J, et al. Nicotinamide Phosphoribosyl Transferase a Reliable Marker of Progression in Cervical Dysplasia. Anticancer Res 2017;37:4821-5. [Crossref] [PubMed]
37. Bi TQ, Che XM, Liao XH, et al. Overexpression of Nampt in gastric cancer and chemopotentiating effects of the Nampt inhibitor FK866 in combination with fluorouracil. Oncol Rep 2011;26:1251-7. [Crossref] [PubMed]
38. Liu B, Zheng X, Li J, et al. Atovaquone inhibits colorectal cancer metastasis by regulating PDGFRβ/NF-κB signaling pathway. BMC Cancer 2023;23:1070. [Crossref] [PubMed]
39. Wen Y, Zhu Y, Zhang C, et al. Chronic inflammation, cancer development and immunotherapy. Front Pharmacol 2022;13:1040163. [Crossref] [PubMed]
40. Capece D, Verzella D, Flati I, et al. NF-κB: blending metabolism, immunity, and inflammation. Trends Immunol 2022;43:757-75. [Crossref] [PubMed]
41. Wang Y, Mao Y, Wang C, et al. RNA methylation-related genes of m6A, m5C, and m1A predict prognosis and immunotherapy response in cervical cancer. Ann Med 2023;55:2190618. [Crossref] [PubMed]
42. Sommermann TG, O’Neill K, Plas DR, et al. IKKβ and NF-κB transcription govern lymphoma cell survival through AKT-induced plasma membrane trafficking of GLUT1. Cancer Res 2011;71:7291-300. [Crossref] [PubMed]
43. Khan AUH, Salehi H, Alexia C, et al. Glucose Starvation or Pyruvate Dehydrogenase Activation Induce a Broad, ERK5-Mediated, Metabolic Remodeling Leading to Fatty Acid Oxidation. Cells 2022;11:1392. [Crossref] [PubMed]
44. Chen J, Zou L, Lu G, et al. PFKP alleviates glucose starvation-induced metabolic stress in lung cancer cells via AMPK-ACC2 dependent fatty acid oxidation. Cell Discov 2022;8:52. [Crossref] [PubMed]
45. Sun MW, Lee JY, de Bakker PI, et al. Haplotype structures and large-scale association testing of the 5’ AMP-activated protein kinase genes PRKAA2, PRKAB1, and PRKAB2 [corrected] with type 2 diabetes. Diabetes 2006;55:849-55. Erratum in: Diabetes 2006;55:1904.
46. Rana R, Shearer AM, Fletcher EK, et al. PAR2 controls cholesterol homeostasis and lipid metabolism in nonalcoholic fatty liver disease. Mol Metab 2019;29:99-113. [Crossref] [PubMed]
47. Cui Z, Zheng C, Lin Y, et al. Lactate Dehydrogenase C4 Accelerates Triple-Negative Breast Cancer Progression by Promoting Acetyl-CoA Acyltransferase 2 Lactylation to Increase Free Fatty Acid Accumulation. Adv Sci (Weinh) 2025;12:e11849. [Crossref] [PubMed]
48. Sun J, Ding J, Yue H, et al. Hypoxia-induced BNIP3 facilitates the progression and metastasis of uveal melanoma by driving metabolic reprogramming. Autophagy 2025;21:191-209. [Crossref] [PubMed]
49. Raimundo N, Song L, Shutt TE, et al. Mitochondrial stress engages E2F1 apoptotic signaling to cause deafness. Cell 2012;148:716-26. [Crossref] [PubMed]
50. Wu H, Wang Y, Lin Y, et al. Potential treatment approaches for malignant peritoneal mesothelioma: in vivo and in vitro experimental study of natural killer cell immunotherapy. Cancer Biol Med 2024;21:1078-94. [Crossref] [PubMed]
51. Shafer-Weaver KA, Sayers T, Kuhns DB, et al. Evaluating the cytotoxicity of innate immune effector cells using the GrB ELISPOT assay. J Transl Med 2004;2:31. [Crossref] [PubMed]
52. Kazim M, Ganguly A, Malespini SM, et al. Granzyme B-Targeting Quenched Activity-Based Probes for Assessing Tumor Response to Immunotherapy. J Am Chem Soc 2025;147:37926-39. [Crossref] [PubMed]
53. Zhang Y, Zou J, Chen R. An M0 macrophage-related prognostic model for hepatocellular carcinoma. BMC Cancer 2022;22:791. [Crossref] [PubMed]
54. Liu J, Chen H, Qiao G, et al. PLEK2 and IFI6, representing mesenchymal and immune-suppressive microenvironment, predicts resistance to neoadjuvant immunotherapy in esophageal squamous cell carcinoma. Cancer Immunol Immunother 2023;72:881-93. [Crossref] [PubMed]
55. Cao R, Zhang W, Zhang H, et al. Comprehensive Analysis of the PRDXs Family in Head and Neck Squamous Cell Carcinoma. Front Oncol 2022;12:798483. [Crossref] [PubMed]
56. Prochazkova K, Finke A, Tomaštíková ED, et al. Zebularine induces enzymatic DNA-protein crosslinks in 45S rDNA heterochromatin of Arabidopsis nuclei. Nucleic Acids Res 2022;50:244-58. [Crossref] [PubMed]
57. Lu S, Gao C, Wang Y, et al. Phylogenetic Analysis of the Plant U2 snRNP Auxiliary Factor Large Subunit A Gene Family in Response to Developmental Cues and Environmental Stimuli. Front Plant Sci 2021;12:739671. [Crossref] [PubMed]
58. Cui Y, Chen H, Xi R, et al. Whole-genome sequencing of 508 patients identifies key molecular features associated with poor prognosis in esophageal squamous cell carcinoma. Cell Res 2020;30:902-13. [Crossref] [PubMed]
59. Otvos B, Alban TJ, Grabowski MM, et al. Preclinical Modeling of Surgery and Steroid Therapy for Glioblastoma Reveals Changes in Immunophenotype that are Associated with Tumor Growth and Outcome. Clin Cancer Res 2021;27:2038-49. [Crossref] [PubMed]
60. Jin X, Yan J, Chen C, et al. Integrated Analysis of Copy Number Variation, Microsatellite Instability, and Tumor Mutation Burden Identifies an 11-Gene Signature Predicting Survival in Breast Cancer. Front Cell Dev Biol 2021;9:721505. [Crossref] [PubMed]