Integrated analysis of prognostic immune-related genes in the tumor microenvironment of ovarian cancer.

Background: Ovarian cancer (OC) is a major cause of most gynecological cancer deaths, and the rates of incidence and mortality are increasing worldwide. However, factors in the tumor microenvironment (TME) related to OC and certain prognostic markers of OC are still unknown. We aimed to identify biomarkers connected to prognostic immunity based on clinical patients' data from The Cancer Genome Atlas (TCGA).
Methods: We used the ESTIMATE algorithm to compute the immune and matrix scores of OC patients from TCGA. Next, differentially expressed genes (DEGs) according to the immune and matrix scores were obtained. Subsequently, genes (GZMB, C2orf37, CXCL13, and UBD) connected with prognostic immunity were determined. Moreover, functional enrichment analysis and the protein-protein interaction network showed that these genes were enriched in many biological processes related to immune function. The Tumor Immune Estimation Resource (TIMER) algorithm was also used to analyze the immune prognostic genes according to six immuno-infiltrating cells.
Results: According to high/low immune-scores and matrix-score groups, 682 common genes were identified, within 420 upregulated genes and 262 downregulated genes. Gene ontology (GO) analysis of biological process primarily enriched in T cell activation, regulation of lymphocyte activation and lymphocyte differentiation. OS analysis showed 45 genes (6.6%) were relevant in the final results. The Kaplan-Meier plotter database verified the top 10 genes related to prognosis, but only GZMB, C2orf37, CXCL13 and UBD were related to overall survival (OS). Conclusions: GZMB, CXCL13, and UBD may influence prognosis via their effects on the infiltration of immune cells and therefore represent potential targets for OC immunotherapy. 2022 Annals of Translational Medicine. All rights reserved.

Introduction

According to global cancer data in 2020, the incidence and mortality rates for ovarian cancer (OC) ranked 8th among female malignant tumors (1). OC is considered the deadliest malignant tumor in the female reproductive system because of its hidden onset and difficult early diagnosis. The ovarian epithelial tumor is the most common pathological type of OC, accounting for 80–95% of ovarian malignancies. The latest research showed the 5-year survival rate of epithelial OC as <40% (2). Despite active cytoreductive surgery, 80% of patients still relapse within 2 years, and recurrent OC lacks a valid treatment. The failure to effectively eradicate OC can be attributed to the complex interconnected signal network coupled with the unique peritoneal tumor microenvironment (TME). Previous studies have shown that high levels of immune-cell infiltration are associated with favorable outcomes (3,4), so exploring the TME to find prognostic immune-related biomarkers of OC should yield beneficial therapeutic targets. We present the following article in accordance with the STREGA reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-21-7014/rc).

Methods

Gene expression datasets

Data of gene expression profiles and clinical characteristics including grade, staging, and survival time, as well as the outcome of OC patients were obtained from The Cancer Genome Atlas (TCGA) portal (https://tcga-data.nci.nih.gov/tcga/). Immune and matrix scores were computed by ESTIMATE, which provided the matrix cells and immune cell infiltration levels from the tumor tissues (5). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Correlation analysis of expressed genes

On the basis of the ESTIMATE results, all specimens were divided into high and low immune/matrix score groups to determine the overlapping genes. |Log (FC)| >1, P value <0.05, and false discovery rate (FDR) <0.05 (6) were defined as the limits of significant differentially expressed genes (DEGs). Heatmaps were generated by the heat map package in R.

Survival analysis

Kaplan-Meier (K-M) survival estimation was applied to determine the relevance between immune/matrix scores and overall survival (OS) of patients, and then to recognize possible prognostic genes.

Functional analysis

In this study, to make certain the feasible function of genes, we used gene ontology (GO) analysis and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis using the cluster Profiler package in R, and a P value <0.05 was defined as a meaningful result. In addition, in order to analyze interacting genes, STRING was used to establish a protein-protein interaction (PPI) network, which was reconstructed in Cytoscape. Finally, using Molecular Complex Detection (MCODE) in Cytoscape, clustering was obtained based on the densely connected region with topology.

TME analysis

The Tumor Immune Estimation Resource (TIMER), including more than 10,000 samples of 32 cancer types from TCGA, is a comprehensive resource for systematic analysis of different immuno-infiltrating cancer types (https://cistrome.shinyapps.io/timer/). We analyzed the prognostic level of OC genes and their correlations with immune cell infiltration, including B cells, macrophages, CD4+ T cells, neutrophils, CD8+ T cells, and dendritic cells.

Statistical analysis

Differential analysis of expressed genes, function analysis and unsupervised clustering analysis were performed using R version 3.6.

Results

Relationship between immune and matrix scores with OC prognosis

We discovered that the immune score was significantly associated with OC prognosis. In this research, 371 OC datasets were obtained from TCGA database containing gene expressions and clinical characteristics. According to the ESTIMATE algorithm, matrix scores ranged from −2,067.71 to 1,568.51, and immune scores were distributed between −1,781.66 and 2,529.21. This research showed that the distribution of matrix scores (P=0.225, P=0.331) and immune scores (P=0.628, P=0.378) did not vary across the different stage levels and different clinical grade levels (data not shown). In order to explore the possible connection between the prognosis of OC and the immune and matrix scores, the patients were divided into two groups with high/low-scores, and K-M survival prognostic curves were constructed. The immune scores were markedly positively correlated with OS, indicating a positive factor in patient prognosis ().

Figure 1

Relationship between immune and stromal scores and ovarian cancer prognosis. Patients’ OS based on (A) immune scores and (B) stromal scores. OS, overall survival.

Analysis of DEGs

In this study, 629 upregulated and 520 downregulated immune-related genes were selected by comparing the different immune-score groups. Similarly, 734 upregulated genes and 398 downregulated genes were in the stromal-score groups (). Moreover, as shown in , there were 420 up- and 262 downregulated overlapping genes in the immune- and matrix-score groups.

Figure 2

Comparison of gene expression and ovarian cancer patients’ immune and stromal scores. Heatmap of significant differentially expressed genes based on immune and stromal scores (A,B); red indicates higher expressed genes and green indicates lower ones. Venn diagram analysis of aberrantly expressed genes based on immune and stromal scores (C,D).

Functional analysis of prognostic genes

The functions of the overlapping genes were predicted, and 40 pathways and 635 GO results were produced by KEGG and GO analysis. GO analysis includes three aspects: biological process, cellular component and molecular function. The results indicated that 0042110: T cell activation, 0009897: external side of the plasma membrane, 0051249: regulation of lymphocyte activation, 007159: leukocyte cell-cell adhesion, and 0030098: lymphocyte differentiation were the top five from the GO analysis (). KEGG results showed that these genes were concentrated in hsa04060 (chemokine signaling pathway) and hsa04060 (cytokine-cytokine receptor interaction) ().

Figure 3

Enrichment analysis of immune-related genes. (A) GO analysis and (B) KEGG analysis. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process; CC, cellular component; MF, molecular function. X-axis represents gene number.

PPI network construction and analysis

The PPI network was used to display the connections and functions of the immune-related prognostic genes (). Module A, including 33 edges and 9 nodes in the network, contained CCL13, CXCL9, CXCL13, CCL19, CXCL5, CCL21, CCR2, CCR5, and CCL5 (). Module B had 6 nodes and 15 edges: CD3G, LCP2, CD3D, CD247, CD3E, and ITK were the main genes (). Module C had 7 nodes and 15 edges (), containing ITGAM, ITGAX, FCER1G, ITGB2, CD53, CYBB, and LILRB2.

Figure 4

The functional analysis of immune-related genes. (A) Functional analysis of immune-related genes. (B-D) Top three modules in the protein-protein interaction networks.

Survival analysis of overlapping genes

The survival package in R showed 45 genes related to OS selected from 682 overlapping genes, but only the top 10 were showed (P<0.05, ). Next, we predicted their functions. KEGG and GO analysis distinguished 20 pathways and 249 GO terms, showing that the 45 genes mainly participated in regulation of T cell activation (GO: 0050863), leukocyte cell-cell adhesion (GO: 0007159), in cytokine-cytokine receptor interaction (hsa04060) and viral protein interaction with cytokines and cytokine receptors (hsa04061) ().

Figure 5

Kaplan-Meier survival curves for immune-related genes associated with overall survival. Horizontal axis: overall survival time, years, Vertical axis: survival function.

Figure 6

Functional analysis of prognosis-related immune genes. (A) GO analysis and (B) KEGG analysis. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process; CC, cellular component; MF, molecular function.

Validation in the K-M plotter database

Because it is important to test and verify prognostic genes in other databases, we used the K-M plotter database to validate the results. Only 4 genes had a relationship with OC survival. In addition, the functions of GZMB, UBD (GABBR1), and C2orf73 in OC have not been reported in previous research (). Through verification in GEPIA database, we found that the expression levels of CXCL13 and BUD in patients with advanced ovarian cancer were lower than those in early ovarian cancer. Data is not shown in the text.

Figure 7

OS curves of the 10 prognostic genes. (A) CCR1 (SCYAR1), (B) C2orf73, (C) CD80, (D) CXCL13, (E) GBP5, (F) GZMB, (G) HPSE2, (H) IL7R, (I) SLAM1, and (J) UBD (GABBR1). OS, overall survival; HR, hazard ratio.

Correlation of prognostic genes’ expressions with immune cell infiltration level

We investigated whether the prognostic genes’ expressions correlated with immune cell infiltration levels in OC and our results showed that GZMB, CXCL13, and UBD expression levels negatively correlated with tumor purity and significantly correlated with immune cell infiltration levels (). The relation between immune cell (B cells, CD8+ cells, CD4+ cells, macrophages, neutrophils and dendritic cells) infiltration and clinical characteristics was further examined by the survival module of TIMER and the results showed that dendritic cell infiltration could predict better survival ().

Figure 8

Correlation of different gene expressions with levels of immune cell infiltration.

Figure 9

Relationship between clinical characteristics and levels of immune cell infiltration. OV, ovarian cancer.

Discussion

OC has the highest mortality of all gynecological malignancies. During the past 40 years, the OC survival rate has improved but two-thirds of women still die within 10 years of diagnosis (7). Less than 25% of patients are diagnosed in the early stages, because of the lack of specific diagnostic symptoms and/or biomarkers. OC preferentially metastasizes into the peritoneal cavity (8), and the interaction between the transformed cells and the unique peritoneal TME affects not only the outcome of disease, but also the treatment response (2). The lack of effective eradication of OC can be thought to due to the complex interconnected signal networks coupled with the unique peritoneal TME. Thus, it is important to understand the TME of OC patients by identifying immune-related prognostic genes that predict OS. Previous study shown that the more immune cell number predicts the better outcomes, and using the ESTAMATE algorithm could assess the tumor purity and get the matrix and immune scores (3).

The TME plays a key role in tumor progression, metastasis, and recurrence, causing a poor prognosis (9-11), because of the highly immunosuppressive in OC tumor microenvironment (OC-TME), allowing evasion of immune surveillance and unhindered tumor growth. In this research, immune-related DEGs were comprehensively analyzed, and the relationship among DEGs, prognosis, and immune cells was explored. Patients with higher immune scores had a longer OS (P=0.037).

In order to identify the immune-related genes, we classified the samples into high/low immune-score and matrix-score groups. In these, 682 common genes were screened out, comprising 420 up- and 262 downregulated genes. GO analysis revealed the biological processes were primarily enrichment of T cell activation, regulation of lymphocyte activation, and lymphocyte differentiation. The results indicated that immune-related DEGs may promote the evolution of OC by influencing these biological processes.

OS analysis was performed to further determine the prognostic-related immune genes and among 682 relevant immune genes, 45 (6.6%) were relevant in the final result, but we only show the top 10 genes related to prognosis. In addition, PPI modules were established to show the influence and function of the prognostic immune-related genes. Nodes with high connectivity within the modules included GZMB, CD80, CXCL13, SLAMF7, CCR1, and CCR7. Next, we used the K-M plotter database to verify the top 10 genes related to prognosis, and found that only GZMB, C2orf37, CXCL13, and UBD were related to OS. Granzyme B (GZMB) is a serine protease with the highest content secreted by cytotoxic T lymphocytes (CTLs) and natural killer (NK) cells (12). Many studies have shown that GZMB is detected in lung carcinomas (13), primary human breast carcinomas (14), urothelial carcinomas (15), and nasal-type NK/T-cell lymphoma (16). However, no study of this gene has been carried out in OC. The CXC subfamily chemokine 13 (CXCL13) is a member of the CXC chemokine family. CXCL13 can bind to its specific receptor CXCR5 and is secreted by B cells, follicular helper T cells, dendritic cells, and macrophages, and plays an important role in immune cell recruitment, activation, and adaptive immune response regulation (17,18). According to one study, it was been confirmed that the combination of CXCL13, CD8, and CXCR5 is an independent factor for OS and progression-free survival in patients with high-grade OC (19), and high CXCL13 expression is associated with prolonged survival. Our results are consistent with this previous study. C2orf37, also named DCAF17, is involved in gonad development and functions in both sexes (20,21). It is a key protein in the pathogenesis of the multisystem disorder Woodhouse-Sakati syndrome. However, there have been no studies linking DCAF17 to OC. The ubiquitin protein D (UBD) is a member of the ubiquitin protein family. UBD is reported to be a tumor oncogene that promotes tumor recurrence and metastasis, via multiple signaling pathways including PI3K/Akt, Wnt/beta-catenin, and nuclear factor kappa-B signals (22-24). However, our research showed that the higher the UBD expression, the better the outcome for patients. Therefore, we need to conduct more experiments to verify this.

More and more studies indicate the importance of B cells, dendritic cells, CD4+ T cells, neutrophils, CD8+ T cells, macrophages, and other immune cell types in OC (2). Tumor-associated macrophages are considered to be involved in the immunosuppressive process in OC, and play a key role in tumor metastasis, angiogenesis, cell invasion, and early recurrence (25,26). In high-grade OC, high-level infiltration of CD8+ T cells is normally related to a better prognosis (27). A previous study found that dendritic cells play a vital function in initiating and regulating immune responses, with a connection to patients’ final results. Also, mature dendritic cells are associated with increased OC immune soaking and can affect the prognosis (28). High levels of dendritic cell infiltration predict favorable survival outcomes in OC, which is consistent with the results of previous studies. The formation of neutrophil extracellular traps is also related to improvement of OS in high-grade OC (29).

From our study, we found that GZMB, UBD, and CXCL13 had a relationship with tumor purity, but the specific mechanism is unknown. Therefore, our research has the following limitations. First, we urgently need biological experiments to verify our results, because our research is based on data analysis. Second, we lack the molecular mechanisms for these genes, which will require further exploration.

The article’s supplementary files as