Identification and analysis of survival-associated ceRNA triplets in prostate adenocarcinoma.

Prostate adenocarcinoma (PRAD) is one of the most common causes of cancer-associated mortality worldwide. Recent evidence has emphasized the role of competitive endogenous RNAs (ceRNA) in prostate cancer. However, the current understanding of the roles that ceRNAs play in survival-associated PRAD remains in its infancy. In the present study, a PRAD-specific ceRNA network was constructed by integrating long non-coding RNA (lncRNA)-microRNA (miRNA)-gene interactions using experimental and computational methods, as well as expression correlations from The Cancer Genome Atlas database. The topological features of the ceRNA network were then analyzed and the PRAD-risk lncRNAs were compared with non-risk lncRNAs within this network. It was revealed that PRAD-risk lncRNAs had a higher degree, closeness and betweenness centrality, but also had the shortest path length. Finally, 42 significant PRAD-survival-associated triplets were identified. Notably, these triplets may form a compacted subnetwork composed of only 25 nodes (5 miRNAs, 4 lncRNAs and 16 genes) and 32 edges, indicating that some nodes were involved in many triplets. Among this subnetwork, mir-21 indicated the highest degree centrality and was demonstrated to exert its oncogenic effects in prostate tumors by downregulating transforming growth factor β receptor 2 (TGFBR2). Two triplets (MIR22HG_hsa-mir-21_TGFBR2 and MIR22HG_hsa-mir-21_BCL2) were finally identified; not only were they significantly associated with PRAD survival but they also had the highest average degree in the identified subnetwork. The results from the present study provide further insights into the understanding of the potential roles and interactions of ceRNA triplets and potential prognosis markers for PRAD. Copyright: © Li et al.

Introduction

Prostate adenocarcinoma (PRAD) is one of the most common causes of cancer-associated mortality worldwide (1,2). PRAD is the most frequent prostate cancer histological subtype worldwide, and is derived from basal cells differentiating into glandular cells (3). The annual morbidity rate of PRAD has grown steadily, and the incidence of this cancer has increased by 14% over the last two decades worldwide (4,5). Approximately 30% of men with the disease develop clinical recurrence, and this highly aggressive form of PRAD can lead to mortality (6). One method of curing PRAD is gene therapy (7). Thus, the identification and analysis of survival-associated biomarkers provides key opportunities for improving the prognosis of patients.

Long non-coding RNAs (lncRNAs) are defined as non-protein coding transcripts >200 nucleotides in length (8), and they play an important role in various types of cancer (9), such as colon cancer, breast cancer and PRAD (10,11). Prostate cancer antigen 3, a well-characterized lncRNA, has been approved by the US Food and Drug Administration for clinical decisions regarding repeat biopsies for prostate cancer (12). lncRNAs, and their associated genes, can be negatively regulated by microRNAs (miRNAs) (13–15). Accumulating evidence has suggested that lncRNAs and protein coding genes could serve as competitive endogenous RNAs (ceRNAs), which exert their decoy activity by recruiting miRNA molecules via base-pairing with miRNA-recognition elements (MREs), which subsequently compete with common miRNAs, thus contributing to tumor development, progression and metastasis (16,17). For example, in multiple myeloma, lncRNA metastasis-associated lung adenocarcinoma transcript 1 acts as an oncogene by sponging miR-509-5p to modulate forkhead box P1 expression (18). Based on the ceRNA mechanism, miRNAs and ceRNAs could influence the expression of one another (19). For instance, phosphatase and tensin homolog pseudogene 1 functions as a ceRNA in order to modulate levels of PTEN by decoying miRNA (miR)-106b and miR-93 in gastric cancer (20). Recently, ceRNA-based computational methods have been used in cancer-associated studies (21,22). Wu et al (23) investigated the miR-133b-mediated lncRNA-mRNA ceRNA network and laid the foundation for future investigation into the role of lncRNAs in colorectal cancer. Another study demonstrated that indirect interactions among networks of ceRNAs may aid in improving responses to cancer therapy and the development of new therapeutic interventions (24). However, the role of ceRNAs in PRAD remains to be fully investigated.

The aim of the present study was to provide further insights into the current understanding of the potential interactions of ceRNA triplets, and provide potential prognostic markers for PRAD.

Materials and methods

Expression of genes, lncRNAs and miRNAs

First, the RNA-sequencing version 2 datasets for PRAD were downloaded from The Cancer Genome Atlas (TCGA) database (http://tcga-data.nci.nih.gov/) (25,26). Raw read counts for each exon were derived from exon quantification files provided by TCGA level three dataset (the calculated expression value of a particular composite exon of a gene, per sample) (27). Then, the reads per kilobases per million reads (RPKM) value of one gene or lncRNA could be obtained using the following formula: RPKM=ECx109/SCxl, in which EC represents mapped read counts on all exons of one gene or lncRNA, SC represents mapped read counts on all exons of a sample, and 1 represents the total length of all exons on a gene or lncRNA. The exon structures of genes and lncRNAs were downloaded from GENCODE (version 14; http://www.gencodegenes.org/). The expression dataset of miRNA (level three) was also obtained from TCGA database (25). Finally, the matched RPKMs of genes, lncRNAs and miRNA expression profiles were extracted from 494 tumor samples and 52 normal samples (25).

Construction of the ceRNA network of PRAD

The miRNA-target gene interactions were downloaded from the mirTarBase (28) and TarBase (29) databases. Experimental and computational methods were used to predict miRNA-lncRNA interactions. First, miRanda (30), TargetScan (31), PITA (32) and RNAhybrid (33) were used to predict candidate miRNA-lncRNA interactions with default parameters. Experimental interactions in starBase (34) and DIANA-LncBase (35) were used to filter candidate miRNA-lncRNA interactions, which provided high-throughput and photoactivatable ribonucleoside-enhanced cross-linking and immunoprecipiation experimental data. Next, the correlations between miRNA and target coding genes/lncRNAs were calculated using Pearson's correlation coefficient (PCC) using R software. Due to the negative regulatory relationship between miRNAs and their targets (genes/lncRNAs), only interactions with a PCC <-0.1 and P<0.01 were used to construct the PRAD ceRNA network.

Comparison of PRAD-risk lncRNAs and non-disease lncRNAs

PRAD-risk lncRNAs were downloaded from LncRNADisease (36), which stores manually corrected associations between diseases and lncRNAs. To compare PRAD-risk lncRNAs and non-disease lncRNAs, the degree, betweenness centrality, closeness centrality and shortest path were analyzed for PRAD-risk lncRNAs and non-disease lncRNAs in the ceRNA network.

Identifying differentially expressed triplets in PRAD

First, the fold-change method was applied to identify differentially expressed miRNAs, mRNAs and lncRNAs between PRAD samples and normal samples with their corresponding RPKM profile. A molecule was considered to be differentially expressed with |log2fold change| >1. Subsequently, differentially expressed triplets (lncRNA-miRNA-mRNA) were extracted from the PRAD ceRNA network; each element of these triplets was differentially expressed in PRAD.

Identifying survival-associated triplets in PRAD

All PRAD-survival-associated triplets were extracted from the identified candidates, and a weighted expression (WE) score was calculated for each differentially expressed triplet for identification. For a triplet i, explnci, expmiri and expgenei were defined as representing the normalized expression value of lncRNA, miRNA and gene in this triplet, respectively. wlnc and wmir represent the weighted score of lncRNA and miRNA expression, respectively. The WE scores of triplets (i) were calculated using the following formula: WEi=wlnc × explnci + expgenei + wmir × expmiri. The expression values of lncRNAs were normalized with the following formula: explnc=[x-min (x)]/[max(x)-min(x)], in which min(x) and max(x) represent the minimum and maximum expression values of the lncRNA, respectively. Using the same method, expmiri and expgenei could be obtained. wlnc was calculated using the mean expression of all genes divided by the mean expression of all lncRNAs. Furthermore, wmir could be calculated by the mean expression of all genes divided by the mean expression of all miRNAs. For differentially expressed triplets, the median WE score was used as a threshold to classify all the PRAD tumor samples into high-risk (WE score >31.69) and low-risk groups (WE score <31.69). A Kaplan-Meier analysis, a common measurement of the fraction of patients living for a certain amount of time between two groups (37), was performed to test for significance between the two groups. The log-rank test was used to evaluate the significance of differences in survival time. Triplets with P<0.05 were identified as significant prognostic biomarkers.

Topological measurement

In a given graph (G), G=(V, E), in which V represents a set of nodes, and E represents a set of edges. The degree (D) is defined as the number of edges that connect to a node. If it is assumed that there are k edges linked to node v, then the degree of node v could be described as: D(v)=k.

Betweenness centrality measures the centrality of the node in a network and is defined as the number of shortest paths from each node to all others that pass through this node; it reflects the amount of control that a node exerts over the interactions of other nodes in the network. The betweenness centrality (BC) of node V is described as: BC(v)=[Σs≠v≠t (σst(v)/σst)], where σst represents the number of shortest paths from node s to node t, and σst(v) represents the number of those paths that pass through node v.

The closeness centrality is how close a node is to other nodes in the network and is defined as the average mean path from this node to all other nodes. The closeness centrality (C) of node v is defined as:

Where d(u,v) is the shortest distance between node u and node v, and n is the number of nodes of the network.

These topological measures were analyzed and visualized using Cytoscape (version 3.2.1) (38).

Results

ceRNA network for PRAD

To construct a PRAD ceRNA network, interactions between miRNAs and target genes were downloaded from TarBase and mirTarBase with experimentally supported datasets, and miRNA-lncRNA interactions were predicted using a computational program and filtered using experimental data. For miRNA-lncRNA interactions, interactions predicted by more than two programs (from either TargetScan, miRanda, PITA or RNAhybrid) were kept as candidates, and then filtered by experimental interactions in starBase and DIANA-LncBase. Following de-redundancy, 43,497 validated miRNA-target pairs and 314,729 miRNA-lncRNA interactions were kept for further analysis. Next, a PCC <-0.1 and P<0.01 were used as a cut-off to filter negatively regulated associations between the miRNAs-genes and lncRNAs. Finally, a ceRNA network for PRAD was constructed, including 2,159 interactions (516 miRNA-gen interactions and 1,643 miRNA-lncRNA interactions) between 210 lncRNAs, 169 miRNAs and 391 genes (Fig. 1). To dissect the topological characteristics of the ceRNA network for PRAD, the degree of node distribution was analyzed within the network. The ceRNA network for PRAD followed the power-law distribution with an R2 of 0.893 and a slope of −1.398 (Fig. 2A). When comparing the three types of node, it was revealed that the degrees of the miRNA nodes were higher than those of the genes and lncRNAs, indicating that miRNAs accounted for the major proportion, whereas genes and lncRNAs accounted for only a small proportion of the ceRNA network (Fig. 2B and C).

Figure 1.

Prostate adenocarcinoma-specific competitive endogenous RNA network. Triangular nodes represent miRNAs, square nodes represent lncRNAs, and circular nodes represent genes. The black circles, with dotted lines represent the lncRNAs, miRNAs and genes that were studied further. miRNA, microRNA; lncRNA, long non-coding RNA.

Figure 2.

Analysis of the topological characteristics of the prostate adenocarcinoma-specific competitive endogenous RNA network. (A) The degree distributions of all nodes. (B) The proportions of gene, lncRNA and miRNA nodes. (C) The degree distributions of gene, lncRNA and miRNA nodes. miRNA, microRNA; lncRNA, long non-coding RNA.

PRAD-risk lncRNAs were obtained from LncRNA Disease. There were six PRAD-risk lncRNAs in the PRAD ceRNA network [H19, LINC00963, maternally expressed gene 3, differentiation antagonizing non-protein coding RNA, PVT1 and taurine upregulated gene 1 (TUG1)]. In order to compare PRAD-risk lncRNAs with non-disease lncRNA, the topological characteristics of degree, betweenness centrality, closeness centrality and shortest path were calculated within the PRAD ceRNA network. It was revealed that PRAD-risk lncRNAs had a higher degree, closeness and betweenness centrality, but a lower shortest path length (Fig. 3). For example, the degree of TUG1, which functions as a tumor suppressor by regulating PTEN expression in prostate cancer (39), ranked fourth among all of the lncRNAs (degree=45; Fig. 2A).

Figure 3.

Topological analysis of PRAD-risk lncRNAs and non-disease lncRNAs. The red and blue bars represent PRAD-risk lncRNAs and non-disease lncRNAs, respectively. A comparison of the (A) shortest path, (B) betweenness centrality, (C) closeness centrality and (D) degree. *P<0.05; **P<0.01. PRAD, prostate adenocarcinoma; lncRNA, long non-coding RNA.

Survival-associated triplets in PRAD

To identify survival-associated triplets, differentially expressed triplets (lncRNA-miRNA-mRNA) were first extracted from the PRAD ceRNA network, meaning that each element of the triplet was differentially expressed. There were 159 differentially expressed triplets. WE scores were calculated for each differentially expressed triplet. All patients with PRAD were classified into a high-risk group or a low-risk group according to the median WE score for each differentially expressed triplet. The Kaplan-Meier analysis was then performed to test the significance of the two groups. Finally, a total of 42 survival-associated triplets were identified (P<0.05; Table I and Fig. 4).

Table I.

A total of 42 survival-associated triplets.

Triplets	P-value
MIR22HG_hsa-mir-210_NCAM1	0.006
MIR22HG_hsa-mir-210_BDNF	0.007
MIR22HG_hsa-mir-20a_TGFBR2	0.031
RP11-444D3.1_hsa-mir-21_TIMP3	0.033
LINC00958_hsa-mir-21_TIMP3	0.033
MIR22HG_hsa-mir-21_TIMP3	0.033
MIR22HG_hsa-mir-20a_CCND2	0.034
RP11-444D3.1_hsa-mir-21_TPM1	0.034
LINC00958_hsa-mir-21_TPM1	0.034
LINC00958_hsa-mir-210_GPD1L	0.034
MIR22HG_hsa-mir-17_TGFBR2	0.034
RP11-444D3.1_hsa-mir-21_TGFBR3	0.035
LINC00958_hsa-mir-21_TGFBR3	0.035
RP11-444D3.1_hsa-mir-21_TP63	0.035
LINC00958_hsa-mir-21_TP63	0.035
LINC00958_hsa-mir-21_SOX5	0.035
LINC00958_hsa-mir-21_PRRG4	0.035
MIR22HG_hsa-mir-21_TPM1	0.035
RP11-444D3.1_hsa-mir-21_BCL2	0.035
LINC00958_hsa-mir-21_BCL2	0.035
MIR22HG_hsa-mir-20a_BCL2	0.035
LINC00958_hsa-mir-21_SERPINB5	0.035
RP11-444D3.1_hsa-mir-21_SERPINB5	0.035
MIR22HG_hsa-mir-17_SELE	0.036
MIR22HG_hsa-mir-17_CCND2	0.037
MIR22HG_hsa-mir-21_TGFBR3	0.038
MIR22HG_hsa-mir-21_TP63	0.038
RP11-444D3.1_hsa-mir-21_SOX5	0.038
MIR22HG_hsa-mir-21_SOX5	0.038
MIR22HG_hsa-mir-21_PRRG4	0.038
RP11-444D3.1_hsa-mir-21_PRRG4	0.038
MIR22HG_hsa-mir-21_TGFBR2	0.038
RP11-444D3.1_hsa-mir-21_TGFBR2	0.038
LINC00958_hsa-mir-21_TGFBR2	0.038
MIR22HG_hsa-mir-21_CFL2	0.038
RP11-444D3.1_hsa-mir-21_CFL2	0.038
LINC00958_hsa-mir-21_CFL2	0.038
MIR22HG_hsa-mir-21_BCL2	0.038
MIR22HG_hsa-mir-21_SERPINB5	0.038
MIR22HG_hsa-mir-17_BCL2	0.042
LINC00968_hsa-mir-197_ABCC3	0.044
LINC00958_hsa-mir-210_NCAM1	0.048

miR, microRNA.

Figure 4.

A total of 42 survival-associated triplets. The X-axis represents the-log(P-value) from the Kaplan-Meier analysis.

Of these triplets, it was revealed that some genes, miRNAs or lncRNAs were involved in a number of the triplets identified. These triplets were then extracted and a survival-associated sub-network was constructed. Notably, these 42 survival-associated triplets were composed of only 25 nodes (5 miRNAs, 4 lncRNAs and 16 genes) and 32 edges (Fig. 5A). The miRNA with the highest degree was miR-21. The lncRNAs and genes with most connections were MIR22HG, and transforming growth factor-β receptor 2 (TGFBR2) and BCL2, respectively. Thus, two triplets were identified (MIR22HG_hsa-mir-21_TGFBR2 and MIR22HG_hsa-mir-21_BCL2). The survival curves of these triplets are presented in Fig. 5B and C, respectively. The red curve represents samples with a higher WE score and the blue curve represents samples with a lower WE score.

Figure 5.

Survival-associated subnetwork and Kaplan-Meier curves. (A) Survival-associated subnetwork constructed from the 42 triplets. Triangular nodes represent microRNAs, square nodes represent long non-coding RNAs and circular nodes represent genes. The Kaplan-Meier survival curves of (B) MIR22HG_hsa-mir-21_TGFBR2 and (C) MIR22HG_hsa-mir-21_BCL2 are presented. The red curve represents samples with a higher WE score. The blue curve represents the samples with a lower WE score. The numbers of patients included in the high and low WE score groups are both 247. TGFBR2, transforming growth factor β receptor 2. miR, microRNA; WE, weighted expression.

Discussion

The ceRNA hypothesis has been reported to represent a novel post-transcriptional layer of gene regulation operating via miRNA competition (16). With the crosstalk of ceRNAs, it has been demonstrated that miRNAs, and their ceRNA targets, can connect directly or indirectly to form a complex ceRNA network (40). In the present study, a PRAD-specific ceRNA network was constructed by integrating experimentally validated and computationally predicted miRNA-targeted gene/lncRNA interactions, as well as the negative expression correlations between miRNA-target genes and lncRNAs. The ceRNA network for PRAD contained 2,159 interactions (516 miRNA-gen interactions and 1,643 miRNA-lncRNA interactions) between 210 lncRNAs, 169 miRNAs and 391 genes. The ceRNA network followed the power-law distribution with an R2 of 0.893 and a slope of −1.398 (Fig. 2A). lncRNAs with the largest degree value indicates it has an important role in PRAD. In the present study, the average degree, shortest path, betweenness centrality and closeness centrality between PRAD-risk lncRNAs and non-disease lncRNAs were compared. PRAD-risk lncRNAs demonstrated a higher degree, closeness centrality and betweenness centrality, but a lower shortest path length, suggesting that they were more important in communication and the diffusion of information than non-disease nodes in the ceRNA network. Next, differentially expressed lncRNA-miRNA-gene triplets were identified as candidates and PRAD survival-associated triplets were extracted. Finally, 42 PRAD-survival-associated triplets were obtained. Notably, these triplets constructed a compact network composed of only 25 nodes and 32 edges, indicating that some nodes were included in a number of triplets. It was demonstrated that there were only four lncRNAs (LINC00968, MIR22HG, RP11-444D3.1 and LINC00958) within this subnetwork. To the best of our knowledge, there have been no previous reports of these lncRNAs and their association with PRAD, although they have been previously associated with different types of cancer. LINC00968 and MIR22HG were revealed to be differentially expressed in lung squamous cell carcinoma (41,42), while LINC00958 was previously identified as a candidate oncogene in bladder cancer (43).

The miRNA with the highest degree was miR-21, which may target and inhibit the tumor suppressor gene PTEN to promote the proliferation and invasion of prostate cancer cells (44). The genes with the highest degrees were TGFBR2 and BCL2, which are apoptosis regulators associated with many different types of cancer. The TGF-β family serves a fundamental role in a number of different cellular functions in a developmental, context-dependent and cell type-specific manner (e.g., cell migration, survival, proliferation and differentiation) (45). TGFBR2 has previously been indicated to act as a tumor suppressor gene (46,47).

Finally, two triplets (MIR22HG_hsa-mir-21_TGFBR2 and MIR22HG_hsa-mir-21_BCL2) were identified in the present study that were not only associated with PRAD survival but also had the highest average degree in the sub-network (Fig. 5A). Notably, it was revealed that although these two triplets were significant as a whole, the nodes in these triplets were not significantly associated with PRAD survival alone. It has been reported that upregulation of miR-21 may serve as an independent predictor of progress-free survival in patients with advanced prostate cancer (48,49). In addition, miR-21 could exert its oncogenic effects in prostate tumors by downregulating TGFBR2, thus inhibiting the tumor suppressive activity of the TGF-β pathway (50). The triplet results from the present study complement these previous studies, suggesting that miR-21 may be outcompeted by ceRNAs (MIR22HG, TGFBR2 and BCL2), and play an important role in PRAD. miRNAs have been proposed as promising anticancer therapeutic targets (51). As miRNAs are located in the center of the ceRNA network, it is reasonable to envision the potential of these cancer-associated ceRNAs as therapeutic targets.

The success of the present study can be attributed to two aspects: First, the interactions used to construct the PRAD-specific ceRNA network were supported by both experimental and computational approaches. Furthermore, the negative expression regulatory mechanisms between miRNAs and their targets (lncRNAs/genes) in PRAD were also considered. This provided more accuracy in characterizing the ceRNA regulatory associations with PRAD. Secondly, it is more reasonable to identify survival-associated triplets by integrating information from the lncRNAs, genes, miRNAs and topological features, which could therefore help us to understand the ceRNA interactions and the role of lncRNAs in PRAD.

In summary, survival-associated ceRNA triplets involved in PRAD were identified in the present study by constructing a PRAD-specific ceRNA network. These ceRNA triplets may serve as potential therapeutic targets and prognostic biomarkers for PRAD, and the results provide important insights into the understanding of the potential interactions of ceRNA in PRAD.