Introduction
============

Multiple myeloma (MM) is a malignant tumor of plasma cells occurring in bone marrow (BM), characterized by clonal proliferation of malignant plasma cells and over-production of monoclonal immunoglobulin (M-protein), and leads to severe osteolytic lesions ([@b1-ol-0-0-10100],[@b2-ol-0-0-10100]). MM has become the second most frequent hematological malignancy in the world, and each year, there are over 14,000 new cases in the US alone ([@b3-ol-0-0-10100]). Due to the obscure symptoms and unknown pathogenesis, identifying the molecular biomarkers of early diagnosis, preventive and personalized therapy, is very important and urgently needed.

MicroRNAs (miRNAs), a class of small non-coding RNAs with approximately 22 nucleotides, can play a regulatory role in biological processes, such as cell differentiation, proliferation and apoptosis ([@b4-ol-0-0-10100],[@b5-ol-0-0-10100]). Researches have demonstrated that miRNAs which can be a potential biomarker for the diagnosis, prognosis and treatment of myeloma, play a vital role in the onset, development, relapse and drug resistance of MM ([@b6-ol-0-0-10100]). It is highly demanded to identify the target genes of miRNAs, since the intricate functional mechanism of miRNAs in diseases and the regulatory role of miRNAs on the expression of target genes can be displayed by suppressing protein translation or inducing mRNAs recession ([@b7-ol-0-0-10100],[@b8-ol-0-0-10100]). Therefore, the identification of target genes contributes to further understanding the pathogenesis of miRNAs in MM.

Although a large number of the interactions between the miRNAs and mRNAs have been confirmed, the accurate prediction of miRNA target genes is still a challenge due to the limitation of prediction algorithms ([@b9-ol-0-0-10100]). For example, the prediction approach based on evolutionary conservation is only applied to the conservative functional target sites ([@b10-ol-0-0-10100]). Therefore, the comparison of differential gene expression based on genome-wide becomes a promising candidate approach to identify the specific miRNA regulation, which depends on the combination of expression profiling and sequence information, rather than solely relying on evolutionary conservation. In this study, a new prediction algorithm called targetScore was proposed, which was specifically developed to detect the target genes of a specific miRNA in exceptional cell-circumstance and operates on the whole gene set to more closely model the overall possibility instead of only on a pre-filtered set of genes. Consequently, target genes of four miRNAs (i.e. miR-148a, miR-99b, miR-20a, miR-21) associated with MM were predicted with the targetScore approach, which contributes to enhance the understanding of MM pathobiology and provide a novel way for future study in MM.

Materials and methods
=====================

### Microarray data information

The expression profile of genes associated with MM can be gained from the public functional genomics data repository named Gene Expression Omnibus (GEO; <http://www.ncbi.nlm.nih.gov/geo/>), and microarray data of GSE6477 based on the GPL96 Platform was imported. Eventually, 162 genechips about the expression data from different stages of plasma cell neoplasm can be used for further analysis, containing 73 genechips from newly diagnosed myeloma of human bone marrow plasma cells, 28 genechips from relapsed myeloma, 24 genechips from smoldering myeloma 22 genechips from monoclonal gammopathy of undetermined significance, and 15 genechips from normal donors were adopted as a control group. Furthermore, expression profile data of 12,437 genes were detected by mapping between the probe and the gene ([@b11-ol-0-0-10100],[@b12-ol-0-0-10100]).

### Differentially expressed genes (DEGs) analysis

For the gene expression profile obtained from the GSE6477, an R software package called limma was used to analyze the expression of known genes and detect the genes with larger differential expression. In order to identify the DEGs, the expression matrix of genes was subjected to t-test and F-test using the limma package, and the lmFit function was used to linearly fit the data. Furthermore, using empirical eBayes command in limma package to calculate the consensus pooled variance of genes, and adjust the associated p-value. Ultimately, 117 DEGs with logFC absolute values ≥2 and P-value \<0.05 were acquired.

### The calculation of targetScore values

TargetScore value is an integrative probabilistic score for target gene prediction of specific miRNA. Assuming there are N known genes, and *x* = *(x~1~, x~2~,......., x~N~)^T^* is regarded as the log fold-change *(x~f~)* or sequence scores *(x~l~)* value of known genes. Accordingly, *x* ∈ *(x~f~, x~1~, x~2~,......., x~L~)* is defined as the sequence scores of genes in L sets. Consequently, the targetScore value of known genes can be computed from the following formula:

t

arg

et Score

=

σ

(

\-

log

F

C

)

\[

1

K

\+

1

∑

x

ε

(

x

f

,

x

1

,

x

2

,

......

x

l

)

p

(

t

\|

x

)

\]

![](ol-17-05-4375-g00.jpg)

where

σ

(

\-

log

F

C

)

=

1

1

\+

exp

(

log

F

C

)

,

and

p

(

t

\|

x

)

![](ol-17-05-4375-g01.jpg)

is assigned as the posterior distribution of the latent variable to predict miRNA target genes for a given x. The value of posterior distribution can be calculated by integrating the prior probabilities and likelihood functions. It is worth mentioning that the targetScanCS and targetScanPCT values of the given genes as the parameter of prior probabilities can be acquired from experimental data. Additionally, the Variational Bayesian Formula and Gaussian Mixture Model (VB-GMM) can be applied to calculate the maximum likelihood functions of genes, then avoids the overlap of data ([@b13-ol-0-0-10100]). Besides, a Variational Bayesian Expectation-Maximization (VB-EM) algorithm is used to optimize the parameters of the VB-GMM. A 'target component' can be represented by the mixture component related to miRNAs in the negative fold-change and sequence score. Therefore, interaction between miRNAs and target mRNAs can be inferred by analyzing the posterior distribution of the target component of the latent variables. As a sigmoid-transformed fold-change, the targetScore value of latent variables can be computed by weighting the averaged posterior values of all diagnostic target components ([@b14-ol-0-0-10100]). This algorithm is available at [www.genelibs.com/gb](www.genelibs.com/gb).

Results
=======

### Identifications of DEGs in MM

One hundred and sixty-two samples containing MM disease and control group were obtained from the GEO database for further analysis. The limma package was used to screen the genes with differential expression. Eventually, 117 expression genes with larger difference that satisfying P-value \<0.05 and logFC absolute value \>2 were obtained for subsequent study. The distribution of DEGs between the P-values and logFC values is shown in [Fig. 1](#f1-ol-0-0-10100){ref-type="fig"}. It is worth mentioning that a gene with logFC value lower than minus 2 is a downregulated differentially expression gene, and a gene with logFC value greater than positive 2 is an upregulated DEGs. Consequently, 101 downregulated DEGs and 16 upregulated DEGs were confirmed. The relevant parameters of top 10 DEGs with the smallest P-value were selected and listed in [Table I](#tI-ol-0-0-10100){ref-type="table"}. It can be seen from [Table I](#tI-ol-0-0-10100){ref-type="table"} that the logFC values of the DEGs are negative values and the absolute values are \>2. Moreover, P-values and adjusted P-values of the DEGs are far lower than 0.05, indicating that the DEGs obtained by the GEO and limma package have a very important statistical significance.

### The acquisition of targetScore values

Acquiring the logFC, targetScanCS and targetScanPCT values of the miRNAs, miR-148a, miR-99b, miR-20a and miR-21 from the real data, and then import these results into equation ([@b1-ol-0-0-10100]) to compute the targetScore values. Consequently, a total of 312 latent target genes of miR-148a were selected with a targetScore values \>0.5, 247 latent target genes of miR-99b were selected with a targetScore values \>0.3, 154 latent target genes of miR-20a were selected with a targetScore values \>0.3, and 121 latent target genes of miR-21 were selected with a targetScore values \>0.4. Theoretically, the greater the values of targetScore, the greater the probability of the miRNA targets. Therefore, the optimal target genes can be confirmed by further screening.

### Screening of the optimal target genes

In order to identify the optimal target gene, the intersection of latent target genes from the four miRNAs was plotted and shown in [Fig. 2A](#f2-ol-0-0-10100){ref-type="fig"}. It can be seen from [Fig. 2A](#f2-ol-0-0-10100){ref-type="fig"} that a total of 147 potential target genes were identified from the genes of the four miRNAs by integrating the bioinformatical analysis. Noting that 1 and 23 consistently expressed genes with intersection above four and three miRNAs were confirmed, respectively. Then taking the intersection of the genes screened by targetScore value for the four miRNAs and genes selected by differential expression. Consequently, 34 optimal target genes associated with MM were identified (as showed in [Fig. 2B](#f2-ol-0-0-10100){ref-type="fig"}). Besides, inspired by the above results, we used the Cytoscape software to construct the interaction network between miRNAs and target mRNAs. Consequently, the interaction network between mRNAs and miRNAs, including hsa-miR-99b, hsa-miR-20a, hsa-miR-148a and hsa-miR-21, was constructed (data not shown). Subsequently, we integrated mRNAs associated with these four miRNAs. As a result, the co-expressed genes between these miRNAs were displayed in [Fig. 3](#f3-ol-0-0-10100){ref-type="fig"}, which showed that a miRNA regulated by multiple miRNAs may be closely associated with certain diseases. Furthermore, the optimal target genes (blue node) were screened out. It is worth mentioning that latent target genes among these interactions were validated (red node). There are consistently expressed genes between the optimal target genes and the validated genes, such as TGFBI, HIF1A, TGFBR2. Furthermore, the relevant biological parameter information, including P-value, targetScore, targetScanCS, targetScanPCT, and logFC, was filtered and listed in [Table II](#tII-ol-0-0-10100){ref-type="table"}.

### Gene mutations in the MM

Identifying the mutated genes in MM is great promise for personalized treatment, whereby the patients who are associated with MM and occurred specific gene mutations would be considered using the appropriate targeted for treatment ([@b15-ol-0-0-10100],[@b16-ol-0-0-10100]). In order to confirm the meaningfully mutated genes, Lohr *et al* ([@b17-ol-0-0-10100]) have selected 203 patients with multiple myeloma to identify parallel sequencing of paired tumor/normal samples. Consequently, 879 mutated genes associated with MM were verified ([@b17-ol-0-0-10100]). Then integrating the genes obtained by targetScore and the mutated genes, total of 41 consistently expressed genes were extracted (as shown in [Table III](#tIII-ol-0-0-10100){ref-type="table"}), including 1 target gene with 36 mutations, 2 target genes with 5 mutations, 5 target genes with 4 mutations, 5 target genes with 3 mutations, 10 target genes with 2 mutations and 18 target genes with 1 mutation. Furthermore, 5 optimal target genes (SYK, LCP1, HIF1A, ALDH1A1, MAFB) were confirmed by intersecting the genes in [Tables II](#tII-ol-0-0-10100){ref-type="table"} and [III](#tIII-ol-0-0-10100){ref-type="table"}, the venn plot of gene intersection was shown in [Fig. 4](#f4-ol-0-0-10100){ref-type="fig"}.

Discussion
==========

In order to seek for the causes and possible mechanisms of multiple myeloma formation and development, a large number of basic and clinical studies have been performed in the past few decades. However, the pathogenesis of MM is still unclear, and individual differences in treatment outcomes of MM are greater. In addition, the prognosis evaluation index cannot accurately reflect the complexity of the disease, causing a still high mortality of MM in the world ([@b18-ol-0-0-10100]). Therefore, it may be the key to effectively solve the clinical treatment problem by deeply studying the pathogenesis and prognostic factors of MM disease. Numerous studies have suggested that miRNAs have the potential for diagnosis and prognosis in cancers, and on the basis of cell properties and developmental stage of tumors, the expression profile of miRNAs can be used to distinguish cancer with higher accuracy ([@b19-ol-0-0-10100],[@b20-ol-0-0-10100]). As a molecular biomarker of disease diagnosis, the abnormal expression of the miRNAs have been demonstrated to cause the drug resistance of MM cells ([@b21-ol-0-0-10100]). Therefore, mRNAs as the downstream part of miRNAs, the target genes associated with MM are considered for studying the mechanism of miRNAs in MM.

miR-21 as a member of oncogenic miRNAs, is aberrantly expressed in hematological malignancies, like myeloma ([@b22-ol-0-0-10100]). Similarly, the altered expression levels of other miRNAs (miR-148a, miR-99b and miR-20a) in MM patients have been investigated by Huang *et al* ([@b23-ol-0-0-10100]). Demonstrating that miR-21, miR-148a, miR-99b and miR-20a are potential prognostic biomarkers of MM disease, and their aberrant expression provides the possible potential for identification, classification and prognosis determination of the tumor. For example, the expression levels of miR-20a and miR-148a were associated with the prognosis of MM and the expression of miR-99b was affected by the karyotype of MM disease ([@b23-ol-0-0-10100]).

In this study, we used microarray technology to improve the demand for the genome-wide, and gene expression levels of 162 samples, including 73 new MM patients, 24 smoldering MM patients, relapsed MM patients and 15 normal donor controls, were adopted. Furthermore, a new target prediction algorithm called targetScore with high accuracy, that does not rely on the genetic conservatism was used. Eventually, approximately 22% (i.e. 147) predicted target genes regulated by at least 2 miRNAs were identified. Particularly, gene HLA-DPB1 is generally regulated by 4 miRNAs, which provides instructions for the preparation of proteins, and plays an important role in the immune system. Additionally, integrating the analysis of miRNAs-mRNAs expression data to predict optimal target genes, namely HLA-DPB1, P2RY13, DPYSL2, TGFBI, MS4A6A, KCTD12, TMEM156 and PDK4 with larger differential expression and higher targetScore values (more than 0.5) are selected. As the optimal target genes, the interaction of the 8 genes with relevant miRNAs was shown in [Fig. 3](#f3-ol-0-0-10100){ref-type="fig"}. It is worth mentioning that the inhibitory effect of TGFBI gene as the target gene for oncogenic miR-21 has been verified by Lu *et al* ([@b24-ol-0-0-10100]).

For MM patients, the mutations of genes are commonly found in subclonal populations. Therefore, the mutated genes could become the partial strategy of the targeted treatment, and the importance of clinical treatment decisions has been reported ([@b25-ol-0-0-10100]). In this study, through the integration of the mutated genes in MM patients and genes predicted by targetScore algorithm, 41 possible target genes were selected ([Table III](#tIII-ol-0-0-10100){ref-type="table"}). Moreover, 13 optimal target genes with more than mutations were selected from the 41 possible target genes, including NRAS, ZFHX3, NIDI, BCL7A, DICER1, RBI, SHYN1, CAMTA1, SYNE1, CDKN1B, CYLD, MAF and MAFB. Specifically, mutations occur frequently in NRAS gene, and the number of mutations even reached 36, which could cause the cell malignant proliferation and metastasis. In addition, studies have shown that the low expression of CYLD with a mutation number of 3 is related to disease development from monoclonal gammopathy of undetermined significance to MM, and the functional detection results indicated that the expression of CYLD suppresses the cell proliferation and survival of MM disease ([@b26-ol-0-0-10100]). It is worth mentioning that the function mechanism of RBI gene on the tumor has been studied, and the results revealed that RBI can be considered as the anti-proliferative target of miR-20a ([@b27-ol-0-0-10100]). Besides, in these 41 possible target genes, the targeted functions of HIF1A and TGFBR2 on the miR-20a and miR-21 have been studied ([@b28-ol-0-0-10100],[@b29-ol-0-0-10100]). Therefore, the above analysis shows that the target genes predicted by the targetScore algorithm have a higher accuracy.

Our present study revealed the potential target genes of miRNAs associated with MM. The biological information of relevant target genes was obtained by the GEO database and a novel target prediction approach called targetScore. We integrated the sequence information and expression profile of genes and then we further computed the targetScore value of the obtained genes using VB-GMM method. Eventually, a total of 121,18 and 154 possible target genes of the miRNAs, miR-21, miR-99b and miR-20a with targetScore value \>0.4, respectively, were detected. Particularly, the number of predicted target genes for miRNAs with targetScore \>0.5 reached 312. And 147 target genes were targeted by at least 2 miRNAs, where 34 of these 147 target genes with larger differential expression were selected. Furthermore, 41 target genes were screened out by the intersection of the mutated genes and genes obtained by targetScore algorithm, and three genes, RBI, HIF1A and TGFBR2 as the targets of miR-20a, miR-21 have been validated. However, the functional effect of target genes obtained in the present study on the pathogenesis of miRNAs in MM is not clear, so the pathway and GO enrichment analysis of targets in MM is necessary.

We thank Dr Xin Yu for the help and advice for this project.

Funding
=======

This study was supported by Heilongjiang Harbin Science and Technology Bureau City Youth Reserve Talent Project (2014RFQGJ091), Science and Technology Commission of Heilongjiang Province (2016--235), and Application Technology Research and Development Project of Harbin City, Heilongjiang Province (2013AA3BS020).

Availability of data and materials
==================================

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors\' contributions
=======================

YX performed the experiments and analyzed the data, and was also a major contributor in writing the manuscript. HL and GN made a substantial contribution to the study conception and design. YX, HL and GN performed the statistical analysis. JZ was involved in the conception and design of the study as well as in the drafting of the manuscript. All authors read and approved the manuscript.

Ethics approval and consent to participate
==========================================

Not applicable.

Patient consent for publication
===============================

Not applicable.

Competing interests
===================

The authors declare that they have no competing interests.

![Volcano plot of the aberrantly expressed genes screened from 162 MM samples.](ol-17-05-4375-g02){#f1-ol-0-0-10100}

![(Α) Identification of 147 possible target genes with at least 2 miRNAs intersection, and (B) 34 consistently expressed genes by the interaction of DEGs and targetScore value genes.](ol-17-05-4375-g03){#f2-ol-0-0-10100}

![Interaction network of co-expressed genes between mRNs and miRNAs. miRNAs are represented in yellow, mRNAs are represented in other colours (green, red and blue). Notably, the red nodes represent the genes that have been validated, and blue nodes are the optimal target genes.](ol-17-05-4375-g04){#f3-ol-0-0-10100}

![Optimal target genes obtained by the intersection of the 34 DEGs and the 41 mutated genes.](ol-17-05-4375-g05){#f4-ol-0-0-10100}

###### 

The relevant parameters of top 10 DEGs with the smallest P-value.

  DEGs     logFC       AveExpr    t            P-value    adj.P.Val   B
  -------- ----------- ---------- ------------ ---------- ----------- -----------
  RNASE2   −3.919468   8.845879   −17.63099    1.07E-39   1.33E-35    77.958873
  CLC      −4.696261   7.664386   −13.052750   3.72E-27   2.31E-23    50.463704
  PRG3     −3.546961   8.184811   −12.19434    9.35E-25   3.88E-21    45.159977
  RNASE3   −3.870937   8.008265   −10.85179    5.10E-21   1.59E-17    36.885963
  PRG2     −5.654262   8.129576   −10.70904    1.27E-20   3.15E-17    36.011896
  ELANE    −4.391511   5.948464   −9.631113    1.12E-17   1.98E-14    29.480459
  DEFA4    −2.035331   8.836430   −9.494378    2.61E-17   4.05E-14    28.663198
  MPO      −2.329628   7.083587   −9.196154    1.64E-16   2.26E-13    26.891907
  LYVE1    −3.327286   6.777089   −8.601285    6.07E-15   7.55E-12    23.411790
  EPX      −2.80708    8.27019    −8.414076    1.86E-14   2.10E-11    22.333413

###### 

Analysis results of 34 optimal target genes selected by intersecting the targetScore and differential expression values.

  Target genes   logFC       targetScanCS   targetScanPCT   targetScore   P-value
  -------------- ----------- -------------- --------------- ------------- ----------
  SYK            −2.396184   −0.03            0             0.547351      6.72E-06
  KCTD12         −2.762166     0            −0.22           0.650777      4.82E-11
  C1QA           −2.595153     0            −0.22           0.643825      7.46E-08
  P2RY13         −2.924695   −0.171           0             0.660584      1.17E-12
  PRG2           −5.654262     0              0             0.402927      1.27E-20
  TMEM156        −2.149213     0            −0.22           0.619641      7.75E-06
  LUM            −2.235238   −0.341           0             0.594835      0.000323
  HLA-DPB1       −2.169956     0            −0.22           0.620972      2.44E-07
  ANXA1          −2.117414   −0.248           0             0.564145      1.03E-05
  HMOX1          −2.602141   −0.12            0             0.635262      5.02E-10
  CXCL12         −2.248040   −0.079         −0.22           0.904481      3.30E-07
  HLA-DPA1       −2.789554     0            −0.22           0.651823      5.23E-10
  DPYSL2         −2.023163   −0.189         −0.9            0.800525      0.000448
  PID1           −2.417350   −0.154         −0.06           0.904022      1.93E-07
  ALDH1A1        −2.596302   −0.337           0             0.643168      1.34E-08
  NRP1           −2.108773   −0.425         −0.8            0.891752      2.72E-12
  TGFBI          −3.190161   −0.346         −0.6            0.960319      5.87E-10
  RNASE6         −3.864246     0            −0.04           0.657398      1.42E-06
  EPX            −2.80708    −0.107           0             0.627697      1.86E-14
  PTPN22         −2.278062     0            −0.06           0.613351      1.55E-08
  CD27           −2.420251     0            −0.05           0.634639      2.60E-07
  LCP1           −2.049416     0            −0.07           0.612925      2.50E-06
  GPR137B        −2.67012    −0.126         −0.51           0.931729      6.99E-08
  PDK4           −2.134344   −0.197         −0.57           0.894196      1.04E-06
  HIF1A          −2.091514   −0.179         −0.62           0.825931      1.12E-08
  CLC            −4.696261     0              0             0.400681      3.72E-27
  CD5L           −2.070227   −0.054           0             0.522844      5.17E-06
  HIF1A          −2.091514   −0.179         −0.62           0.825932      1.12E-08
  RGS1           −2.269436   −0.112         −0.09           0.876758      0.000206
  IGLV1-44       −4.651009     0              0             0.400512      2.34E-10
  AHR            −2.440098   −0.199         −0.03           0.779391      6.36E-06
  MS4A6A         −2.461868   −0.147           0             0.630938      4.39E-11
  MAFB           −3.830404   −0.376         −0.89           0.978761      6.50E-09
  KLF4           −2.130986   −0.267         −0.38           0.893878      1.85E-06

###### 

Characteristics of 41 consistently expressed genes obtained by integrating the targetScore values and mutations of genes.

  miRNA      Genes      logFC       Target-ScanCS   Target-ScanPCT   Target-Score   Mutation
  ---------- ---------- ----------- --------------- ---------------- -------------- ----------
  miR-148a   NRAS       −0.520168   −0.335          −0.88            0.555780       36
  miR-148a   ZFHX3      −0.757232   −0.118          −0.28            0.655262       5
  miR-148a   NID1       −0.729052   −0.157          −0.22            0.644748       5
  miR-21     BCL7A      −0.814101   −0.216          −0.34            0.465047       4
  miR-148a   DICER1     −0.500383   −0.15           −0.76            0.547017       4
  miR-20a    RB1        −1.164287   −0.036          −0.57            0.509868       4
  miR-99b    SHTN1      −1.687627     0               0              0.312160       3
  miR-99b    CAMTA1     −0.241022   −0.227          −0.11            0.445970       3
  miR-20a    SYNE1      −0.671450   −0.189          −0.72            0.442826       3
  miR-148a   CDKN1B     −0.963022   −0.201          −0.83            0.717540       3
  miR-20a    CYLD       −0.448714   −0.045          −0.37            0.406729       3
  miR-21     ALDH1A1    −2.596302   −0.337            0              0.643168       2
  miR-99b    MMP8       −1.486319     0               0              0.301628       2
  miR-148a   MMP13      −1.390973   −0.365          −0.05            0.799989       2
  miR-20a    TNFRSF21   −0.767567   −0.283          −0.69            0.457822       2
  miR-21     ZFP36L2    −0.738853   −0.122          −0.16            0.453354       2
  miR-21     CD36       −1.595103   −0.066          −0.06            0.623486       2
  miR-99b    HIF1A      −2.091514     0               0              0.329239       2
  miR-20a    PDGFRA     −0.556430   −0.192          −0.9             0.424678       2
  miR-99b    MAF        −1.499598     0               0              0.302369       4
  miR-21     MEF2A      −0.538136   −0.078          −0.06            0.421767       1
  miR-148a   MAP3K4     −0.472509   −0.366          −0.51            0.534801       1
  miR-99b    SGK1       −1.491325     0               0              0.301908       1
  miR-148a   IRS2       −1.239373     0             −0.24            0.536011       1
  miR-21     SESN1      −0.851813   −0.097          −0.2             0.470898       1
  miR-21     STAT3      −0.756632   −0.222          −0.29            0.456121       1
  miR-148a   ARL6IP1    −0.604213   −0.403          −0.85            0.593083       2
  miR-99b    SYK        −2.396184     0               0              0.339027       1
  miR-20a    ANKRD12    −0.751340   −0.089          −0.9             0.455298       1
  miR-21     TGFBR2     −0.634521   −0.115          −0.47            0.437027       1
  miR-148a   PBX1       −0.581344   −0.073          −0.05            0.582475       1
  miR-148a   CDH1       −1.284355     0             −0.05            0.540909       1
  miR-21     KLF6       −1.196904   −0.105          −0.35            0.526831       1
  miR-99b    MAFB       −3.830404     0               0              0.362044       4
  miR-148a   MED12      −0.485766   −0.163          −0.11            0.540587       1
  miR-20a    DUSP2      −1.228701   −0.364          −0.96            0.532599       2
  miR-99b    LCP1       −2.049416     0               0              0.327691       1
  miR-20a    SHOC2      −0.403998   −0.089          −0.56            0.400203       1
  miR-148a   TBL1XR1    −0.447336   −0.121          −0.58            0.523952       1
  miR-99b    LMO2       −1.828608     0               0              0.318704       1
  miR-148a   EPAS1      −0.983346   −0.152          −0.65            0.722513       1

[^1]: Contributed equally
