Construction of lncRNA-mediated ceRNA network to reveal clinically relevant lncRNA biomarkers in glioblastomas.

Cross-talk between competing endogenous RNAs (ceRNAs) play key roles in tumor development. In this study, we performed exon-level expression profiling on 26 glioblastomas (GBMs) and 6 controls to identify long non-coding RNAs (lncRNAs) of GBM initiation and progression using lncRNA-mediated ceRNA network (LMCN). The mRNA and lncRNA expression data, as well as miRNA-target interactions were firstly collected. Then, we used hypergeometric test to detect the lncRNA-mRNA interactions, followed by the construction of LMCN based on Pearson correlation coefficient. With the goal of investigation of the network organization, degree distribution of LMCN was performed. Next, the synergistic, competing lncRNA modules were identified using jActiveModule plug-in of Cytoscape. Moreover, we implemented the pathway analysis for its mRNAs in the module to explore the functions of significant lncRNAs. Using the criteria of degrees >50, 8 hub genes were identified, including EPB41L4A-AS1, ZRANB2-AS2, XIST, HOTAIR, TRAF3IP2-AS1, TPT1-AS1, PVT1 and DLG1-AS1. Furthermore, 1 synergistic, competitive module was identified. In this module, lncRNAs XIST and PVT1 were also the hubs in the synergistic, competing lncRNA module. Functional annotation demonstrated that 5 pathways were identified, including cytokine-cytokine receptor interaction, neuroactive ligand-receptor interaction, and mTOR signaling pathway. We have successfully identified several hubs (such as XIST and PVT1) and significant pathways (for instance, cytokine-cytokine receptor interaction, and neuroactive ligand-receptor interactions) for GBM via establishing the LMCN. These findings might offer potential biomarkers to early diagnose, and predict GBM prognosis in the future.

Introduction

Glioblastomas (GBM), one of the most aggressive and lethal types of brain tumor, accounts for more than half of all gliomas (1). The GBM prognosis is still very poor, although the standard of care has improved survival (2,3). The dismal outcome makes GBM an urgent subject of cancer research. Thus, detecting suitable biomarkers is urgently needed for GBM diagnosis as well as prognosis.

Well-defined risks for GBM are radiation exposure and certain genetic syndromes (4). For example, mutation of PTEN, EGFR, TP53, PIK3CA, IDH1, NF1, RB1 and PDGFRA have been identified in GBM genome, whereas the functional consequences of many of these alterations are still unknown (5–7). With the development of high-througput technology, it has been well accepted that at least 90% of genome sequences is transcribed into noncoding RNAs (ncRNAs), but <2% of the genome encode proteins (8,9). It is increasingly recognized that the roles of ncRNAs constitute an important component in the most biological processes, in part due to the findings that most of the transcripts do not encode for proteins (10,11). According to the size of ncRNA, they are classfied into small ncRNAs as well as long ncRNAs (lncRNAs). Growing evidence has demonstrated that lncRNAs participate in many kinds of biological pathways, for example, neural lineage commitment, and immune response (12–14). Furthermore, lncRNAs have been demonstrated to be important regulators of disease processes (15,16). Significantly, lncRNAs have been revealed to regulate tumor growth and metastasis, including GBM. For example, Pastori et al indicated that lncRNA HOTAIR plays important roles in glioblastoma through promoting cell proliferation (17). Another study has indicated that the oncogenic lncRNA MALAT1 is related to survival and metastasis in GBM (18). However, the overall study of lncRNAs in brain tumors has lagged behind (19,20). Therefore, the extraction of lncRNA signatures is challenging.

Significantly, competing endogenous RNAs (ceRNAs) have been reported to be crucial in illuminating how lncRNAs regulate the expression of coding genes (21,22). For instance, a previous study has reported that the ceRNA interaction network for GBM reveals oncogenic pathways (23). Moreover, Zhang et al (24) have proposed the biological functions of lncRNA associated-ceRNAs in GBM. Collectively, the dysregulated ceRNA network might offer new hope for revealing the pathogenesis underlying GBM.

Herein, in order to better understand the pathogenic processes of GBM, we sought to identify the candidate lncRNA signatures through establishing a functional lncRNA-mediated ceRNA network (LMCN) for GBM.

Materials and methods

The specific steps of our study were divided into the following procedures: Firstly, mRNA as well as lncRNA expression data of GBM were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo/) based on the platform of exon-array data on the Affymetrix Human 1.0 ST Array. Simultaneously, miRNA-target interactions were downloaded from starBase v2.0 (25). Then, a hypergeometric test was applied to identify the competing lncRNA-mRNA interactions, following by co-expression analysis to extract the co-expressed lncRNA-mRNA pairs to further establish the LMCN. Functional analyses for the mRNAs in the LMCN to reveal the biological roles of lncRNAs.

Collecting miRNA-target interactions and affymetrix microarray

The mRNA as well as lncRNA data of GBM (accession number GSE9385) (26) were received from the GEO database, which is on the basis of the (GPL5188) platform of GeneChip Human Exon 1.0 ST Array. The GSE9385 included 26 GBMs, 22 oligodendrogliomas (ODs) as well as 6 control samples. In our study, with the goal of better understanding the molecular mechanisms of GBM, we only selected 26 GBMs and 6 control samples for subsequent analysis. After the probes were mapped to the gene symbols, we obtained a total of 35,367 genes for the following study.

Identification of miRNA-target interactions

starBase v2.0 provides high-quality and experimentally validated miRNA-target interactions. In the present study, the experimentally confirmed lncRNA-miRNA as well as miRNA-mRNA interactions were obtained from starBase 2.0. Then, we aligned the obtained 35,367 genes to the lncRNA-miRNA as well as miRNA-mRNA interactions obtained from starBase 2.0, and we derived the new expression profiles of 11,423 genes (including 212 lncRNAs and 11,211 mRNAs). Moreover, we obtained the interactions covering any genes of 212 lncRNAs and 11,211 mRNAs from the lncRNA-miRNA and miRNA-mRNA interactions, respectively. Ultimately, 2551 lncRNA-miRNA interactions and 352,647 miRNA-mRNA interactions were obtained.

Identifying potential ceRNA interactions

Herein, to get the competing lncRNA-mRNA interactions, the hypergeometric test was used. In an attempt to measure the enrichment significance of the common miRNAs, we computed the P-values. Subsequently, false discovery rate (FDR) was utilized to correct the original P-values based on Benjamini & Hochberg method (27). The interactions with FDR <0.01 was considered as the significant and competing interactions.

Co-expression analysis

With an attempt of measuring the co-expression probability of lncRNA-mRNA interactions, we calculated the PCC relying on the expression of the potential lncRNA-mRNAs pair. In this study, the PCC absolute value of an edge was determined as the weight values, and only edges having weight value >1.15 were reserved to establish the LMCN based on Pearson correlation coefficient. Cytoscape (28) was employed to illustrate the highly competitive LMCN.

Topological analysis for LMCN

As documented, topological centralities (degree, closeness as well as betweenness) are broadly utilized to reveal the properties of network (29). Among these centralities, degree is the simplest index. Herein, we analyzed the degree distribution of all nodes in LMCN and we determined the nodes having degrees greater than 50 as hubs.

Mining synergistic, competing lncRNA modules from the LMCN

Although LMCN can supply an overall view of all lncRNA-mRNA interactions, the sub-networks or clusters reveal a more detailed picture of how lncRNAs synergized with the competing mRNAs. Thus, in this study, the jActiveModule plug-in of Cytoscape was employed to further screened out the synergistic competing modules from the LMCN as previously described (30).

Analysis of functional annotation of synergistic, competing lncRNA modules

DAVID is widely used to analyze genes obtained from the genomic experiments. In the present study, we utilized the DAVID to extract significant pathway terms enriched in genes in the module based on the ‘guilt by association’ strategy, thereby revealing the underlying biological process of lncRNAs. We used Fisher's exact test to detect the significant pathways using the criteria of FDR <0.001.

Results

Constructing a highly competitive LMCN

To establish a highly competitive LMCN, the hypergeometric test was first utilized to detect the candidate lncRNA-mRNA interactions via measuring the significance of the common miRNAs among each lncRNA-mRNA interaction. Eventually, using the criteria of FDR <0.01, we identified 212 lncRNAs, 11,199 mRNAs, as well as 218,696 ceRNA interactions. Based on a weight value >1.15, a highly competive LMCN was built, which included 199 lncRNAs, 2052 mRNAs, and 2715 ceRNA interactions. This LMCN is shown in Fig. 1. We observed that the lncRNAs were distributed in the central area of the LMCN, while the mRNAs were distributed in the outside layer (Fig. 1).

Figure 1.

Establishment of a highly competitive LMCN. The LMCN included 199 lncRNAs, 2052 mRNAs, 2715 ceRNA interactions. Pink represents lncRNAs, blue stands for mRNAs. lncRNA, long non-coding RNA; LMCN, lncRNA-mediated competing endogenous RNA network.

Degree distribution of the LMCN

Then, we investigated the degree properties of LMCN to reveal the organization of the LMCN in GBM patients. The degree distribution of the LMCN (R2=0.99964) demonstrated power law distributions (Fig. 2A). In general, the node with a higher degree demonstrated a hub, which regulated more ceRNA interactions. As displayed in Fig. 2B, the degree distribution of the mRNAs was smaller than that of the lncRNAs, which suggested that though lncRNAs did not encode for proteins, these lncRNAs showed more degree distribution than mRNAs within the LMCN. On the basis of the traits of the lncRNAs of the LMCN, some hub nodes were considered as risk lncRNAs for GBM. With the goal of identifying the hub nodes of the LMCN, the degrees of all nodes in the entire network were ordered in a descending rank on the basis of their degree distribution. Using the degrees >50 as the cut-off threshold, there were 8 hub genes, including EPB41L4A-AS1 (degree=164), ZRANB2-AS2 (degree=122), XIST (degree=102), HOTAIR (degree=93), TRAF3IP2-AS1 (degree=61), TPT1-AS1 (degree=53), PVT1 (degree=53), and DLG1-AS1 (degree=52).

Figure 2.

Degree analysis of the LMCN. (A) Degree distributions of the LMCN reflected a power-law distribution. (B) Box-type scatterplot.

Identifying synergistic, competing lncRNA modules from the LMCN

To explore the modularity features of the LMCN, the jActiveModule plug-in of Cytoscape was employed. One synergistic, competitive module was identified, which contained 83 nodes, as listed in Fig. 3. Similar to the LMCN, synergistic, competitive module possessed a scale-free characteristics with power law degree distribution (Fig. 4, R2=0.99431). Based on the network organization, we found that lncRNA XIST competed with 26 mRNAs and 6 lncRNAs (LINC00667, TTTY15, TPT1-AS1, DLG1-AS1, DLEU2 and LINC00839) in the module (Fig. 3), suggesting its significant roles in GBM. Of note, three lncRNAs XIST, TPT1-AS1 and PVT1 in the LMCN were also the hubs in the synergistic, competing lncRNA module.

Figure 3.

A synergistic, competing lncRNA module identified using jActiveModule plug-in of Cytoscape. LMCN, lncRNA-mediated ceRNA network.

Figure 4.

Degree analysis revealing the specific properties of the synergistic, competing lncRNA module. lncRNA, long non-coding RNA.

Analysis of functional annotation

To further investigate the biological significance of lncRNAs in GBM, pathway enrichment analyses of mRNAs in the synergistic, competing lncRNA module were conducted. According to the FDR <0.001, 5 pathways enriched by the mRNAs of the module were identified, including cytokine-cytokine receptor interaction, neuroactive ligand-receptor interaction, and mTOR signaling pathway. Specific information is shown in Table I.

Table I.

The significant pathways relying on FDR<0.001.

Pathway terms	FDR values
Neuroactive ligand-receptor interaction	4.28×10−12
Cytokine-cytokine receptor interaction	4.38×10−05
mTOR signaling pathway	1.99×10−04
Oocyte meiosis	2.40×10−04
Spliceosome	3.63×10−04

FDR, false discovery rate.

Discussion

The ceRNAs have been demonstrated to regulate the expression based on the competing mechanisms, which play important roles in diverse tumor pathological as well as physiological processes. The ‘guilt by association’ method is typically, utilized to explore lncRNA roles (31). Thus, we decided to use a more comprehensive method to establish a functional LMCN, module extraction, and functional analysis. The results suggested that the GBM-related LMCN can be applied to accelerate biosignature identification and therapeutic development for GBM.

In the present study, hub lncRNAs XIST, and PVT1 in the LMCN were also the hubs in the synergistic, competing lncRNA module. The lncRNA XIST is the main regulator of X inactivation in mammals (32). Previous studies have indicated that XIST exerts critical functions in cell proliferation, cell differentiation, and genome maintenance (33). Significantly, a former study has implicated that XIST is overexpressed in glioma tissues and glioblastoma stem cells, and knockdown of XIST plays tumor-suppressive roles via decreasing cell proliferation, and migration as well as inducing apoptosis (34). PVT1 is an oncogene as well as a Myc protein target known to be upregulated in transformed cells (35). As reported, lncRNA PVT1 is related to cell proliferation, lymph node invasion, cell apoptosis, and metastasis (36,37). The aberrant expression of PVT1 has been identified in cancers, including gastric cancer, prostate cancer, and breast cancer (38). However, no knowledge on the relationship between PVT1 and GBM has been reported. Demonstrated here, we infer that regulating the expression of XIST and PVT1 might be of therapeutic benefit in GBM patients.

Neuroactive ligand-receptor interaction was the most significant pathway in this study. Neuroactive ligand-receptor interactions involve in signaling molecules and interactions which exert critical functions in many cellular processes, for example, apoptosis and cell proliferation. Moreover, apoptosis and cell proliferation are key processes for cancer onset and progression (39,40). Cytokine-cytokine receptor interaction is another significant pathway. Cytokines are crucial regulators involved in cell growth, and angiogenesis (41). As reported, cytokines can induce inflammatory responses which are induced through binding to the corresponding receptors of cytokines (42). Significantly, cell growth, angiogenesis, and inflammation contribute to cancer development and progression (43,44). Importantly, cytokine-cytokine receptor interaction has been indicated to play key roles in cancer (45). Remarkably, neuroactive ligand-receptor interactions and cytokine-cytokine receptor interactions have been indicated to be highly associated with GBM (41).

Therefore, we suggested that altered expression of neuroactive ligand-receptor interactions and cytokine-cytokine receptor interactions caused GBM occurrence.

In the present study, computational method of LMCN network was used to identify the potential lncRNAs and significant pathways. This approach potentially helps us move closer to understanding the biological phenomena of GBM. However, some limitations still remained. To begin with, samples were limited. In addition, study was performed using bioinformatics method, but was not validated using experiments. Moreover, the data used in our study were retrieved from the public database, not produced by us. Therefore, it is urgent to develop chip data on GBM and to do animal experiments to confirm our obtained findings.

In short, we identified several hub lncRNAs (such as XIST and PVT1) and significant pathways (for instance neuroactive ligand-receptor interactions and cytokine-cytokine receptor interactions) for GBM based on the establishment of LMCN. These findings might offer potential biomarkers to early diagnose, and predict GBM prognosis.