A Novel Prognostic Risk Model for Cervical Cancer Based on Immune Checkpoint HLA-G-Driven Differentially Expressed Genes

Human leukocyte antigen G (HLA-G) is a potential checkpoint molecule that plays a key role in cervical carcinogenesis. The purpose of this study was to construct and validate a prognostic risk model to predict the overall survival (OS) of cervical cancer patients, providing a reference for individualized clinical treatment that may lead to better clinical outcomes. HLA-G-driven differentially expressed genes (DEGs) were obtained from two cervical carcinoma cell lines, namely, SiHa and HeLa, with stable overexpression of HLA-G by RNA sequencing (RNA-seq). The biological functions of these HLA-G-driven DEGs were analysed by GO enrichment and KEGG pathway using the “clusterProfiler” package. The protein-protein interactions (PPIs) were assessed using the STRING database. The prognostic relevance of each DEG was evaluated by univariate Cox regression using the TCGA-CESC dataset. After the TCGA-CESC cohort was randomly divided into training set and testing set, and a prognostic risk model was constructed by LASSO and stepwise multivariate Cox regression analysis in training set and validated in testing set or in different types of cervical cancer set. The predictive ability of the prognostic risk model or nomogram was evaluated by a series of bioinformatics methods. A total of 1108 candidate HLA-G-driven DEGs, including 391 upregulated and 717 downregulated genes, were obtained and were enriched mostly in the ErbB pathway, steroid biosynthesis, and MAPK pathway. Then, an HLA-G-driven DEG signature consisting of the eight most important prognostic genes CD46, LGALS9, PGM1, SPRY4, CACNB3, PLIN2, MSMO1, and DAGLB was identified as a key predictor of cervical cancer. Multivariate Cox regression analysis showed that this signature is an independent risk factor for the overall survival of CESC patients. Kaplan-Meier survival analysis showed that the 5-year overall survival rate is 23.0% and 84.6% for the high-risk and low-risk patients, respectively (P<0.001). The receiver operating characteristic (ROC) curve of this prognostic model with an area under the curve (AUC) was 0.896 for 5 years, which was better than that of other clinical traits. This prognostic risk model was also successfully validated in different subtypes of cervical cancer, including the keratinizing squamous cell carcinoma, non-keratinizing squamous cell carcinoma, squamous cell neoplasms, non-squamous cell neoplasms set. Single-sample gene set enrichment (ssGSEA) algorithm and Tumor Immune Dysfunction and Exclusion (TIDE) analysis confirmed that this signature influence tumour microenvironment and immune checkpoint blockade. A nomogram that integrated risk score, age, clinical stage, histological grade, and pathological type was then built to predict the overall survival of CESC patients and evaluated by calibration curves, AUC, concordance index (C-index) and decision curve analysis (DCA). To summarize, we developed and validated a novel prognostic risk model for cervical cancer based on HLA-G-driven DEGs, and the prognostic signature showed great ability in predicting the overall survival of patients with cervical cancer.

Introduction

Cervical cancer ranks as the fourth most prevalent cancer in the world (1). Persistent infection with human papillomavirus (HPV) has been identified as a necessary but not sufficient factor that leads to cervical cancer. More than 85% of the global cervical cancer burden occurs in less developed countries that lack organized screening and HPV vaccination programmes (2). Treatment depends on the disease extent at diagnosis and might involve surgical management for early-stage disease or adjuvant treatment or combination therapy for locally advanced disease. Survival prediction mainly depends on International Federation of Gynecology and Obstetrics (FIGO) stages. The expected 5-year survival rate remains poor (15% or less) for patients with recurrence or metastasis (stage IV) (3). Unfortunately, there are few effective therapeutic strategies specifically for recurrent or metastatic cervical cancer. Strikingly, immunotherapies hold promise for these patients (4). Thus, to identify novel biological markers, potential therapeutic targets and develop effective immunotherapy strategies are urgently needed in clinical settings.

Human leukocyte antigen G (HLA-G) is a tolerogenic nonclassical HLA-class I antigen with immunoinhibitory function. HLA-G gene polymorphism likely determines women at higher risk for cervical HPV infections and/or viral persistence (5). Ectopic expression of HLA-G molecule in cervical lesions helps tumour cells escape immunosurveillance by generating inhibitory signals (6). Increasing evidence suggests that HLA-G expression is strongly associated with high-grade tumour and poor prognosis in patients with cancer, especially lung cancer, colorectal cancer, cervical cancer (7). Based on the restricted tissue expression pattern and immunosuppressive functions of HLA-G molecule, many recent studies have consistently claimed that HLA-G could be a new immune checkpoint molecule in tumours (8–13). Our previous studies support that upregulation of HLA-G expression in the tumour microenvironment plays a key role in early cervical carcinogenesis (5, 6, 14). These results were in good agreement with results from other studies, and HLA-G may have significance in the early diagnosis of cervical malignant lesions (15–17). Notably, biomarkers concerning immune checkpoint blockade (ICB) therapy has become a promising strategy in cancer treatment. Therefore, it is urgent to understand the downstream signal pathway induced by HLA-G expression and the dynamic indicators to monitor the effect of immunotherapy after anti-HLA-G. However, its role in HLA-G-driven differentially expressed genes (DEGs) in tumour growth and advanced-stage cervical cancer remain unknown.

This study aimed to develop and validate a prognostic risk model by screening the immune checkpoint molecule HLA-G and its DEGs in cervical cancer based on RNA sequencing (RNA-seq) and the TCGA dataset by a series of bioinformatics methods.

Materials and Methods

Data Acquisition and Processing

The RNA-seq FPKM (fragments per kilobase of exon model per million mapped reads) data, Methylation450k data of patients with cervical cancer, and corresponding clinical data were obtained from the TCGA (https://portal.gdc.cancer.gov/). After removing normal samples, follow-up data and clinicopathological information were matched with the expression profiles. In the TCGA cervical squamous cell carcinoma (CESC) dataset, a total of 275 patients with a survival status and follow-up of more than 30 days were included (Table 1). The study design was illustrated in Figure 1.

TABLE 1

Table 1 Clinical information of included CESC patients based on the overall survival in TCGA dataset.

FIGURE 1

Figure 1 Flow diagram of this study.

Cell Transfection Assays

The human cervical carcinoma cell lines SiHa and HeLa (ATCC, MD, USA) were cultured in DMEM supplemented with 10% heat-inactivated FBS (Gibco, Grand Island, USA). Cells were cultured at 37°C in a humidified 5% CO₂ atmosphere. We constructed two cervical carcinoma cell lines that stably overexpressed HLA-G. Both SiHa and HeLa were transfected with the recombinant pVITRO2-hygro-mcs vector (In vivoGen, USA) containing HLA-G cDNA or empty vector by Lipofectamine 3000 (InvitroGen, USA). The details of the transfection performance have been described previously (18).

Cell Sequencing Data Analysis of the Transfectants

Total RNA was extracted from SiHa-pVITRO2-mcs, SiHa-HLA-G-pVITRO2-mcs, HeLa-pVITRO2-mcs and HeLa-HLA-G-pVITRO2-mcs cells. RNA sequencing was performed using the BGISEQ-500 platform (BGI, Wuhan, China). The details of the experimental process have been described in the previous study (19). The experimental process is based on the complete genomics sequencing technology and DNA NanoBalls (DNBs) technology for sequencing library construction and combined primer anchor synthesis (cPAS) for sequencing. Before downstream analysis, raw reads containing the adaptors sequences and low-quality reads were filtered out, and clean reads were mapped to reference genes using Bowtie2. Gene expression levels were quantified by the RSEM software package. Uniform manifold approximation and projection (UMAP) were performed to reduce the high-dimensional dataset to a three-dimensional (3D) space using the “umap” R package.

Identification of HLA-G-Driven DEGs

The “DESeq2” R package was performed to screen differentially expressed genes induced by the HLA-G molecule in cervical carcinoma cells by comparing the whole transcriptome expression of SiHa-HLA-G-pVITRO2-mcs versus SiHa-pVITRO2-mcs cells and HeLa-HLA-G-pVITRO2-mcs versus HeLa-pVITRO2-mcs cells. The “ggplot2” and “pheatmap” packages were used to generate volcano and heat maps of HLA-G-driven DEGs. Meanwhile, an UpSet plot was drawn to screen out the upregulated and downregulated overlapping DEGs between the SiHa group and HeLa group (|Log₂ Fold Change (FC) | > 0.5, p-value < 0.001) using the “UpSetR” R package.

Gene Functional Enrichment Analyses

To investigate the biological functions of HLA-G-driven DEGs, functional enrichment was comprehensively detected by gene ontology (GO) terms enrichment and KEGG pathway using the “clusterProfiler” R package. GO enrichment was classified into three categories: biological processes (BP), cellular components (CC), and molecular functions (MF). GO enrichment and KEGG pathway were based on the threshold of p-value < 0.05, and the top ten enrichment items are illustrated in figures.

Protein-Protein Interaction (PPI) Network Construction

Associations between these coding products of HLA-G-driven DEGs were constructed by the STRING database (https://string-db.org/). All p-values < 0.05 were considered significant. A PPI network was constructed based on these HLA-G-driven DEGs under the threshold of a minimum required interaction score > 0.4. Node genes that interact with HLA-G or prognostic model-related genes were retained and visualized using the Cytoscape open-source bioinformatics software platform (20).

Evaluation of Methylation Prognostic Indicators

Methylation data of cervical cancer tissues matched with survival information of patients were used to assess the relationship between DNA methylation levels and overall survival (OS) status by univariate Cox regression analysis. In addition, eight signature-related gene DNA methylation statuses and their expression levels, clinical information and survival status were investigated through MEX. PRESS web sites (http://mexpress.be).

Prognostic Risk Model Construction

The overlapping HLA-G-driven DEGs were analyzed by univariate Cox regression analysis using the “survival” package. Least absolute shrinkage and selection operator (LASSO) analysis was performed to select candidate DEGs with penalty parameter tuning adjusted by 10-fold cross validation using the “glmnet” package. To verify candidate DEGs and the risk model, the 275 enrolled patients in the TCGA-CESC dataset were randomly divided into training set (150 cases) and testing set (125 cases) using the “caret” package. The optimal prognostic signature was constructed based on the LASSO filtered candidate DEGs by the following stepwise multivariate Cox analysis in the training set using the “My.stepwise” package. Then, the testing set and entire set were used for validation. The risk score was calculated according to the following formula:

$$Risk score={\displaystyle \sum _{i=1}^N\text{βi}\left(\text{mRNAi}\right)\ast \text{expression}\left(\text{mRNAi}\right)}$$

(N represents the number of gene signatures; βi represents the coefficients of each gene).

Prognostic Signature for Cervical Cancer

CESC patients were divided into low-risk group and high-risk group according to the risk score using the surv_cutpoint function in the “survminer” package. Then, the risk score distribution, the survival status and survival time for each patient, as well as the heat map of the risk-related genes expression levels in each patient were plotted in the training set, testing set and entire set. Kaplan-Meier (KM) curves depicting overall survival were performed to investigate the survival difference between the above mentioned low-risk group and high-risk group. The predictive performance of the m6A risk signature for 1-, 3-, and 5-year OS was assessed by receiver operating characteristic (ROC) curves carried out using the “timeROC” package and concordance index (C-index) analyses were performed and validated in the training set, testing set, and entire set. The area under the curves (AUC) on the ROC were calculated.

Validation of the Risk Score Signature for Subtypes of Cervical Cancer

The clinical data of cervical cancer subtypes were obtained from TCGA. Keratinizing squamous cell carcinoma (SCC) or non-keratinizing SCC, squamous cell neoplasms or non-squamous cell neoplasms were grouped and used to verify the corresponding risk model. Similarly, the prediction efficiency of the risk score system was evaluated in different subtypes of cervical cancer.

Nomogram Construction and Validation

Prognostic factors such as age, clinical stage, histological grade, pathological types, and risk score were used for the selection of the prognostic parameters to construct the nomogram through LASSO. The nomogram was constructed to assess the probability of 1-, 3-, 5- and 10-year overall survival with calibration curves for CESC patients using the “rms” R package. Then, the performance of the nomogram was visualized by AUC, concordance index, and DCA results.

Evaluation of Immune Status Between Low- and High-Risk Groups

To investigate the relationships between the immune system and the prognostic signature, we further analysed the immune status between patients in low-risk group and high-risk group. The immune activity between the two groups was quantified based on single-sample gene set enrichment analysis (ssGSEA) by exploring the signatures of immune-related gene. Then, the relationships between the risk score and immune cells or relationships between the risk group and immune gene set expression level were further analysed.

Immunotherapy Responsiveness Evaluation of Risk Signature and Related Genes

The Tumor Immune Dysfunction and Exclusion web (TIDE, http://tide.dfci.harvard.edu/login/) was used to predict the immunotherapy responsiveness of risk signature and related genes. TIDE is based on tumor pre-treatment expression profiles and can be used to estimate multiple published transcriptomic biomarkers to predict patient response. The published biomarkers (TIDE, MSI.score, TMB, CD274, CD8, IFGN, T.clonality, B.clonality, Merk18) based on their predictive power of response outcome were used for comparison.

Clinical Benefits Evaluation of CESC OS Signatures

Recent years the net reclassification index (NRI), integrated discrimination improvement (IDI) have become the two popularity alternatives assessment in risk evaluation or prediction between different models (21, 22). We searched ten recently published CESC prognostic signatures, and used NRI, IDI method in “survIDINRI” package to compare the clinical benefits of our signature and published signatures (23–32).

Results

Acquisition of HLA-G-Driven DEGs in Cervical Cancer

In this study, HLA-G-driven DEGs were screened out by comparing different gene expression profiles of SiHa-HLA-G-pVITRO2-mcs versus SiHa-pVITRO2-mcs cells and HeLa-HLA-G-pVITRO2-mcs versus HeLa-pVITRO2-mcs cells (Supplementary Figure 1). Read quality statistics after filtering by RNA-seq analysis are shown in Figure 2. After normalization (Figure 2A), the expression level of 1108 HLA-G-driven DEGs, as shown in heatmaps, had clear clusters within different experimental groups (Figure 2B). UMAP a non-linear dimension reduction tool uses a “force-based layout” algorithm to built a weighted graph represent data at low dimensions (33). The UMAP analysis of the resultant total expression of different cell groups revealed clear spatial partitioning of each cell group (Figure 2C). The density plot also showed the difference expression value distribution of SiHa and HeLa cells with or without overexpression of HLA-G (Figure 2D). We identified 1900 upregulated and 2774 downregulated genes in the SiHa group (Figures 3A, C) and 3231 upregulated and 3436 downregulated genes in the HeLa group (Figures 3B, C), with thresholds of |Log₂ FC| > 0.5 and p-value < 0.001. Finally, we selected the overlapping genes in the upregulated group and downregulated group, including 391 upregulated genes and 717 downregulated genes, which were identified as HLA-G-driven DEGs in cervical cancer (Figures 3A–C).

FIGURE 2

Figure 2 Read quality statistics after filtering by RNA-seq analysis. (A) The gene expression normalization of SiHa-pVITRO2-mcs, SiHa-HLA-G-pVITRO2-mcs, HeLa-pVITRO2-mcs and HeLa-HLA-G-pVITRO2-mcs cells after being sequenced by the BGISEQ-500 platform and quantified by the RSEM software package. (B)Heatmaps draw by “pheatmap” of 1083 genes that significantly co-expressed with HLA-G in SiHa and HeLa cells. (C) UMAP plotting of SiHa-pVITRO2-mcs, SiHa-HLA-G-pVITRO2-mcs, HeLa-pVITRO2-mcs and HeLa-HLA-G-pVITRO2-mcs cells with 1083 genes that significantly co-expressed with HLA-G. (D) The expression values density of SiHa-pVITRO2-mcs, SiHa-HLA-G-pVITRO2-mcs, HeLa-pVITRO2-mcs and HeLa-HLA-G-pVITRO2-mcs by “plotDensities” function in “limma” package.

FIGURE 3

Figure 3 Comparison of gene expression profile which HLA-G induced of cervical cancer. Volcano of significantly DEGs in SiHa (A) and HeLa (B), (C) Venn diagram, (D) the bubble plot of KEGG analysis, (E–G) the bubble plot of GO functional enrichment analysis. BP: Biological Precess. CC, Cellular Components. MF, Molecular Function.

Functional Annotations of HLA-G-Driven DEGs

The 1108 HLA-G-driven DEGs were mainly enriched in 559 GO terms and 68 KEGG pathways according to the functional enrichment analysis (p-adjust < 0.05). The KEGG results showed that the top five significantly enriched pathways were ErbB signalling pathway, steroid biosynthesis, MAPK signalling pathway, longevity regulating pathway-multiple species, and EGFR tyrosine kinase inhibitor resistance (Figure 3D). Meanwhile, the results indicated by GO enrichment showed that the top five significantly enriched terms in the biological process (BP) group were regulation of small molecule metabolic process, regulation of lipid metabolic process, connective tissue development, regulation of lipid biosynthetic process, and regulation of steroid metabolic process. In the cellular component (CC) group, the top five significantly enriched terms were coated vesicle membrane, coated vesicle, cell-substrate junction, cell- substrate adherens junction, and clathrin-coated vesicle membrane. In the molecular function (MF) group, the top five significantly enriched terms were protein serine/threonine kinase activity, DNA-binding transcription activator activity, RNA polymerase II-specific, epidermal growth factor receptor binding, growth factor binding, and protein phosphorylated amino acid binding (Figures 3E–G). These functional annotations suggested possible pathways or functional imbalances in which HLA-G-driven DEGs might be involved in the process of cervical cancer.

Construction of PPI Network

We constructed a PPI network that incorporated 1097 nodes and 4072 edges based on the data from the STRING database (PPI enrichment p-value < 0.0001). The co-expression network was processed via Cytoscape software to identify possible key modules, and 375 genes that were directly or indirectly related to HLA-G and model-related genes were retained. A denser PPI network containing 375 nodes (34.2% of total) and 3062 edges (75.2% of total) is shown in Figure 4. HLA-G, CD46, PGM1, CACNB3, PLIN2, and MSMO1 directly interacted with several genes, which indicated that the model-building genes may involve various biological functions.

FIGURE 4

Figure 4 PPI network based on the data from STRING database.

Construction of a Prognostic Model Based on HLA-G-Driven DEGs

In a biological context, LASSO is widely used for gene expression and genomic data variable selection to identify a parsimonious model—a minimal subset of variables required to explain the data. There are relatively a few detailed reports of prognostic models based on gene expression profiles in cervical cancer (26, 28–30, 32). In this study, we selected the CESC prognostic markers from HLA-G driven DEGs buy LASSO regression. First, we integrated clinical data such as survival status, follow-up data, and clinicopathological information and matched them with the mRNA expression profiles of patients in the TCGA-CESC dataset. A total of 275 patients in the TCGA-CESC dataset with complete overall survival data. A total of 171 genes related to overall survival were screened in the TCGA-CESC cohort by univariate Cox regression analysis. Twenty-eight genes selected from LASSO (Supplementary Figure 2) were then used to identify key predictors of HLA-G-driven DEGs in the prognostic model for predicting the overall survival of cervical cancer patients by stepwise multivariate Cox regression. Then, eight signature-related genes are shown in the forest plots of Figure 5, connecting HLA-G-driven DEGs to CD46, LGALS9, PGM1, SPRY4, CACNB3, PLIN2, MSMO1 and DAGLB in the OS model. These eight genes all showed a significant relationship with overall survival in both univariate (Figure 5A) and multivariate (Figure 5B) Cox regression analyses. Among them, only high expression of LGALS9 was related to a favourable prognosis in patients of cervical cancer (P=0.0072). According to the coefficient of each gene in the stepwise multivariate Cox (Figure 5C Beta), the risk score for this prognostic OS model=(1.10*expression level of CD46) + (-0.30*expression level of LGALS9) + (0.73*expression level of PGM1) + (0.41*expression level of SPRY4) + (0.47*expression level of CACNB3) + (0.27*expression level of PLIN2) + (0.52*expression level of MSMO1) + (0.63*expression level of DAGLB). As depicted in Figure 5D, we observed a significant correlation between the risk score signature and overall survival (P<0.001), which indicated that the risk score of this prognostic model could be used as one of the independent risk factor for cervical cancer patients. The KM curve and ROC curve of risk-related genes were also analysed. As shown in Supplementary Figure 3, the p-value of KM for all genes was less than 0.05. CD46, PGM1, DAGLB, and SPRY4 had relatively large AUC values for the ROC curves.

FIGURE 5

Figure 5 Forest plots presenting the univariate and multivariate Cox regression analysis of prognostic HLA-G-driven DEGs in model for overall survival (OS). (A) Forest plot for the univariate Cox analysis of risk score, risk signature related genes(CD46, LGALS9, PGM1, SPRY4, CACNB3, PLIN2, MSMO1, and DAGLB) and other clinicopathological factors: Pathologic_M, Pathologic_N, Pathologic_T, Stage and age in overall survival (OS) of cervical Cancer. (B) Forest plot for the relationship between risk signature related genes (CD46, LGALS9, PGM1, SPRY4, CACNB3, PLIN2, MSMO1, DAGLB) and other clinicopathological factors: Pathologic_M, Pathologic_N, Pathologic_T, Stage age in overall survival (OS) prediction of cervical Cancer. (C) Forest plot for the multivariate Cox analysis of signature related related genes in overall survival (OS) prediction of cervical Cancer. (D) Forest plot for the multivariate Cox analysis of risk signature and other clinicopathological factors: Pathologic_M, Pathologic_N, Pathologic_T, Stage age in overall survival (OS) prediction of cervical Cancer. (CI, confidence interval; HR, hazard ratio).

Application and Validation of Prognostic Risk Model

After we obtained a relatively fixed cut-off value of the risk score by the “Surv_cutpoint” function through the training set to evaluate the testing set and the entire set. According to the fixed risk score cut-off value of 1.844096, patients were divided into low-risk scores group and high-risk scores group. The risk score distribution of each patient showed a high correlation with each patient’s survival status in the TCGA-CESC dataset, as shown in the dot plot figure (Figures 6A, B). Kaplan-Meier survival analysis showed that the OS time of the low-risk scores group was significantly longer than that of the high-risk scores group. The 5-year overall survival rates were 84.6% and 23.0% for the low- and high-risk scores patients, respectively (Figure 6C). The relationship between the relative expression level of eight signature-related genes and the relative risk score groups was shown in the heatmap (Figure 6D). In addition, the ROC curves based on the risk score were plotted, and the AUC signatures were 0.824 at 1 year, 0.837 at 3 years, 0.896 at 5 years, and 0.866 at 10 years (Figure 7A). We observed that the risk score signature had better predictive ability than other clinical factors (age, stage, pathological types, etc.) at 1, 3, 5, and 10 years (Figures 7D–H). In addition to the training set, we used the testing set and entire set of data to validate the prognostic value of this risk model. The results of the ROC curves in both the testing set and entire set were similar to those in the training set (Figures 7B, C), which indicated that the risk score signature has better sensitivity and specificity. In addition to the ROC curve, Kaplan-Meier survival subgroup analysis further validated the value of prognostic risk models in age (<45, 45-69 and >69), stage (stage I, II, III, IV) or pathological type (T, N, M).

FIGURE 6

Figure 6 Risk score analysis, Kaplan-Meier analysis for the validation of prognostic model in training set, testing set and entire set. (A) Rank of risk score and distribution of groups. Patients were divided into low-risk group and high-risk group based on the median value of the risk score calculated. (B) The survival status and survival time of patients ranked by risk score. (C) Kaplan-Meier suggested that high-risk score group had shorter overall survival than low-risk score group in TCGA-CESC cohort. (D) The heatmap of the expression of eight HLA-G-driven DEGs in the signature in low- and high-risk score groups.

FIGURE 7

Figure 7 ROC Curves of OS for this prognostic model. (A) training set, (B) testing set, (C) entire set, (D) age, (E) stage, (F) pathologic_T, (G) pathologic_N, (H) pathologic_M.

Cervical cancer staging classification can stratify patient survival and help guide management more accurately. Multiple validation trials have demonstrated that stratified management of patients improves 5-year overall and disease-free survival. Given that stage T, N, M and ages were predictors of prognosis in cervical cancer patients, we compared the prognostic performance of these clinicopathological features with risk groups basted on our signature by survival curves. As depicted in Figure 8, the low risk patients owned relative long OS compared with high risk patients in the same clinicopathological grouping. In stage, Stage IV patients in low risk groups showed little difference OS with Stage I patients in high risk (Figure 8A). Similar significant different results of high and low risk groups were also found in different pathological M, pathological N, stage and age (Figures 8B-D). In pathological T, only T4 showed no relative significance, it may be because the patient’s T4 stage has too much influence on the patient, or the number of patients in the included T4 stage is too small (Figure 8E). Both these Kaplan-Meier curves and ROC curves suggested that this prognostic risk model constructed in this study could be used as a valuable prognostic tool to predict the OS of cervical cancer patients.

FIGURE 8

Figure 8 Kaplan-Meier Curves of OS for this prognostic model. (A) age, (B) stage, (C) pathologic_T, (D) pathologic_N, (E) pathologic_M.

Validation of the Prognostic Model in Subtypes of Cervical Cancer

To validate the prediction efficiency of the risk score signature of this corresponding risk model in different subtypes of cervical cancer. This prognostic model was also successfully validated in squamous cell neoplasms, non-squamous cell neoplasms, keratinizing SCC, and non-keratinizing SCC (Figure 9). The results of Kaplan-Meier analysis showed that the overall survival time of CESC patients with low-risk scores was significantly longer than that of patients with high-risk scores in those subtypes (P<0.001) (Figures 9A–D). The AUC values at 5 years were 0.827 for squamous cell neoplasms (Figure 9E), 0.92 for non-squamous cell neoplasms (Figure 9F), 0.954 for keratinizing squamous cell carcinoma (Figure 9G), and 0.604 for non-keratinizing squamous cell carcinoma (Figure 9H). These results indicate that this prediction model has a certain predictive effect on most CESC tumour types.

FIGURE 9

Figure 9 Validate the prediction efficiency of the risk score signature in different subtypes of cervical cancer. (A, E) squamous cell neoplasms, (B, F) non-squamous cell neoplasms, (C, G) keratinizing squamous cell carcinoma, (D, H) non-keratinizing squamous cell carcinoma.

Construction of a Nomogram Based on the Eight Genes

After LASSO regression was used to test the significance of those traditional clinical prognostic parameters (age, sex, stage, grade, and pathological types) with the risk score signature (Figures 10B, C). A quantitative nomogram for CESC prognosis integrated the eight-gene signature, and the clinical prognostic parameters were established (Figure 10A), integrating the pathological TNM, stage, and risk score based on the LASSO results. The calibration curves after patients grouped into intervals containing 80 subjects on average were plotted to validate the nomogram performance. In the calibration plots, the predictive curves were close to the ideal curve for 1-, 3-, 5- and 10-year OS, indicating good performance (Figure 10D). Moreover, according to the C-index plot, at all inspection times, the nomogram curve was above the curve of tumour pathological TNM stage at all times. During most of the testing time, the nomogram curve was above the curve of the risk score in the total TCGA cohort (Figure 10E). These results combined with the ROC results (Figure 10F) indicated that the nomogram predicted survival better than a single prognostic factor. The DCA results showed that in a large risk threshold range, the nomogram and the risk score curve had higher benefits than the extreme curve, so their optional risk threshold ranges were relatively large and relatively safe (Figure 10G). This means that the the constructed nomogram has more clinical usefulness in predicting survival probability of cervical cancer patients. If the threshold probability was above 10%, using this nomogram to predict prognosis added more benefit than either the treat-none scheme or the treat-all scheme.

FIGURE 10

Figure 10 The nomogram model. (A) Prognostic nomogram for patients with cervical cancer, (B, C) Construction of prognostic nomogram models, (D) Calibration curves, (E) The calculation of the C-indexes in nine models, (F) ROC curve of overall survival for eight genes signature and clinical parameters, (G) Decision curve analysis.

Evaluation of the Relationship Between the HLA-G-Driven DEGs Signature and Immune Status in CESC

Immune infiltration in tumour microenvironment plays a key role in tumour progression and disease prognosis. We investigated the expression profiles of the 29 immune-related gene sets in the TCGA-CESC cohort. As depicted in Figure 11A, the low-risk group showed more immune activities, which could be caused by the higher infiltrated immune cells. Furthermore, the immune status of each patient in the low-risk group and high-risk group showed different degree of heterogeneity (Figure 11B). To understand the relationship between the risk model and the immune microenvironment more specifically, we analysed the risk score and the corresponding immune microenvironment score directly from each patient (Figures 11C–K). We found that the infiltration of CD8+ T cells, HLA, APC co-inhibition, and T cell co-stimulation was less in the high-risk group than in the low-risk group in the OS model. These results suggested that the eight genes in this prognostic model may play an important role in the immune regulation of the tumour microenvironment.

FIGURE 11

Figure 11 The low-risk group and high-risk group showed different immune status. (A, B) the expression profiles of the 29 immune-related gene sets, (C–K) the immune status of each patient in the low-risk group and high-risk group.

Response Assessment of Signature Related Genes to Immune Checkpoint Blockade

The accuracy of the risk signature and eight signature related genes in comparison with other nine published immune checkpoint blockade response biomarkers were compared in different immunotherapy cohorts in TIDE web (Figure 12). Athough there is no data on cervical cancer, the signature gave an AUC greater than 0.5 in 10 out of the 21 cohorts (Figure 12A). The AUC of risk signature for PD1 response in head and neck squamous cell carcinoma is 0.94, which is larger than that of all other published markers. The AUC of the risk signature predicted the OS of melanoma treated with CTLA4 is 0.74. In addition, all these signature related genes also showed significant immune checkpoint blockade response in more than 9 cohorts (Figure 12B). The p-value of the signature prediction for the OS of melanoma treated with CTLA4 is less than 0.05. Above results suggesting that this signature and related genes may be robust predictive biomarkers for immune checkpoint blockade response.

FIGURE 12

Figure 12 Comparison of immunotherapy response. (A) The immunotherapy response of risk signature and related genes were compared with the immune checkpoint blockade response status biomarkers by AUC in TIDE web. (B) The eight gene signature is association with overall survival of melamoma patients treated in CTLA4 (by two-sided Wald test).

Prognostic-Related Gene Methylation Status of the Eight-Gene Signature

We analysed DNA methylation data containing 309 CESC tissues and 3 normal tissues. The methylation β value of each methylation site is displayed in the scatter chart based on the arrangement of methylation sites. Through univariate Cox regression analysis, 10 beneficial sites and 9 harmful sites were collected from a total of 131 methylation sites of the eight related genes (Figure 13). PGM1 contained 3 beneficial sites, while LGALS9 had 3 harmful sites. PLIN2 and MSMO1 showed no significant sites for OS. In order to show the prognostic impact of these methylation sites, we performed survival status and clinicopathological information analysis. As depicted in Figure 14, the methylation sites of the PGM1, SPRY4, and LGALS9 genes significantly correlated with the overall survival.

FIGURE 13

Figure 13 The methylation annotation file of the eight prognostic-related genes. Mean gene methylation level and the OS significance of eight risk signature genes (A) CACNB3, (B) SPRY4, (C) PGM1, (D) CD46, (E) LGALS9, (F) DAGLB, (G) PLIN2, and (H) MSMO1. (green: beneficial, red: harmful, turquoise: mean methylation level in normal, purple: mean methylation level in tumor sample).

FIGURE 14

Figure 14 The relationship between the risk signature and DNA methylation status. The visualization of the relationships between TCGA expression, DNA methylation and clinical data by MEXPRESS for (A) CACNB3, (B) CD46, (C) PGM1, (D) SPRY4, (E) LGALS9, (F) PLIN2, (G) MSMO1 and (H) DAGLB.

The Eight-Gene Signature may Lead to Better Treatment Decisions for CESC Patients

Risk prediction have played important roles in making treatment decisions and improve treatment guidelines. Therefore, the accuracy of model predictions is crucial for beneficial treatment decisions. To further evaluate the predictive level of our model, we searched and compared 10 published cervical cancer prognostic signatures with our eight gene signature by NRI and IDI analysis. As results showed in Supplementary Figure 4, our eight gene signature performed better than the other then signature in 1, 3, 5 years OS prediction. but in ten year OS prediction, our signature showed less efficient than a seven gene signature (PLOD2, DSG2, SPP1, CXCL8, MCM5, HLTF, and KLF4) (32) and a six-gene signature (APOC1, GLTP, ISG20, SPP1, SLC24A3 and UPP1) (24). These results reveals that our signature may lead to a better treatment decisions comparison with other published prognostic signatures

Discussions

According to the fact that HLA-G is a normal immune signal that can block the immune response, many recent studies have consistently claimed that HLA-G could be a new immune checkpoint molecule in tumours (8, 10, 11, 13, 34). HLA-G could inhibit almost all stages of the antitumour response by interacting with its inhibitory receptors ILT2 and ILT4 to transmit inhibitory signals. In addition, as a potential tumour marker, HLA-G can predict the diagnosis and prognosis of various types of cancer (7, 9, 11, 35). Therefore, exploring the relationship between HLA-G-positive tumours and their driving DEGs or their impact on the tumour microenvironment is very important for developing effective immunotherapy strategies. In our study, we established and validated, for the first time, a novel prognostic risk model for cervical cancer based on HLA-G-driven DEGs by using a series of bioinformatics methods. The eight HLA-G-driven DEG signatures of our prognostic model suggested great ability in predicting the overall survival of patients with cervical cancer. The prognostic value of the risk score signature was also evaluated in different subtypes of cervical cancer, which has been successfully verified in keratinizing SCC, non-keratinizing SCC, squamous cell neoplasms and non-squamous cell neoplasms. In addition, we established a nomogram model that integrated traditional clinical prognostic factors with the risk score signature and developed a quantitative method for CESC prognosis.

In this study, we obtained 1108 overlapping HLA-G-driven DEGs, which were subjected to GO and KEGG. The GO analysis suggested that HLA-G-driven DEGs were significantly enriched in coated vesicle membrane and cell-substrate adherens junction (in terms of CC group), regulation of lipid biosynthetic and metabolic process (in terms of BP group), protein serine/threonine kinase activity and epidermal growth factor receptor binding (in terms of MF group). These results suggested that HLA-G-driven DEGs were highly correlated with the cancer progression process (36–38). Additionally, KEGG suggested that HLA-G-driven DEGs were significantly enriched pathways in the ErbB, steroid biosynthesis, and MAPK, which are involved in cervical carcinogenesis (38–40). Previous studies have identified that somatic mutations of ErbB2 or ErbB3 receptors are significantly increased in cervical cancer, which provides potential therapeutic targets for patients (40, 41). After multivariate Cox regression analysis, eight genes were identified as key predictors of HLA-G-driven DEGs in the prognostic risk model for predicting the overall survival of patients with cervical cancer, including CD46, LGALS9, PGM1, SPRY4, CACNB3, PLIN2, MSMO1, and DAGLB.

Among eight genes in the prognostic model, CD46 was reported to be a key costimulatory molecule on CD4+ T lymphocytes and could facilitate the differentiation of CD4+ T cells into T regulatory 1 cells (Tr1 cells) secreting interleukin-10 (IL-10) and regulate the tumour immune microenvironment (TIME) (42, 43). The LGALS9 expression level was significantly correlated with HPV16 or HPV18 infection of cervical cells and showed a trend towards improved survival of patients (44, 45). SPRY4 is an inhibitor of the receptor-transduced MAPK signalling pathway, which plays a key role during carcinogenesis (46). CACNB3 may play a regulatory role in transcription factors or/and calcium transport (47). PGM1 was reported to be associated with glycogen metabolism (48). PLIN2 (perilipin 2) is a lipid storage-associated cell protein that is associated with the metabolism of intracellular lipid droplets (49). MSMO1 is believed to function in cholesterol biosynthesis (50). DAGLB is involved in proinflammatory signalling and the regulation of lipid pathways of the immune system (51). In the present study, the HLA-G-driven DEG signature suggested a great predictive ability in the OS of cancer patients.

Similar to HLA-G, CD46 was found to be upregulated in a variety of cancer cells (52). Consistent with previous data (53), our results suggested that high expression of CD46 was related to poor prognosis of cervical cancer. The AUC value at 5 years OS for CD46 was 0.729 (Supplementary Figures S3A and 3I). Notably, CD46 expression in normal tissue is low, and it would be a suitable target for cancer immunotherapy (54). LGALS9 (also known as Galectin-9) is a β-galactoside binding protein that is the ligand for T cell immunoglobulin mucin-3 (TIM-3). Checkpoint molecule TIM-3 exerts immunosuppressive function by interacting with LGALS9, leading to immune escape during carcinogenesis (55, 56). Immunotherapy targeting the PD-1/PD-L1 has been successful in advanced cervical cancer, but the potential clinical value of other checkpoint molecules remains unknown, such as TIM-3/LGALS9 or HLA-G/ILT2. In this study, LGALS9 was a protective factor affecting the prognosis of patients, and its expression level was negatively associated with prognosis (P<0.001) (Supplementary Figures S3B, J). Notably, long noncoding RNA SPRY4-IT1 (SPRY4 intronic transcript 1) is a transcript variant derived from the second intron of the SPRY4 gene, and it could be a good candidate biomarker for discriminating patients with cervical cancer in female population (57, 58). In addition, TIDE analysis showed that these genes also showed significant response in different immunotherapy treatments. Though the essential functions of the other six genes in cervical cancer remain unclear. Thus, further experimental validations of their functions in cervical cancer could strengthen our conclusion.

DNA methylation, the best characterized epigenetic alteration, the results in the exclusion of certain genes and plays an essential role in the progression of cancer (59). In this study, the lower expression of LGALS9 was related to the poor prognosis of cervical cancer patients in our prognostic OS model. A previous study found that this biological effect was mediated through the aberrant epigenetics of LGALS9, which was facilitated by the recruitment of DNMT3A to its promoter region (60). Using 5-aza-2’-deoxycytidine treatment, LGALS9 or SPRY4 expression can be affected by DNA methylation (60, 61). In recent years, DNA methylation changes have been utilized as biomarker for the early detection, diagnosis and prognosis of cervical cancer.

In addition, we explored the correlation between the risk model of cervical cancer and tumour-infiltrating immune cells. We compared the expression profiles of 29 immune-related gene sets between low-risk group and high-risk group, and explored their correlation with this prognostic model. The results suggested that there was less infiltration of CD8+ T cells and HLA in the high-risk group. A previous study showed that the overexpression of HLA-G molecules on the surface of tumour cells affects tumour-specific T cell immunity, and upregulate the expression of ILT2 on tumour-infiltrating CD8+ T cells in the cancer microenvironment (10). Therefore, blocking checkpoint molecule HLA-G could be used as a promising immunotherapeutic strategy for the treatment of cervical cancer.

Conclusions

In summary, we established a novel prognostic risk model using HLA-G-driven DEGs through comprehensive bioinformatics analysis for the first time. These findings provide new insight into the biological markers, DNA methylation, and immunotherapy strategies in cervical cancer.

Data Availability Statement

RNA sequencing data can be accessed via the GEO repository with accesison number GSE208119, access here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE208119.

Author Contributions

H-HX wrote the first draft of the manuscript. All authors listed have made a substantial, direct and intellectual contribution to the work, and approved the submitted version.

Funding

This work was supported by grants from National Natural Science Foundation of China (81901625), Natural Science Foundation of Zhejiang province (LY20H100004, LGF22H160058), Health Bureau of Zhejiang Province (2020KY1037).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.851622/full#supplementary-material

Supplementary Figure 1 | The protein expression of HLA-G in transfected and un-transfected cells by Flow Cytometry using HLA-G specific antibody MEM-G/09-FITC. (A) SiHa-pVITRO2-mcs. (B) SiHa-HLA-G-pVITRO2-mcs. (C) HeLa-pVITRO2-mcs. (D) HeLa-HLA-G-pVITRO2-mcs. MEM-G/09-FITC is an IgG1 antibody detecting native form of HLA-G on the cell surface.

Supplementary Figure 2 | Establishment of prognostic risk model. (A) Tuning parameter (λ) selection in LASSO regression analysis. (B) LASSO coefficient profiles of 1108 HLA-G-driven DEGs.

Supplementary Figure 3 | Kaplan-Meier and ROC of OS for genes in the signature. (A, I) CD46, (B, J) LGALS9, (C, K) PGM1, (D, L) SPRY4, (E, M) CACNB3, (F, N) PLIN2, (G, O) MSMO1 and (H, P) DAGLB.

Supplementary Figure 4 | Clinical benefits evaluation by NRI and NRI analysis. Graphically display IDI, continuous NRI, and median improvement of eight gene signature compared with other 10 published CESC prediction signatures by IDI.INF.GRAPH in “survIDINRI” package.

References

1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2018) 68(6):394–424. doi: 10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

2. WHO. Cervical Cancer . World Health Organization. Available at: https://www.who.int/news-room/fact-sheets/detail/human-papillomavirus-(hpv)-and-cervical-cancer (Accessed September 21, 2020).

Google Scholar

3. Cohen PA, Jhingran A, Oaknin A, Denny L. Cervical Cancer. Lancet (2019) 393(10167):169–82. doi: 10.1016/S0140-6736(18)32470-X

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Yang Y. Cancer Immunotherapy: Harnessing the Immune System to Battle Cancer. J Clin Invest (2015) 125(9):3335–7. doi: 10.1172/JCI83871

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Xu HH, Zhang X, Zheng HH, Han QY, Lin AF, Yan WH. Association of HLA-G 3' UTR Polymorphism and Expression With the Progression of Cervical Lesions in Human Papillomavirus 18 Infections. Infect Agent Cancer (2018) 13:42. doi: 10.1186/s13027-018-0217-2

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Xu HH, Yan WH, Lin A. The Role of HLA-G in Human Papillomavirus Infections and Cervical Carcinogenesis. Front Immunol (2020) 11:1349. doi: 10.3389/fimmu.2020.01349

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Liu L, Wang L, Zhao L, He C, Wang G. The Role of HLA-G in Tumor Escape: Manipulating the Phenotype and Function of Immune Cells. Front Oncol (2020) 10:597468. doi: 10.3389/fonc.2020.597468

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Carosella ED, Ploussard G, LeMaoult J, Desgrandchamps F. A Systematic Review of Immunotherapy in Urologic Cancer: Evolving Roles for Targeting of CTLA-4, PD-1/PD-L1, and HLA-G. Eur Urol (2015) 68(2):267–79. doi: 10.1016/j.eururo.2015.02.032

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Własiuk P, Putowski M, Giannopoulos K. PD1/PD1L Pathway, HLA-G and T Regulatory Cells as New Markers of Immunosuppression in Cancers. Postepy Hig Med Dosw (Online) (2016) 70(0):1044–58. doi: 10.5604/17322693.1220994

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Dumont C, Jacquier A, Verine J, Noel F, Goujon A, Wu CL, et al. CD8+PD-1-ILT2+ T Cells Are an Intratumoral Cytotoxic Population Selectively Inhibited by the Immune-Checkpoint HLA-G. Cancer Immunol Res (2019) 7(10):1619–32. doi: 10.1158/2326-6066.CIR-18-0764

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Rouas-Freiss N, LeMaoult J, Verine J, Tronik-Le Roux D, Culine S, Hennequin C, et al. Intratumor Heterogeneity of Immune Checkpoints in Primary Renal Cell Cancer: Focus on HLA-G/Ilt2/Ilt4. Oncoimmunology (2017) 6(9):e1342023. doi: 10.1080/2162402X.2017.1342023

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Kuroki K, Matsubara H, Kanda R, Miyashita N, Shiroishi M, Fukunaga Y, et al. Structural and Functional Basis for LILRB Immune Checkpoint Receptor Recognition of HLA-G Isoforms. J Immunol (2019) 203(12):3386–94. doi: 10.4049/jimmunol.1900562

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Carosella ED, Rouas-Freiss N, Tronik-Le Roux D, Moreau P, LeMaoult J. HLA-G: An Immune Checkpoint Molecule. Adv Immunol (2015) 127:33–144. doi: 10.1016/bs.ai.2015.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Li XJ, Zhang X, Lin A, Ruan YY, Yan WH. Human Leukocyte Antigen-G (HLA-G) Expression in Cervical Cancer Lesions is Associated With Disease Progression. Hum Immunol (2012) 73(9):946–9. doi: 10.1016/j.humimm.2012.07.041

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Zheng N, Wang CX, Zhang X, Du LT, Zhang J, Kan SF, et al. Up-Regulation of HLA-G Expression in Cervical Premalignant and Malignant Lesions. Tissue Antigens (2011) 77(3):218–24. doi: 10.1111/j.1399-0039.2010.01607.x

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Yoon BS, Kim YT, Kim JW, Kim SH, Kim JH, Kim SW. Expression of Human Leukocyte Antigen-G and its Correlation With Interleukin-10 Expression in Cervical Carcinoma. Int J Gynaecol Obstet (2007) 98(1):48–53. doi: 10.1016/j.ijgo.2007.03.041

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Rodríguez JA, Galeano L, Palacios DM, Gómez C, Serrano ML, Bravo MM, et al. Altered HLA Class I and HLA-G Expression is Associated With IL-10 Expression in Patients With Cervical Cancer. Pathobiology (2012) 79(2):72–83. doi: 10.1159/000334089

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Lin A, Yan WH, Xu HH, Gan MF, Cai JF, Zhu M, et al. HLA-G Expression in Human Ovarian Carcinoma Counteracts NK Cell Function. Ann Oncol (2007) 18(11):1804–9. doi: 10.1093/annonc/mdm356

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Huang J, Liang X, Xuan Y, Geng C, Li Y, Lu H, et al. A Reference Human Genome Dataset of the BGISEQ-500 Sequencer. Gigascience (2017) 6(5):1–9. doi: 10.1093/gigascience/gix024

CrossRef Full Text | Google Scholar

20. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Pencina MJ, D'Agostino RB Sr, D'Agostino RB Jr, Vasan RS. Evaluating the Added Predictive Ability of a New Marker: From Area Under the ROC Curve to Reclassification and Beyond. Stat Med (2008) 27(2):157–72 discussion 207-12. doi: 10.1002/sim.2929

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Pencina MJ, D'Agostino RB Sr, Steyerberg EW. Extensions of Net Reclassification Improvement Calculations to Measure Usefulness of New Biomarkers. Stat Med (2011) 30(1):11–21. doi: 10.1002/sim.4085

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Ju M, Qi A, Bi J, Zhao L, Jiang L, Zhang Q, et al. A five-mRNA Signature Associated With Post-Translational Modifications can Better Predict Recurrence and Survival in Cervical Cancer. J Cell Mol Med (2020) 24(11):6283–97. doi: 10.1111/jcmm.15270

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Qu X, Shi Z, Guo J, Guo C, Qiu J, Hua K. Identification of a Novel Six-Gene Signature With Potential Prognostic and Therapeutic Value in Cervical Cancer. Cancer Med (2021) 10(19):6881–96. doi: 10.1002/cam4.4054

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Yu S, Li X, Zhang J, Wu S. Development of a Novel Immune Infiltration-Based Gene Signature to Predict Prognosis and Immunotherapy Response of Patients With Cervical Cancer. Front Immunol (2021) 12:709493. doi: 10.3389/fimmu.2021.709493

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Pan J, Xu L, Pan H. Development and Validation of an M6a RNA Methylation Regulator-Based Signature for Prognostic Prediction in Cervical Squamous Cell Carcinoma. Front Oncol (2020) 10:1444. doi: 10.3389/fonc.2020.01444

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Qin R, Cao L, Ye C, Wang J, Sun Z. A Novel Prognostic Prediction Model Based on Seven Immune-Related RNAs for Predicting Overall Survival of Patients in Early Cervical Squamous Cell Carcinoma. BMC Med Genomics (2021) 14(1):49. doi: 10.1186/s12920-021-00885-3

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Zuo Z, Xiong J, Zeng C, Jiang Y, Xiong K, Tao H, et al. Exploration of a Robust and Prognostic Immune Related Gene Signature for Cervical Squamous Cell Carcinoma. Front Mol Biosci (2021) 8:625470. doi: 10.3389/fmolb.2021.625470

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Yang Y, Li Y, Qi R, Zhang L. Constructe a Novel 5 Hypoxia Genes Signature for Cervical Cancer. Cancer Cell Int (2021) 21(1):345. doi: 10.1186/s12935-021-02050-3

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Shang C, Huang J, Guo H. Identification of an Metabolic Related Risk Signature Predicts Prognosis in Cervical Cancer and Correlates With Immune Infiltration. Front Cell Dev Biol (2021) 9:677831. doi: 10.3389/fcell.2021.677831

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Feng Y, Wang Z, Yang N, Liu S, Yan J, Song J, et al. Identification of Biomarkers for Cervical Cancer Radiotherapy Resistance Based on RNA Sequencing Data. Front Cell Dev Biol (2021) 9:724172. doi: 10.3389/fcell.2021.724172

PubMed Abstract | CrossRef Full Text | Google Scholar

32. An R, Meng S, Qian H. Identification of Key Pathways and Establishment of a Seven-Gene Prognostic Signature in Cervical Cancer. J Oncol (2022) 2022:4748796. doi: 10.1155/2022/4748796

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Becht E, McInnes L, Healy J, Dutertre CA, Kwok IWH, Ng LG, et al. Dimensionality Reduction for Visualizing Single-Cell Data Using UMAP. Nat Biotechnol (2019) 37:38–44. doi: 10.1038/nbt.4314

CrossRef Full Text | Google Scholar

34. Krijgsman D, Roelands J, Hendrickx W, Bedognetti D, Kuppen PJK. HLA-G: A New Immune Checkpoint in Cancer? Int J Mol Sci (2020) 21(12):4528. doi: 10.3390/ijms21124528

CrossRef Full Text | Google Scholar

35. Hiraoka N, Ino Y, Hori S, Yamazaki-Itoh R, Naito C, Shimasaki M, et al. Expression of Classical Human Leukocyte Antigen Class I Antigens, HLA-E and HLA-G, is Adversely Prognostic in Pancreatic Cancer Patients. Cancer Sci (2020) 111(8):3057–70. doi: 10.1111/cas.14514

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Urabe F, Kosaka N, Ito K, Kimura T, Egawa S, Ochiya T. Extracellular Vesicles as Biomarkers and Therapeutic Targets for Cancer. Am J Physiol Cell Physiol (2020) 318(1):C29–39. doi: 10.1152/ajpcell.00280.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Ma Y, Zhang S, Jin Z, Shi M. Lipid-Mediated Regulation of the Cancer-Immune Crosstalk. Pharmacol Res (2020) 161:105131. doi: 10.1016/j.phrs.2020.105131

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Guo YJ, Pan WW, Liu SB, Shen ZF, Xu Y, Hu LL. ERK/MAPK Signalling Pathway and Tumorigenesis. Exp Ther Med (2020) 19(3):1997–2007. doi: 10.3892/etm.2020.8454

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Manzo-Merino J, Contreras-Paredes A, Vázquez-Ulloa E, Rocha-Zavaleta L, Fuentes-Gonzalez AM, Lizano M. The Role of Signaling Pathways in Cervical Cancer and Molecular Therapeutic Targets. Arch Med Res (2014) 45(7):525–39. doi: 10.1016/j.arcmed.2014.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zammataro L, Lopez S, Bellone S, Pettinella F, Bonazzoli E, Perrone E, et al. Whole-Exome Sequencing of Cervical Carcinomas Identifies Activating ERBB2 and PIK3CA Mutations as Targets for Combination Therapy. Proc Natl Acad Sci U.S.A. (2019) 116(45):22730–6. doi: 10.1073/pnas.1911385116

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Cancer Genome Atlas Research Network, Albert Einstein College of Medicine, Analytical Biological Services, Barretos Cancer Hospital, Baylor College of Medicine, Beckman Research Institute of City of Hope, et al. Integrated Genomic and Molecular Characterization of Cervical Cancer. Nature (2017) 543(7645):378–84. doi: 10.1038/nature21386

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Ni Choileain S, Hay J, Thomas J, Williams A, Vermeren MM, Benezech C, et al. TCR-Stimulated Changes in Cell Surface CD46 Expression Generate Type 1 Regulatory T Cells. Sci Signal (2017) 10(502):eaah6163. doi: 10.1126/scisignal.aah6163

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Kemper C, Verbsky JW, Price JD. Atkinson JP. T-Cell Stimulation and Regulation: With Complements From CD46. Immunol Res (2005) 32(1-3):31–43. doi: 10.1385/IR:32:1-3:031

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Punt S, Thijssen VL, Vrolijk J, de Kroon CD, Gorter A, Jordanova ES. Galectin-1, -3 and -9 Expression and Clinical Significance in Squamous Cervical Cancer. PloS One (2015) 10(6):e0129119. doi: 10.1371/journal.pone.0129119

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wang L, Zhao Y, Wang Y, Wu X. The Role of Galectins in Cervical Cancer Biology and Progression. BioMed Res Int (2018) 2018:2175927. doi: 10.1155/2018/2175927

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Sun M, Huang F, Yu D, Zhang Y, Xu H, Zhang L, et al. Autoregulatory Loop Between TGF-β1/miR-411-5p/SPRY4 and MAPK Pathway in Rhabdomyosarcoma Modulates Proliferation and Differentiation. Cell Death Dis (2015) 6(8):e1859. doi: 10.1038/cddis.2015.225

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Cui C, Merritt R, Fu L, Pan Z. Targeting Calcium Signaling in Cancer Therapy. Acta Pharm Sin B (2017) 7(1):3–17. doi: 10.1016/j.apsb.2016.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Jin GZ, Zhang Y, Cong WM, Wu X, Wang X, Wu S, et al. Phosphoglucomutase 1 Inhibits Hepatocellular Carcinoma Progression by Regulating Glucose Trafficking. PloS Biol (2018) 16(10):e2006483. doi: 10.1371/journal.pbio.2006483

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Conte M, Franceschi C, Sandri M, Salvioli S. Perilipin 2 and Age-Related Metabolic Diseases: A New Perspective. Trends Endocrinol Metab (2016) 27(12):893–903. doi: 10.1016/j.tem.2016.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Gabitova L, Restifo D, Gorin A, Manocha K, Handorf E, Yang DH, et al. Endogenous Sterol Metabolites Regulate Growth of EGFR/KRAS-Dependent Tumors via LXR. Cell Rep (2015) 12(11):1927–38. doi: 10.1016/j.celrep.2015.08.023

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Hsu KL, Tsuboi K, Adibekian A, Pugh H, Masuda K, Cravatt BF. Daglβ Inhibition Perturbs a Lipid Network Involved in Macrophage Inflammatory Responses. Nat Chem Biol (2012) 8(12):999–1007. doi: 10.1038/nchembio.1105

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Geller A, Yan J. The Role of Membrane Bound Complement Regulatory Proteins in Tumor Development and Cancer Immunotherapy. Front Immunol (2019) 10:1074. doi: 10.3389/fimmu.2019.01074

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Chen H, Deng Q, Wang W, Tao H, Gao Y. Identification of an Autophagy-Related Gene Signature for Survival Prediction in Patients With Cervical Cancer. J Ovarian Res (2020) 13(1):131. doi: 10.1186/s13048-020-00730-8

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Elvington M, Liszewski MK, Atkinson JP. CD46 and Oncologic Interactions: Friendly Fire Against Cancer. Antibod (Basel) (2020) 9(4):59. doi: 10.3390/antib9040059

CrossRef Full Text | Google Scholar

55. Curley J, Conaway MR, Chinn Z, Duska L, Stoler M, Mills AM. Looking Past PD-L1: Expression of Immune Checkpoint TIM-3 and its Ligand Galectin-9 in Cervical and Vulvar Squamous Neoplasia. Mod Pathol (2020) 33(6):1182–92. doi: 10.1038/s41379-019-0433-3

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Zhang L, Tian S, Pei M, Zhao M, Wang L, Jiang Y, et al. Crosstalk Between Histone Modification and DNA Methylation Orchestrates the Epigenetic Regulation of the Costimulatory Factors, Tim-3 and Galectin-9, in Cervical Cancer. Oncol Rep (2019) 42(6):2655–69. doi: 10.3892/or.2019.7388

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Dong J, Su M, Chang W, Zhang K, Wu S, Xu T. Long non-Coding RNAs on the Stage of Cervical Cancer (Review). Oncol Rep (2017) 38(4):1923–31. doi: 10.3892/or.2017.5905

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Cao Y, Liu Y, Lu X, Wang Y, Qiao H, Liu M. Upregulation of Long Noncoding RNA SPRY4-IT1 Correlates With Tumor Progression and Poor Prognosis in Cervical Cancer. FEBS Open Bio (2016) 6(9):954–60. doi: 10.1002/2211-5463.12102

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Zhu H, Zhu H, Tian M, Wang D, He J, Xu T, et al. and Hydroxymethylation in Cervical Cancer: Diagnosis, Prognosis and Treatment. Front Genet (2020) 11:347. doi: 10.3389/fgene.2020.00347

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Zhang L, Tian S, Zhao M, Yang T, Quan S, Yang Q, et al. SUV39H1-DNMT3A-Mediated Epigenetic Regulation of Tim-3 and Galectin-9 in the Cervical Cancer. Cancer Cell Int (2020) 20:325. doi: 10.1186/s12935-020-01380-y

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Zhou X, Xie S, Yuan C, Jiang L, Huang X, Li L, et al. Lower Expression of SPRY4 Predicts a Poor Prognosis and Regulates Cell Proliferation in Colorectal Cancer. Cell Physiol Biochem (2016) 40(6):1433–42. doi: 10.1159/000453195

PubMed Abstract | CrossRef Full Text | Google Scholar