Comprehensive Landscape of Modules-Dysregulated Functions Reveals a Profound Role of ceRNAs in Coronary Heart Disease

Coronary heart disease (CHD) is one of the most common severe cardiovascular diseases. Competitive endogenous RNAs (ceRNA) play critical roles in complex diseases. However, our understanding of the dysregulated functions of ceRNAs in CHD remains limited. Here, we systematically analyzed the alterations of ceRNAs and identified the specific functions based on dysregulated modules from the ceRNA network. A total of 2457 significantly differential expressed genes and 212 differential expressed lncRNAs were identified. We got 76679 regulator relationship between different expression genes and miRNAs and 336 regulator relationship between differential expressed lncRNAs and miRNAs. We constructed the ceRNA network and selected five dysregulated modules. Furthermore, CHD specific functions based on dysregulated modules from the ceRNA network were identified, including histone acetylation, platelet degranulation, cAMP-dependent protein kinase complex, xenobiotic transport and so on. Our results will provide novel insight for a better understanding of the mechanism of ceRNAs and facilitate the identification of novel diagnostic and therapeutic biomarkers in CHD.


Introduction
Coronary heart disease (CHD) is one of the most common severe cardiovascular diseases [1]. CHD results from the buildup of plaques in the coronary arteries, and the ruptured plaques may induce thrombosis in coronary atherosclerosis [2,3]. Despite remarkable success in the surgical and medical management of CHD, many interventions are palliative rather than curative, and some survivors still have signifcant residual hemodynamic and electrical conduction abnormalities and experience cardiovascular complications over the long term [4,5]. Te genetic susceptibility to CHD has aroused much attention. Studies have shown that more than 400 genes are associated with the incidence, pathogenesis, and progression of CHD [6,7]. MicroRNAs (miRNAs) are a type of small noncoding RNA composed of 21-22 nucleotides [8,9]. Tey exert their biological efects by silencing genes posttranscriptionally via binding to miRNA response elements (MREs) in the target mRNA [10,11]. Te report suggested that miR-451 is upregulated in CHD and is associated with the PI3K-Akt-mTOR pathway; the data indicated that miR-451 might be a novel biomarker for CHD [12]. Long noncoding RNAs (lncRNAs) defned as nonprotein-coding RNAs longer than 200 nucleotides [13]. Some studies verifed that the abnormal expression of lncRNA in CHD [14][15][16]. Liu et al. reported that lncRNA-ANRIL is expressed low in the serum of patients with CHD, and it has high predictive value both for efective treatment and for poor prognosis of them. lncRNA-ANRIL is also an independent risk factor for patients poor prognosis [17,18]. Competing endogenous RNAs (ceRNAs) are transcripts that can regulate each other at the posttranscription level by competing for shared miRNAs [19,20]. Numerous studies have highlighted that lncRNAs can bind to miRNA sites as ceRNAs, thereby afecting and regulating the expression of mRNAs and target genes [21,22]. One study observed that the direct binding between lncRNA-MEG3 and miR-26a was confrmed via dual-luciferase reporter assay, which indicated that lncRNA-MEG3 could sponge miR-26a as a ceRNA in CHD [23]. Te other study provided evidence that lncRNA HIF1A-AS1 could act as ceRNA to adsorb miR-204 to suppress miR-204 expression and elevate SOCS2 expression in a mouse model, which was established by left coronary artery occlusion [24].
In the present study, based on the ceRNA hypothesis, we constructed a ceRNA network focused on CHD mechanisms and therapy opportunities by utilizing associated miRNA, lncRNA, and mRNA expression profles. We constructed a ceRNA regulator network among lncRNA-miRNA-mRNA after acquiring the interaction relationship between miRNA and dysregulated genes and miRNA and dysregulated lncRNA. Drug treatment could infuence the molecular network of disease mechanisms. Terefore, we constructed a ceRNA and drug-target interaction network based on the interaction network. Five functionally-related modules were identifed from the ceRNA network. Dysregulated modules based on the ceRNA mechanism play an important role of controlling CHD-specifc biological functions. Drugs involved in modules, which target multiple genes and are regulated by more than one miRNA and lncRNA, may become promising therapeutic targets. We presented a landscape of dysregulated ceRNAs acting as a valuable resource for the investigation of the relationship between ceRNA and CHD. Flowchart of the strategy to identify modules-based CHD specifc functions is shown in Figure 1.

Expression Dataset of Coronary Heart
Disease. Gene expression and lncRNA expression data were obtained from the exoRBase database (http://www exoRBase org) [25], which is a repository of lncRNA and mRNA derived from RNA-seq data analyses of human blood exosomes. All collected datasets were processed through a consistent bioinformatics pipeline. Te expression profles included a total of 38 samples, consisting of 6 disease samples and 32 normal samples. RNA-seq data standardized by TPM, expression profle transformed by log(x+1) for analysis.

Identifcation of Epigenetically Dysregulated lncRNAs and
PCGs. Based on the expression data from the database, t-test and fold change [26] were applied to identify the diferent expression genes (DEGs) and diferent expression lncRNAs (DELs) among CHD and normal samples. Te fold change value was calculated by the average expression of genes or lncRNAS in the disease sample divided by the average gene of genes or lncRNAS in the normal sample. After signifcance analysis and false discovery rate (FDR) adjusting, FDR <0.05 was considered as signifcant diference between DEGs and DELs.

Identifying ceRNA Regulatory
Relationships. Tis report sets up ceRNA regulatory relationships in two steps. Firstly, we downloaded experimentally confrmed interaction relationship between miRNA and lncRNA, miRNA and target gene from starBase database (http://starbase sysu.edu cn/ [27]), selected regulatory relationship between miRNA and DELs and miRNA and DEGs. Tese two types of relations performed hypergeometric to identify shared miRNA. Te specifc calculation formula is as follows: Where m represents the total number of human miRNAs, t represents the number of miRNAs that interact with mRNA, n represents the number of miRNAs that interact with lncRNA, r represents the number of miRNAs shared mRNA and circRNA. By performing the hypergeometric test on the pairing relationship between mRNA and miRNA and the relationship between lncRNA and miRNA, we restricted FDR <0.05 and received a preliminary ceRNA relationship.
Secondly, studies have found that overexpression of lncRNA reduces the number of free miRNA molecules, which in turn leads to gene high expression which are targeted by miRNA. Coexpression is a common method for determining regulatory relationships. We calculated the Pearson's correlation, thus getting the corresponding correlation coefcient cor values of mRNA and lncRNA in normal samples and disease samples, respectively. Considering the statistical signifcance of the algorithm and the size of the network, we believed that the coefcient cor between mRNA and lncRNA in disease and normal samples is greater than 0.5 and the p value smaller than 0.05, which corresponds to the coexpression relationship. Te mathematical formula is as follows: |cor(case) − cor(control)|0.5. (2)

Identifying Dysregulated Modules of ceRNA Network.
We constructed a ceRNA regulatory network based on the interaction data between miRNA and lncRNA, miRNA, and mRNA. Te Cytoscape software was used for the construction of networks. Based on the ceRNA network, we identifed modules using GraphWeb, a web server tool for identifying the network-based biomarkers that best represent the properties of the network [28]. Te GraphWeb tools consisted of three component processes: network datasets (to input human protein-protein interaction pairs; ceRNA network in this present study); network algorithm (we used the betweenness centrality clustering method and the default values were set); and network settings (including default edge settings, node settings, and module settings with less than 3 nodes and insignifcant modules hidden). It has recently been reported that the emerging ceRNA networks are involved in drug resistance of disease [29]. In order to explore the interaction between drugs and ceRNA. Furthermore, we downloaded all drug-targeted data from the DrugBank database (https://www drugbank ca/) [30], matched the genes in the ceRNA network and modules with the drug-target genes, and added the drugs that target the genes in the ceRNA network and modules, obtaining a ceRNA network and modules that include drug information. Step 2 Characterization of ceRNA relationship

CHD Specifc Functions Annotation
Step 3 dysregulated modules based on ceRNA network Step 4 CHD specifc functions

Identifying Dysregulated lncRNAs and PCGs in CHD.
We obtained 16225 genes and 1735 lncRNA expression profles from the exoRBase database. Firstly, we used T-test and fold change to identify 2457 DEGs and 212 DELs compared to normal samples. Te ratio of DEGs to all genes is 15% (Figure 2(a)), and the ratio of DELs to all lncRNAs is 12% (Figure 2(b)). For DEGs and DELs, we found that there are diferences in expression levels between CHD and normal samples. Some genes are upregulated in normal samples and downregulated in CHD samples, but other genes are opposed (Figure 2(c) and Figure 2(d)). Signifcantly downregulated DEGs TRIM69 has been afrmed that is associated with heart failure (Figure 2(e)). Tey used robust rank aggregation to derive an integrative transomic score for each annotated gene associated with a heart failure trait, TRIM69 for incident HF with reduced ejection fraction [32]. We also identifed signifcantly upregulated and downregulated DELs (Figure 2(f )). Yan et al. suggested that lncRNA necrosis-related factor (NRF) was shown to be increased in acute myocardial infarction (AMI) patients with heart failure compared with AMI patients without heart failure and had predictive value for the diagnosis of heart failure [33]. Te report also suggested that lncRNA SNHG14 (small nucleolar RNA host gene 14) was validated to sponge both miR-322-5p and miR-384-5p to elevate PCDH17 levels in vivo and in vitro cardiac hypertrophy models. Te subsequent rescue assays revealed that miR-322-5p and miR-384-5p restored SNHG14 depletion-mediated suppression on hypertrophy in Ang-II-induced cardiomyocytes [34].

Characterizing ceRNA Relationship and Constructing
Network in CHD. We obtained 10212 pairs interactions between lncRNA and miRNA, 423975 pairs interactions between genes and miRNA from the starBase database, and selected 336 pairs interactions between DELs and miRNA and 76679 pairs interaction between DEGs and miRNA. We identifed 6334 lncRNA and gene pairs that shared miRNA after applied hypergeometric. Ten, we used coexpression identifed expression-regulated lncRNA and gene pairs. We fnally got an lncRNA-miRNA-mRNA regulator relationship as a candidate ceRNA and constructed a ceRNA network based on the regulator relationship among lncRNA, miRNA, and mRNA. Furthermore, we added drug information into the network based on the DrugBank (Figure 3).
By integrating lncRNAs-miRNA interactions, mRNAs-miRNA interactions, and drug-target relationships, the ceRNA network was constructed. Te network included 1407 edges and 290 nodes, respectively, 13 lncRNAs, 139 miRNAs, 90 mRNAs, and 48 drugs. Some genes are targeted by multiple drugs, such as PDE4D. lncRNA linc00867, with a high degree, may play an important role in the ceRNA network. We identifed that miRNA hsa-miR-224-5p targeted gene MAPKAPK2, lncRNA LINC00665 regulated miRNA hsa-miR-224-5p in CHD patients, and hsa-miR-224-5p was shared between LINC00665 and MAPKAPK2. At the same time, the nodes of LINC00665 and MAPKAPK2 have a relatively high degree in networks. Terefore, ceRNA regulator relationship among LINC00665, hsa-miR-224-5p, and MAPKAPK2 was identifed. From the DrugBank database, we got the small molecule drug staurosporine, that is, a potent protein kinase C inhibitor which enhances cAMPmediated responses in human neuroblastoma cells [35]. Tis present study suggested that staurosporine infuenced ceRNA relationship LINC00665-hsa-miR-224-5p-MAPKAPK2 through targeted gene MAPKAPK2.

Identifying Dysregulated Modules Based on ceRNA Network in CHD.
We identifed fve modules based on the ceRNA network by GraphWeb (Figure 4). In each module, we have identifed many key nodes. Such as in module2 (Figure 4(b)), gene TMX3 downregulated in CHD samples compared to normal samples, FC value equal to 0.41, lncRNA NUTM2A-AS1 FC value equal to 0.47. Downregulated lncRNA LINC00665 FC value equal to 0.40, meanwhile gene PRKG1 FC value is 0.49 in module4 (Figure 4(d)). In module5, linRNA TAPT1-AS1 FC value is 0.25 and gene SOCS5 FC value is 0.41 (Figure 4(e)). Li et al. used whole-exome sequencing (WES) to explore pathogenic variants present in a Zhuang family with coronary artery disease, they discovered two pathogenic mutations, PPP2R3A and TMX3 in four members of the CHD group and two members of the high-risk group [36]. Te elements in each module have a close regulatory relationship. In module1, we found the regulated relationship among hub nodes OIP5-AS1-has-mir-153-3p-CREBBP, and in module2 we identifed ceRNA relationship NUTM2A-AS1-hsa-mir-107-PDE4D, iloprost is a small molecule drug that was used to treat pulmonary arterial hypertension (PAH) through regulated target gene PDE4D [37].
Our study suggested that the drug iloprost impacted ceRNA among NUTM2A-AS1, hsa-mir-107 and PDE4D. LINC00667 is hub node in ceRNA network, we found some ceRNA relationships including LINC00667 in module3 (Figure 4(c)). In our work, to explored the relationship between ceRNA network and drug, we identifed functionally modules (ceRNA-drug) to explored the relationship between ceRNA network and drug. For instance, module4 contains various ceRNA interactions, involving the hub gene MAPKAPK2. Many genes, such as CMPK1, PDE4D, and MAPKAPK2, are targeted by multidrugs. Moreover, other studies have reported that the PDE4D gene interacts with serum triglyceride, and the haplotypic association found in the present population is indicative of the populationspecifc risk associated with CHD [38]. Tus, we considered using drugs that targeted the gene PED2D for the treatment of CHD patients.

Dysregulated Modules of the ceRNAs Network
Contributing to CHD Specifc Functions. We performed function enrichment analysis based on the whole ceRNA network and fve dysregulated modules in CHD ( Figure 5). Myeloid cell diferentiation and histone acetyltransferase have been discerned in the ceRNA network ( Figure 5(a)), meanwhile other reports also suggested ceRNA regulated myeloid cell diferentiation. Cheng et al. reported that ceRNA regulated relationship infuenced histone acetyltransferase binding function [38], confrming the reliability of dysregulated modules in our results. We obtained 10 biological functions from module 1, such as histone modifcation and acetylation, other researches have reported that aberrant histone modifcations, including acetylation and trimethylation, were found both in global histone and specifc MCP-1 gene locos in monocytes from patients with CHD ( Figure 5(b)). Aberrant epigenetic modifcation enzymes expression may be the regulatory mechanism responsible for aberrant histone modifcations [39,40]. In module2, platelet degranulation pathway was identifed, at the same time, other reports suggested anticoagulated peripheral venous blood from 19 patients with stable CHD and 19 normal control subjects was incubated with or without various platelet agonists and analyzed by whole blood fow cytometry, and it turns out that circulating degranulated platelets were increased in patients with CHD compared with control subjects (Figure 5(c)) [41]. Some studies on account of high-throughput gene expression data have also shown that platelet degranulation is associated with CHD [42].
More interestingly, we also identifed several KEGG pathways related to CHD, vascular smooth muscle contraction identifed in network (Figure 6(a)), ovarian steroidogenesis identifed in module2 (Figure 6(b)  ( Figure 6(c)). Some genes are linked with multiple genes in the ceRNA network. Tese hub genes participated in many ceRNA regulator relationships. For instance, the hub gene CREBBP from module1 plays an important role in CHD, and studies have indicated that CITED2, which encodes a CREBBP/EP300 interacting transcriptional modulator of HIF1A and TFAP2, has a causative impact in the development of CHD in humans [43]. In module2, hub gene GNAS interacted with many miRNAs and lnRNAs, GNAS-EDN3 gene loci is associated with hypertension, left ventricular wall thickness, stroke and coronary artery disease based on 29 genome-wide signifcant variants [44].
Terefore, hub genes which are involved in lots of ceRNA relationships may afect many CHD-related specifc biological functions, in turn infuencing the progression of coronary heart disease. Our work conducted ceRAN network-related drugs in CHD, and identifed fve functionallyrelated modules from the ceRNA network. Drugs involved in functionally-related modules, which target multiple genes and are regulated by more than one miRNA and lncRNA, may become promising therapeutic targets.

Discussion
Te majority of patients with CHD are diagnosed at advanced stages and have poor overall survival. Tanks to highthroughput sequencing technology and the rapid development of bioinformatics, we are able to discover the various aberrant expressions of RNAs in disease samples. Unlike classic molecular or cellular biology studies that focus on a specifc molecular interaction, the ceRNA network is constructed to provide a more comprehensive view of the RNA regulatory mechanism during CHD. Current databases have the apparent advantage that they contain multifarious highthroughput data types. Tis enables researchers to discover changes between diferent kinds of RNA and to uncover the regulatory mechanisms between the ceRNAs. Terefore, our study used expression and regulation data from the exoRBase and starBase databases to identify ceRNA relationship pairs related to coronary heart disease, and applied drug-target information from DrugBank to build a ceRNA drug network. Mined modules from the network and identifed functional pathways related to the network and modules.

Conclusions
In conclusion, the present study constructed a novel ceRNA regulator network including some drug information. Te network may be associated in the progression of CHD. We also mined CHD specifc functions and pathways from dysregulated fve modules of the ceRNAs network which may be associated in the progression of CHD. However, the sample size of CHD is not large, and the data types are not sufcient, such as clinical data and prognostic information.     In the future, if more types of data are available, it will be benefcial to the study of the mechanism of coronary heart disease.

Data Availability
Gene expression and lncRNA expression data were obtained from the exoRBase database (http://www exoRBase org).

Conflicts of Interest
Te authors declare no conficts of interest.