Network based multifactorial modelling of miRNA-target interactions

Construction of miRNA:target network

Four separate networks were used throughout our study which either vary in size (Sample vs Real) or lack additional interaction factors (Sample vs Sample+, Real vs Real+), summarized in Table 1 and illustrated in Fig. 1A. For Sample network, arbitrary miRNA:target interactions were generated along with arbitrary expression values. The Sample+ network had additional information, on top of Sample network, regarding miRNA:target interaction (Table 1). miRNA:target pairs from miRTaRBase, expression values of miRNA, and genes from TCGA were fed into SPONGE for sparse partial correlation analysis (List et al., 2019). Resulting miRNA:gene pairs were used to construct the Real network. The Real+ network had additional Cross-linking, Ligation, and Sequencing of Hybrids (CLASH) and Covalent Ligation of Endogenous Argonaute-bound RNAs-Cross-linking Immunoprecipitation (CLEAR-CLiP) (Helwak et al., 2013; Moore et al., 2015) data on top of Real network. Additional parameters for the Real+ network were curated from literature (see Section 2 in Supplemental Materials and Methods). Complex network analysis and visualization could be performed on constructed networks by tidygraph R package (Pedersen, 2020). Moreover, the networks could be exported to Cytoscape by RCy3 package for visualization and further analysis (Gustavsen et al., 2019).

Triggering perturbation and subsequent calculations

In order to demonstrate calculation steps after a perturbation event, Sample network at Table 1 was used. Initially, the network is assumed in steady-state condition (Fig. 2A; Fig. S2) and needs at least one trigger for initiating calculations. The trigger can be a change in the expression level of one or more genes (Fig. 2B; Fig. S3). After the trigger, the network undergoes the iterative cycle of calculations at each of which distribution of miRNA in the local neighborhood is recalculated (Fig. 2C). Based on new miRNA distribution, the expression level of each node (i.e., ceRNA) is updated (Fig. S4). Due to the common targeted nodes, the change in one neighborhood spreads to other neighborhoods (Fig. 2D), consequently has potential to affect the whole network due to the “ripple effect”.

During calculations, the following assumptions were adopted; (1) Transcription and degradation rates of miRNAs are steady and equal. (2) All available miRNAs are recycled as in miRNA:ceRNA interactions, targets are degraded and miRNAs are unaffected. (3) ceRNA targets also have stable transcription and degradation rates, and these rates are equal.

The repression level (referred to as Repression Count in formulas) of a miRNA on the individual target (RC_ij) is calculated according to Eq. (1). Briefly, the expression level of a given miRNA, namely Exp(miRNA)_j, is distributed among targets proportional to the expression of an individual target, Exp(Target)_i, over total expression levels of targets in the local neighborhood of miRNA_j. If an mRNA is targeted by more than one miRNA, the cooperative repression of miRNAs on the target is calculated by summing the repression level of each miRNA (see Fig. 1B, total repression on Node 3, on the right panel).

Multifactorial calculations in miRNA:target network

Our model integrates multiple factors when calculating overall miRNA activity. We classified factors into two categories. Affinity factors determine binding interaction between miRNA and target and they alter the amount of miRNA sequestered to the target. Degradation factors dictate the degradation efficiency of sequestered miRNA on its target. In other words, affinity factors exert their influence before or during binding, degradation factors exert their influence after binding.

From the literature, binding free energy (Cao & Chen, 2012; Helwak et al., 2013), seed type (Werfel et al., 2017) and binding region (e.g., 5′UTR, 3′UTR, CDS) (Hausser et al., 2013; Helwak et al., 2013) data were retrieved and plugged into Real+ network. For Sample+ at Table 1 network we used arbitrary interaction values.

The normalized values of factors are used to determine affinity and degradation of miRNA on its each target (Fig. S5) using Eqs. (2)–(4). Affinity factor of a miRNA:target interaction (AF_ij) is merely multiplication of normalized energy E^′_ij with normalized seed type effect STE^′_ij (Eq. 2). Amount of miRNA per its target (RC_ij, see Fig. S5C) is calculated by taking affinity factor into account (Eq. 3). Since, not all miRNA:target binding events result in degradation of target, cooperative degradation of an mRNA (R_i) by targeting miRNAs is calculated by aggregating multiplication of sequestered miRNA (RC_ij) with degradation factors such as binding region (RE^′_ij, normalized region effect) (Eq. 4 and Figs. S5D–S5E).

Analysis of real breast cancer networks and determining perturbation efficiencies

Expression levels of miRNA and genes in tumor and normal tissues of 87 patients were retrieved from TCGA; https://www.cancer.gov/tcga. We used sparse partial correlation method, SPONGE, to filter out genes with weak association (List et al., 2019). Log transformed miRNA and gene expression values were used and genes with p-adj less than 0.2 were considered, as suggested by authors (List et al., 2017). For Real network, predicted miRNA:target interactions from miRTarBase (Chou et al., 2017) and for Real+ network, experimental miRNA:target dataset (i.e., CLEAR-CLIP and CLASH datasets) were used as binary matrix for correlation analysis (Helwak et al., 2013; Moore et al., 2015).

Our package has functions to calculate perturbed count and perturbation efficiency for each node after using that node as the trigger. The perturbed count is the number of affected nodes, perturbation efficiency is the average of percent change in the expression level of perturbed nodes. Figures 1B and 1E illustrates perturbation efficiency in a tiny network for two different gene perturbations.

Through the parallel processing, we evaluated the perturbation efficiencies in Real and Real+ networks of all nodes in breast cancer tissue samples (Tange, 2011). Throughout the study, we used three-fold increase in expression level as the trigger. Perturbing nodes that show high perturbation efficiency in tumor and normal tissues were subjected to enrichment analysis in KEGG or GO using DAVID (Huang, Sherman & Lempicki, 2009; Huang da, Sherman & Lempicki, 2009).

Calculating threshold value for significant perturbation

For a perturbation to be considered significant we calculated the minimum required number of affected nodes and the number of samples (e.g., patient and normal). The distribution of perturbation efficiencies calculated in simulated networks was fitted to finite mixture model (Trang et al., 2015). After simulating Real and Real+ networks, we calculated that a node causing perturbation in at least 10 samples with 78 perturbed nodes and nodes that perturbed at least 216 nodes in 18 samples were considered as significantly perturbing, respectively.

Basic perturbation calculations using Sample network

To assess the effects of expression level changes in ceRNA regulation based on miRNA and target abundance we constructed the Sample network given in Fig. 2 and Table 1. After an increase in the expression level of a gene (G2), expression values of other genes also changed due to the redistribution of miRNA among its targets. Step by step calculations in gene and miRNA levels were described in Fig. 2. After the increase of expression level of G2, the miRNA that is found in the same neighborhood (M1) became less repressive on its remaining targets (G1, G3, and G4) resulting in increased expression in G1, G3, and G4. The similar finding was shown in an earlier ceRNA hypothesis model (Ala et al., 2013). More importantly, expression levels of G5 and G6 changed even though they don’t share common miRNA with the initial trigger gene (G2) (Fig. 2D). Previous studies showed that if the abundance of a gene increases remaining targets were affected due to shared targets of multiple miRNAs (Lai, Wolkenhauer & Vera, 2016; Salmena et al., 2011; Tay, Rinn & Pandolfi, 2014). Accordingly, the effect of perturbation in a single gene could spread through the whole miRNA:target network. It’s important to note that the changes in gene expression levels will have more pronounced effect if miRNA:target ratio is high which was reported in previous findings (Arvey et al., 2010; Bosson, Zamudio & Sharp, 2014; Denzler et al., 2014).

Perturbation calculations in Sample+ network: effects of binding energy, seed type and region

Our approach could calculate miRNA repression activity more accurately, by integrating seed type, region, and energy parameters. To demonstrate this capability, we generated Sample+ network, summarized in Table 1 and Table S1. Expression values of Sample and Sample+ networks were identical however Sample+ contains additional interaction parameters (Table S1) because of which steady-state miRNA distribution in Sample+ (Fig. 3A; Fig. S6) was different than that of Sample network (Fig. 2A). For instance, proportional distribution of G1:M1 interaction in Fig. 2A was higher than the G1:M1 interaction in Sample+ network (Fig. 3A; Fig. S5C) (286 vs 6.58). When the expression of Gene2 (G2) increased (Fig. 3B; Fig. S7), expression values of all genes also changed at various levels because of the contribution of interaction parameters (Figs. 3B–3D; Figs. S7 and S8).

In Sample+ network, G4 was targeted by two miRNAs, M1 and M2, hence expression level changes in G4 have more prominent effect in the whole network when compared to G2 (Fig. S9). The fact that each node had different expression levels and topological features suggested each node had different perturbation efficiency. Our package could screen all nodes for their perturbation efficiencies. When Sample and Sample+ networks were screened for perturbation efficiencies, the results were drastically different (Table S6 and Section 1 in Supplemental Materials and Methods) although Sample and Sample+ networks had the same expression values. This finding suggested additional parameters for miRNA:target interactions should be integrated if available for more accurate calculations.

Compiling realistic network and identifying nodes causing widespread perturbation

For large-scale analysis, we constructed the Real network with miRNA:target pairs from miRTarBase comprising of 3,265 genes and 581 miRNAs (Table 1). Many instances of the Real network were constructed by overlaying the expression data of each patient. We constructed 174 networks from 87 patients who had both normal and cancer expression data available. Subsequently, the perturbation efficiency of each node in all networks was calculated. The finite mixture model was used to find significantly perturbing nodes as described in “Materials and Methods”. As a result, 70 of 3,265 genes and 27 of 581 miRNAs were found to have significant perturbations in both tissues (Fig. 4A). It was observed that 29 genes and 1 miRNA had tumor tissue-specific perturbing activity. Additionally, in normal tissue samples of these 87 patients, 46 genes and 4 miRNA had showed robust perturbation efficiency. Please note that not all tumor-specific perturbing genes exhibit differential expression between normal and cancer tissues (shown in Fig. S11). Thus, our approach had potential to spot a gene that has insignificant change in expression but becomes effective in relation to expression levels of remaining one or more nodes in miRNA:target network.

In order to assess the functional annotation of highly perturbing genes, enrichment and disease associations were examined. The considerable number of these genes, exactly 55, were enriched (p-value < 0.05) in critical pathways in cancer such as PI3K-Akt signaling pathway (4.69 fold enrichment, p-value 2.04e−08), proteoglycans in cancer (6.07 fold, p-value 1.03e−07), FoxO signaling pathway (6.64 fold, p-value 4.55e−06). It was observed that 126 of these genes had enriched in biological processes and molecular functions, including negative regulation of apoptotic process (GO:0043066, 4.75 fold, and p-value 5.15e−07), cellular response to epidermal growth factor stimulus (GO:0071364, 23.12 fold, and p-value 5.33e−06), cadherin binding involved in cell-cell adhesion (GO:0098641, 8.25 fold, p-value 1.35e−11), transcription factor binding (GO:0008134, 5.32 fold, p-value 1.54e−05). The full list of enriched GO terms and KEGG pathways is presented in Table S7. Interestingly, tumor-specific perturbing genes showed enrichment in cancer-associated pathways and biological processes (Fig. 4B) while normal tissue-specific perturbing genes were not enriched significantly in the same pathways or processes. While four of the top five enriched KEGG pathways were evidently associated with cancer, the remaining pathway (focal adhesion) was shown to be associated with adhesion, migration, and invasion in breast cancer (Luo & Guan, 2010). Additionally, tumor-specific perturbing genes were enriched in several GO terms such as regulation of amino acid metabolism and apoptotic processes which play a significant role in tumor progression (Fig. 4B) (Vettore, Westbrook & Tennant, 2019). Also, the negative regulation of transcription regulation from RNA polymerase II (GO:0000122) may be important in manipulating tumor microenvironment and communications of cancer cells (Venkatraman et al., 2020). In utero embryonic development process (GO:0001701) plays role in breast cancer metastasis and mammary development in utero (Howard & Veltmaat, 2013). The full list of GO terms and KEGG pathways enriched by tumor-only perturbing genes is presented in Table S8.

DisGeNet platform (Piñero et al., 2020), a database listing critical genes for numerous diseases, was used to analyze critical genes in our findings. We found that 29 of 99 perturbing genes in tumor tissue had breast cancer disease association score greater than 0.1, with 6.95e−08 hypergeometric p-value and 2.9 fold enrichment.

In addition to highly perturbing genes, our analysis revealed highly perturbing miRNAs. In order to associate highly perturbing miRNAs with diseases, The Human microRNA Disease Database (HMDD, Version 3.2) was used (Huang et al., 2018). Almost all highly perturbing miRNAs (31 of 32, Fig. 4A) in our study were denoted as breast cancer-associated (2.42 fold enrichment with p-value 7.12e−11) in HMDD, 20 of which were attributed with causality. The remaining non-breast cancer-associated miRNA, miR-99b, may play a regulatory role in proliferation and migration processes in breast cancer by affecting the TGF-β signaling pathway (Turcatel et al., 2012).

It’s worth mentioning that the analysis thus far evaluated 87 normal and tumor samples as a whole. However, per sample basis analysis is also important since miRNAs or genes may show unanticipated behavior at the individual sample level. Figures 4C and 4D summarize perturbed node counts for each miRNA across all samples. For example, miR-30a-5p isoform which suppresses proliferation, migration, and tumor growth (Yang et al., 2017) had been observed as highly effective in almost all normal tissue samples but not in all tumor tissue samples. As another example, miR-183 which was commonly up-regulated in tumor samples and one of the key regulators of the metastatic process in breast cancer tissues (Cao et al., 2020), exhibited diverse perturbing efficiencies between normal and tumor tissues in our analysis. miR-182-5p and miR-21-5p had been found as highly efficient in tumor tissues, but not in all normal tissues. These results suggested that our ceRNA network-dependent approach can process and interpret many individual networks.

Detecting perturbation efficiency of nodes in presence of interaction parameters

For more accurate analysis we generated a network with experimental data. The interactions in the Real+ network was originated from CLEAR-CLiP and CLASH datasets as opposed to the Real network which was comprised of predicted miRNA:target interactions. A total of 1,348 genes and associated 284 miRNAs were compiled as the Real+ network (Table 1; Fig. 1A). While significant factors such as energy, seed type, and region on the target sequence in miRNA:target interaction in the template network were the same for each patient data set, expression levels of miRNA and gene were different in each data set. We calculated the perturbation efficiencies for each node in each patient.

By fitting distribution of node perturbations of simulated networks to finite mixture model, a node was defined as a “perturbing node” if it affected at least 216 nodes in at least 18 samples. It had been observed that 475 of 1,348 genes showed high perturbation efficiency in normal and tumor tissue samples. While 12 of these genes were highly effective in tumor tissues specifically, no genes were detected showing normal tissue-specific effect. On the other hand, 83 of 284 miRNAs were found to be a highly perturbing node in both tissues. KEGG or GO term enrichment analysis of 475 perturbing genes by DAVID tool (Huang, Sherman & Lempicki, 2009; Huang da, Sherman & Lempicki, 2009) (Fig. S10) showed that 60 of these genes were enriched in signal pathways like FoxO signaling pathway (p-value 0.0058 and 2.78 fold enrichment). Also, 399 of perturbing genes were significantly enriched in cancer-associated biological process and molecular functions such as cadherin binding involved in cell-cell adhesion (GO:0098641 with p-value 1.19e−08 and 3.68 fold enrichment), cell-cell adhesion (GO:0098609 with p-value 5.30e−08 and 3.65 fold enrichment), microtubule cytoskeleton organization (GO:0000226 with 5.23e−04 p-value and 4.82 fold enrichment), negative regulation of Ras protein signal transduction (GO:0046580 with p-value 3.80e−03 and 7.61 fold enrichment). The full list of enriched GO terms and KEGG pathways is presented in Table S9. When we analyzed the critical genes associated with breast cancer from the DisGeNET database (Piñero et al., 2020), 32 of the disease-related genes with a score greater than 0.1 were observed (1.50 fold enrichment with p-value 0.0051).

HMDD database analysis showed that 67 of 79 miRNAs (2.13 fold enrichment with p-value 2.05e−17) had been associated with breast cancer and 40 of these miRNAs reported as having the potential to cause breast cancer (Huang et al., 2018). The remaining 12 miRNAs not associated with breast cancer in HMDD were explored further in the literature. There are numerous studies that associate breast cancer with miR-28, miR-1287, miR-3065, miR-500a, and miR-99b (Yang et al., 2011; Schwarzenbacher et al., 2019, Martinez-Gutierrez et al., 2020; Aushev et al., 2016; Turcatel et al., 2012). They have the potential to be used for biomarker, therapeutic or regulatory purposes. Additionally, miR-501, miR-532, and miR-589-3p, which had perturbation efficiency in tumor tissue, were not associated with breast cancer in HMDD. Several researches about miR-589-5p association with drug response and abnormal regulation of miR-532-5p in breast cancer were reported (Uhr et al., 2019; Huang et al., 2020). miR-577 had been reported as a potential target for breast cancer therapy in an earlier study due to its function as a suppressor of metastasis that is induced by epithelial-mesenchymal transition (Yin et al., 2018). On the other hand, at the time of writing this study, there were no reports that directly associate miR-2116, miR-3127, and miR-501-3p with breast cancer.

Compared to simulations made with networks constructed with miRTarBase data set, the Real network, results from Real+ network showed that the high number of genes/miRNAs perturbed the large number of nodes (Chou et al., 2017). This finding emphasized the importance of integrating miRNA:target interaction parameters as they profoundly affect the proportional distribution of miRNA expression. As more miRNA:gene interaction parameters become available in future studies by novel experimental methods, our method has the potential to integrate them and thus providing more accurate depiction of ceRNA network perturbations and their consequences.

Evaluation of performance and runtime

Since human miRNA:target networks comprise thousands of nodes and interaction we evaluated the runtime of our package. To benchmark runtime, 500, 1,000, 5,000, and 10,000 miRNA:target pairs were sampled from the network constructed using miRTarBase data set. Using these networks, the runtime results of simulations originating from a single node and perturbation efficiency analysis of all nodes in networks were benchmarked. It had been determined that runtime increased significantly with increasing network size, for instance, perturbation efficiency analysis for 1,000 and 10,000 interactions takes 84 and 245 s, respectively (Figs. S13 and S14).