Characterization of Dysregulated lncRNA-Associated ceRNA Network Reveals Novel lncRNAs With ceRNA Activity as Epigenetic Diagnostic Biomarkers for Osteoporosis Risk.

The altered expression of long non-coding RNAs (lncRNAs) has been implicated in the development and human diseases. However, functional roles and regulatory mechanisms of lncRNA as competing endogenous RNAs (ceRNAs) in osteoporosis and their potential clinical implication for osteoporosis risk are largely unexplored. In this study, we performed integrated analysis for paired expression profiles and regulatory relationships of dysregulated lncRNAs, mRNAs, and miRNAs based on "ceRNA hypothesis," and constructed an osteoporosis-related dysregulated miRNA-mediated lncRNA-mRNA ceRNA network (DysCeNet) composed of 105 nodes (including eight miRNAs, 24 mRNAs, and 73 lncRNAs) and 515 edges. Functional analysis suggested that the DysCeNet was involved in known osteoporosis or bone metabolism-related biological processes and pathways. Then, we performed random forest-based feature selection for 73 lncRNAs with ceRNA activity and identified 25 of 73 lncRNAs as potential diagnostic biomarkers. A random forest-based classifier composed of 25 lncRNA biomarkers (RF-25lncRNA) was developed for predicting osteoporosis risk. Performance evaluation with the leave-one-out cross-validation (LOOCV) procedure showed that the RF-25lncRNA achieved a good performance in distinguishing high- and low-bone mineral density (BMD) subjects in different osteoporosis datasets. Our study for the first time revealed a global view of lncRNA-associated ceRNA regulation in osteoporosis and provided novel lncRNAs with ceRNA activity as candidate epigenetic diagnostic biomarkers for early detection of osteoporosis risk.

Introduction

Osteoporosis is a progressive systemic skeletal disease with low bone density and deterioration of bone architecture (Garnero, 2017; Tu et al., 2018). It is estimated that over 14 million men and women in the United States will have osteoporosis by 2020 (Burge et al., 2007). People with osteoporosis substantially have an increased risk of bone fragility and fracture, leading to increase pain, disability, nursing home placement, total health care costs, and mortality (Tu et al., 2018). Early detection and intervention for osteoporosis risk are an effective way to delay the development of disease and improve the quality of life of patients. Therefore, there is an urgent need to identify novel molecular biomarkers in the clinical assessment of osteoporosis risk.

With the development and application of high-throughput omics technologies, it has been shown that osteoporosis has genetic and molecular heterogeneity like many other common complex diseases, and disruption in some molecular pathways can disturb the equilibrium of bone turnover and thereby contribute to osteoporosis (Al Anouti et al., 2019). Long non-coding RNAs (lncRNAs) are a newly discovered class of non-coding RNAs (ncRNAs) that are longer than 200 bp (Pelechano and Steinmetz, 2013) and have attracted much attention in recent medical studies. A large number of studies have demonstrated that lncRNAs is an important player of the genomic regulatory network and is involved in a wide variety of biological progress (Mercer and Mattick, 2013; Guo et al., 2016; Quinn and Chang, 2016; Bunch, 2018; Sun et al., 2020). lncRNAs have been attributed to various functions in development, differentiation, and human disease by negatively or positively regulating gene expression at the transcriptional, post-transcriptional, and epigenetic levels (Kornienko et al., 2013; Maass et al., 2014). Growing functional roles of lncRNAs has been highlighted in osteoporosis during recent studies. For example, XIST, a well-known major effector of the X-inactivation process, was recently reported to promote osteoporosis through inhibiting bone marrow mesenchymal stem cell differentiation (Chen et al., 2019). Mei et al. (2019) found that lncRNA ZBTB40-IT1 played a critical role in bone metabolism and can be modulated by osteoporosis GWAS risk SNPs suppresses osteogenesis. It has become increasingly clear that lncRNAs can act as competing endogenous RNAs (ceRNAs) to interact other RNA molecules by competing for binding to shared microRNAs, which has been implicated in development and human disease including osteogenesis (Tay et al., 2014; Huang et al., 2019a, b; Silva et al., 2019). However, genome-wide exploration for miRNA-mediated lncRNA-associated ceRNA mechanism in osteoporosis and their potential clinical implication for osteoporosis risk remained largely unknown.

In this study, we tried to construct a global miRNA-mediated lncRNA–mRNA ceRNA network in the development of osteoporosis by integrating paired expression profiles and regulatory relationship of lncRNAs, mRNAs, and miRNAs based on “ceRNA hypothesis,” and further to uncover novel lncRNAs with ceRNA activity as epigenetic diagnostic biomarkers for identifying people at high risk for developing osteoporosis.

Materials and Methods

Patients and Samples

Ten samples [including five high-bone mineral density (BMD) subjects and five low-BMD subjects] with corresponding transcriptome gene expression microarray data (Affymetrix Human Exon 1.0 ST Array) and epigenomic miRNA microarray data (Affymetrix Multispecies miRNA-2 Array) were obtained from the Gene Expression Omnibus (GEO) database (the accession number is GSE62589[1]). Two other osteoporosis datasets with transcriptome gene expression microarray data were obtained from the GEO database, including GSE56814 dataset (16 high-BMD subjects and 15 low-BMD subjects) (the accession number is GSE56814[2]) and GSE13850 dataset (10 high-BMD subjects and 10 low-BMD subjects) (the accession number is GSE13850[3]).

Acquisition and Analysis of Expression Profiles

Raw transcriptome gene expression microarray data (CEL files) profiled on Affymetrix GeneChip Human Exon 1.0 ST Array (HuEx-1_0-st) and Affymetrix Human Genome U133A Array (HG-U133A) were obtained from the GSE63402 and GSE13850. These raw data were processed and normalized using the Robust Multichip Average (RMA) algorithm of R package “oligo” for background subtraction, quantile normalization, and summarization. Then, all probes of microarray were mapped into protein-coding genes using the R package “biomaRt” (Durinck et al., 2009). LncRNA expression profiles were obtained using repurposing strategy by mapping array probes into the human genome (GRCh 38) and lncRNA annotations from the GENCODE database[4] (Zhou et al., 2016) using SeqMap tool (Jiang and Wong, 2008). Human mature miRNAs were retrieved from the miRNA microarray. Finally, expression profiles of 20,068 mRNA, 7821 lncRNAs, and 1100 miRNA were obtained for further analysis.

Differential expression analysis of lncRNAs, miRNAs, and mRNAs between high- and low-BMD subjects was performed using the R package “limma” (Ritchie et al., 2015). Those lncRNAs, miRNAs, and mRNAs with p < 0.05 were considered as differentially expressed genes. Hierarchical clustering was performed to investigate the expression patterns between high- and low-BMD subjects.

Construction of Dysregulated lncRNA-Associated ceRNA Network

The experimentally validated miRNA–mRNA and miRNA–lncRNA interaction data were collected from the TarBase database[5] (Li et al., 2014). The dysregulated lncRNA-associated ceRNA network (DysCeNet) in osteoporosis was constructed based on the “ceRNA hypothesis” as follows: (i) Pearson correlation coefficient (PCC) was calculated to measure the expression correlation between differentially expressed mRNAs and lncRNAs from matched mRNA and lncRNA expression profiles. Those dysregulated lncRNA–mRNA pairs with PCC > 0.5 were considered as candidate co-dysregulated lncRNA–mRNA ceRNA crosstalk; (ii) expression correlation between differentially expressed miRNAs and differentially expressed mRNAs, and between differentially expressed miRNAs and differentially expressed lncRNAs was evaluated using PCC; (iii) for a candidate co-dysregulated lncRNA–mRNA ceRNA crosstalk, both mRNAs and lncRNAs in this lncRNA–mRNA ceRNA crosstalk are targeted and co-expressed negatively with a common miRNA; this candidate co-dysregulated lncRNA–mRNA ceRNA crosstalk was selected as dysregulated miRNA-mediated lncRNA–mRNA ceRNA crosstalk; (iv) all miRNA-mediated lncRNA–mRNA ceRNA crosstalks were integrated to form a global DysCeNet.

Identification of lncRNA Biomarkers With ceRNA Activity Using a Machine Learning Method

To identify potential lncRNA biomarkers with ceRNA activity, a random forest approach and leave-one-out cross-validation (LOOCV) were used to select optimal lncRNAs biomarkers using the R package “randomForest” and out-of-bag (OOB) error, which measure the performance of the model on the training set (Lv et al., 2019; Tan et al., 2019). The OOB error will produce an unbiased estimate for the classification error, while the bagging method will decrease the chance of overfitting (Toth et al., 2019). Then, a random forest-based classifier was built using the optimal lncRNA biomarkers, and a receiver operating characteristic (ROC) curve and the area under ROC curve (AUC) was used to measure the diagnostic performance of the lncRNA classifier (Lai et al., 2019).

Functional Enrichment Analysis

Functional enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for mRNAs in the DysCeNet was conducted to predict potential biological functions of lncRNAs in the DysCeNet using the R package “clusterprofiler” (Yu et al., 2012). Those significantly enriched GO terms with p < 0.05 with mutually overlapping gene sets were clustered together using the Enrichment Map plugin in Cytoscape environment (Merico et al., 2010).

Results

Identification of Differentially Expressed mRNAs, miRNAs, and lncRNAs Associated With Osteoporosis

To identify potential risk mRNAs, miRNAs, and lncRNAs in osteoporosis, we performed a comparative analysis for expression profiles of mRNAs, miRNAs, and lncRNAs between high- and low-BMD subjects. A total of 68 mRNAs, 11 miRNAs, and 95 lncRNAs were identified as differentially expressed (p < 0.05) in high-BMD subjects compared with low-BMD subjects (Supplementary Table S1). Among them, there were 73 unregulated genes (including six mRNAs, 10 miRNAs, and 57 lncRNAs) and 101 downregulated genes (including 62 mRNAs, one miRNA, and 38 lncRNAs) in high-BMD subjects compared with low-BMD subjects (Figure 1A). Hierarchical clustering analysis showed that expression patterns of these differentially expressed mRNAs, miRNAs, and lncRNAs were capable of distinguishing high-BMD subjects from low-BMD subjects (Figure 1B).

FIGURE 1

Differential expressed analysis between high- and low-BMD subjects. (A) Volcano plot of the distribution of differentially expressed mRNAs, miRNAs, and lncRNAs. (B) Hierarchical clustering heatmap and dendrogram of 10 samples based on differentially expressed mRNAs, miRNAs, and lncRNAs.

Construction and Analysis of Dysregulated lncRNA-Associated ceRNA Network

To construct an osteoporosis-related dysregulated miRNA-mediated lncRNA–mRNA ceRNA network, we performed an integrated analysis for paired expression profiles and regulatory relationships of dysregulated lncRNAs, mRNAs, and miRNAs of 10 samples in the GSE62589 as described in Section “Materials and Methods.” Finally, an osteoporosis-related dysregulated miRNA-mediated lncRNA–mRNA ceRNA network was constructed and was composed of 105 nodes and 515 edges (including eight miRNAs, 24 mRNAs, and 73 lncRNAs) (Figure 2A) (Supplementary Table S2). To further explore the functional implication of the DysCeNet, we performed functional enrichment analysis for mRNAs in the DysCeNet and found that mRNAs in the DysCeNet were significantly enriched in blood vessel development, NIK/NF-kappaB signaling, bone mineralization involved in bone maturation, osteoclast differentiation, and cell aging.

FIGURE 2

Construction and characterization of dysregulated lncRNA-associated ceRNA network (DysCeNet). (A) A global view of the DysCeNet in osteoporosis. (B) The functional enrichment map of GO terms for mRNAs as ceRNA counterparts of lncRNA biomarkers. (C) Enriched KEGG pathways for mRNAs as ceRNA counterparts of lncRNA biomarkers.

Identification of Potential lncRNA Biomarkers for the Osteoporosis Risk

To identify potential lncRNA biomarkers for the osteoporosis risk, we performed feature selection for lncRNAs in the DysCeNet using a random forest model. Finally, 25 lncRNAs of 73 lncRNAs in the DysCeNet were identified as potential biomarkers for the osteoporosis risk according to their discriminative power using OOB error (Table 1). Of them, 10 lncRNA biomarkers are up-regulated and 15 lncRNA biomarkers are down-regulated in high-BMD subjects compared with low-BMD subjects (Figure 3A). To test whether these 25 lncRNA biomarkers could efficiently distinguish high- and low-BMD subjects, we performed hierarchical clustering for 10 samples in the GSE62589 according to the expression pattern of 25 lncRNA biomarkers. As shown in Figure 3B, all 10 samples were classified into two clusters according to the expression pattern of seven lncRNA biomarkers with 100% accuracy. The results of hierarchical clustering demonstrated the potential of 25 lncRNAs as diagnostic biomarkers for osteoporosis risk. Therefore, a random forest-based classifier composed of 25 lncRNA biomarkers was developed. The performance of the RF-25lncRNA for predicting osteoporosis risk was evaluated in the GSE62589 dataset using the LOOCV procedure, in which nine samples were used as the training set and the remaining one was served as the test sample. Results of performance evaluation showed that the RF-25lncRNA achieves a perfect predictive performance in distinguishing high- and low-BMD subjects with an AUC of 1.0 (Figure 3C).

TABLE 1

Genomic information of 25 lncRNA biomarkers with ceRNA activity for osteoporosis risk.

Ensembl id	Gene name	Location	p-valuea
ENSG00000228290	RP11-132M7.3	Chr 6: 85,399,143-85,419,252 (+)	0.0008
ENSG00000258610	RP11-488C13.7	Chr 14: 77,245,295-77,248,592 (+)	0.0004
ENSG00000242009	RP11-359D24.1	Chr 3: 81,889,698-81,972,012 (+)	0.0078
ENSG00000254562	RP11-64I17.1	Chr 11: 38,668,026-38,695,791 (+)	0.0011
ENSG00000224239	AC090044.2	Chr 3: 898,807-899,774 (+)	0.0063
ENSG00000279785	AC008267.8	Chr 7: 66,474,556-66,475,749 (+)	0.0007
ENSG00000264108	RP11-357H3.1	Chr 18: 73,720,073-73,727,816 (−)	0.0027
ENSG00000227885	RP11-79N23.1	Chr 6: 53,605,547-53,611,820 (−)	0.0088
ENSG00000272599	RP11-152N13.16	Chr 10: 74,884,331-74,885,163 (−)	0.0053
ENSG00000236145	AC079988.3	Chr 2: 121,935,903-121,940,156 (+)	0.0061
ENSG00000255471	RP11-736K20.5	Chr 11: 86,603,256-86,636,079 (−)	0.0296
ENSG00000203496	AL133216.1	Chr 10: 38,453,181-38,466,176 (+)	0.0444
ENSG00000250893	AC098869.2	Chr 4: 40,426,119-40,427,585 (+)	0.0299
ENSG00000278990	AL132796.3	Chr 14: 98,916,902-98,917,801 (−)	0.0098
ENSG00000235890	TSPEAR-AS1	Chr 21: 44,506,807-44,516,575 (+)	0.0414
ENSG00000235823	OLMALINC	Chr 10: 100,373,099-100,454,043 (+)	0.0207
ENSG00000253210	AC040970.1	Chr 8: 141,124,717-141,130,064 (−)	0.0281
ENSG00000228705	LINC00659	Chr 20: 62,774,121-62,776,996 (−)	0.0193
ENSG00000236054	Z82246.1	Chr 22: 32,583,300-32,584,204 (+)	0.0097
ENSG00000251329	LINC02353	Chr 4: 32,351,038-32,353,220 (+)	0.0171
ENSG00000234584	AC019186.1	Chr 2: 157,725,708-157,736,005 (−)	0.0104
ENSG00000237017	AC245052.4	Chr 19: 54,119,511-54,125,343 (−)	0.0191
ENSG00000230753	ZNF341-AS1	Chr 20: 33,787,373-33,811,109 (−)	0.0077
ENSG00000240591	AL096701.3	Chr 22: 31,452,881-31,464,304 (+)	0.0179
ENSG00000226649	AC019118.1	Chr 2: 3,156,756-3,157,797 (+)	0.0087

FIGURE 3

Performance evaluation of identified 25 lncRNA biomarkers in the discovery dataset. (A) Expression pattern of 2525 lncRNA biomarkers in high-BMD and low-BMD subjects. (B) Hierarchical clustering heatmap and dendrogram of 10 samples based on expression patterns of 25 lncRNA biomarkers. (C) The ROC curves of RF-25lncRNA in the GSE62589 dataset.

Further Validation of lncRNA Biomarkers for the Osteoporosis Risk in Two Other Independent Datasets

To test the robustness of the lncRNA biomarkers for the osteoporosis risk, 25 lncRNA biomarkers were applied to the independent GSE56814 dataset. The RF-25lncRNA correctly classified 11 of 15 low-BMD subjects and nine of 16 high-BMD subjects, achieving an AUC of 0.733 (95% CI: 0.553–0.913) with an accuracy of 64.5% (Figure 4A).

FIGURE 4

Independent validation of identified 25 lncRNA biomarkers in two independent testing datasets. (A) The ROC curves of RF-25lncRNA in the GSE56814 dataset. (B) The ROC curves of two lncRNA biomarkers (RP11-488C13.7 and RP11-152N13.16) in the GSE13850 dataset.

The performance of the lncRNA biomarkers was further tested in another independent GSE13850 dataset. However, transcriptome gene expression data of the GSE13850 was profiled on the HG-U133A platform, and only two lncRNAs (RP11-488C13.7 and RP11-152N13.16) of 25 lncRNA biomarkers in the RF-25lncRNA were covered on the HG-U133A platform. When two lncRNA biomarkers (RP11-488C13.7 and RP11-152N13.16) were applied to the 20 samples of the GSE13850, two lncRNA biomarkers (RP11-488C13.7 and RP11-152N13.16) correctly classified seven of 10 low-BMD subjects and seven of 10 high-BMD subjects, achieving an AUC of 0.66 (95% CI: 0.396–0.924) with an accuracy of 70% (Figure 4B).

Discussion

Early detection and intervention for osteoporosis are crucial to prevent fragility fractures and delay disease progression. Traditional markers for osteoporosis risk are BMD, vitamin D, alkaline phosphatase, and so on (Parveen et al., 2019). Increasing evidence has suggested that altered molecular profiles contributed to the osteoporosis and outcome, which provided novel insights into molecular basis of osteoporosis and also highlighted the potential of molecular factors as markers for osteoporosis diagnosis and prognosis (Yu et al., 2013; Makitie et al., 2018; Gong et al., 2019; Ukon et al., 2019). Biomarker identification has been proven to be an effective way to recognize people at high risk for developing osteoporosis and have attracted much attention in the clinical decision-making for osteoporosis management. Previous studies have focused on mRNA profiles and miRNA profiles and identified some candidate mRNA- or miRNA-based biomarkers. For example, circulating miRNAs, hsa-miR-122-5p and hsa-miR-4516, have been found to be diagnostic biomarkers for osteoporosis risk (Mandourah et al., 2018). Another study by Seeliger et al. (2014) also identified five freely circulating miRNAs associated with osteoporosis fractures. There is increasing evidence that lncRNAs also play important roles in the pathogenesis of osteoporosis (Peng et al., 2018; Feng et al., 2019; Liu et al., 2019). In this study, we obtained lncRNA profiles of osteoporosis patients by repurposing array probes on publicly available microarray data and performed genome-wide analysis for expression patterns of lncRNAs, miRNAs, and mRNAs between high- and low-BMD subjects. A total of 68 mRNAs, 11 miRNAs, and 95 lncRNAs were found to be associated with osteoporosis, which provided a candidate resource for experimental scientist for studying the molecular mechanism of osteoporosis.

Through many candidate osteoporosis-related mRNAs, miRNAs, and lncRNAs were identified in our study, regulatory relationships and mechanisms among these different types of RNA molecules in the development of osteoporosis are still unknown. It is well known that different RNA molecules can act as ceRNAs to communicate with and co-regulate each other by competing for binding to shared miRNAs (Tay et al., 2014). LncRNAs have been reported as key components of the ceRNA-mediated regulatory network, and aberrant miRNA-mediated lncRNA–mRNA ceRNA crosstalk has been implicated in many human complex diseases, including cancers. However, functional roles and regulatory mechanisms of lncRNA as ceRNAs in the development of osteoporosis and their potential implication for osteoporosis are largely unexplored. To explore the ceRNA activity of lncRNAs in the development of osteoporosis, we performed integrated analysis for paired expression profiles of dysregulated 68 mRNAs, 11 miRNAs, and 95 lncRNAs and miRNA–target regulatory relationship based on “ceRNA hypothesis,” and constructed an osteoporosis-related miRNA-mediated dysregulated lncRNA–mRNA ceRNA network (DysCeNet) composed of 105 nodes and 515 edges (including eight miRNAs, 24 mRNAs, and 73 lncRNAs). Network analysis suggested that a large proportion of deregulated lncRNAs (76.8%, 73/95) in osteoporosis function as ceRNA and communicated with 24 mRNAs by competing for eight common miRNAs (Figure 2A), which implied that extensive variation in miRNAs and lncRNAs disrupted the miRNA-mediated lncRNA–mRNA ceRNA regulatory network contributing to osteoporosis at the post-transcriptional level. Functional analysis through functional enrichment analysis for mRNAs in the DysCeNet found that mRNAs as ceRNA counterparts of lncRNA biomarkers were involved in known osteoporosis or bone metabolism-related biological progression and pathways, including blood vessel development, NIK/NF-kappaB signaling, bone mineralization involved in bone maturation, osteoclast differentiation, and cell aging (Figures 2B,C). For example, blood vessels in the bone play vital roles for the formation of new bone and promoting blood vessel growth could reverse the weakening of bones and treat osteoporosis (Sivaraj and Adams, 2016). NIK/NF-kappaB signaling has been shown to play an important role in the positive and negative regulation of cytokine-mediated osteoclast formation and activation (Boyce et al., 2015). These results suggested that dysregulated expression of lncRNA ceRNAs and the resultant perturbation in miRNA-mediated lncRNA–mRNA crosstalk in the DysCeNet are involved in osteoporosis-biological processes and contributed to the osteoporosis.

A large number of studies have indicated the superior potential of lncRNAs as diagnostic and prognostic biomarkers compared with protein-coding genes due to the fact that lncRNAs were expressed in much more cell- type-, tissue-, and disease-specific patterns that are closely more associated with their function (Hauptman and Glavaè, 2013; Chen et al., 2018). lncRNA biomarkers have been widely investigated and identified in various cancers during the past years (Zhou et al., 2018a, b, 2019; Bao et al., 2019). To explore the potential application of lncRNAs with ceRNA activity as diagnostic biomarkers for osteoporosis risk, we performed random forest-based feature selection for 73 lncRNAs with ceRNA activity and identified 25 of 73 lncRNAs as potential diagnostic biomarkers. To accelerate the clinical application, we also developed a random forest-based classifier (RF-25lncRNA) composed of seven lncRNA biomarkers and test the performance of the RF-25lncRNA in different osteoporosis datasets. Performance evaluation of the LOOCV procedure showed that the RF-25lncRNA achieved a good performance in distinguishing high- and low-BMD subjects in three osteoporosis datasets. These results demonstrated that seven lncRNA with ceRNA activity may become reliable and powerful epigenetic diagnostic biomarkers for early detection of osteoporosis risk.

There are some limitations to this study. First, the performance of newly identified lncRNA biomarkers was validated only in three osteoporosis datasets because of no other publicly available osteoporosis datasets. Second, the biological functions of newly identified lncRNA biomarkers are unknown, although they were found to have ceRNA activity in our study. Therefore, more laboratory and clinical researches were needed.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: Three osteoporosis datasets were obtained from the Gene Expression Omnibus (GEO) database (GSE62589, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62589, GSE56814, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE56814, and GSE13850, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE13850).

Author Contributions

YZ conceived and designed the experiments. MZ and LC analyzed the data. MZ wrote the manuscript. All authors read and approved the final manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.