Identifying the optimal target genes associated with multiple myeloma by a novel bioinformatical analysis
- Authors:
- Published online on: March 4, 2019 https://doi.org/10.3892/ol.2019.10100
- Pages: 4375-4382
-
Copyright: © Xue et al. This is an open access article distributed under the terms of Creative Commons Attribution License.
Abstract
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 (1,2). 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 (3). 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 (4,5). 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 (6). 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 (7,8). 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 (9). For example, the prediction approach based on evolutionary conservation is only applied to the conservative functional target sites (10). 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 (11,12).
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 = (x1, x2,……., xN)T is regarded as the log fold-change (xf) or sequence scores (xl) value of known genes. Accordingly, x ∈ (xf, x1, x2,……., xL) 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:
target Score=σ(-logFC)[1K+1∑xε(xf,x1,x2,……xl)p(t|x)] where σ(-logFC)=11+exp(logFC), and p(t|x)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 (13). 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 (14). This algorithm is available at 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. 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. It can be seen from Table I 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 (1) 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. It can be seen from Fig. 2A 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). 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, 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.
Table II.Analysis results of 34 optimal target genes selected by intersecting the targetScore and differential expression values. |
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 (15,16). In order to confirm the meaningfully mutated genes, Lohr et al (17) 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 (17). Then integrating the genes obtained by targetScore and the mutated genes, total of 41 consistently expressed genes were extracted (as shown in Table III), 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 and III, the venn plot of gene intersection was shown in Fig. 4.
Table III.Characteristics of 41 consistently expressed genes obtained by integrating the targetScore values and mutations of genes. |
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 (18). 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 (19,20). As a molecular biomarker of disease diagnosis, the abnormal expression of the miRNAs have been demonstrated to cause the drug resistance of MM cells (21). 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 (22). 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 (23). 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 (23).
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. 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 (24).
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 (25). 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). 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 (26). 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 (27). Besides, in these 41 possible target genes, the targeted functions of HIF1A and TGFBR2 on the miR-20a and miR-21 have been studied (28,29). 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.
Acknowledgements
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.
References
Eslick R and Talaulikar D: Multiple myeloma: From diagnosis to treatment. Aust Fam Physician. 42:684–688. 2013.PubMed/NCBI | |
Falank C, Fairfield H and Reagan MR: Signaling interplay between bone marrow adipose tissue and multiple myeloma cells. Front Endocrinol (Lausanne). 7:672016. View Article : Google Scholar : PubMed/NCBI | |
de Mel S, Lim SH, Tung ML and Chng WJ: Implications of heterogeneity in multiple myeloma. BioMed Res Int. 2014:2325462014. View Article : Google Scholar : PubMed/NCBI | |
Xie S, Zhang Y, Qu L and Xu H: A Helm model for microRNA regulation in cell fate decision and conversion. Sci China Life Sci. 56:897–906. 2013. View Article : Google Scholar : PubMed/NCBI | |
Peng L, Li Y, Zhang L and Yu W: Moving RNA moves RNA forward. Sci China Life Sci. 56:914–920. 2013. View Article : Google Scholar : PubMed/NCBI | |
Zhang J, Xiao X and Liu J: The role of circulating miRNAs in multiple myeloma. Sci China Life Sci. 58:1262–1269. 2015. View Article : Google Scholar : PubMed/NCBI | |
Chen J, Deng S, Zhang S, Chen Z, Wu S, Cai X, Yang X, Guo B and Peng Q: The role of miRNAs in the differentiation of adipose-derived stem cells. Curr Stem Cell Res Ther. 9:268–279. 2014. View Article : Google Scholar : PubMed/NCBI | |
Yang Y, Li F, Saha MN, Abdi J, Qiu L and Chang H: miR-137 and miR-197 induce apoptosis and suppress tumorigenicity by targeting MCL-1 in multiple myeloma. Clin Cancer Res. 21:2399–2411. 2015. View Article : Google Scholar : PubMed/NCBI | |
Alexiou P, Maragkakis M, Papadopoulos GL, Reczko M and Hatzigeorgiou AG: Lost in translation: An assessment and perspective for computational microRNA target identification. Bioinformatics. 25:3049–3055. 2009. View Article : Google Scholar : PubMed/NCBI | |
Friedman RC, Farh KK, Burge CB and Bartel DP: Most mammalian mRNAs are conserved targets of microRNAs. Genome Res. 19:92–105. 2009. View Article : Google Scholar : PubMed/NCBI | |
Chng WJ, Kumar S, Vanwier S, Ahmann G, Price-Troska T, Henderson K, Chung TH, Kim S, Mulligan G, Bryant B, et al: Molecular dissection of hyperdiploid multiple myeloma by gene expression profiling. Cancer Res. 67:2982–2989. 2007. View Article : Google Scholar : PubMed/NCBI | |
Tiedemann RE, Zhu YX, Schmidt J, Yin H, Shi CX, Que Q, Basu G, Azorsa D, Perkins LM, Braggio E, et al: Kinome-wide RNAi studies in human multiple myeloma identify vulnerable kinase targets, including a lymphoid-restricted kinase, GRK6. Blood. 115:1594–1604. 2010. View Article : Google Scholar : PubMed/NCBI | |
Khan AA, Betel D, Miller ML, Sander C, Leslie CS and Marks DS: Transfection of small RNAs globally perturbs gene regulation by endogenous microRNAs. Nat Biotechnol. 27:549–555. 2009. View Article : Google Scholar : PubMed/NCBI | |
Li Y, Goldenberg A, Wong KC and Zhang Z: A probabilistic approach to explore human miRNA targetome by integrating miRNA-overexpression data and sequence information. Bioinformatics. 30:621–628. 2014. View Article : Google Scholar : PubMed/NCBI | |
Palumbo A and Anderson K: Multiple myeloma. N Engl J Med. 364:1046–1060. 2011. View Article : Google Scholar : PubMed/NCBI | |
Mahindra A, Laubach J, Raje N, Munshi N, Richardson PG and Anderson K: Latest advances and current challenges in the treatment of multiple myeloma. Nat Rev Clin Oncol. 9:135–143. 2012. View Article : Google Scholar : PubMed/NCBI | |
Lohr JG, Stojanov P, Carter SL, Cruz-Gordillo P, Lawrence MS, Auclair D, Sougnez C, Knoechel B, Gould J, Saksena G, et al Multiple Myeloma Research Consortium, : Widespread genetic heterogeneity in multiple myeloma: Implications for targeted therapy. Cancer Cell. 25:91–101. 2014. View Article : Google Scholar : PubMed/NCBI | |
Augustson BM, Begum G, Dunn JA, Barth NJ, Davies F, Morgan G, Behrens J, Smith A, Child JA and Drayson MT: Early mortality after diagnosis of multiple myeloma: Analysis of patients entered onto the United Kingdom Medical Research Council trials between 1980 and 2002 - Medical Research Council Adult Leukaemia Working Party. J Clin Oncol. 23:9219–9226. 2005. View Article : Google Scholar : PubMed/NCBI | |
Lu J, Getz G, Miska EA, Alvarez-Saavedra E, Lamb J, Peck D, Sweet-Cordero A, Ebert BL, Mak RH, Ferrando AA, et al: MicroRNA expression profiles classify human cancers. Nature. 435:834–838. 2005. View Article : Google Scholar : PubMed/NCBI | |
Lawrie CH, Soneji S, Marafioti T, Cooper CD, Palazzo S, Paterson JC, Cattan H, Enver T, Mager R, Boultwood J, et al: MicroRNA expression distinguishes between germinal center B cell-like and activated B cell-like subtypes of diffuse large B cell lymphoma. Int J Cancer. 121:1156–1161. 2007. View Article : Google Scholar : PubMed/NCBI | |
Munker R, Liu CG, Taccioli C, Alder H and Heerema N: MicroRNA profiles of drug-resistant myeloma cell lines. Acta Haematol. 123:201–204. 2010. View Article : Google Scholar : PubMed/NCBI | |
Ma J, Liu S and Wang Y: MicroRNA-21 and multiple myeloma: Small molecule and big function. Med Oncol. 31:942014. View Article : Google Scholar : PubMed/NCBI | |
Huang JJ, Yu J, Li JY, Liu YT and Zhong RQ: Circulating microRNA expression is associated with genetic subtype and survival of multiple myeloma. Med Oncol. 29:2402–2408. 2012. View Article : Google Scholar : PubMed/NCBI | |
Lu Y, Xiao J, Lin H, Bai Y, Luo X, Wang Z and Yang B: A single anti-microRNA antisense oligodeoxyribonucleotide (AMO) targeting multiple microRNAs offers an improved approach for microRNA interference. Nucleic Acids Res. 37:e242009. View Article : Google Scholar : PubMed/NCBI | |
Mulligan G, Lichter DI, Di Bacco A, Blakemore SJ, Berger A, Koenig E, Bernard H, Trepicchio W, Li B, Neuwirth R, et al: Mutation of NRAS but not KRAS significantly reduces myeloma sensitivity to single-agent bortezomib therapy. Blood. 123:632–639. 2014. View Article : Google Scholar : PubMed/NCBI | |
van Andel H, Kocemba KA, de Haan-Kramer A, Mellink CH, Piwowar M, Broijl A, van Duin M, Sonneveld P, Maurice MM, Kersten MJ, et al: Loss of CYLD expression unleashes Wnt signaling in multiple myeloma and is associated with aggressive disease. Oncogene. 36:2105–2115. 2017. View Article : Google Scholar : PubMed/NCBI | |
Trompeter HI, Abbad H, Iwaniuk KM, Hafner M, Renwick N, Tuschl T, Schira J, Müller HW and Wernet P: MicroRNAs MiR-17, MiR-20a, and MiR-106b act in concert to modulate E2F activity on cell cycle arrest during neuronal lineage differentiation of USSC. PLoS One. 6:e161382011. View Article : Google Scholar : PubMed/NCBI | |
Kim YJ, Hwang SJ, Bae YC and Jung JS: MiR-21 regulates adipogenic differentiation through the modulation of TGF-β signaling in mesenchymal stem cells derived from human adipose tissue. Stem cells. 27:3093–3102. 2009.PubMed/NCBI | |
He M, Wang QY, Yin QQ, Tang J, Lu Y, Zhou CX, Duan CW, Hong DL, Tanaka T, Chen GQ, et al: HIF-1α downregulates miR-17/20a directly targeting p21 and STAT3: a role in myeloid leukemic cell differentiation. Cell Death Differ. 20:408–418. 2013. View Article : Google Scholar : PubMed/NCBI |