SPMLMI: predicting lncRNA–miRNA interactions in humans using a structural perturbation method

Data collection

Our goal was to develop a method that is able to predict candidate interactions between lncRNAs and miRNAs. Three kinds of data were collected from various databases and a bilayer network of lncRNAs and miRNAs was constructed to represent the complex relationships. First, known lncRNA–miRNA interaction data were downloaded from the lncRNASNP database version 2.0 (http://bioinfo.life.hust.edu.cn/lncRNASNP) (Miao et al., 2018) and used to construct an lncRNA–miRNA network. Second, expression data for lncRNAs were downloaded from the NONCODE database (http://www.noncode.org/) (Fang et al., 2018) and used to calculate the expression similarity between lncRNAs and to construct the lncRNA similarity network. Additionally, miRNA expression data were downloaded from the miRmine database (https://guanfiles.dcmb.med.umich.edu/mirmine/index.html) (Panwar, Omenn & Guan, 2017) and used to construct the miRNA similarity network. Table 1 shows the data obtained from publicly available databases.

Construction of an lncRNA–miRNA bilayer network

An lncRNA–miRNA bilayer network containing two kinds of nodes (miRNAs and lncRNAs) and three types of relationships (miRNA–miRNA similarity relationship, lncRNA–lncRNA similarity relationship, as well as the known lncRNA–miRNA interaction relationship) was constructed.

To construct the lncRNA–miRNA interaction network, data were downloaded from the version 2.0 of the lncRNASNP database (Miao et al., 2018). This database contains known lncRNA–miRNA interactions with experimental validation collected from version 2.0 of the starBase (Li et al., 2014). Then, we define matrix LMnet as the known lncRNA–miRNA interaction network, where if lncRNA i and miRNA j are connected, the element LMnet = 1; otherwise, LMnet = 0. A total of 39,366 associations between 3,150 lncRNAs and 262 miRNAs were retained after removing duplicates. Figure 1B presents a simple example of lncRNA–miRNA interaction network construction.

The lncRNA similarity network and miRNA similarity network were constructed based on expression profiles. First, expression profile data for miRNAs were downloaded from the miRmine database (Panwar, Omenn & Guan, 2017). The database contains data for 2588 miRNAs in 135 different human tissues and cell lines. After removing missing items and matching the miRNA IDs in the lncRNA–miRNA interaction network (LMnet),expression data for 262 miRNAs in 124 different human tissues and cell lines were retained. Second, lncRNA expression profiles were downloaded from the NONCODE database (Fang et al., 2018). Expression profiles of 171,814 transcripts in 16 different human tissues and eight cell lines were obtained. After the removal of duplicates and matching lncRNA IDs in LMnet, 3,150 lncRNA expression records for 24 different human tissues and cell lines were retained.

To describe the expression characteristics of ncRNAs in different tissues or cell lines in an unbiased fashion, the expression profiles were standardized before calculating pairwise similarity values. The Pearson correlation coefficient was used to measure the similarity of ncRNAs and the expression similarity matrix of lncRNAs and miRNAs was obtained. The similarity matrix represented the similarity network of these two types of ncRNAs, referred to as LSnet and MSnet, respectively, as shown in Figs. 1A and 1C. For a pair of RNAs, a higher correlation indicates greater similarity in expression, in general.

Finally, the above three networks are integrated to build a bilayer network. Specifically, the lncRNA–miRNA bilayer network can be represented by a block adjacency matrix A ∈ R^N×N, (1)${\displaystyle A=\left[\begin{array}{cc}LSent\hfill & LMnet\hfill \\ LMne{t}^T\hfill & MSnet\hfill \end{array}\right],}$

where N = 3,412 (the total number of nodes in the network, 3,150 lncRNAs and 262 miRNAs) and LMnet^T is the transpose of LMnet.

Overall, the lncRNA–miRNA bilayer network was composed of the miRNA similarity network, lncRNA similarity network, and edges connecting the two networks. Figure 1 demonstrates the workflow for constructing the bilayer network and its matrix representation.

Structural consistency and structural perturbation method

Structural consistency was first proposed by Lü et al. (2015) and can be used to evaluate the likelihood of links prediction in complex network. It is defined as the consistency of network structural features before and after the elimination of partial associations at random. For the rationale and the specific derivation, refer to the supplementary material (File S1). In this study, this method was applied to the lncRNA–miRNA bilayer network A and its sub-networks (LSnet and MSnet) and the structural consistency of each network was evaluated separately.

Generally, the link prediction problem refers to the problem of estimating the probability of the existence of unobserved links according to known topological information. The network structure perturbation method involved in the structure consistency calculation process can be used to predict missing links (Lü et al., 2015). After obtaining the lncRNA-miRNA bilayer network (adjacency matrix, A) in ‘Construction of an lncRNA-miRNA bilayer network’, we derive the perturbation matrix A′ by the structural perturbation method. The prediction matrix $\widehat{A^{\prime }}$ is then obtained by averaging t independent perturbation calculations. See the supplementary materials for more details (File S1). In this way, elements in the prediction matrix $\widehat{A^{\prime }}$ can be interpreted as scores for each pair of nodes in the lncRNA–miRNA bilayer network. The scores in $\widehat{A^{\prime }}$ determine the extent of all unobserved lncRNA–miRNA interactions, and we assume that a higher score corresponds to a greater likelihood of a potential interaction.

Evaluation metrics

To evaluate the performance of SPMLMI, k-fold cross-validation was performed based on verified lncRNA–miRNA interactions downloaded from the lncRNASNP database (Miao et al., 2018). In the cross-validation process, the original set was randomly divided into k equally sized subsets. Among these k subsets, one was used as the validation data for testing the model, and the others were used as the training data. Specifically, we randomly reset 1∕k known interaction entries of the training dataset to unknown, and prioritize the candidate interaction entries in the prediction matrix by SPMLMI scores. The known interactions are considered as positive, and the rest as negative. True positive rate is the proportion of overlap between the top 1∕k ranked candidate interactions with the known interactions. Then, the cross-validation was repeated k times, taking each of the k subsets as the validation data in turn. The cross-validation results are averaged to produce a single estimation.

The AUC (area under the receiver operating characteristic curve) was estimated to evaluate the performance of lncRNA–miRNA interaction prediction results. In particular, the TPR (true positive rate) and FPR (false positive rate) were calculated by varying the threshold and the ROC (receiver operating characteristic) curve and AUC value were obtained, where TPR and FPR represent the percentage of test samples above or below the given thresholds, respectively.

In addition, we applied SPMLMI to datasets analyzed by two additional methods, EPLMI (Huang, Chan & You, 2018) and SNFHGILMI (Fan, Cui & Zhu, 2020), for an objective comparison of performance.

Structural consistency of different networks

We calculated the structural consistency of three datasets (i.e., the EPLMI dataset provided by Huang, Chan & You (2018), SNFHGILMI dataset provided by Fan, Cui & Zhu (2020), and SPMLMI dataset obtained in this study). Each group of datasets was composed of an miRNA similarity network (MSnet), an lncRNA similarity network (LSnet), and a bilayer network (Bilayer − net). We randomly selected 10% of the total links from each network as the perturbation set (note that the bilayer network contains two kinds of links). As shown in Fig. 2A, in three different datasets, the structural consistency values for the lncRNA–miRNA bilayer network (Bilayer − net, yellow bar) were 0.2413 ± 0.0015, 0.3036 ± 0.0011, and 0.2168 ± 0.0005, respectively. The structural consistency of these three bilayer networks (Bilayer − net) were all greater than that of the similarity network (where the blue bar is MSnet and the orange bar is LSnet in Fig. 2), indicating that considering more information can improve the structural consistency of the network. These results show that the construction of the bilayer network improves the inherent link predictability of the network and will be helpful for the prediction of “missing links” (i.e., newly established or novel interactions).

Performance evaluation and parameter optimization

First, to evaluate the performance of our proposed SPMLMI for the prediction of lncRNA–miRNA interactions, we adopted 2-fold, 5-fold and 10-fold cross-validation frameworks. To avoid the bias of sample division in cross-validation, we repeated each experiment 10 times and took the average value. Figure 2C shows the corresponding ROC curve for each k-fold cross-validation. The average areas under the curves were 0.8617 ± 0.0013, 0.9512  ± 0.0034, and 0.9517 ± 0.0016 for 2-fold, 5-fold, and 10-fold cross-validation, respectively. These results indicate that k-fold cross-validation increase improves the performance of SPMLMI. There was no significant difference (paired t-test, p-value = 0.1997) between 5-fold and 10-fold cross-validation. Therefore, in subsequent experiments, we use 5-fold cross-validation to save time.

Second, to optimize t (i.e., the number of perturbations), we tuned the parameter from 1 to 16. Figure 2B shows boxplots of the variation of AUC values with increasing t. The AUC values increased as t increased; however, for values of t greater than or equal to 8, AUC was relatively stable. For simplicity, we set the parameter t to 8 for subsequent experiments.

Comparison with previous methods

Furthermore, we compared SPMLMI with the previously established EPLMI (Huang, Chan & You, 2018) based on the two-way diffusion algorithm and SNFHGILMI (Fan, Cui & Zhu, 2020) based on the heterogeneous graph inference method. Initially, we evaluated the same dataset (i.e., the SPMLMI dataset obtained in this study, see ‘Materials & Methods’ for details) by these three methods. Figure 3A shows the ROC curves for the prediction performance of these three models. The average AUC for SPMLMI was 0.9512 ± 0.0034, which was higher than the AUC values of 0.7999 ± 0.0014 for EPLMI and 0.5163 ± 0.0013 for SNFHGILMI, respectively. We also compared the methods using the original datasets. Using the EPLMI dataset provided by Huang, Chan & You (2018), the average AUC values for the three methods (SPMLMI, EPLMI, and SNFHGILMI) were 0.8767 ± 0.0033, 0.8417 ± 0.0017, and 0.5184 ± 0.0023, respectively (Fig. 3B). For the SNFHGILMI dataset provided by Fan, Cui & Zhu (2020), the average AUC values were 0.8653 ± 0.0021, 0.8570 ± 0.0012, and 0.9426 ± 0.0035, respectively (Fig. 3C).

These results prove that the newly established SPMLMI has good performance and generalizability across datasets.

Case studies

Two lncRNAs, SNHG14 (UBE3A-ATS/NONHSAT041137) and MALAT1 (NONHSAT022132), with important biological functions were selected to validate the prediction performance of SPMLMI. SNHG14, also known as UBE3A-ATS, overlaps with the entire UBE3A gene and its paternal expression is thought to negatively regulate paternal UBE3A expression in neurons, thereby determining tissue-specific imprinting of UBE3A in the brain (Galiveti et al., 2014; Stanurova et al., 2016). MALA1 is located less than 70 kb from NCRNA00084 on chromosome 11q13.1 (Hutchinson et al., 2007) and is dysregulated in many diseases, such as alcohol-related diseases, endometrioid endometrial carcinoma, and non-small cell lung cancer (Ji et al., 2003; Kryger et al., 2012; Li et al., 2016). For each of this two lncRNAs in the interaction network, known interaction entries were randomly reset to unknown. We implemented SPMLMI using the reset data set and prioritized all candidate miRNAs according to their scores in the prediction matrix.

As shown in Table 2, among the top 10% (26) predicted SNHG14-related miRNAs, 18 lncRNA–miRNA interactions were confirmed by experimental validation (containing 17 of the 21 known interactions we collected from the lncRNASNP2 database Miao et al., 2018). It should be noted that the miRNA hsa-miR-126-3p (marked with an asterisk in Table 2), ranked 16th, was not collected from the lncRNASNP2 database (Miao et al., 2018) but was reported by an experimental study of transcriptome-wide Ago2:RNA interactions in human brain tissues using HITS-CLIP (i.e., CLIP-seq) (Boudreau et al., 2014). This information was retrieved from DIANA-LncBase v3 (Karagkouni et al., 2020). Several studies have shown that miR-126 directly targets the 3′-UTR of insulin receptor substrate-1 (IRS-1) (Ryu et al., 2011; Zhang et al., 2008; Zhou et al., 2013). Over-expression of miR-126 leads to downregulation of IRS-1 expression, and a consequent impairment in insulin signaling (Ryu et al., 2011). Moreover, the research results of Zhou et al. (2013) showed over-expression of miR-126 down-regulated IRS-1, suppressed AKT and ERK1/2 activation, CRC cells proliferation, migration, invasion, and caused cell cycle arrest. Therefore, lncRNA SNHG14, as an interacting molecule of miR-126, may play a potential role in regulating cell metabolism and cell cycle.

Similarly, for MALAT1, we examined the candidate miRNAs ranked by prediction scores. The results showed that the top 10% of predicted MALAT1-related miRNAs were all verified in published studies (Table S1). The reason is most likely that lncRNA MALAT1 has more known interaction information (containing 113 known interactions collected from the lncRNASNP2 database (Miao et al., 2018)). Therefore, we examined more predicted candidate miRNAs and found that a total of 106 were verifiable. Additionally, in a genome-wide analysis of miRNA–mRNA interactions in marrow stromal cells, the interaction between lncRNA MALAT1 and miRNA hsa-miR-487b-3p was confirmed (Balakrishnan et al., 2014); this candidate miRNA was ranked 91st (Table S1) based on prediction scores. The experimental study by Xi et al. (2013) demonstrated that miR-487b directly targets gene Wnt5a. The protein Wnt5a is an important activator of the Wnt/Ca²⁺ signaling pathway, which plays an important role in embryonic development, cell differentiation, inflammation and other physiological processes (De, 2011). Because of sharing miRNA miR-487b with Wnt5a gene, lncRNA MALAT1 may also be involved in this complex signaling network.

The above results suggest that SPMLMI has good prediction ability and can be used to screen candidate target miRNAs for a given lncRNA.