This article is mentioned in:

Introduction

Ovarian cancer (OC), which frequently occurs in postmenopausal women (1), is a cancer with no apparent symptoms until it reaches an advanced stage. The symptoms include bloating, abdominal swelling, pelvic pain and loss of appetite (2). The three most common OC subtypes are high-grade serous carcinomas, sex cord stromal tumors and germ cell tumors (3), which may metastasize to the peritoneum, liver, lungs or lymph nodes (4). OC is difficult to detect and metastasizes early in disease progression; therefore, patients with OC frequently have a poor prognosis (5). Globally, OC affects 1.2 million women and led to 161,100 mortalities in 2015 (6,7). Thus, obtaining an improved understanding of OC progression and recurrence is of great importance for improving its prognosis.

Previously, genes affecting OC were identified, including p21-activated kinase 4 (Pak4), cyclin E1 (CCNE1), RNA binding motif protein 3 (RBM3), YY1 associated protein 1 (YAP) and prominin-1 (CD133). Pak4 overexpression has been reported to contribute to OC cell migration, invasion and proliferation, thus making it a promising prognostic indicator and therapeutic target (8). CCNE1 amplification has been demonstrated to markedly reduce disease-free survival and overall survival, thus indicating that CCNE1-targeted treatment may benefit patients with OC who have upregulated CCNE1 expression (9,10). In epithelial ovarian cancer, RBM3 expression has been associated with cisplatin sensitivity and correlated with a positive patient prognosis (11). Furthermore, YAP has been associated with cell growth and tumorigenesis, and its coexpression with TEA domain transcription factor 4 serves as a predictor of a poor outcome (12,13). CD133 also serves as a predictor of poor OC patient survival, thus suggesting that it may serve as a biomarker of cancer stem cells during disease (14). However, the further identification of prognostic indicators and their potential uses is required.

In recent years, bioinformatics analysis of expression profile data has been gradually used to examine the pathogenesis of human diseases (15). It is known that the progress of a disease is usually mediated by multiple relevant genes and not by a single gene (16,17). Therefore, the present study was designed to mine subnetwork features and build a model to assess OC recurrence risk. In the present study, OC expression profiles were downloaded from a public database, and differentially expressed genes (DEGs) were analyzed and functionally enriched. Following identification of a functional subnetwork, a random forest classifier was constructed and validated. It was considered that this constructed classifier may provide an improved approach for predicting the prognoses of patients with OC.

Materials and methods

Data source

RNA sequence data from 307 OC samples and their corresponding clinical data (including patient vital status and overall survival) were downloaded from The Cancer Genome Atlas (TCGA; http://cancergenome.nih.gov) database. Based on patient vital status, 180 recurrent samples and 72 disease-free samples were identified. Additionally, the GSE44104 dataset [including 20 recurrent and 40 non-recurrent OC samples; platform, (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA] and the GSE49997 dataset (including 204 OC samples; platform, GPL2986 ABI Human Genome Survey Microarray version 2, Applied Biosystems; Thermo Fisher Scientific, Inc.) were obtained from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo) database and utilized as validation sets.

Data preprocessing and DEG screening

Using a z-score algorithm (18), the expression value of each gene was normalized to a normal distribution (mean=0; variance=1). The samples were subsequently analyzed using the differential expression via distance synthesis (DEDS) algorithm, which may be applied to obtain differential expression levels via the distance synthesis of relevant data (19). This approach was used to screen DEGs in the recurrent samples relative to the disease-free samples.

Functional enrichment analysis

Gene ontology (GO; http://www.geneontology.org) analysis may be used to predict the potential functions of gene products (20). Using GO Term Finder (http://search.cpan.org/dist/GO-TermFinder/) as previously described (21), upregulated and downregulated genes were separately enriched. Functional terms with P<0.05 and an association with at least three genes were selected as significant terms.

Weighted gene coexpression network analysis (WGCNA)

Genes may jointly influence alterations in functional terms through their interactions, with functional consistencies between genes also confirmed by significant expression level correlations. To systematically analyze how DEGs with similar expression profiles co-affect OC prognosis, WGCNA (22) was performed. Based on the obtained weighted gene coexpression levels, DEG coexpression was suggested to be significantly associated with OC prognosis. To further identify the genes that were able to differentiate between patients with OC with different prognoses, the verified DEG protein-protein interaction (PPI) pairs were used to construct a PPI network using BioNet 1.24.1 package (http://www.bioconductor.org/packages/release/bioc/html/BioNet.html).

Subnetwork analysis and classifier construction

The BioNet package (http://bionet.bioapps.biozentrum.uni-wuerzburg.de) provides an extensive framework that enables functional subnetworks to be isolated from biological networks. Using the BioNet package in R 3.1.0 (23), as previously described (24), subnetwork analysis was conducted for the PPI network with the P-value/false discovery rate set to 0.01. With the subnetwork nodes as feature genes, a random forest classifier (25) was constructed. For sample labels (recurrence/no recurrence), true and false positive rates were calculated and combined with leave-one-out cross validation (26). Additionally, a receiver operating characteristic (ROC) curve was constructed (27) to evaluate the classification efficiency of the random forest classifier.

Validation using other independent datasets. To confirm that the subnetwork nodes were able to effectively differentiate between patients with OC with different prognoses, the classification efficiency of the random forest classifier for the validation set (GSE44104) was analyzed and presented using a confusion matrix. Additionally, a Kaplan-Meier (KM) survival analysis (28) was performed and combined with the clinical information belonging to the TCGA dataset.

Results

DEG screening

Following implementation of the DEDS algorithm, a total of 44 upregulated and 117 downregulated genes were identified in the recurrent samples relative to the non-recurrent samples, with more downregulated genes identified compared with upregulated genes. Additionally, a volcano plot was constructed to examine DEG expression distributions (Fig. 1).

Figure 1. Volcano plot illustrating differentially expressed gene expression distributions. Upregulated genes are depicted in red and downregulated in green.

Functional enrichment analysis

Using the GO Term Finder, significant GO terms were enriched for the upregulated and downregulated genes separately. For the upregulated genes, the enriched GO terms were primarily associated with the ‘regulation of synapse assembly’ (P=1.04×10−7), ‘regulation of synapse organization’ (P=5.26×10−7), and ‘regulation of synapse structure or activity’ (P=5.71×10−7; Table IA). The downregulated genes were associated with ‘single-multicellular organism process’ (P=4.11×10−5), ‘multicellular organismal process’ (P=5.14×10−5) and ‘positive regulation of steroid hormone biosynthetic process’ (P=3.83×10−3; Table IB), which included cytochrome P450 family 17 subfamily A member 1 (CYP17A1).

Table I. Significant functional terms enriched for the upregulated and downregulated genes.

Table I.

Significant functional terms enriched for the upregulated and downregulated genes.

A, Upregulated genes
Term	Corrected P-value	Count	Genes
Regulation of synapse assembly	1.04×10−7	5	DKK1, LINGO2, TPBG, SLITRK5, SLITRK3
regulation of synapse organization	5.26×10−7	5	DKK1, LINGO2, TPBG, SLITRK5, SLITRK3
regulation of synapse structure or activity	5.71×10−7	5	DKK1, LINGO2, TPBG, SLITRK5, SLITRK3
synapse assembly	2.08×10−6	5	DKK1, LINGO2, TPBG, SLITRK5, SLITRK3
positive regulation of synapse assembly	6.78×10−6	4	LINGO2, TPBG, SLITRK5, SLITRK3
synapse organization	2.23×10−5	5	DKK1, LINGO2, TPBG, SLITRK5, SLITRK3
regulation of nervous system development	2.72×10−4	6	LRRC7, DKK1, LINGO2, TPBG, SLITRK5, SLITRK3
positive regulation of nervous system development	4.44×10−4	5	LRRC7, LINGO2, TPBG, SLITRK5, SLITRK3
regulation of developmental process	1.131×10−3	8	LRRC7, DKK1, SLITRK3, BVES, LINGO2, MEGF10, SLITRK5, TPBG
regulation of multicellular organismal development	1.96×10−3	7	LRRC7, DKK1, LINGO2, MEGF10, TPBG, SLITRK5, SLITRK3
positive regulation of developmental process	2.64×10−3	6	LRRC7, DKK1, LINGO2, TPBG, SLITRK5, SLITRK3
system development	2.84×10−3	10	LRRC7, DKK1, SLITRK3, WDR69, BVES, C8orf85, LINGO2, MEGF10, TPBG, SLITRK5
regulation of multicellular organismal process	3.63×10−3	8	LRRC7, DKK1, SLITRK3, BVES, LINGO2, MEGF10, SLITRK5, TPBG
single-organism developmental process	5.75×10−3	11	LRRC7, DKK1, SLITRK3, WDR69, BVES, C8orf85, LINGO2, MEGF10, MMP7, TPBG, SLITRK5
developmental process	6.65×10−3	11	LRRC7, DKK1, SLITRK3, WDR69, BVES, C8orf85, LINGO2, MEGF10, MMP7, TPBG, SLITRK5
single-multicellular organism process	8.10×10−3	11	LRRC7, DKK1, SLITRK3, WDR69, BVES, C8orf85, LINGO2, MEGF10, MMP7, TPBG, SLITRK5
positive regulation of multicellular organismal process	8.16×10−3	6	LRRC7, DKK1, LINGO2, TPBG, SLITRK5, SLITRK3
multicellular organism development	9.65×10−3	10	LRRC7, DKK1, SLITRK3, WDR69, BVES, C8orf85, LINGO2, MEGF10, TPBG, SLITRK5
regulation of cellular component biogenesis	9.92×10−3	5	DKK1, LINGO2, TPBG, SLITRK5, SLITRK3
B, Downregulated genes
Terms	Corrected P-value	Count	Genes
single-multicellular organism process	4.11×10−5	23	GJA8, TRIM71, TCF23, FOXL2, SCNN1G, CPNE5, ARX, CACNA2D2, SLC32A1, AQP5, CCR4, SPTB, HOXC5, RD3, FGF9, PDZD7, PROK1, IL33, WNT6, CYP17A1, BMP6, LGI1, COLEC11
multicellular organismal process	5.14×10−5	25	GJA8, TRIM71, TCF23, FOXL2, SCNN1G, CPNE5, ARX, CACNA2D2, SLC32A1, AQP5, CCR4, TAS1R3, SPTB, HOXC5, RD3, FGF9, PDZD7, PROK1, IL33, WNT6, ACCN3, CYP17A1, BMP6, LGI1, COLEC11
eye development	5.32×10−5	7	RD3, GJA8, FGF9, BMP6, AQP5, FOXL2, WNT6
sensory organ development	7.88×10−5	8	RD3, GJA8, FGF9, BMP6, PDZD7, AQP5, FOXL2, WNT6
multicellular organism development	1.86×10−4	20	GJA8, TRIM71, TCF23, FOXL2, CPNE5, ARX, SLC32A1, AQP5, CCR4, SPTB, HOXC5, FGF9, RD3, PDZD7, PROK1, WNT6, CYP17A1, BMP6, LGI1, COLEC11
anatomical structure development	1.87×10−4	21	GJA8, TRIM71, TCF23, FOXL2, CPNE5, ARX, CACNA2D2, SLC32A1, AQP5, CCR4, SPTB, HOXC5, FGF9, RD3, PDZD7, PROK1, WNT6, CYP17A1, BMP6, LGI1, COLEC11
single-organism developmental process	5.26×10−4	21	GJA8, TRIM71, TCF23, FOXL2, CPNE5, ARX, CACNA2D2, SLC32A1, AQP5, CCR4, SPTB, HOXC5, FGF9, RD3, PDZD7, PROK1, WNT6, CYP17A1, BMP6, LGI1, COLEC11
system development	5.74×10−4	18	GJA8, TRIM71, TCF23, FOXL2, CPNE5, ARX, SLC32A1, AQP5, CCR4, SPTB, FGF9, RD3, PDZD7, PROK1, WNT6, CYP17A1, BMP6, LGI1
developmental process	6.76×10−4	21	GJA8, TRIM71, TCF23, FOXL2, CPNE5, ARX, CACNA2D2, SLC32A1, AQP5, CCR4, SPTB, HOXC5, FGF9, RD3, PDZD7, PROK1, WNT6, CYP17A1, BMP6, LGI1, COLEC11
positive regulation of steroid hormone biosynthetic process	3.83×10−3	2	BMP6, CYP17A1
positive regulation of growth	4.05×10−3	5	CACNA2D2, FGF9, LGI1, CPNE5, ARX
nervous system development	4.48×10−3	12	FGF9, PDZD7, TRIM71, WNT6, CPNE5, CYP17A1, ARX, BMP6, LGI1, SLC32A1, CCR4, SPTB
positive regulation of organ growth	7.93×10−3	3	CACNA2D2, FGF9, ARX
camera-type eye development	8.92×10−3	5	RD3, GJA8, AQP5, FOXL2, WNT6
anatomical structure morphogenesis	9.15×10−3	12	FGF9, PROK1, PDZD7, TRIM71, FOXL2, WNT6, CPNE5, ARX, BMP6, LGI1, AQP5, SPTB

WGCNA

WGCNA was performed and turquoise and grey modules were identified within the cluster dendrogram (Fig. 2). Within the turquoise module, 87 genes were identified (14 upregulated and 73 downregulated), while in the grey module, 74 were identified (30 upregulated and 44 downregulated). Correlations between the two modules regarding status/prognosis were also analyzed. The results demonstrated that the two modules had significant correlations with status and prognosis (P<0.05; Fig. 3). Furthermore, sample correlations were further examined via heatmap analysis, and it was indicated that the samples in the turquoise module were more strongly correlated when compared with those in the grey module (Fig. 4).

Figure 2. Weighted gene coexpression network analysis cluster dendrogram. Grey means grey module identified in WGCNA with 30 upregulated and 44 downregulated; turquoise means turquoise module identified in WGCNA with 14 upregulated and 73 downregulated genes. WGCNA, Weighted gene coexpression network analysis.

Figure 3. Module-trait relationship graph illustrating that turquoise and grey modules have significant correlations with status/prognosis. Red color indicates a positive correlation, while green means a negative correlation. Grey means grey module identified in WGCNA with 30 upregulated and 44 downregulated; turquoise means turquoise module identified in WGCNA with 14 upregulated and 73 downregulated genes.

Figure 4. Heatmap displaying gene correlations for the turquoise and grey modules. A deeper color indicates a stronger correlation, while a lighter color suggests a weaker correlation. Grey means grey module identified in WGCNA; turquoise means turquoise module identified in WGCNA. WGCNA, Weighted gene coexpression network analysis.

Subnetwork analysis and classifier construction

Following building of the PPI network, subnetwork analysis was performed and a significant subnetwork was identified (Fig. 5). The importance scores for the 16 subnetwork nodes [transcription factor GATA-4 (GATA4); fibroblast growth factor 9 (FGF9); aromatase (CYP19A1); 3β-hydroxysteroid dehydrogenase/δ5-4-isomerase type 2 (HSD3B2); corticosteroid 11β-dehydrogenase isozyme 1 (HSD11B1); CYP17A1; pituitary homeobox 2 (PITX2); left-right determination factor 1 (LEFTY1); homeobox protein ARX (ARX); estrogen receptor β (ESR2); steroidogenic factor 1 (NR5A1); forkhead box protein L2 (FOXL2); myocardin (MYOCD); steroidogenic acute regulatory protein mitochondrial (STAR); vesicular inhibitory amino acid transporter (SLC32A1); and twist-related protein 1 (TWIST1)] are listed in Table II. There were multiple interactions among these subnetwork nodes, including HSD3B2-NR5A1, HSD11B1-HSD3B2, CYP17A1-GATA4, ARX-FOXL2, MYOCD-GATA4, STAR-FGF9 and SLC32A1-PITX2. The subnetwork nodes were taken as feature genes and a random forest classifier was constructed. The generated ROC curve demonstrated that the true and false positive rates separately were 92 and 23% when classifying the recurrent and non-recurrent OC samples (Fig. 6A).

Figure 5. Significant subnetwork identified within the protein-protein interaction network. Upregulated genes are denoted in red and downregulated in green. GATA4, transcription factor GATA-4; FGF9, fibroblast growth factor 9; CYP19A1, aromatase; HSD3B2, 3β-hydroxysteroid dehydrogenase/δ5-4-isomerase type 2; HSD11B1, corticosteroid 11β-dehydrogenase isozyme 1; CYP17A1, cytochrome P450 family 17 subfamily A member 1; PITX2, pituitary homeobox 2; LEFTY1, left-right determination factor 1; ARX, homeobox protein ARX; ESR2, estrogen receptor β; NR5A1, steroidogenic factor 1; FOXL2, forkhead box protein L2; MYOCD, myocardin; STAR, steroidogenic acute regulatory protein mitochondrial; SLC32A1, vesicular inhibitory amino acid transporter; TWIST1, twist-related protein 1.

Figure 6. ROC curve analysis. ROC curves for the (A) Cancer Genome Atlas and (B) GSE44104 datasets. FPR, false positive rate; TPR, true positive rate; AUC, area under the curve; ROC, receiver operating characteristic.

Table II. Importance scores of the subnetwork nodes.

Table II.

Importance scores of the subnetwork nodes.

Node	Score
ARX	15.203
CYP17A1	10.956
CYP19A1	10.902
ESR2	10.537
FGF9	10.148
FOXL2	9.827
GATA4	9.786
HSD11B1	9.223
HSD3B2	9.196
LEFTY1	8.259
MYOCD	7.867
NR5A1	7.634
SLC32A1	6.367
STAR	6.167
PITX2	5.094
TWIST1	4.708

[i] GATA4, transcription factor GATA-4; FGF9, fibroblast growth factor 9; CYP19A1, aromatase; HSD3B2, 3β-hydroxysteroid dehydrogenase/δ5-4-isomerase type 2; HSD11B1, corticosteroid 11β-dehydrogenase isozyme 1; CYP17A1, cytochrome P450 family 17 subfamily A member 1; PITX2, pituitary homeobox 2; LEFTY1, left-right determination factor 1; ARX, homeobox protein ARX; ESR2, estrogen receptor β; NR5A1, steroidogenic factor 1; FOXL2, forkhead box protein L2; MYOCD, myocardin; STAR, steroidogenic acute regulatory protein mitochondrial; SLC32A1, vesicular inhibitory amino acid transporter; TWIST1, twist-related protein 1.

Validation using other independent datasets

The random forest classifier was used to differentiate between samples in the GSE44104 validation set, and the prediction accuracies for the non-recurrence and recurrence groups were 87.5 and 85%, respectively (Fig. 6B). These findings indicated that the subnetwork nodes were important in predicting OC prognosis.

Of the 307 OC samples in the TCGA dataset, only 262 samples remained upon removal of samples without follow-up or survival time information. These 262 samples were divided into high- and low-risk groups using the random forest classifier, and further examined using KM survival analysis. The results demonstrated that the low risk group had a significantly longer survival time compared with the high-risk group (P=0.0166; Fig. 7A). Subsequently, the samples from the second validation set (GSE49997) were divided into high- and low-risk groups using the random forest classifier. Following analysis using a KM survival curve, a significant difference was noted in survival time between the low- and high-risk groups (P=0.0165; Fig. 7B). These findings suggested that the random forest classifier had portability and repeatability.

Figure 7. KM survival analysis. KM survival curves for the (A) Cancer Genome Atlas and (B) GSE49997 datasets. Predicted high risk (cluster 1) groups are indicated in red and low risk (cluster 2) groups in black. KM, Kaplan-Meier.

Discussion

In the present study, 44 upregulated and 117 downregulated genes were identified in the recurrent samples relative to the non-recurrent samples. When performing WGCNA, turquoise and grey modules were identified that had significant correlations with status and prognosis. Furthermore, a significant subnetwork was identified from the PPI network, with the subnetwork nodes (including GATA4, FGF9, CYP19A1, HSD3B2, HSD11B1, CYP17A1, PITX2, LEFTY1, ARX, ESR2, NR5A1, FOXL2, MYOCD, STAR, SLC32A1 and TWIST1) being utilized as feature genes for constructing a random forest classifier. Moreover, the classification efficiency of the random forest classifier was validated and confirmed.

Of the identified nodes, previous studies have reported that GATA4 overexpression, in conjunction with human epidermal growth factor receptor 2, may predict a shorter disease-free survival time and may be utilized as a prognostic marker in patients with ovarian granulosa cell tumor to optimize follow-up management in the early stages (29,30). In OC, GATA4 and transcription factor GATA-6 expression is frequently lost, with this feature specifying the histological subtype prior to tumorigenic transition of the ovarian surface epithelium (31). The upregulation of FGF9 has been detected in primary ovarian endometrioid adenocarcinomas carrying a defective Wnt/β-catenin pathway and serves an important role in promoting the cancer phenotype (32,33). The aromatase enzyme encoded by the CYP19A1 gene acts in the conversion of androgen to estrogen, with CYP19A1 variants potentially able to influence OC susceptibility (34). Thus, GATA4, FGF9 and CYP19A1 appear to be implicated in the mechanisms of OC pathogenesis. Enrichment analysis demonstrated that CYP17A1, which is involved in ‘positive regulation of steroid hormone biosynthetic processes’, was downregulated. A previous study demonstrated that steroid hormones serve a role in OC pathogenesis, and their receptors influence OC patient survival (35). Therefore, CYP17A1 may affect OC patient prognoses through the positive regulation of steroid hormones.

The overexpression of PITX2 has been implicated in OC progression by facilitating cell growth, migration and invasion, and thus may potentially serve as a therapeutic target for patients with high-grade OC (36,37). PITX2 has also been demonstrated to promote OC cell proliferation through the Wnt pathway, which is closely associated with ovarian development and OC (38). In ovarian clear cell carcinoma, the overexpression of LEFTY, a transforming growth factor-β superfamily member, exhibits an anti-tumor effect by affecting cell proliferation and cellular susceptibility to apoptotic signals (39). Estrogen receptor β, which is encoded by the ESR2 gene, has been suggested to serve as a critical factor during OC carcinogenesis (40). Furthermore, within the ESR2 promoter region, the genotypic and allelic frequencies of the single nucleotide polymorphism (SNP) rs3020449 have been demonstrated to exhibit significant differences based on OC stage, thus indicating that SNP rs3020449 may be associated with OC progression (41). These findings suggest that PITX2, LEFTY1 and ESR2 may also serve roles in OC pathogenesis.

NR5A1 serves an important role in ovarian function and development, with an NR5A1 mutation reported to induce 46 XY disorders of sex development (42). FOXL2 is critical for GC (granulosa cell) differentiation during the process of folliculogenesis, with its downregulation potentially serving as an ovarian granulosa cells tumor prognostic factor (43). TWIST1, which may induce epithelial-mesenchymal transition and contribute to tumor metastasis, has been demonstrated to be associated with poor survival in patients with cancer (44,45). Thus, NR5A1, FOXL2 and TWIST1 may be associated with patient survival in OC. Moreover, the multiple interactions within the PPI subnetwork (including HSD3B2-NR5A1, HSD11B1-HSD3B2, CYP17A1-GATA4, ARX-FOXL2, MYOCD-GATA4, STAR-FGF9 and SLC32A1-PITX2) indicate that HSD3B2, HSD11B1, CYP17A1, ARX, MYOCD, STAR and SLC32A1 may also function in OC by interacting with other genes.

However, the present study has several limitations to note. First, the datasets used in the present study had sample size differences, platform differences and data heterogeneities that may affect the prediction accuracy of the random forest classifier. Second, the smaller patient numbers and analytical methods may limit the predictive capability of the present model. Finally, only bioinformatics analyses were conducted in the present study, and no direct experimental validation was performed. Therefore, further analyses are required to validate the obtained results.

In conclusion, 44 upregulated and 117 downregulated genes associated with OC recurrence were identified. Furthermore, the 16 subnetwork node genes that were identified may be critical molecules associated with OC recurrence.

Acknowledgements

Not applicable.

Funding

The present study was supported by grants from the Hubei Province's Outstanding Medical Academic Leader Program and the Hubei Province Health and Family Planning Scientific Research Project (grant no. WJ2015MA024) and the Natural Science Foundation of Hubei Province (grant no. 2017CFB335).

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

LC and LL performed data analyses and wrote the manuscript. LW, XL and HX contributed significantly in data analyses and manuscript revision. JZ conceived and designed the study. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References