Discovering MicroRNA-Regulatory Modules in Multi-Dimensional Cancer Genomic Data: A Survey of Computational Methods

Sequence-Based miR–Target Prediction.

A number of bioinformatics algorithms for predicting miR recognition sites within transcripts have been developed using knowledge from experimentally validated target sites. Early experimental studies found that a primary determinant of target specificity was perfect complementarity (canonical site) at the 5′ end of the miR “seed region” at positions 2–7.^1,2,62 Therefore, these algorithms initially were primarily focused on sequence complementarity between the seed region of miR and the 3′ untranslated region (UTR) of the putative target.^62,63 However, given the large number of randomly occurring six nucleotide sequences in a 3′ UTR of a gene, perfect seed match itself is a poor predictor of miR regulation.^64–66

Several studies have proposed that target sites in which the pairing between miR seed and mRNA does not completely match (termed non-canonical sites) also confer regulatory effect.^67,68 A recent study by Agarwal et al. ⁶⁹ found that while miRs bind to non-canonical sites, there was no detectable repression based on mRNA stability or translation using multiple cell types. This finding has supported the focus on canonical binding sites by sequence-based target prediction programs. Other algorithms consider additional data, including mRNA secondary structure and target accessibility^70,71; however, this also results in a large number of predictions with many false positives (FPs). ² Many algorithms, including TargetScan ⁶⁹ and TargetRank, ⁶² use evolutionary conservation of the target site to select predicted targets based on conservation to reduce FP predictions. However, approximately 20% of functional target sites are not conserved between mammals, and conservation is further decreased in a step-wise manner in larger taxonomic groups, ⁷² indicating that sensitivity of target prediction decreases with higher conservation thresholds. ⁶⁵

The growing number of experimental MTI data has prompted more recent use of machine learning (ML) algorithms to train classifiers directly on the experimental data. For example, miRSVR exploits mRNA expression data from miR transfection experiments to train an ML algorithm to predict MTIs. ⁶⁴ The miRSVR scoring model for predicted MTIs is calibrated to correspond linearly with the probability of downregulation of the target, providing a meaningful guide for selecting a score cutoff. Agarwal et al. ⁶⁹ recently compared 17 sequence-based target prediction algorithms and found that the number of potential MTIs varied greatly, reflecting the varied strategies of these algorithms. Their analysis found that TargetScan version 7 performed significantly better than the existing models and was as good as recent high-throughput experimental approaches to identify effective target sites. Table 1 summarizes the target prediction algorithms discussed earlier. Sequence-based miR target prediction algorithms have been comprehensively reviewed elsewhere. ⁷³

OPEN IN VIEWER

Table 1.

List of software, websites, and references to methods for sequence-based miR target prediction and method for inferring miR–target relationships using paired expression profiles of miRs and genes in single-cancer datasets.

SEQUENCE BASED METHODS FOR miR TARGET PREDICTION (PREDICTING WHETHER A GIVEN mRNA IS TARGETED BY A miR)
METHOD/REFERENCE/SOFTWARE	DATA TYPES	COMMENTS
TargetRank Nielsen et al. 62 http://hollywood.mit.edu/targetrank/	*mRNA sequence	*scoring system based on sequence complementarity and conservation
miRanda Enright et al. 114	*mRNA sequence	*sequence complementarity and estimated minimum free energy *source code in C freely available
TargetScan (Version 7) Agarwal et al. 69 http://www.targetscan.org/vert_71/	*mRNA sequence	*scoring system based on 14 features found to be informative of binding efficacy using a regression model
STarMir Rennie et al. 70 http://sfold.wadsworth.org/cgi-bin/starmir.pl	*mRNA sequence	*model for binding site predictions trained on mir binding sites from CLIP data
miRanda-miRSVR Betel et al. 64 http://www.microrna.org/microrna/home.do	*mRNA sequence	*regression model trained on miranda predicted target site features and mir transfection data to predict target site binding efficacy

METHODS FOR INFERRING miR–TARGET RELATIONSHIPS USING PAIRED MIR AND GENE EXPRESSION PROFILES IN SINGLE-CANCER DATASETS
METHOD	DATA TYPES	COMMENTS
Correlation coefficient based methods Peng et al. 47	*miR & gene expression *sequence predicted targets	*Proposed permutation based method to estimate FDR of MTIs
MLR Yang et al. 17	*miR & gene expression *sequence predicted targets *CNA *PM	*models gene expression by a linear combination of all miR expression profiles (adjusting for epigenetic and genomic effects)
LASSO Lu et al. 52	*miR & gene expression	*models gene expression given multiple potentially competing miRs
Elastic net regression Sass et al. 82	*miR & gene expression	*found superior performance for identification of experimentally validated MTIs versus Lasso and PCC
Causal inference (IDA) Le et al. 88	*miR & gene expression	*ensemble of LASSO, PCC, and IDA detected more MTIs than any single method
Maximal information content (MIC) Le et al. 88	*miR & gene expression	*mutual information based method to detect linear and non-linear associations between two variables
GenmiR++ Huang et al. 51	*miR & gene expression *sequence predicted targets	*Bayesian inference method for MTI prediction

Abbreviations: FDR, false discovery rate; MTI, miR–target interactions; UTR, untranslated region; LASSO, least absolute shrinkage and selection operator; PM, DNa promoter methylation; CNA, copy number abnormalities; PCC, Pearson's correlation coefficient; IDA, interventional calculus when the directed acyclic graph is absent; CLIP, crosslinking and immunoprecipitation.

High-Throughput Experimental Methods for miR–Target Identification.

High-quality experimentally derived training data are generally required to improve sequence-based target prediction performance. ⁷⁴ High-throughput methods such as those employing crosslinking and immunoprecipitation (CLIP) are an important class of capture-based methods for detection of direct miR–target binding events associated with the Argonaute protein (Ago).^75,76 Argonaute high-throughput sequencing of RNAs isolated by CLIP simultaneously sequences Ago-miR and Ago-mRNA binding sites to identify interaction sites between miR–target pairs. ⁷⁶ One limitation of this approach is that miR–target complexes are dissociated prior to sequencing, requiring the target sequence in each miR–target pair to be inferred computationally, which is prone to error.

Recently a method for producing ligation of the miR–target pair called crosslinking, ligation, and sequencing of hybrids (CLASH) has shown to be more robust than CLIP for identification of miR target sites. ⁶⁷ The former method is similar to CLIP, but adds a ligation step between the miR and target, allowing direct characterization of the chimeric or hybrid miR–target to unambiguously identify the miR bound at a specific target site. A novel finding from the CLASH analysis was the detection of strongly overrepresented motifs in the interaction sites of several miRNAs, suggesting that individual miRs systematically differ in their binding site modes. Although this likely affects the response of RISC to miR–target binding, it is unclear how it impacts in vivo function of MTIs. ⁶⁷

While CLASH holds much promise, at the present time, this method has a very low yield with only ~2% of the reads obtained in an experiment corresponding to miR–target chimeras. Thus, further improvements to CLASH will be needed before comprehensive mapping of MTIs will be possible. ⁷⁷ As each cell line has a different miR expression profile, the cell line used in an experimental analysis will yield different sets of MTIs than other cell lines or disease conditions. For example, an miR with low expression in a cancer (sub)type profiled using CLASH may not be detected, whereas this miR may have high expression in another (sub) type of cancer. Therefore, studies have integrated knowledge of MTIs from one cell line or condition identified by CLIP and CLASH with miR and gene expression profiles from the cell line or condition of interest, yielding condition-specific MTIs. For example, StarBase, a database of 108 CLIP- and CLASH-based datasets from 37 studies, integrates TCGA pan-cancer expression data with CLIP-seq data to provide MTI pan-cancer predictions. ⁷⁸ Recent comprehensive reviews of experimental target prediction methods are available elsewhere.^75,77

While high-throughput experimental (in vitro) approaches are being increasingly performed, they currently suffer from several disadvantages compared to in silico predicted MTIs (based on sequence based and/or structural stability): 1) the number of experimentally derived MTIs remains far fewer than the predicted MTIs ⁷⁹ and 2) the tissue specificity of the experimental method may exclude availability of cancer (sub) types of interest to the investigator. Therefore, the most common approach to facilitate inference of cancer-specific MTIs involves profiling paired miR and gene expression datasets. Integrating condition-specific (dynamic) expression profiles with in silico predicted MTI datasets (static) datasets (as discussed in the next section) has been shown to improve MTI prediction. ⁸⁰ Several recent studies have demonstrated better MTI prediction in cancer datasets by integrating CLIP data with sequence-based target predictions.^11,81 In the next section, we discuss commonly applied methods to integrate paired miR and gene expression with in silico and in vitro data types for multi-dimensional analysis to identify significant MTIs.

Correlation Coefficient-Based Methods.

The principle of assuming that the expression levels of miRs and target mRNAs are negatively correlated is commonly used to detect MTIs.^47,48,84 These methods typically select potential miR–target pairs that (i) are negatively correlated above some statistical significance threshold and (ii) have been identified to interact using sequence-based target prediction or experimental methods (Fig. 2). The large number of miR–target correlations calculated necessitates estimation of the false discovery rate (FDR), defined as the number of FP divided by the number of FP and true positives (TPs). Peng et al. ⁴⁷ proposed a permutation-based method to estimate the FDR of miR–target correlations at a given statistical threshold (Fig. 3A). The FDR was defined as the ratio of the number of correlated miR–mRNA pairs above a given threshold (eg, Pearson's r < −0.5 and P < 0.01) in a randomly permuted dataset (ie, FP) to the number of pairs above the threshold in the original dataset (ie, TP and FP). To generate the randomly permutated miR–mRNA datasets, the sample labels for miR and mRNA were randomly swapped such that the samples in the random miR datasets did not correspond to the samples in the random mRNA datasets. This process was repeated 100 times, and a median value of FDR was selected. Using this approach, the authors were able to ascertain that Pearson's r < −0.55 with P-value < 0.01 was associated with an approximately 5% FDR. An alternative method to estimate FDR using empirical Bayes (EB) is discussed in the “Bayesian models” section (Fig. 3B).

OPEN IN VIEWER

Figure 2.

Inferring MTIs by integrating match miR–gene expression profiles and sequence-based target prediction data. Sequence predicted targets from a pre-selected database (orange matrix) depicted as a binary matrix (indicating the presence or absence of miR–target pairs, as dark orange or light orange boxes, respectively). Expression profiles from matched gene (green matrix) and miR (red matrix) microarrays are correlated using the PCC and input into the purple matrix. MTIs with PCC above a selected threshold and present in the sequence-based target prediction database are indicated as dark purple boxes, while all other pairs as light purple boxes.

OPEN IN VIEWER

Figure 3.

Methods for estimation of FDR from MTI data. (A) Permutation-based method to estimate FDR at a given statistical threshold. Randomly permuted datasets are generated to estimate the proportion of FPs above a given threshold. To generate randomly permuted datasets, the sample labels for miR and gene expression matrices are randomly swapped such that samples in the random miR datasets do not correspond to samples in the random gene datasets (swapping of sample labels indicated by red arching arrows above miR and above gene expression matrices). PCC is calculated for each miR–gene pair in the sequence-based target prediction matrix (dark orange boxes above). The distribution of PCC values is depicted as a gray histogram (below). PCC values below a threshold of −0.4 (vertical dashed red line) represent FP in the permuted data (gray). PCC values below this threshold using the original (unpermuted dataset) represent TPs (yellow region in the histogram). The FDR is calculated as the number of FP divided by the number of FP and TP. (B) EB estimation of FDR using local FDR. PCC values for all miR–gene pairs within the sequence-based target prediction matrix are calculated and transformed into z-scores. The z-scores corresponding to false (null) interactions are depicted as a histogram in gray (below), and TPs are shown in yellow within the histogram. The empirical null distribution (orange curve) and theoretical null distribution (blue curve) are shown overlapping the histogram. Note that the theoretical null is too narrow for the data, whereas the empirical null (determined using local FDR) is a better fit for the data.

For miR and gene expression profiles from heterogeneous conditions such as multiple (sub)types of cancer samples or cancer and healthy controls, correlations are often calculated for each condition separately rather than over the data as a whole.^81,85,86 The rationale for this approach is that each cancer type may have differences in dysregulation of MTIs. Farazi et al. ⁸¹ found significant differences in the medians of the correlation distributions for miRs and their predicted targets comparing between three subtypes of breast cancer. They also detected distinct correlations between the expression of specific miR families and their predicted targets among the three cancer subtypes, providing a rationale for performing separate MTI analysis based on the subtype. Importantly, the top-ranked miRs according to regulatory activity did not necessarily overlap with the top-ranked miRs according to association with tumor subphenotype. This disconnection becomes especially relevant for validation and development of future miR-based therapies.

MLR and Regularized Least-Squares Models.

MLR models. Whereas correlation-based approaches consider only the pairwise miR–gene expression, MLR models gene expression by a linear combination of all miR expression profiles targeting the gene. Furthermore, other epigenetic or genomic factors can also be modeled such that gene expression is the response (dependent) variable and the transcriptional and epigenetic regulators are the independent variables in the models (Fig. 4). Integrating gene expression with associated alterations in genomic, epigenetic, and miR expression has been undertaken in several studies to identify molecular drivers of cancer (sub) types. Yang et al. ¹⁷ applied an MLR model to analyze gene expression of each gene in a mesenchymal signature of ovarian cancer. The putative regulatory factors of each gene selected for analysis were the associated DNA copy number abnormalities (CNA), promoter methylation (PM), and miR expression level, which were used as independent variables in the regression model. Based on this analysis, the investigators detected a set of 219 genes predicted to be targeted by 19 miRs in an miR–mRNA network. These genes could be used to distinguish the mesenchymal subgroup of ovarian cancer from other ovarian cancer subtypes.

OPEN IN VIEWER

Figure 4.

Factors controlling gene expression in cancer. The expression of the oncogene MYC is regulated by miRs gene expression (pink square, left), the availability of the miR-RISC (shown in blue), DNA-methylation promoter region (purple hexagon), and gene copy number abnormalities (duplicated copies of MYC shown as purple squares). MiR loading to RISC leads to targeting the 3′ UTR of the MYC mRNA, where it is degraded and translation by the ribosome is inhibited.

Regularized Least Squares.

Least absolute shrinkage and selection operator. The least absolute shrinkage and selection operator (LASSO) algorithm performs variable selection procedure by estimating linear regression coefficients through L1-constrained least squares. This L1 constraint (penalty) tends to result in regression coefficients that are exactly zero, thereby imposing sparsity on feature selection, making the models more interpretable. Several investigators have argued that LASSO avoids overfitting in the presence of noisy expression data and a large number of explanatory variables, leading to better prediction accuracy.^53,87 Li et al. ⁵⁴ found LASSO with negative coefficient constraint as the best performing method for the joint analysis of miR and mRNA expression data compared to Pearson's correlation-based methods and Genmir++ (a Bayesian learning method to infer the probabilities of targets) ⁵¹ using a database of experimentally verified MTIs as a reference.

MLR with Elastic Net Penalty.

While LASSO models multiple potentially competing miRs, ⁵² several investigators^82,86 have noted that it fails to account for co-expression between miRs targeting the same gene. The sparse solution (L1 penalty) implemented in LASSO selects one representative miR from each correlated group of miRs, disregarding other potentially biologically relevant co-expressed miRs within the group. In contrast, the ridge regression model maintains predictors in the model by using an L2 penalty in which the coefficients are small but non-zero. However, this approach does not facilitate feature selection. The elastic net penalty overcomes the drawbacks of LASSO and ridge regression by combining L1 and L2 penalties in order to account for co-expression among miRs while simultaneously performing feature selection. Compared to Pearson's correlation and LASSO, elastic net regression in combination with negativity constraint coefficients was applied to a head and neck cancer dataset and was found to identify MTIs with greater enrichment for experimentally validated MTIs. ⁸²

Ensemble of Computational Methods.

Le et al. ⁸⁸ investigated the performance of eight methods for joint analysis of miR–target expression and prediction methods using datasets from multiple cancer types. The computational methods compared were (i) Pearson's correlation coefficient, (ii) MIC (maximal information coefficient, a mutual information based method to detect linear and non-linear associations between two variables), ⁸⁹ (iii) z-score (a network inference method to estimate the effects of gene knockout experiments), ⁹⁰ (iv) IDA (interventional calculus when the directed acyclic graph is absent, a causal inference method), ⁹¹ (v) ProMISe (probabilistic MiR–mRNA interaction signature, a method that estimates the probability of a gene to be a target of an miR by taking competition between genes and between miRs into account), ⁹² (vi) GenMiR++, ⁵¹ (vii) LASSO, and (viii) elastic net regression methods.

An ensemble of three methods (Pearson, IDA, and LASSO) obtained more targets than any single method and identified targets with enhanced functional enrichment. From this, we can infer that different classes of computational methods tend to identify unique sets of validated targets and therefore each has its own merits. ⁹¹ It is important to note that this analysis did not incorporate sequence-based target predictions. Thus, the investigators hypothesized that the poor performance by GenMiR++ and ProMISe, two methods originally designed to incorporate sequence-based predicted targets, may have been due to the use of miR and gene expression data alone. Le et al. ⁹³ provided a software pipeline integrating the above computational methods in the R package miRLAB that is freely available on Bioconductor at http://bioconductor.org/packages/miRLAB/.

The statistical approaches described in this section quantify the association between any given single miR and single gene (refer “Correlation coefficient-based methods” section) or associations between multiple (co-expressed) miRs and a single gene (refer “MLR and regularized least-squares models” section). These approaches are summarized in Table 1. In the next section, we explore methods for identifying associations between co-expressed miRs and groups of gene targets in a network of “many-to-many” miRs and genes termed an MRM. As is discussed in the following, a greater understanding of the (patho)biological roles of miRs and their targets can be gained by identifying MRM compared to single MTI analysis.

Bipartite Graphs and Biclique Enumeration Approach.

A bipartite graph (network) consists of two sets of nodes (in this case, miR and target) and a set of edges connecting the nodes (in this case, edges represent association strength between miR and target). ³⁶ A putative module (in this case, MRM) in the bipartite network has been postulated to correspond to a biclique, a special type of bipartite network where every node in the first set (miRs) is connected to every node in the second set (target genes). ⁴⁷ A biclique is called a maximal (complete) biclique if it is not contained in a larger biclique. To perform enumeration of the maximal bipartite cliques (EBC) within the bipartite network representing putative MRMs, the module input consensus algorithm (MICA), ⁹⁴ the most efficient algorithm for finding bicliques, is often used. ⁴⁷ Zhang et al. ⁹⁵ argued that MICA is designed for general graphs and unable to take advantage of the bipartite structure. Therefore, the investigators developed the maximal biclique enumeration algorithm, which outperformed MICA. ⁹⁵

Peng et al. ⁴⁷ developed one of the earliest approaches to identify MRM using a maximal biclique method on MTIs discovered in a multi-step approach. In the first step, pairwise Pearson's correlation between the differentially expressed (DE) miRs and genes across all samples was preformed to identify putative MTIs. The MTIs were selected if they exceeded a predefined FDR threshold (discussed in the “Correlation coefficient-based methods” section). Then the MTIs were further selected if present in a set of sequence-based target predictions, resulting in a binary matrix of MTIs. This matrix of selected MTIs was graphically represented as a bipartite network. Using MICA, maximal bicliques (ie, putative MRMs) were identified, which comprised between one and three miRs per MRM.

Maximal biclique-based methods to discover putative MRM have been argued to be too sensitive to noise in the data and frequently produce MRM with a high level of redundancy and only a single miR that cannot be used to explore miR combinatorial regulation.^95,96 Missing subsets (false negatives) or erroneous (FP) MTIs may have an adverse effect on the quality of the maximal bicliques detected. Furthermore, searching only for maximal bicliques may be too restrictive as they are defined by an all-to-all relation between miRs and targets within the MRM. ⁹⁷ To add flexibility to MRM detection, several studies^98,99 have applied quasi-bicliques, which exhibit a most-to-most interaction between miRs and genes within the MRM (Fig. 5). ⁹⁷ Quasi-bicliques allow all nodes in the bipartite network (miRs and genes) to accommodate missing interactions up to some user-determined level. Veksler-Lublinsky et al. ⁹⁹ found that MRMs discovered using the quasi-biclique method more significantly identified MTIs than a maximal biclique approach.

OPEN IN VIEWER

Figure 5.

Comparison of maximal biclique and quasi-biclique enumeration of bipartite graph for identification of MRMs. A matrix of miRs and genes (purple) indicating known MTIs (dark purple boxes) is represented as a bipartite graph with miRs shown as circles on the left (labeled 1–5) and genes shown on the right (labeled 1–8) in a) and b). MTIs are depicted with an edge (straight line) connecting the miR and gene. (A) Maximal biclique enumeration (all-to-all interactions): miRs and genes within a maximal biclique (representing an MRM) are shown in blue, whereas those excluded from a maximal biclique are shown in gray. Unlike miR-1 and miR-2, miR-3 is not connected to all genes within the biclique and is therefore excluded. Likewise, miR 4 and miR 5 do not form a maximal biclique and are excluded. (B) Quasi-biclique enumeration (most-to-most interactions): quasi-biclique enumeration allows some missing interactions within the biclique; therefore, miR-3 is included with miR-1 and miR-2 (blue, above); miR-4 and miR-5 form a biclique (red, below). Thus, in this toy example, quasi-biclique method is more sensitive to detecting MRM than maximal biclique enumeration.

Kim et al. ¹⁰⁰ have argued that the bi-relationships modeled between miRs and targets using bicliques are unsuitable for complex genetic interactions because information is lost. Instead, they applied hypergraphs to generalize the concept of an edge between nodes to a hyperedge by which more than two variables could be connected simultaneously. As the weight of a hyperedge reflects the strength of higher order dependency among variables, it was hypothesized that each hyperedge potentially behaves as a gene module. To model prostate cancer stage-specific MRM networks, the investigators developed a hypergraph model with each hyperedge represented by miR and gene expression for each cancer stage. Using a learning hypergraph model, hyperedges having high discriminative capacity between cancer stages were selected The investigators identified putative prostate cancer stage-specific MRMs; however, it is unclear whether the hypergraph structure improves MRM discovery over the maximal biclique approach utilized by Peng et al. ⁴⁷

Matrix Factorization Approach.

An alternative strategy to the biclique approach was implemented by Zhang et al. ⁹⁶ by integrating miR and target expression profiles using non-negative matrix factorization (NMF) to identify putative MRM. The NMF method decomposes a matrix to find two smaller (lower rank) non-negative matrices, allowing substructures to be readily identified within the data. This method is similar to principal component analysis, another unsupervised matrix decomposition technique, except that it employs a constraint of non-negativity instead of orthogonality. ¹⁰¹ Zhang et al. ⁹⁶ extended the NMF method by simultaneously factoring two variable types (expression profiles for miRs and genes) into a common basis matrix and two lower rank matrices. The decomposed matrices were then used to identify MRM with each miR or gene permitted to be assigned to multiple MRMs. Sparsity was induced via L1 penalty, and they named their algorithm the sparse network-regularized multiple NMF (SNMNMF) technique. To guide the optimization process, the investigators incorporated sequence-based predicted target data and GGI using a semi-supervised learning framework to define constraints for MRM. Using a multiplicative interactive procedure, the model is learned until convergence to a local optimum.

Comparing NMF to EBC methods using the TCGA ovarian cancer data, only 19% of the modules identified by EBC were enriched for at least one gene ontology (GO) biological term compared to 53% of the modules identified by NMF, indicating that this method performed better than EBC. The limitations of the NMF approach include the requirement for a pre-defined number of modules in order to perform the matrix factorization (which may be difficult to estimate empirically), the solution is often not unique, and it does not enforce a negative correlation constraint. A negativity constraint is implemented by other MRM methods to account for the fact that since miRs mainly repress mRNA target expression levels, positively correlated MTIs are likely to be FP.^3,82

An additional limitation of the conventional NMF method is that it only permits two types of genomic data (eg, miR and gene expression) to be analyzed. Therefore, Zhang et al. ⁴⁹ implemented a joint NMF (jNMF) semi-supervised framework for the integrative analysis of more than two types of genomic data to identify groups of methylomic markers, miRs, and genes in putative modules in ovarian cancer datasets. Yang and Michailidis ¹⁰¹ have noted that the jNMF method is not different from NMF, and as a consequence, jNMF does not distinguish between different data sources in the integrative analysis, which can be problematic for analyzing heterogeneous data. Thus, the investigators developed an integrative NMF (iNMF) method that applied a novel tuning selection procedure that allows the model to adapt to the level of heterogeneity among the datasets. The iNMF method was found to be more robust to heterogeneous noise across the data sources than jNMF for the detection of true modules. These matrix factorization approaches exemplify transformation-based integration of multi-dimensional datatypes (refer “Framework for elucidating MiR-regulatory modules (MRM) from multi-dimensional omics data” section).

BN Models.

A BN is a graphical model based on probabilistic relationships among the measured variables, ¹⁰² which can be applied to miR and mRNA expression datasets to discover MRM.^103,104 Inference of the BN structure (topology) involves searching among the possible networks and then scoring these structures. Two nodes are expected to be connected to each other if one node (ie, miR expression) affects another (eg, gene target expression). If the search space (number of potentially interacting variables) is not sufficiently restricted, the process of learning a BN can become very computationally time consuming, as all possible networks will be formed. ¹⁰³ Therefore, most BN approaches rely on constraints to decrease the search space to provide a more compact representation of probability distributions. Prior information such as predicted target data can be incorporated to improve network construction.

After the topology of the BN is determined, the strength of the relationship between any two nodes is quantified using conditional probabilities (network parameters). Thus, the expression values (activities) of miRs and genes are represented as nodes, and the dependence between nodes (regulatory interactions) is represented by the edges within the network.

Jin and Lee ¹⁰⁵ proposed a method to integrate miR–gene expression data with GGI data using the TCGA ovarian cancer dataset as an alternative to the NMF method developed by Zhang et al. ⁴⁹ To constrain the search space, DE genes along with labeled samples were input into a biclustering algorithm. Clustering has been previously demonstrated to aid in the inference of BN by reducing the parameter space and avoiding highly correlated gene profiles from inhibiting interaction inference. ¹⁰⁶ The interaction of sample and gene expression profiles was modeled as a bipartite graph and generated subgraphs, termed gene-sample modules (GSMs). Within each GSM, the expression trend was similar for most genes among the subset of selected samples.

Since previous studies have shown that not all the genes in cancer-related pathways undergo expression changes, the investigators expanded the GSMs using GGI data by including genes that are significantly correlated with at least one gene in a GSM. Next, to discover MRMs, a BN model was used to estimate dependencies between expression values of miRs and genes in the GSM using Bayes information criterion (BIC). The BIC is a score function used to assess the degree to which the BN explains the data (ie, whether it provides a “good” structure). ¹⁰² The search space is constrained to the miRs whose correlation coefficient values for genes in a given GSM are in the top T%. The investigators found that the average number of enriched pathways in modules using this method was larger than that in the NMF approach by Zhang et al. ⁹⁶ described in the “Matrix factorization approach” section.

Liu et al. ¹⁰⁴ were the first to model separate BN structures in two different conditions (eg, cancer and normal) in order to effectively identify both strong and subtle MTIs. For each condition, the BN learning was performed using the miR and gene expression as input, initially generating separate BN models for each condition. The BN learning iteratively evaluated the MTIs, removing MTIs from the initial structure using the expression data and retaining only high-confidence MTIs, with the goal of selecting the structure that best fits the data based on a scoring function. To avoid statistically insignificant results and overfitting with small sample sizes, a bootstrapping step (sampling with replacement) was added to achieve reliable inference and integration of the BN models from the two conditions. To significantly reduce the search space, the BN was assumed to have a bipartite graph topology, and sequence-based predicted targets were further used to constraint the initial network. Furthermore, miR and gene expression profiles were discretized into a binary status: “up” or “down” regulated. This condition known as “splitting” approach, initially used to generate separate models, was shown to capture complex MTIs from the cancer and control samples. Using datasets from epithelial and mesenchymal cell lines, their method named Bayesian networks splitting and averaging was found to identify more co-regulated targets by multiple miRs compared to conventional BN that did not separate the conditions. This example of intermediate integration ⁴⁰ learned after combining the sample types was also found to result in improved performance over concatenation-based integration by Gevaert et al. ¹⁰⁷ for integration of clinical and gene expression data.

Zacher et al. ³³ formulated a BN integrating miR and gene expression to infer miR activities in a specific condition versus a reference condition (eg, controls). In their model, if the activity of certain miRs changes between conditions, a shift in the expression value of the targets is expected. For example, an miR is considered to be active in a condition if it is upregulated and its gene targets are downregulated compared to the reference condition. While TF-target activity can also be inferred using this approach, we will focus on the use of this application for miR–target discovery. The investigators applied a Bayesian linear regression, where the expression level of each gene in each condition was determined by a linear combination of miR activities. Then a score was used to rank the degree to which miRs explained the observed differential gene expression, where a score close to 1 indicates that the corresponding miR is essential (ie, high miR activity) for explaining the differential gene expression in the condition versus reference. The miR activity states are inferred using Markov Chain Monte Carlo sampling. The model swaps the activity states of any two miRs with opposite activity states sampled from the posterior distributions. To reduce the size of the potential network, the investigators (i) input only DE miRs and genes between the two conditions, (ii) limited MTIs to those present in experimental or sequence-based predicted MTIs datasets, and (iii) required that the target genes for two miRs exhibit a minimal overlapping similarity (ie, 0.8). While this method identifies miRs (and their targets) with condition-specific activity, it does not directly identify MRMs in contrast to the other BN methods described earlier.

The investigators applied their method to ovarian cancer data from TCGA. They separated patients into “early” or “late” relapse (>1 year) and discovered 12 miRs active in the “early” but not in the “late” relapse patients, with target genes of 11 of the 12 miRs directly associated with relapse-free survival times. This method named Bayesian influence of regulation of transcriptional activity is freely available through the R package birta on Bioconductor at https://www.bioconductor.org/packages/release/bioc/html/birta.html.

These BN methods represent examples of transformation-based (intermediate) integration.^40,57 While the methods described earlier facilitate MTI or MRM discovery using datasets from a single-cancer type (or subtype), other recent works have investigated methods for inferring MTIs and MRMs across many types of cancers, which are discussed in the next section.

MLR and Regularized Least-Squares Models.

Multiple linear regression. Jacobsen et al. ⁸ analyzed associations between miR and gene expression for individual cancer types using an MLR model that was fit to the gene expression, taking into account biases due to copy number abnormalities (CNA) and PM at the gene locus (Fig. 6). The linear coefficients for each miR in the MLR model indicate the magnitude of association between the miR–target pairs for each cancer type, where the greater (negative) value of the coefficient infers greater miR activity on the target gene. This method more accurately evaluated miR–target expression association in the presence of CNA and PM abnormalities that influence gene expression. To evaluate the relative strength of association between given miR and gene in at least half of the cancer types while avoiding potential biases in cancer datasets (such as sample size), the authors developed a rank-based, statistic, association RECurrence (REC), under the null hypothesis that no association exists between the miR–target pair in any of the cancer types. Jacobsen et al. ⁸ applied a model-based (late integration) approach^40,57 (refer “Framework for elucidating MiR-regulatory modules from multi-dimensional omics data” section) that applies an MLR model to each cancer type individually and combines the results (ranks) into a pan-cancer REC score.

OPEN IN VIEWER

Figure 6.

MLR of mRNA expression in a single-cancer type to identify miR–mRNA associations. mRNA expression from each sample is modeled as the response variable (left), and miR expression, DNA copy number abnormalities (CNA), and level of PM are the independent variables. The MLR evaluates pairwise associations between mRNA j and miR k while controlling for the confounding variables CNA and PM. The greater the negative value of β _jk in the MLR, the greater the miR–mRNA association. The degrees of PM and mRNA and miR expression are shown as a heat plot (with increased, unchanged, and decreased levels shown in red, gray, and blue, respectively). CNA are expressed as an increase in the number of copies (gain), no change in copy number (wild type), or loss of copies of the gene, shown in red, gray, and blue, respectively. Modified from Jacobsen et al. ⁸

Using this approach, the investigators identified 143 miR–target pairs (40 miRs and 72 mRNA targets) having both strong negative REC score and exceeding sequence-based predicted target thresholds (using TargetScan and miRanda scores). The 40 miRs in the pan-cancer network were more likely to be dysregulated by genetic and epigenetic alterations – a common property of cancer driver genes – than 180 other miRs expressed across all cancer types. The results of this analysis are available in the CancerMiner website (http://cancerminer.org).

Regularized Regression (LASSO).

Setty et al. ⁸⁷ applied a LASSO model using CNA, miR binding site counts in mRNA 3′ UTR, and TF binding site counts in the gene promoter as covariates to explain gene expression changes (tumor versus normal) from subtypes of glioblastoma multiforme. First, LASSO was applied on a sample-by-sample basis to identify key direct regulators (miR and TF) that account for DE genes in glioblastoma multiforme relative to normal brain samples. A second LASSO model used a multi-task learning framework (group LASSO) to learn regression models of all samples simultaneously. Multi-task learning seeks to perform individual tasks (eg, selecting MTIs in subtypes), while exploiting the relationships between several learning tasks to increase the power of the search, thereby improving model performance. Group LASSO, an extension of LASSO, was developed for multi-task learning to select a common subset of features among tasks. ¹⁰⁹ For both LASSO models, Setty et al. ⁸⁷ used a feature-scoring scheme to determine significant regulators for common and subtype-specific gene expression. This scoring scheme estimated increase in total residual error in the models when an miR was excluded in order to determine its degree of influence in predicting gene expression changes. Using this method, common and subtype-specific miR and TFs were identified. This software is available in the R package RegulatorInference (https://bitbucket.org/leslielab/regulatorinference/).

Le and Bar-Joseph ⁵³ developed a probabilistic regression model called protein interaction-based miR modules (PIMiM) to discover MRMs in specific cancers (conditions) and across multiple cancers. The goal of this approach was to determine MRMs that explain gene expression as a function of the miR expression and the set of proteins the genes interact with. This module-based method assigned miRs and predicted targets to one of K modules (where K was a predetermined number). A regularized probabilistic regression model was learned in which gene expression data were regressed to the expression data of the miRs in the assigned modules. Their algorithm used sequence-based target prediction and GGI data as constraints. The cancer-specific approach used an L1-norm to encourage sparser modules and additional terms to reward sequence-predicted MTIs and GGI pairs to the same network.

PMiM was also found to detect MRMs with higher functional enrichment than the matrix factorization method by Zhang et al. ⁹⁶ when PIMiM was applied to the ovarian cancer dataset from TCGA. The investigators also found that using the prior knowledge from GGI datasets greatly increased precision and recall of 115 miRs known to participate in ovarian cancer. A potential disadvantage of supervised MRM identification methods such as PMiM is that the MRMs identified naturally tend toward the reference (in this case, GGI data). Yang and Michaelidis ¹⁰¹ have argued that publically available GGI datasets are prone to accumulated errors and oversimplification, and therefore, supervised MRM methods depend on the accuracy of these databases. Indeed, a recent evaluation of Mirsynergy, ⁵⁴ a LASSO and clustering-based MRM algorithm that is also dependent on GGI data, was found to produce module structures that were highly dependent on initial clustering of miRs and the GGI data. ¹¹⁰ Moreover, supervised methods are less likely to select new candidates based on the existing data and therefore are at a disadvantage for novel discovery. ¹⁰¹

To integrate multiple types of cancers, PMiM uses a group LASSO to regularize the models over multiple cancer (sub)types using an L1/L2 penalty to encourage miRs and genes to be assigned to the same modules across cancer types. One limitation of this approach is that the total number of modules must be determined in advance and is fixed for each cancer (sub)type. The group selection approaches implemented by Setty et al. ⁸⁷ and Le and Bar-Joseph ⁵³ are examples of model-based (late) integration methods (refer “Framework for elucidating MiR-regulatory modules from multidimensional omics data” section) in which different models are learned simultaneously across cancer types, increasing the power to detect pan-cancer MTIs compared to traditional approaches that integrate parameters after each cancer has been separately modeled.

Bayesian Models.

Chen et al. ⁸⁶ demonstrated the advantage of joint analysis of multiple cancers using an EB method over a single-cancer analysis for identification of MTIs. The underlying goal of their method was to integrate MTI commonalities among cancers by explicitly borrowing information among the cancer types to increase the power of detecting TP and reducing FP MTIs, respectively. It was recognized that different cancer types have distinct sets of MTIs, and therefore, the investigators expected that only a fraction of the MTIs would be shared between cancers. Additionally, cancers that were more closely related (eg, ovarian and breast cancers) were expected to share a higher degree of common MTIs than distantly related cancers. Thus, similarity between two cancers was quantified by the fraction of shared MTIs.

In the first stage, Pearson's correlation coefficients (PCCs) for all miR–mRNA pairs were expressed as z-scores (derived from Fisher transformation), which approximately followed a normal distribution, where increasingly negative z-scores were more likely to represent MTIs. Next, to correct for multiple hypothesis testing, a variant of the Benjamini–Hochberg FDR, termed local FDR, was applied to the set of all z-scores. ¹¹¹ The local FDR method estimated the empirical null distribution (histogram of z-scores) using maximum likelihood, where the central peak of the distribution mainly consisted of null cases (non-interacting pairs) and the (negative) tail tended to contain non-null cases (interacting pairs with negative z-scores). Thus, the FDR could be determined at any given z-score threshold; a lower (absolute) z-score threshold resulted in a greater sensitivity at the cost of a higher FDR and vice-versa (Fig. 3B). Using samples from glioblastoma multiforme and ovarian cancer, the investigators estimated that ~10% of miR–mRNA pairs may interact.

In the next stage, the estimated similarity between cancers based on the fraction of shared MTIs was calculated. Then the estimated similarity among cancers and the estimated MTI probability in individual cancers (obtained from local FDR) were combined to derive the posterior marginal probability of MTIs between miR–target pairs. The prior probabilities that measure cancer similarity were estimated from the observed data using an iterative updating procedure. Unlike the local FDR approach for single-cancer analysis, the pan-cancer approach was shown to change the order of z-scores, thereby re-prioritizing candidate MTIs by sharing information between cancers. This method called joint analysis of multiple cancers for miR–gene interactions (MCMG) was shown to have a higher precision than Pearson's correlation-based and LASSO approaches for identification of MTIs in a dataset of multiple cancer types. The MCMG method is implemented in R and available for download (http://bioinformatics.med.yale.edu/group).

Li and Zhang ¹¹² developed a pan-cancer analysis method called pan-cancer miRNA–target associations (PanMiRa) that directly infers the posterior distribution of the pan-cancer MTIs and accounts for potential genomic confounders. Unlike the MCMG approach that predicts cancer-specific MTIs by borrowing information from other cancers, PanMiRa aims to infer recurrent MTIs across all cancers while taking into account biases due to CNA and PM. Similar to Jacobsen et al. ⁸ , PanMiRa applies an MLR model for individual cancers with gene expression as the response (dependent) variable and miR expression, CNA, and PM as the independent variables (Fig. 6). The coefficient in the MLR model for miR expression (denoting relationship between miR k gene i in cancer type d) is converted to a t-statistic, which is subsequently transformed into z-scores. To reduce FP, only the interactions with negative z-scores in at least 75% of the cancer types were retained. Next, the investigators exploited the empirical distribution approach of the local FDR method applied to the z-scores to estimate the joint posterior of the true interactions across the pan-cancer types via an EB algorithm. The miRs and targets involved in pan-cancer interactions detected using PanMiRa were significantly enriched for known oncomiRs and oncogenes, respectively. The investigators also found a significantly higher enrichment number of MTIs detected by CLASH when comparing PanMiRa to the recurrence association method by Jacobsen et al. ⁸ and randomly shuffled posteriors. The source code is freely available in R (http://www.cs.utoronto.ca/~yueli/PanMiRa.html).