The miRNA expression profiles of KIRC used in this study are obtained from the TCGA database. There are 611 samples, in which cancer and adjacent normal samples are 539 and 72, respectively. The miRNAs with zero variance in expression are excluded and a set with 1343 miRNAs is obtained. A normalization is made on each miRNA using FPKM.

Feature selection is a common method to eliminate redundant and irrelevant variables with the aim of reducing feature dimensionality and discovering a valid subset of genes as a diagnostic biomarker to make data fitting simpler and prediction more accurate [29, 30].

The cause of KIRC is generally due to mutation of some genes. However, the probability of mutation is small, which means that only a small number of mutations in cancer and adjacent normal cells lead to differential expressions and even cellular carcinogenesis. Therefore, we speculate that the number of the most differential miRNAs is countable.

In addition, there are often multiple miRNAs that act in combination in a biochemical reaction. That is, there are single miRNAs that do not express differently between cancer and normal tissues, but two or more of these miRNAs that together show differential expressions. In that case, the miRNAs that have a difference in combination cannot be found by conventional methods. To treat this situation, full enumeration is considered. To discover the optimal miRNAs for discrimination between KIRC and normal samples, each miRNA tuple in each dimensional combination is tried. For miRNA expression profiles, the number of miRNAs is usually much larger than the number of samples, so the time complexity of the full permutation is O(n!), which is too computationally intensive.

Thus, a random miRNA selection way is considered. Firstly, m miRNAs are randomly selected every time, and a corresponding classification result such as accuracy is calculated on the expression profiles corresponding to cancer and normal tissues. If one of the n miRNAs is viewed to be important for discrimination between KIRC and normal samples, interference with the miRNA will significantly impact the classification result.

We propose a method to quantify this impact. First of all, 70% of the data set is randomly divided as a training set. Classifiers such as logistic regression (LR), decision tree classifier (DTC), support vector machine (SVM), k-nearest-neighbor (KNN), multinomial naïve Bayes (MNB), and Fisher’s linear discriminant (LDA) are trained by using the Python package Scikit-learn [31]. Then, predictions are made on the 30% of the left data. That is, the classification accuracy $Scor{e}_1$is obtained. After that, a one-time permutation on expressions of miRNA i is made. The classifiers are then used to make predictions on the processed data with classification accuracy to be $Scor{e}_2$. The importance of the $i^{th}$ miRNA for the accurate classification of KIRC and adjacent normal samples is calculated as follows,

(1)$$impac{t}_i=Score1-Score2$$

where $impac{t}_i$ is the importance score of the $i^{th}$ miRNA.

Due to the uneven distribution of positive and negative samples, the classification accuracy is modified and expressed as follows,

(2)$$Accuracy=\frac{\frac{TP}{TP+FN}+\frac{TN}{TN+FP}}{2}$$

where FN, TP, FP and TN represent false negative, true positive, false positive and true negative, respectively. Positive and negative samples refer to KIRC and normal tissues.

The above method is repeated for n rounds in order to count as many miRNA combinations as possible. Then the average importance scores corresponding to all the listed miRNAs are obtained.

Considering the small sample size with unknown distribution, six classifiers (i.e., DTC, MNB, KNN, LDA, LR, and SVM) are utilized for training and validation.

Once the importance score of each miRNA is obtained, how to select the important genes becomes a new problem. If the importance score threshold is set manually, subjective factors will be inevitably introduced. Here we choose to use an unsupervised clustering approach to select miRNAs according to the importance scores.

A clustering algorithm based on probability distribution was used in this study, which is implemented based on the mixture module of the Python Scikit-learn package [31]. It is assumed that the sample points can be divided into k clusters, also known as k components, and each cluster obeys a different Gaussian distribution. The probability of each sample belonging to each distribution is calculated, and the sample is classified into the cluster with the highest probability. Then, based on the expectation maximization (EM) algorithm, the parameters of the Gaussian distribution are updated after several iterations until the model converges to a locally optimal solution.

Finally, the cluster with the highest importance scores is used as the important feature subset. That is

(3)$$p(\mathrm{X})=\sum _{\mathrm{i}=1}^{\mathrm{k}}{\alpha }_{\mathrm{i}}*p\left(X\mid {\mu }_i,{\mathrm{\Sigma }}_i\right)$$

where k is the number of components. ${\alpha }_i$ is the probability that the sample point belongs to the i${}^{th}$ Gaussian distribution, and $p(X{\mu }_i,{\mathrm{\Sigma }}_i)$ is the probability density function of the i${}^{th}$ Gaussian distribution, where ${\mu }_i$ is the mean vector and ${\mathrm{\Sigma }}_i$ is the covariance matrix.

After clustering on the importance scores, important feature subsets are obtained from different classifiers. Because the decision boundaries of these classifiers are different, the obtained important feature subsets are different.

To obtain the optimal feature subset, an ensemble strategy is proposed based on the idea of ensemble learning to combine all important feature subsets. The simplest way is to select the important feature subset by the best classifier. But it may lead to model overfitting. The union among these important feature subsets is calculated to get the optimal feature subset.

The combination of the results from individual classifiers can reduce the possibility of feature missed selection caused by a large hypothesis space, reduce the learning algorithm into local optimum, and thus improve the generalization effect of the model. It can also expand the hypothesis space and learn better approximations.

Therefore, this strategy of training the classifiers individually to find the important feature subsets and then combination can find the optimal features that are closer to the real important features and perform better on other data sets.

Evaluations are performed for the selected genes. One is on the testing set, and the other is on the new data sets built from the TCGA data set with only the selected features. In addition, we made cross-validation to test the validity of the obtained significant genes.

First, the important feature subsets derived from six different classifiers are evaluated. Each time, an important feature subset and all genes are used to train the classifier separately. We use the following evaluation metrics for accuracy evaluation based on the confusion matrix, i.e., TP rate, FP rate, Precision, Recall, F1 measure, and the accuracy after balancing positive and negative samples, as are expressed as follows,

(4)$$T{P}_{rate}={\displaystyle \frac{TP}{TP+FN}},$$ $$F{P}_{rate}={\displaystyle \frac{FP}{FP+TN}},$$ $$Precision={\displaystyle \frac{TP}{TP+FP}},$$ $$Recall={\displaystyle \frac{TP}{TP+FN}}$$ $$F1measure={\displaystyle \frac{2*precision*Recall}{Precision+Recall}}.$$

Second, the ensemble feature subset and the differential genes are evaluated. Three traditional differential analysis methods are used to search for the differential genes. The first two are performed by the R package limma and edgeR, and the last one is performed by a simple two-sample t-test. The p-value and logFoldChange value of each gene are calculated after the differential analysis to screen out differential genes. The threshold for p-value is set to 0.05, and the threshold for logFoldchange is calculated as follows,

(5)$$\log F{C}_{cutoff}=\overline{|\log FC|}+2*\sigma (\overline{\mid \log FC}\mid ).$$

Then new data sets are built on the selected genes. The same evaluation metrics for accuracy evaluation based on the confusion matrix are used.

Apart from that, one 5-fold cross-validation is made to verify the stability and reliability of the results. The TCGA dataset is divided into five mutually exclusive subsets, and the gene importance is calculated separately after multiple training and testing. Five sets of selected genes are obtained.

To validate the effectiveness of the obtained optimal feature subset, independent testing sets of KIRC are selected from the GEO database.

Combinations of any pair of miRNAs are enumerated to see the effect of different combinations on the sample grouping, and then all the important miRNAs are also considered as a combination.

Due to the high dimension of the data, it’s difficult to observe the distribution. To see the effect of the optimal feature subset on the sample grouping, the principal component analysis method (PCA) is used. PCA retains most of the effective information while reducing the number of features. Here, only two principal components are finally retained for data visualization.

The data were equally divided into halves, each of which corresponded to the training set and testing set. In the training set, 70% of the samples were randomly selected to train the six classifiers and the remaining out-of-bag samples were used to calculate the importance scores of each miRNA according to Eqn. 1. The procedure was repeated 700,000 times. The average importance score of each miRNA was obtained, as shown in Fig. 1. It can be seen that the miRNAs with high important scores account for only a fraction of all the miRNAs, about 1%, with miR-621, miR-210, and miR4456 ranking high in most of the classifiers. Finding common miRNAs using different classification methods indicates the validity of our method. The phenomenon that different miRNAs appear with high important scores indicates that sample distribution affects the classification results according to different classifiers with different classification boundaries.

For miRNAs that are not involved in the cancer process and are expressed at similar levels to non-cancerous cells, they are not important in the prediction of KIRC. According to the previous hypothesis, only a small number of variants of miRNAs cause KIRC. Therefore, there will be a large number of miRNAs with low important scores. In contrast, the important miRNAs are viewed as outliers. This can be verified in Fig. 2, which shows the effectiveness of the clustering method. The important genes are selected based on the Gaussian boundaries presented by short blue horizontal lines shown in Fig. 2.

To verify the validity of the found miRNAs, six classifiers were separately trained on the training data, and the classification accuracies were calculated on the independent test set. The classification results are shown in Table 1. It can be seen that LR and LDA using all the miRNAs for classification are more efficient than those using the selected miRNAs. And the higher accuracy of LR also means that there may be overfitting. However, DTC, SVM, KNN, and MNB achieve better classification results considering the selected miRNAs. Especially, there is no good hyperplane using SVM with all the miRNAs considered, while it shows much better classification results under the selected miRNAs. The corresponding confusion matrices are shown in Figs. 3,4.

To compare the optimal feature subset with the differential genes found by traditional differential analysis methods, four new datasets were constructed and used for evaluation.

Volcano plots of the differential genes are shown in Fig. 5. The differential genes are divided into up-regulated genes and down-regulated genes. Blue points in Fig. 5 represent down-regulated genes, and red points in Fig. 5 represent down-regulated genes. The top 10 differential genes are labeled with gene names.

By making a comparison between the proposed method with the prevailing feature selection methods such as edgeR, limma, and t-test, the quantitative assessment metrics of the classification results on the four data sets are shown in Table 2, and the corresponding ROC curves and AUCs are shown in Fig. 6.

From Table 2 we can see that the significant genes found by this method achieved the same or even better classification results than traditional methods, using only 15 selected genes, which is much smaller than the number of genes found by traditional differential analysis.

Each pair of important miRNAs was enumerated and GSE independent testing sets GSE109368, GSE151423, and GSE 151419 were used to verify the effectiveness of each gene combination. After enumeration, it is found that some pair of miR-210, miR-140 and miR-1270 are more effective in dividing the sample group as shown in Fig. 7.

All the important miRNAs were also treated as a combination, the features were compressed into two dimensions by PCA as shown in Fig. 8.

From Figs. 7,8, it can be found that the cancer samples show some degree of separation from the normal samples, which verifies the validity of the important miRNAs.