The investigation of the role of small bio-molecules (called microRNAs) in biological functions is a very popular topic in bioinformatics, since microRNAs have been shown to present novel therapeutic methods for diseases like cancer or Hepatitis C. In order to predict the involvement of microRNAs in biological functions many statistical approaches have been used that involve p-value calculations, with the most popular one being Fisher's exact test. However, it has been shown that data distribution does not match with any of the theoretical distributions used by the aforementioned approaches. Thus, an empirical randomization approach is preferred. Nevertheless, such analyses are computationally intensive. In this paper, we present a novel approach for microRNA enrichment analysis using Machine Learning techniques, in order to predict p-values instead of calculating them using randomization experiments. This simplifies the work for bioinformatics data analysts, helping them to eficiently perform multiple enrichment analysis tasks.
Transcriptomics, a popular sub-field of bioinformatics, deals with
RNA biomolecules that are produced by DNA transcription and
either contain instructions on how to produce proteins or
perform other functions vital to life. One such class of RNAs are
microRNAs (or miRNAs) which are short biomolecules that
“target” specific genes and block their expression. This means that
they afect biological functions by preventing protein production.
Dysregulation of certain miRNAs has been linked to diseases
like cancer, neuronal diseases, inflammatory diseases, etc [
The most popular statistical test used in the relevant
literature is Fisher’s exact test [
Indeed, it was shown in 2015 [
Generally, such approaches measure the strength of association based on the real data, without making assumptions about the structure of the data like theoretical statistical tests. In order to test whether a group of miRNAs is associated to a specific biological process, the method calculates a statistical measure relevant to the miRNA group and compares it to the statistical measure calculated for a large number of randomly assembled miRNA groups. Thus, if the number of random groups presenting a better behaviour (in terms of the statistical measure) is small, then the significance of the association between the query and the biological function under examination is large.
However, this approach is computationally intensive, because
it involves a very large number of operations between sets (union
and intersection) leading to large execution times (e.g., in the
order of hours, on 8 cores of a Intel Xeon CPU). For this reason,
attempts have been made to alleviate this issue, either by making
the set operations faster [
In this paper, we present a novel approach for miRNA enrichment analysis using machine learning techniques, to predict p-values using features of miRNA groups, relevant to the problem, instead of calculating them using randomization experiments. This simplifies the work for bioinformatics data analysts, helping them to eficiently perform multiple enrichment analysis tasks. Our contributions are: (a) framing the problem, (b) data set creation and feature engineering, (c) determining a shortlist of promising machine learning models, using cross-validation, and (d) fine-tuning to determine the best models for our case. Our approach shows that the best model demonstrates an 2 score above 90%, and = 0.048. 2.1
Given a miRNA group as a query, the biological function
permutation test introduced in [
Based on the previous, we can now describe one popular permutation test: Let be the set of genes containing the union of targets from all miRNAs in a query group, and also let the set of genes participating in a biological function. The statistical measure used to compare the query to the randomly assembled miRNA groups is the left-sided overlap and it is defined as follows: ( ∩ ) - (, ) = ()
Essentially, the left-sided overlap is defined as the proportion of targeted genes that also participate in the biological function. Then, we create 1 million random miRNA groups with the same size as the query and calculate the left-sided overlap for each of them. The empirical p-value is defined as (where overlap is the left-sided overlap):
({ : ( , ) ≥ (, )}) − = which is the proportion of randomly assembled groups presenting a larger left-sided overlap than the query. 2.2
The above analysis relies on a very large number of union and intersection set operations. Given that this analysis is performed for more than one biological functions, and that more than 20K biological functions exist, a few million union and about 20 billion intersection operations are performed in the span of the analysis.
The software implemented by the authors in [
Then in [
In order to train a machine learning model to predict p-values (p_value) as accurately as possible, we selected features from our dataset based on their biological meaning and their relevance to the analysis. The features are summarized below: miRNA group size (mirna_group_size): The number of miRNAs in a miRNA group.
that identifies each biological function. This string consists of 2 letters (same for all biological functions) and a numerical part. We turned this ID into a numerical value, by stripping the letters from the string.
Number of common genes (number_of_common_genes): The number of genes targeted by a miRNA group that also belong to the biological function.
lap as defined in Section 2.1.
from Section 2.1, the right-sided overlap is defined as: ( ∩ ) ℎ - (, ) = () Two-sided overlap (overlap): Given and from Section 2.1, the two-sided overlap (or Jaccard coeficient) is defined as: ( ∩ ) - (, ) = ( ∪ )
(common_genes_proportion_to_total): The number of common genes between M and G as the proportion of the total number of genes in the universe of genes.
Common gene list (common_genes): The list of common genes, sorted by alphabetical order. We used label encoding to turn the values into categorical values.
chromosomes on which the common genes are located, sorted by alphabetical order. We used label encoding to turn the values into categorical values.
(number_of_common_chr): The number of chromosomes on which the common genes are located.
number of chromosomes on which the common genes are located divided by the number of chromosomes of the genes targeted by the miRNA group.
The number of chromosomes on which the common genes are located, divided by the number of chromosomes of the genes belonging in the biological function.
of chromosomes on which the common genes are located, divided by the number of chromosomes for the union of the genes contained in and .
Using these features, we produced a dataset for groups of
miRNAs of size 10, 25, 50, 100 (containing 100 groups of each size)
for all biological functions as described in the Gene Ontology [
In order to find the best method to use for our case, we are going to use the following algorithms: • Linear/Ridge/Lasso Regressors: traditional regression methods, fitting a linear equation that uses the least squares method (with several variants of regularization). • Decision Tree Regressor: uses decision trees to make predictions. • Random Forest Regressor: estimates target value by combining average estimation values of several individual prediction models based on classifying decision trees for a number of subsets of the data set. • Adaboost Regressor: combines multiple weak decision trees into one. • Gradient Boost (XGBoost, LightGBM, CatBoost): in XGBoost, the estimation of the target is done by combining estimates of many individual prediction models based on decision trees. LightGBM is similar to XGBoost, but prediction is much faster than XGBoost. Regarding CatBoost, the key diference is that it builds decision symmetric trees. Both LightGBM and CatBoost can inherently handle categorical features. • MLPerceptron Regressor: typical Neural Net configuration. 4
In this section, we present preliminary results on a dataset that has been created by only using 1000 random miRNA groups to calculate p-values. This was done due to time constraints and to reach some preliminary conclusions about the dataset, in order to compare the algorithms presented in Section 3.2. Given the definition of an empirical p-value in Section 2.1, we expect that the p-values in the dataset will present a small accuracy, since the number of random miRNA groups is three orders of magnitude smaller than the required number. This is expected to create larger errors than the ones that the permutation test produces. However, these preliminary results are a first indicator on which algorithms present the worst performance in terms of accuracy and thus, which algorithms will be featured in the final work. 4.1
In Figure 1 we can see a heat map of the linear correlation between the features. Each cell represents the correlation measured by Pearson correlation coeficient between features in row x and column y. The values of range between -1 and 1. The larger the deviation from 0 the stronger the positive (for values closer to 1) or negative (for values closer to -1) correlation is. It is interesting to note here that the p-value, which is the feature we want to predict, does not seem to have a strong linear correlation (either positive or negative) with any of the features individually. Nevertheless, this does not mean that the p-value is not related with the features in a non-linear way. More specifically, it is expected to be related with the left-sided overlap (see p-value definition in Section 2.1). 4.2
In this section we are going to perform a preliminary evaluation
of the algorithms we outlined in Section 3.2. In order to evaluate
and compare the algorithms, we are going to use the following
statistical measures [
Model selection and implementation was performed with the method of k-Fold cross-validation. We first split the dataset into the training and testing dataset. The training dataset is split into k folds (groups). Then k-1 of those groups are used to train the model and the group left is used for validation. The score of the validation is recorded. Then, the process is repeated, but this time another group is used for validation and the rest for retraining, until all folds have been exhausted. The average of all validation scores is the final validation score. Finally, the test score is retrieved by testing the model against the testing data set. Final validation scores for each algorithm are shown in Table 1 and the testing scores in Table 2. Model parameters were selected via Grid Search for the Linear Regressors and the rest via Random Search, due to the large number of parameters and memory constraints.
It should be noted that the algorithms based on linear regression models do not perform as well as the other algorithms. Our approach shows that the best model is LightGBM, demonstrating a = 0.048. 5
In this paper, we presented an approach for miRNA enrichment analysis using machine learning techniques, in order to predict p-values instead of calculating them using randomization experiments. The goal is to simplify the work for bioinformatics data analysts, facilitating multiple enrichment analysis. Preliminary results showed that the best model demonstrates an 2 score above 90%, and = 0.048. As next steps, we plan to expand our dataset to include p-values estimated using 1 million random groups in order increase their accuracy as well as the accuracy of the prediction.
We acknowledge support of this work by the project “Moving from Big Data Management to Data Science" (MIS 5002437/3) which is implemented under the Action“Reinforcement of the Research and Innovation Infrastructure”, funded by the Operational Programme “Competitiveness, Entrepreneurship and Innovation” (NSRF 2014-2020) and co-financed by Greece and the European Union (European Regional Development Fund). 2 2 N/A N/A
N/A Management (SSDBM 2020). Association for Computing Machinery (ACM).