Identification of functional lncRNAs through constructing a lncRNA-associated ceRNA network in myocardial infarction

1.   Introduction

Myocardial infarction (MI) is a serious kind of coronary heart disease, refers to the ischemic necrosis of the heart muscle ^[1]. On the basis of coronary artery disease, the blood flow in coronary arteries is sharply reduced or interrupted, resulting in severe and persistent acute ischemia of the corresponding myocardium, which eventually leads to ischemic necrosis of the myocardium ^[2]. Common risk factors of MI include age, high blood pressure, smoking, diabetes, high blood lipids, obesity, and inactivity ^[3]. MI can be life-threatening if it is not detected in time or in severe cases. Therefore, it is very important to improve the detection and diagnosis of MI through identifying disease biomarkers and revealing the relevant molecular mechanism and biological function.

In previous studies, long noncoding RNAs (lncRNAs) have been suggested to play critical roles in the physiological and pathological processes of heart. For example, Liu et al. showed that lncRNA CAIF could suppress cardiac autophagy and attenuate MI by targeting p53-mediated myocardin transcription ^[4]. Wang et al. suggested that lncRNA APF could regulate autophagic cell death and MI by targeting miR-188-3p, which may serve as a potential diagnostic target ^[5]. Zhang et al. found that lncRNA ZFAS1 could function by limiting the intracellular level of SERCA2a protein and contribute to the impairment of cardiac contractile function in MI ^[6]. In addition, more researches have shown that lncRNAs could function as microRNA (miRNA) sponges to combine with the miRNA response elements and mediate gene silencing, which is known as competing endogenous RNA (ceRNA) hypothesis ^[7]. The ceRNA relationship between lncRNAs and mRNAs has been considered to play an important role in the research of MI-related regulatory mechanism. For example, based on ceRNA hypothesis, Zhang et al. constructed dysregulated lncRNA-mRNA networks for both MI occurrence and recurrence. Four lncRNAs (RP1-239B22.5, AC135048.13, RP11-4O1.2, RP11-285F7.2) and three lncRNAs (RP11-363E7.4, CTA-29F11.1, RP5-894A10.6) were respectively identified from these two networks for distinguishing controls from patients of MI and considered as potential biomarkers ^[8]. LncRNA KLF3-AS1 acted as ceRNA to sponge miR-138-5p and then regulated Sirt1 for inhibiting cell apoptosis and attenuating MI progression ^[9]. It follows that the critical lncRNAs could be effectively identified and analyzed based on the lncRNA-associated ceRNA network. However, only a few studies focused on identifying lncRNA regulators via biological network, especially via ceRNA network. Song et.al firstly demonstrated that ceRNA network was participated in the pathological processes of cardiovascular diseases, such as cardiac hypertrophy ^[10]. But the study only emphasized the miRNA mediated regulatory mechanisms and ignored the gene expression information. Furthermore, module analysis of ceRNA network in this study was simple. But this study demonstrated the potential value of network analysis in cardiovascular diseases.

In network biology, it is possible to explain topological properties and define functional modules. Such as hub nodes, which are considered as nodes with the biggest degree. Hub nodes represent the important measure of network centrality. A module represents a group of genes taking part in specific function. Both network topological properties and functional modules contribute to comprehensively understand the biological mechanisms of human diseases ^[11]. Several graph-based algorithms have been developed to perform module identification, such as MCODE ^[12], CFinder ^[13] and Markov Cluster ^[14]. These approaches all focus on the detection of densely connected subgraphs of the network, which could be applied to aid in analysis and solve biological problems. In the field of network based investigation of MI, some studies recently have identified crucial genes based on multiple network analysis algorithms, such as WGCA analysis and microRNA-Gene analysis ^[8,15,16]. However, all these methods could identify some novel regulators, the network construction protocol ignored some important information, such as gene expression and upstream-downstream regulation relationships, which could yield the false positive results. Song et.al firstly constructed the mRNA-mRNA ceRNA pairs in cardiovascular disease based on gene expression information and miRNA-target information. Thus, in this study, we integrated miRNA-mRNA/lncRNA regulatory relationships and lncRNA-mRNA co-expression information to construct MI related ceRNA network and reveal the potential disease genes.

In this study, we firstly obtained differentially expressed (DE) mRNAs and lncRNAs based on MI-related gene expression profile. Then, on the basis of ceRNA hypothesis, we obtained significant lncRNA-mRNA ceRNA pairs and constructed a lncRNA-associated ceRNA network in MI. Through extracting hub subnetwork and identifying functional modules from the lncRNA-associated ceRNA network, we performed an in-depth analysis of the network and focused on the network centrality and important local nodes. Finally, random walk algorithm was used for prioritizing potential MI-related crucial lncRNAs. In this manuscript, we have 4 sections, above contents were Section 1, the rest of the paper is organized as follows: Section 2 contains the methodology (method). Section3 contains the results. Section 4 contains the conclusions and policy implications.

2.   Materials and methods

2.1.   Dataset and differentially expressed genes

We downloaded MI-related gene expression profile GSE97320 from Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/). This dataset contained 6 human samples referring to 3 healthy people and 3 patients with acute MI, which were sufficient for differential expression analysis.

For obtaining both mRNA and lncRNA expression profiles with the same samples, we performed probe re-annotation to the dataset. Dataset GSE97320 was conducted on the platform of GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array). Thus, we downloaded the corresponding fasta file of probe sequences from the annotation file (http://www.affymetrix.com). We also downloaded the fasta file of human genome (GRCh38) and the gtf file of annotation file from GENCODE database (https://www.gencodegenes.org). Using BLASTn, we identified probe-matched human lncRNA transcript sequence and mRNA transcript sequence. The alignment results were filtered according to the following criteria: (I) perfect match; (II) specific match; and (III) identified by ≥ 4 probes. If multiple probes mapped to the same lncRNA/mRNA, we assigned the mean value to the lncRNA/mRNA.

Then, based on the mRNA and lncRNA expression profiles, we used SAM to obtain DE mRNAs and lncRNAs with the threshold of 2-fold change (FC) and P-value < 0.05. This result was also considered as the step1 to construct MI ceRNA network (Step1).

2.2.   Construct a lncRNA-associated ceRNA network in MI

On the basis of the ceRNA hypothesis and the above DE mRNAs and DE lncRNAs (Step1), we constructed a lncRNA-associated ceRNA network in MI via following two steps (Figure 1).

Figure 1.  The pipeline of ceRNA network construction for MI.

DownLoad: Full-Size Img PowerPoint

First, we downloaded all the curated 423,975 miRNA-mRNA interactions from starBase, containing 386 miRNAs and 13,861 mRNAs ^[17]. We mapped MI-related DE mRNAs to these miRNA-mRNA interactions for extracting miRNA-DE mRNA interactions. With the miRNA and DE lncRNA sequences as input data, we used the miRanda tools with default parameters for identifying significant miRNA-DE lncRNA interactions.

Then we used hypergeometric test for calculating the statistical significance of the number of the shared miRNAs based on miRNA-DE mRNA interactions and miRNA-DE lncRNA interactions. Specifically, if the number of the shared miRNAs between a lncRNA and a mRNA met the threshold of hypergeometric test p-value < 0.05, we could get an initial lncRNA-mRNA pair. The formula of hypergeometric test could be seen as follows:

where, m is the number of miRNAs in starBase, n is the number of miRNAs in the miRNA-DE lncRNA interactions, t is the number of miRNAs in the miRNA-DE mRNA interactions, and r is the number of the shared miRNAs between a lncRNA and a mRNA (Step2).

Secondly, according to the ceRNA hypothesis, in addition to the number of the shared miRNAs met the threshold of hypergeometric test p-value < 0.05, the lncRNAs and mRNAs that formed ceRNA relationships should be expression correlated. So we further computed Pearson correlation coefficients (PCC) for the above initial lncRNA-mRNA pairs. Those lncRNA-mRNA pairs met the threshold of PCC > 0.6 were considered as significant lncRNA-mRNA ceRNA pairs.

We integrated all the significant lncRNA-mRNA ceRNA pairs and MI-related DE mRNA-DE mRNA pairs from HPRD for constructing a lncRNA-associated ceRNA network in MI (Step3).

2.3.   Extract a hub subnetwork from the lncRNA-associated ceRNA network

In a network, hub nodes are the nodes whose degrees exceed the specified threshold, which measure the centrality of network nodes. Hub nodes are not vulnerable because they have the biggest degree connecting with most other nodes. However, once hub nodes are destroyed, the entire network would be easily destroyed. That is, in our lncRNA-associated ceRNA network, if the lncRNAs/mRNAs that are hub nodes show MI-related abnormal expression, these nodes and the ceRNA network may play important roles in the occurrence and development of MI. Therefore, we chose top 10% lncRNAs and top 5% mRNAs with the bigger degrees together as hub nodes. And then, the interactions between these hub lncRNAs and hub mRNAs were extracted from our lncRNA-associated ceRNA network. Simultaneously, we mapped these hub lncRNAs and hub mRNAs into the above miRNA-DE lncRNA interactions and miRNA-DE mRNA interactions for obtaining miRNA-hub lncRNA pairs and miRNA-hub mRNA pairs, respectively. If the number of the shared miRNAs between a hub lncRNA and a hub mRNA exceeded 25, the miRNAs were retained. Thus, we constructed a hub subnetwork, which contained hub lncRNAs, hub mRNAs, some crucial miRNAs and the interactions between them.

2.4.   Identify MI-related functional modules

Functional modules represent a group of genes that exert specific function, which may contribute to comprehensively understand the biological mechanisms of MI. Thus, in this study, we used the Molecular Complex Detection (MCODE) plug-in in software Cytoscape to identify MI-related functional modules from the lncRNA-associated ceRNA network. The MCODE algorithm is based on graph-theoretical analysis, which clusters a given network by topology for finding densely connected regions ^[12]. The criteria that we used for identifying MI-related functional modules were as follows: MCODE scores > 5, degree cut-off = 2, node score cut-off = 0.2, max depth = 100, and k-score = 2.

2.5.   Detect MI-related crucial lncRNAs using random walk

For detecting MI-related crucial lncRNAs and prioritizing potential MI-related risk factors, we performed random walk with restart (RWR) to our lncRNA-associated ceRNA network. First, we downloaded MI-related disease genes from DisGeNET, which contains the largest available genes associated with various human diseases ^[18]. Then, these known MI-related disease genes were mapped to the lncRNA-associated ceRNA network as seed nodes for RWR, in order to prioritize crucial lncRNAs associated with MI.

The RWR is an iterative algorithm that starts with a source node for given (such as MI-related disease genes) and transfers from a certain node to a random neighbor node, which is as follows:

where, p⁰ is the initial probability vector, p^t (p^t+1) is an iterative vector in which the ith (i + 1th) element represents the probability of the random walker appearing at a certain node in step t (t + 1), r is the restart probability of the random walk at the initial nodes in each step, and W is the probability transition matrix that was obtained from the column-normalized adjacency matrix of the lncRNA-associated ceRNA network.

In this study, we Set 1 to the known MI-related disease genes and 0 to the other genes for the initial probability vector p⁰, and then performed the iteration steps according to the above formula. The iteration would finish when the change between p^t and p^t+1 was lower than 10^-6. After that, we obtained the real score for each lncRNA of the lncRNA-associated ceRNA network. For computing statistical significance of the scores of these lncRNAs, we further performed random sampling for 3000 times with maintaining the degree of network nodes invariable. In every random experiment, we obtained the random score for each lncRNA of the lncRNA-associated ceRNA network. We counted the times that the random score of each lncRNA was larger than the real one. The p-value of the lncRNA was computed by the ratio of the times and 3000. Thus, we computed statistical significance of all the lncRNAs and considered those lncRNAs with p-value < 0.05 as MI-related crucial lncRNAs.

3.   Results

3.1.   Differential expression analysis

In this study, MI-related gene expression profile GSE97320 was downloaded from GEO. Using probe re-annotation technique, we obtained mRNA and lncRNA expression profiles with the same samples (3 healthy people and 3 patients with acute MI), respectively. In order to performing the following in-depth analysis, we computed the correlation between samples by PCC. We found that most of the samples showed strongly correlated in mRNA expression profile or lncRNA expression profile (Figure 2A, D). It demonstrated that our samples were stable and suitable for the subsequent analysis. Then, SAM test with |FC| > 2 & P-value < 0.05 was used to identify 840 DE mRNAs (Figure 2B) and 137 DE lncRNAs (Figure 2E) between 3 healthy people and 3 patients with acute MI, respectively. For further visualizing up-regulated and down-regulated DE mRNAs/lncRNAs, we also provided Volcano Plot (Figure 2C, F). We conducted pathway annotations for the 840 DE mRNAs and found that they were enriched in MI-related pathways, such as Jak-STAT signaling pathway, MAPK signaling pathway and Wnt signaling pathway (Figure 2G). Actually, previous study has shown that JAK/STAT signaling pathway was involved in the onset of MI and could affect left ventricular remodeling after MI ^[19]. Down-regulation of LINC-PINT promoted the expression of miR-208a-3p and suppressed the expression of JUN, which led to the inactivation of the MAPK signaling pathway in MI and thus played a protective role ^[20]. Wnt/β-catenin signaling pathway could be regulated by MiR-34a and further promoted the apoptosis of myocardial cells in the rat model of MI ^[21]. These results showed that our DE mRNAs played important roles in MI-related biological processes and were worth studying better.

Figure 2.  Analysis of gene expression profile and differentially expressed (DE) genes. (A) Most of the samples showed strong correlation between samples by Pearson correlation coefficients (PCC) in mRNA expression profile; (B) SAM test with |FC| > 2 & P-value < 0.05 was used to identify 840 DE mRNAs; (C) Volcano Plot visualized up-regulated and down-regulated DE mRNAs; (D) Most of the samples showed strong correlation between samples by PCC in lncRNA expression profile; (E) SAM test with |FC| > 2 & P-value < 0.05 was used to identify 137 DE lncRNAs; (F) Volcano Plot visualized up-regulated and down-regulated DE lncRNAs; (G) The results of pathway enrichment for the 840 DE mRNAs. The node size represents the number of overlapped genes between DE mRNAs and pathways. The node colour represents the statistical significance of enrichment.

DownLoad: Full-Size Img PowerPoint

3.2.   Construction of the lncRNA-associated ceRNA network in MI

MI-related DE mRNAs and 423,975 miRNA-mRNA interactions from starBase were used to extract miRNA-DE mRNA interactions. The miRanda tools were used to identify significant 181,310 miRNA-DE lncRNA interactions, containing 2589 miRNAs and 110 DE lncRNAs. Then, integrating hypergeometric test p-value < 0.05 and PCC > 0.6, we identified significant lncRNA-mRNA ceRNA pairs. By further integrating MI-related DE mRNA-DE mRNA pairs from HPRD, we finally constructed a lncRNA-associated ceRNA network in MI. This network consisted of 78 lncRNAs, 608 mRNAs and 2433 edges between them (Figure 3A). We then analyzed the degree distribution of the network, which revealed that node degree distribution followed the power law distribution (Figure 3B, R² = 0.83), suggesting the network was a scale free network and some high-degree nodes (defined as hubs) maintained the network architecture. Furthermore, we also extracted the cluster coefficients and average short path length and of the network and compared with 1000 times random networks (Figure 3C, D). Results indicated that cluster coefficient of true network was smaller than random networks. This results implied that function modules were functioned in this network. Average short path length of the true network was larger than random networks, indicating that the network had reduced global efficiency.

Figure 3.  Construction of the lncRNA-associated ceRNA network in MI. (A) Network visualization. Yellow triangle node represents lncRNA. Blue circular node represents mRNA; (B) Degree distribution of the lncRNA-associated ceRNA network; (C) Cluster coefficients of the lncRNA-associated ceRNA network and random networks; (D) Average short path lengths of the lncRNA-associated ceRNA network and random networks.

DownLoad: Full-Size Img PowerPoint

Additionally, we computed the intersection of the 608 mRNAs and 965 MI-related disease genes from DisGeNET and found 137 intersection genes under the threshold of hypergeometric test p-value < 0.01 (Figure 4A). That is, over 20% mRNAs in the network were MI-related disease genes, which suggested that the mRNA nodes of the lncRNA-associated ceRNA network were reliable. We further calculated the degrees of these intersection genes, which were also MI-related disease genes. Results showed that the degrees of these disease gene nodes were significantly larger than the other ones of our lncRNA-associated ceRNA network (Figure 4B). It demonstrated that the MI-related disease genes were the core of the network and played a crucial role in maintaining the stability of the network.

Figure 4.  Network analysis of lncRNA-associated ceRNA network. (A) The intersection of 608 mRNAs in the network and 965 MI-related disease genes from DisGeNET with hypergeometric test p-value < 0.01; (B) The degrees of the given disease genes were significantly larger than the other ones of the lncRNA-associated ceRNA network. (C) Intersection of Top 10 genes with high degree, betweenness and closeness. (D) Visualization of the hub subnetwork. Yellow triangle node represents lncRNA. Green diamond node represents miRNA. Blue circular node represents mRNA; (E) The results of pathway enrichment for the mRNAs of the hub subnetwork. The node size represents the number of overlapped genes between the mRNAs in subnetwork and pathways. The node colour represents the statistical significance of enrichment.

DownLoad: Full-Size Img PowerPoint

3.3.   Extraction of a hub subnetwork

For understanding the overall structural characteristics of the network, three most robust topological features were calculated, including degree, betweenness and closeness. Table 1 and Figure 4C displayed the top 10 nodes with the best topological features. Two lncRNAs, RP11-546D6.3 and RP11-488L18 were considered as the crucial lncRNAs that might control the physiological processes of MI. To reveal the more functional lncRNAs from network analysis, 10 Top 10% lncRNAs and top 5% mRNAs with the bigger degrees were chosen together as hub nodes. All these hub lncRNAs, hub mRNAs, some crucial miRNAs and the interactions between them were extracted to construct a hub subnetwork. The subnetwork was consisted of 7 lncRNAs, 7 miRNAs and 27 mRNAs (Figure 4D). These miRNAs and mRNAs have been validated to be associated with the biological process and molecular function of MI. Such as miR-373-3p, it allowed proliferation and tumorigenesis of primary human cells, which maintained both oncogenic RAS and active wild-type p53 ^[22]. The miR-373-3p was shown to be implicated in the regulation of cell growth, proliferation, apoptosis, migration and invasion, deregulation of which could occur in many cancers ^[23]. The expression of miR-93-5p was validated to have a cardio-protective effect after MI. And miR-93-5p-enhanced adipose-derived stromal cells-induced exosomes could play a protective effect on infarction-induced myocardial damage and prevent cardiac injury by inhibiting inflammatory response and autophagy ^[24]. The ceRNA relationship between miR-302a-3p and SNHG16 resulted in a reduce role in promoting cell proliferation induced by overexpression of SNHG16 ^[25]. Oxidative stress has been related to aging and tumorigenesis via modulating cell cycle. Then, down-regulation of miR-20b-5p and miR-106a-5p under oxidative stress contributed to regulation of G1/S-phase transition of the cell cycle in multipotent stromal cells ^[26]. Constrained expression of miR-20a-5p in SH-SY5Y cells could suppress the expression of ATG7, and thus inhibiting autophagy and cellular proliferation while promoting apoptosis ^[27]. In addition, gene E2F2 was identified as the targets of let-7i-5p, which mediated its effect in regulating cell cycle of cardiomyocyte ^[28]. Other genes including CEBPB, SOCS3, BCL6, LYN, RAC1, and HSPA6 were also reported to be closely associated with the pathology of MI ^{[29,30,31,32,33]}. More interestingly, the 7 lncRNAs we found in the hub subnetwork may become new potential biomarkers for the research of MI-related molecular mechanism. Furthermore, we performed pathway enrichment analysis for the mRNAs of the hub subnetwork. Results showed that they were enriched in some MI-related important pathways, such as ErbB signaling pathway, T cell receptor signaling pathway, natural killer cell mediated cytotoxicity and VEGF signaling pathway and so on (Figure 4E). For instance, several target genes of miR-373 tend to reveal crosstalks between miR-373 and NF-κB, TGF-β, Wnt/β-catenin and JAK/STAT signaling pathways, which were verified as classical tumor-related signaling pathways ^[34,35,36]. These results suggested that hub nodes that located in the core of the network could display very important information from a biological standpoint.

Table 1.  The top 10 nodes with the most important topological features.

Gene	Degree	Gene	Betweenness	Gene	Closeness
RP11-546D6.3	132	RP11-546D6.3	32118.88	RP11-546D6.3	0.0000650660
AC092192.1	128	AC011330.13	28570.82	IKBIP	0.0000650491
RP11-488L18.10	124	RP11-488L18.10	28076.46	JAK1	0.0000650280
AC011330.13	116	AC092192.1	22461.34	YME1L1	0.0000649857
RP11-127B20.2	101	RP11-127B20.2	19713.15	BCL10	0.0000649604
RP11-639B1.1	100	RP11-639B1.1	16209.76	LYN	0.0000649351
AP001171.1	99	CTC-471J1.2	15097.64	MAPK1	0.0000649224
RP11-77H9.8	93	AP001171.1	15015.73	RP11-488L18.10	0.0000649224
CTC-471J1.2	85	RP11-260M2.1	13593.16	DOCK4	0.0000649140
CTD-2528L19.6	85	RP11-77H9.8	9529.205	EMB	0.0000649013

---

 | Show Table

DownLoad: CSV

3.4.   Identification of MI-related functional modules

As we know, modules correspond to predefined gene sets, which were easier to analyze and interpret. Biological processes often demonstrate themselves as coordinated changes across functional modules ^[37]. In this study, MCODE plug-in was used to identify 2 MI-related functional modules from the lncRNA-associated ceRNA network. Module 1 contained 4 lncRNAs and 179 mRNAs (Figure 5A). Module 2 contained 8 lncRNAs and 278 mRNAs (Figure 5C). In a network containing interacting genes, modules are typically identified as the focus of attention. And then, in a module, a group of genes with the greatest connectivity are the core genes that investigate the characteristics of the module. We could see that the 4 lncRNAs in module 1 and 8 lncRNAs in module 2 were the hub nodes at the center of the modules with the greatest degree and connectivity, which may play important roles in MI-related specific function. And the mRNAs as the leaf nodes were the direct targets of these hub lncRNAs, some of which were validated to be related to MI ^[38,39,40]. We further performed pathway enrichment for the mRNAs of the two functional modules, respectively. Results showed that they were enriched in some of the same pathways, such as FoxO signaling pathway, ErbB signaling pathway, Ubiquitin mediated proteolysis, T cell receptor signaling pathway (Figure 5B, D). In addition, module 1 was involved in some specific pathways including VEGF signaling pathway, Jak-STAT signaling pathway, MAPK signaling pathway, WNT signaling pathway, p53 signaling pathway, TGF-beta signaling pathway, GnRH signaling pathway and Focal adhesion (Figure 5B). Module 2 was also referred to some other specific pathways, such as Endocytosis, Axon guidance, Fc epsilon R1 signaling pathway, Insulin signaling pathway, RNA transport, AMPK signaling pathway, mRNA surveillance pathway and Phagosome (Figure 5D). These results showed that different functional modules had both similar and different features. However, all the pathways were associated with the biological process and molecular function of cardiovascular disease related to MI.

Figure 5.  Identification of the functional modules from the lncRNA-associated ceRNA network. (A) Visualization of module 1. Yellow triangle node represents lncRNA. Blue circular node represents mRNA; (B) The results of pathway enrichment for the mRNAs of module 1; (C) Visualization of module 2. Yellow triangle node represents lncRNA. Blue circular node represents mRNA; (D) The results of pathway enrichment for the mRNAs of module 2.

DownLoad: Full-Size Img PowerPoint

3.5.   Detection of MI-related crucial lncRNAs

In our lncRNA-associated ceRNA network, we used random walk algorithm for detecting and prioritizing lncRNAs. After 3,000 times random sampling, we computed statistical significance of the lncRNAs. Thus, 3 MI-related crucial lncRNAs including CTD-3092A11.2, RP5-821D11.7 and CTC-523E23.1 were detected with the threshold of p-value < 0.05 (Figure 6A). Actually, these 3 crucial lncRNAs were also MI-related DE lncRNAs with |FC| > 2 (Figure 6A). CTD-3092A11.2 was up-regulated expression. RP5-821D11.7 and CTC-523E23.1 were down-regulated expression. By mapping these 3 lncRNAs into the lncRNA-associated ceRNA network and further extracting their corresponding miRNAs, we obtained a new network (Figure 6B). This network had 3 lncRNAs, 8 mRNAs, 16 miRNAs and the ceRNA relationships between them. Pathway enrichment analysis from miEAA revealed that this network participated in myocardial infarction via multiple pathways, such as ferroptosis, PPAR signaling, VEGF signaling and so on. Additionally, the results demonstrated that these 3 crucial lncRNAs could act as ceRNAs to play key roles in the biological processes of MI, which could be used as the potential biomarkers in MI.

Figure 6.  Detection of crucial lncRNAs. (A) 3 crucial lncRNAs including CTD-3092A11.2, RP5-821D11.7 and CTC-523E23.1 were detected with p-value < 0.05. They were also MI-related DE lncRNAs with |FC| > 2; (B) Network visualization. Yellow triangle node represents lncRNA. Blue circular node represents mRNA. (C) miEAA enrichment analysis for the miRNAs in (B).

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

MI is a common acute critical disease in clinical practice. With the aging of the society, the acceleration of the pace of modern life, the change of diet and the influence of social and psychological factors, the incidence of MI in China has been presenting an increasing trend year by year. MI occurs due to intense spasm or occlusion of the coronary artery, resulting in severe and prolonged myocardial ischemia or necrosis ^[41]. LncRNAs could act as miRNA sponges by harboring miRNA response elements (MRE) to compete with mRNAs for the shared miRNAs. Dysregulated expression of key lncRNAs could disturb the balance of ceRNA network, thus leading to the occurrence and development of various cancers ^[42]. Many studies have explored the pathogenesis and molecular functions of MI from the perspective of lncRNA-related ceRNAs. For example, Zhang et al. showed that lncRNA H19 could harbor miR-22-3p and regulate the expression of KDM3A gene for improving MI-induced myocardial injury. They revealed the importance of the H19/miR-22-3p/KDM3A ceRNA relationship in MI ^[43]. Liang et al. found lncRNA 2810403D21Rik/Mirf could act as a ceRNA of Mir26a. Down-regulation of 2810403D21Rik/Mirf caused up-regulation of Mir26a, and then elevated autophagy and relieved cardiac injury in MI mice. This study demonstrated that lncRNA 2810403D21Rik/Mirf could function as an anti-autophagic molecule and become a novel therapeutic approach for MI ^[44]. These researches have achieved remarkable results, but the mechanisms underlying the function of lncRNAs in MI have still not been explicitly delineated.

Thus, the present study constructed a lncRNA-associated ceRNA network based on ceRNA hypothesis. Interestingly, we found that over 20% mRNAs in the network were the given disease genes for MI. And the degrees of these given disease genes were significantly larger than the other ones of the network. Then, a hub subnetwork was constructed by integrating hub lncRNAs, hub mRNAs, some crucial miRNAs and the interactions between them. The miRNAs and mRNAs were shown to play important roles in MI-related biological function. Furthermore, 2 MI-related functional modules were identified from the lncRNA-associated ceRNA network. Pathway enrichment results showed that some similar and different pathways were enriched in these functional modules, which were related to cardiovascular disease related to MI. Finally, 3 MI-related crucial lncRNAs, CTD-3092A11.2, RP5-821D11.7 and CTC-523E23.1 were detected as potential biomarkers based on network random walk, which were enriched to the cell death related pathways. These lncRNAs may be the important player of the regulatory network and be involved in a wide variety of MI-related biological progress. We also performed a benchmark analysis for random walk algorithm on three types of network, including ceRNA network, co-expression network and hypergeometric network. Briefly, ceRNA network strategy was used in this study; co-expression network strategy was performed to identify positive correlated lncRNA-mRNA pairs with Pearson correlation coefficient > 0.6; hypergeometric network strategy was performed to identify significant lncRNA-mRNA pairs with hypergeometric p-value < 0.05. We performed 1000, 2000, 3000, 4000, 5000 and 10000 times network permutations to identify potential disease lncRNAs in 3 networks. Results showed that our proposed ceRNA network strategy identified more lncRNAs than other networks in each network permutations (Figure 8). This result also demonstrated that our protocol to identify novel lncRNAs in MI was reliable.

Figure 7.  Expression of 20 functional lncRNAs that identified from this study.

DownLoad: Full-Size Img PowerPoint

Figure 8.  Benchmark analysis of random walk on ceRNA network, co-expression network and hypergeometric network.

DownLoad: Full-Size Img PowerPoint

In a word, we identified MI-related potential biomarkers and functional modules under the background of ceRNA network. Totally, 20 functional lncRNAs were identified and considered as the potential regulators of MI (Figure 6). However, there were some limitations in our study. Firstly, limited data are available in the study. If a large number of MI-related mRNA and lncRNA expression profile data are released, we may discover more valuable information. Secondly, the questions about molecular mechanism of MI are extremely complicated. We only provide new research and diagnosis targets for MI from the ceRNA network perspective. If combined with the experimental research, we will understand the pathogenesis and biological function in depth.

Acknowledgements

We thanks GEO database for data supporting. ZBB designed this project, ZBB, MY and LM processed the data and wrote the manuscript.

Conflict of interest

The authors declare that they have no conflicts of interest to disclose.