Development and Validation of a Prognostic Model of RNA-Binding Proteins in Colon Adenocarcinoma: A Study Based on TCGA and GEO Databases

Introduction

Colon adenocarcinoma (COAD) is one of the leading causes of cancer morbidity and mortality worldwide. It is the third and second most prevalent cancer among men and women, respectively.¹ Although the survival rate of patients with COAD has improved due to advances in surgical techniques,² the prognosis of elderly colon cancer patients remains poor due to the increased risk of postoperative complications and postoperative mortality.³ Therefore, to elucidate more effective therapeutic strategies, it is critical to clarify the molecular mechanisms of COAD pathogenesis and progression.

RNA-binding proteins (RBPs) are proteins that bind to various RNAs, such as rRNA, miRNA, snRNA, ncRNA, mRNA, snoRNA and tRNA. To date, more than 1500 RBP-encoding genes have been experimentally validated in the human genome.⁴ RBPs play a key role in RNA processing by mediating posttranscriptional regulation; for example, these proteins can regulate mRNA stability, localization, variable splicing, polyadenylation and translation efficiency.⁵ Posttranscriptional regulation is essential in life processes. Dysregulated RBP expression is associated with many human diseases.^5,6 RBPs have been reported to contribute to the development and progression of cardiovascular disease through their involvement in posttranscriptional regulation.⁷ Some studies have reported that RBPs are expressed in various cancers, which affects mRNA-to-protein translation and is involved in carcinogenesis.^8,9 Currently, there are only a few reports on the critical role of RBPs in cancer development and progression.^10–13 Therefore, a comprehensive analysis is necessary to further understand the function and role of RBPs in cancer.

In the present study, we performed a detailed analysis based on RNA expression sequences and the corresponding clinical data of COAD downloaded from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. We also screened prognosis-related gene signatures by applying consensus clustering analysis, least absolute shrinkage selection operator (LASSO) regression analysis and Cox regression analysis. Moreover, we developed a prognostic model to predict the prognosis of colon adenocarcinoma and its therapeutic targets.

Materials and Methods

Data Collection and Preprocessing

We obtained RNA sequence information for 473 COAD and 41 normal colon tissues and the clinical information of the corresponding patients from the TCGA database as a training set. We performed differential analysis using the R package “limma” and the Wilcox test to identify differentially expressed (DE) RBPs between COAD and normal colon tissues. The threshold values were| log2 fold change (FC)| ≥ 0.5 and adj P value < 0.05.

GO and KEGG Pathway Enrichment Analysis of the DE RBPs

GO analysis consists of 3 main components: biological processes (BPs), cellular components (CCs), and molecular functions (MFs).¹⁴ KEGG analysis provides a biological interpretation of genomic sequences and other high-throughput data.¹⁵ We used the R package “clusterProfiler” for GO and KEGG enrichment analysis of the DE RBPs, and set significance thresholds of P and false discovery rate (FDR) less than 0.05.

Protein-Protein Interaction (PPI) Network Construction and Module Acquisition

We investigated the interactions between the associated proteins of these DE RBPs in the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database.¹⁶ Then, Cytoscape (version 3.7.2) was used to construct the PPI networks.¹⁷ Moreover, we used the Cytoscape plugin Molecular Complexity Detection (MCODE) to obtain the more important modules in the PPI networks.

Identification of Prognosis-Related RBPs

We performed univariate Cox regression analysis to identify prognosis-related RBPs among DE RBPs using the R package “survival”, followed by LASSO regression analysis with a P value < 0.01. We used the R package “glmnet” to screen for prognosis-related RBP-encoding genes. After that, we performed multivariate Cox regression analysis to screen the best RBP-encoding genes associated with prognosis and obtained their standardized regression coefficients. In addition, risk scores were calculated by the following formula: risk score = expression of gene 1 × coefficient of gene 1 + expression of gene 2 × coefficient of gene 2 + … expression of gene N × coefficient of gene N.^18,19 Based on median risk scores, we divided all COAD patients into high-risk and low-risk groups. The difference in survival between the two groups was then analyzed by the R package “survival” using Kaplan-Meier curve analysis and Log rank test. In addition, receiver operating characteristic (ROC) curves were plotted and the area under the curve (AUC) values were calculated for 1, 3, and 5 years using the R package “survival ROC” to assess predictive ability.²⁰

Internal and External Validation of Prognosis-Related RBPs

We first performed internal validation by first randomly dividing the COAD samples within the TCGA database into two equal numbers, one as a training cohort and the other as a validation cohort. In addition, we further validated the prognosis-related RBPs by combining the GSE17538, GSE39084, GSE39582, and GSE103479 datasets from the GEO database into one dataset as an external validation cohort. Inclusion criteria for the GEO database dataset were a sample size greater than 50 and complete survival and clinicopathological data, including overall time, survival status, age, sex, and tumor stage. We then used the same formula as above to calculate the risk score for each COAD patient.

Establishment of the RBP Nomogram

According to the results of the multivariate Cox analysis, we used the R package “rms” to predict overall survival (OS) in the TCGA cohort of COAD patients at 1, 3 and 5 years, and based on this, we generated a prognostic nomogram of prognosis-related RBP-encoding genes and plotted calibration curves to assess the prognostic performance of the nomogram.

Verification of Expression Level and Prognostic Significance

We examined the expression of hub RBPs at the translational level with the Human Protein Atlas (HPA) online data.²¹ To demonstrate the relationship between DERBPs and OS, we plotted Kaplan-Meier survival curves. Finally, we performed the Log rank test to test the significance of the difference between the two groups.

Cell Lines and Cell Culture

COAD cell lines SW480, RKO, HCT116, HCT-15, DLD-1, HT-29, HCT-8 cells and human colonic epithelial cell (NCM460) were purchased from ATCC (American Type Culture Collection, Manassas, VA, USA). All COAD cell lines were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM, Gibco) supplemented with 10% fetal bovine serum (FBS, Gibco), 100 U/mL penicillin, and 100 μg/mL streptomycin at 37 °C in 5% CO2.

RNA Extraction and Quantitative Real-Time PCR (qRT-PCR)

Total RNAs were extracted from cell lines with TRizol Reagent (Invitrogen). Total RNAs were reversely transcribed into cDNAs with PrimeScript RT Master Mix (Takara) and then used to perform quantitative real-time PCR (qRT-PCR) with SYBR qPCR Master Mix (Vazyme). GAPDH was used as an internal control for gene quantification. The 2−ΔCT was calculated for every sample and normalized to GAPDH. The primer sequences used were shown in Table 1.

Table 1 Primer Sequences Used in the qRT-PCR Assay

Results

DE RBPs Screening

We screened 496 DE RBPs, including 346 upregulated RBPs and 150 downregulated RBPs, between COAD tissues and normal colon tissues by using the R package “limma” and the Wilcox test. The heat map and volcano plot are shown in Figure 1.

Figure 1 Identification of differentially expressed RNA-binding proteins (DE RBPs) in The Cancer Genome Atlas (TCGA) dataset. (A) Heat map of the DE RBPs based on their log2-transformed fragments per kilobase of transcript per million mapped reads (FPKM) expression values. Red represents high expression, and green represents low expression. (B) Volcano plot of DE RBPs between colon adenocarcinoma (COAD) and normal colon tissue; 208 were upregulated, and 122 were downregulated. Red: upregulated RBPs; black: unchanged RBPs; green: downregulated RBPs.

GO and KEGG Pathway Enrichment Analysis of the DE RBPs

With the R package “clusterProfiler”, we performed GO and KEGG pathway enrichment analysis of the DE RBPs in COAD. The DE RBPs were significantly enriched in ncRNA processing, ribosome biogenesis, RNA catabolic processes, and RNA splicing in the BP analysis; cytoplasmic ribonucleoprotein granule, ribonucleoprotein granule, ribosome and spliceosomal complex in the CC analysis; and catalytic activity, acting on RNA, nuclease activity, ribonuclease activity and mRNA 3ʹ-untranslated region (-UTR) binding in the MF analysis (Figure 2A). According to the KEGG analysis, DE RBPs were significantly enriched in pathways related to RNA transport and degradation, mRNA monitoring, spliceosomes, ribosomes, ribosome biogenesis in eukaryotes, and coronavirus disease-COVID-19 (Figure 2B).

Figure 2 Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed RBPs (DE RBPs). (A) Bubble plot of enriched GO terms. (B) Bubble plot of enriched KEGG pathways.

PPI Network Construction and Key Modules Selection

We established a PPI network using Cytoscape software based on the data in the STRING database; the network, included 481 nodes and 8223 edges (Figure 3A). We used the MODE tool to screen the key modules and obtained the most significant modules, which consisted of 69 nodes and 2279 edges (Figure 3B).

Figure 3 Protein-protein interaction (PPI) network and modules analysis. (A) PPI network for differentially expressed RBPs (DE RBPs). (B) Key module 1 in the PPI network.

Selection of Prognosis-Related RBPs

We identified 481 DE RBPs from the PPI network in the TCGA cohort. We performed univariate Cox regression analysis to investigate the prognostic significance of these DE RBPs and obtained 11 hub RBPs associated with prognosis (Figure 4A). Subsequently, we analyzed these 11 prognostically relevant pivotal candidate RBPs by LASSO regression analysis (Figure 4B and C) and multivariate Cox regression analysis and found that seven pivotal RBPs were independent predictors of prognosis in COAD patients (Figure 4D).

Figure 4 Selection of prognosis-related RBPs in the training cohort. (A) Univariate Cox regression analysis. (B and C) Least absolute shrinkage and selection operator (LASSO) regression analysis. (D) Multivariate Cox regression analysis to screen out the key RNA-binding proteins (RBPs) most relevant to prognosis.

Construction of the RBP-Related Risk Score Model

We developed a risk score model based on seven previously obtained key RBP-encoding genes that are associated with prognosis and calculated the risk score for each COAD patient using the following formula: risk score = (−0.066 × RPL10L Exp) + (−0.39811 × ERI1 Exp) + (−0.52662 × POP1 Exp) + (0.247129 × CAPRIN2 Exp) + (−0.36382 ×TDRD7 Exp) + (−0.44221× SNIP1 Exp) + (−0.19892×PPARGC1A Exp). Of the 7 RBPs screened for prognostic relevance, CAPRIN2 was a risk factor (HR>1), while RPL10L, ERI1, POP1, TDRD7, SNIP1 and PPARGC1A were protective factors (HR < 1) (Table 2). We divided the 446 COAD patients into a high-risk group (n = 223) and a low-risk group (n = 223) according to the median risk score of the training cohort. Survival analysis revealed significantly longer OS in the low-risk group than in the high-risk group (p = 6.673e-08) (Figure 5A). We performed a time-dependent ROC analysis to assess the prognostic ability of seven RBP biomarkers. The AUCs of this RBP risk score model for 1, 3 and 5 years were 0.694, 0.709 and 0.665, respectively. (Figure 5B), indicating moderate diagnostic performance. For the low-risk and high-risk groups, we also show expression heat maps, patient survival status, and risk scores for the signature consisting of seven RBPs (Figure 5C). In addition, we used the same formula for the TCGA internal dataset and the GEO external validation dataset to build predictive models to assess whether these 7 RBPs have prognostic value in the TCGA database (Supplementary Figures 1 and 2) and other COAD patient cohorts (Figure 6). The results of the above analysis are similar to those of the TCGA external training cohort analysis. The above results indicate that the prognostic model we developed has good sensitivity and specificity.

Table 2 Seven Prognosis-Associated Hub RNA-Binding Proteins (RBPs) Identified by Multivariate Cox Regression Analysis

Table 3 The Prognostic Value of Different Clinical Parameters in the Cancer Genome Atlas (TCGA) Cohort

Figure 5 Risk score analysis of the seven-gene prognostic model in The Cancer Genome Atlas (TCGA) cohort. (A) Survival curve for low-risk and high-risk groups. (B) Receiver operating characteristic (ROC) curves for predicting overall survival (OS) based on the risk score. (C) Expression heat map, risk score distribution, and survival status.

Figure 6 Risk score analysis of the eight-gene prognostic model in the GEO database combined cohort. (A) Survival curve for low-risk and high-risk groups. (B) Receiver operating characteristic (ROC) curves for predicting overall survival (OS) based on the risk score. (C) Expression heat map, risk score distribution, and survival status.

Relationship Between Clinical Features and the Seven Prognosis-Related RBP-Encoding Genes

We performed univariate, LASSO regression and multivariate Cox regression analyses of clinical characteristics in the training and testing sets. Age (p < 0.001), tumor stage (p = 0.013) and risk score (p < 0.001) were found to be independent prognostic indicators in the TCGA cohort (Table 3), while age (p < 0.001), tumor stage (p < 0.001) and risk score (p < 0.001) were independent prognostic indicators in the GEO database combined cohort (Table 4). Moreover, results from the TCGA internal training and validation cohorts also showed that age, tumor stage, and risk score were independent prognostic factors (Supplementary Tables 1 and 2).

Table 4 The Prognostic Value of Different Clinical Parameters in the GEO Database Combined Cohort

Construction of the RBP Nomogram

We integrated the 7-RBP signature with clinical features to create a nomogram to develop quantitative methods for assessing COAD prognosis (Figure 7A). The ability to plot a vertical line between the total point axis and each prognostic axis to calculate the survival of patients with COAD at 1, 3, and 5 years can be helpful for practitioners involved in making clinical decisions about patients with COAD. The nomogram calibration curves show that that we plotted the nomogram with better accuracy (Figure 7B).

Figure 7 Nomogram (A) for predicting the 1-, 3-, and 5-year overall survival (OS) of colon adenocarcinoma (COAD) patients in The Cancer Genome Atlas (TCGA) cohort. (B) Calibration curves for Nomogram.

Expression of the Seven Prognosis-Related RBP-Encoding Genes and Their Survival Impact

The expression of seven key RBPs was further analyzed with the HPA database, which contains immunohistochemical results for five RBPs (not including PPARGC1A and POP1) in COAD tissues and normal colon tissues (Figure 8). To validate the prognostic value of key RBP-encoding genes in the TCGA cohort, we analyzed the seven RBP-encoding genes by Kaplan-Meier curve analysis and showed that patients with COAD with low CAPRIN2 expression and high POP1 and PPARGC1A expression had higher OS rates than those with other expression patterns (Figure 9).

Figure 8 Verification of hub RNA-binding protein (RBP) expression in colon adenocarcinoma (COAD) tissues and normal colon tissues using the Human Protein Atlas (HPA) database. (A) CAPRIN2, (B) ERI1, (C) RPL10L, (D) SNIP1 and (E) TDRD7.

Figure 9 Kaplan-Meier curve analysis of the three prognostic RNA-binding proteins (RBPs), including CAPRIN2, POP1 and PPARGC1A, in The Cancer Genome Atlas (TCGA) cohort.

qRT-PCR Validation of the Seven Prognosis-Related RBP-Encoding Genes

qRT-PCR was used to validate the expressions of the seven prognosis-related RBP-encoding genes, including CAPRIN2, RPL10L, ERI1, POP1, TDRD7, SNIP1 and PPARGC1A in seven COAD cells (SW480, RKO, HCT116, HCT-15, DLD-1, HT-29 and HCT-8) and one human colonic epithelial cell (NCM460). The results showed that RPL10L, ERI1, POP1, TDRD7, SNIP1 and PPARGC1A were lowly expressed in most COAD cell lines compared to the control cell (NCM460), while CAPRIN2 was highly expressed (Figure 10). This is consistent with the prognostic significance.

Figure 10 qRT-PCR validation of seven prognosis-related RBP-encoding genes in colon adenocarcinoma (COAD) cell models compared to control cells. *p < 0.05. **p < 0.01. ***p < 0.001. Abbreviation: ns, not significant.

Discussion

Several studies have suggested that RBPs are closely associated with the development and progression of many cancers.^8,22 However, there are few comprehensive studies on the function and prognosis of RBP in patients with COAD. The main purpose of this study was to analyze the relationship between RBPs and the prognosis of patients with COAD. First, we identified 496 RBPs that differed between COAD tissues and normal colon tissues in the TCGA database. Then, we further analyzed the related biological pathways and constructed PPI networks for these RBPs. In addition, we identified a total of 11 survival-related RBPs significantly associated with COAD survival by univariate Cox regression analysis. Moreover, gene signatures were constructed by LASSO and multiple Cox regression analysis, and we identified seven hub RBP-encoding genes (CAPRIN2, RPL10L, ERI1, POP1, TDRD7, SNIP1 and PPARGC1A). We later performed survival analysis and ROC curve analysis of the hub RBPs to further investigate their clinical significance. Finally, we constructed a risk model based on seven hub RBPs associated with prognosis to predict the prognosis of COAD. Our findings are useful as they revealed novel biomarkers for the diagnosis and prediction of prognosis in COAD patients.

The GO functional pathway enrichment analysis indicated that the DE RBPs were greatly enriched in ncRNA processing, ribosome biogenesis, rRNA metabolic processes, RNA catabolic processes, RNA splicing, mRNA catabolic process, regulation of translation, and catalytic activity acting on RNA. In recent years, many studies have demonstrated the role of abnormal RNA metabolism and RNA processing in various diseases.^23–25 Posttranscriptional regulation of RNA stability is critical in regulating the expression of target genes. RBPs interact with RNA to form ribonucleoprotein complexes to increase the stability of mRNAs and thereby promote the expression of target genes, a process that is important in a variety of diseases. LncRNA OCC-1 regulates the levels of a large number of mRNAs directly related to colorectal cancer cell growth by regulating RBP HuR stability at the posttranscriptional level.²⁶ The RBP IGF2BP3/ELAVL1 complex is critical for regulating the stability of oncogenic mRNA in colorectal cancer.²⁷ RBP Msi1 inhibits the p21cip1 expression by directly binding to the shared sequence of p21cip1 3-’UTR in colon cancer cells.²⁸ In addition, ribonucleoprotein granules are key factors in protein biosynthesis. Alterations in ribonucleoproteins affect the translation process and are associated with tumor progression.²⁹ Analysis of the enriched KEGG pathways revealed that the aberrantly expressed RBPs regulate colon carcinogenesis and progression by affecting RNA transport, ribosome biogenesis, spliceosomes, mRNA monitoring pathways and RNA degradation.

We screened seven key prognosis-related RBPs in COAD patients, including CAPRIN2, RPL10L, ERI1, POP1, TDRD7, SNIP1 and PPARGC1A, by Cox regression analysis. Ribosomal protein L10 (RPL10) is an oncogene that encodes a protein known as oncoprotein QM.³⁰ The expression of RPL10 was reported to be negatively correlated with prostate adenocarcinoma grade.³¹ A previous study reported that RPL10 is specifically expressed in ovarian cancer cells and is considered a novel biomarker for diagnostic and/or therapeutic targeting of ovarian cancer.³² Exo-ribonuclease 1 (ERL1) is a ribonuclease involved in RNA metabolism, and it has been reported that ERL1 plays an important role in defining HOXC8 expression for the correct patterning of the skeleton.³³ However, there have been no studies of ERL1 in human cancer. “Pyrin-only” protein 1 (POP1, also known as POPDC1 and BVES) is a membrane protein essential for regulating tight junction (TJ) formation and blocking epithelial-mesenchymal transition (EMT), and through its epigenetic silencing, it may promote EMT in colon cancer.³⁴ Cytoplasmic activation/proliferation-associated protein 2 (CAPRIN2) is an LRP5/6 binding protein.³⁵ It has been reported that CAPRIN2 promotes the development of hepatoblastoma³⁶ and oral squamous cell carcinoma³⁷ progression by activating the classical WNT signaling pathway. Tudor domain-containing 7 (TDRD7) is a novel causative gene that can cause congenital cataracts such as azoospermia in human males.³⁸ Nevertheless, its role in cancer has also not been reported. Smad nuclear interacting protein 1 (SNIP1) is a major gene that regulates cyclin D1 mRNA stability during cotranscription or posttranscription stages, and its overexpression alone is sufficient to cause tumorigenesis.³⁹ Increasing research has indicated that SNIP1 might be a significant prognostic predictor in patients with cervical cancer,⁴⁰ non-small cell lung cancer,⁴¹ osteosarcoma,⁴² and tongue squamous cell carcinoma.⁴³ The above studies suggest that SNIP1 may serve as an important prognostic marker in cancer. Peroxisome proliferator-activated receptor coactivator 1A (PPARGC1A) might be implicated in the development of colorectal cancer.⁴⁴ Previous studies showed that PPARGC1A induces the expression of oxidative phosphorylation genes and promotes tumor growth in colon cancer.^45,46

We applied multivariate Cox regression analysis to construct a risk model consisting of these seven key RBPs that is valuable for predicting the prognosis of patients with COAD. Analysis of ROC curves over time showed that these seven genes had good diagnostic ability, and that the risk model we developed could be used effectively to identify COAD patients with a poor prognosis. However, to date, there have been few studies on the mechanisms by which these hub RBPs affect COAD pathogenesis. Therefore, further studies are very important. In addition, we constructed a nomogram capable of predicting OS in patients with COAD at 1, 3 and 5 years. We used Kaplan-Meier curves to evaluate the prognostic value of these seven pivotal RBPs and found that three of them (CAPRIN2, POP1, PPARGC1A) were associated with patient prognosis.

Nevertheless, there are still a few limitations of this study. First, the prediction model we developed was based on TCGA data only and was not clinically validated or prospectively studied; moreover, the limited clinical information in the TCGA dataset may reduce the accuracy of Cox regression analysis. Nevertheless, our predictive model constructed based on 7 RBPs showed great potential to predict the prognosis of COAD patients, which is crucial for clinical decision making.

Conclusion

In conclusion, we screened prognosis-related genes in COAD patients based on RBP-related genes. In addition, we validated the expressions of the seven prognosis-related RBP-encoding genes by qRT-PCR and constructed a prognosis-related nomogram for patients with COAD. The gene signatures we screened and the nomogram we created bring some enlightenment to the clinical. Nevertheless, the RBPs we screened still need to be further explored to determine if they might be useful for molecularly targeted therapy of COAD patients.