Identification of the miRNA signature associated with survival in patients with ovarian cancer

Identification of miRNA signature associated with survival time

First, we attempted to estimate the survival time of patients with ovarian cancer using miRNA expression profiles. We used a dataset containing 209 patients with ovarian cancer along with the expression profiles of 415 miRNA and survival time data. OV-SURV included the optimal feature selection algorithm IBCGA and identified a miRNA signature consisting of 39 miRNAs that could be used to predict the survival time of patients with ovarian cancer.

The OV-SURV-average yielded a correlation coefficient of 0.77±0.01 and mean absolute error of 0.69±0.02 years between real and estimated survival time. The data strongly suggested that the identified miRNAs were effective at estimating the survival time of patients with ovarian cancer. The system flow chart of this work is depicted in Figure 1.

We compared this OV-SURV method with standard multiple linear regression [36], LASSO [37], and Elastic net [38]. The LASSO method resulted in a correlation coefficient and a mean absolute error of 0.48 and 1.00 years, respectively. The Elastic net method resulted in a correlation coefficient and a mean absolute error of 0.55 and 1.00 years, respectively. The multiple linear regression method resulted in a correlation coefficient and a mean absolute error of 0.66 and 0.84 years, respectively. A comparison of these results is shown in Table 1. The correlation plots between the estimated and real survival time obtained using OV-SURV, LASSO, Elastic net, and multiple linear regression are shown in Figure 2.

Further, OV-SURV was validated on an independent dataset from the TCGA database to evaluate the prediction performance. The dataset consisting of 160 living patients who suffer from ovarian cancer with follow-up time was used for the validation. The follow-up times of 160 patients were up to 5 years and the mean and standard deviation of the follow-up times were 20.42 ± 16.30 months. The mean estimated survival time of OV-SURV was 33.85 ± 5.24 months. There were 123 patients whose estimated survival time was longer than their actual follow-up time. The result revealed that OV-SURV achieved the classification accuracy of 77% for estimating the survival of patients with ovarian cancer. The mean absolute error of predicting the remaining 37 patients was 0.93 years which is slightly longer than the 0.69 years in Table 1. The prediction performance of OV-SURV on an independent test cohort was measured based on our previous study [18]. The prediction result of OV-SURV for each of the 160 patients is shown in Supplementary Figure 1.

Ranking of the identified miRNA signature

To prioritize the miRNAs of the identified miRNA signature, miRNAs were ranked based on their contributions to the survival estimation using main effect difference (MED) analysis. There are 39 miRNAs in the identified miRNA signature. The 39 miRNAs of the signature, miRNA IDs, and MED scores are shown in Table 2. The top 10 miRNAs ranked by MED analysis were hsa-miR-19b, hsa-let-7f, hsa-miR-323, hsa-miR-1978, hsa-miR-128, hsa-miR-1237, hsa-miR-486, hsa-miR-98, hsa-miR-933, and hsa-miR-889, and these miRNAs were then analyzed further. Among the miRNAs in the signature identified by OV-SURV, these 10 miRNAs were significantly associated with survival time. We then assessed the involvement of these miRNAs and their corresponding genes in cancer and various biological pathways.

Validating the prognostic potential of top 10 ranked miRNAs

We validated the prognostic power of the top 10 ranked miRNAs using Kaplan-Meier (KM) survival curves. The TCGA dataset contains 486 tumor samples were considered for this analysis. Survival analysis was performed using Kaplan-Meir plotter [39]. Five of the top 10 ranked miRNAs, hsa-let-7f, hsa-miR-1237, hsa-miR-98, hsa-miR-933, and hsa-miR-889, were significantly associated with survival in patients with ovarian cancer. The survival groups were distributed using the log-rank test, a p-value<0.05 was considered the cut-off to describe the statistical significance in the survival analysis. These five miRNAs, hsa-let-7f, hsa-miR-1237, hsa-miR-98, hsa-miR-933, and hsa-miR-889, have p-values of 0.01, 0.00033, 0.007, 4.9e-0.5, and 0.011, respectively; and hazard ratios of 0.73, 0.61, 0.73, 1.62, and 1.4, respectively.

Pathway analysis of the identified miRNA signature

The biological relevance of the selected miRNAs was analyzed using the DIANA tool. The pathway analysis helps to determine whether these miRNAs were statistically significant in Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways, including signalling pathways, cell processes and cancer pathways. The analysis revealed that the top 10 selected miRNAs were involved in several KEGG pathways, the most significant pathways including the transforming growth factor-β signalling pathway (hsa04350; p=3.29e-06), proteoglycans in cancer (hsa05205; p=7.93e-05), sphingolipid metabolism (hsa00600; p=0.0001), hippo signalling pathway (hsa04390; p=0.0005), and adherens junction (hsa04520; p=0.009). The details of the KEGG pathway analysis for the top 10 ranked miRNAs and their number of target genes are shown in Supplementary Table 1.

Additionally, the KEGG analysis of all 39 miRNAs in the signature revealed that few miRNAs were highly enriched in fatty acid biosynthesis, fatty acid metabolism, ECM/receptor interactions, and the hippo signalling pathway. The most highly significant pathway found in the KEGG pathway analysis was fatty acid biosynthesis (p<1E-325). There were four miRNAs (hsa-miR-182, hsa-miR-34a, hsa-miR-342, and hsa-miR-1304) that were highly enriched in fatty acid biosynthesis (hsa00061) and targeted genes such as FASN, ACSL4, and ACACA. Additionally, five miRNAs, namely, hsa-let-7f, hsa-miR-34a, hsa-miR-342, hsa-miR-1304, and hsa-miR-24, targeted the genes FASN, PTPLB, SCD, ACSL4, ECHS1, and TECR, which were all involved in fatty acid metabolism (p=1.02E-14). Fatty acid synthase (FAS) is an enzyme responsible for the synthesis of fatty acids and has been identified in most human carcinomas. FAS protein expression is associated with poor prognosis in both prostate and breast cancers [40, 41]. An in vivo microarray study revealed that there was a significant relationship between phosphorylated AKT and the expression of FAS; inhibition of FAS activity resulted in down regulation of phosphorylated AKT, which initiates apoptosis in ovarian cancer cells [42]. Inhibition of FAS and activation of AMP-activated protein kinase is selectively cytotoxic to SKOV3 human ovarian cancer cells [43]. Additionally, inhibition of FAS is significantly associated with survival in xenograft models [44]. This analysis suggested that the selected miRNAs were relevant to cancer and related biological signalling pathways. The heatmap showing KEGG pathway enrichment analysis of the miRNA signature is shown in Figure 3A.

Thus, this analysis indicated that fatty acid synthesis in ovarian cancer cells played an essential role in ovarian cancer progression. The identified miRNAs involved in fatty acid synthesis and fatty acid metabolism may possess important roles in ovarian cancer survival. The results of the KEGG pathway analysis for the miRNA signature and its target genes (various pathways, numbers of miRNAs involved, and the corresponding p-values) are shown in Supplementary Table 2. Fatty acid biosynthesis and fatty acid metabolism flow diagrams are shown in Supplementary Figures 2, 3.

Functional annotations of the miRNA signature

Functional annotations of the identified miRNA signature were analyzed using GO enrichment analysis. GO annotation analysis revealed that the top 10 ranked miRNAs were involved in several biological processes (BPs), cellular component (CCs), and molecular functions (MFs). The top 10 ranked miRNAs were involved in several BPs, such as gene expression (GO:0010467; p<1e-325), biosynthetic process (GO:0009058; p<1e-325), cellular nitrogen compound metabolic process (GO:0044403; p<1e-325), symbiosis, encompassing mutualism through parasitism (GO:0044403; p=2.78e-14), and viral process (GO:0016032; p=3.14e-13) by targeting 82, 360, 419, 68, and 62 genes, respectively, to name a few.

The top 10 ranked miRNAs were involved in various CCs, includes organelle (GO:0043226; p<1e-325), cytosol (GO:0005829; p=4.46e-13), protein complex (GO:0043234; p=2.55e-11), nucleoplasm (GO:0005654; p=9.93e-11), and microtubule organizing center (GO:0005815; p=0.004). The top 10 ranked miRNs involved in MFs, such as ion binding (GO:0043167; p<1e-325), RNA binding (GO:0003723; p=2.22e-15), enzyme binding (GO:0019899; p=2.09e-08), poly(A) RNA binding (GO:0044822; p=1.16e-07), protein binding transcription factor activity (GO:0000988; p=4.34e-05), and nucleic acid binding transcription factor activity (GO:0001071;p=0.006). In GO annotation analysis, among the BPs most genes involved in cellular nitrogen compound metabolic process (419 genes), in CCs most genes involved in organelle (834 genes), and in MFs, most genes are involved in ion binding (447 genes). The details of top 10 ranked miRNAs involvement in GO annotations is listed in Supplementary Table 3.

Go enrichment analysis revealed that miRNA signature is highly enriched in gene expression, catabolic process, RNA binding, cytosol, cellular nitrogen compound metabolic process, viral process cellular protein modification process, organelles, and enzyme binding. GO analysis of the miRNA signature is shown in Figure 3B. The results of the GO analysis for the miRNA signature and its target genes are shown in Supplementary Table 4.

Expression difference of the top ranked miRNAs across stages

Further, we wished to enquire about the relative expression levels of the top ranked miRNAs across different stages of ovarian cancer. We employed relative miRNA expression analysis across stages of ovarian cancer using UALCAN web portal [45]. The expression analysis results provide the significant expression difference of top 10 ranked miRNAs across different stages of patients with ovarian cancer. Of the top 10 ranked miRNAs, except the miRNA, hsa-miR-1978, remaining miRNAs are significantly expressed across stage 2, 3, and 4 of patients with ovarian cancer. Statistical significance of the relative miRNA expression across stages of ovarian cancer patients is provided in Supplementary Table 5. Box-whisker plot representation of relative expression difference across stages of ovarian cancer patients is given for the top ranked miRNAs in Supplementary Figure 4.

Roles of the top ranked miRNAs in ovarian cancer

The roles of the top 10 ranked miRNAs in ovarian cancer were analyzed using experimentally validated literature.

hsa-miR-19b: A quantitative real-time PCR study involving ovarian cancer cells showed that hsa-miR-19b was significantly expressed between cancer and control group with a p-value of 0.035 [46]. Over expression of hsa-miR-19b facilitated the invasion and migration of ovarian cancer cells (CAOV-3) and serves as an oncogenic role in ovarian cancer progression [47].

hsa-let-7f: Zheng et al. identified lower expression levels of hsa-let-7f in ovarian cancer samples when compared to controls, further lower levels of hsa-let-7f associated with poor prognosis of patients with ovarian cancer [48]. Elevated expression level of let-7f was found in primary ovarian carcinomas and suggesting a role of this miRNA in tumor progression [49].

hsa-miR-323: In silico miRNA expression analysis in ovarian cancer reported that dysregulation of hsa-miR-323 in ovarian cancer cells [50]. Differential expression of miR-323 was observed in epithelial ovarian cancer cells and miR-323 was significantly down-regulated in epithelial ovarian cancer cell lines compared to ovarian surface epithelium cells [14].

hsa-miR-128: hsa-miR-128 downregulates the colony stimulating factor-1, a key regulator of ovarian cancer, resulted into reduction of cell motility and adhesion in ovarian cancer cells [51].

hsa-miR-486: hsa-miR-486 downregulated the estrogen receptor alpha (ERα)-mediated olfactomedin 4, which contributes to ovarian cancer progression in cancer cells [52].

hsa-miR-98: A qRT-PCR study epithelial ovarian cancer revealed that hsa-miR-98 regulates the cisplatin resistance of epithelial ovarian cancer cells and was associated with poor outcome of patients with ovarian cancer [53].

A summary of the top 10 ranked miRNAs and their roles in ovarian cancer are shown in Table 3. In summary, six miRNAs, namely, hsa-miR-19b, hsa-let-7f, hsa-miR-323, hsa-miR-128, hsa-miR-486, and hsa-miR-98, had experimentally validated support to confirm their associations in ovarian cancer. The remaining four miRNAs, such as hsa-miR-1978, hsa-miR-1237, hsa-miR-933, and hsa-miR-889, were not previously reported to have roles in ovarian cancer. However, they were actively involved in other types of cancer, such as kidney cancer [54], breast cancer [55], gastric cancer [56], and squamous cell carcinoma [57]. Additionally, we verified the importance of these four miRNAs, hsa-miR-1978, hsa-miR-1237, hsa-miR-933, and hsa-miR-889 in ovarian cancer by performing overall survival analysis using Kaplan-Meir survival curves. The miRNAs, hsa-miR-1978, hsa-miR-1237, hsa-miR-933, and hsa-miR-889 obtained p-values of 0.22, 0.01, <0.001, and <001, respectively. The three miRNAs, hsa-miR-1237, hsa-miR-933, and hsa-miR-889 were significantly associated with the overall survival in patients with ovarian cancer. The KM survival plots for hsa-miR-889, hsa-miR-1237, and hsa-miR-933 are shown in Supplementary Figure 5.

Furthermore, we validated the involvement of top 10 ranked miRNAs in various cancers including ovarian cancer using experimentally validated studies. All of the top 10 ranked miRNAs were involved in major types of cancers, in which 6 miRNAs, hsa-miR-19b, hsa-let-7f, hsa-miR-323, hsa-miR-128, hsa-miR-486, and hsa-miR-98 were involved in ovarian cancer progression. The top 10 ranked miRNAs and experimentally validated evidences across various cancers is listed in Supplementary Table 6.

Differential expression of miRNA signature in ovarian cancer

Furthermore, to examine the involvement of the miRNA signature and the expression levels of the miRNAs of the signature in ovarian cancer, we predicted the association of the miRNA signature with ovarian cancer using dbDEMC 2.0 and HMDD v3.0. The prediction results showed that 29 miRNAs were differentially expressed in ovarian cancer. The expression levels of the remaining 10 miRNAs had not been previously reported in ovarian cancer, although these miRNAs were known to be actively involved in other major cancer types. These results suggested that these 10 miRNAs, including hsa-miR-323, hsa-miR-1978, hsa-miR-933, hsa-miR-616, hsa-miR-1245, hsa-miR-19b-2, hsa-miR-342, hsa-miR-513a, hsa-miR-577, and hsa-miR-1295, should be further evaluated to determine their roles in ovarian cancer survival. The prediction list of the identified miRNA signature associated with ovarian cancer is shown in Supplementary Table 7.

Additionally, we predicted the miRNA gene target network for the top ranked miRNAs using Cytoscape version 3.6. There are total 5009 predicted gene targets, in which 3893 interactions are predicted using MicroCosm version 5, and 116 interactions are predicted using TargetScan version 6.2. The miRNA-gene target interaction network is visualized in Supplementary Figure 6, for better visualization, only gene targets predicted by TargetScan are shown.

Dataset

We obtained miRNA expression data from patients with ovarian cancer using the TCGA database. The number of patients with ovarian cancer in the TCGA database was 586. We downloaded level-3 miRNA expression data from the TCGA portal, and the miRNA profiling was performed using the Illumina HiSeq 2000 miRNA sequencing platform. We filtered the dataset using the following criteria: (i) only patients with clinical data and survival information were included; (ii) to statistically remove the bias, the patients with a survival period of less than 30 days were excluded; and (iii) all patient lists and the corresponding survival periods were merged into a single data file to eliminate duplicate entries. After this filtering process, we identified 209 patients with the expression profiles of 415 miRNAs along with corresponding clinical data and days to death. Additionally, we used the dataset consisting of 160 living patients who suffer from ovarian cancer with follow-up time for the validation. Clinical information of the independent test cohort is listed in Supplementary Table 8.

OV-SURV

The proposed method OV-SURV was designed to estimate the survival time of patients with ovarian cancer by identifying a miRNA signature and then to discover their potential roles of miRNA biomarkers in ovarian cancer. SVM is effective for solving classification and regression problems [60] and can interpret the property values of a large number of samples in a multidimensional space. Owing to its effective regression abilities, SVR has been used to solve many biological prediction problems. The optimisation problem of ν-SVR for finding a model w can be defined as follows:

$$\min [\{\frac{1}{2} w^T(\varnothing (x_i)+b)+C(\nu \epsilon +\frac{1}{m}{\displaystyle {\sum }_{i=1}^m({\xi }_i+{\xi }_i^{*}))}\}]$$ (1)

Here, 0 ≤ ν ≤ 1, ξ_i ≥ 0, (x₁, y₁)…(x_m, y_m) F0CE R^d ×R are training data, y_i is the target values, x_i is the feature vector. We used d is the number of features in each instance and m is the number of training instance where d=415 and m=209 in this study. Let C>0 be the regularisation parameter, ε≥0 be an error sensitivity parameter, and b is a constant. The parameters C and ε are used to avoid overfitting the training data. The parameters C, ε, a Gaussian kernel γ and the feature selection are simultaneously optimized by OV-SURV while maximizing the estimation accuracy.

IBCGA

The feature selection algorithm IBCGA uses an intelligent evolutionary algorithm (IEA) [34], which can efficiently solve large parameter optimisation problems, and is good at deriving an optimised SVR model with feature selection. IBCGA can automatically select a small set of informative features from a large number of candidate features. OV-SURV identifies a signature of informative miRNAs associated with ovarian cancer survival based on the cooperation of IBCGA and ν-SVR using the objective function of maximising the correlation coefficient (CC) in terms of 10-CV.

$CC=\frac{{{\displaystyle \sum }}_{i=1}^N(x_i-\overline{x}) (z_i-\overline{z})}{\sqrt{\left[{{\displaystyle \sum }}_{i=1}^N{(x_i-\overline{x})}^2\right]\left[{{\displaystyle \sum }}_{i=1}^N{(z_i-\overline{z})}^2\right]}}$ (2)

where x_i and z_i are real and predicted survival time of the i-th miRNA, $\overline{x}$ and $\overline{z}$ are the corresponding means, and N is the total number of instances. Mean absolute error (MAE) is also used for estimating the prediction performance. The MAE is defined as follows

$MAE= \frac{1}{N}{\displaystyle \sum }_{i=1}^N{\left|x-z\right|}^2$ (3)

In this work, the LibSVM package [61] was used for implementation of ν-SVR.

The chromosome of IEA comprises 415 binary genes representing miRNAs for feature selection and three 4-bit genes for encoding parameters γ, C, and ν of ν-SVR. The chromosome encoding was designed as described in previous studies [17, 34]. The IBCGA can simultaneously obtain a set of solutions, Xr, where r = r_end, r_{end +} 1,…, r_start in a single run. The feature selection algorithm IBCGA used can be described as follows:

Step 1: (Initialisation) Randomly generate an initial population of Npop individuals. In this work, Npop=50, r_start = 50, r_end = 10, r = r_start.

Step 2: (Evaluation) Evaluate the fitness value of all individuals using the fitness function (CC) of 10-CV.

Step 3: (Selection) Use a conventional method of tournament selection that selects a winner from two randomly selected individuals to generate a mating pool.

Step 4: (Crossover) Select two parents from the mating pool to perform an orthogonal array crossover operation of IEA.

Step 5: (Mutation) Apply a conventional bit mutation operator to parameter genes and a swap mutation to the binary genes for keeping r selected features. The best individual was not mutated for the elite strategy.

Step 6: (Termination test) If the stopping condition for obtaining the solution Xr is satisfied, output the best individual as the solution Xr. Otherwise, go to Step 2.

Step 7: (Inheritance) If r > r_end, randomly change one bit in the binary genes for each individual from 1 to 0; decrease the number r by one and go to Step 2. Otherwise, stop the algorithm.

Step 8: (Output) Obtain a set of m miRNAs from the chromosome of the best solution Xm among the solutions Xr, where r = r_end, r_{end +} 1,…, r_start.

Multiple linear regression analysis

We employed multiple regression to estimate the survival time in patients with ovarian cancer [36]. A general multiple linear regression can be defined as

$y_i={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\cdots +{\beta }_n{x}_n+\epsilon$ (4)

KEGG pathway and GO term analysis

We used the DIANA-mirPath version 3 web-based server to analyze the miRNA profiles [62]. The DIANA-micro-T-CDS algorithm provided the predicted miRNA targets for the pathway analysis. The p-value threshold was set to 0.05 and Fishers’s exact test (hypergeometric distribution) was used for the enrichment analysis. We utilised the Network of Cancer Genes (NCG) [63], miRTarBase [64], and GeneCards [65] databases to obtain the cancer-related genes and predict gene targets.

Availability of data and materials

All the data used in this analysis can be found at TCGA data portal [https://cancergenome.nih.gov/].