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

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 1A,1B).

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 (Figure 2A,2B). Moreover, as shown in Figure 2C,2D, 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 (Figure 3A). KEGG results showed that these genes were concentrated in hsa04060 (chemokine signaling pathway) and hsa04060 (cytokine-cytokine receptor interaction) (Figure 3B).

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 (Figure 4A). Module A, including 33 edges and 9 nodes in the network, contained CCL13, CXCL9, CXCL13, CCL19, CXCL5, CCL21, CCR2, CCR5, and CCL5 (Figure 4B). Module B had 6 nodes and 15 edges: CD3G, LCP2, CD3D, CD247, CD3E, and ITK were the main genes (Figure 4C). Module C had 7 nodes and 15 edges (Figure 4D), 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, Figure 5). 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 6A,6B).

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 (Figure 7). 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 (Figure 8). 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 9).

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.

Acknowledgments

Funding: This work was supported by the “234 Discipline Climbing Program” of the First Affiliated Hospital of Naval Military Medical University (2019YXK014).

Footnote

Reporting Checklist: The authors have completed the STREGA reporting checklist. Available at https://atm.amegroups.com/article/view/10.21037/atm-21-7014/rc

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://atm.amegroups.com/article/view/10.21037/atm-21-7014/coif). The authors report that this article was supported by the “234 Discipline Climbing Program” of the First Affiliated Hospital of Naval Military Medical University (2019YXK014). The authors have no other conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Jiang Y, Wang C, Zhou S. Targeting tumor microenvironment in ovarian cancer: Premise and promise. Biochim Biophys Acta Rev Cancer 2020;1873:188361. [Crossref] [PubMed]
3. Barnes TA, Amir E. HYPE or HOPE: the prognostic value of infiltrating immune cells in cancer. Br J Cancer 2018;118:e5. [Crossref] [PubMed]
4. Tas F, Erturk K. Tumor Infiltrating Lymphocytes (TILs) May be Only an Independent Predictor of Nodal Involvement but not for Recurrence and Survival in Cutaneous Melanoma Patients. Cancer Invest 2017;35:501-5. [Crossref] [PubMed]
5. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
6. Sui J, Xu SY, Han J, et al. Integrated analysis of competing endogenous RNA network revealing lncRNAs as potential prognostic biomarkers in human lung squamous cell carcinoma. Oncotarget 2017;8:65997-6018. [Crossref] [PubMed]
7. Rizzuto I, Stavraka C, Chatterjee J, et al. Risk of Ovarian Cancer Relapse score: a prognostic algorithm to predict relapse following treatment for advanced ovarian cancer. Int J Gynecol Cancer 2015;25:416-22. [Crossref] [PubMed]
8. Ahmed N, Stenvers KL. Getting to know ovarian cancer ascites: opportunities for targeted therapy-based translational research. Front Oncol 2013;3:256. [Crossref] [PubMed]
9. Chen B, Chen W, Jin J, et al. Data Mining of Prognostic Microenvironment-Related Genes in Clear Cell Renal Cell Carcinoma: A Study with TCGA Database. Dis Markers 2019;2019:8901649. [Crossref] [PubMed]
10. Huang S, Zhang B, Fan W, et al. Identification of prognostic genes in the acute myeloid leukemia microenvironment. Aging (Albany NY) 2019;11:10557-80. [Crossref] [PubMed]
11. Zhang FP, Huang YP, Luo WX, et al. Construction of a risk score prognosis model based on hepatocellular carcinoma microenvironment. World J Gastroenterol 2020;26:134-53. [Crossref] [PubMed]
12. Rousalova I, Krepela E. Granzyme B-induced apoptosis in cancer cells and its regulation Int J Oncol 2010;37:1361-78. (review). [PubMed]
13. Kontani K, Sawai S, Hanaoka J, et al. Involvement of granzyme B and perforin in suppressing nodal metastasis of cancer cells in breast and lung cancers. Eur J Surg Oncol 2001;27:180-6. [Crossref] [PubMed]
14. Hu SX, Wang S, Wang JP, et al. Expression of endogenous granzyme B in a subset of human primary breast carcinomas. Br J Cancer 2003;89:135-9. [Crossref] [PubMed]
15. D'Eliseo D, Pisu P, Romano C, et al. Granzyme B is expressed in urothelial carcinoma and promotes cancer cell invasion. Int J Cancer 2010;127:1283-94. [Crossref] [PubMed]
16. Ko YH, Park S, Jin H, et al. Granzyme B leakage-induced apoptosis is a crucial mechanism of cell death in nasal-type NK/T-cell lymphoma. Lab Invest 2007;87:241-50. [Crossref] [PubMed]
17. Im SJ, Hashimoto M, Gerner MY, et al. Defining CD8+ T cells that provide the proliferative burst after PD-1 therapy. Nature 2016;537:417-21. [Crossref] [PubMed]
18. Havenar-Daughton C, Lindqvist M, Heit A, et al. CXCL13 is a plasma biomarker of germinal center activity. Proc Natl Acad Sci U S A 2016;113:2702-7. [Crossref] [PubMed]
19. Yang M, Lu J, Zhang G, et al. CXCL13 shapes immunoactive tumor microenvironment and enhances the efficacy of PD-1 checkpoint blockade in high-grade serous ovarian cancer. J Immunother Cancer 2021;9:e001136. [Crossref] [PubMed]
20. Gurbuz F, Desai S, Diao F, et al. Novel inactivating mutations of the DCAF17 gene in American and Turkish families cause male infertility and female subfertility in the mouse model. Clin Genet 2018;93:853-9. [Crossref] [PubMed]
21. Ali A, Mistry BV, Ahmed HA, et al. Deletion of DDB1- and CUL4- associated factor-17 (Dcaf17) gene causes spermatogenesis defects and male infertility in mice. Sci Rep 2018;8:9202. [Crossref] [PubMed]
22. Lee CG, Ren J, Cheong IS, et al. Expression of the FAT10 gene is highly upregulated in hepatocellular carcinoma and other gastrointestinal and gynecological cancers. Oncogene 2003;22:2592-603. [Crossref] [PubMed]
23. Gong P, Canaan A, Wang B, et al. The ubiquitin-like protein FAT10 mediates NF-kappaB activation. J Am Soc Nephrol 2010;21:316-26. [Crossref] [PubMed]
24. Yuan R, Wang K, Hu J, et al. Ubiquitin-like protein FAT10 promotes the invasion and metastasis of hepatocellular carcinoma by modifying β-catenin degradation. Cancer Res 2014;74:5287-300. [Crossref] [PubMed]
25. Yin M, Li X, Tan S, et al. Tumor-associated macrophages drive spheroid formation during early transcoelomic metastasis of ovarian cancer. J Clin Invest 2016;126:4157-73. [Crossref] [PubMed]
26. Drakes ML, Stiff PJ. Regulation of Ovarian Cancer Prognosis by Immune Cells in the Tumor Microenvironment. Cancers (Basel) 2018;10:302. [Crossref] [PubMed]
27. Ovarian Tumor Tissue Analysis (OTTA) Consortium. Dose-Response Association of CD8+ Tumor-Infiltrating Lymphocytes and Survival Time in High-Grade Serous Ovarian Cancer. JAMA Oncol 2017;3:e173290. [Crossref] [PubMed]
28. Truxova I, Kasikova L, Hensler M, et al. Mature dendritic cells correlate with favorable immune infiltrate and improved prognosis in ovarian carcinoma patients. J Immunother Cancer 2018;6:139. [Crossref] [PubMed]
29. Muqaku B, Pils D, Mader JC, et al. Neutrophil Extracellular Trap Formation Correlates with Favorable Overall Survival in High Grade Ovarian Cancer. Cancers (Basel) 2020;12:505. [Crossref] [PubMed]

(English Language Editor: K. Brown)