2.1 Identification of differential expression of circRNAs, miRNAs, and mRNAs
The volcano plots for the DECs, DEMs, and DEGs and the heat maps for the microarray data were shown in figure 2. The basic information regarding the expression of circRNAs, miRNAs, and mRNAs in the microarray datasets were listed in table 1.
Three microarray datasets (GSE116925, GSE35956, and GSE156508) were included in our study to obtain the differently expressed mRNAs. After the analysis with GEO2R tool, a total of 988 mRNAs were found in gene chip GSE35956 (adj p<0.05 and logFC>|2|); 992 mRNAs were determined in gene chip GSE156508 (adj p<0.05 and logFC>|0.5|); and 908 mRNAs were identified in gene chip GSE116925 (adj p<0.05 and logFC>|1|). Subsequently, we integrated the filtered mRNAs of the three datasets and integrated them with a Venn diagram, and finally, the 104 mRNAs which could be found in at least two data sets were defined as the DEGs (figure 3A).
As for the miRNAs, 320 miRNAs were found in gene chip GSE74209 (adj p<0.05 and logFC>|1|) and 186 miRNAs were found in gene chip GSE93883 (adj p<0.05 and logFC>|1|); then finally the 41 overlapped miRNAs were defined as the DEMs (figure 3B).
As for the circRNAs, we got 10 differentially expressed circRNAs (DECs) from the gene chip GSE161361 (adj p<0.01 and logFC>|4|) and treated them as the target circRNAs (figure 3C). The basic information regarding the target circRNAs was listed in table 2.
2.2 Identification of the target miRNAs and target genes and construction of the circRNA-miRNA-mRNA network
By intersecting predicted target miRNAs and DEMs obtained from GEO, we confirmed 5 target miRNAs finally (figure. 4A).
The predicted mRNAs which could bind to the target miRNAs were obtained by using the miRWalk and the TargetScan website. Then by intersecting the predicted mRNAs and DEGs obtained from GEO, 38 target genes were confirmed after the duplicates were removed (figure. 4B-F).
Finally, the circRNA-miRNA-mRNA network containing circRNA-miRNA pairs and miRNA-mRNA pairs was generated by cytoscape software (figure. 4G), and the heatmap of the 38 target genes included in the network was constructed (figure. 5).
2.3 RNA-seq verification of target genes
To further analyze the relationship between target genes and osteoporosis, we established a model of osteoporosis in ovariectomized rats, and the RNA-seqwas proceeded to analyze the differentially expressed genes (figure. 6 A and B). As shown in figure 6C, the gene expression of MFAP5, CAMK2A, and RGS4 was significantly different between the control and OVX rats. Besides, GO and KEGG enrichment analysis of thedifferentiallyexpressedgenes (DEGs) indicated that the down-regulated genes in the OVX rats were primarily enriched in the PI3K-Akt signaling pathway, endocrine resistance pathway, osteoblast differentiation, and bone development, et.al (Supplement figure 2A and B), while the up-regulated genes in the OVX rats were primarily enriched in the thyroid hormone signaling pathway, and HIF-1 signaling pathway, et.al (Supplement figure 2C and D). The above GO and KEGG items were closely associated withosteoporosis.
2.4 PPI network establishment and hub genes identification of the target genes
In total, after removing nodes not linked to any other nodes, 32 nodes and 74 edges were mapped in the PPI network (figure. 7A). The Cytohubba app in Cytoscape was used to identify hub genes in the PPI network and a significant module containing 5 nodes and 10 edges was identified. These highest-scoring nodes were screened as hub genes: DIRAS2, CAMK2A, MAPK4, CDC42BPA, and RGS4 (figure. 7B). Besides, it's worth noting that the CAMK2A and RGS4 in these hub genes were also significant differences in the RNA sequencing results.
Then, a circRNA-miRNA-hub genesubnetwork based on the five hub genes was constructed (figure. 7C). The binding sites of circRNA-miRNA were shown in supplement figure 1 and table 4, and the binding sites of the miRNA-hub gene were shown in table 5.
2.5 qRT-PCR analysisof relative RNAs
To further examine the relationship between relative RNAs and osteogenesis, we analyzed the correlation between the circRNAs, miRNAs, mRNAs in the ceRNA sub-network and the osteogenesis indicator OPG, in the same 12 cases of human bone tissues. As shown in figure 8, the expression of hsa_circ_0028877, hsa_circ_0082916, hsa_circ_0030712, DIRAS2, CAMK2A, and MAPK4 was significantly correlated with osteogenesis.
Besides, with the correlation analysis of miRNA-circRNA and miRNA-mRNA, we judged that the expression of hsa_circR_0028877-hsa-miR-1273f - CAMK2A and hsa_circR_0028877-hsa-miR-1273f - DIRAS2 conformed to the common regulatory characteristics of the ceRNA mechanism.
2.6 Function enrichment analysis of target genes, target miRNAs, and DECs
The function enrichment analysis of target genes was performed by the Metascape website and ggplot2 package (version 3.3.3). As shown in figure 9A and B,the enrichment analysis of target genes was mainly enriched in the‘MAPK cascade’, ‘cellular response to hormone stimulus’, and ‘regulation of oxidoreductase activity’, et.al. which revealed that the target genes might affect the development of osteoporosis in these signal pathways.
Because the circRNAs belong to a type of novel non-coding RNA, their function was still not being well understood and the functional analysis could not be directly performed on the DECs. Thus, we hypothesized that the host genes of DECs play significant roles in the ceRNA network and the analysis of the host genes could reveal the potential mechanisms of DECs. Then the ggplot2 package was applied for further enrichment analysis. As shown in figure 9C, the enrichment analysis of host genes was mainly enriched in the items of ‘cadherin binding’, ‘rRNA methyltransferase activity’, and ‘cellular response to oxidized low-density lipoprotein’, et.al.
The enrichment analysis of the target miRNAs was performed bythe mirPath v.3 tool in the DIANA Tools website. MiRPath is a miRNA pathway analysis web-server in which the enrichment analysis of miRNAs was performed indirectly through the corresponding target genes. As shown in figure 9D and E, the enrichment analysis was mainly enriched in the items of ‘PI3K-Akt signaling pathway’, ‘Vitamin digestion and absorption’, ‘Thyroid hormone signaling pathway’, and ‘cell death’, et.al.
2.7 Identification of the potential drugs and construction of drug-target gene interaction network
We uploaded 38 target genes from the ceRNA network to the Drugbank, TTD, and CLUE database and obtained 62 kinds of candidate compounds (drugs). Then the drug-target gene interaction network was constructed to visualize their interactions (figure. 10A). Besides, as shown in fig. 10B, to further screen the candidate compounds, the sub-network which contains the drugs that simultaneously target multiple target genes was constructed and 4 drugs were screened out: RKI-1447, FRAX486, Hyaluronic-acid, and Fostamatinib (figure. 10C-F).