An Immune-Related Gene-Based Signature as Prognostic Tool in Ovarian Serous Cystadenocarcinoma

Introduction

Ovarian cancer, a malignancy that arises in ovarian tissue, is classified by the World Health Organization as serous carcinoma, mucinous carcinoma, endometrioid carcinoma, or clear-cell carcinoma, or as Brenner, mixed, or undifferentiated types.¹ Ovarian serous cystadenocarcinoma (OSCC) is an aggressive type of ovarian cancer with high recurrence rate, metastatic ability, mortality, and resistance to conventional therapy.² Although treatment methods for ovarian cancer are gradually improving, the prognosis of advanced serous ovarian neoplasms is still poor, indicating the need for new biomarkers.

Tumor-associated macrophages (TAMS) are the major lineage cells in the ovarian tumor microenvironment (TME) and have been shown to play an important role in immunosuppression through tumor cell invasion, angiogenesis, metastasis, and early relapse.^3,4 Immune cells, mesenchymal cells, endothelial cells, inflammatory mediators, and extracellular matrix (ECM) molecules have also been found in the TME of serous ovarian neoplasms,^5,6 with immune cells playing a defensive role through chemotaxis.^7,8 Based on these findings and given that immune escape is a major cause of tumorigenesis,⁹ we hypothesized that classifying immune cells and estimating their abundance in the ovarian TME would help clarify their involvement in cancer. Therefore, in this study, we used weighted gene co-expression network analysis (WGCNA), univariate Cox regression analysis, and least absolute shrinkage and selection operator (LASSO) analysis to identify immune-related modules and calculate a “signature index”. Time-dependent receiver operator characteristics curve (tROC) analysis was also performed, suggesting that the signature index is a potential immune-related prognostic biomarker for serous ovarian neoplasms.

Materials and Methods

Data Procession

The gene expression profiles of the GSE53963 dataset were extracted from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/).¹⁰ GSE53963 is based on GPL6480 and contains 174 serous ovarian cancer samples, which were used as the training set to create an immune-related gene-based signature. The extracted gene expression profiles were normalized using the normalizeBetweenArrays function in the limma R package.¹¹ RNA sequencing (displayed as read count) and clinical information from patients with ovarian serous cancer were acquired from The Cancer Genome Atlas (TCGA; https://www.cancer.gov/) to validate the prognostic ability of the immune-related gene-based signature.¹² RNA-sequencing data were normalized using the VOOM function in the limma R package.¹¹ If a gene was measured using multiple probes, its expression was considered as the mean expression across all the probes.

WGCNA

Immune-related genes (IRGs) were obtained from the InnateDB database (http://www.innatedb.com)¹³ and their expression profiles in the GSE53963 dataset were extracted by WGCNA to identify phenotype-related gene modules. A module eigengene was defined as the first principal component of the module’s expression matrix; module membership (MM), as the association between a gene and its module; and gene significance (GS), as the correlation of a gene with a phenotype. Genes significantly associated with MM (P < 0.05) were defined as hub genes. Clinicopathological features including death (event), tumor grade and overall survival (OS) were considered as phenotypes in WGCNA.

Univariate Cox Regression and LASSO Analysis

Univariate Cox regression was performed to assess the potential relation of hub IRGs with prognosis. P values < 0.05 were considered statistically significant. The LASSO method was then applied using the glmnet package (https://CRAN.R-project.org/package=glmnet) in R for feature selection to identify IRGs with non-zero coefficients. The IRG-based signature to predict prognosis was generated based on the following formula:

where Expr represents the expression value of the gene and Coef the coefficient of the corresponding gene provided by LASSO.

The survivalROC package was used to draw the tROC curves and assess the prognostic value of the IRG-based signature.¹⁴ OSCC patients were divided into high- and low-signature index groups according to the corresponding cut-off, and their Kaplan–Meier survival curves were compared using the log-rank method. Differences with P < 0.05 were considered statistically significant. Univariate/multivariate Cox proportional hazards regression was also applied to compare the relative prognostic value of the signature index with that of routine clinicopathological features.

Immunohistochemical Analysis

The expression of signature IRGs in normal ovarian and ovarian cancer tissues was explored using the Human Protein Atlas database (http://www.proteinatlas.org), which contains information on a wide variety of human protein-coding genes.¹⁵ The expression of a given gene was investigated in control and cancer tissues using the same antibody, and expression level was semi-quantitated as “high” (at least 75% of cells staining positive), “medium” (25–75%) or “low” (no more than 25%).

Gene Set Enrichment Analysis (GSEA)

GSEA was performed using normalized mRNA expression profiles of a JAVA-based software (http://software.broadinstitute.org/gsea/index.jsp)¹⁶ to explore the differences in biological processes (BPs) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways between the high- and low-signature index groups. The gene sets c5.bp.v6.2.symbols.gmt and c2.cp.kegg.v6.2.symbols.gmt from the Molecular Signatures Database (MSigDB) were used as a reference.¹⁷ Nominal P values < 0.05 were considered statistically significant.

Tumor Microenvironment Cell (TMC) Abundance Between High- and Low-Signature Index Groups

A total of 24 types of TMCs were identified for each sample using single-sample GSEA (ssGSEA) based on a previous study.¹⁸ The marker genes of all TMCs were also obtained from previous studies.^19,20 The ssGSEA analysis was performed using the gsva function of the GSVA package in R¹⁶ and the abundance of TMCs was compared between the high- and low-signature index groups in the GSE53963 and TCGA datasets.

Statistical Analysis

Student’s t-test was used to compare continuous variables. Univariate/multivariate Cox proportional hazards analysis and Kaplan–Meier survival analysis were used to compare survival between low- and high-signature index groups. All analyses were conducted in R version 3.6 (https://www.r-project.org/). All tests were two-sided, and values with P < 0.05 were considered statistically significant.

Results

Identification of Survival-Related Gene Modules

To identify the survival-related gene modules of OSCC, the expression profiles of 3517 overlapped genes were extracted by WGCNA (Figure 1A). The power value was set to 3, which was the lowest value with an independence degree of up to 0.95 (Figure 1B). Nine modules were identified (Figure 1C), of which the red module was positively associated with the stage of OSCC, and the yellow with the overall survival (OS) of OSCC patients (Figure 1D and E).

Figure 1 Weighted gene co-expression network analysis. (A) Venn diagram of immune-related and all genes in the GSE53963 dataset. (B) Scale-free fit index for soft-thresholding powers. (C) Dendrogram of all differentially expressed genes, clustered based on different metrics. (D) Heatmap of the correlation between module eigengenes and clinical traits. (E) Correlation and significance of the red and yellow modules with stage of OSCC or overall survival (OS) of OSCC patients, respectively.

Construction of a Prognostic Signature Based on IRGs

A total of 114 hub genes from the red and yellow modules were subjected to univariate Cox regression analysis, identifying 10 IRGs as significant (P < 0.05) (Table 1). Among them, six genes (ARHGEF18, PLEKHA7, MTOR, VPS45, BRCA1, and HINT2) were identified as model characteristic genes (Figure 2A) and their characteristic coefficients were calculated to obtain the signature index. Analysis of tROCs revealed that the area under the curve (AUC) tended to 1 for two- and five-year survival in the GSE53963 training set, suggesting that the signature index is an effective prognostic marker for OSCC (Figure 2B). Patients in the high-signature index group had a significantly lower survival rate than patients in the low-signature index group (Figure 2C), and similar survival results were obtained in the TCGA validation dataset (Figure 2D). Furthermore, analysis of the TCGA dataset based on another optimized cut-off to distinguish between high- and low-signature index groups showed that the signature was an independent predictor of survival after adjustment for routine clinicopathological features (Figure 2E).

Table 1 Immune-Related Genes in the GSE53963 Dataset

Figure 2 Construction of the immune-related prognostic signature for ovarian serous cystadenocarcinoma. (A) Determination of the six model characteristic genes by least absolute shrinkage and selection operator analysis. (B) Time-dependent receiver operating characteristic (tROC) curves to evaluate the ability of the signature index to predict 2-year and 5-year survival in the GSE53963 dataset. (C and D) Survival probability in the GSE53963 and The Cancer Genome Atlas datasets. (E) Forest plot of the signature index. Abbreviation: AUC, area under the tROC curve.

Model Characteristic Genes Highly Expressed at the Protein Level in OSCC

The six identified model characteristic genes were present in The Human Protein Atlas, where they were expressed at higher levels in OSCC samples than in normal ovary samples (Figure 3), strongly supporting our results. The immunohistochemical image credit: Human Protein Atlas. The image available from Human Protein Atlas, the URL of each gene was shown in Supplementary Materials.

Figure 3 Immunohistochemical analysis. Expression of immune-related signature genes in normal ovarian and ovarian cancer tissues explored using the Human Protein Atlas database. Abbreviation: OSCC, ovarian serous cystadenocarcinoma.

Molecular Function in the High- and Low-Signature Index Groups

GSEA indicated that several BPs, including “Golgi organization” and “autophagosome organization” (Figure 4A), as well as several KEGG pathways, including “inositol phosphate metabolism” and “oxidative phosphorylation”, were significantly enriched in the GSE53963 training set (Figure 4B).

Figure 4 Gene set enrichment analysis. (A and C) Biological processes and (B and D) Kyoto Encyclopedia of Genes and Genomes pathways enriched in the (A and B) GSE53963 and (C and D) The Cancer Genome Atlas datasets.

In the TCGA validation set, enriched BPs included “activation of GTPase activity” and “activation of Jun kinase activity” (Figure 4C), while enriched KEGG pathways included “adherens junctions” and the “phosphatidylinositol signaling system” (Figure 4D).

Immune Infiltration Score in the High- and Low-Signature Index Groups

A total of 24 types of immune cells were analyzed: memory and naive B cells; activated and resting dendritic cells; endothelial cells; eosinophils; fibroblasts; M0-, M1-, and M2-type macrophages; activated and resting mast cells; monocytes; neutrophils; activated and resting natural killer (NK) cells; plasma cells; memory activated, memory resting, and naive CD4 T cells; CD8 T cells; follicular helper T cells; gamma delta T cells; and transferred regulatory T cells. In the GSE53963 dataset, abundance of fibroblasts, activated NK cells, and memory-activated CD4 T cells differed significantly between the low- and high-signature index groups (P < 0.05). In the TCGA dataset, in contrast, abundance of naive B cells, endothelial cells, eosinophils, fibroblasts, macrophages M0, and memory resting and naive CD4 T cells differed significantly between the two groups (P < 0.05) (Figure 5).

Figure 5 Abundance of immune cells in ovarian serous cystadenocarcinoma patients in the (A) GSE53963 and (B) The Cancer Genome Atlas datasets. Relative abundance was determined indirectly from the expression level of marker genes corresponding to each type of immune cell.

Discussion

OSCC is a deadly gynecological cancer associated with high rates of recurrence and mortality. Immunotherapy has recently been applied for personalized treatment of serous ovarian neoplasms²¹ based on the association of immune cells in the TME with OSCC.⁷ The combination of poly(ADP-ribose) polymerase (PARP) and immune checkpoint inhibition have also been proposed as a therapeutic approach against ovarian cancer,²² while immune checkpoint blockade therapy, PARP inhibitors, neoantigen vaccines, and personalized T-cell therapy have been applied in the clinic.²³ Several studies have also reported gene expression characteristics for OSCC, though the clinical relevance of disease models based on such expression profiling remains unclear.

Here, in order to improve the prognosis of patients with OSCC, we developed a novel gene-based signature index to identify new immune-related prognostic biomarkers. IRGs were extracted from the InnateDB database, which includes a comprehensive list of genes related to innate immunity,¹³ and their expression profiles in GSE53963 were extracted by WGCNA to identify gene modules associated with OSCC. Subsequent LASSO analysis identified ARHGEF18, PLEKHA7, MTOR, VPS45, BRCA1, and HINT2 as characteristic genes, and their coefficients were used to calculate the signature index.

PLEKHA7 can alter E-cadherin-containing cell–cell contacts, while high PLEKHA7 expression has been associated with good prognosis in epithelial ovarian cancer.²⁴ MTOR and VPS45 have also been related to prognosis in OSCC,^25,26 while ARHGEF18 has been identified as an indicator of genetic susceptibility to non-idiopathic pulmonary arterial hypertension.²⁷ BRCA1 mutations have been found in invasive ovarian cancer patients,²⁸ while we are unaware of studies reporting associations between HINT2 and OSCC. By modulating HINT2 mRNA translation, dynamic N6-methyladenosine (m6A) RNA modification can suppress ocular melanoma.²⁹ Thus, the signature developed in the present study using the six characteristic genes may serve as a new prognostic biomarker for OSCC. Using the median value of this signature as a cut-off, we were able to stratify patients into groups showing different survival.

Our GSEA analysis revealed that oxidative phosphorylation, an active metabolic pathway in serous ovarian neoplasms,³⁰ is significantly enriched in the GSE53963 training set. Moreover, we identified adherens junctions, which are involved in the Mical-L2 regulated epithelial–mesenchymal transition (EMT) in ovarian cancer cells,³¹ as an enriched KEGG pathway in the TCGA validation set. The loss of E-cadherin during the EMT disrupts adherens junctions, thus preventing anoikis and promoting metastasis and chemoresistance.³²

Furthermore, 24 immune cell types were identified in the TME of OSCC by ssGSEA analysis. Among them, dendritic and NK cells, which are an integral part of the innate immune response to prevent the development of OSCC,³³ had higher expression in the low-signature index group than in the high-signature index group. Moreover, the survival rate in the high-signature index group was lower than that in the low-signature index group, consistent with the results of a survival analysis of this study. Taken together, our results suggest that the signature index may serve as an immune marker for OSCC, while also providing a foundation for further mechanistic and therapeutic studies.

Our study was limited to the size of the samples in the public databases from which we obtained our data, and we performed analyses based only on bioinformatics. Our results should be verified and extended in biological experiments, particularly the hypothesis that the genes in our novel signature are involved in immune responses to ovarian cancer.

Conclusion

In summary, we developed an immune-related gene-based signature that may be useful as a prognostic biomarker for OSCC patients. Moreover, the insights from our research may provide leads for further experiments to understand the onset, progression and treatment of OSCC.