Exploring prognostic potential of long noncoding RNAs in colorectal cancer based on a competing endogenous RNA network

Zhi-Dong Yang, Hui Kang

Abstract BACKGROUND Colorectal cancer (CRC) is one of the most prevalent tumors worldwide.Recently, long noncoding RNAs (lncRNAs) have been shown to influence tumorigenesis and tumor progression by acting as competing endogenous RNAs(ceRNAs). It is difficult to extract prognostic lncRNAs and useful bioinformation from most ceRNA networks constructed previously.AIM To construct a prognostic related ceRNA regulatory network and lncRNA related signature based on risk score in CRC.METHODS RNA transcriptome profile and clinical information of 506 CRC patients were downloaded from the Cancer Genome Atlas database. R packages and Perl program were used for data processing. Cox regression analysis was used for prognostic model construction. Quantitative real-time polymerase chain reaction was used to detect the expression of lncRNAs.RESULTS A prognostic-related ceRNA network was constructed, including 9 lncRNAs, 44 mRNAs, and 30 miRNAs. In addition, a four-lncRNA model was constructed using multivariate Cox regression analysis, which could be an independent prognostic model in CRC. The risk score for each patient was calculated, and the 506 patients were divided into high and low-risk groups (253 for each group)based on the median risk score. The results of the survival analysis showed that patients with a high-risk score had a poor survival rate. Furthermore, the predictive value of the four-lncRNA model was evaluated in GSE38832. Patient survival probabilities could be better predicted when combing the risk score and clinical features. Gene Set Enrichment Analysis results verified that a number of cancer-related signaling pathways were enriched with a high-risk score in CRC.Finally, we validated a novel lncRNA (LINC00488) using quantitative real-time polymerase chain reaction in 22 paired CRC patient tumor tissues compared to adjacent non-tumor tissues.CONCLUSION The four-lncRNA model could give better predictive value for CRC patients. Our understanding of the lncRNA-related ceRNA regulatory mechanism could provide a potential diagnostic indicator for CRC patients.

Key words: Colorectal cancer; Long noncoding RNA; miRNA; mRNA; Transcription factor; Competing endogenous RNA; Survival

INTRODUCTION

Colorectal cancer (CRC) is one of the three most common cancers worldwide[1].Although there have been improvements in the treatment of this disease, numerous patients are diagnosed at an advanced stage of the disease because of a shortage of useful biomarkers[2]. Advanced stage disease with distant metastasis is the main reason for the poor therapeutic efficacy of surgery in CRC. Therefore, discovering highly effective diagnostic biomarkers and understanding the molecular mechanisms underlying CRC tumorigenesis are needed to significantly increase the survival rate of CRC patients.

Long noncoding RNAs (lncRNAs) are non-coding RNAs of over 200 nucleotides that are located in the nucleus or cytoplasm. Several studies have demonstrated that lncRNAs play a significant role in the development of malignant tumors[3,4]. Salmena et al[5]hypothesized that lncRNAs and mRNAs could interact with each other via microRNAs (miRNAs). This hypothesis was built on the fact that lncRNAs contain miRNA response elements (MREs) and, thus, could competitively bind to the core seed regions of miRNAs to further regulate the expression of target genes. Increasing evidence in recent years has supported this theory by demonstrating that lncRNAs can act as competing endogenous RNAs (ceRNAs), causing decreased miRNA binding to mRNA targets[6]. Indeed, an increasing number of papers have implicated ceRNAs in cancer development[7].

Recently, studies have confirmed that the regulatory networks of lncRNA, miRNA,mRNA, and ceRNA play significant roles in the development of gastric cancer[8],cholangiocarcinoma[9], and CRC[10-12]. However, most of these generated ceRNA networks contain a large number of lncRNAs, and it is difficult to extract the lncRNAs closely related to survival prognosis. In this paper, univariate Cox regression analysis was performed on all the differentially expressed lncRNAs before a ceRNA network was constructed. The RNA expression profiles were combined with patient clinical features to generate the ceRNA network. Through the multivariate Cox regression model, a four-lncRNA model based on the risk score of the ceRNA network was developed. Furthermore, the association between this model, clinical information, and survival prognosis was analyzed. Key lncRNAs in the ceRNA network were validated in clinical CRC tissues by quantitative real-time polymerase chain reaction (qRT-PCR).Our study provided a deep understanding of ceRNAs, and the four-lncRNA model could identify novel candidate biomarkers for CRC.

MATERIALS AND METHODS

Raw data and preprocessing

The raw RNA sequences of 506 colorectal cancer patient samples (534 tumor and 41 non-tumor tissues) were retrieved from the Cancer Genome Atlas (TCGA) database(https://portal.gdc.cancer.gov/). The criteria of inclusion were: (1) Clinical information (age, gender, pathological stage, survival time, and status) should be complete; (2) Sequencing data should be complete, including RNA-Seq data for both mRNAs and lncRNAs; and (3) The follow-up survival time should be ≥ 30 d. Our study followed the TCGA publication guidelines.

Expression profile analysis

We used the Ensemble annotation file (http://www.ensembl.org/index.html) to define and annotate lncRNAs and mRNAs. We identified differentially expressed mRNAs (DEmRNAs), lncRNAs (DElncRNAs), and miRNAs (DEmiRNAs) using program codes written with R and Perl software. A log2 absolute of the fold change ≥2 and P ＜ 0.01 were considered significant. The “heatmap” and “gplots” functions of the R software package were used to generate a volcano map of the differentially expressed RNAs. Functional enrichment analysis was carried out for the DEmRNAs using the “clusterProfiler” package in the R software. Transcription factors (TFs) were identified amongst the clustered gene sets according to the ENCODE CHIP-seq data at the University of California Santa Cruz Genome Browser.

Univariate Cox regression analysis

Univariate Cox regression was performed for all the aberrantly expressed lncRNAs to evaluate the correlation between expression and overall survival (OS). Only DElncRNAs correlated with prognosis were used for the analysis.

ceRNA regulatory network construction and functional annotation

The ceRNA network of CRC was constructed following the steps listed below: (1) The miRcode database (http://www.mircode.org) was used to scientifically retrieve DElncRNA-miRNA interactions and matched DEmiRNA interactions were performed; (2) DEmiRNA-targeted mRNAs were predicted using miRTarBase(http://www.targetscan.org/), miRDB (http://www.mirdb.org/), and TargetScan(http://www.microrna.org/microrna/home.do) databases. Only the mRNAs presented in all three databases and overlapping with the DEmRNAs were used to construct the ceRNA network; and (3) The resulting map was generated using Cytoscape 3.5.1.

Four-lncRNA predictive model construction

A prognostic predictive model based on the risk score was applied using multivariate Cox regression analysis for all the DElncRNAs involved in the ceRNA network. As shown the formula in the picture, where βi stands for the coefficient for each gene and Xi represents the z-score. Kaplan-Meier (K-M) survival analysis and receiver operating characteristic (ROC) curve analysis were used to assess the predictive value of this four-lncRNA model. Univariate and multivariable Cox regression analyses were used to evaluate whether the risk score was an independent predictor when compared with the clinical information.

Four-lncRNA model validation

The raw data was downloaded from the Gene Expression Omnibus (GEO) database.Expression profiles were processed with R software. The risk scores of the patients were calculated using the Four-lncRNA model. K-M and ROC curve analyses were carried out to assess the predictive value of the four-lncRNA model.

qRT-PCR

Total RNA was isolated from 22 paired CRC patient tissues using TRIzol reagent(Invitrogen, CA), cDNA was synthesized using the Transcriptor First Strand cDNA Synthesis Kit (Roche) following the manufacturer’s instructions. qRT-PCR was performed using FastStart Universal SYBR Green Master (ROX) (Roche) on the LightCycler 480. The primers used for target amplification are presented in Table 1.Relative expression was determined using the 2-ΔΔCtmethod.

Risk score combined with clinical feature for CRC prognostic prediction

The prognostic value of American Joint Committee on Cancer (AJCC) tumor invasion depth (T), lymph node metastasis (N), distant metastasis (M), and TNM stage was assessed through ROC curve and survival analyses. In addition, risk score and clinical features were combined, and the prediction value was evaluated using a K-M estimator.

Gene Set Enrichment Analysis

Gene Set Enrichment Analysis (GSEA) v4.0.1 was carried out using the JAVA program with enrichment analysis performed on normalized mRNA expression profiles. Random sample permutations (n = 1000) were performed to calculate the P value. A false discovery rate ＜ 0.01 and a nominal P ＜ 0.05 were considered significant.

Nomogram construction

The risk score and clinical features were combined, and a nomogram was generated using the “rms” package of the R software.

Statistical analysis

R software and relevant R packages were used for all statistical analyses. Survival analysis was performed using the “survival” package. ROC curves were plotted using the “pROC” package. Pearson correlation analysis was carried out to determine correlations among DERNAs.

RESULTS

Cancer-specific RNAs in CRC

According to the cutoff threshold, 823 lncRNAs were differentially expressed (622 upregulated and 201 downregulated) (Figure 1A). We identified 95 lncRNAs that correlated with OS using univariate Cox regression analysis (Supplementary Table S1). In addition, we identified 869 upregulated and 756 downregulated mRNAs.Furthermore, 258 miRNAs (186 upregulated and 72 downregulated) were differentially expressed in the CRC tumors compared to the adjacent non-tumor tissues. The volcano maps of the RNAs are presented in Figure 1B and C. To further understand the molecular mechanisms involved in CRC development, functional enrichment analyses of the Kyoto Encyclopedia of Genes and Genomes pathway and Gene Ontology were performed for the DEmRNAs. The results are presented in Figures 2 and 3. TF-mRNA analysis was performed using a false discovery rate ＜ 0.05.The most significant TFs and their links are shown in Figure 4.

The ceRNA network in CRC

To understand the interactions between the DElncRNAs, DEmRNAs, and DEmiRNAs, 72 paired lncRNA-miRNA interactions (9 lncRNAs and 30 miRNAs)were constructed between the aberrantly expressed lncRNAs and miRNAs through the miRcode database. We found that 1419 mRNAs were targeted by these 30 miRNAs in the miRTarBase, miRDB, and TargetScan databases. The possible interactions between these mRNAs and 1625 DEmRNAs were cross-checked. As a result, 44 mRNAs were involved in the network. The final CRC ceRNA network was constructed by the incorporation of nine DElncRNAs, 44 DEmRNAs, and 30 DEmiRNAs (Figure 5).

The ceRNA network could indicate possible interactions between the differentially expressed RNAs in CRC. To validate this hypothesis, Pearson correlation analysis was carried out on the mRNA expression levels of the RNAs included in the network. The results showed that only lncRNA AL139147.1 had a direct positive linear correlation with FOXQ1 and CBX2, whereas four miRNAs (has-miR-32, has-miR-98, has-miR-141,has-miR-192) had negative linear correlations with the DEmRNAs (R ＞ 0.35, P ＜ 0.01,Figure 6).

Four-lncRNA model for prognosis and survival

In order to identify potential prognostic markers, the expression levels of nine lncRNAs involved in the constructed ceRNA network were measured. The results showed that LINC00488 was downregulated whereas AL139147.1, DSCR4, DSCR8,H19, KCNQ1OT1, LMO7-AS1, MIR31HG, and MYO16-AS1 were upregulated in theCRC tumors compared to the non-tumor tissues (P ＜ 0.05, Figure 7A).

Table 1 Sequences of primers used for polymerase chain reaction

A four-lncRNA multivariate Cox regression model was made for these nine lncRNAs of the ceRNA network. The risk scores were identified for each patient using the following formula: Risk score = (0.20671 ＊ DSCR4) + (-0.27339 ＊ LINC00488) +(0.09205 ＊ MIR31HG) + (0.17811 ＊ LMO7-AS1). The 506 patients were categorized into a high- or low-risk group (253 patients in each group) using the median risk score as the cutoff. Furthermore, the K-M method with the log-rank statistical test was used to evaluate survival between the two groups. We found that the high-risk group had a worse OS (P = 0.003, Figure 7B). Moreover, the area under the ROC curve (AUC) for the model was 0.628 (P ＜ 0.05, Figure 7C), indicating the prognostic value of the fourlncRNA model based on the risk score. The heatmap, risk score distribution, and survival status for the expression profiles are presented in Figure 7D.

Four-lncRNA model validation using the GEO database

To demonstrate the predictive value of the four-lncRNA model for CRC patients, we downloaded GSE38832 from the GEO database, including the gene profiles and clinical features. The risk scores were calculated for 122 patients in this cohort.Interestingly, the survival analysis demonstrated that the patients with a high-risk score had a shorter survival time than those with a low-risk score (P = 0.024; Figure 7E), and the area under the ROC curve was 0.649 (Figure 7F).

Prognostic prediction of the risk score combined with clinical significance for CRC

Correlation analysis between the risk score and patient clinical features revealed that the risk score was associated with tumor invasion depth (T), lymph node metastasis(N), and distant metastasis (M). However, no association was observed with other parameters (e.g., gender and age) (Table 2). Furthermore, univariate and multivariate Cox regression analyses were carried out to assess whether the four-lncRNA model and clinical features correlated with the OS (P ＜ 0.05, Figure 8). Univariate analysis showed that TNM stage (P ＜ 0.001), T (P = 0.01), N (P ＜ 0.001), M (P ＜ 0.001), and risk score of the four-lncRNA model (P ＜ 0.001) were significantly associated with the OS.In addition, we found that the four-lncRNA model could be an independent prognostic factor for CRC patients using multivariate analysis (P = 0.003). The fourlncRNA model risk score with clinical features combined with mRNA expression levels is comprehensively presented in the heatmap (P ＜ 0.05, Figure 9).

Interestingly, the ROC and survival analyses for the AJCC stage could effectively predict the prognosis of CRC patients. For the integrated analysis of risk score and AJCC stage, the results showed that patients with a high-risk score and tumor grade had a significantly worse prognosis (P ＜ 0.001, Figure 10).

Four-lncRNA model-associated biological pathways

GSEA was carried out between the high- and low-risk groups (253 for each group) in the TCGA cohort to identify the high-risk-associated pathways. We found that“epithelial-mesenchymal transition (EMT),” “MAPK signaling pathway,” and some other cancer-related biological pathways were significantly enriched (P ＜ 0.05, Figure 11). The results suggest that the four-lncRNA model may represent the status of these biological processes to further predict the OS of CRC patients.

Nomogram construction

Based on all the analyses, a comprehensive nomogram that integrated the risk scores and clinical information was constructed to predict the OS of CRC patients. Using the nomogram results, we could better predict the 1-, 3-, and 5-year survival probabilities of CRC patients (P ＜ 0.05, Figure 12).

qRT-PCR validation

Figure 1 Volcano map of differentially expressed RNAs in the Cancer Genome Atlas database.

Based on the results of the bioinformatics analysis, we found that lncRNA LINC00488 involved in the ceRNA network was a novel lncRNA. Therefore, we evaluated the mRNA expression levels of LINC00488 by qRT-PCR in 22 paired patient samples(CRC and normal tissues). The results showed that LINC00488 was downregulated in the CRC tumors compared to the corresponding non-tumor tissues (P ＜ 0.05, Figure 13). These results demonstrated that our bioinformatics analysis was reliable.

DISCUSSION

CRC is the third leading cause of cancer mortality[1,13]. Despite significant advances in the diagnosis and treatment of cancer, progress in improving the OS of CRC patients has been slow. LncRNAs have been the focus of the search for new biomarkers to improve the prognosis of these patients due to the accumulating evidence that lncRNAs have an essential regulatory function in tumorigenesis[14]. LncRNAs can act as ceRNAs in both tumorigenesis and biological processes. Unlike other methods for ceRNA network construction, we performed univariate regression analysis on all the aberrantly expressed lncRNAs and screened those closely related to survival prognosis to subsequently construct a ceRNA network and prognostic model. This method effectively filtered out a large number of lncRNAs that were unrelated to survival, which provides an effective basis for survival prognosis evaluation for CRC patients.

In recent years, several lncRNAs have been identified as diagnostic and prognostic biomarkers[15,16]. Our ceRNA network was constructed from 9 lncRNAs, 30 miRNAs,and 44 mRNAs. We found that LINC00488 was downregulated while the other eight lncRNAs (AL139147.1, DSCR4, DSCR8, H19, KCNQ1OT1, LMO7-AS1, MIR31HG, and MYO16-AS1) were upregulated in CRC tumors compared to non-tumor tissues. Ding et al[17]found that H19 was upregulated in tumors from 185 CRC patients compared to non-tumor tissues and promoted EMT through the H19/miR-29b/PGRN axis in CRC.Wang et al[18]demonstrated that high H19 expression was associated with 5-fluorouracil resistance. Eide et al[19]investigated the relationship between miR-31 and its host transcript lncRNA MIR31HG and found a strong correlation. In particular,their results showed that high lncRNA MIR31HG expression correlated with a shorter relapse-free survival in CRC patients. Numerous studies have demonstrated that lncRNAs DSCR8 and DSCR4 are involved in the progression of multiple cancers[20].Moreover, Li et al[21]demonstrated that lncRNA KCNQ1OT1 was upregulated in colon cancer and regulates oxaliplatin chemoresistance through the miR-34a/ATG4B pathway. Enrico et al[22]also demonstrated that lncRNA KCNQ1OT1 could be a predictive clinical biomarker for CRC patients. All these results are consistent with our data presented in this study.

We found a novel lncRNA (LINC00488) in our ceRNA network. Expression analysis of LINC00488 in 22 paired CRC tumor tissues and adjacent non-tumor tissues confirmed the bioinformatics results. Correlation regression analysis revealed that AL139147.1 had direct linear correlations with forkhead box protein Q1 (FOXQ1) (R =0.383, P ＜ 0.001) and protein chromobox 2 (CBX2) (R = 0.354, P ＜ 0.001). Several studies have demonstrated that many lncRNAs can interact with these two proteins to regulate cancer prognosis[23,24].

Gene expression profile-based prognostic lncRNAs, miRNAs, and mRNA signatures for prognosis prediction in patients with cancer have been investigated[10,25].Seeking more reliable prognostic biomarkers for patient survival seems necessary to provide more personalized cancer management strategies for CRC patients. The contribution of mRNAs to CRC clinical prognosis has been demonstrated in the past decade. Several studies have also suggested the potential utilization of lncRNAs as prognostic markers and therapeutic targets[26-28]. However, the lncRNAs involved in the ceRNA network have not yet been reported for prognosis. Through multivariate Cox regression analysis, we identified four putatively informative genes (LINC00488,LMO7-AS1, MIR31HG, and DSCR4), which were used for the establishment of the four-lncRNA model. Based on this model, patients in the high-risk group survived significantly shorter than those in the low-risk group.

Figure 2 Functional enrichment analysis of Gene Ontology.

The AJCC stage provides important clinical information and is preferred for the management of CRC[29]. Interestingly, the prognostic prediction was more effective by combining the risk score and AJCC stage, which indicates that the four-lncRNA model based on risk score could improve the accuracy of prognostic prediction in CRC. Furthermore, nomogram, a statistical model based on the combination of risk score and clinical features, provides a more extensive prediction for each patient of CRC. Based on the nomogram, we could predict the 1-, 3-, and 5-year survival probabilities.

EMT is a major player in modulating tumor metastasis[30]. Many signaling pathways contribute to the regulation of EMT, such as Wnt/β-catenin and Notch[31,32]. The GSEA results showed that the “EMT” was significantly enriched. In this study, we found that COL1A1 and TPM2, which are key genes in EMT, were significantly expressed in the ceRNA network. COL1A1 encodes alpha 1 type I collagen, which is a target protein for TGFβ1. Pudova et al[33]found that COL1A1 is positively correlated with hexokinase 3 (HK3), and HK3 upregulation is associated with EMT. The study showed that TPM2 K152Ac was the most significantly downregulated protein in liver metastases in CRC[34]. We also found other signaling pathways (e.g., MAPK signaling pathway[35], regulation of binding[36], and focal adhesion[37]) that correlated with cancer.These results revealed a more extensive perspective in understanding the signaling pathways involved in CRC.

The major findings reported in this study could be considered a preliminary attempt to explore the prognostic lncRNAs, miRNAs, and mRNAs in CRC. However,this study has some limitations: All CRC patient data used to construct the prognostic signature were retrieved from the TCGA database. We only used one dataset in GEO to validate the predictive value of our model. Further validation based on other independent datasets is required to evaluate the general applicability of the signature in CRC. In addition, only LINC00488 related to OS was investigated for its expression pattern using qRT-PCR in 22 paired CRC patient samples. Additional studies are needed to confirm the biological functions of the lncRNAs, miRNAs, and mRNAs of the ceRNA network in CRC.

Table 2 Clinical and pathological characteristics of 506 colorectal cancer patients

Figure 4 Transcription factor-differentially expressed mRNA regulatory network.

Figure 5 mRNA-long noncoding RNA-miRNA regulatory network.

Figure 6 Linear regression of differentially expressed RNA expression levels (R ＞ 0.35, P ＜ 0.01).

Figure 7 Predictive value assessment of the four-long noncoding RNA model.

Figure 8 Cox regression analysis of clinical features and risk score in colorectal cancer patients.

Figure 9 Heatmap of the four-long noncoding RNA model combined with clinical features and long noncoding RNA expression.aP ＜ 0.05.

Figure 10 Prognostic value of clinical features and risk score for the overall survival of colorectal cancer patients (P ＜ 0.01).

Figure 11 Gene Set Enrichment Analysis for high-risk patients (P ＜ 0.05).

Figure 12 Nomogram construction.

Figure 13 Expression levels of two long noncoding RNAs in the competing endogenous RNA network related to overall survival.

ARTICLE HIGHLIGHTS