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

Abstract

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.

Keywords: HLA-G; cervical cancer; immune checkpoint; mRNA signature; prediction; prognosis.

Figure 1

Flow diagram of this study.

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

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.

Figure 4

PPI network based on the data from STRING database.

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).

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

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.

Figure 8

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

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.

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.

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.

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).

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

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.