Identification of an energy metabolism‑related gene signature in ovarian cancer prognosis.

Changes in energy metabolism may be potential biomarkers and therapeutic targets for cancer as they frequently occur within cancer cells. However, basic cancer research has failed to reach a consistent conclusion on the function(s) of mitochondria in energy metabolism. The significance of energy metabolism in the prognosis of ovarian cancer remains unclear; thus, there remains an urgent need to systematically analyze the characteristics and clinical value of energy metabolism in ovarian cancer. Based on gene expression patterns, the present study aimed to analyze energy metabolism‑associated characteristics to evaluate the prognosis of patients with ovarian cancer. A total of 39 energy metabolism‑related genes significantly associated with prognosis were obtained, and three molecular subtypes were identified by nonnegative matrix factorization clustering, among which the C1 subtype was associated with poor clinical outcomes of ovarian cancer. The immune response was enhanced in the tumor microenvironment. A total of 888 differentially expressed genes were identified in C1 compared with the other subtypes, and the results of the pathway enrichment analysis demonstrated that they were enriched in the 'PI3K‑Akt signaling pathway', 'cAMP signaling pathway', 'ECM‑receptor interaction' and other pathways associated with the development and progression of tumors. Finally, eight characteristic genes (tolloid‑like 1 gene, type XVI collagen, prostaglandin F2α, cartilage intermediate layer protein 2, kinesin family member 26b, interferon inducible protein 27, growth arrest‑specific gene 1 and chemokine receptor 7) were obtained through LASSO feature selection; and a number of them have been demonstrated to be associated with ovarian cancer progression. In addition, Cox regression analysis was performed to establish an 8‑gene signature, which was determined to be an independent prognostic factor for patients with ovarian cancer and could stratify sample risk in the training, test and external validation datasets (P<0.01; AUC >0.8). Gene Set Enrichment Analysis results revealed that the 8‑gene signature was involved in important biological processes and pathways of ovarian cancer. In conclusion, the present study established an 8‑gene signature associated with metabolic genes, which may provide new insights into the effects of energy metabolism on ovarian cancer. The 8‑gene signature may serve as an independent prognostic factor for ovarian cancer patients.

Introduction

Ovarian cancer is a frequently diagnosed gynecologic cancer with a high mortality rate; the estimated annual incidence of this disease worldwide is >200,000 individuals, with ~125,000 deaths (1,2), and as it is difficult to detect at an early stage, the majority of patients are diagnosed at an advanced stage (3). Although advances in chemotherapy, radiotherapy, surgery and targeted therapy have achieved progress in the treatment of ovarian cancer (4), the 5-year OS rate of patients with advanced disease is ~30% (5,6). Clinicopathological features such as histological grade and the International Federation of Gynecology and Obstetrics (FIGO) staging system are widely used prognostic indicators for ovarian cancer (7); however, they are not effective in predicting the survival and chemotherapy response in patients with ovarian cancer (8,9). Currently, no fully validated and clinically applicable tests are available for guiding ovarian cancer treatment decisions. Previous studies have demonstrated that the progression of ovarian cancer is associated with a variety of pathways involved in energy metabolism, including galactose metabolism, which is associated with the risk of developing ovarian cancer (10,11), and that adipocytes promote ovarian cancer metastasis and provide energy for rapid tumor growth (12). Therefore, studying local energy metabolism status and its prognostic value for patients with ovarian cancer may help improve the prediction of clinical outcomes of ovarian cancer and provide reference for personalized medical treatment.

The early observations by Dr Otto Warburg revealed that fundamental metabolic differences existed between malignant tumor cells and adjacent normal cells (13). Other studies also reported the association between altered cellular metabolism and therapeutic outcomes (14–17). Energy metabolism reprogramming, which is a hallmark of cancer, allows tumor cells to produce ATP to maintain the reduction-oxidation balance and biomolecular synthesis required for cell growth, proliferation and migration (18). The metabolic phenotype of cancer cells is heterogeneous; a number of tumor cells are mainly dependent on glycolysis, whereas others exhibit a metabolic phenotype characterized by oxidative phosphorylation (OXPHOS) (19,20). Increasing evidence has demonstrated that glycolysis and oxidation have a symbiotic metabolic relationship in tumor cells; for example, lactic acid and pyruvate produced by glycolysis can be transferred to and used as a substrate for tricarboxylic acid (TCA) intermediates and ATP production in adjacent cancer cells (21). Similarly, malignant tumor cells utilize free fatty acids and ketones released by neighboring catabolic cells for energy production (12,22), which can also be realized by glutamine metabolization through TCA cycle (23). Mitochondrial OXPHOS, driven by glutamine, is a major source for ATP synthesis under hypoxic conditions (24). These observations suggest that targeting particular metabolic pathways in cancer may be an effective strategy for cancer therapy. Therefore, an in-depth understanding of energy metabolism in tumors may contribute to the development of new therapies.

Previous studies have demonstrated that metabolic abnormalities lead to different prognosis of patients, and metabolism-related genes can be used as prognostic markers of tumors. For example, Wu et al (25) have demonstrated that lipid metabolism-related genes can predict the prognosis of patients with glioma. Zhou et al (26) identified a 29 energy metabolism-related gene signature, including interleukin-4, carbohydrate sulfotransferases and branched chain amino acid transaminase 1 (BCAT1), to evaluate the prognosis of diffuse glioma. Genes related to amino acid metabolism such as BCAT2, glutamate-cysteine ligase catalytic subunit and aminoadipate aminotransferase can also predict the prognosis of glioma (27). Ma et al (28) have reported that metabolic deregulations mediate the dedifferentiation of papillary thyroid carcinoma and developed a metabolic gene signature that may be used as a biomarker for dedifferentiated thyroid cancer. Disorders in the metabolic pathway of sputum may affect the progression of breast cancer (29). Liu et al (30) developed a signature of four metabolic genes to predict the overall survival (OS) of patients with liver cancer. However, the expression patterns of metabolism-related genes in ovarian cancer are still unclear, and it is thus necessary to study metabolism-related gene characteristics in ovarian cancer.

The aim of the present study was to identify ovarian cancer molecular subtypes based on energy metabolism-related genes and gene signatures of energy metabolism markers to improve the current understanding of the molecular mechanisms in ovarian cancer energy metabolism and clinical prognosis.

Materials and methods

Data collection and processing

The latest clinical follow-up information of 587 ovarian cancer cases and RNA-seq data from 379 cases were downloaded from The Cancer Genome Atlas (TCGA; http://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga). Genomic Data Commons Application Programming Interface was used to retrieve the data on 29 April 2019. The follow-up information and RNA-seq samples were matched, and 362 cases were selected as they were followed up for >30 days. The samples were randomly divided into two groups (ratio, 3:1), one of which served as the training set (N=271), whereas the other served as the test set (N=91). Similarly, the Affymetrix Human Genome U133 Plus 2.0 Array (http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus) was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), from which the GSE26193 (31) gene expression profile was obtained. A total of 107 samples with clinicopathological characteristics were selected for the test set, and the sample information of each group is presented in Table I. A total of 11 human metabolism-associated pathways were obtained from the Reactome database (https://reactome.org/) (32) (Table II), and a total of 594 genes involved in energy metabolism were identified from these pathways [Tables SI, SII (https://github.com/biocn/OD_data/blob/master/Table%20S2.docx) and SIII]. For the RNA-seq data, >50% genes with an expression level of 0 in each sample were removed. For the chip data, normalization of microarray data was obtained from Affymetrix platform using Robust Multi-Array Average method (33). Probes were matched to genes, and those matched to multiple genes were removed; multiple probes were matched to the median of a gene to obtain a gene expression profile. Workflow is presented in Fig. 1. Ovarian cancer specimens included in TCGA datasets were surgically resected prior to systemic treatment. Samples selected for TCGA analysis had >70% tumor cell nuclei and <20% necrosis. Data on normal ovary tissues from healthy subjects were obtained from the GTEx database (https://www.gtexportal.org/home/index.html).

Table I.

Clinical information of each data set.

Characteristic	TCGA training datasets (n=271)	TCGA validation datasets (n=91)	GSE26193 (n=107)
Age, years
≤60	149	50	NA
>60	122	41	NA
Survival status
Alive	110	31	31
Dead	161	60	76
Tumor stage
I	1	0	21
II	17	3	10
III	207	77	59
IV	43	11	17
Grade
G1	0	1	7
G2	34	8	33
G3	229	80	67
G4	1	0	0

TCGA, The Cancer Genome Atlas; NA, no data available.

Table II.

Pathways involved in energy metabolism in the Reactome database.

Metabolic pathway	Pathway ID	Gene count
Biological oxidations	R-HSA-211859	221
Metabolism of carbohydrates	R-HSA-71387	292
Mitochondrial Fatty Acid Beta-Oxidation	R-HSA-77289	38
Glycogen synthesis	R-HSA-3322077	16
Glycogen metabolism	R-HSA-8982491	27
Glucose metabolism	R-HSA-70326	92
Glycogen breakdown (glycogenolysis)	R-HSA-70221	15
Glycolysis	R-HSA-70171	72
Pyruvate metabolism	R-HSA-70268	31
Pyruvate metabolism and Citric Acid (TCA) cycle	R-HSA-71406	55
Citric acid cycle (TCA cycle)	R-HSA-71403	22
Sum	881 (unique, 594)

Figure 1.

Flow chart of the present study.

Univariate Cox proportional hazard regression analysis

As previously described (34), univariate Cox proportional hazard analysis was conducted to determine the impact of each energy metabolism gene in order to select genes significantly related to patient OS in the training data set. P<0.05 was selected as the threshold.

Identification of molecular types associated with metabolic genes

Nonnegative matrix factorization (NMF) is an unsupervised clustering method widely used in identifying genomics-based tumor molecular subtypes (35,36). To further determine the association between expression levels of energy metabolism genes and phenotypes, the NMF method was used to cluster samples according to the expression profiles of energy metabolism-related genes associated with ovarian cancer prognosis, and standard ‘brunet’ for 50 iterations was selected by NMF. The number of clusters ‘k’ was set between 2 and 10, the average profile width of the common member matrix was calculated using the R package ‘NMF’ (37), and the minimum member of each subclass was set to 10.

Association between molecular subtypes and tumor microenvironment

TIMER is a network resource for systematically assessing the clinical impact of different immune cells in different cancer types (38), and was used to estimate the abundance of six immune cell types, including B, CD4T and CD8T cells, neutrophils, macrophages and dendritic cells, in the tumor microenvironment. The abundance of immune cells in the tumor microenvironment was analyzed in the three different molecular subtypes.

Analysis of genetic differences in molecular subtypes

To detect biological function differences among the three molecular subtypes, differentially expressed genes among the molecular subtypes were analyzed by DESeq2 (39), and the threshold was set to false discovery rate (FDR) <0.05 and |log2(fold-change)|>1. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) enrichment analysis of altered genes categorized by biological process were conducted using the R software packages ‘GSVA’ (40) and ‘clusterProfiler’ (41). The functions of differentially expressed genes were analyzed using the threshold P<0.05.

Screening of robust energy metabolism-related prognostic feature genes

LASSO is a regression modeling with a large number of potential prognostic features; it can perform automatic feature selection that results in signatures with generally effective performance for predicting prognosis (42). The LASSO method has been used in combination with the Cox model for survival analysis, and has been successfully applied for building sparse signatures for survival prognosis in a range of application areas, including oncology (43–45). In the current study, genes with differential expression among the three molecular subtypes were selected for univariate survival analysis, and prognostic genes were screened using the threshold value of 0.05. Robust prognostic characteristic genes were screened using the R software package ‘glmnet’ (46), and the optimal characteristics were evaluated by 10-fold cross validation.

Construction of energy metabolism-related prognostic gene signatures

To obtain a robust prognostic gene, multivariate Cox regression analysis was performed, and the following risk scoring model was constructed:

Where n is the number of prognostic genes, Expk is the expression value of a prognostic gene, and eHRk is the estimated regression coefficient of a gene in the multivariate Cox regression analysis. Specific genes were selected by LASSO analysis, and the 8-gene signature was constructed using the calculated risk scores.

Correlation analysis between gene signatures and KEGG pathways

Gene Set Enrichment Analysis (GSEA; http://software.broadinstitute.org/gsea/downloads.jsp) (47) was performed using the MSigDB 6.2 (48). The C2 Canonical pathway gene set collection containing 1,320 gene sets was selected for analysis. Gene sets with FDR <0.05 after 1,000 permutations were considered as significantly enriched.

Comparison with existing prognostic features of ovarian cancer

To assess the survival classification and predictive power of the 8-gene signature, four published prognostic features were retrospectively reviewed, including a 3-gene (49), 6-gene (50), 8-gene (51) and 5-gene signature (52). To ensure that the models were comparable, according to the corresponding genes in the four models, the risk scores of each ovarian cancer sample in TCGA were calculated using the same method, the ROC of each model was evaluated, and the samples were divided according to the median risk score. OS differences between the high- and low-risk groups were calculated.

Statistical analysis

Kaplan-Meier (KM) curves were plotted when Youden's index in each data set was used as a cutoff to compare the survival risk between the high- and low-risk groups, and the data were analyzed by the log-rank test followed by Bonferroni correction. Multivariate Cox regression analysis was performed to evaluate whether gene markers were independent prognostic factors. Statistical analysis between multiple groups was performed by ANOVA, followed by Dunnett's test. ROC analysis was performed using R package ‘pROC’ (https://cran.r-project.org/web/packages/pROC/index.html). The χ2 test was used for mutation frequency detection. All statistical analyses were processed in R 3.4.3 (https://mirrors.tuna.tsinghua.edu.cn/CRAN/) with default software parameters, unless otherwise stated. P<0.05 was considered to indicate a statistically significant difference.

Results

Identification and molecular classification of genes involved in energy metabolism

For the TCGA training set, univariate regression analysis was performed to establish an association between patient OS and energy metabolism-related gene expression. A total of 39 energy metabolism-related genes with prognostic significance (Table SIV) were identified. With the minimum member of each subclass set to 10 and the optimal cluster number set to 3, the average profile width of the common member matrix was determined by the R package ‘NMF’ according to the cophenetic, dispersion and silhouette indicators (Fig. 2A). The results exhibited different expression patterns with prognostic differences in energy metabolism-related gene expression profiles in the three subtypes (Fig. 2B). In addition, significant differences in OS time among the three subtypes were detected (Figs. 2C and S1); of note, the C1 group was associated with the worst prognosis. Further analysis was conducted on the association among gene mutations in the three subtypes, and 20 genes with the highest mutation rate in each subtype were selected to obtain a total of 39 genes. The genes with the highest mutation rates among the three subtypes were not identical (Figs. 2D and S2). The mutation frequencies of transformation/transcription domain-associated protein (TRRAP), zinc finger protein 551, myosin heavy chain 13 (MYH13) and EvC ciliary complex subunit 2 in C1 were significantly higher compared with those in C2 and C3 (χ2 test P<0.05), the mutation frequency of the von Willerbrand factor in C2 was higher compared with that in the other two groups (P<0.05), and the mutation frequency of filamin B in C3 was higher compared with that in the other two groups (χ2 test p<0.05).

Figure 2.

(A) Consensus map of NMF clustering. (B) Heat map of gene metabolism-related gene expression levels in different molecular subtypes. (C) Overall survival curves of the three molecular subtypes. (D) Heatmaps of mutation frequencies of the top 20 genes with the highest number of mutations in each subtype in each sample. (E) Heat map of KEGG pathway scores for three molecular subtypes. (F) Comparison with published subtypes; different colors represent different subtypes, and the ordinate represents the percentage of samples. NMF, nonnegative matrix factorization; KEGG. Kyoto Encyclopedia of Genes and Genomes.

Further analysis by the R software package ‘GSVA’ revealed that 102, 60 and 21 significantly different KEGG pathways were present in C1 vs. C2/C3, C2 vs. C1/C3 and C3 vs. C1/C2, respectively. In the C1 group, the scores of pathways associated with tumorigenesis and tumor development, such as the ‘TGF BETA SIGNALING PATHWAY’ and ‘ECM RECEPTOR INTERACTION’, were significantly higher compared those of the other two groups. In the C2 group, the scores of major diseases such as ‘ALZHEIMERS DISEASE’ and ‘PARKINSONS DISEASE’ were significantly higher compared with those of the other two groups, whereas the scores of pathways such as ‘PATHWAYS IN CANCER’ and ‘PROSTATE CANCER’ were significantly lower compared with those of the other two groups. In the C3 group, the overall pathway score was low. These result suggested that the C1 subtype may be associated with a poor prognosis (Fig. 2E). Comparison between the subtypes identified in the preset study and previously published molecular subtypes of ovarian cancer (53) demonstrated that the C1 subtype was comparable with the published mesenchymal subtype (78.26% match), which is also associated with a poor prognosis (Fig. 2F).

Microenvironmental characteristics and metabolic pathway differences of tumor molecular subtypes

No significant differences were detected after comparing the clinical features such as stage, grade and age in the three subtypes (Table III), suggesting that independent factors may be affecting the different clinical outcomes of the three subtypes. The tumor microenvironments were analyzed, and the immune cell contents of the three subtypes were compared; the B-cell score was significantly lower, whereas CD4 and CD8 cell, neutrophil, macrophage and dendritic cell scores were higher in the C1 subtype compared with those in the C2 and C3 subtypes (Fig. 3A), which suggested that the immune cells of the tumor microenvironment in the C1 subtype was more active; in addition, a gene expression matrix was extracted based on the gene set in the 11 energy metabolism-related pathways, the values of which were defined by the mean value of the corresponding gene in each sample, and the expression levels of the genes in these pathways in the different subtypes were determined. The results demonstrated that the expression of genes in ‘Biological oxidations’ in C3 was significantly higher compared with that in C1 and C2, whereas the expression of genes in ‘Metabolism of carbohydrates’ was higher in C1 compared with that in C2 and C3. The expression levels of genes in the pathways ‘Mitochondrial fatty acid beta-oxidation’, ‘Pyruvate metabolism’, ‘Citric acid cycle (TCA cycle)’ and ‘Pyruvate metabolism and Citric Acid (TCA) cycle’ was significantly lower in C1 compared with those in C2 and C3, whereas the expression of genes in ‘Glycogen synthesis’ and ‘Glycogen metabolism’ in C2 were higher compared with those in C1 and C3; in addition, the expression of genes in ‘Glycogen breakdown (glycogenolysis)’ in C2 was lower compared with that in C1 and C3, and no significant differences in the genes in ‘Glucose metabolism’ and ‘Glycolysis’ were observed among the three subtypes (Fig. 3B). These results suggested that most of the energy metabolism pathways in the three subtypes were different, and these subtypes exhibit different metabolic patterns.

Table III.

Clinical information statistics of three molecular subtypes.

Clinical features	C1	C2	C3
Status
Alive	17	46	78
Dead	45	64	112
Stage
I	0	0	1
II	0	8	12
III	51	80	153
IV	10	21	23
Grade
G1	1	0	0
G2	5	10	27
G3	56	96	157
G4	0	1	4
Age, years
≤60	31	66	102
>60	31	44	88

Figure 3.

(A) Immune cell scores among the three molecular subtypes. (B) Differences in expression levels of 11 energy metabolism pathways among the three molecular subtypes. TPM, transcripts per million.

Identification of differentially expressed genes in the C1 subtype

To analyze the molecular mechanism underlying the poor prognosis of patients with the C1 subtype, the gene expression differences between C1 and C2 or C3 were evaluated. Compared with C2, 342 genes were upregulated and 331 were downregulated in C1 (Fig. 4A); similarly, compared with C3, 316 genes were upregulated and 257 were downregulated in C1 (Fig. 4B). The differentially expressed genes contained a total of 888 genes (Table SV), of which 359 were shared between C2 and C3 (Fig. 4C). GO biological process and KEGG functional enrichment analysis was performed on the 888 DEGs, which identified 28 enriched KEGG pathways (Table SVI), suggesting that the 28 shared pathways may be involved cancer progression; the enriched KEGG pathways included the ‘PI3K-Akt signaling pathway’, ‘cAMP signaling pathway’ and ‘ECM-receptor interaction’ (Fig. 4E). In addition, 515 GO terms were enriched in the C1 subtype in the biological process category (Table SVII), including ‘angiogenesis’ and ‘regulation of trans-synaptic signaling’ (Fig. 4D), which were also associated with the development of cancer.

Figure 4.

(A) Volcano plot of differentially expressed genes between the C1 and C2 subtypes. (B) Volcano plot of differentially expressed genes between the C1 and C3 subtypes. (C) The intersection between differentially expressed genes. (D) Top 20 enriched Gene Ontology terms in the ‘biological process’ category. (E) Top 20 Kyoto Encyclopedia of Genes and Genomes biological pathways. Colors indicate saliency; size indicates the number of genes in the pathway. FDR, false discovery rate.

Identification of 8-gene signature for ovarian cancer survival

The association between the differentially expressed genes in C1 and prognosis was analyzed, and 82 significant prognostic factors were selected as candidate genes (Table SVIII). Dimensional-reduction analysis was performed by LASSO, with the choice of 10-fold cross-validation and minimized error rate when λ=0.0992708 (Fig. 5A and B); as a result, Tolloid-like 1 gene (TLL1), Type XVI collagen (COL16A1), prostaglandin F2 alpha (PTGFR), cartilage intermediate layer protein 2 (CILP2), kinesin family member 26b (KIF26B), interferon inducible protein 27 (IFI27), growth arrest-specific gene 1 (GAS1) and chemokine receptor 7 (CCR7) were selected, a number of which were associated with the development of ovarian cancer. For example, COL16A1 is associated with the development of ovarian cancer (54), PTGFR is a potential serum marker for early diagnosis of ovarian cancer (55), upregulation of KIF26B enhances the proliferation and migration of ovarian cancer cells (56), and IFI27 promotes epithelial-mesenchymal transformation and induces ovarian tumorigenicity (57). The 8-gene signature was established by a multi-factor COX regression analysis using the following model: Risk Score=0.0603 × expTLL1 + 0.006 × expCOL16A1 + 0.1139 × expPTGFR + 0.0011 × expCILP2 + 0.0142 × expKIF26B - 0.0003 × expIFI27 + 0.0011 × expGAS1−0.0898 × expCCR7. The risk score of each sample was calculated. The results demonstrated that as the risk score increased, the survival time gradually reduced (Fig. 5C). In addition, when the expression levels s of IFI27 and CCR7 decreased, the expression levels of TLL1, COL16A1, PTGFR, CILP2, KIF26B and GAS1 increased (Fig. 5C). The ROC analysis demonstrated that the AUC was 0.83 (Fig. 5D), and when the samples were grouped according to the Youden's index of the ROC (cut-off, −0.1147472), a highly significant difference was observed in the prognosis between the high- and low-risk groups (Fig. 5E).

Figure 5.

(A) Trajectory change of each independent variable. The × axis represents the log value of the independent variable lambda; the y axis represents the coefficient of the independent variable. (B) Confidence intervals of lambda. (C) Risk score, survival time, survival status and expression of the 8-gene signature in the training set. (D) ROC curve of the 8-gene signature classification. (E) Kaplan-Meier survival curve based on the 8-gene signature in the training set. ROC, receiver operating characteristic; AUC, area under the curve; HR, hazard ratio; CI, confidence interval; L1 Norm, lambda.

The robustness of the 8-gene signature model

In order to verify the robustness of the 8-gene signature model, the risk score of each sample in the TCGA validation dataset was first calculated. The association between the risk score and gene expression was consistent with the training set (Fig. 6A), and the AUC in the TCGA validation dataset was 0.67 (Fig. 6B). The TCGA validation dataset samples were divided into high- and low-risk groups according to the threshold of the training set, and significant prognostic differences were identified between the two groups (Fig. 6C). The model was verified in the external validation set (GSE26193), and the ROC analysis demonstrated that the AUC was 0.63 (Fig. 6D), and the low-risk group exhibited a significantly better prognostic result compared with the high-risk group (Fig. 6E). Therefore, the model exhibited effective prognostic classification performance in the test and external validation sets.

Figure 6.

(A) Risk score, survival time, survival status and 8-gene expression in TCGA validation set. (B) ROC curve of the 8-gene signature classification in TCGA validation set. (C) KM OS curve based on the 8-gene signature in TCGA validation set. (D) ROC curves of OS based on the 8-gene signature classification in the GSE26193 validation set. (E) KM OS curve based on the 8-gene signature classification in the GSE26193 validation set. ROC, receiver operating characteristic; AUC, area under the curve; KM, Kaplan-Meier; OS, overall survival; HR, hazard ratio; CI, confidence interval.

Clinical independence of the 8-gene signature model

To identify the independence of the 8-gene signature model in clinical applications, relevant hazard ratio (HR), 95% confidence interval (CI) of HR and P-value were analyzed by univariate and multivariate COX regression in TCGA training set, TCGA test set and the GEO validation set. The clinical information of TCGA, GSE44001 patient records, including age, differentiation, clinical stage, and grouping information of the 8-gene signature were systematically processed (Table IV). In TCGA training set, univariate COX regression analysis demonstrated that the high-risk group and age were significantly associated with survival; however, the corresponding multivariate COX regression analysis identified that only the high-risk group (HR, 2.56; 95% CI, 1.82–3.59; P=5.66×10−8) exhibited clinical independence. In TCGA test set, univariate and multivariate COX regression analysis demonstrated that the high-risk group was significantly associated with survival (HR, 1.682; 95% CI, 0.729–3.882; P=0.022). In GSE44001, univariate COX regression analysis demonstrated that the high-risk group and stage were associated with survival; corresponding multivariate COX regression analysis revealed that the high-risk group (HR, 1.604, 95% CI, 0.494–50.041; P=0.017) and grade (HR, 2.203; 95% CI, 1.628–2.982; P=3.12×10−7) exhibited significant differences in predicting ovarian cancer prognosis. The expression trends of the eight genes were also analyzed in different groups based on age (>60 and ≤60 years), stage, grade and lymphatic invasion. Among them, IFI27 was significantly upregulated in the ≤60 years group, G3 and G4; COL16A1, PTGF and KIF26B were significantly upregulated in patients with lymphatic invasion (Fig. S3). TLL1, COL16A1, PTGFR, CILP2, KIF26B, IFI27, GAS1 and CCR7 were significantly upregulated in tumor samples compared with normal samples (Fig. S4). These results suggested that the 8-gene signature was a prognostic indicator independent of other clinical factors and had clinical value.

Table IV.

Univariate and multivariate Cox regression analysis of prognostic clinical factors and clinical independence.

A, TCGA training datasets
	Univariate analysis			Multivariate analysis
Variable	HR	95% CI	P-value	HR	95% CI	P-value
8-gene risk score, High vs. low	3.233	2.163–4.833	<0.001[a]	2.56	1.8227–3.592	<0.001[a]
Age, >60 vs. ≤60 years	1.016	1.002–1.031	0.020[a]	1.014	0.9994–1.028	0.986
Grade, G3/G4 vs. G1/G2	1.201	0.757–1.904	0.436	1.112	0.7914–1.563	0.539
Stage, III vs. I/II	1.991	0.812–4.877	0.131	1.216	0.7877–1.877	0.377
Stage, IV vs. I/II	1.402	0.865–2.272	0.170
B, TCGA test datasets
	Univariate analysis			Multivariate analysis
Variable	HR	95% CI	P-value	HR	95% CI	P-value
8-gene risk score, High vs. low	1.958	1.103–3.473	0.021[a]	1.682	0.729–3.882	0.022[a]
Age, >60 vs. ≤60 years	1.024	0.996–1.052	0.083	1.026	0.997–1.056	0.075
Grade, G3/G4 vs. G1/G2	1.118	0.477–2.618	0.796	1.875	0.880–3.997	0.103
Stage, III vs. I/II	1.615	0.222–11.760	0.635	0.956	0.399–2.294	0.921
Stage, IV vs. I/II	1.379	0.479–3.972	0.551
C, GSE26193
	Univariate analysis			Multivariate analysis
Variable	HR	95% CI	P-value	HR	95% CI	P-value
8-gene risk score, High vs. low	2.341	1.374–3.988	0.002[a]	1.604	0.494–50.041	0.017[a]
Grade, G3/G4 vs. G1/G2	1.062	0.661–1.709	0.801	2.203	1.628–2.982	<0.001[a]
Stage, III vs. I/II	3.799	1.959–7.366	<0.001[a]	0.729	0.483–1.099	0.132
Stage, IV vs. I/II	2.442	1.582–3.768	<0.001[a]

P<0.05. TCGA, The Cancer Genome Atlas; HR, hazard ratio; CI, confidence interval.

GSEA analysis of enriched pathway differences between the high- and low-risk groups

In TCGA dataset, GSEA was performed to determine the significantly enriched pathways in the high- and low-risk groups, and a total of 26 pathways were identified (Table SIX). In the high-risk group, the enriched pathways were mainly associated with the occurrence, invasion and metastasis of ovarian cancer, including ‘basal cell carcinoma’, ‘focal differentiation’, ‘pathways in cancer’ and ‘gap junction’ (Fig. 7A-D). In the low-risk group, mainly immune-related pathways were enriched, such as ‘antigen processing and presentation’, ‘intestinal immune network for IGA PRODUCTIO’, ‘primary immunodeficiency’ and ‘natural killer cell-mediated cytotoxicity’ (Fig. 7E-H). Thus, the 8-gene signature may be involved in important biological processes in the development and progression of ovarian cancer.

Figure 7.

(A) Basal cell carcinoma, (B) focal adhesion, (C) pathways in cancer, (D) gap junction, (E) antigen processing and presentation, (F) intestinal immune network for IgA production, (G) primary immunodeficiency and (H) natural killer cell-mediated cytotoxicity pathways in Gene Set Enrichment Analysis in The Cancer Genome Atlas dataset.

The performance of the 8-gene signature

By referring to previous studies, four robust prognostic risk models were identified, including a 3-gene (49), a 6-gene (50), an 8-gene (51) and a 5-gene (52) signature. To improve the comparability of the models, the risk score of each ovarian cancer sample in TCGA was calculated according to the corresponding genes in the four models by applying the same method. The ROC of each model was evaluated, the samples were divided into high- and low-risk groups according to the median risk score, and the OS differences between the two groups were calculated (Fig. 8). Significant differences were observed in the four models in predicting OS prognosis between the two groups (P<0.05). However, the AUCs of the four gene models were lower compared with the 8-gene signature developed in the present study, indicating that the 8-gene signature in the present study exhibited a better predictive performance.

Figure 8.

(A) ROC curve of the 3-gene signature. (B) KM curve of prognostic differences between high- and low-risk groups based on the 3-gene signature. (C) ROC curve of the 6-gene prognostic signature. (D) KM curve of prognostic differences between high- and low-risk groups based on the 6-gene signature. (E) ROC curve of the 8-gene signature. (F) KM curve of prognostic differences between high- and low-risk groups based on the 8-gene signature. (G) ROC curve of the 5-gene signature. (H) KM curve of prognostic differences between high- and low-risk groups based on the 5-gene signature. ROC, receiver operating characteristic; AUC, area under the curve; KM, Kaplan-Meier; HR, hazard ratio; CI, confidence interval.

Discussion

Previous studies have demonstrated that abnormal metabolism is one of the markers of cancer cells, and that there are differences between healthy and tumor cells in energy metabolism (58–60); in addition, various catabolic pathways, for example, glycolysis, OXPHOS and fatty acid metabolism, are involved in energy metabolism (61). In the current study, RNA-seq data were used to detect local energy metabolism-related gene expression status and their prognostic value for patients with ovarian cancer. The results identified three energy metabolism-related molecular subtypes, of which C1 was associated with a poor prognosis, and the profile of the immune cells in the tumor microenvironment in this subtype was different from the other two subtypes. Previous studies have demonstrated that the interaction between cancer cells and the tumor microenvironment affects cancer proliferation, energy metabolism, metastasis and recurrence (62), and energy metabolism serves an important role in immune regulation (63). The occurrence of aerobic glycolysis in tumor cells can shape the immune system by increasing the transcription of cytokines and inhibiting the differentiation of monocytes into dendritic cells (64,65). Therefore, abnormal energy metabolism may lead to different prognostic outcomes by altering the state of immune cells in the tumor microenvironment.

Ovarian cancer is a highly heterogeneous disease, as patients with similar Tumor-Node-Metastasis (TNM) staging have different survival times (66). At present, traditional clinicopathological indicators such as tumor size, vascular invasion, portal vein tumor thrombus and TNM staging do not satisfy the current needs in predicting individual outcomes, especially risk stratification, and the ‘one-size-fits-all’ treatment strategy has been demonstrated to be ineffective (67,68). Screening prognostic molecular markers that fully reflect the biological characteristics of tumors has value in individualized prevention and treatment of patients with ovarian cancer (69,70). In the current study, abnormal tumor energy metabolism was associated with poor prognosis. Energy metabolism-related gene expression levels were associated with clinical and molecular characteristics of patients with ovarian cancer. Thus, the present study developed a signature based on subtypes of abnormal energy metabolism that could stratify high- and low-risk patients. Of note, such energy metabolism-related genes may serve as a powerful prognostic indicator and help decide targeted therapies for patients with ovarian cancer based on energy metabolism.

In the 8-gene signature developed in the present study, TLL1, COL16A1, PTGFR, CILP2, KIF26B and GAS1 were identified as risk factors, whereas IFI27 and CCR7 were protective factors. Previous studies have reported that COL16A1 is involved in the progression of ovarian cancer (53), PTGFR is a potential serum marker for early diagnosis of ovarian cancer (54), upregulation of KIF26B promotes the proliferation and migration of ovarian cancer cells (55), and IFI27 promotes epithelial-mesenchymal transition and induces ovarian tumorigenicity (56). Thus, these genes are associated with the prognosis of ovarian cancer. The GSEA results of the present study revealed that the eight genes were enriched in the pathways and biological processes of ovarian cancer development. These results suggested that the signature may have clinical value and may provide a potential target for the diagnosis of ovarian cancer.

Although the association between the expression levels of energy metabolism-related genes and the prognosis of ovarian cancer were analyzed by bioinformatics and the characteristics related to energy metabolism were explored, the current study had limitations; for example, a number of samples lacked clinical follow-up information, and factors such as the presence of other diseases were not considered to distinguish their effects from those of prognostic biomarkers. In addition, the results were obtained only through bioinformatics analysis; other experiments should be performed to ensure the accuracy of the current results. Finally, key proteins in all metabolic pathways are under the control of post-translational modifications, and not necessarily by up/downregulation of their absolute levels.

In conclusion, the current study determined the expression levels of energy metabolism-related genes and their predictive values in ovarian cancer prognosis and established an 8-gene signature related to energy metabolism, which was be able to determine the risk levels in patients with ovarian cancer. The test and validation datasets exhibited high AUCs, and the results were independent of clinical features. Compared with clinical features, the 8-gene signature exhibited improved survival risk prediction for patients with ovarian cancer. Thus, the 8-gene signature developed in the present study may be used as a molecular diagnostic test in assessing the prognosis of patients with ovarian cancer.