Abstract

Background. Cervical cancer (CC) is one of the most frequent female malignancy. Cancer stem cells (CSCs) positively affect survival outcomes in cancer patients, but in cervical cancer, the mechanism of tumor stem cells is still uncertain. Methods. RNA-seq data and related clinical follow-up of patients suffering from CC were from TCGA. Consensus clustering screened prognostic mRNAsi-related genes and identified molecular subtypes for CC. Based on the overlapping differentially expressed genes (DEGs) in subtypes, we employed LASSO and multivariate Cox regression to screen prognostic-related genes and established the RiskScore system. The patients were grouped by RiskScore, the prognosis was analyzed by the Kaplan-Meier (K-M) curve among the various groups, and the precision of the RiskScore was assessed by the ROC curve. Finally, the potential worth of RiskScore in immunotherapy/chemotherapy response was assessed by evaluating TIDE scores and chemotherapy drug values. Results. We noticed that patients with low mRNAsi had a shorter survival and then identified three molecular subtypes (C1-3), with the C1 having the worst prognosis and the lowest mRNAsi. Finally, we identified 7 prognostic-related genes (SPRY4, PPP1R14A, MT1A, DES, SEZ6L2, SLC22A3, and CXCL8) via LASSO and Cox regression analysis. We established a 7-gene model defined RiskScore to predict the prognosis of CC patients. K-M curve indicated that low RiskScore patients had improved prognosis, and ROC curves indicated that RiskScore could precisely direct the prognostic evaluation for those suffering from the cancer. This was also confirmed in the GSE44001 and GSE52903 external cohorts. Patients were more sensitive to immunotherapy if with low RiskScore, and RiskScore exhibited precise assessment ability in predicting response to immunological therapy in CC patients. Conclusion. CC stemness is associated with patient prognosis, and the RiskScore constructed based on stemness characteristics is an independent prognostic index, which is expected to be a guide for immunotherapy, providing a new idea for CC clinical practice.

1. Introduction

Although cervical cancer (CC) was preventable and treatable, it remains a global burden on female’s health. In Africa, South America, Southeast Asia, and other countries, the incidence and mortality rates are the highest, reaching 15.3-40.1% and 7.7-28.6%, respectively [1]. It was believed that effective implementation of early screening, HPV vaccination, and other means can reduce the incidence of CC and that patients who receive treatment early had higher overall survival (OS) rate [2, 3]. However, in some developing countries, most patients were diagnosed with stage III/IV [46]. Some research data showed that the clinical outcomes in late-stage CC patients were still not optimistic [7]. Consequently, discovering serviceable therapeutic targets and prognostic biomarkers for CC is highly necessary.

Growing evidence suggested that tumor-infiltrating immune cells (TLIs) from the tumor microenvironment (TME) were associated with the progressiveness, aggressiveness, and treatment responsiveness of tumor [8, 9]. Recent studies confirmed that immune cell infiltration could predict CC recurrence, prognosis, and treatment response. For example, Kawachi et al. [10] demonstrated that tumor-associated CD204+ M2 macrophages are predictive metrics for cervical adenocarcinoma. Zhang et al. [11] assessed tumor-infiltrating biomarkers (PD-1, PD-L1, CD3, CD4, CD8, CD20, CD56, and CD68) in CC using immunohistochemistry and found that increased CD4, CD8, CD20, and CD56 signal predicts favorable neoadjuvant chemotherapy response. Someya et al. [12] also confirmed that CD8 or FoxP3 combined with EQD2 may help predict the treatment effect of CC radiotherapy, which was found to help optimize the selection of individualized treatment for CC. In addition, immunotherapies, which blocked immune checkpoints like CTLA-4 and PD-1/PD-L1, gradually have become second-line therapies in malignancies [13, 14]. Some immune components of TME can suppress antitumor immune responses, not all patients will benefit from immunotherapy [15]. Therefore, exploring useful therapeutic targets and prognostic biomarkers for CC is highly critical.

Cancer stem cells (CSCs) were associated with adverse clinical outcomes of cancer. CSCs were a subpopulation of tumors, mainly responsible for tumor maintenance and spread. CSCs would produce many differentiated cells with high value-added and self-renewal capacity and were the main components of tumor populations [16]. Activation of CSCs was known to lead to tumor recurrence and chemotherapy resistance [17]. In recent years, stemness indices were conducted to profile stemness characteristics. For example, Malta et al. [18] used a machine learning algorithm to calculate the stemness index of a pluripotent stem cell sample, which in turn assessed the degree of oncogenic dedifferentiation. Currently, some mRNA expression-based stemness indices (mRNAsi) were associated with low-grade glioma and gastric cancer prognosis, providing new ideas for prognosis, recurrence, and metastasis prediction of gastric cancer and low-grade glioma [19, 20]. However, in CC, the association between stemness and prognosis has not been elucidated. Therefore, this study utilized the mRNA expression profile of CC in TCGA database to focus on the potential prognostic value of mRNAsi and stemness index-related genes affecting its prognosis using bioinformatics algorithms. We attempt to establish a clinical prediction model for CC prediction of prognosis and biomarkers of immunotherapy response to provide a basis for clinical guidance.

2. Materials and Methods

2.1. Data Source and Processing

The high-throughput RNA-seq data (TCGA-CESC) and clinical information of CC were from TCGA database. In this study, we followed the methods of Li et al. [21] to process RNA data. Samples without survival time, status, or clinical follow-up information were removed. After screening, 288 primary tumor samples were retained for subsequent studies in TCGA-CESC cohort. In addition, the GSE44001 and GSE52903 datasets were from the GEO database as validation cohorts. Similarly, the same treatment was done in the GSE44001 and GSE52903 cohorts, and 300 and 55 samples were retained, respectively.

2.2. Acquisition of Stemness Index Based on RNA-Seq

Stemness index model trained from the Progenitor Cell Biology Consortium database with OCLR algorithm was applied for tumor stemness calculation [18, 22]. With stemness index being a value between 0 (the lowest) and 1 (the highest), how similar tumor cells are to stem cells can be measured by the index. A stronger stem cell property was related to a stemness index close to 1. From a previously reported study, we acquired stem cell indices based on the transcriptome of each CESC sample [23]. In the following sections, the index referred to the mRNA expression-based stemness index (mRNAsi).

2.3. Identification of mRNAsi-Related Genes

The Spearman correlation between each mRNAsi and the corresponding gene was calculated in TCGA-CESC cohort, and genes whose correlation index and were seen as mRNAsi-related genes (mRNAsiRGs). Then, univariate Cox regression analysis combined with the clinical information of CESC patients was conducted on the screened mRNAsiRGs to identify which could affect the overall survival (OS) of CESC.

2.4. Consensus Clustering

These genes, which had been found to be highly correlated with the mRNAsi, were used to perform the consensus clustering [24]. Patients were clustered into different stemness subtypes. The cumulative distribution function (CDF) and consensus matrices determined the optimal subgroup number.

2.5. Stemness Subtype Comparison

Chi-squared test was conducted to compare different clinical variables combined TNM stage, stage, grade, age, and stemness in the subtypes for difference exploration among the different stemness subtypes of TCGA datasets. Next, molecular characterization information of TCGA-CESC was obtained from Thorsson et al.’s pan-cancer landscape study [23]. Differences in the genomic landscape among the three subtypes were analyzed. Finally, we counted the gene mutations in different subtypes.

2.6. DEG Analysis and Functional Enrichment Analysis

R package limma [25] was used to calculate DEGs between the C1 vs. C2&C3 (C1) group, the C2 vs. C1&C3 (C2) group, and the C3 vs. C1&C2 (C3) group ( and ). Subsequently, GO and KEGG functional enrichment analyses were performed on DEGs in groups C1, C2, and C3 using the R package clusterProfiler [26], and the pathway with was selected as the most significantly enriched pathway.

2.7. Construction and Verification of the Stemness-Related Prognosis Model

Firstly, the DEGs in groups C1, C2, and C3 were taken as a concurrent set, and then, univariate Cox regression analysis was performed with R package survival to select the DEGs that had an impact on the OS of CC patients. Next, in order to reduce the model complexity, the R package glmnet [27] was used for LASSO Cox regression analysis, and tenfold cross-validation selected appropriate penalty parameter lambda to remove the genes with strong correlation in the model to reduce the dimension. Finally, multivariate Cox regression analysis was performed using R package survival to screen DEGs affecting OS in CC patients and construct a risk assessment model named RiskScore calculated by the following formula: where was the Cox regression coefficient normalized by the -score of the gene and was the expression data of the gene. We used the best cutoff value found by the surv_cutpoint function of the R package survminer for sample division into different groups. The surv_cutpoint function could calculate the statistical value of all potential subgroups in the sample, and in this study, we chose the minimum value as the best subgroup cutoff value [28]. K-M curves were employed to compare the OS difference. The timeROC package [29] was used to perform ROC analysis and calculate the AUC value to evaluate the prognostic ability of RiskScore. Finally, to verify the robustness of RiskScore, it was tested in the GSE44001 and GSE52903 cohorts.

2.8. Correlation between the RiskScore and Clinical Characteristics

We evaluated the distribution of RiskScore in different clinical features including TNM stage, grade, stage, and age. To assess the robustness and accuracy of the prognostic model, survival differences were compared for the two groups stratified by clinical characteristics.

2.9. GSEA Enrichment Analysis

To explore the biological processes associated with RiskScore, we performed the GSEA enrichment analysis using the fgsea package, and clusterProfiler package was utilized for functional annotation.

2.10. Determination of Tumor-Infiltrating Immune Cells

The CIBERSORT algorithm [30] was utilized to compute the relative abundance of 22 common tumor-infiltrating immune cells in TCGA-CESC cohort. Next, immune cell infiltration was assessed using the ESTIMATE algorithm [31] to compute the StromalScore, ImmuneScore, and ESTIMATEScore in different subtype samples. To further elucidate their TME differences, we downloaded 29 TME signatures from previous studies [32] that could predict immune therapy response using R package GSVA [33] to obtain the ssGSEA score.

2.11. Prediction of Treatment Response

We from cancer drugs sensitivity genomics (GDSC, http://www.cancerrxgene.org) downloaded about 1000 cancer drug sensitivity data, using the cancer cells from the CESC-related department. The half-maximal inhibitory concentration () is an important index to evaluate the effect of drugs or the treatment response of samples. The lower the , the stronger the antitumor ability [34]. Using the pRRophetic package, of the chemotherapy drugs was calculated by ridge regression [35]. Next, tumor immune dysfunction and exclusion (TIDE) scores of CC patients were downloaded. TIDE score is an index to assess the response to immunotherapy response. The higher the TIDE score, the lower the benefit of immunotherapy [36]. Finally, to verify the predictive value of RiskScore in immunotherapy response, we performed validation in the IMvigor210 cohort.

3. Results

3.1. mRNAsi Is Associated with CC Prognosis

Initially, we evaluated the correlation of clinicopathological features between mRNAsi in CESC patients in TCGA-CESC cohort. There was no significant correlation between mRNAsi and age, T stage, grade, stage, and N stage (Figure 1(a)). Next, CC patients in TCGA-CESC cohort were clustered in the high-mRNAsi group and the low-mRNAsi group. To evaluate the effect of mRNAsi on OS in CESC patients, K-M analyses showed a vitally worse OS of the low-mRNAsi group than the high-mRNAsi group (Figure 1(b)). Besides, we also compared whether there were differences in mRNAsi between different groups of clinicopathological characteristics. Violin plot of mRNAsi differences among different grade, age, TNM stage, and stage subgroups was exhibited (Figure 1(c)). Combined with the results, no significant difference in mRNAsi was found between different subgroups of TNM stage and age, and the mRNAsi of CESC patients in the G2 and G3 subgroups was significantly higher than that of G1. Further, we found that higher mRNAsi meant poorer OS in both early- and advanced-stage patients (Figure 1(d)). The results suggested that mRNAsi based on CSCs was closely related to their prognosis in CESC.

3.2. Molecular Subtypes of CESC Based on mRNAsiRGs

In total, the expression profiling of 667 mRNAsiRGs in 288 samples from RNA-seq data in TCGA-CESC cohort was obtained by the Spearman correlation analysis. Potential prognostic significance of each mRNAsiRGs was assessed using univariate Cox regression analysis, which screened a sum of 31 significant OS-related mRNAsiRGs from TCGA-CESC cohort (Figure 2(a)). The functional enrichment analysis is shown in Figure S1. The 288 samples were subjected to consensus clustering. Patients were categorized into three significantly distinct subtypes (C1, C2, and C3 group; Figure S2A-B, Figure 2(b)). Those in the C1 demonstrated a significantly worse OS and PFS (Figures 2(c) and 2(d)). The results of PCA analysis for the three subgroups of patients are presented in Figure 2(e). In addition, significant differences were detected in mRNAsi among the three subtypes of patients, with the highest mRNAsi being observed in the C3 subtype (Figure 2(f)). Finally, the heat map demonstrated the difference in the expression levels of 31 mRNAsiRGs in C1-3 subtypes which was shown in Figure 2(g). Finally, in TCGA-CESC cohort, distribution of various clinical features was analyzed in C1-3 subtypes defined in this study to see if there were any differences. We found that TNM stage and age did not differ significantly in the subtypes, but higher stemness features could be clearly observed in patients with C1 and C2 subtypes (Figures 3(a)3(c)).

3.3. Associations between Molecular Subtypes and Genomic Landscapes

Further, to explore the differences in genomic landscape between stemness subtypes. We found that the C3 subtype showed higher levels of TMB and purity, while the C1 subtype showed higher levels of ploidy (Figure 4(a)). Meanwhile, we found that the majority of patients with the C3 subtype belonged to the C2 subtype with IFN-γ dominant characteristics (Figure 4(b)). In addition, we also analyzed the genomic mutation landscape among the molecular subtypes. We found a significant correlation between molecular subtypes and mutations. The waterfall plot showed the integration status of somatic mutations in TCGA-CESC, where Missense_Mutation and Nonsense_Mutation were the most common mutation types, while TTN, MUC4, and EP300 were the genes with high mutation frequency in CESC (Figure 4(c)).

3.4. Immunological Microenvironment in Stemness Subtypes of CESC

To further elucidate the differences between the immune microenvironment with molecular subtypes of patients, we calculated the immune cell infiltration level in the TME of CESC patients using the CIBERSORT and ESTIMATE algorithms. First, according to the CIBERSORT results we observed that the level of partial immune cell infiltration differed between subtypes, where Tregs, T_cells_follicular_helper, T_cells_CD8, and Dendritic_cells_resting were highly infiltrated in the C3 subtype, while Macrophages_M0 and T_cells_CD4_memory_resting were highly infiltrated in the C1 subtype (Figure 5(a)). Meanwhile, ESTIMATE results showed that the StromalScore, ImmuneScore, and ESTIMATEScore of the C1 subtype were significantly higher than those of other subtypes, with higher immune cell infiltration (Figure 5(b)). In addition, to further evaluate TME differences, 29 gene signatures representing tumor main functional components, immune cells, stromal cells, and other cell populations in the TME were downloaded from previous studies, and their TME gene signature ssGSEA scores were calculated. We observed higher ssGSEA scores for features such as angiogenesis fibroblasts (angiogenesis, endothelium, cancer-associated fibroblasts, matrix, matrix remodeling, and granulocyte traffic), protumor immune infiltrate (myeloid cells traffic and immune suppression by myeloid cells), and EMT signature proliferation rate (EMT signature) in the C1 subtype (Figures 5(c) and 5(d)). Finally, we calculated 14 cancer-related pathway activities using the R package PROGENy, and the EGFR, hypoxia, MAPK, TGFb, TNFa, and NFκB pathways were activated in the C1 subtype (Figures 5(e) and 5(f)). Overall, immune cell activity and pathway activity related to tumorigenesis and progression were higher in the TME of C1 subtype patients, which may also be the reason for the shorter OS in C1 subtype patients.

3.5. Differences in Immunotherapy between Molecular Subtypes

Immunotherapy is a promising emerging option in the treatment of CC, supported by several findings [14, 37]. T-cell-inflamed gene expression profile (T-cell-inflamed GEP) was a novel immunotherapeutic biomarker that predicts clinical response to PD-1-targeted immune checkpoint blockade. GEP included genes associated with antigen presentation, and GEP includes IFNγ-responsive genes associated with antigen presentation, chemokine expression, cytotoxic activity, and adaptive immune resistance [38]. Cytolytic activity (CYT) score was also a novel immunotherapy biomarker that could characterize the antitumor immune activity of CD8+ cytotoxic T cells and macrophages [39]. Therefore, we first evaluated the potential efficacy of immunotherapy in patients with different CC subtypes using biomarkers such as T-cell-Inflamed GEP score, Th1/IFNγ gene signature ssGSEA score, and CYT score. In this study, we found that the T-cell-inflamed GEP score, Th1/IFNγ gene signature ssGSEA score, and CYT score of the C3 subtype were higher than those of the C1 and C2 subtypes (Figures 6(a)6(c)). Considering that cancer immunotherapy for immune checkpoint blockade (ICB) was based on the inhibition of key immune checkpoints, we evaluated a representative number of molecules whose expression levels are shown in Figure 6(d). Finally, we found higher MDSC, CAF, exclusion, and TIDE scores in the C1 subtype (Figure 6(e)). These researches indicated that patients with the C3 subtype might respond more favourably to immunotherapy.

3.6. Identification of DEGs Based on Different Subtypes

To identify DEGs in different subtypes, 439 DEGs were identified in group C1 compared to C2&C3, of which 373 were upregulated and 66 were downregulated (Figure 7(a)). 301 DEGs were identified in group C3 compared to C1&C2, of which 280 were downregulated and 21 were upregulated (Figure 7(b)). DEGs could not be identified in the C2 group compared to C1&C3. Further, we analyzed the biological functions of these DEGs. Interestingly, we found that DEGs in the C1 group were mainly enriched in CAF and EMT-related pathways, a result consistent with the TME trend in this study (Figures 7(c) and 7(d)).

3.7. Development and Validation of the Prognostic RiskScore

The prognostic RiskScore was constructed based on DEGs in the subtypes, and 478 DEGs were obtained by taking a concatenation of genes in the C1 and C3 groups. In TCGA-CESC cohort, univariate Cox regression analysis combined with the survival information of patients identified 127 genes affecting patients’ OS (), including 117 risk factors and 10 protective factors. Then, the model was further compressed by LASSO Cox regression analysis to remove the genes with a high goodness of fit. The model reached optimality at , for which we selected 12 genes (ITGA5, SPRY4, P4HA3, PPP1R14A, CXCL8, SLC22A3, TFPI, SEZ6L2, MT1A, ARMCX1, DES, and PLOD2) for subsequent analysis (Figures 8(a) and 8(b)). Finally, based on multivariate COX regression analysis and Akaike information criterion [40], the model showing a minimum AIC was seen as the optimal model. We finally identified seven genes, including SPRY4, PPP1R14A, MT1A, DES, SEZ6L2, SLC22A3, and CXCL8, as mRNAsi-related genes affecting OS in CC (Figure 8(c)). The final prognostic prediction model for CC was constructed: .

The RiskScore of patients in TCGA-CESC cohort was calculated based on the model, and the best cutoff was selected using the R package survminer to divide the sample into a low-risk group () and a high-risk group (). We observed that as the RiskScore increased, more patients died. In addition, risk factors such as SPRY4, PPP1R14A, MT1A, SEZ6L2, SLC22A3, and CXCL8 were significantly highly expressed in the high-risk group, and protective factor DES was obviously highly expressed in the low-risk group (Figure 8(d)). K-M analysis showed that patients with low risk had a better OS than those with high risk (Figure 8(e)), and the ROC curve revealed that AUC of 1-, 3-, and 5-year OS was 0.82, 0.81, and 0.82, respectively (Figure 8(f)). To assess the predictive robustness of the RiskScore, we obtained RiskScores for patients in the independent validation sets GSE44001 and GSE52903 cohorts. Similarly, survival analysis showed that patients with low risk had a better OS than those with high risk. The prediction of 1-, 3-, and 5-year survival probabilities indicates that RiskScore still has excellent AUC values, demonstrating that RiskScore has excellent performance in assessing the prognosis of patients suffering from CC (Figures 8(g)8(j)).

3.8. Clinical Correlation Analysis of the Prognostic RiskScore

In order to determine the relationship between RiskScore and clinicopathological characteristics, interaction between RiskScore and different clinical parameters (N stage, age, and T stage) was discussed. We found increased RiskScore in patients in higher T stage, M stage, and mRNAsi-low groups (Figure 9(a)). Further, we observed a significantly higher RiskScore in patients with C1 subtype than C2 and C3 subtypes (Figure 9(b)). Also, we evaluated the association between RiskScore and molecular subtypes and could observe that the high-risk group was mainly C1 subtype and mRNAsi-low patients (Figure 9(c)). In addition, we performed a survival analysis of high- and low-risk patients in early- and advanced-stage patients. We noted that RiskScore had promising predictive value in different clinical groups (Figure 9(d)). Overall, RiskScore constructed based on prognosis-related mRNAsi features of DEGs in subtypes had good robustness.

3.9. Characteristics of Immune Infiltration/Pathways between the RiskScore Groups

In TCGA-CESC cohort, GSEA enrichment analysis was performed for patients in different risk groups. We found that the high-risk group was significantly enriched in CAF and synthesis and metabolism-related pathways such as KEGG_N_GLYCAN_BIOSYNTHESIS and KEGG_O_GLYCAN_BIOSYNTHESIS and EMT-related pathways such as KEGG_ECM_RECEPTOR_INTERACTION, KEGG_FOCAL_ADHESION, and KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION (Figure 10(a)). We also found that RiskScore was significantly negatively correlated with mRNAsi (Figure 10(b)). TME differences in the two risk groups were also analyzed. First, immune cell infiltration was assessed by ESTIMATE, and the results showed a higher ImmuneScore in the low-risk group (Figure 10(c)). Additionally, relative abundance of 22 immune cells was compared, and we observed that the level of Macrophages_M1, T_cells_CD8, Mast_cells_resting immune cell infiltration, and Dendritic_cells_resting was higher in the low-risk group (Figure 10(d)). The correlation between RiskScore and 22 immune cell infiltration winds is shown in Figure 10(e). To further reveal the TME differences, we found that RiskScore was significantly positively correlated with gene signatures such as angiogenesis fibroblasts and EMT signature proliferation rate (Figure 10(f)). Finally, we also analyzed the relationship between RiskScore and pathway activity scores and found a significant positive correlation between RiskScore and angiogenesis-related pathways (Figure 10(g)).

3.10. Explorations of Clinical Applications for RiskScore

In this study, we examined the validity of the RiskScore for predicting response to immunotherapy by calculating the T-cell-inflamed GEP score, the Th1/IFNγ gene signature ssGSEA score, and the CYT score in the high- and low-risk groups. The results showed that T-cell-inflamed GEP score and Th1/IFNγ gene signature ssGSEA score were significantly increased in the low-risk group (Figures 11(a)11(c)). Subsequently, the expression levels of immune checkpoint molecules in patients in different risk groups were counted. We found that CTLA4, BTLA, HAVCR2, and TIGIT were significantly increased in the low-risk group (Figure 11(d)). Further, we observed that patients in the high-risk group had higher scores of MDSC, CAF, exclusion, and TIDE, indicating that the higher-risk patients were not sensitive to immunotherapy (Figure 11(e)). The correlation between RiskScore and T-cell-inflamed GEP score, Th1/IFNγ gene signature ssGSEA score, immune checkpoint, and TIDE score is shown in Figure 11(f). Furthermore, we examined the ability of RiskScore to predict patient response to ICB therapy to further validate the association of RiskScore with immunotherapy. In the anti-PD-L1 cohort (IMvigor210 cohort), patients in the low-risk group showed significantly prolonged overall survival and clinical benefit (Figures 11(g)11(i)). Notably, patients with low RiskScore also exhibited remarkable survival advantages in the early-stage groups (stage I+II and stage III+IV) (Figures 11(j) and 11(k)). Finally, we examined the validity of RiskScore for chemotherapy response prediction and calculated for each sample. Based on the R package pRRophetic of the pRRophetic algorithm, we analyzed the s of 6 chemotherapeutics (Paclitaxel, Gemcitabine, Cisplatin, Gefitinib, Mitomycin C, and Sunitinib) used in CC treatment and predicted 3 significant correlations between RiskScore and drug sensitivity (Figures 11(l) and 11(m)). The results suggest that RiskScore may also serve as a predictor of immunotherapy.

4. Discussion

CC remains the leading cause of global cancer-related deaths in women, though advances have been made in screening, diagnosis, and treatment. CSCs are highly heterogeneous cell populations, and an in-depth investigation of the interactions between highly heterogeneous CSCs and their TME will help to explore CSC-targeted therapeutic strategies and improve the sensitivity of current CC immunotherapies. This study first explored the relationship between stemness characteristics and tumor prognosis and TME in CC; identified three molecular subtypes including C1, C2, and C3; and constructed a clinical prediction model for predicting prognosis and immunotherapy response.

First, the relationship was analyzed between mRNAsi scores and clinical outcomes, and we found that patients with high mRNAsi had higher OS in both early and more advanced patients in this finding, which has been confirmed in several studies. Zhang et al. [20] found that patients with low mRNAsi had worse OS in low-grade gliomas. These results suggested that mRNAsi was intimately connected to the prognosis of CC patients. Next, based on the mRNAsi-related genes influencing OS of CC, we found that CC could be classified into C1-3 molecular subtypes. Based on stemness characteristics, C3 subtype has the highest mRNAsi and TMB. TMB is an alternative biomarker for neoantigens [41]. Marabelle et al. [42] found that high TMB solid tumor patients were more sensitive to pembrolizumab monotherapy. Although this study did not exactly analyze TMB differences in CC and their impact on prognosis, combined with the results of subtype survival analysis and mRNAsi scores, we speculate that low TMB leads to poor prognosis in CC patients. In addition, we found that TTN, MUC4, and EP300 were the most frequently mutated genes in CC. These genes have been shown to be closely associated with cancer. MUC4 is a classic biomarker that helps tumor cells evade immune cell recognition [43]. Recent studies by Rowson-Hodel et al. [44] show that MUC4 further enhances the survival and metastasis of disseminated tumor cells through the physical interaction of platelets and macrophages. Zhu et al. [45] showed that higher TMB was related to EP300 mutations, which promoted bladder cancer’s antitumor immunity. TTN is a prognostic marker of colon cancer immune microenvironment [46]. And generally speaking, we found for the first time three stem-related subtypes of CC, which provides new insights into the pathological study of CC.

In this study, a clinical prediction model composed of SPRY4, PPP1R14A, MT1A, DES, SEZ6L2, SLC22A3, and CXCL8 genes was constructed. SPRY4 may function as an antitumor gene in CRC, and SPRY4 suppresses CRC progression by inhibiting EZH2 [47]. The pathogenesis of human melanoma involves PPP1R14A, which drives tumorigenesis and Ras activity through the activation of growth-promoting ERM family factor and inhibition of tumor suppressor merlin [48]. DES has been shown to inhibit telomerase activity in prostate cancer cells [49]. Chen et al. [50] found that SEZ6L2 expression level was significantly related to lymph node metastasis, TNM stage, and HER-2 status of breast cancer, and knocking down SEZ6L2 could significantly inhibit breast cancer cell line proliferation. Fu et al. [51] confirmed that SLC22A3 could inhibit ESCC metastasis, whereas dysregulated SLC22A3 promoted filopodia formation and cell invasion through lowering its direct relation to ACTN4, increasing ACTN4 actin-binding activity in normal esophageal cells and ultimately leading to early development and progression of ESCC. It has been reported that increased CXCL8 levels indicate poor clinical prognosis and tumor progression in gastric cancer patients. CXCL8 is largely secreted via macrophages and forms an immunosuppressive microenvironment in gastric cancer by inducing PD-L1+ macrophages [52]. However, the function of these genes has not been elucidated in CC, and also, further experiment is required to verify their function.

In recent years, due to its complex composition and mechanism, TME is critical in tumor progression and tumorigenesis and has become a research hotspot in the field of cancer [53]. Increasing evidence suggested that TICs and stromal components are closely related to breast cancer development [54, 55]. From this report, we found that there were differences in the infiltration levels of T_cells_follicular_helper, Macrophages_M1, and T_cells_CD8 in patients with different RiskScore groups. CD8+ T cells and CD4+ T cells are considered to be the major components of TLIs and are involved in antitumor immunity [56, 57]. Liu et al. suggested that high expression of TIGIT in CD8+ T cells of CC patients promoted their depletion and that blocking TIGIT in activated CD8+ T cells attenuated the inhibitory effect of SHIP-1 on CD8+ T cells and enhanced the activation of NF-κB and ERK [58]. In pancreatic cancer, Sehgal et al. found that CD8+ T cells blocking PD-1 and CXCR4 perform antitumor effects [59]. Tumor-associated macrophages [60] also have a major role in TME, and their different forms contribute to OS in tumor patients. TAM consists of multiple macrophage subpopulations, including M0, M1, and M2 macrophages [61]. M2 macrophages drive immunosuppression and angiogenesis by producing immunosuppressive factors such as IL-10 [62]. In this study, lower RiskScore CC patients possessed improved immune activity and OS, and it was worth being aware that higher RiskScore patients had lower TIDE scores. We speculated that patients with lower RiskScore with TME in an activated state might be more sensitive to immunotherapy. To confirm this conjecture, we validated the predictive performance of RiskScore for immunotherapy response in the IMvigor210 cohort receiving anti-PD-L1 therapy, and the results confirmed our conjecture. These results suggested that higher immunosuppressive cell infiltration in patients with higher RiskScore promoted their resistance to immunotherapy, suggesting that poor prognosis in high-risk patients was connected to the induction of an immunosuppressive microenvironment.

In addition, we found that the high-risk group was mainly enriched in CAF and EMT-related pathways, for example, KEGG_ECM_RECEPTOR_INTERACTION and KEGG_FOCAL_ADHESION. In cancer research, it is always known that most cancer-related deaths resulted from metastasis of cancer cells [63, 64]. ECM-receptor interactions play an essential role in tumor shedding, adhesion, degradation, motility, and proliferation. ECM is found to be upregulated in prostate cancer tissues [65] and to also participate in gastric cancer invasion and metastasis [66]. Furthermore, ECM in colorectal cancer stimulates epithelial-mesenchymal transition (EMT) in cancer cells [66]. And the link between these pathways and stemness in CC needs to be demonstrated by more advanced studies.

In a word, our study identified stemness-associated genes associated with prognosis in TCGA public database and further constructed a 7-gene prognostic signature with a promising ability in predicting survival outcomes and immunotherapy response in CC patients. However, there are some limitations. Firstly, our model was developed based on a public database, lacking clinical samples for validation, and subsequent multicenter, multidata clinical trials are needed. Secondly, we identified some differential pathways in two RiskScore groups, the association in these pathways and tumor stemness was not designed for further exploring the molecular mechanisms in a wet trial, and more systematic and in-depth studies in wet trials are needed to elucidate the association.

5. Conclusion

We utilized mRNAsi-related genes from CSC features to develop new molecular subtypes of stemness signatures in CC, which might expect to complement the existing pathological typing and provide insights into the mechanisms of metastasis and drug resistance in CC. Our research also developed a 7-gene prognostic signature for CC, and there were promising clinical benefits of RiskScore in predicting prognosis and treatment benefits in CC.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no competing interest.

Authors’ Contributions

All authors contributed to this present work: [9] and LY designed the study, and [48] acquired the data. JHZ drafted the manuscript, YYH and [37] revised the manuscript. All authors read and approved the manuscript.

Acknowledgments

The study was funded by the General Program of Chongqing Natural Science Foundation (cstc2021jcyj-msxmX0219): Study on the Role of CircRNA-0008193 Targeting miR-1226-5p in Regulating PD-L1 in Immune Escape of Cervical Cancer.

Supplementary Materials

Supplementary 1. Figure S1: 31 significant OS-related mRNAsiRGs functional enrichment analysis.

Supplementary 2. Figure S2: consistency cluster analysis (A) CDF for to . (B) Relative change in area under the CDF curve according to different values.