A Novel Immune-Related Signature to Predict Prognosis and Immune Infiltration of Cervical Cancer

Background

Cervical cancer is the fourth most common cancer and the leading cause of cancer-related death among women globally [1]. Ascribed to improvement of socio-economic levels and diminishing prevalence of sexually transmitted disease, the incidence rate of cervical cancer has decreased in many countries for the past few decades [2,3]. However, the mortality rate remains high and the median survival time of advanced cervical cancer patients is only 16.8 months [4,5]. According to the estimates produced by International Agency for Research on Cancer, there were 342 000 deaths caused by cervical cancer in 2020 worldwide [1]. Consequently, novel and effective therapeutic strategies are urgently needed.

In recent years, immunotherapy has been applied in cervical cancer. It can recruit the host immune system and kill cancer cells directly, thereby attenuating the recurrence and metastasis rates of cancer [6]. Nonetheless, there is significant heterogeneity in the therapeutic response and prognoses of cervical cancer patients with immunotherapy [7]. Only a subset of patients can derive benefit from immunotherapy. Thus, it is of great clinical significance to screen the sensitive patients. Extensive publications have reported that immune cell infiltration is closely related to the clinical therapeutic effect in various tumors [8,9]. Clinical efficacy of immunotherapies depends on a high density of preexistent tumor infiltrating immune cells [10]. Therefore, it is imperative for developing the suitable biomarkers to predict the immune cells infiltration and prognosis of cervical cancer patients.

In this research, we extracted and established a novel immune-related signature based on the analysis of cervical cancer genome and immunome. This signature can predict the immune cells infiltration and prognosis of cervical cancer simultaneously. It may contribute to better clinical personalized treatment and improve cervical cancer patients’ prognoses.

Material and Methods

DATA ACQUISITION:

The transcriptomic data of 306 cervical cancer samples from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) was downloaded by using UCSC Xena (https://xenabrowser.net/datapages/). After integrating clinical survival information, 293 cervical cancer patients were enrolled. In addition, the data of 10 normal cervical samples from Genotype-Tissue Expression (GTEx) databases and 3 normal paracancerous samples from TCGA was evaluated and integrated. A microarray dataset GSE44001 based on the platform GPL14951 (Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip) was obtained from GEO database (http://www.ncbi.nlm.nih.gov/geo). The expression values were normalized using robust multi-array average (RMA) method and transformed into the log2 scales. The TCGA cohort was randomly categorized into a training cohort and a validation cohort. The entire TCGA cohort and GSE44001 cohort were regarded as the internal testing set and external testing set, respectively.

CONSTRUCTION OF A NOVEL IMMUNE-RELATED SIGNATURE:

The R package “DESeq2” was used to screen the differentially expressed genes (DEGs) between cervical cancer tissues and normal cervical tissues [11]. The threshold criteria were set as |log2 (Fold Change) |>2 and adjusted P value <0.05. Volcano plots were created by using R package “ggplot2” (version 3.3.3). The list of immune-related genes (IRGs) was downloaded from Immunology Database and Analysis Portal (https://www.immport.org/home) [12]. A Venn diagram was generated to obtain the overlapping genes between differentially expressed genes (DEGs) and immune-related genes.

Univariate Cox proportional hazards regression analyses were performed to investigate the correlation between the expression of those overlapping genes and the overall survival (OS) of cervical cancer patients [13]. Genes with P value <0.05 were selected as candidate genes. Furthermore, Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis was conducted with 10-fold cross-validation by using R package “glmnet” and “survival”. The genes with non-zero coefficients were selected as the hub immune-related genes to construct the novel signature. The risk score of the novel signature for each cervical cancer patient was calculated according to the expression level of each gene and corresponding LASSO regression coefficient. The cervical cancer patients in the training cohort, validation cohort, internal testing cohort, and external testing cohort were categorized into low-risk and high-risk groups based on the median value of risk scores.

VALIDATION OF THE NOVEL SIGNATURE:

To validate the predictive power of this novel immune-related signature, a Kaplan-Meier (K-M) survival curve was applied to estimate and compare the overall survival rates of low-risk and high-risk groups. The time-dependent receiver operating characteristic (ROC) curve illustrated the predictive sensitivity and specificity of this signature at 1, 3, and 5 years. All analyses described above were performed in 4 datasets: the training cohort, validation cohort, internal testing cohort, and external testing cohort. Univariate and multivariate Cox proportional hazards regression analyses were used to determine whether the novel signature was an independent prognostic factors for cervical cancer.

CONSTRUCTION AND VALIDATION OF A NOMOGRAM:

A nomogram was constructed with both the risk score of the novel signature and traditional clinical characteristics, including age, body mass index (BMI), clinical stage, M stage, N stage, and histological type. The calibration plots visually revealed the nomogram-predicted survival probability at 1, 3, and 5 years, which were achieved using R package “rms” (version 6.2-0) and “survival” (version 3.2-10) [13]. Decision curve analysis (DCA) was used to evaluate the clinical net benefit of different prediction models.

FUNCTIONAL ENRICHMENT ANALYSIS:

Functional enrichment analyses were performed to examine the biological functions of hub immune-related genes. Gene Ontology (GO) terms, including biological processes (BP), cellular components (CC), and molecular function (MF), were enriched and visualized via R package “GOplot”, “clusterProfiler”, and “ggplot2” [14,15].

ESTIMATION OF IMMUNE CELLS INFILTRATION:

The proportion of immune cells infiltration in each cervical cancer sample of TCGA datasets was estimated though a computational analytical tool CIBERSORT (http://cibersort.stanford.edu/) [16]. The relationship between the risk scores of novel signature and immune cells infiltration was investigated using single-sample GSEA (ssGSEA) method in R package “GSVA” (version 1.34.0) [17,18]. The fractions of immune cells between high-risk and low-risk groups were also compared.

VALIDATION OF THE HUB IMMUNE-RELATED GENES EXPRESSION:

To validate the expression of the 10 hub immune-related genes, we detected the mRNA and protein expression levels of these genes in cervical cancer tissues and normal tissues. The differential mRNA expression of 10 hub immune-related genes between cervical cancer tissues and normal tissues was investigated by RNA sequencing. The protein expression of these immune-related genes in cervical cancer tissues and normal tissues was detected by immunohistochemistry. The staining results were acquired from the Human Protein Atlas database (https://www.proteinatlas.org/).

STATISTICAL ANALYSIS:

All statistical analyses were processed with R software. Univariate and multivariate Cox regression analyses were conducted to determine the hazard ratios (HR) and 95% confidence intervals (CI) for events. The Mann-Whitney U test or t test was used to compare the gene expression between cervical normal and cancer tissues and the immune cells infiltrating fraction between low-risk and high-risk groups. Correlations between the risk scores of novel signature and immune cells infiltration were calculated using Spearman’s rank correlation analyses. A two-tailed P<0.05 was regarded as statistically significant (* P<0.05, ** P<0.01, *** P<0.001).

Results

IDENTIFICATION OF HUB IMMUNE-RELATED GENES:

Based on the differential expression analysis, there were 4338 differentially expressed genes presenting 2-fold or greater change between cervical cancer tissues and normal tissues. Therein, 2801 genes displayed up-regulated expression and 1537 genes displayed down-regulated expression (Figure 1A). Among 4338 differentially expressed genes, there were 298 genes having immunological-related functions (Figure 1B). Further univariate Cox proportional hazards regression analysis evaluated the prognostic values of these 298 genes. As shown in Supplementary Table 1, there were 38 genes significantly related with the overall survival of cervical cancer patients (P<0.05). Subsequently, LASSO regression modeling further minimized overfitting and selected 10 genes as the hub immune-related genes: TFRC, ZAP70, SPP1, CHIT1, CAMP, CSF2, TUBB3, LEPR, DLL4, and DES (Figure 1C). The forest plot showed that these hub immune-related genes have promise for use as prognostic biomarkers in cervical cancer (all P<0.05) (Figure 1D).

CONSTRUCTION OF A NOVEL IMMUNE-RELATED SIGNATURE:

The novel immune-related signature was constructed based on those 10 hub immune-related genes. The risk score of this novel immune-related signature was calculated using the expression level and corresponding correlation coefficient of each gene. Risk score=0.03749*TFRC-0.14105*ZAP70+0.05229*SPP1-0.06349*CHIT1-0.05564*CAMP+0.00302*CSF2+0.05782*TUBB3+0.06949*LEPR +0.13824*DLL4-0.03868*DES. The median value of risk scores for novel signature was identified as the cutoff point to categorize patients into a low-risk group and a high-risk group. In the training set, the cutoff value was 0.508 (Figure 2A). According to the cutoff, the patients in the high-risk group had significantly worse overall survival than those in the low-risk group (P<0.001). We evaluated the sensitivity and specificity of this novel immune-related signature for predicting the prognosis of patients with cervical cancer by using an ROC curve. The AUC of this signature reached 0.776 at 1 year, 0.747 at 3 years, and 0.778 at 5 years, which indicated a relatively high reliability of this novel signature. Taken together, this novel signature could predict the overall survival of cervical cancer patients.

VALIDATION OF THE NOVEL IMMUNE-RELATED SIGNATURE:

To validate the predictive accuracy and stability of this novel signature, the similar workflow was employed in another 3 datasets (Figure 2B–2D). The distribution of risk scores, survival statuses, and hub immune-related genes expression for these 3 datasets was exhibited. All cervical cancer patients were allocated into a high-risk or a low-risk group based on the cutoff point of the training cohort. For the validation cohort, Kaplan-Meier survival analysis revealed that patients in the high-risk group had worse prognosis than patients in the low-risk group (P<0.001). AUC values of 0.854, 0.857, and 0.853 were obtained for the 1-, 3-, and 5-year survival rates, respectively. Similar results were obtained in the internal and external testing sets. Significantly worse overall survival rates were observed in the cervical cancer patients with high-risk scores than those with low-risk scores (P<0.001 for the entire TCGA dataset and P=0.023 for GSE44001). The AUC values in 1, 3, and 5 years were all above 0.8 for the entire TCGA dataset and above 0.6 for GSE44001. Moreover, ROC curves demonstrated that the novel immune-related signature had better predictive ability than other independent genes, with the area under curve (AUC) of 0.739 (Supplementary Figure 1). Altogether, the novel signature was a highly shown to be a reliable biomarker for predicting the prognosis of cervical cancer patients.

PROGNOSTIC VALUE OF THE NOVEL IMMUNE-RELATED GENE SIGNATURE:

To determine whether the novel immune-related gene signature was an independent prognostic factor for cervical cancer patients, univariate and multivariate Cox regression analyses were performed in the entire TCGA cohort. As shown in Table 1, the immune-related gene signature, body mass index (BMI), clinical stage, M stage, and N stage were significantly associated with the overall survival of cervical cancer patients (all P<0.05). Further multivariate analysis revealed that the novel immune-related signature remained an independent prognostic indicator (HR=3.291, 95% CI=1.609–6.729, P=0.001). Then, a clinically quantitative prognostic model nomogram was developed based on the risk scores of this novel signature and some clinicopathological factors, including age, BMI, clinical stage, M stage, N stage, and histological type (Figure 3). The concordance index of this nomogram was 0.813 (0.782–0.843). Moreover, the calibration curves of this nomogram showed the excellent concordance between the predictive and actual 1-, 3-, and 5-year overall survival rates in patients with cervical cancer (Figure 4A–4C). The decision curve analysis (DCA) also indicated that the nomogram model, constructed by both novel signature and clinical characteristics, achieved the highest clinical net benefit across the entire range of reasonable threshold probabilities (Figure 4D). In addition, the risk score of the novel signature provided greater clinical net benefit than traditional clinical characteristics.

GENE ENRICHMENT ANALYSIS:

To further explore the biological function of the 10 hub immune-related genes, biological processes (BP), cellular components (CC), and molecular function (MF) enrichment analyses were carried out. The top 5 ranking pathways for each method are shown in Figure 5. The enriched pathways were mainly related to immunology, including T cell differentiation, T cell activation, negative thymic T cell selection, and immunological synapse.

CORRELATION BETWEEN THE IMMUNE-RELATED GENE SIGNATURE AND IMMUNE CELL INFILTRATION:

Since the most enriched pathways of the 10 hub immune-related genes were correlated with immunology, the association between the risk scores of this novel signature and immune cells infiltration was further explored. Spearman correlation tests indicated that the risk scores of this novel signature was negatively correlated with the infiltration of various immune cells, including T cells (r=−0.470, P<0.001), CD8 T cells (r=−0.294, P<0.001), T helper cells (r=−0.2331, P<0.001), Th1 cells (r=−0.306, P<0.001), TFH cells (r=−0.261, P<0.001), Treg cells (r=−0.398, P<0.001), B cells (r=−0.397, P<0.001), and DC cells (r=−0.279, P<0.001) (Figure 6A). Comparison of immune cells infiltrated fraction between low-risk and high-risk groups obtained similar results. The fraction of many immune cell subtypes in the high-risk group was significantly lower than that in the low-risk group (all P<0.001) (Figure 6B).

VALIDATION OF THE HUB GENES EXPRESSION IN CERVICAL NORMAL AND CANCER TISSUES:

All 10 hub immune-related genes were significantly differentially expressed in cervical normal and cancer tissues. Among them, the mRNA expression of TFRC, ZAP70, SPP1, CHIT1, CAMP, CSF2, and TUBB3 in cervical cancer tissues were significantly higher than that in normal cervical tissues, while LEPR, DLL4, and DES presented the opposite expression (Figure 7A). We also investigated the protein expression of TFRC, CAMP, ZAP70, TUBB3, and CHIT1 in cervical cancer tissues and normal tissues (Figure 7B). TFRC, CAMP, ZAP70, TUBB3, and CHIT1 proteins were highly expressed in cervical cancer tissues. These results were in accordance with the previous bioinformatics analysis based in TCGA database.

Discussion

Although the HPV vaccination and cervical cancer screening have been popularized, cervical cancer remains the second most common malignant cancer among females [19]. The therapeutic effect of metastatic and recurrent cervical cancer is very limited and the 5-year survival rate is still poor [20]. Thus, identifying effective predictive and therapeutic targets is a hot spot for current research. In recent years, many studies have shown that the expression of immune-related genes and the infiltration of immune cells were strongly associated with the clinical outcomes of a variety of cancers [21]. These findings provided new insights for the development of innovative prognostic markers and precision therapeutic targets for cervical cancer. Although previous studies reported various biomarkers, some studies only explored a single gene and some studies were only based on a single microarray dataset or database, thus limiting the credibility and comprehensiveness of the research.

In this study, through univariate Cox and LASSO regression analyses, we screened out 10 immune-related genes: TFRC, SPP1, CAMP, CSF2, TUBB3, ZAP70, CHIT1, LEPR, DLL4, and DES. With the rapid development of high-throughput sequencing, the direct relationship between some of these candidate genes and cervical cancer was explored in previous studies. Transferrin receptor (TFRC), a type II transmembrane glycoprotein, serves as an important iron transporter and participates in the onset and progression of many cancers, including cervical cancer [22]. Secreted phosphoprotein 1 (SPP1) is a secreted phosphoprotein that plays a critical role in the proliferation, migration, and chemotaxis of cancer cells [23]. Higher expression of SPP1 was strongly related to worse overall survival of cervical cancer patients [24]. Cathelicidin antimicrobial protein (CAMP) is mainly expressed in different host defense cells, including macrophages, neutrophils, and endothelial cells. A number of studies have found that CAMP conferred a tumorigenic effect in various cancers, which could promote tumor growth and invasion through recruiting immune cells [25,26]. Colony-stimulating factor 2 (CSF2) is an immune modulatory cytokine stimulating stem cells to produce granulocytes and monocytes. The roles of CSF2 in tumor progression vary in different cancers [27]. In cervical cancer, CSF2 upregulation may induce tumor-related leukocytosis and result in poor prognosis [28]. Tubulin-beta3 (TUBB3) is overexpressed in the more aggressive and resistant cancer phenotypes, which was closely related with the sensitivity of antimicrotubular chemotherapeutic agents [29–31]. Our study produced similar results. The expression of TFRC, SPP1, CAMP, CSF2, TUBB3, ZAP70, and CHIT1 in cervical cancer tissues were significantly higher than that in normal cervical tissue. The opposite results were observed in the expression levels of LEPR, DLL4, and DES. All of these 10 immune-related genes showed diagnostic and prognostic value in cervical cancer.

Nevertheless, a single biomarker may omit large amounts of information and it is difficult to evaluate the prognosis comprehensively. Therefore, a novel predictive signature incorporating these 10 immune-related genes was constructed in this study. Compared with those single biomarker, this signature provided a more reliable and more comprehensive assessment in cervical cancer. It was proved that the risk score of this novel signature could be an independent prognostic predictor with high sensitivity and specificity. The signature successfully classified patients into high-risk and low-risk subgroups with different survival outcomes. Patients with high risk scores exhibited worse overall survival than those with low risk scores. In addition, the risk score of this novel signature presented higher net benefit than traditional clinical characteristics, including age, BMI, clinical stage, distant metastatic status, lymph node status, and histological type. The model generated by both clinical characteristics and this novel signature was the best in predicting the overall survival of cervical cancer patients. Furthermore, the model had the ability to provide a quantifiable assessment of the survival risk per individual patient. This is of great significance in clinical application in terms of individual treatment of cervical cancer.

Further gene ontology enrichment analyses indicated that these 10 genes were significantly enriched in immune-related pathways, including T cell differentiation, T cell activation, and immunological synapse. Therefore, the association between the novel signature and immune cells infiltration in cervical cancer was explored in our study. The results showed that the risk scores of this novel signature had significant negative correlations with the infiltrated statuses of most immune cell types, especially T cell subsets. Compared with the low-risk group, the proportions of infiltrated immune cells were significantly lower in the high-risk group. Compelling evidence has demonstrated that the infiltration of immune cells has an important effect on patients’ prognoses in a number of human malignancies [32–34]. T cells in particular are considered a major element of the tumor microenvironment. Intratumoral T cells display dysfunctional states in mouse models and cancer patients, which is shaped by multifaceted distinct immunosuppressive mechanisms [35]. A major mechanism is the increased expression of inhibitory receptors on T cell surfaces, such as programmed cell death-1 (PD-1), T cell immunoglobulin mucin receptor 3 (TIM3), lymphocyte activation gene 3 (LAG3), and cytotoxic T lymphocyte-associated protein 4 (CTLA4) [36–38]. Another is the reduced production of cytokines for T cells activation, such as interleukin-2 (IL-2), tumor necrosis factor (TNF), β-chemokines, and interferon-γ (IFNγ) [39,40]. Low infiltration of intratumoral T cells predicted poor overall survival and recurrence-free survival [41]. Our findings are consistent with these previous studies. Hence, it is reasonable to infer that our novel signature has potential for use as a predictive and therapeutic target for immune escape in cervical cancer.

There were several limitations in the present study. Firstly, this was a retrospective study. More multicenter prospective clinical cohort studies are needed to validate the prognostic value of our novel signature. Secondly, a small portion of data in the TCGA cohort was missing, which potentially hindered a more comprehensive analysis. Last but not least, the underlying biological mechanisms of the novel signature and immune cell infiltration need to be investigated by in vivo and in vitro experiments.

Conclusions

This novel signature incorporating 10 immune-related genes could provide a more comprehensive assessment in cervical cancer. It would be a valuable biomarker to individualize the risk assessment and predict the prognosis and immune infiltration for each patient. In the era of precision medicine, this novel signature may hold great promise to assist treatment selection and prognostic stratification in cervical cancer.