{'fraction_non_alphanumeric': 0.0385621704303023, 'fraction_numerical': 0.019643360302700964, 'mean_word_length': 5.3345232024477305, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 0, 'https://': 5, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Molecular characterization of colorectal cancer: A five-gene prognostic signature based on RNA-binding proteins

Background: Colorectal cancer (CRC) is one of the most common cancers worldwide. RNA-binding proteins (RBPs) regulate essential biological processes and play essential roles in a variety of cancers. The present study screened differentially expressed RBPs, analyzed their function and constructed a prognostic model to predict the overall survival of patients with CRC. Methods: We downloaded CRC RNA-sequencing data from the Cancer Genome Atlas (TCGA) portal and screened differentially expressed RBPs. Then, functional analyses of these genes were performed, and a risk model was established by multivariate Cox regression. Results: We obtained 132 differentially expressed RBPs, including 66 upregulated and 66 downregulated RBPs. Functional analysis revealed that these genes were significantly enriched in RNA processing, modification and binding, ribosome biogenesis, post-transcriptional regulation, ribonuclease and nuclease activity. Additionally, some RBPs were significantly related to interferon (IFN)-alpha and IFN-beta biosynthetic processes and the Toll-like receptor signaling pathway. A prognostic model was constructed and included insulin like growth factor 2 messenger ribonucleic acid binding protein 3 (IGF2BP3), poly (A) binding protein cytoplasmic 1 like (PABPC1L), peroxisome proliferator activated receptor gamma coactivator 1 alpha (PPARGC1A), peptidyl- transfer ribonucleic acid hydrolase 1 homolog (PTRH1) and tudor domain containing 7 (TDRD7). The model is an independent risk factor for clinicopathological characteristics. Conclusion: Our study provided novel insights into the pathogenesis of CRC and constructed a prognostic gene model, which may be helpful for determining the prognosis of CRC.


INTRODUCTION
Colorectal cancer (CRC) is the third most common cancer with more than 2.2 million new cases and 1.1 million deaths estimated to occur by 2030. [1] It is predicted that the incidence of CRC in young individuals will increase by as much as 90% by 2030. [2] There are no obvious symptoms in the early stage of CRC, but changes in stool habits, emaciation, hematochezia, anemia, colonic fistula or intestinal perforation can be found during the progression and distant metastasis of CRC. [3] Nearly 5% of CRC cases are genetic syndromes, including familial adenomatous polyposis (FAP), hereditary nonpolyposis colorectal cancer (HNPCC), or mutY deoxyribonucleic acid glycosylase-associated polyposis (MAP), and 15-20% of the cases have a family history of hereditary diseases. However, the majority of cases are sporadic diseases without a family history of CRC, which is closely associated with microsatellite instability and lifestyle factors. [4][5][6] Tumor recurrence is the primary barrier to improving overall survival. [7] To reduce the recurrence rate and mortality of patients with CRC, it is crucial to improve the surveillance ability following surgery. The traditional prognostic criterion for CRC is the American Joint Committee on Cancer/Union Internationale Contre le Cancer (AJCC/UICC) tumor-node-metastasis (TNM) staging system, but the system provides limited prognostic information and cannot predict the response to therapy due to significantly different clinical outcomes possibly emerging in the same TNM stages. [1] Therefore, some new methods to classify cancer, including molecular pathways, mutation status, cell origin, gene expression-based stratification and the TNM immune staging system, have been used. [1,8,9] In recent years, biomarkers have also been applied for early diagnosis and prognosis, and therefore, it is vital to find new biomarkers that reflect the occurrence, progression and prognosis of CRC.
RNA-binding proteins (RBPs) play essential roles in post-transcriptional events and regulate essential biological processes, including RNA localization, stability, transport, editing, degradation and translation, by forming ribonucleoprotein (RNP) complexes with target RNA. [10] There are 1,542 RBPs that have been experimentally validated, accounting for approximately 7.5% of all protein-coding genes. [11] RBPs are ubiquitous in cells and are evolutionarily conserved so that they can play a central role in essential biological functions. [11] Because the regulation of RBPs to RNAs occurs in rapid and efficient ways, small changes in RBPs can alter gene expression and cell fate. [12,13] Abnormalities in RBPs can induce diseases, including cancer, and correlate with patient prognosis. [10,[14][15][16] Previous studies have found that several RBPs are significantly related to CRC, such as lin-28 homolog B (LIN28B), human antigen R (HuR) and ribonucleic acid binding motif protein 3 (RBM3). LIN28B is overexpressed in nearly 30% of CRC, which indicates worse survival and recurrence. [17,18] The expression of HuR is significantly enhanced in colon tumors by regulating target transcribed mRNAs and is vital for neoplastic transformation and cancer development. [19] Over-expressed RBM3 in colorectal cancer cells increases proliferation and engenders chemotherapy resistance. [20] RBPs are closely related to CRC, and therefore, it is useful to perform a systematic analysis of RBPs in colon tumors.
With the development of high-throughput sequencing technology, a vast number of data sets have been generated and conserved in public portals, such as TCGA and Gene Expression Omnibus (GEO). However, the mining of these data is limited as it requires professional software and advanced expertise for researchers to process and analyze. [21] In the present study, we downloaded CRC data from the TCGA database and used bioinformatics to screen differentially expressed RBPs in tumor and normal tissue and further analyze their function in CRC. We also constructed a prognostic model to predict the overall survival of patients with CRC, which may have prognostic value in the clinic in future.

Data acquisition and analysis
The corresponding RNA-sequencing data containing 473 colon tumor samples and 41 normal colon samples were downloaded from the Genomic Data Commons data portal (https://portal.gdc.cancer.gov/). The R package "edgeR" was used to normalize the above data and analyze the expression of each RBP between tumor and normal samples. A false discovery rate (FDR) <0.05 and | log 2 fold change (FC)| ≥1 were set as filter conditions for differentially expressed genes (DEGs). The Gene Expression Omnibus series GSE17538 dataset was obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) and used as a validation cohort.

Gene Ontology, Kyoto Encyclopedia of Genes and Genomes analysis, and protein-protein interaction analysis
The differentially expressed RBPs were uploaded to the Database for Annotation Visualization and Integrated Discovery (DAVID) version 6.7 for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Both P < 0.05 and FDR < 0.05 were considered statistically significant. GO analysis included biological processes, cellular components and molecular functions. Subsequently, the STRING online tool (http://www.string-db.org/) was used to investigate the protein-protein interactions (PPIs) of differentially expressed RBPs. Visualization of the interaction network was performed by Cytoscape 3.8.0 software. Module analyses were performed by Molecular Complex Detection (MCODE). Modules with scores and node numbers greater than 5 were selected as key modules. The DAVID online tool and R package "ggplot2" were used to analyze and visualize GO terms and KEGG pathways of hub genes and key module genes.

Establishment of the prognostic gene signature for predicting the survival rate
Multivariate Cox regression analysis was performed to construct a prognostic gene signature. The result was shown as risk score = (Coefficient mRNA1 × expression of mRNA 1 ) + (Coefficient mRNA2 × expression of mRNA 2 ) + ... + (Coefficient mRNAn × expression of mRNA n ). The R packages "survival" and "survminer" were used to draw the Kaplan-Meier survival curve. To test the independence of the prognostic gene signature of other clinical parameters in TCGA, Cox regression analyses of clinical parameters were performed.

Validation of the gene change
The expression of prognostic genes was validated at the mRNA level by the Tumor Immune Estimation Resource database (TIMER) database (https://cistrome. shinyapps.io/timer/) and at the protein level by the Human Protein Atlas database (HPA) (http://www. proteinatlas.org). Additionally, genetic alterations of prognostic genes were determined by the cBioPortal online tool (http://www.cbioportal.org/).

Statistical analysis
In this study, we used SPSS 19.0 (IBM, Armonk, New York) and R software (R Foundation for Statistical Computing, Vienna, Austria) to process and analyze the data. For the Kaplan-Meier estimates, a significant difference in two-group survival curves was assessed by a log-rank test. P < 0.05 was considered statistically significant.

Identification of differentially expressed RBPs
To compare the differentially expressed RBPs (1,542 selected) between 41 normal colon samples and 473 tumor colon samples from the TCGA database, we used an R package to analyze the data. We obtained 132 differentially expressed RBPs (P < 0.05, |log 2 FC)| ≥ 1.0), including 66 upregulated RBPs and 66 downregulated RBPs [Supplementary Table S1]. We also constructed an expression heat map for all differentially expressed RBPs [ Figure 1].

Gene Ontology and KEGG pathway analysis of differentially expressed RBPs
To explore the functions and mechanisms of the above RBPs, the upregulated and downregulated RBPs were uploaded to DAVID for GO and KEGG pathway analysis, respectively. As shown in Table 1, for biological process (BP) analysis, upregulated RBPs were significantly enriched in RNA processing (ncRNA metabolic process, ncRNA processing, rRNA processing, rRNA metabolic process, tRNA processing and tRNA metabolic process), ribosome biogenesis, RNP complex biogenesis, RNA modification, post-transcriptional regulation of gene expression, RNA localization and regulation of translation. The downregulated RBPs were significantly enriched in the post-transcriptional regulation of gene expression, mRNA metabolic process and mRNA processing. For GO cellular component (CC) analysis [ Table 1], upregulated RBPs were significantly enriched in the nuclear lumen, nucleolus, intracellular organelle lumen, membrane-enclosed lumen, RNP complex, intracellular nonmembrane-bound organelle, nuclear pore and pore complex; but the downregulated RBPs were not significantly enriched in any cellular component. The GO molecular function (MF) analysis showed that the upregulated RBPs were significantly related to RNA binding, ribonuclease activity and nuclease activity, and downregulated RBPs were significantly related to RNA binding, including double-stranded RNA binding and single-stranded RNA binding [ Table 1]. KEGG pathway results showed that the upregulated RBPs were related to RNA degradation, and downregulated RBPs were related to the Toll-like receptor signaling pathway, but neither was significant [ Table 1].

Protein-protein interaction network construction
We uploaded these differentially expressed RBPs to the STRING database and constructed a coexpression network with 131 nodes and 411 edges. Then, the degree of each protein was calculated, and 10 genes were screened as hub genes according to the top eight scores of the degree, including nucleolar protein 56 (NOP56), ribosomal ribonucleic acid processing 9 (RRP9), bystin like (BYSL), dyskeratosis congenita 1 (DKC1), ribosome biogenesis regulator 1 homolog (RRS1), U3 small nucleolar ribonucleic acid-associated protein 14 homolog A (UTP14A), block of proliferation 1 (BOP1) and damage specific deoxyribonucleic acid binding protein 1 and cullin 4 associated factor 13 (DCAF13) [ Figure 2a]. We further selected three possible key modules with the MCODE plug-in with scores and node numbers greater than 5 [ Figure 2b-2d]. Then, GO and KEGG pathway analyses of 10 hub genes and key module genes were performed. As Figure 3 shows, for biological processes, these genes were significantly enriched in the ncRNA and rRNA processes, ribosome complex biogenesis, and regulation of IFN-alpha and IFN-beta biosynthetic processes. For cellular components, hub genes and selected module genes were significantly enriched in the nuclear lumen, nucleolus, intracellular organelle lumen, membrane-enclosed lumen and intracellular nonmembrane-bound organelle. For MF, these genes were significantly related to RNA binding and  Patients were divided into high-risk and low-risk groups with a cutoff of 0.99980 (the median of the risk score). The GSE17538 dataset was obtained to validate the risk score model. Each patient's risk score was calculated, and the cutoff of the high-risk and low-risk groups was − 0.63054 (the median of the risk score). The overall survival and expression heat maps of the high-risk and low-risk groups were plotted with R software [ Figure 4]. Compared with the low-risk group, the overall survival was significantly poorer in the high-risk group in TCGA (P < 0.001) and the GEO data set (P < 0.05).

Independent risk factors for clinicopathological characteristics to predict overall survival
To identify independent clinical risk factors, including age, sex, TNM pathological stage, primary site (left side and right side) and our risk score model, univariate and multivariate Cox regression analyses were performed with the TCGA and GEO data sets. The right side of the colon includes the cecum, hepatic flexure of colon, ascending colon and transverse colon, and the left side colon extends from the splenic flexure, descending colon and sigmoid. [22] As shown in Figure 5, our risk score model and TNM pathological stage were both independent prognostic factors for the overall survival.

External validation using online databases
To validate the expression levels and investigate genetic alterations of five prognostic model genes, HPA (http://www.proteinatlas.org/), TIMER (https://cistrome. shinyapps.io/timer/) and cBioPortal database (https:// www.cbioportal.org/) were used in this study. We selected 392 patients with colon cancer to investigate the genetic alterations of the five genes in the cBioPortal database. As shown in Figure 6a, PABPC1L possessed the most frequent genetic alterations (10%), and amplification mutation was the most common. In accordance with our results, the expression of IGF2BP3 and PTRH1 was significantly enhanced, and TDRD7 was decreased in the tumor samples compared with the normal samples in the HPA and TIMER databases [ Figures 6b and 7]. There was no immunohistochemistry information about the PABPC1L and PPARGC1A proteins in HPA. The TIMER database showed that the PABPC1L mRNA level was significantly increased and that PPARGC1A was decreased in the tumor group compared with the normal group [ Figure 7].

DISCUSSION
CRC is one of the most commonly diagnosed types of cancers but also the most preventable and treatable cancer. [23] Therefore, it is vital to screen for early CRC. Previous studies have shown that incidence and mortality steadily declined in countries with programmatic screening. [24,25] Colonoscopy may be the most effective strategy to screen for CRC, with 0.022 life-years gained and 1,068 CRCs prevented, but it is associated with rare and severe complications. [26,27] However, colonoscopy may be suitable as the second procedure following a positive first test. [28] Therefore, noninvasive tests are needed as supplements. Fecal deoxyribonucleic acid (DNA) testing is a new non-invasive strategy for screening, and Cologuard ® -targeted multiple DNA sequences, including mutant kirsten rat sarcoma viral oncogene homolog (KRAS), actin, fecal immunochemical test (FIT), aberrantly methylated bone morphogenetic protein 3 (BMP3) and the N-myc downstream regulated gene family member 4 (NDRG4) promoter regions, have been approved for CRC screening. [29] Compared with the FIT, a noninvasive screening method, Cologuard ® had a higher sensitivity (74% vs 92%) but lower specificity (95-96% vs 87-90%) for detection of CRC. [29] Therefore, finding more sensitive and specific biomarkers is essential for a fecal DNA test to diagnose CRC. Similar to molecular mutation and methylation, gene expression-based methods have also been used as tools for diagnosis, classification, treatment and prognosis. [8,[30][31][32][33] In this study, we used the gene expression data of RBPs from TCGA to identify the differentially expressed RBPs and analyze their function as well as construct a multiple gene signature predicting overall survival in patients with CRC.
We obtained 132 differentially expressed RBPs, including 66 upregulated and 66 downregulated RBPs. The biological processes of these DEGs were enriched in RNA processing and RNA metabolic process, ribosome biogenesis, RNP complex biogenesis, RNA modification, post-transcriptional regulation of gene expression and RNA localization. The gene mutations of RBPs changed their function, potentially inducing cancer. For example, dicer 1, ribonuclease III (DICER1), an endoribonuclease, processes precursor RNA molecules into mature forms, which is vital for miRNA and RNA interference biogenesis pathways. [34] The mutation of DICER1 impairs miRNA-mediated gene suppression, which can induce tumor susceptibility. [35] Post-transcriptional modifications frequently occur in RNA-binding elements within RBPs, which can change RBP binding properties, function and subcellular localization, and therefore, may be among the most critical mechanisms of RBP dysfunction in cancer. [36] For GO CC analysis, significantly differentially expressed RBP enrichment occurred in the nuclear lumen, nucleolus, intracellular organelle lumen, membrane-enclosed lumen, RNP complex, intracellular non-membrane-bound organelle, nuclear pore and pore complex. RBPs are ubiquitous in cells, and therefore, regulate a variety of biological functions. Nuclear maturation of transcripts, nuclear export of RNA, and RNA processing require nuclear RBPs. [37] In addition, several small nuclear ribonucleoproteins (U1, U2, U4, U5, and U6 snRNP) form a spliceosome to participate in pre-mRNA splicing, and abnormalities in RNA splicing frequently occur at the origin of many disorders, including cancer. [38,39] The RNP mutation regulates p53 expression, which induces many diseases, including cancer. [40] Furthermore, heterogeneous nuclear RNP D-like is aberrantly expressed in CRC, which promotes the growth of these cells by activating cell cycle progression. [41] RBPs can mediate organelle-coupled translation as post-transcriptional regulators. [42] The GO MF analysis showed that the differentially expressed RBPs were significantly related to RNA binding, ribonuclease activity and nuclease activity. RBPs can bind targeted RNAs by RNA-binding domains composed of 60-100 amino acids. As stated above, ribonuclease processes precursor RNA molecules into mature forms, and mutation of special ribonucleases, such as DICER1, induces cancer. [35] We screened eight hub genes, including NOP56, RRP9, BYSL, DKC1, RRS1, UTP14A, BOP1 and DCAF13, by constructing a protein-protein interaction (PPI) network and calculating the degree score. NOP56 is required for rRNA methylation, and abnormalities in NOP56 induce altered ribosome biogenesis and cell cycle progression. [43] Additionally, NOP56 is hyperactivated and essential for Myc-induced cell transformation in Burkitt's lymphoma. [44] Deacetylation of RRP9 is a prerequisite for pre-rRNA processing, [45] but there is no definite report about the relevance between RRP9 and cancer. BYSL participates in ribosomal processing of 18S rRNA and is significantly enhanced in hepatocellular carcinoma specimens and is essential for nucleogenesis in cancer cell proliferation. [46] DKC1 is remarkably increased in CRC samples, and DKC1 enhances angiogenesis by promoting hypoxia inducible factor-1 (HIF-1α) transcription and facilitates metastasis in CRC. [47] RRS1 is necessary for ribosome biogenesis and is significantly increased in CRC samples to promote cancer development. [48] UTP14A is essential for 18S rRNA processing and promotes p53 degradation by binding p53. [49] Nucleolar UTP14A is significantly higher in CRC samples than in normal samples. [50] BOP1 is a direct Wnt/β-catenin target gene and can induce epithelial-mesenchymal transition (EMT), cell migration and experimental metastasis of CRC cells. [51] DCAF13 is upregulated in breast cancer, lung cancer and hepatocellular carcinoma and indicates poor overall survival. [52][53][54] As mentioned above, DKC1, RRS1, UTP14A and BOP1 were related to CRC, which was reported in previous studies, but there has been no detailed study about the relationship between the NOP56, RRP9, BYSL and DCAF13 genes and CRC. Further research on the mechanism of the NOP56, RRP9, BYSL or DCAF13 genes in CRC may be helpful.
To investigate the functions of these hub genes and three key module genes, GO and KEGG analyses were performed. The results indicated that these genes possibly function by regulating RNA processes, RNA binding, ribosome complex biogenesis, IFN-alpha and IFN-beta biosynthetic processes and the Toll-like receptor signaling pathway. TLR3, TLR7 and TLR8 genes in the third module [ Figure 2d] were enriched in IFN-alpha and IFN-beta biosynthetic processes and the Toll-like receptor signaling pathway, and these three mRNAs were significantly decreased in CRC tissues compared with normal tissues. Toll-like receptors (TLRs) are widely expressed in tumor cells and may play a significant role in cancer biology. [55] Clinically, TLR agonists have been used as anticancer agents that activate immune cells in the tumor microenvironment. [56] IFN-α, -β, and -γ induce growth inhibition or cancer cell death. [57] Genetic variation of IFN is associated with the risk of colon cancer development. [58] Additionally, we constructed a multiple gene risk model predicting overall survival, including IGF2BP3, PABPC1L, PPARGC1A, PTRH1 and TDRD7. The model is an independent risk factor for clinicopathological characteristics to predict overall survival. IGF2BP3 regulates IGF2 translation and function and is expressed in testis tissues. [59] IGF2BP3 is correlated with several cancers and is significantly upregulated in CRC tissues. [59] PABPC1L is an important paralog of PABPC1 that regulates and stabilizes mRNA translation. [60] A study reported that PABPC1L is overexpressed in CRC and promotes CRC cell proliferation and migration by regulating the protein kinase B (AKT)-signaling pathway. [60] A recent study validated that PABPC1L might promote colon tumorigenesis by regulating mRNA splicing. [61] Mutation of PPARGC1A is related to the risk of CRC, and PPARGC1A protects against tumorigenesis by regulating the fate of the enterocyte cells. [62,63] PTRH1 is peptidyl-tRNA hydrolase 1 and salvages tRNA from peptidyl-tRNA, which is toxic to cells when accumulating. [64] TDRD7 is a scaffold protein and can induce congenital cataracts. [64] However, there is no definitive study about the relationship between PTRH1 or TDRD7 and cancer.
Finally, the expression of risk model genes was validated using online databases. IGF2BP3, PABPC1L and PTRH1 were significantly upregulated in CRC tissues, but PPARGC1A and TDRD7 were downregulated compared with the normal group. Our results were consistent with those in the online databases.
One limitation of the study is the lack of experimental data validating the variations in the risk model genes. Although, all of the model gene mRNA levels in CRC were validated by the TIMER database, only IGF2BP3, PTRH1 and TDRD7 protein levels were validated by the HPA database. Therefore, it would be more convincing if a related experiment was performed to investigate the protein expression levels of the model genes, especially PABPC1L and PPARGC1A. Another limitation is the lack of more clinical data validating the prognostic signature established in the study. GSE17538 was used as the validation cohort, but it is necessary to collect more data and samples in the clinic to validate the signature before clinical application.
In conclusion, we screened differentially expressed RBPs, identified hub genes and analyzed the functions and related pathways in CRC. Subsequently, we constructed a risk score model including the IGF2BP3, PABPC1L, PPARGC1A, PTRH1 and TDRD7 genes. The prognostic model may help predict the overall survival of patients with CRC and be applied in the clinic in the future. Contd...