Small Nucleolar RNAs (snoRNAs)-Based Risk Score Classifier Predicts Overall Survival in Bladder Carcinoma

Background

Bladder cancer is a heterogeneous malignancy that is responsible for an estimated 81 400 new cases and 17 980 deaths in the United states in 2020 [1]. Although progress has been made in surgical resection technology and neoadjuvant chemotherapy, survival after surgical resection differs considerably among BLCA patients [2–4]. The predictive accuracy of the current clinical staging system is insufficient for this medically precise era and for the display of the molecular characteristics of BLCA. Hence, cancer researchers have been interested in developing an accurate stratification method that can predict the prognosis and direct treatment of BLCA patients to improve their survival [5–10]. It is also imperative to determine the biological characteristics of BLCA.

Genome researchers have long regarded the non-coding RNAs of the human genome as ‘junk’ DNA [11]. Specifically, the rapid development of high-throughput technology helped scientists to identify a large number of non-coding RNAs, which account for almost 60% of the transcriptional output in human cells [12–16]. This huge number of non-coding RNAs makes it difficult to dismiss them as “junk”. Hence, recent extensive research has identified that non-coding RNAs, especially long non-coding RNAs and micro RNAs, are directly linked to the tumorigenesis and development of cancers [17,18]. Notably, small nucleolar RNA (snoRNAs), a class of small (60–300 nucleotides) non-coding RNAs, has been well documented in rRNA biogenesis. Recently, however, some studies have highlighted that snoRNAs are involved in cancer development and progression [19]. Because these studies dismissed the non-coding RNAs of the human genome as ‘junk’ DNA, the relationships between snoRNAs and BLCA have not been well characterized.

The Cancer Genome Atlas (TCGA) consortium and other public genomic datasets offer high-throughput data for the subclassification of BLCA, as well as predicting diagnosis and prognosis [20–22]. Gong et al. developed a database of snoRNA in cancers to systematically quantify and deposit snoRNA expression profiles in more than 10 000 samples across 31 cancer types based on TCGA, including BLCA [23]. This work greatly promoted the development of snoRNAs analysis in cancers. However, very few other studies have explored the role of snoRNAs in BLCA. Furthermore, no investigative reports on models that predict the survival status of BLCA based on snoRNAs are available.

Based on TCGA data portal, we analyzed the expression profiles and clinical outcomes of BLCA patients. Survival-associated snoRNAs were identified and submitted to develop a novel 12-snoRNAs-based risk score classifier. To leverage the complementary value of molecular and clinical parameters, clinical factors were integrated to build a nomogram, which allowed improved prediction of BLCA patient survival.

Material and Methods

DATA ACQUISITION:

Expression profiles of snoRNAs were downloaded from the SNORic database [23]. Corresponding clinical data of BLCA patients were acquired from the TCGA database (https://cancergenome.nih.gov/).

SURVIVAL ANALYSIS:

To generate a prognostic classifier to predict the outcome of BLCA, we used univariate Cox analysis to explore the relationships between snoRNAs expression levels and the overall survival (OS) of BLCA patients. Patients with OS less than 90 days were removed to obtain more accurate results [24]. Then, multivariate Cox regression analysis was conducted to screen out independent prognostic factors. We calculated a signature involving the single prognostic parameters achieved from the previous step to assess the OS risk based on the individual expression of the prognostic snoRNAs, weighted by the regression coefficient. The data were then divided into high- or low-risk groups by using the median risk score as the threshold value. To leverage the complementary value of the snoRNAs and several indispensable clinical parameters, including age, sex, tumor T stage, tumor N stage, tumor M stage, and histological grade, we integrated them by using a nomogram graph.

FUNCTIONAL CHARACTERIZATION OF THE SNORNAS-BASED RISK SCORE:

To further gain a biological understanding of the snoRNAs-based risk score, we collected genes in the whole genome which were significantly correlated with the included prognostic snoRNAs that were obtained from SNORic. Significantly correlated genes were identified with the |Spearman coefficient| ≥0.4 and FDR <0.05. These significantly correlated genes were submitted to a “clusterProfiler” package in R software for gene functional enrichment analysis. The background of the gene list was measured by “Homo sapiens.” Statistically significant terms were identified when the p-value and the q-value were both less than 0.05.

STATISTICAL ANALYSIS:

All statistical analyses were performed using SPSS (version 22.0, SPSS Inc., Chicago, IL, USA) and R (version 3.3.1; https://www.r-project.org/). snoRNAs expression values were converted to log2 (RMPK+1) and the average expression value across all samples less than 1 were identified as undetected snoRNAs. The difference for snoRNAs with OS was evaluated using the log-rank test. Univariate survival Cox analysis was performed by using “survival” packages of R software. Time-dependent (receiver operating characteristic) ROC curve and corresponding area under the curve (AUC) values, which could assess the performance of prognostic signatures, were calculated by using the “survivalROC” package. The nomogram graph was derived by using the “RMS” package. Statistical significance was defined as P<0.05, unless specified otherwise.

Results

IDENTIFICATION OF SURVIVAL-ASSOCIATED SNORNAS:

Expression profiles of 412 snoRNAs were acquired from the SNORic database. Finally, a total of 366 BLCA patients with OS more than 90 days were included in the present study. Table 1 shows the clinical information of the patients with BLCA from TCGA. Univariate Cox hazard analyses were performed to assess the relationship between clinical parameters and clinical outcomes of the BLCA patients. Age (HR=1.792, 95% CI: 1.262–2.543; P=0.001) and AJCC stage (HR=2.161, 95% CI: 1.459–3.202; P<0.001) were significantly correlated with OS of BLCA patients. However, no significant correlations were observed between OS and sex or histological grade. All 368 snoRNAs were submitted to a univariate Cox analysis and 58 snoRNAs were identified as survival-associated snoRNAs (Figure 1, Table 2).

CONSTRUCTION OF SNORNAS-BASED PROGNOSTIC SIGNATURE:

Significant factors from univariate selection were kept in the multivariate analysis by using backward selection (Table 3). A multivariate Cox regression analysis was used to build a prognostic signature that selected 12 out of the 58 snoRNAs. Kaplan-Meier (K-M) curves were used to display the prognostic value of each snoRNA (Figure 2). We then derived a prognostic signature for each patient based on the individual expression levels of the 12 survival-associated snoRNAs multiplied by their coefficients in the multivariate Cox analysis: prognostic signature=(SNORD114–11 * (−0.168)+SNORD114–14 * (0.201)–SNORD114–15 * (0.229)+SNORD114–9 * (0.543)–SNORA55 * (0.198)–SNORA60 * (0.192)–SNORD88A * (0.167)–SNORD69 * (0.314)–SNORD20 * (0.226)+U49A * (0.432)–SNORD51 * (0.335)+U74 * (0.514) (Figure 3). Patients were separated into high- and low-risk groups by the median value of the prognostic signature. K-M survival plots indicated that patients in the high- and low-risk groups had distinct clinical outcomes (HR=2.500, 95% CI: 1.828–3.420, P<0.001; Figure 4A). The AUC value of the ROC curve was 0.719. The threshold was 2000 days (Figure 4B). The AJCC TNM stage appears to be a moderate prognostic indicator for clinical outcome predicting (HR=2.155, 95% CI: 1.553–2.993, P<0.001; Figure 4C). The AUC value of the TNM stage was 0.636 (Figure 4D). We integrated the clinical factors and risk score to create a composite nomogram based on the results of the multivariate Cox regression analyses to predict 1-year, 2-year, and 3-year OS. (Figure 5). We then submitted clinical parameters and risk scores in multivariate Cox analysis. Notably, the prognostic signature remained an independent significant prognostic risk factor (HR=3.300, 95% CI: 2.203–4.943, P<0.001) (Table 4).

MOLECULAR FUNCTION OF PROGNOSTIC SIGNATURE:

Genes that were significantly correlated with 12 survival-associated snoRNAs were obtained from SNORic based on Spearman correlation analysis. The network was conducted to display the relationships between snoRNAs and genes (Figure 6). Red lines indicate positive correlation relationships while blue lines indicate negative correlation relationships. SnoRNAs-related genes were further submitted to gene functional enrichment analysis. For the biological process, the extracellular structure organization, the extracellular matrix organization, and skeletal system development were the commonly enriched categories (Figure 7A). For the cellular component ontology, the enriched categories were correlated with proteinaceous extracellular matrix, endoplasmic reticulum lumen, and contractile fiber (Figure 7B). With regards to the molecular function, the snoRNA-related genes mainly showed enrichment in cell adhesion molecule binding, actin binding, and sulfur compound binding (Figure 7C). The disease ontology suggested that these genes were enriched in several types of cancer (Figure 8A). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed significant pathways with these genes (Figure 8B). “Focal adhesion” was the most significant of the enriched terms.

Discussion

Here, we performed a systematic analysis of the snoRNAs and identified a risk score based on the expression profiles of 12 survival-associated snoRNAs in BLCA patients based on 366 clinical cases. Our study resulted in the following: (1) the identification of 58 prognostic relevant snoRNAs; (2) the development of a twelve-snoRNAs-based risk score classifier that predicts OS in BLCA; (3) a pathway analyses that revealed the molecular characteristics of the risk score; and (4) the construction of a nomogram to leverage the complementary value of molecular and clinical factors.

We first identified several snoRNAs that were correlated with the clinical outcomes of BLCA. Given the stable nature of snoRNAs in the human body, these snoRNAs have inherent advantages for use as molecular biomarkers [25–27]. Regrettably, few studies have assessed snoRNAs as diagnostic and prognostic tools for BLCA. This is mainly due to the conventional prejudice that snoRNAs mainly function to modify, mature, and stabilize rRNAs [28]. Furthermore, high-throughput DNA sequencing techniques help to identify cancer-specific snoRNAs [29]. Hence, our group comprehensively analyzed the clinical significance and potential molecular characteristics of snoRNAs to identify the diagnostic and prognostic biomarkers in BLCA.

An important goal of the present study was the construction of risk scores based on the expression pattern of snoRNAs to create a risk stratification model for the practice of precision medicine. Previously, studies involved several molecules, including mRNA expression, copy number variation, DNA methylation, lncRNAs, and miRNAs. For example, Liu et al. proposed a clinical multidimensional transcriptome signature for survival predictions of patients with BLCA [30]. Aberrant DNA methylation can provide reliable biomarkers in the prediction of clinical outcomes of common urological cancers [31]. Several studies have explored the prognostic value of non-coding RNAs in developing prognostic signatures for BLCA [32,33]. We constructed the risk score model with the expectation that it would be applicable to clinical management from the perspectives of snoRNAs. On the one hand, studies that refer to the clinical significance of snoRNAs are limited. The present study provides novel insights into the clinical value and molecular mechanisms of snoRNAs in BLCA. On the other hand, we optimized the prognostic model and considered both molecular and clinical features. The snoRNAs possess the possibility and feasibility as biomarkers and therapeutic targets in clinical use. Aberrant snoRNAs expression was found in many cancers, and the expression level was correlated with diagnosis, classification of subtypes, and patient survival. Moreover, snoRNAs were stably expressed and detectable in body fluids, including blood plasma, serum, and urine of cancer patients, indicating that snoRNAs have the potential to serve as biomarkers in clinical use. Urine testing is very important for early diagnosis of BCLA. Since snoRNAs are stably expressed and detectable in urine, snoRNAs detection in urine might be useful for BLCA diagnosis and prognosis prediction [34–36].

It is interesting to explore the biological characteristics reflected by the risk score. Nowadays, the cellular regulatory roles of snoRNAs in cancers are widely understood. Hence, we selected snoRNAs-related genes to explore the molecular characteristics of snoRNAs in BLCA. We found that prognostic snoRNAs are mainly involved in several signal transduction pathways, such as focal adhesion, extracellular matrix (ECM)-receptor interaction, and peroxisome proliferator-activated receptor (PPAR) signaling pathway. In these findings, the risk score could reflect the cell-cell interaction status of BLCA. Notably, focal adhesion was the most significant KEGG pathway in the pathway functional enrichment analysis. In biological activity, focal adhesion is a sub-cellular regulatory structure that mediates mechanical force and regulatory signals transmitted between the ECM and an interacting cell [37]. The focal adhesion pathway interacts closely with other indispensable oncogenic pathways and is actively involved in the progression of cancers [38,39]. Interestingly, inhibitors focused on the focal adhesion pathway could be effective anti-tumor targets [40]. The biological characteristics of the risk score are correlated with cell-cell interactions, which indicated that these snoRNAs actively participate in the progression and/or metastasis of BLCA, although the results need to be further explored.

Conclusions

In summary, the development of a prognostic signature, as defined by expression profiles of 12 survival-associated snoRNAs, could be an excellent predictor of the clinical outcome of BLCA patients. However, studies of other independent cohorts are required to validate these findings. There is no evidence that the prognostic signature can predict prognosis in patients who received adjuvant therapies after surgical resection. In the present study, the process involves mechanisms that should be validated by in vitro or in vivo experiments. Clinical information was integrated into snoRNAs expression profiles for the first time to construct a snoRNAs-based risk score by our group.