A 25 Immune-Related Gene Pair Signature Predicts Overall Survival in Cervical Cancer.

Mounting evidence suggests that the tumor microenvironment plays an important role in the occurrence and development of cancer, with immune system dysfunction being closely related to malignant cancers. We aimed to screen immune-related genes (IRGs) to generate an IRG pair (IRGP)-based prognostic signature for cervical cancer (CC). Datasets were obtained from The Cancer Genome Atlas and Gene Expression Omnibus databases and used as training and validation cohorts, respectively. Using the ImmPort database, IRGs in control and CC samples were compared, and differentially expressed genes were identified to construct an IRGP prognostic signature. Based on this analysis, 25 IRGPs were identified as important factors for the prognosis of CC. Univariate and multivariate Cox regression analyses further showed that the IRGP signature was an independent prognostic factor of overall survival. In summary, we successfully constructed an IRGP prognostic signature of CC, providing insights into immunotherapy for CC.

Introduction

Cervical cancer (CC) is among the commonest malignant tumors and the deadliest malignancy in females, with about 604 000 new and 342 000 death cases reported worldwide in 2020. In recent years, CC incidence has declined due to increased awareness; however, its mortality rate remains high in developing countries. CC is difficult to diagnose at an early stage, limiting treatment only to options available in the advanced stage. Currently, the standard treatment options for CC include surgery, chemotherapy, and radiotherapy. However, the effectiveness of chemotherapy in CC is not remarkable owing to drug resistance. Therefore, developing better prognostic tools for predicting and improving the outcome of CC cases is highly needed.

Immune-related genes (IRGs), with important regulatory effects on the immune system, are reportedly associated with the occurrence and development of several cancers.[1,3] The key characteristic of IRGs is that they regulate complex modulatory networks in cancers, thereby making them therapeutically significant targets that could constitute biomarkers for predicting prognosis in cancer. Specifically, tumor-infiltrating immune cells are highly relevant to CC development.

Recently, immune-related prognostic signatures using microarray and RNA-sequencing data have been developed for several cancers. For example, a robust IRG pair (IRGP) signature was identified for predicting prognosis and immune heterogeneity in patients with glioblastoma. Moreover, a 14 IRG-based prognostic model was constructed for overall survival (OS) prediction in gastric cancer using The Cancer Genome Atlas (TCGA) and ImmPort databases. Other immune-related signatures have also been identified for several malignancies, including hepatocellular carcinoma, osteosarcoma, ovarian carcinoma and colorectal cancer.

In the present study, we generated an IRGP signature by combining RNA-sequencing (RNA-seq) data and clinical findings of CC cases from TCGA and verifying them using Gene Expression Omnibus (GEO). The constructed prognostic model may help predict CC outcomes with improved accuracy.

Materials and Methods

Data collection and processing

The clinical and RNA-seq data of CC cases (n = 309; comprising 3 control and 306 tumor specimens) were retrieved on April 20, 2021 from TCGA database (https://portal.gdc.cancer.gov/) and were utilized as the training cohort for constructing a prognostic signature. The GSE44001 (n = 300) dataset was downloaded on April 26, 2021 from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) and served as the validation cohort for verifying the efficacy of the IRGP signature. IRGs were obtained from ImmPort (https://www.immport.org/shared/genelists), a large public repository providing immunological data for humans. Gene expression profiling from RNA-seq data underwent normalization to Fragments per Kilobase Million, and samples with a value of 0 (zero) were deleted. Accordingly, gene expression profiling data in the GSE44001 dataset underwent conversion with probes changed to related gene symbols. In both training and validation cohorts, individuals with incomplete survival data were excluded from the analysis.

Construction of an IRGP prognostic signature

In total, 2483 IRGs were acquired from ImmPort; subsequently, differentially expressed IRGs in CC were identified for constructing the IRGP prognostic signature. With median absolute deviation >0.5, IRGs expressed in both training and validation cohorts showing a large variation were further analyzed; otherwise, they were excluded. IRGs were paired to form IRGPs, which were scored according to the relative expression of the 2 paired IRGs, including IRG1 and IRG2. In case IRG1 showed elevated expression compared with IRG2 in a particular IRGP, the output was scored as 1; otherwise, the output was considered 0. Each IRGP with small variation and unbalanced distribution, for example, fluctuation range showing a score of 0 or 1, was removed, considering the remaining IRGPs as the initial candidates. Univariate, Lasso, and multivariate regression analyses were performed for determining the correlations between the expression of IRGPs and OS in the training cohort, and those with P < .05 were selected. The log-rank test was used to identify associations of OS with IRGP in the training cohort. Then, the Lasso penalized Cox regression analysis (1000 iterations) was performed to avoid overfitting of the data, to obtain IRGPs and construct a prognostic signature. Finally, a total of 25 IRGPs were identified for prognostic signature building. A 1-year time-dependent receiver operating characteristic (ROC) curve for the training cohort was designed. The optimal cut-off, where sensitivity + specificity was largest, was determined, and the CC cases were classified as high- or low-risk based on this value. Using these coefficients for prognostic IRGPs, a model was established to calculate the IRGP risk score for every specimen. The risk score signature was derived as:

where Expr and Coef are the relative expression and Lasso regression-derived coefficient of gene pairs, respectively.

Validation of the IRGP signature

For predicting the reliability and stability of the developed IRGP prognostic signature, R (survival package) software was used for Kaplan-Meier curve generation based on high- and low-risk cases in both cohorts. Then, univariable and multivariable Cox proportional-hazard analyses of other clinical factors were performed for evaluating the accuracy of the prognostic signature.

Assessment of the immune microenvironment

CIBERSORT analysis, an instrumental learning method that estimates cell type abundances from a large tissue transcriptome, was performed for the enrichment analysis of immune cells in high- and low-risk patients and the enrichment of 22 different immune cells in CC samples was analyzed. The normal samples and low-expression genes were removed, and the remaining data were corrected. The presence of immune cells in high- and low-risk cases in the training cohort was analyzed using CIBERSPORT with approximately 1000 permutations.

Gene Ontology (GO) and gene set enrichment analysis (GSEA)

To elucidate the biological mechanism of action of IRGPs in CC, GO analysis and GSEA were performed. R (gProfiler) was used for the GO analysis. GSEA was performed for the high- and low-risk cases using R (fgsea package), with 1000 permutations. Genes with a false discovery rate of <0.05 were retained.

Statistical analysis

R (v3.6.3) and Perl (v5.30.0) interfaces were used to perform the abovementioned statistical analyses. The glmnet, survival, survival ROC, and limma packages of R were used to perform the univariable and multivariable risk regression model analyses, to draw ROC curves, and to plot survival curves, respectively. For all results, P < .05 indicated statistical significance.

Results

Samples from the TCGA database (n = 306, excluding the normal samples) were analyzed as the training cohort, and those from the GEO (n = 300) database were analyzed as part of the validation cohort. A total of 2483 IRGs were retrieved from ImmPort. Overall, 379 IRGs were expressed in both training and validation sets. Subsequently, 13 033 IRGPs were used for constructing a prognostic signature. The associations of IRGs with OS were examined in the training group by univariable regression analysis. A total of 25 IRGPs (Table 1) had significant associations with OS in CC. The risk scores for various samples were calculated in the training group according to the developed risk score formula. A 1-year time-dependent receiver operating characteristic (ROC) curve was plotted for the training cohort dataset, and an optimal cut-off of 1.274 was obtained (Figure 1), according to which the training cohort was assigned to high- and low-risk cases. High-risk cases showed remarkably decreased OS rates compared with low-risk cases (P < .01) (Figure 2(a)). Next, univariable and multivariable Cox regression analyses were performed for evaluating the effects of IRGP risk score, age, and tumor grade and stage on the time-independent ROC curve of the risk score in the training cohort. Univariable analysis showed a significant association of the IRGP prognostic signature with OS (hazard ratio [HR] = 4.466, 95% confidence interval [CI]: 2.540-7.852; P < .01) (Figure 3(a)). Multivariable analysis further indicated that the IRGP prognostic signature was a prognostic factor, independently of age, and tumor grade or stage (HR = 4.958, 95% CI: 2.655-9.258; P < .01) (Figure 3(b)).

Table 1.

Detailed information on about the 25 immune-related gene pairs (IRGPs).

IRGP1	Immune processes	IRGP2	Immune processes	Coef
ADRM1	Antigen_Processing_and_Presentation	MIF	Antimicrobials	−0.179492135
APOBEC3H	Antimicrobials	BTC	Cytokines	−0.463445295
C5AR1	Chemokine_Receptors	STC1	Cytokines	−0.181534534
CXCL14	Cytokines	ANGPTL2	Cytokine_Receptors	−0.000291734
CXCL14	Cytokines	PPP3CB	NaturalKiller_Cell_Cytotoxicity	−0.005971396
CXCL2	Cytokines	RAF1	NaturalKiller_Cell_Cytotoxicity	0.086843698
DES	Antimicrobials	EPOR	Cytokine_Receptors	−0.085941485
DLL4	Antimicrobials	DES	Antimicrobials	0.449944957
DUOX1	Antimicrobials	NRP1	Cytokine_Receptors	−0.332649118
FLT3LG	Cytokines	INHBA	Cytokines	−0.199692472
HLA-DQA2	Antigen_Processing_and_Presentation	CCL3	Antimicrobials	−0.213899733
IL1B	Antimicrobials	CD3D	TCRsignalingPathway	−0.017464468
IL1B	Antimicrobials	DUOX1	Antimicrobials	0.394694083
IL1B	Antimicrobials	EDN1	Chemokines	0.737385178
IL34	Cytokines	OSM	Cytokines	−0.213905919
INHBA	Cytokines	PRKCQ	TCRsignalingPathway	0.059166996
IRF5	Antimicrobials	LIF	Cytokines	−0.047922745
JAK1	Antimicrobials	APOBEC3C	Antimicrobials	0.401495096
LMBR1	Antimicrobials	MAP3K14	TCRsignalingPathway	0.064208499
LTBP3	Cytokines	MAP3K14	TCRsignalingPathway	0.5265958
MICA	NaturalKiller_Cell_Cytotoxicity	IFIH1	Antimicrobials	0.354710175
NRP1	Cytokine_Receptors	THRA	Cytokine_Receptors	0.204863519
PLXNB3	Cytokine_Receptors	FGFR2	Cytokine_Receptors	0.386753047
TLR3	Antimicrobials	CXCR6	Antimicrobials	0.682446531
VAV3	BCRSignalingPathway	NRP1	Cytokine_Receptors	−0.41334932

Figure 1.

Time-dependent receiver operating characteristics (ROC) curve for the immune-related gene pair (IRGP) risk score in the training set. An IRGP risk score of 1.247 was considered the optimal cut-off for classifying cases into the high- and low-risk groups.

Figure 2.

Overall survival (OS) for different risk cases according to the optimal cut-off 1.247 in training and validation. Kaplan-Meier survival curves for high- and low-risk cases in the training (a) and validation (b) cohorts.

Figure 3.

Associations within the immune-related gene pair (IRGP) prognostic signature and clinical data with OS in the training and validation cohorts. Univariable and multivariable Cox analyses of clinical parameters and the IRGP prognostic signature in the training (a b) and validation (c d) cohorts.

IRGP signature is a valid independent prognostic factor

An independent validation set, GSE44001 (n = 300), was used for verifying the consistency of the prognostic value of the developed IRGP signature. Risk scores were calculated for various patients of the validation set, who were assigned to the high- and low-risk groups based on the optimal cut-off value obtained in the training cohort. In the validation cohort, high-risk cases had decreased OS compared with their low-risk counterparts (P = .02). Univariable and multivariable Cox regression analyses demonstrated that the IRGP prognostic signature indeed constituted an independent prognostic factor (univariable: HR = 1.542, 95% CI: 1.125-2.111; P = .007; multivariable: HR = 1.580, 95% CI: 1.141-2.189; P = .006) (Figure 3(c) and (d)).

Infiltration of immune cells in high- and low-risk cases

Based on the CIBERSORT algorithm, 22 types of immune cells identified within patients with CC were systematically analyzed in both risk groups for proportions analysis. Seven immune cell types showed markedly different proportions between the high- and low-risk groups (P < .05) (Figure 4(a)). In particular, resting dendritic cells, resting mast cells, activated CD4+ memory T cells, and CD8+ T cells had elevated infiltration levels in low-risk cases compared with those in high-risk cases, whereas M0 macrophages, activated mast cells and neutrophils had lower proportions. M0 macrophages and CD8+ T cells were remarkably enriched in the high- and low-risk groups, respectively (both P < .01 Figure 4(b)).

Figure 4.

Immune cell infiltration levels in low- and high-risk cervical cancer cases: (a) In total, 22 distinct immune cells were examined for abundance based on CIBERSORT algorithm in low- and high-risk cases, and (b) Abundance levels of 7 immune cells subsets showed significant differences between the 2 patient groups.

GO and GSEA findings

Gene Ontology (GO) enrichment analysis showed that the IRGPs were significantly enriched mainly in 17 signaling pathways, including immune response regulating, adaptive immune response, and lymphocyte-mediated pathways (Figure 5(a)). Gene set enrichment analysis (GSEA) was used for exploring significantly altered signaling pathways in high-risk and low-risk CC cases of the training cohort. The most significant signaling pathways were identified according to normalized enrichment scores (Figure 5(b)). Several related pathways, including immune response regulating signaling pathway, adaptive response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains, lymphocyte-mediated immunity and antigen receptor-mediated signaling pathway, were significantly enriched in low-risk cases, which suggested that these signaling pathways may have key functions in CC development.

Figure 5.

Functional enrichment of 25 immune-related gene pairs: (a) results of Gene Ontology (GO) enrichment analysis; the immune-related gene pairs (IRGPs) were significantly enriched in 17 signaling pathways and (b) Results of gene set enrichment analysis (GSEA); immune response regulating signaling pathway, adaptive response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains, and lymphocyte-mediated immunity and antigen receptor-mediated signaling pathway were significantly enriched in low-risk cases.

Discussion

Considering the important role of tumor immunity in CC, the TCGA training cohort was analyzed to establish a prognostic signature involving 25 IRGPs containing 46 unique IRGs for OS prediction in CC cases. The reliability of the prognostic signature was confirmed with an independent external validation cohort from the GEO database. Univariable and multivariable analyses in the training and validation sets showed that age and the tumor pathological stage were not independent OS predictors in CC; however, the independent prognostic ability of the constructed IRGP signature was proven. Hence, this IRGP prognostic signature may aid in predicting the survival outcomes of patients with CC.

IRGs have different roles in tumors. For instance, CXCL14, a highly conserved homeostatic chemokine, is responsible for immune cell recruitment and maturation and drives epithelial cell movement helping to establish immune monitoring in the normal epithelium. IL-1B gene polymorphisms are associated with CC risk in the Chinese Uygur population. Moreover, DLL4 inhibits radiation resistance and metastasis in CC and can potentially be used as a biomarker for predicting radiosensitivity and prognosis in CC. VAV3 promotes gastric cell proliferation and significantly contributes to prostate cancer growth and malignancy. STC1 is a glycoprotein hormone involved in calcium/phosphorus homeostasis, being tightly associated with tumor occurrence and development. Indeed, RNA interference-mediated inhibition of STC1 expression in CaSki cells (a human papillomavirus type 16-positive cell line) cell growth, migration and invasion were greatly enhanced. Furthermore, reduced levels of d Nrp1+ Treg cells in patients with CC are directly related to reduced tumor mass. High CXCL12/CXCR4 and CXCL16/CXCR6 ratios in cervical intraepithelial neoplasia and CC suggest that the occurrence of CC is a continuous process. Indeed CXCR6 is considered a molecular marker and a potent prognostic factor of CC. In addition, Nrp1 is reported as a new tumor marker in hepatocellular carcinoma ; it also regulates gastric cancer progression and enhances the invasive and migratory properties of lung adenocarcinoma cells. The INHBA gene plays an important role in many tumors, with high INHBA levels in esophageal squamous cell carcinoma and other solid tumors reflecting poor prognosis. Moreover, INHBA overexpression results in poor clinical outcomes in bladder urothelial carcinoma, indicating its potential as a prognostic marker and a new therapeutic target. INHBA expression was associated with colon cancer prognosis and recurrence. High CD3D amounts are tightly associated with low survival in breast cancer.

Human papilloma virus (HPV)-related persistent infection is the main cause of CC. The application of the HPV vaccine involves the alteration of the tumor immune microenvironment to prevent CC development. APOBEC3 family genes are translated into cytidine deaminases that counter viral infection and retrotransposition ; the haplotype II APOBEC3H, with higher frequency in individuals of African descent, is translated into a protein exerting the greatest antiviral effects on cells. Conventional CC treatment options include chemotherapy, radiotherapy, and surgery. However, patients with advanced-stage tumors are susceptible to radiotherapy and chemotherapy resistance. Many CC cases are detected at an advanced stage, usually with a high invasion rate; thus, the 3-year and 5-year mortality rates range between 52% and 79%.[30-32] Although some diagnostic and prognostic markers of CC have been identified, the overall expression of these markers remains to be further elucidated owing to CC heterogeneity.

Previous studies have shown that the balance of the host tumor microenvironment is critical to tumor occurrence and development. The host tumor microenvironment consists of fibroblasts, vascular endothelial cells, immune cells, cytokines, growth factors, hormones, and extracellular matrix components, among other elements. In this study, resting dendritic cells, resting mast cells, activated CD4+ memory T cells, and CD8+ T cells showed elevated infiltration levels in low-risk cases compared with those in high-risk cases, whereas M0 macrophages, activated mast cells, and neutrophils showed lower infiltration. Most of the retrieved IRGs were enriched in the immune response regulating and lymphocyte- and antigen receptor-mediated signaling pathways. Moreover, there were significant differences in the proportions of specific immune cells between high- and low-risk cases. Typically, CD8+ T cells recognize tumor cells and substantially contribute to the immune response, and elevated CD8+ T cell infiltration could lead to improved prognosis. As shown above, CD8+ T cell infiltration was significantly lower in high-risk cases, in disagreement with previously reported findings in papillary carcinoma. Although CD8+ T cells are one of the main anticancer immune cell types, they are usually in a dysfunctional state (T cell exhaustion) when they infiltrate the cancer tissue. The characteristics of CD8+ T cell exhaustion are decreased activity and proliferation, increased apoptosis rate, and decreased production of effector cytokines. The distribution of CD8+ cells may be influenced by inhibitory receptors, small immunosuppressive molecules, immune regulatory cells, transcriptomic changes, and metabolic reprogramming of exhausted CD8+ T cells. Studies have also shown anomalous dendritic cells contribute to CD8+ T cell generation, resulting in a high CD8+ T cell abundance. However, these cells do not exert the corresponding antitumor effects. Many tumor immunotherapies, such as inhibitory antibodies, target CD8+ T cells, achieving good therapeutic effects. Understanding the molecular mechanism of CD8+ T cell exhaustion is very important for establishing reasonable immunotherapeutic interventions.

An IRGP prognostic signature is generated based on a relative expression ranking of genes and does not need data standardization, offering a great advantage over previously determined prognostic profiles. In addition, it is easy to be applied to clinical practice. Recently, other researchers used a 29-IRGP signature for predicting OS in CC cases; however, the current study included more cases and used more stringent filtering conditions to develop a 25-IRGP signature for predicting patient prognosis in CC. Notably, the IRGPs screened and retained for model building in the present study also differed from those examined in previous studies.

If IRGs are detected in a patient, the 25 IRGP-based prognosis signature may be employed for assessing their risk score based on the risk score formula. Moreover, the risk score can be used to determine high- or low-risk cases. Finally, a Kaplan–Meier curve can be plotted for patient prognosis. The IRGP prognosis signature has been widely used in many cancer studies, in particular in hepatocellular carcinoma and osteosarcoma. Moreover, all of them were verified using data in the databases; and the only shortcoming was that these studies were not experimentally verified.

This study had several limitations. The initial features of the prognostic signature are gene expression files based on RNA-seq or microarray data, which have a high cost, a lengthy transformation cycle, and requirements of professional knowledge of biological information; hence, popularizing the model for clinical use would be rather challenging. Furthermore, the biological functions and molecular mechanisms of the 46 IRGs in patients with CC need to be evaluated. Finally, standard genetic analysis by polymerase chain reaction is required to further verify the expression of the identified IRGs in CC.

Conclusion

This study presents a 25-IRGP prognostic signature as a reliable and independent biomarker of OS prognosis in CC. It also provides insights into the potential mechanisms of immunotherapy and the identification of related targets, which may contribute to the development of personalized immunotherapeutic strategies, thereby improving the outcome of CC.