MicroRNA molecules can act as an important tool for the detection of ovarian cancer. Chung et al. (15) reported let-7b, miR-26a, miR-132, and miR-145 as potential biomarkers in ovarian cancer patients. Among the results of Yuan et al.’s (16) study, has-miR-6784-5p, has-miR-6800-5p, and has-miR-5100 are indicating ovarian-associated cancer signature. Jeon et al. (17) reported that the serum and tissue miR-1290 was significantly elevated in patients with epithelial ovarian cancer compared with patients with benign ovarian neoplasm. Chen et al. (18) reported a total of 19 miRNAs, which were identified by random forest models, that were important in cancer diagnosis. In this study, the top five miRNAs with the highest frequency were chosen to be the biomarker candidates for cancer screening, which has-miR-3184-5p achieved a high rank. Yaghoobi et al. (19) proposed a method called EBST that has identified 11 serum miRNAs as potential biomarkers associated with ovarian cancer; among the miRNAs set, has-miR-1228-5p and has-miR-6784-5p were also reported. Zhang et al. (20) reported the four miRNA models that showed very strong performances with AUCs > 0.95 in the biliary tract, bladder, colorectal, esophageal, gastric, glioma, liver, ovarian, pancreatic, and prostate cancers. This study provides proof-of-concept data in demonstrating that the four miRNA (hsa-miR-5100, hsa-miR-1343-3p, hsa-miR-1290, and hsa-miR-4787-3p) model has the potential to be developed into a simple, inexpensive, and non-invasive blood test for the early detection of multiple cancers with high accuracy. Using statistical approaches, Hamidi et al. (21) identified 10 miRNAs regulated in ovarian serum cancer samples compared with non-cancer samples in the publicly available data set GSE106817: hsa-miR-5100, hsa-miR-6800-5p, hsa-miR-1233-5p, hsa-miR-4532, hsa-miR-4783-3p, hsa-miR-4787-3p, hsa-miR-1228-5p, hsa-miR-1290, hsa-miR-3184-5p, and hsa-miR-320b. However, the approach of the previous study (21) failed to take into account the non-linearity structure in big data; therefore, in this paper, we are implementing a new machine-learning variable selection approach called Boruta to address this problem. We will observe that the new miRNAs will be explored by the new method that has not been recognized in the traditional methods.

It is important to note that the choice of feature selection (FS) method should be tailored to the specific characteristics of the data set and research question at hand. Gene expression data are the representation of non-linear interactions among genes (22). By computing analysis of these data, it is expected to gain knowledge of gene functions and disease mechanisms. Statistical methods can only identify linear patterns, while non-linear patterns of relationships remain hidden. As mentioned in many research (23–29), Boruta has superior advantages in terms of feature selection accuracy, stability, and classification performance across different domains such as protein subcellular localization and credit risk assessment, however, especially in microarray data sets of ovarian cancer that have been rarely used before. This is based on some studies on the stability of Boruta (30–32) as a machine-learning method that can more accurately discover new miRNAs that were hidden in statistical methods. Therefore, this work attempts an innovation in two important issues: the identification of new miRNAs based on complex non-linear structures and the comparison of new results with the previous ones, which will be described in the results and discussion section.

We have used the GSE106817, GSE113486, and GSE113740 data sets from the GEO database, which is available at https://www.ncbi.nlm.nih.gov/geo/. The GSE106817 data set started on 13 November 2017 in Kanagawa, Japan, which is serum miRNA profiles of 4,046 women specimens, and which consists of 333 ovarian cancer and 2,759 non-cancer controls and 976 other types of cancer. The GSE106817 data set consists of ovarian cancer patients who were of mean age 57(±12) years, 25% stage I, 10% stage II, 55% serous, 19% clear cell, and 13% endometrioid histology (33). Three microarray data sets totaling to 6,835 unique participants including 728 ovarian cancer patients and 3,892 non-cancer controls were included in the current analysis, all derived from studies originating from a Japanese nationwide research project “Development and Diagnostic Technology for Detection of miRNA in Body Fluids” that is designed to characterize serum miRNAs in over 5,000 participants across several types of cancer using a standardized microarray platform. Supplementary Figure S1 clearly shows the stages of data pre-processing, identification of significant features or predictors, the model building of classifier algorithms, and performance evaluation, which are the four main phases of this analysis.

In cancer prediction models, statistical and machine-learning algorithms have been widely used, providing more accurate prognoses and lower per-patient costs. The high dimensionality of the gene expression profiles is a crucial issue when building cancer-predictive models (36). As a result, we used a machine-learning algorithm based on the random forest classifier, which is easily implemented in the Boruta package in R (37). In many studies involving miRNAs expression data, Boruta has been used to identify important features (38); this could help in the development of biomarkers for cancer diagnosis and prognosis. On the other hand, we used these techniques to characterize miRNAs with biomarker potential that may be useful in the diagnosis and/or prognosis of this disease, potentially assisting public health (39).

We cleaned and normalized the data using the min-max normalization method (40). Since gene expression data sets had too many irrelevant features for classification, feature selection was inevitable. Feature selection techniques can be used in data pre-processing to perform successful data reduction, which is beneficial for finding accurate data models (41). As noted, feature selection techniques have the benefits of reducing over-fitting and reducing model complexity with ease of understanding, as well as training models more quickly.

We split the data using the CARET package into two parts: two-thirds of the data were used for model development or training, while the remaining one-third of the data were used to evaluate or validate the model.

The validation technique is widely used to avoid over-fitting and to check the validity of the models. We evaluated our outcomes employing two external data sets, as shown in the Supplementary Figure S1. The metrics utilized to assess the results of the classification models are expressed below:Accuracy:ACC=TP+TNTP+FP+TN+FN,Sensitivity:SEN=TPTP+FN,Specificity:SPC=TNTN+FP,Kappa:k=Pr(a)−Pr(e)1−Pr(e)where: 1.

TP (true positive) is the number of people who suffer from “cancer” among those who were diagnosed with “cancer.”

As noted in the previous section, we find 10 miRNAs that are has-miR-1290, has-miR-1233-5p, has-miR-1914-5p, has-miR-1469, has-miR-4675, has-miR-1228-5p, has-miR-3184-5p, has-miR-6784-5p, has-miR-6800-5p, and has-miR-5100. We implemented each miRNA separately in models to get their power of prediction individually in classification between cancer and non-cancerous samples. The AUC of each of these miRNAs is listed in Supplementary Table S1A. We observe that in the internal validation, all miRNAs have high AUC (minimum AUC: 86.0%; maximum AUC is 96.8%). The performance measures for LR, DT, RF, ANN, and XGB models are shown in Supplementary Table S3A. We observe that the AUC of LR, RF, ANN, and XGB is 99.9%. Supplementary Table S3A shows the accuracy, sensitivity, specificity, NPV, PPV, and Kappa for LR, DT, RF, ANN, and XGB models in the classification and prediction of ovarian cancer. Four models obtained an AUC of 99.9%; however, DT obtained 98% AUC. In detail, RF has the highest value of accuracy (99.13), specificity (99.51), PPV (95.83), and Kappa (95.35), and LR have high sensitivity (98.96) and NPV (99.88). Figure 1A illustrates the ROC curve for the proposed models of 10 candidate miRNAs that are shown in Supplementary Table S1A. All models except DT have over 99.9% of AUC. Figure 1B shows the individual AUCs of 10 miRNAs in internal data set: has-mir-5100 (93.7%), has-mir-6800-5p (97%), has-mir-6784-5p (94.2%), has-mir-3184-5p (94.2%), has-mir-1228-5p (95.6%), has-mir-4675 (95.4%), has-mir-1469 (96.7%), has-mir-1914-5p (96%), has-mir-1233-5p (97.7%), and has-mir-1290(95.4%). In Supplementary Figure S2, we used a boxplot to display the expression levels of these 10 candidate miRNAs in the cancer and non-cancer groups. In the boxplots, it is clear that four of the miRNAs has-miR-1233-5p, has-miR-1914-5p, has-miR-4675, and has-miR-5100 have higher expression level with various cut-off for cancerous samples, and on average, four of them (has-miR-1228-5p, has-miR-3184-5p, has-miR-6784-5p, and has-miR-6800-5p) have lower expression level for cancerous samples. We used heatmap plots by implementing the “heatmaply” package to underpin the potential relationships between features and the hierarchical clustering analysis using the selected features to recognize different samples in the internal discovery data sets. Supplementary Figure S3 shows a promising result of the hierarchical clustering analysis (heatmap) using the 10 identified miRNAs to differentiate between cancerous and non-cancerous samples in GSE106817. The selected microRNAs are differently expressed in the non-cancer and cancerous classes. This is well illustrated by drawing the heatmap (Supplementary Figure S3).

Supplementary Table S1B,C demonstrate the performance of miRNAs individually in proposed classification models in external validation data sets, as seen by the fact that almost all miRNAs have higher AUC in GSE113486 than in GSE113740. But in detail of each data set, Supplementary Table S3 shows the results of the external data set models based on the Boruta feature selection algorithm. As seen in Supplementary Table S3B LR and ANN have a value of 100 in all seven criteria and XGB has an AUC, specificity, and PPV of 100. Supplementary Figure S4A shows that all models for GSE113486 yielded 100% AUC except DT and RF. Supplementary Figure S4B illustrates how biomarkers perform individually has-mir-5100 (95.8%), has-mir-6800-5p (99.7%), has-mir-6784-5p (97.5%), has-mir-3184-5p (93.9%), has-mir-1228-5p (99.8%), has-mir-4675 (99.1%), has-mir-1469 (100%), has-mir-1914-5p (99.8%), has-mir-1233-5p (99.4%), and has-mir-1290 (96.4%) in GSE113486. Boxplots show us that six of miRNAs (has-mir-1290, has-mir-1233-5p, has-mir-1914-5p, has-mir-1469, has-mir-4675, and has-mir-5100) have upregulated to ovarian cancer samples in GSE113486 (Supplementary Figure S5). In GSE113740, as the second external validation data set, we can see the result of AUC of LR, RF, ANN, and XGB over 94% in Supplementary Figure S6A. We also found AUC for these 10 miRNAs (individually) in external data sets that included individually has-mir-5100 (90.6%), has-mir-6800-5p (89.7%), has-mir-6784-5p (74.4%), has-mir-3184-5p (74.4%), has-mir-1228-5p (85.2%), has-mir-4675 (79.7%), has-mir-1469 (84.4%), has-mir-1914-5p (81.5%), has-mir-1233-5p (86.5%), and has-mir-1290 (91%) as shown in Supplementary Figure S6B. Supplementary Table S3C shows us that RF and XGB have the highest value in Kappa (72.96 and 71.96) in AUC and accuracy (97.2, 93.75), as seen ANN has 100 of sensitivity and NPV. Boxplots (Supplementary Figure S7) show us that six of miRNAs (has-mir-1290, has-mir-1233-5p, has-mir-1914-5p, has-mir-1469, has-mir-4675, and has-mir-5100) have high expression level in ovarian cancer samples. Supplementary Figures S8, S9 show a promising result of the hierarchical clustering analysis (heatmap) using the 10 identified miRNAs to differentiate between the cancerous and non-cancerous samples in GSE113486 and GSE113740, respectively.

The results of the biological insight section tell us about cell analysis for miRNAs that were found in this study based on the findings of the previous studies. The DIANA tool miRPath v.4 was used to perform the pathway enrichment analysis, based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The target genes of miRNA were identified using TargetScan v8.0 databases. The settings of the software were P-value threshold = 0.005 and the FDR correction filter were ticked. It should be mentioned that we used two methods to find the target genes: the first one is the genes union and the second is the pathway union. To investigate the efficiency of the set of biomarkers selected by Boruta and their superiority over the previous similar work done by Hamidi et al. (21), three groups of miRNAs were analyzed by miRPath v.4: (A) common biomarkers of the current study and the previous study by Hamidi et al. (21); (B) biomarkers selected by Boruta in the present study and not identified in the previous work; and (C) biomarkers that were selected in the previous study and were not identified in the current study. The list of genes of these three groups and their analysis results by miRPath v.4 tool are shown in Figure 2. As shown in Figure 2A, among the six common genes between the present and previous work, four genes are involved in at least one known cancer pathway (axon guidance). Among those four genes, hsa-miR-5100 and hsa-miR-1290 are involved in several well-known and important pathways in cancer. Figure 2B shows that among the four specific genes identified by the Boruta technique, three genes are involved in at least two well-known pathways in cancer, among which hsa-miR-4675 is involved in several pathways. However, in Figure 2C, among the four specific genes identified in the previous work of Hamidi et al. (21), only the hsa-miR-320b gene is involved in several important cancer pathways. It should be noted that there are six common paths between Groups A and B, while there are four common paths between A and C. This means that there are more correlation between genes of Group A and B than of Group A and C. This interpretation shows the biological superiority of Boruta's technique over the previous work. A comparison between routes of Group B and C also provides interesting results. Eight pathways are common between the two groups, which are proteoglycans in cancer, ErbB signaling, colorectal cancer, hepatocellular cancer, pathways in cancer, pancreatic cancer, axon guidance, and Hippo signaling. Axon guidance pathway is common among all the three groups. Many axon guidance molecules regulate cell migration and apoptosis in normal and tumorigenic tissues (63). Supplementary Table S4 shows the target genes of the selected microRNAs and the associated KEGG pathways from the genes union method, which indicates the significance of the relationship between the microRNAs and the corresponding pathways under the specified threshold values. Figure 3 shows the network of miRNAs and identified target genes. In this figure, transcription factors and LNC-RNAs have also been added through some studies. References for these interactions are described in Supplementary Table S4.

In Supplementary Table S5, we only selected seven pathways because only the pathways that had very high correlation with miRNAs were selected (including a P-value of < 0.002). Among the top seven pathways identified, based on P-value, were pathways associated with fatty acid biosynthesis, prion diseases, axon guidance, glioma, ErbB signaling pathway, proteoglycans in cancer, and endometrial cancer. All signaling pathways related to miRNAs were used from known pathways, and in general, they play an important role in all types of cancer. According to the KEGG database, some of the published articles confirm the role of some of the selected miRNAs in cancer directly. A number of these documents are summarized in Table 2. Figure 4 shows the predicted pathways of the effect of some of the selected microRNAs that have been taken from the https://targetexplorer.ingenuity.com/index.htm.

Figure 5 presents the common miRNAs between two related studies (18, 21) and miRNAs that were obtained in our study. There is some evidence in the literature for the biomarkers included in our study. Hamidi et al. (21) showed that hsa-miR-5100, hsa-miR-1233-5p, hsa-miR-4532, hsa-miR-1290, has-miR-3184-5p, and hsa-miR-320b could potentially be employed as important biomarkers in ovarian cancer. Jeon et al. (17) investigated that miRNA-1290 in the epithelial ovarian cancer group was significantly overexpressed in serum exosomes and tissues as compared with the benign ovarian neoplasm group. Ying et al. (90) expressed that microarray data analysis showed that hsa-miR-1290 was differentially expressed between COC1 (DDP-sensitive) and COC1/DDP (DDP-resistant) tumor cell lines. Chen et al. (18) showed that only five balanced miRNAs were determined to be important in cancer diagnosis: hsa-miR-663a, hsa-miR-6802-5p, hsa-miR-6784-5p, hsa-miR-3184-5p, and hsa-miR-8073. Furthermore, Chen et al. (18) found that hsa-miR-3184-5p can act as an early biomarker of bladder cancer and as a key regulator of breast cancer. Also, hsa-miR-6784-5p has been reported to be a sensitive serum biomarker for ovarian cancer diagnosis and a key regulator for breast cancer.

In the end, we note that although there are fundamental differences between microarray and RNA-Seq methods for obtaining gene expression data, the data matrix obtained from both methods is completely similar after performing the necessary pre-processing. Therefore, our method is also applicable to RNA-Seq data.