A novel prognostic nomogram based on 5 long non-coding RNAs in clear cell renal cell carcinoma.

Clear cell renal cell carcinoma (ccRCC) is the most common and invasive histological subtype of all kidney malignancies, with high levels of incidence and mortality. In the present study, long non-coding (lnc)RNA expression profiles of patients with ccRCC from The Cancer Genome Atlas database were comprehensively analyzed to identify differentially expressed lncRNAs (DElncRNAs). The patients with ccRCC were then divided into training and validation cohorts. Univariate and LASSO regression analyses were performed to select the most significant survival-associated candidate DElncRNAs in the training cohort. Multivariate Cox regression analysis was then performed to develop a risk score formula and a prognostic nomogram for predicting 3- and 5-year overall survival (OS). The accuracies of the nomogram predictions were evaluated by determining the area under the receiver operating characteristic curve (AUC) and a calibration plot. Finally, functional enrichment analysis and protein-protein interaction network prediction were implemented to predict the functions and molecular mechanisms of the candidate DElncRNAs in ccRCC. A total of 1,553 DElncRNAs were identified, and 5 candidate DElncRNAs (AC026992.2, AC245041.2, LINC00524, LINC01956 and LINC02080) were included in the nomogram. The AUC values for 3- and 5-year overall survival in the training cohort were 0.768 and 0.814, respectively, which were increased compared with that based on the clinical index (0.760 and 0.694, respectively). Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analyses revealed that the 521 mRNAs highly associated with 5 DElncRNAs were primarily involved in 17 terms and 25 pathways, respectively. Based on the 5 DElncRNAs, a novel and convenient prognostic nomogram for predicting 3- and 5-year OS for patients with ccRCC was developed. The results of the present study may be conducive to the development of a precise predictive tool for the prognosis of ccRCC and may provide information regarding the molecular mechanisms of ccRCC. However, additional experimental in vitro and in vivo studies investigating lncRNAs may be required. Copyright: © Wang et al.

Introduction

Renal cell carcinoma (RCC) is the most prevalent type of malignant tumor of the adult kidney worldwide (1). In the United States of America, an estimated ~64,000 individuals were diagnosed with RCC in 2017 (2). Clear cell (cc)RCC is the most common and invasive histological subtype of RCC and accounts for 70–80% of all RCC cases (3,4). The 5-year survival of patients with ccRCC diagnosed in the early stage is >90% (5). However, for patients diagnosed at the advanced stage, the 5-year survival is as low as 12% (5). Furthermore, a large proportion of patients have been diagnosed at the advanced stage (6). The high mortality rate of patients with ccRCC in the advanced stage may be due to lack of typical symptoms and biomarkers with high sensitivity and accuracy for diagnosis in the early stage, and an absence of a reliable risk stratification method for assessing prognosis. Therefore, there is an urgent requirement to identify tumor-specific biomarkers and to develop a nomogram for the precise prediction of prognosis, which may lead to the development of an accurate risk stratification method and guide the clinical diagnosis and treatment of ccRCC.

Long non-coding RNAs (lncRNAs), a class of transcripts, are >200 nt long and lack any protein-coding capacity (7). Previous studies have indicated that lncRNAs are aberrantly expressed in diverse malignant tumor types and serve as pivotal regulators of different biological processes in tumors, including cell proliferation, invasion, apoptosis and metastasis (8). For instance, Ning et al (9) demonstrated that lncRNA nuclear paraspeckle assembly transcript 1 (NEAT1) was overexpressed in ccRCC tissue, which was significantly associated with poor prognosis in patients with ccRCC and silencing of NEAT1 was able to inhibit ccRCC cell proliferation and invasion. In addition, Yang et al (10) demonstrated that lncRNA PVT1 positively regulated ccRCC cell proliferation and invasion by interacting with microRNA (miR)-200s through increasing the expression of zinc finger E-box binding homeobox 1 (ZEB1), ZEB2 and polycomb complex protein BMI-1.

Advances in omics technology have created opportunities to mechanistically elucidate, diagnose and treat cancer in a systematic manner (11). RNA-sequencing profiling has been developed, which may be utilized for identifying novel molecular markers and mechanisms in numerous tumor types (12). Furthermore, with their advantage of the combination of independent prognostic factors to consider their multiple effects on probability of outcome, prognostic nomograms have been widely used as a powerful model for risk evaluation in cancer (13).

In the present study, all original matrix files for ccRCC and corresponding clinical information were downloaded from The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov) to detect the differentially expressed (DE)lncRNAs between ccRCC tissues and normal renal tissues. The best survival-associated candidate DElncRNAs were then selected to develop a prognostic nomogram. Furthermore, functional enrichment analysis was performed to predict the biological functions of the candidate DElncRNAs. The present study may contribute to the determination of independent prognostic lncRNAs in ccRCC and provide additional information regarding the molecular mechanisms of ccRCC.

Materials and methods

Data retrieval and processing

The lncRNA expression profiles were obtained from TCGA, an open database to identify novel biomarkers in cancer research (version: April 5, 2018), and were then subjected to background correction and normalization with Perl 5.0 (http://www.perl.org/). Patients with a follow-up time of <30 days or lack of pathologic diagnosis and corresponding clinical information were removed. The data of 574 tissue samples were included in the present study, comprising 70 adjacent non-tumor kidney tissues and 504 ccRCC tissue samples. Relevant clinical characteristics of the 504 cancer cases were also obtained and collated.

Patient cohort

The 504 patients with ccRCC were randomly divided into two cohorts: The training cohort and the validation cohort. The training cohort comprised 380 cancer cases and the remaining cases were in the validation cohort. For categorical variables, data were expressed as numbers and compared using χ2 tests or Fisher's test, whereas for continuous variables, data were expressed as the mean ± standard deviation and compared using Student's t-tests in SPSS 20.0 (IBM Corp.). P<0.05 was considered to indicate a statistically significant difference.

Processing of lncRNA expression data

The statistical software R (version 3.5.2; http://www.R-project.org) and the Bioconductor package ‘edgeR’ (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) (14) were used to identify DElncRNAs with the criterion of |log fold change (FC)|>2 and adjusted P<0.01, as described previously (15).

Construction and evaluation of DElncRNA-based prognostic nomogram

In the training cohort, a univariate regression analysis was performed to select DElncRNAs that were highly associated with the overall survival (OS) of patients with ccRCC. Subsequently, a LASSO regression analysis was performed to additionally screen out the set of independent prognostic candidate DElncRNAs with the strongest predictive power. Next, the most significant survival-associated candidate DElncRNAs were subjected to a multivariate Cox regression analysis to develop a risk score formula. On the basis of this risk score formula, the risk score of each patient was calculated and patients were then stratified into low- and high-risk groups (cut-off =0.8863). Ultimately, a prognostic nomogram of candidate DElncRNAs was developed for predicting 3- or 5-year OS probabilities. Furthermore, a calibration plot with a bootstrapping set of 1,000 resamples and the receiver operating characteristic curve (ROC) were generated to appraise the predictive capacity of the prognostic nomogram by calculating the area under the curve (AUC). In the same manner, a ROC curve in the validation cohort and a ROC curve based on clinical information were also generated to validate the prognostic nomogram.

Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses

To date, the functions of the majority of lncRNAs remain unknown. Therefore, to investigate the roles of candidate DElncRNAs, Pearson correlation coefficients between candidate DElncRNAs and mRNAs in the expression matrix were calculated using R software. mRNAs with a Pearson correlation coefficient of >0.45 and P<0.01 were selected and subjected to KEGG and GO analysis with an adjusted P<0.01 set as the threshold.

Protein-protein interaction (PPI) network construction and module selection

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (version 10.5; http://string-db.org) (16), comprising 9643,763 proteins from 2,031 organisms, was adopted to predict PPIs for ‘Homo sapiens’, with a confidence score of >0.9 set as the cut-off criterion. Data from the PPI network were processed by Cytoscape and the hub genes (top 10 degree genes in the PPI network), which were considered to be the most important targets of candidate DElncRNAs, were selected. Subsequently, an online tool, UALCAN (http://ualcan.path.uab.edu/index.html) (17), was used to assess the expression of the hub mRNAs in ccRCC.

Survival analysis

A Kaplan-Meier plot and the log-rank test were used to construct survival curves and assess significant differences in OS between the low- and high-risk groups and the associations between the expression levels of candidate DElncRNAs and OS in patients with ccRCC. P<0.05 was set as a cutoff value.

Results

Patient characteristics

A total of 504 patients with ccRCC (380 cases in the training cohort and 124 cases in the validation cohort) were included in the present study. All patients had been pathologically diagnosed with ccRCC. The detailed demographic and baseline characteristics of the two cohorts are summarized in Table I. The mean age was 60.46±12.05 and 61.3±12.77 years in the training and validation cohorts, respectively. There was no significant difference in any of the clinicopathological parameters, including age, sex, ethnicity, tumor stage and survival status between the two cohorts.

Table I.

Baseline characteristics of all patients with clear cell renal cell carcinoma.

Demographic characteristics	Training cohort (n=380)	Validation cohort (n=124)	Total (n=504)	P-value
Age, years	60.46±12.05	62.12±12.76	60.87±12.08	0.272
Sex				0.914
Male	247	133	329
Female	82	42	175
Ethnicity				0.668
Caucasian	338	108	446
African descent	37	13	50
Asian	5	3	8
TNM stage				0.264
I	180	68	254
II	40	14	54
III	94	20	114
IV	60	22	82
Survival status				0.661
Alive	251	85	336
Dead	129	39	168

Identification of DElncRNAs in ccRCC

With the selection criteria set as |logFC|>2 and P<0.01, 1,064 upregulated and 489 downregulated DElncRNAs were identified between ccRCC tissues and normal renal tissues (Fig. 1A). The heatmap of the top 30 DElncRNAs is provided in Fig. 1B.

Figure 1.

Significant DElncRNAs identified in the present study. (A) Volcano plot of DElncRNAs. Red represents 1,064 upregulated lncRNAs and green represents 553 downregulated lncRNAs with |log (FC)|>2 and P<0.01. (B) Heat maps of the top 30 DElncRNAs. Red represents upregulated lncRNAs and green represents downregulated lncRNAs. DElncRNA, differentially expressed long non-coding RNA; FC, fold change; FDR, false discovery rate.

By using univariate regression analysis, 135 DElncRNAs significantly associated with OS of patients were screened out and then additionally analyzed by LASSO regression analysis in the training cohort (Fig. 2A and B). A total of 5 candidate DElncRNAs (AC026992.2, AC245041.2, LINC00524, LINC01956 and LINC02080) were selected to develop a prognostic nomogram (Fig. 2D) and a risk score formula using multivariate Cox regression analysis (Fig. 2C) as follows: Risk score=(expression level of AC026992.2 ×-0.13804) + (expression level of AC245041.2 ×0.10658) + (expression level of LINC00524 ×0.10723) + (expression level of LINC01956 ×0.09903) + (expression level of LINC02080 ×0.31598). The AUC of the prognostic nomogram for 3- and 5-year OS in the training cohort was 0.768 and 0.814, respectively, with a Harrell's concordance index (C-index) of 0.729 (Fig. 3A). The calibration plots for 3- or 5-year survival probabilities in the training cohort are presented in Fig. 3D and E, respectively. The distribution of risk score, survival status and expression profile of the 5 prognostic DElncRNAs for the training cohort are depicted in Fig. 4A.

Figure 2.

Prognostic nomogram based on 5 candidate differentially expressed lncRNAs. (A) LASSO regression coefficient profiles of survival-associated lncRNAs. (B) ‘Leave-one-out’ cross-validation for parameter selection during LASSO regression. (C) Forest plot of multivariate Cox regression analysis. (D) Prognostic nomogram for the prediction of 3- and 5-year overall survival of patients with clear cell renal cell carcinoma. lncRNA, long non-coding RNA.

Figure 3.

Evaluation of nomogram. (A) Area under the receiver operating characteristic curve for 3- and 5-year OS of patients with clear cell renal cell carcinoma in the training cohort based on 7 candidate DElncRNAs. (B) Area under the receiver operating characteristic curve for 3- and 5-year OS of patients with clear cell renal cell carcinoma in the validation cohort based on 7 candidate DElncRNAs. (C) 3- and 5-year survival based on the clinical index. (D) Calibration curve of the nomogram model for 3-year OS in the training cohort. (E) Calibration curve of the nomogram model for 5-year OS in the training cohort. DElncRNA, differentially expressed long non-coding RNA; OS, overall survival.

Figure 4.

Distribution of risk score, survival status and expression profile of 7 prognostic differentially expressed long non-coding RNAs in the (A) training cohort and in the (B) validation cohort.

Verification of the candidate DElncRNA nomogram

To confirm the accuracy of the prediction of 3- or 5-year OS by the nomogram established in the training cohort, it was validated in the 124 patients with ccRCC of the validation cohort. Similar to the procedure in the training cohort, the 124 patients in the validation cohort were stratified into high- and low-risk groups according to their risk score. The AUC for predicting 3- and 5-year OS in the validation cohort was 0.899 and 0.869, respectively, with a C-index of 0.88 (Fig. 3B). The AUC for predicting 3- and 5-year OS based on the TNM staging (AJCC 7th edition, 2010) (18) was 0.760 and 0.694, respectively (Fig. 3C). The distribution of the risk score, survival status and expression profile of the 5 candidate DElncRNAs in the validation cohort are demonstrated in Fig. 4B.

Functional enrichment analysis

A total of 521 mRNAs, the expression levels of which were associated with the 5 candidate DElncRNAs (Pearson correlation coefficient >0.45 and P<0.01), were identified. The GO analysis revealed that the 521 mRNAs were enriched in 17 terms, including receptor ligand activity and channel activity (Fig. 5A). The KEGG analysis indicated that the 521 mRNAs were primarily involved in 25 pathways, among which a few pathways were highly associated with oncogenesis, progression and metastasis of neoplasm, including the Wnt, p53 and mTOR signaling pathways (Fig. 5B).

Figure 5.

Prediction of functions of candidate differentially expressed long non-coding RNAs. (A) Gene Ontology analysis. (B) Kyoto Encyclopedia of Genes and Genomes analysis. P.adjust, adjusted P-value.

PPI network construction and module selection

A PPI network was constructed using the online database STRING and then processed with Cytoscape software (Fig. 6A). 5 hub mRNAs were selected as the most crucial targets of the candidate DElncRNAs, comprising cyclin D kinase 1 (CDK1; degree of connectivity =27), aurora kinase B (AURKB; degree of connectivity =21), cyclin B1 (CCNB2; degree of connectivity =20), BUB1 mitotic checkpoint serine/threonine kinase (BUB1; degree of connectivity =20) and ubiquitin conjugating enzyme E2 C (UBE2C; degree of connectivity =20) (Fig. 6B). The result obtained with the UALCAN tool demonstrated that the 5 mRNAs were expressed at increased levels in ccRCC compared with that in normal tissues (Fig. 6C).

Figure 6.

PPI network and hub genes. (A) PPI network. (B) Five hub mRNAs in the PPI network. (C) Expression of top five degree mRNAs in clear cell renal cell carcinoma. miRNA, microRNA. KIRC, clear cell renal cell carcinoma. PPI, protein-protein interaction; CDK1, cyclin D kinase 1; AURKB, aurora kinase B; CCNB2, cyclin B1; BUB1, BUB1 mitotic checkpoint serine/threonine kinase; UBE2C, ubiquitin conjugating enzyme E2 C.

Kaplan-Meier analysis was performed to determine the association between OS of patients with ccRCC and risk score, and the expression levels of 5 candidate DElncRNAs. The result of the survival analysis demonstrated that the prognosis of the patients with ccRCC in the high-risk group was poorer compared with that in the low-risk group (Fig. 7A). Furthermore, high expression of AC245041.2 (P<0.001), LINC00524 (P<0.001), LINC01956 (P<0.001) and LINC02080 (P<0.001), along with low expression level of AC026992.2 (P<0.001), was associated with a poor outcome for patients with ccRCC (Fig. 7B-F).

Figure 7.

Prognostic value of the selected long non-coding RNAs regarding overall survival. (A) Low- and high-risk group. (B) AC026992.2 (P<0.001). (C) AC245041.2 (P<0.001). (D) LINC00524 (P<0.001). (E) LINC01956 (P<0.001). (F) LINC02080 (P<0.001).

Discussion

ccRCC ranks first in incidence amongst all histological types of kidney malignancies (19), and accounts for nearly 3% of all types of adult malignancies (20). Although the 5-year survival rate of patients with ccRCC at the early stage is >90%, it decreases to 10% in patients with advanced ccRCC (5), and >100,000 patients succumb to ccRCC per year worldwide (21). Therefore, it is imperative to identify tumor-specific markers for risk stratification that may be utilized for assessing the prognosis of patients, and which may facilitate the development of novel strategies for the diagnosis and therapy of ccRCC. lncRNAs, a class of non-coding RNAs of >200 nt in length, have been demonstrated to have a significant role in transcriptional and post-transcriptional regulation, and deregulation of certain lncRNAs is involved in the initiation and progression of various cancer types (22). The roles of lncRNAs have become active areas of research, which will undoubtedly be propitious for the elucidation of the functions of lncRNAs in cancer. However, to date, the functions of the majority of lncRNAs remain elusive (23).

In the present study, a comprehensive analysis of lncRNA expression matrix files and corresponding clinical information of patients with ccRCC was performed. The patients were randomly assigned to training or validation cohorts. In the training cohort, a total of 1,553 DElncRNAs (1,064 upregulated and 489 downregulated DElncRNAs) were identified. Following screening by univariate regression analysis and LASSO regression analysis, the top 5 candidate DElncRNAs (AC026992.2, AC245041.2, LINC00524, LINC01956 and LINC02080) associated with survival were selected and used to develop a prognostic nomogram. In the training cohort, the AUC of the prognostic nomogram for 3- and 5-year OS was 0.768 and 0.814, respectively, and in the validation cohort, it was 0.899 and 0.869, respectively, demonstrating an excellent predictive accuracy for the probability of survival at 3 or 5 years. The AUC for predicting 3-year OS based on the clinical index (0.760) was similar to that in the training cohort. However, the AUC for predicting 5-year OS based on the clinical index (0.694) was markedly decreased compared with that in training cohort, indicating that the accuracy of the prognostic nomogram was improved compared with that of the clinical index for predicting 5-year OS. Previously, Shi et al (24) used 5 lncRNAs (ENSG00000229178, ENSG00000236453, ENSG00000245060, ENSG00000258789 and ENSG00000272558) for predicting 3-year OS of patients with ccRCC. In addition, Qu et al (25) also built a prognostic lncRNA signature for predicting 5-year OS in localized ccRCC with 4 lncRNAs (ENSG00000255774, ENSG00000248323, ENSG00000260911 and ENSG00000231666). In the present study, 5 lncRNAs, which were completely different from lncRNAs used in the studies of Shi et al (24) and Qu et al (25), were used to develop a prognostic nomogram. The AUC for 3-year OS in the training cohort and validation cohort, 0.768 and 0.899, respectively, was superior to that of Shi et al (24), 0.703 and 0.630, respectively. Furthermore, in the study by Qu et al (25), the AUC for 5-year OS in the training and validation cohort was 0.690 and 0.663, respectively, which was lower compared with those of the present study. In addition, the calibration plot for 3- or 5-year OS demonstrated an outstanding consistency between the prediction by the prognostic nomogram and the actual outcome. All of the results suggested that the prognostic nomogram established is suitable for estimating the probability of OS of patients with ccRCC at 3 and 5 years.

Compared with normal renal tissues, the 5 candidate DElncRNAs were aberrantly expressed in ccRCC tissues. Survival analysis for low- and high-risk groups indicated that the high-risk group exhibited a poorer prognosis compared with the low-risk group. Furthermore, high expression levels of AC245041.2, LINC00524, LINC01956 and LINC02080 (P<0.001 for each) along with low expression levels of AC026992.2 (P<0.001) were highly negatively associated with OS of patients with ccRCC.

However, at present, little is known regarding the biological functions of the 5 candidate DElncRNAs. Therefore, to additionally explore the biological roles of the 5 candidate DElncRNAs, 521 mRNAs, the expression levels of which were highly associated with the expression of the 5 candidate DElncRNAs, were selected to perform functional enrichment analysis and a PPI network was constructed. According to the GO analysis, the mRNAs were primarily involved in 17 terms, including ‘receptor ligand activity’ and ‘channel activity’, and the KEGG analysis revealed enrichment in 25 pathways. Several of these pathways are known to be associated with oncogenesis, progression and metastasis of cancer. For example, the Wnt signaling pathway is generally involved in cell proliferation and division via controlling the β-catenin degradation complex (26). Through activation of the Wnt signaling pathway, the lncRNA colon cancer-associated transcript 2 improved the proliferation and invasion of ccRCC cells (27). In addition, the mTOR signaling pathway serves vital roles in modulating diverse biological behaviors, including cell growth, metabolism, protein synthesis and autophagy (28). Liu et al (29) identified that inactivation of the mTOR signaling pathway inhibited apoptosis and promoted cell proliferation in ccRCC. The results of the GO and KEGG analyses conducted in the present study predicted that the 5 candidate DElncRNAs had an important effect on the oncogenesis and progression of ccRCC by affecting a series of biological pathways and processes. The PPI network was constructed to determine the interaction among 521 mRNAs, and 5 hub mRNAs (CDK1, AURKB, CCNB, BUB1 and UBE2C) were identified as the most important targets. All of the 5 top degree mRNAs were overexpressed in ccRCC. CDK1 is essential for the eukaryotic cell cycle by regulating the onset of mitosis and the centrosome cycle (30). A previous study indicated that CDK1 expression was highly associated with the prognosis of patients with RCC (31). Furthermore, Li et al (32) indicated that through targeting CDK1, miR-31-5p inhibited the proliferation, migration and invasion of ccRCC cells. At present, little is known about the functions of AURKB, CCNB, BUB1 and UBE2C in ccRCC. Therefore, future studies investigating the functions of AURKB, CCNB, BUB1 and UBE2C in ccRCC re required.

Although the prognostic nomogram established in the present study demonstrated good predictive accuracy for patients with ccRCC, there are a few limitations that should be addressed. Firstly, as all cases were retrieved from TCGA database, the risk of selection bias could not be excluded. Furthermore, to date, no experimental studies have been performed to examine the functions of 4 of the identified lncRNAs in cancer. Therefore, further in vitro and in vivo studies are required to confirm the results of the present study.

In conclusion, 5 candidate DElncRNAs (AC026992.2, AC245041.2, LINC00524, LINC01956 and LINC02080) were identified in the present study, which were independent prognostic factors for patients with ccRCC, and exhibited potential utility as powerful molecular biomarkers for prognosis and risk assessment. A novel and convenient prognostic nomogram was then developed for predicting 3- and 5-year OS for patients with ccRCC based on lncRNAs. The results of the present study may contribute to an improved understanding of ccRCC at the molecular level. However, additional experimental research concerning lncRNAs in vitro and in vivo is required to verify the results of the present study.