A Comprehensive Risk Assessment Model for Ovarian Cancer Patients with Phospho-STAT3 and IL-31 as Immune Infiltration Relevant Genes

Introduction

Ovarian carcinoma is a malignant tumor that has a high mortality rate.¹ With the recent increase in the world’s aging population, the incidence of ovarian carcinoma has also increased;² the mortality rate for ovarian cancer is higher for women over the age of 65 due to the difficulty in early diagnosis and the survival rate is less than 30%.³ Additionally, there is a lack of treatment options available for patients at advanced stages of the disease. To help patients with a poor prognosis get better nursing, an excellent risk assessment model to identify patients with poor prognosis after surgery needs to be constructed. This can contribute to the understanding of ovarian cancer pathogenesis and management of its patients. Another is to find new cures for patients with poor prognosis.

Firstly, clinical features are commonly used as important factors for predicting the patient’s outcome, but sometimes it becomes very difficult to make a credible prediction. Numerous related genes^4–6 and transcripts^7,8 have been identified as cancer biomarkers, that could benefit the malignancy diagnosis. However, a single gene is not a competent model as it ignores the inter-gene associations and pathways. Recently, co-expressed gene modules⁹ and pathway models have also emerged as potential tools to improve the prognosis of ovarian cancer.¹⁰ Models that combine genes with pathways and clinical features are considered to be useful; they can be exploited and analyzed using various bioinformatics methods.^11–13 Secondly, the identification of cancer biomarkers is key, because they usually act as hub genes in tumorigenesis and metastasis. In recent years, immunotherapy has captured the interest of researchers, clinicians, and pharmaceutical companies, particularly due to its promise in treating various forms of cancer.^14,15 The resulting immune infiltration in ovarian cancer tissues is a major prognostic factor.¹⁶ Therefore, there is a biological rationale for the development of immunotherapy in ovarian cancer based on identifying genes related to immune infiltration, which may benefit the patients with poor prognosis.

In this research, a comprehensive ovarian cancer risk assessment model consisting of prognosis-related pathways and clinic factors was studied, using large-scale bioinformatics analysis and based on the Pathway Dysregulation Score (PDS) matrix, Cox-Proportional Hazards (Cox-PH) regression, and L1-least absolute shrinkage and selection operator (L1-LASSO) penalization. We also generated a gene-based assessment model based on prognosis-related gene biomarkers and a pathway-based assessment model on prognosis-related pathways. When tested on the training set and five validation sets using Kaplan–Meier survival analysis and receiver operating characteristic (ROC) curves, the comprehensive risk assessment model with clinical factors and prognosis-related pathways showed a sound prediction than the gene model and the pathway model. Phospho-STAT3 (p-STAT3) and IL-31 were also found to be related to immune infiltration as indicators to patient survival times and could serve as potential therapeutic targets for patients with poor prognosis.

Materials and Methods

Data and Grouping

The ovarian cancer gene expression profiling data were obtained from the databases, the Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO). The clinical characteristics of the training and validation sets are listed in Table S1.

GSE32062,¹⁷ the training set used for establishing the primary model, has 260 expression profiles downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). The platform of GPL6480, Agilent-014850 Whole Human Genome Microarray 4×44K was used for the expression array of GSE32062.

Five independent datasets were used as validation sets. One was obtained from the TCGA database (https://gdc-portal.nci.nih.gov/). It contained 419 mRNA-seq expression profiles of ovarian tumors; Illumina HiSeq 2000 RNA Sequencing was used. The other 4 datasets that were obtained from GEO were GSE49997, GSE8842, GSE26712, and GSE31245, and their platforms were GPL2986 ABI Human Genome Survey Microarray Version 2, GPL5689 Agilent Human 1 cDNA Microarray, GPL96 [HG-U133A] Affymetrix Human Genome U133A Array, and GPL8300 [HG_U95Av2] Affymetrix Human Genome U95 Version 2 Array, respectively.

Data Pre-Processing and Differentially Expressed Genes (DEGs) Selection

For the data that was downloaded from the Affymetrix platform (GSE26712 and GSE31245), the background was corrected, quantiles were normalized and complementarity of missing values was found, which were conducted using the oligo package in R language (version 3.4.1) (http://www.bioconductor.org/packages/release/bioc/html/oligo.html).¹⁸ For GSE32062, GSE49997 and GSE8842, the limma package in R language was utilized (https://bioconductor.org/packages/release/bioc/html/limma.html).¹⁹ Probes were annotated according to the annotation platform, and an approximately normal distribution was yielded using log2 transformation. TCGA data were subjected to quantile normalization using the preprocess core package²⁰ in R language (version 3.4.1) (http://bioconductor.org/packages/release/bioc/html/prepro-cessCore.html).

In the GSE32062 training set, patients with relapse and recurrence within 6 months and death were assigned into the poor prognosis group while those with recurrence after more than 24 months were assigned into well prognosis group. Screening of DEGs was performed between two groups using the limma package with the threshold of FDR < 0.05 and |logFC| > 0.263.

Screening Prognosis-Related Genes and Clinical Factors

Univariate and multivariate Cox regression analyses²¹ were performed to investigate the significant DEGs and the clinical factors of samples using survival package in R language (version 3.4.1) (http://bioconductor.org/packages/survivalr/). A P-value of less than 0.05, obtained by the Log rank test, was set as the significant threshold. A two-way hierarchical clustering analysis based on the centered Pearson correlation algorithm²² was conducted for identifying prognosis-related DEGs and clinical factors using the heatmap package in R (version 3.4.1) (https://bioconductor.org/packages/release/bioc/html/pheatmap.html).²³

Prognostic Gene-Related Risk Pathways Screening

An annotation file (c2.cp.v6.0.symbols.gmt) containing 217 Biocarta and 186 KEGG accessions were downloaded from the Gene Set Enrichment Analysis (GSEA) database. Then, the matrix of prognosis-related gene expression was converted to a PDS matrix based on a principal components analysis (PCA) dimensionality reduction algorithm using the pathifier package²⁴ in R language (version 3.4.1) (http://bioconductor.org/packages/pathifier/). Finally, the identified pathways were used to fit the Cox-PH model based on the L1-LASSO estimation²⁵ for selecting the optimal panel of prognosis-related pathways. The optimal value for penalization coefficient lambda was selected by running cross-validation likelihood (cvl) 1000 times.

Construction and Verification of Risk Assessment Model Based on Pathways

Based on the optimal panel of prognosis-related pathways, a pathway-based risk assessment model was constructed and the Prognosis Index (PI) of each sample was calculated. The median score of PI was used as the cut-off point, and samples in the training set were divided into high and low-risk groups. The survival package was used to draw Kaplan–Meier survival curves for evaluating the relationship between this risk model, and the outcomes of progression-free survival (PFS) and overall survival (OS). ROC curves were drawn for evaluating the classification efficiency of the model. At the same time, the prognostic ability of this risk model was subsequently assessed in validation sets.

Construction and Verification of Risk Assessment Model Based on Gene Sets

Using the matrix of the expression values of prognosis-related genes as input, the Cox-PH model based on L1-LASSO estimation was applied to screen the optimal prognostic genes. A gene-based risk assessment model, similar to the pathway-based risk assessment model, was constructed by using the median PI score as the cut-off point. The Kaplan–Meier survival analysis was performed and ROC curves drawn in the same way.

Construction and Verification of Comprehensive Risk Assessment Model Based on Clinical Risk Factors and Optimal Panel of Prognosis-Related Pathways

A comprehensive risk assessment model was constructed by combining optimized prognosis-related pathways and clinical factors. The PI score of each sample was recalculated based on the gene and the pathway model. The median PI score was selected, and the samples were divided into high and low-risk groups. Then, the association between the comprehensive model and the prognosis was exploited using Kaplan–Meier survival analysis and ROC curves in the training and validation sets, respectively.

Immunohistochemistry and Immune Infiltration Analysis of Ovarian Cancer Tissues

Ten ovarian cancer tissues were collected in Xin Hua Hospital Affiliated to Shanghai Jiao Tong University School of Medicine, Shanghai, China. Five samples were poor prognosis, and 5 were well prognosis. This research was approved by the Ethics Committee of Xin Hua Hospital Affiliated to Shanghai Jiao Tong University School of Medicine. All recruited patients provided a signed, informed consent prior to the collection of ovarian cancer tissues.

Hematoxylin-eosin (HE) staining was performed for evaluating immune infiltration by observing the histomorphology and cell morphology. Among optimized prognosis-related genes and pathways, candidate genes were chosen according to the literature reviewed on PubMed (https://www.ncbi.nlm.nih.gov/pubmed/). Candidate genes that reported of functions related to immunity were briefly selected. Immunohistochemistry was performed for the selected candidate genes, as Zhang described.²⁶ The phospho-STAT3 antibody (#9145S, Cell Signaling Technology, USA) was used at a dilution ratio of 1:200, and the IL-31 antibody (sc-515,415, Santa Cruz Biotechnology, USA) was used at a dilution ratio of 1:100.

Results

Screening of Prognosis Genes and Clinical Factors

All 264 of the DEGs screened from GSE32062 (training set) are shown in Figure 1A. The two-way hierarchical clustering graph based on gene expression values of DEGs in well and poor prognosis groups was drawn (Figure 1B). Using univariate and multivariate Cox regression analysis, 172 prognosis-related genes and one clinical prognostic factor (namely clinical stage) were screened from the DEGs and the clinical information of samples in the training set. The clustering analysis of the expression values of the 172 prognosis-related genes were conducted. Samples were clearly divided into two groups, namely Group 1 and 2 (Figure 1C), which included 109 and 151 ovarian tumor samples, respectively. Group 1 contained 95 recurrent samples and 14 non-recurrent samples, and Group 2 contained 98 recurrent samples and 53 non-recurrent samples (χ² = 15.2484, P = 9.426E–05). Group 1 contained 14 samples of clinic stage III and 95 samples of clinic stage IV, and group 2 contained 42 samples of clinic stage III and 109 samples of clinic stage IV (χ² = 7.5327, P = 0.006059). The association analysis between the two groups and the prognosis are shown in Figure 1D. According to the results, clustering was significantly associated with the recurrence (Log Rank P = 0.0001928). The samples in Group 2 had a longer PFS than that of the samples in Group 1 (29.29 ± 23.88 vs 23.52 ± 21.34, t-test P = 0.0419). Meanwhile, clustering was significantly associated with the OS (Log Rank P = 0.03222). The samples in Group 2 have a longer OS than that of the samples in Group 1 (44.89 ± 22.68 vs 38.71 ± 28.01, t-test P = 0.0496).

Figure 1 The differentially expressed genes (DEGs), prognosis-related genes, and clinical factors screened from the training set. (A) The volcanic map for DEGs. Red triangles and green triangles indicate genes that were upregulated and downregulated in the well prognosis group, respectively. With the threshold of FDR < 0.05 and |logFC| > 0.263. (B) Two-way hierarchical clustering map of DEGs. Top bar in orange color indicates the well prognosis group and blue color indicates poor prognosis group. (C) Two-way hierarchical clustering map based on 172 prognosis-related genes. The first line below the sample tree indicates cluster 1 (black bars) and 2 (white bars) formed by the category. The second line shows the sample recurrence information, and the orange and blue bars represent the recurrent and non-recurrent samples, respectively. The third line shows the clinic stage information, in which the green and red represent the Clinic Stage III and IV samples, respectively. P values were obtained using Chi-square test. (D) The Kaplan–Meier curves of progression-free survival (PFS, left) and overall survival (OS, right) for clusters and clinic stages. Horizontal axis represents the survival period and the vertical axis represents the frequency.

The clinical stage has been identified to be significantly correlated to prognosis (P values < 0.05, Table S2) and the Kaplan–Meier survival curves for the clinical stage and prognosis are shown in Figure 1D.

Screening of Prognosis Gene-Related Risk Pathways

Based on 172 prognosis-related gene expression matrices, 92 ovarian cancer-related pathways (61 KEGG pathways and 31 Biocarta pathways) and path-based PDS matrices were created using pathifier package. Lambda equalled to 0.1753 while cvl reached the max value of −1145.3164 (Figure 2A); the path coefficient is shown in Figure 2B. Finally, 9 optimal combinations of pathways (one Biocarta and 8 KEGG pathways, see Table 1) were obtained using Cox-PH model, which involved 19 prognostic genes.

Table 1 Optimal Prognostic Pathways Screened Using the Cox-PH Model Based on L1-LASSO Estimation

Figure 2 Prognosis analysis for the selected 9 pathways. (A–C) Two-way hierarchical clustering map based on 9 prognosis-related pathways. The first line below the sample tree indicates the group I (black bars) and II (white bars) formed by category. The second line shows the sample recurrence information, and the orange and blue bars represent the recurrent and non-recurrent samples, respectively. The green and red in third line represent the Clinic Stage III and IV samples, respectively. P values were obtained using Chi-square test. (D) The Kaplan–Meier curves of progression-free survival (PFS, left) and overall survival (OS, right) for clusters. Horizontal axis represents the survival period and the vertical axis represents the frequency.

Based on the PDS matrix values of 9 pathways, cluster analysis was performed in the training set of GSE32062 (Figure 2C). As expected, the samples were clearly divided into Group I and Group II, which included 138 and 122 ovarian cancer tumor samples, respectively. There were 117 recurrent samples and 21 non-recurrent samples in Group I, while 76 recurrent samples and 46 non-recurrent samples in Group II (χ² = 15.963, P = 6.459E–05). Group I contained 95 samples of Clinic Stage III and 43 samples of Clinic Stage IV, while Group II contained 109 samples of Clinic Stage III and 13 samples of Clinic Stage IV (χ² = 14.9182, P = 0.0001123). Meanwhile, the clustering of pathways was significantly associated with PFS (Log Rank P = 0.0002581, Figure 2D left) and OS (Log Rank P = 0.02053, Figure 2D right). The samples in Group II had a longer PFS and OS than those of the samples in Group I (31.39 ± 26.58 vs 22.88 ± 18.46, t-test P = 0.00339; 46.24 ± 28.93 vs 43.55 ± 20.95, t-test P = 0.0397).

Prognostic Abilities of Pathway-Based Risk Assessment Model

Using the median PI score as the cut-off point, a pathway-based risk assessment model was constructed. The prognostic ability of the model in a training set and validation sets are shown in Table S3. Kaplan–Meier survival analysis was applied to assess the association between risk groups and PFS as well as OS. As a result, patients in the low-risk group had longer PFS and OS in all sets (all P < 0.05) except for OS in GSE8842 (P > 0.05). Meanwhile, the ROC for PFS and OS based on this pathway-based model were conducted (Supplemental Figure S1).

The Prognostic Ability of Gene-Based Risk Assessment Model

The optimal prognosis-related genes screened by the Cox-PH model based on L1-LASSO estimation are shown in Table 2. The prognostic ability of gene-based risk assessment model in training and validation sets are shown in Table S4. According to the results of Kaplan–Meier survival analysis (Supplemental Figure S2), low-risk patients had a longer PFS than that of high-risk patients in the training set (P = 2.13E–11) as well as validation sets of TCGA (P = 0.03871) and GSE49997 (P = 0.0021). Unfortunately, low-risk patients had longer OS than that in high-risk patients in the training set (47.92 ± 24.41 and 41.71 ± 25.29, P = 2.13E–11); there was no significant correlation with OS prognosis in all validation sets (P > 0.05). Additionally, the ROC for PFS and OS based on this gene-based model were conducted (Supplemental Figure S2).

Table 2 The Optimal Prognosis-Related Genes Screened by Cox-PH Model Based on L1-LASSO Estimation

The Prognostic Ability of Comprehensive Assessment Model Based on Clinical Risk Factors and Pathways

The clinical stage was the only optimal prognosis-related clinical factor screened by the Cox-PH model based on L1-LASSO estimation (P = 0.0176). Kaplan–Meier survival analysis ROC were used to assess the association between risk groups and PFS as well as OS for the training and validation sets based on this comprehensive assessment model (Figure 3). As shown in Table S5, low-risk patients had longer PFS and OS than high-risk patients in training and validation sets of TCGA and GSE4997 (all P < 0.05).

Figure 3 Kaplan–Meier survival analysis and receiver operating characteristic (ROC) curves of comprehensive assessment model of clinical risk factor and pathways. (A) Kaplan–Meier curves of progression-free survival (PFS) for GEO and TCGA clusters. (B) Kaplan–Meier curves of overall survival (OS) for GEO and TCGA clusters. Horizontal axis represents the survival period and the vertical axis represents the frequency. (C) ROC curve of PFS and OS for comprehensive assessment model of clinical risk factor and pathways. The abscissa represents sensitivity and the ordinate represents specificity.

Clinicopathologic Correlation of Immune Infiltration and Expression of Phospho-STAT3 and IL-31 in Ovarian Cancer Tissues

Histological sections from 6 ovarian cancer patients (3 well and 3 poor prognoses) were analyzed with HE staining and immunohistochemistry. Among 9 optimal prognostic pathways and 17 optimal prognostic genes, p-STAT3, BTL8, CD207, IL-31, RELN, LIX1, and POU4F1 were chosen as candidates related to immune infiltration according to the literature review on Pubmed. The candidate genes were chosen because they were related to immunity and have never been reported in correlation to immune infiltration in ovarian cancer.

According to the morphology of HE stained tissue sections and immunohistochemistry results, the position of concentrated immune cells showed brown signals in the cell cytoplasm (Figure 4B), that confirmed the correlation of immune infiltration and expression of p-STAT3 and IL-31. Tissues with poor prognosis are positively (3 strong and 2 moderate) correlated to the expression of p-STAT3, while well prognosis tissues showed negative or weak expression levels. In addition, IL-31 was detected in all samples, and tissues with poor prognosis showed a larger proportion of IL-31 positive regions (Figure 4A).

Figure 4 Hematoxylin-eosin (HE) staining and immunohistochemistry of phospho-STAT3 (p-STAT3) and IL-31 in ovarian cancer tissues. (A) Immunohistochemistry photos showed upregulated p-STAT3 and downregulated IL-31 in ovarian cancer tumor sections with poor prognosis, and brown positions indicate the expression of detected proteins. (B) Morphology of HE and immunohistochemistry results showed p-STAT3 and IL-31 expressed at the position of concentrated immune cells, which were marked by arrows in photos. Scale bars indicate 100 μm.

Discussion

There are 3 types of risk assessment models associated with the prognosis of ovarian cancer in this study: a pathway model, a gene model, and a comprehensive model of clinical risk factors and pathways. Through a series of analyses and verification, the pathway model failed in the OS verification of GSE8842, while the gene model failed in the PFS verification of GSE8842 and the OS verification of all validation sets. Undoubtedly, the pathway model is superior to the gene model. The comprehensive model had a good prognostic role for PFS and OS in the training and all validation sets and hence is the most credible. However, one of the study’s limitations is that there are only two validation sets for the comprehensive model. Research that includes more validation sets based on ovarian cancer expression profiles and clinical studies is required for further confirmation.

Nine signaling pathways have been identified and might play pivotal roles in the prognosis of ovarian cancer. The significantly activated pathway, Janus kinase/signal transducer and activator of transcription (JAK-STAT, coefficient = –0.9615305 and P = 0.007, see Table 1) pathway, contributes to the invasion and metastasis of tumors;²⁷ it plays a crucial role in orchestrating the immune system and can modulate the polarization of T helper cells.²⁸ The dephosphorylation tyrosine residues bound to STATs is an important step in downregulating JAK-STAT signaling. Mechanistically, constitutively activated JAK signaling causes inflammation and subsequently phosphorylates additional targets in the STAT family.²⁹ Inhibiting STAT3 could suppress the invasion and metastasis of ovarian cancer cells;^30,31 a STAT3 inhibitor drug has been considered as a novel therapeutic approach for head and neck squamous cell carcinoma³² and ovarian carcinoma.³³ Based on the mechanism of JAK-STAT signaling, we explored the clinicopathologic correlation of immune infiltration and expression of phosphorylated-STAT3 (p-STAT3). The expression of p-STAT3 was confirmed in immune cell-concentrated microenvironment in ovarian tumors with poor prognosis (Figure 4). In addition, as STAT3 is a latent transcription factor that resides in the cytoplasm until activated,²⁹ it may be a potential target for ovarian cancer treatment.

Interleukin 31 (IL-31) has also been identified to be closely related to the outcome of ovarian cancer (Table 2) and immune infiltration (Figure 4) in this study. IL-31 is mainly produced by activated T-helper 2 cells, which have effects on certain kinds of tumors.³⁴ IL31RA, the receptor of IL-31, is highly expressed in ovarian cancer tumors.³⁵ Signaling through the IL31RA involves immune-pathways like MAPK and JAK-STAT, indicated the possible correlation between IL-31 and the clinical presentation of cancer. Moreover, previous studies also reported that serum level of IL-31 correlated with carcinoma stages, depth of invasion, lymph node, and even distant metastasis.³⁶ Besides, a recent study suggested that IL-31 gene polymorphisms might be associated with the development and progression of epithelial ovarian cancer in the Chinese population.³⁷ However, the roles of IL-31 in the development and progression of human epithelial ovarian cancer are largely unknown. Therefore, future trials focusing on IL-31 gene polymorphisms, IL-31 expression levels, and clinical outcomes are needed, and an animal experiment simulating and observing the immune system is useful for revealing the undiscovered mechanisms and novel pathways.

Other pathways and genes that are found in our models can also be explained. Lysosome, as a membrane-bound organelle,³⁸ has critical functions in regulating cell adhesion dynamics,³⁹ invasion, and metastasis,⁴⁰ and cell death.⁴¹ The lysosome has also been reported as a mediator of drug resistance in cancer,⁴² and abnormal lysosomal trafficking is thought to induce the enhanced exosomal export of cisplatin, leading to drug resistance in human ovarian cancer cells.⁴³ From this view, the lysosome might provide crucial pathways in clinical outcomes of ovarian cancer patients. The tight junction pathway also refers to cell adhesion⁴⁴ and it is involved in ovarian cancer progression and targeted therapy.^45,46 Cell adhesion is a critical process in the migration and invasion of cancer cells⁴⁷ and several studies focusing on inhibiting cancer metastasis have been conducted via targeting focal adhesion pathways.⁴⁸ Hyaluronan synthase (HAS1) and Reelin (RELN) are also two significant DEGs according to our results. Previous reports have consistently revealed that upregulated HAS1 enhances the metastasis of cancer cells.⁴⁹ The overexpression of HAS1 might be related to poor outcomes in ovarian cancer patients as it is expressed higher in the poor prognosis group.⁵⁰ RELN was also reported to be associated with metastasis in several types of cancer,⁵¹ and the knockdown of RELN or its downstream components significantly increase cell motility.⁵²

In this study, we have confirmed that the pathway-based model was more convincing and had greater applicability than the gene-based model, and we have generated a comprehensive risk assessment model consisting of a clinical risk factor and pathways that showed a possible bright foreground. Besides, p-STAT3 and IL-31 showing correlation to immune infiltration of ovarian cancer tissues may be novel potential therapeutic targets for treating ovarian cancer with poor prognosis. For instance, we can explore the curative effect of combination therapy with conventional medicals with inhibitors or monoclonal antibodies of p-STAT3 and IL-31. The strategy for searching superior and applicable assessment models for the study might be appropriate for the prognosis risk assessment and understanding of ovarian cancer pathogenesis. The selected prognostic models from our study might help to improve the nursing of ovarian cancer patients; however, further validations based on a large number of ovarian cancer expression profiles and clinical trials are required. The analysis of clinical samples and animal experiments are also needed to confirm the efficacy and feasibility of p-STAT3 and IL-31.