Identification of Prognosis-Related RNA-Binding Proteins to Reveal the Role of RNA-Binding Proteins in the Progression and Prognosis of Colon Cancer

Introduction

Colon cancer, a major malignancy of the alimentary canal, ranked third among malignant tumors in terms of morbidity worldwide.¹ A relevant study revealed that more than one million people develop colon cancer each year, and the disease-specific mortality rate in developed countries was approximately 33%.² The mortality of colon cancer is on the rise due to changes in diet and lifestyle.³ Although colon cancer treatment options (eg, surgery, chemoradiotherapy, and immunotherapy) have greatly improved, the 5-year survival rate remains in a low level.

RNA binding proteins (RBPs) are a variety of proteins that interact with RNA, which are widely expressed in cells.^4–6 Through high-throughput screening, 1542 RBPs were identified, accounting for 7.5% of all protein-coding genes.⁷ These RBPs affect post-transcriptional events and regulate the physiological events of cells, thus, involved in many biological processes, such as RNA splicing, mRNA stabilization, and protein translation.^4,8 As RBPs play a variety of key functions in post-transcriptional events, the changes of RBPs are related to the occurrence and development of many human diseases. Some RBPs were also aberrantly changed in tumors, which influenced the translation of mRNA to protein level. For instance, due to the silence of miR-34a, dysregulated expression of IGF2BP3 accelerates gastric carcinogenesis.⁹ QKI-5 regulates cell proliferation in lung cancer by changing cancer-associated alternative splicing.¹⁰ However, the role of RBPs in the development of colorectal cancer remains unclear.

The Cancer Genome Atlas (TCGA) is the largest cancer database, containing samples of more than 20,000 primary cancers and normal matched samples of multiple cancer types. Therefore, we can investigate tumor gene data in greater depth with bioinformatics methods. Furthermore, it can be linked to clinical data to obtain more valuable and meaningful results. Here, we described the expression and lineage of RBPs in colon cancer and investigated the mechanisms of RBPs in the development and prognosis of colon cancer.

Methods and Materials

Data Acquisition

We identified and downloaded the transcriptome data of patients with colon cancer from the TCGA database through the R package, including 41 cases of paracancerous normal tissue and 473 cases of tumor tissue. Furthermore, relevant clinical information of the 473 colon cancer patients was obtained, such as age, gender, stage, tumor, lymph node and metastasis stage, survival status, and survival duration (Table 1). Finally, the “Limma” package in R software was utilized to correct the downloaded transcriptome data.

Table 1 Clinical Characteristics of Included Patients in the Study

Gene Function Enrichment Analysis

In order to explore the major biological process of selected hub genes, methods were utilized to conduct the gene functional enrichment analyses, including Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene ontology (GO). We utilized the Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) to identify enriched KEGG and GO themes.

Survival Analysis and Expression Comparison of Hub Genes

Clinical information for TCGA colon cancer, including survival time, survival state, and tumor, nodes, metastases (TNM) staging, were also downloaded from the TCGA database (samples with missing information were excluded). The survival R package was applied in survival analyses for hub genes. For the overall survival rates, the Log rank test was used to detect significant differences. The results were visualized using Kaplan–Meier survival curves, and a P-value <0.05 was considered statistically significant.

Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) (version 3.0, The Broad Institute of MIT and Harvard, http://software.broadinstitute.org/gsea/downloads.jsp) was conducted between colon cancer and paracancerous normal tissues to study the biological characteristics of colon cancer. In detail, the “collapse data set to gene symbols” was set to false, the number of marks was set to 1000, the “permutation type” was set to phenotype, the “enrichment statistic” was set to weighted, and the Signal-to-Noise metric was used for ranking genes. The high expression group was used as the experimental group, and the low expression group was used as the reference group. The “c2.cp.kegg.v7.0.symbols.gmt” gene sets database was used for enrichment analysis. A gene set size >500 and <15, false discovery rate (FDR) <0.25, and nominal P-values <0.05 were regarded as the cut-off criteria.

Statistical Analysis

All analyses were performed using R 3.6.1. All statistical tests were two-sided, and a P-value <0.05 was considered statistically significant. Continuous variables that conformed to the normal distribution were compared with an independent t-test for comparison between groups, while continuous variables with skewed distribution were compared with the Mann–Whitney U-test. The correlation matrix was constructed by R software based on Pearson’s Correlation Coefficient. The relationship between immune cell infiltration and overall survival was analyzed through the Kaplan–Meier curve, which was evaluated by the Log rank test. Time-dependent receiver operating characteristics (ROC) curves were used to analyze the sensitivity and specificity of the recurrence prediction model. The univariate and multivariate regression models were used to analyze the effects of individual variables on survival. The nomogram was constructed with the regression coefficients based on the Cox analysis.

Results

Differential Expression Screening of Colon Cancer

The transcriptome data of 473 colon cancer cases and 41 adjacent normal tissue cases were obtained from the TCGA database for differential expression analysis. A total of 473 RBPs were identified as differentially expressed RBPs (DERBPs) between colon cancer and normal tissues, including 321 upregulated and 152 downregulated (P < 0.05, Figure 1A, Table 2). The heatmap of the top 10 upregulated and top 10 downregulated DERBPs is shown in Figure 1B.

Table 2 Univariate Cox Regression Analysis to Screen RNA Binding Proteins Associated with Overall Survival in Colon Cancer Patients

Figure 1 Identification of DERBPs. (A) volcano plots of 473 DERBPs in colon cancer and normal tissues from TCGA database. (B) Heatmap plots of top 10 up-regulated and top 10 down-regulated DERBPs. The colors in the heatmaps from green to red represent expression level from low to high. The red dots in the volcano plots represent up-regulation, the green dots represent down-regulation and black dots represent genes without differential expression.

Functional Annotation of the 473 DERBPs

In order to fully understand the biological attributes of these 473 DERBPs, we conducted KEGG and GO analysis. Based on the results of DAVID, the top three enriched GO terms for upregulated genes and downregulated genes were ncRNA processing, ribosome biogenesis and rRNA metabolic processing, RNA splicing, regulation of translation and regulation of cellular amide metabolic processing (Figure 2A). The top biological pathways enriched for upregulated genes and downregulated genes were ribosome biogenesis in eukaryotes, RNA transport and spliceosome, spliceosome, and RNA transport and influenza A pathways (Figure 2B).

Figure 2 GO (A) and KEGG (B) enrichment analysis of DERBPs.

Establishment of Immune Prognosis Model

For the purpose of revealing the relationship between the 473 DERBPs and overall survival, 25 prognostic DERBPs were identified by utilizing univariate Cox regression analysis (Figure 3A, Table 2). TCGA colon cancer data were randomly divided into two sets (training set: validation set, 3:2). Then, multivariate Cox regression analysis was applied to select the independent indicators for the overall survival based on the training set and finally got 10 DERBPs (PPARGC1A, ZNF385A, SRP14, RIOK1, ERI1, NOL3, RRS1, TDRD6, AEN, and PNLDC1) for the construction of the prognostic index (Figure 3B, Table 3). After constructing the prognostic index, patients were separated into high-risk and low-risk groups (Figure 3C and D). A heatmap was utilized to visualize the difference between gene expression profiles in low- and high-risk patients in the colon cancer training set (Figure 3E). The results from the K–M analysis indicated that high-risk patients had lower overall survival than low-risk patients in both the training and validation groups (P < 0.001) (Figure 4A–C). The ROC curve revealed that the risk model had good sensitivity and specificity in predicting survival risk (area under the ROC curve (AUC) = 0.782, AUC = 0.744 for five years overall survival in the training and validation groups, respectively) (Figure 4D–F). In order to explore whether the constructed RBP risk scoring model was independent of age, gender, stage, and other clinical-pathological parameters, univariate and multivariate Cox regression analyses were performed for age, gender, stage, TNM, and risk score. In the univariate Cox model, age, pathological stage, pathological T, N, M stage, and high-risk score were associated with poor survival (Figure 5A). In the multivariate Cox model, only age, stage, and risk score worked as independent predicted factors (Figure 5B). To better predict the prognosis of patients with colon cancer at three and five years post-surgery, a new nomogram was constructed from the variables associated with Overall Survival (OS) (age, pathological stage, TNM stage, and risk score) (Figure 5C–E).

Table 3 Multivariate Cox Regression Analysis to Establish RNA Binding Proteins Risk Prediction Model

Figure 3 (A) Univariate survival analysis by cox proportional hazards models to select prognostic key RBPs. (B) Multivariate cox regression model for 10 prognostic RBPs used to construct RBPs risk score model. (C) Distribution of immune risk scores in colon cancer patients. (D) Distribution of survival status in colon cancer patients. (E) Distribution of specific risk factors in the high- and low-risk groups (divided by median value).

Figure 4 (A) Kaplan-Meier curve analysis of high-risk and low-risk patients in the training cohort. (B) Kaplan-Meier curve analysis of high-risk and low-risk patients in the testing cohort. (C) Kaplan-Meier curve analysis of high-risk and low-risk patients in the entire TCGA cohort. (D) Time dependent ROC curve analysis of the training cohort. (E) Time dependent ROC curve analysis of the testing cohort. (F) Time-dependent ROC curve analysis of the entire TCGA cohort.

Figure 5 Cox’s proportional hazard model of correlative factors in colon cancer patients. (A) Univariate COX regression analysis for seven clinicopathological parameters affecting the overall survival. (B) Multivariate COX regression analysis for seven clinicopathological parameters affecting the overall survival. (C) An established nomogram to predict colon cancer survival based on cox model. (D and E) Plots displaying the calibration of each model comparing predicted and actual 3- and 5-year overall survival.

Clinical and Prognostic Correlation of 10 Model Genes and RBP Risk Scores

We further investigated the proportion of each model gene in different pathological stages. We demonstrated that SRP14, PPARGC1A, and ER1 were most significantly associated with the development of colon cancer (Figure 6). In regard to the immune genes risk score, a strong correlation with the pathological, clinical N, and clinical M stages was identified (Figure 7).

Figure 6 Correlation analysis between TNM&Stage and 10 model genes in colon cancer cases. (A) Correlation analysis between tumor stage and 10 model genes expression in colon cancer cases. (B) Correlation analysis between node stage and 10 model genes expression in colon cancer cases. (C) Correlation analysis between metastasis stage and 10 model genes in colon cancer cases. (D) Correlation analysis between pathologic stage and 10 model genes expression in colon cancer cases.

Figure 7 Correlation between RBPs risk scores and various clinical factors. (A) Age. (B) Gender. (C) Stage. (D) T stage. (E) N stage. (F) M stage.

Gene Set Enrichment Analysis of Risk Scores

A GSEA of risk scores based on the TCGA colon cancer cohort was performed to explore the biological relevance of risk scores involved in colon cancer progression. GSEA analysis indicated high-risk scores were associated with MYC_TARGETS_V2, UV_RESPONSE_UP, and Wnt/beta-catenin signaling pathways (Figure 8A–C). In addition, low-risk scores were associated with the protein secretion pathway (Figure 8D).

Figure 8 Gene set enrichment analysis of RBPs risk scores. (A–C) high risk scores. (D) low risk scores.

Discussion

Malignant tumors are characterized by uncontrolled cell growth, mainly due to the dysregulated expression of cancer driver genes that regulate cell proliferation.¹¹ Post-transcriptional mechanisms can greatly influence the gene expression patterns in cancer cells, in which RBPs play key roles. They can interact with target mRNAs in a sequence-dependent and structure-dependent manner and determine cellular behavior by manipulating the processing of these mRNAs.⁴ It has been reported that RBPs show dysregulated expression in various human cancers. However, little is currently known about the expression patterns and roles of RBPs in colon cancer. Therefore, it was crucial to investigate the RBP subsets to evaluate risk and tumor prognosis in colon cancer.

In our study, we conducted a comprehensive and detailed assessment of RBPs in colon cancer based on the data from a large set of samples. All gene expression data and patients’ clinical characteristics information were downloaded from the TCGA dataset. We analyzed the 1542 RBPs between colon cancer and normal tissues. Eventually, we verified 473 differentially expressed RBPs. Moreover, we identified and constructed a 10 hub RBP risk score model for colon cancer via univariate and multivariate Cox regression analyses, including PPARGC1A, ZNF385A, SRP14, RIOK1, ERI1, NOL3, RRS1, TDRD6, AEN, and PNLDC1. Furthermore, to investigate the prognostic value of the model, we performed the ROC curve and investigated the association between the model and clinical features. As expected, the high-risk group was correlated with worse overall survival and was inclined to have advanced stages and higher histological grades, which might manifest a poor outcome.

According to the results of the biological functions and pathway enrichment analysis of these differentially expressed RBPs, the upregulated RBPs were significantly enriched in ribosome biogenesis in eukaryotes, RNA transport, and spliceosome, while downregulated differentially expressed RBPs were enriched in spliceosome, RNA transport, and influenza A pathways. Recently, mechanistic data have emerged, suggesting a broader role for dysregulated ribosome biogenesis in the development and progression of most spontaneous cancers.¹² Moreover, it is reported that RNA polymerases are consistently dysregulated in cancer, which is mostly mediated through upstream oncogenetic and tumor-suppressive signaling pathways rather than through mutations.¹³ These results suggest that RBPs can affect the growth of tumor cells by regulating multiple biological pathways.

Several genes in the RBPs model have been investigated in human cancers. RIOK1 activates NF-KB signal transduction, which promotes cell cycle progression and tumor lung colonization in vivo. It is demonstrated that RIOK1 is overexpressed in different subtypes of human lung cancer and breast cancer, suggesting that RIOK1 is a potential therapeutic target, especially in Ras-driven cancer.¹⁴ Ribosome biogenesis regulator 1 protein homolog (RRS1) is an important factor in ribosome biogenesis. At the molecular level, RRS1 silencing decreased the expression of M-phase inducer phosphatase 3 (CDC25C), cyclin-dependent kinase 1 (CDK1), and antigen Ki-67 (Ki-67), and increased the protein levels of cyclin-dependent kinase inhibitor 1 (CDKN1A) and tumor suppressor p53 (p53). In conclusion, RRS1 may promote the development of colon cancer. Therefore, targeting RRS1 may be a promising treatment strategy for patients with colorectal cancer (CRC).¹⁵

However, there were some limitations in our research. First, the sample size in our study was small, and a larger cohort and more abundant sequencing results are needed. Second, we only focused on the gene expression and gene mutation levels but ignored other events, such as gene methylation and copy number amplification, which are also important in tumor progression. Third, the effect of RBPs on the progression and prognosis of colorectal cancer needs to be verified in vivo and in vitro. Finally, our findings are significant but may not be useful to change the current clinical practice. They can be considered a first step towards a promising treatment strategy for patients with CRC. Further investigation is needed by a prospective study.

In summary, our study sheds light on the utility of RBPs in the prognosis of colon cancer. The constructed RBPs risk scoring model is reliable in predicting the prognosis of colon cancer, and this risk scoring model is an independent influencing factor for the prognosis of colon cancer. With the rapid development of high-throughput technology, we are confident that our risk scoring model has great potential in clinical practice.