Construction and validation of an RNA-binding protein-associated prognostic model for colorectal cancer

Data sources

mRNA sequencing and clinical data for CRC were downloaded from the TCGA database (Data Release 24.0-7 May 2020, https://portal.gdc.cancer.gov/). Detailed data selection criteria described in our previous study: filter criteria of mRNA sequencing: (1) the keywords of cases are (Primary Site) “colon, rectosigmoid junction and rectum”, (Program) “TCGA”, (Project) “TCGA-COAD, TCGA-READ”, (Disease Type) “Adenomas and Adenocarcinomas”. (2) The keywords of files are (Data Category) “Transcriptome Profiling”, (Data Category) “Gene Expression Quantification”, (Experimental Strategy) “RNA-Seq”, (Workflow Type) “HTSeq-FPKM”, filter criteria of clinical data: (Data Category) “Clinical”, (Data Format) “BCR XML”) (Miao et al., 2020). We obtained 692 mRNA expression profiles of CRC. Among them, 51 (7.4%) samples were normal and 641 samples (92.6%) were cancerous. The TCGA-CRC cohort consisted of 624 patients. The clinicopathological features of the patients are listed in Table S1. The GEO dataset was chosen mainly as a validation dataset to prove the accuracy of analysis results based on TCGA dataset. The data access (GSE39582, GSE87211) was obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) for verification studies. The GSE39582 dataset consisted of 585 samples, including 19 colon mucosa samples and 566 stage I–IV colon adenocarcinoma samples. The GSE87211 dataset included 363 samples, including 160 rectal mucosa samples and 203 rectal tumor samples (Hu et al., 2018; Marisa et al., 2013). We merged the two datasets through a Perl script. Raw data from gene chips were normalized using the RMA algorithm provided by the R-package “limma” (Ritchie et al., 2015). Perl script and the R-package “sva” were used to merge microarray data and decrease heterogeneity between the two studies (Zhang et al., 2019).

Identification of differential expression RBPs

RNA binding proteins (RBPs) were obtained from the literature (Gerstberger, Hafner & Tuschl, 2014). Perl software was used to extract and integrate the data. We used the Wilcox Test in R-package “limma” to screen differentially expressed RBPs. False-positive discovery (FDR) < 0.05 and Log2 | (fold change, FC) | > 1 was set as the cutoffs. R-package “pheatmap” was used to analyze Bidirectional hierarchical clustering and draw a heat map.

Functional enrichment analysis of RBPs

Functional annotation and enrichment analysis of up-regulated and down-regulated RBPs were based on the R-packages “clusterProfiler”, “org.Hs.eg.db”, “ggplot2” and “enrichplot”, which classified Gene Ontology (GO) to the Biological Processes (BP), Cellular Components (CC) and Molecular Functions (MF). FDR < 0.05 was set as the cutoff. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis was performed by Gene Set Enrichment Analysis (GSEA). It was considered statistically significant that the number of permutations was set to 1,000, P-value < 0.05 and FDR < 0.25.

Protein–protein interaction network constructed and module selection

The STRING database (http://www.string-db.org/) (Von Mering et al., 2003) was used to evaluate the protein–protein interaction (PPI) information from differently expressed RBPs. Cytoscape software (version 3.7.2) was used to construct and visualize the PPI network. The essential genes and modules were selected in the PPI network using the Molecular Complex Detection (MCODE) plug-in with MCODE and node count numbers more significant than five (Bader & Hogue, 2003).

Construction prognosis model of RBPs and survival analysis

Univariate Cox regression analysis was used to screen out the RBPs that had a significant relationship with the OS of CRC patients (P < 0.05). We used multivariate Cox regression analysis to construct an optimal model of RBPs. We used Kaplan–Meier analysis and the log-rank test (P < 0.05) to assess survival. Each patient’s risk score was calculated using the following formula: Risk score = ${\sum }_{n=1}^n\mathrm{C}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\left(n\right)\times \mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\left(n\right)$, $\mathrm{C}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{i}\mathrm{c}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}\left(n\right)$ representing the regression coefficient and $\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}\left(n\right)$ indicating the relative expression levels of each prognosis-associated RBPs normalized by z-score in the prognostic risk score model, respectively (Liu et al., 2019; Zhang et al., 2020). Survival information was obtained from the clinical data of the aforementioned samples. The median risk score was chosen as the cutoff value for the CRC cohort dichotomy, and the CRC patients were categorized into high and low-risk groups. We used the R-packages “survminer” and “survival” to draw the survival curve according to the high and low-risk values. The Receiver Operating Characteristic (ROC) curve was used to examine the specificity and sensitivity of survival prediction by the protein characteristic risk score and was drawn using the R-package “survival ROC”. The area under curve (AUC) value is an indicator of prognostic accuracy (Sachs, 2017). We drew a heat map, risk curve and survival state diagram based on the different risk scores of our patients. Independent prognostic proteins were identified using univariate and multivariate Cox hazard regression analysis. A total of 758 CRC patient samples with prognostic information from the GSE39582 and GSE87211 datasets were used as verification cohorts to verify the prognostic model’s predictive ability. Protein co-expression was analyzed using the R-packages “ggalluvial”, “ggplot2” and “dplyr”.

Construction the nomogram

We used the expression level of RBPs by the R-package “rms” to construct the nomogram and predict the likelihood of OS. Calibration curves were used to estimate the consistency between predicted and actual survival, and the performance of the model in predicting prognosis was evaluated by the consistency index (C-index). C-index values of 0.5 and 1.0 represent a random probability and an excellent ability to predict survival, respectively.

Identification of differently expressed RBPs in CRC patients

One thousand, five hundred and forty-two RBPs were used in our study. One hundred and seventy-eight differently expressed RBPs were excluded, and 121 up-regulated and 57 down-regulated RBPs were included in our study (Fig. 2).

Identification of involved signaling pathways of the differently expressed RBPs

We divided the selected RBPs into up-or down-regulated expression groups to explore their function and mechanisms. Functional enrichment analysis was conducted using R software. The results showed that 105 GO terms of BP, 21 GO terms of CC and 48 GO terms of MF were significant in the up-regulated group (FDR < 0.05). Seventy-three GO terms of BP, three GO terms of CC and 25 GO terms of MF were significant in the down-regulated group (FDR < 0.05). The up-regulated RBPs were mainly enriched in ncRNA metabolic process during BP analysis, ncRNA processing, nucleic acid phosphodiester bond hydrolysis, RNA phosphodiester bond hydrolysis and ribonucleoprotein complex biogenesis (Fig. 3A). The down-regulated differently expressed RBPs were significantly enriched in mRNA processing, RNA splicing, defense response to a virus, response to a virus and translation regulation (Fig. 3B). During the CC analysis, the up-regulated RBPs were enriched in the nucleolar part, cytoplasmic ribonucleoprotein granule, ribonucleoprotein granule, nuclear envelope and chromatin (Fig. 3A). In contrast, the down-regulated RBPs were enriched in the endolysosome membrane, apical dendrite and endolysosome (Fig. 3B). In terms of MF, the up-regulated differently expressed RBPs were notably enriched during catalytic activity, when acting on RNA, and during nuclease activity, ribonuclease activity, endonuclease activity and tRNA binding (Fig. 3A). The down-regulated differently expressed RBPs were significantly enriched in mRNA binding, AU-rich element-binding, mRNA 3′-UTR AU-rich region binding, mRNA 3′-UTR binding double-stranded RNA binding (Fig. 3B). The top-10 terms of BP, CC and MF are shown in Fig. 3. The detailed information of GO the terms is shown in Table S2. According to the GSEA analysis, KEGG pathways of differently expressed RBPs showed that 66/165 gene sets were up-regulated in phenotype h; 35 gene sets were significant at FDR < 25% and 99/165 gene sets were up-regulated in phenotype l, and 55 gene sets were significant at FDR < 25%. The top-five up-regulated gene sets were basal cell carcinoma (NES = 2.15, P = 0.000), dilated cardiomyopathy (NES = 2.11, P = 0.000), glycosaminoglycan biosynthesis chondroitin sulfate (NES = 2.08, P = 0.000), complement and coagulation cascades (NES = 2.06, P = 0.002), and hypertrophic cardiomyopathy hcm (NES = 1.91, P = 0.002), while the top-five downregulated gene sets were RNA degradation (NES = −2.08, P = 0.008), valine leucine and isoleucine degradation (NES = −2.06, P = 0.004), ubiquitin-mediated proteolysis (NES = −2.02, P = 0.006), propanoate metabolism (NES = −1.95, P = 0.006), and oocyte meiosis (NES = −1.92, P = 0.014) (Figs. 3C–3L). The detailed GSEA analysis is shown in Table S3.

PPI network construction and hub modules screening

We constructed a PPI network with 147 nodes and 543 edges based on the STRING database’s data using Cytoscape software to further explore the roles of differentially expressed RBP in CRC (Fig. 4A). The co-expression network was constructed by MODE Plug-in to the study potential hub modules and the top-four crucial modules. The compulsory modules 1, 2, 3 and 4 consisted of 18 nodes and 144 edges; nine nodes and 28 edges; five nodes and 10 edges; and six nodes and 12 edges, respectively (Fig. 4B). RBPs in the compulsory modules 1–4 were mainly enriched in ribosome biogenesis, ncRNA processing, mRNA processing, mRNA binding, defense response to a virus, response to a virus, DNA alkylation and DNA methylation. The PPI network is shown in Table S4.

Construction prognosis model of RBPs and survival analysis

One hundred forty-seven differently expressed RBPs were identified from the PPI network. To further exploring the prognostic value of these RBPs, we carried out a univariate Cox regression analysis and screened out 15 potential prognostic-associated key RBPs, including six Low-Risk BAPs (HR < 1, P < 0.05) and nine High-Risk RBPs (HR > 1, P < 0.05, Fig. 5A). Nine RBPs were filtered out by multivariate COX regression analysis (Fig. 5B). The prognosis model was constructed based on these nine RBPs. The full name and coefficients of these proteins are shown in Table 1.

The risk score for each patient was calculated based on the following formula:

Risk score = 0.2943 × expression level of PNLDC1 + 0.4835 × expression level of RRS1 + 0.3609 × expression level of HEXIM1 + (−0.1757) × expression level of PPARGC1A + (−0.3660) × expression of PPARGC1B + (−0.4157) × expression level of BRCA1 + 0.5460 × expression level of CELF4 + 0.4375 × expression level of AEN + 0.1639 × expression level of NOVA1.

The risk score was applied to predict prognosis, and 586 CRC patients were divided into high-risk and low-risk subgroups according to the median risk score. Kaplan–Meier cumulative curves showed a significant difference in OS between the high-risk and low-risk groups. Patients with high-risk scores had a worse OS than those with low-risk scores (Fig. 5C). A ROC analysis was used further to assess the prognostic accuracy of the prognosis model. Our results indicated that the AUC value was 0.712 (Fig. 5D), which showed a moderate diagnostic ability. An elevated risk score was correlated with a decreased survival time and an increase in mortality (Fig. 5E). Age, clinical-stage, primary tumor size (T), regional lymph node involvement (N), distant metastasis (M) and risk score were significantly associated with OS in the univariate independent prognostic analysis (P < 0.001, Fig. 5F; Table 2). The multivariate independent predictive analysis demonstrated that age, T, and risk scores were also independent prognostic predictors (P < 0.05, Fig. 5G; Table 2). We assessed whether the nine-RBPs predicting model had an equivalent prognostic value in other CRC cohorts and used the same formula in the validation group in which data come from the GSE39582 and GSE87211 datasets. The result showed that patients with high-risk scores also had a worse OS than those with low-risk scores; the risk score was also an independent prognostic predictor in the GSE39582 and GSE87211 cohorts (Fig. 6; Table 3). These results showed that the prognostic model, based on the nine RBPs, had higher sensitivity and specificity.

Construction of a nomogram based on the nine crucial RBPs

A nomogram can quantitatively identify an individual’s clinical risk by combining multiple risk factors (Liang et al., 2015; Won et al., 2015). We used the nomogram to predict the probability of 1-, 3- and 5-year OS by merging the nine RBPs signatures for TCGA-CRC and verified these results using the GEO data sets (Fig. 7A). The distribution point of each RBPs was proportional to its risk contribution to the survival and was normalized to the distribution of 0–100. We calculated each patient’s total points by summarizing the number of points for all RBPs and estimated the survival rates by drawing a vertical line between each prognosis axis and the total point axis. This tool may help practitioners make clinical decisions for CRC patients. The calibration curve demonstrated that the actual survival rate matched the predicted 1-, 3- and 5-year survival rate. The C-index was 0.77 in the TCGA database (Fig. 7B). The nomogram was verified in the GEO cohorts, the C-index was 0.73, and the 1-, 3- and 5-year calibration curves are shown in Fig. 7C.

Co-expression analysis of RBPs

To explore more potential prognosis proteins, we used the nine RBPs in the model for co-expression analysis. CELF4 and LUZP4 were co-expressed; NOVA1 was co-expressed with RBFOX3, ELAVL4, DZIP1, SAMD4A, ZCCHC24, AFF3, RBM20, ENOX1, RBPMS2 and CPEB1. PPARGC1A was co-expressed with RBM47, RAVER2 and CPEB3. PPARGC1B was co-expressed with RBM47, RAVER2 and LRRFIP2. RRS1 was co-expressed with BOP1 (∣correlation coefficient∣＞0.6, P＜0.001, Fig. 8). The co-expression analysis is shown in Table S5.

Validation of the expression of hub RBPs

The IHC results of the HPA database were used to explore the expression of the nine hubs RBPs in CRC, and the result showed that RRS1 and BRCA1 were highly expressed in the CRC tissues compared with normal tissues. However, PNLDC1 was more highly expressed in the normal colorectal tissues than in the cancer tissues. HEXIM1 was highly expressed in both normal and colorectal tissues. In contrast, CELF4 was not expressed in normal and colorectal tissues. Furthermore, the protein expression of NOVA1 was not significantly different between normal and CRC tissues (Fig. 9).