Mining of immunological and prognostic-related biomarker for cervical cancer based on immune cell signatures

Abstract

Background: Immunotherapy has changed the therapeutic landscape of cervical cancer (CC), but has durable anti-tumor activity only in a subset of patients. This study aims to comprehensively analyze the tumor immune microenvironment (TIME) of CC and to mine biomarkers related to immunotherapy and prognosis.

Methods: The Cancer Genome Atlas (TCGA) data was utilized to identify heterogeneous immune subtypes based on survival-related immune cell signatures (ICSs). ICSs prognostic model was constructed by Cox regression analyses, and immunohistochemistry was conducted to verify the gene with the largest weight coefficient in the model. Meanwhile, the tumor immune infiltration landscape was comprehensively characterized by ESTIMATE, CIBERSORT and MCPcounter algorithms. In addition, we also analyzed the differences in immunotherapy-related biomarkers between high and low-risk groups. IMvigor210 and two gynecologic tumor cohorts were used to validate the reliability and scalability of the Risk score.

Results: A total of 291 TCGA-CC samples were divided into two ICSs clusters with significant differences in immune infiltration landscape and prognosis. ICSs prognostic model was constructed based on eight immune-related genes (IRGs), which showed higher overall survival (OS) rate in the low-risk group (P< 0.001). In the total population, time-dependent receiver operating characteristic (ROC) curves displayed area under the curve (AUC) of 0.870, 0.785 and 0.774 at 1-, 3- and 5-years. Immunohistochemical results showed that the expression of the oncogene (FKBP10) was negatively correlated with the degree of differentiation and positively correlated with tumor stage, while the expression of tumor suppressor genes (S1PR4) was the opposite. In addition, the low-risk group had more favorable immune activation phenotype and higher enrichment of immunotherapy-related biomarkers. The Imvigor210 and two gynecologic tumor cohorts validated a better survival advantage and immune efficacy in the low-risk group.

Conclusion: This study comprehensively assessed the TIME of CC and constructed an ICSs prognostic model, which provides an effective tool for predicting patient's prognosis and accurate immunotherapy.

Keywords: biomarkers; cervical cancer; immunohistochemistry; immunotherapy; tumor immune microenvironment.

Figure 1

Flow chart of the study.

Figure 2

Determination of immune cell signatures (ICSs) subtypes. (A) Consensus matrix of the TCGA-CC cohort with appropriate k values (k = 2). (B) PCA validation of clustering results. (C) Kaplan-Meier curves of OS for CC patients in both ICSs clusters (P = 0.016). (D) Violin plots of 22 tumor-infiltrating immune cell types of two ICSs clusters by CIBERSORT algorithm. (E) Differential analysis of 27 immune marker genes in two ICSs clusters (ns, not significant; *P< 0.050; **P< 0.010; ***P< 0.001). (F) Differences between two ICSs clusters in the degree of enrichment of indicator signals for specific immune functions (ns, not significant; *P< 0.050; **P< 0.010; ***P< 0.001). (G) Biological processes of two ICSs clusters using GSVA analysis. Heatmap colors indicate ICSs infiltrate levels, with red indicating high infiltrate levels and blue indicating low infiltrate levels. (ns, not significant; *P< 0.050; **P< 0.010; ***P< 0.001).

Figure 3

Construction and validation of ICSs prognostic model in the TCGA cohort. (A) Heatmap depicted the expression levels of survival-related DEGs in different ICSs clusters and the distribution of clinical traits of patients. The rows represent survival-related DEGs and the columns represent samples. (B, C) Determination of the number of survival-related DEGs into the multivariate Cox regression model by LASSO analyses. (D–F) Kaplan-Meier curves of OS for the high and low-risk groups in the total population, training set, and validation set (P< 0.001). (G–I) Time-dependent ROC curve in the total population, training set, and validation set.

Figure 4

Independent prognostic analysis of Risk score. (A) Nomogram for predicting the probability of patient mortality at 1-, 3-, or 5- year OS based on two independent prognosis factors (***P< 0.001). (B) Calibration curves of the nomogram for predicting the probability of OS at 1-, 3-, or 5- year. (C) Decision curve analyses (DCAs) of the nomograms for 1-, 3-, or 5- year risk. (D) Time-dependent ROC curve of two independent prognosis factors. (E) Heatmap and table showing the distribution of Stage I-IV between high and low-risk groups (P = 0.036).

Figure 5

Immunohistochemical results of FKBP10 and S1PR4. (A) Immunohistochemical images of FKBP10 in high and low differentiation groups. (B) Immunohistochemical images of S1PR4 in high and low differentiation groups. (C) Immunohistochemical images of FKBP10 in four different tumor stages (Stage I-IV). (D) Immunohistochemical images of S1PR4 in four different tumor stages (Stage I-IV).

Figure 6

Immune infiltration landscape in the high and low-risk groups. (A) Alluvial diagram of the distribution with patients in different ICSs clusters, risk groups, and survival outcomes. (B) Heatmap depicted the infiltration of survival-related ICSs in the high and low-risk groups and the distribution of clinical traits of patients. The rows represent survival-related ICSs and the columns represent samples. (C) Differential analysis of immune score, stromal score and estimate score in the high and low-risk groups (ns, not significant; *P< 0.050; **P< 0.010; ***P< 0.001). (D) Violin plots of 22 tumor-infiltrating immune cell types of high and low-risk groups by CIBERSORT algorithm. (E) Violin plots of 10 tumor-infiltrating immune cell types of high and low-risk groups by MCPcounter algorithm (ns, not significant; *P< 0.050; **P< 0.010; ***P< 0.001). (F) Heatmap and table showing the distribution of pan- SCC immune subtypes (IS1, IS2, IS3, IS4, IS5, and IS6) between high and low-risk groups (P = 0.001). (G) Heatmap and table showing the distribution of immune subtypes (C1, C2, C3, C4, C5, and C6) between high and low-risk groups (P = 0.001). (H) Biological processes of high and low-risk groups using GSVA analysis. Heatmap colors indicate ICSs infiltrate levels, with red indicating high infiltrate levels and blue indicating low infiltrate levels. (ns, not significant; *P< 0.050; **P< 0.010; ***P< 0.001).

Figure 7

Comparison of immunotherapy predictive biomarker. (A) Differential analysis of 27 immune marker genes in the high and low-risk groups (ns, not significant; *P< 0.050; **P< 0.010; ***P< 0.001). (B) Differences between high and low-risk groups in the degree of enrichment of indicator signals for specific immune functions (ns, not significant; *P< 0.050; **P< 0.010; ***P< 0.001). (C, D) The waterfall diagram of the top 20 driver genes between the high (C) and low risk-score (D) of CC patients. (E, F) Differential analysis of MHC molecules (E) and chemokines (F) in the high and low-risk groups (ns, not significant; *P< 0.050; **P< 0.010; ***P< 0.001). (G) CYT difference in the high and low-risk groups (ns, not significant; *P< 0.050; **P< 0.010; ***P< 0.001). (H) STING difference in the high and low-risk groups (ns, not significant; *P< 0.050; **P< 0.010; ***P< 0.001). (I) Scatterplots depicting the correlation between Risk score and CYT (R = -0.380, P<0.001). (J) Scatterplots depicting the correlation between Risk score and STING (R = -0.250, P<0.001).

Figure 8

The role of Risk score in immunotherapeutic response prediction. (A–D) The distribution plot of ips_ctla4_pos_pd1_pos (A), ips_ctla4_pos_pd1_neg (B), ips_ctla4_neg_pd1_pos (C), and ips_ctla4_neg_pd1_neg (D) scores in TCGA-CC cohort. (E, F) Intrinsic connection of Risk score and MHC, EC, SC, IPS score in TCGA-CC and IMvigor210 cohorts, with red indicating positive correlations and blue indicating negative correlations. The asterisks represented the statistical P value (*P< 0.050). (G) Distribution of clinical response rates for anti-PD-L1 immunotherapy in the high and low-Risk score groups in the IMvigor210 cohort (P = 0.020). (H) Risk score in groups with different anti-PD-L1 clinical response status (P = 0.160). (I) Kaplan-Meier curves for OS in the IMvigor210 cohort for the high and low-risk groups (P = 0.001). (J) Kaplan-Meier curves for OS in the TCGA-OC cohort for the high and low-risk groups (P = 0.002). (K) Kaplan-Meier curves for OS in the TCGA-EC cohort for the high and low-risk groups (P = 0.030).