Development and validation of a prediction model for lung adenocarcinoma based on RNA-binding protein.

BACKGROUND: RNA-binding proteins (RBPs) have been found to participate in the development and progression of cancer. This present study aimed to construct a RBP-based prognostic prediction model for lung adenocarcinoma (LUAD).
METHODS: RNA sequencing data and corresponding clinical information were acquired from The Cancer Genome Atlas (TCGA) and served as a training set. The prediction model was validated using the dataset in Gene Expression Omnibus (GEO) databases. Univariate and multivariate Cox regression analyses were conducted to identify the RBPs associated with survival. R software (http://www.r-project.org) was used for analysis in this study.
RESULTS: Nine hub prognostic RBPs (CIRBP, DARS2, DDX24, GAPDH, LARP6, SNRPE, WDR3, ZC3H12C, ZC3H12D) were identified by univariate Cox regression analysis and multivariate Cox regression analysis. Using a risk score based on the nine-hub RBP model, we separated the LUAD patients into a low-risk group and a high-risk group. The outcomes revealed that patients in the high-risk group had poorer survival than those in the low-risk group. This signature was validated in the GEO database. Further study revealed that the risk score can be an independent prognostic biomarker for LUAD. A nomogram based on the nine hub RBPs was built to quantitatively predict the prognosis of LUAD patients.
CONCLUSIONS: Our nine-gene signature model could be used as a marker to predict the prognosis of LUAD and has potential for use in treatment individualization. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Lung cancer has the highest mortality among cancers and the second highest incidence in both men and women (1). There were 228,820 new cases of lung cancer diagnosed in the United States and more than 135,720 related deaths in 2020 (1). An estimated 85% of lung cancers are non-small cell lung cancer (NSCLC), which mainly consists of lung squamous cell carcinoma and lung adenocarcinoma (LUAD) (2). Of these, LUAD is the most prevalent NSCLC subtype, constituting about 40% of all lung cancer cases (3). Despite significant development in the diagnosis and treatment of lung cancer, the prognosis of lung cancer remains unfavorable mostly due to a lack of early diagnoses (4). Therefore, in order to better improve the prognosis of lung cancer patients, new prognostic indicators for lung cancer urgently need to be identified. RNA-binding proteins (RBPs) are proteins that can bind and impact a variety of RNAs (5). RBP dysregulation can lead to alteration in RNA metabolism, resulting in a change in the transcriptome and proteome of the cells, further affecting cell growth, proliferation, invasion, and death (6). RBPs can bind and influence RNA to improve the stability of target messenger RNA (mRNAs) and thus promote the expression of genes. Through this mechanism, RBPs perform a crucial role in the development of various diseases (7,8). Thus far, nearly 1,500 RBPs have been reported in the human genome (9). Accumulating evidence indicates that RBPs play an essential role in processes such as preserving the physiological balance of cells, cell development, and stress response (10). Mechanisms of RBPs regulation have been explored in cancer cells, including stability, translation, alternative splicing, subcellular localization, and polyadenylation (11). It has been demonstrated that RBP dysregulation contributes to cancer progression by altering the expression of oncogenes or tumor suppressor genes (12). LIN-28 homolog B (LIN28B) promotes versatility through interaction with the let-7 family microRNA, and is pivotal in the pathogenesis of colorectal cancer (13,14). Meanwhile, polypyrimidine tract-binding protein 3 (PTBP3) mediates the variable shear of caveolin 1 (CAV1), thus affecting the migration and invasion of gastric cancer (15). In lung cancer, the downregulation of cell proliferation splicing factor Quaking (QKI) is associated with poor survival rate (16). RBPs play a regulatory role in many steps of gene expression, Li et al. have established a predictive model in lung squamous cell carcinoma (LUSC) (17), but few RBPs have been systematically investigated in LUAD thus far.

The purpose of this study was therefore to construct an RBP-based prognosis prediction model for LUAD. We downloaded RNA sequencing data and clinical information for LUAD from The Cancer Genome Atlas, and differentially expressed RNA-binding protein (DERBP) genes were screened. For the purpose of recognizing key RBPs in LUAD that could serve as prognostic markers, we performed a functional analysis of the downloaded data. We present the following article in accordance with the TRIPOD reporting checklist (available at http://dx.doi.org/10.21037/atm-21-452).

Methods

Screening of differentially expressed RBPs in LUAD and data processing

We obtained the RNA sequencing data of 497 LUAD tumor tissue samples and 54 normal lung tissue samples, along with matching clinical information from TCGA. The P values of differentially expressed genes (DEGs) between the LUAD sample of TCGA was analyzed with Wilcoxon test using the “limma” R package (18). To identify RBPs that were differentially expressed in LUAD, the cutoff threshold in TCGA was set as |log2 fold change (FC)| ≥1.0, and the false discovery rate (FDR) was <0.05. All samples were used as the training cohort. In addition, the external validation cohort consisted of the expression data and comparative clinical data that were acquired from the Gene Expression Omnibus (GEO) database (GSE13213, n=117). This research was conducted in line with the Declaration of Helsinki (as revised in 2013).

Kyoto Encyclopedia of Genes and Genomes pathway and Gene Ontology enrichment analyses

In order to identify the biological function of differentially expressed RBPs, R software was used for Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses. The GO analysis terms included biological process (BP), cellular component (CC), and molecular function (MF). The FDR (q value) <0.05 was regarded as the statistical criterion.

Protein-protein interaction network structuring and key module screening

We submitted 349 RBPs to the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (http://www.string-db.org/). A protein-protein interaction (PPI) network was framed with Cytoscape v3.7.1 software. and its subnetworks were visualized. The MCODE tool in Cytoscape was used to manage the PPI networks, recognize the key networks, and construct a network of 349 nodes and 3,782 edges.

Building a prognostic model

The “survival” package from R software was used to conduct univariate Cox regression analysis for all hub RBPs in the key modules of the training data set. We further screened for significant candidate genes using the log-rank test. Subsequently, based on the initially screened significant candidate genes noted above, a multivariate Cox regression model was built to calculate the risk score to evaluate the prognosis of LUAD patients. The risk score of each LUAD specimen was determined by the following formula:

LUAD patients were separated into high-risk and low-risk groups according to the median value of the risk score. Subsequently, a log-rank test was performed to compare the difference of OS between the two subgroups. Furthermore, we used the “survivalROC” package to assess the performance of the prognostic model and drew a receiver operating characteristic (ROC) curve (19). Then, to confirm the predictive power of this prognostic model, 117 LUAD patient samples with credible prognostic data were downloaded from the GSE13213 dataset as a validation cohort.

Nomogram construction and validation of expression level

We used R software to construct a nomogram consisting of nine hub RBPs with the aim of providing a basis for clinicians to evaluate the prognosis of lung cancer patients. The expression level of 9 key RBPs was verified in the Human Protein Atlas (HPA) database (http://www.proteinatlas.org/) (20).

Statistical analysis

Survival curves were constructed by Kaplan-Meier analysis, and the significance of differences was evaluated by log-rank test. Wilcoxon signed-rank test were used to probe into quantitative variables. Significance was defined as P<0.05. ALL statistical analyses were performed using the R software v4.0.3.

Results

Selection and identification of differentially express RBPs in LUAD

The study process is diagrammatized in . We acquired LUAD data from TCGA database, which included 54 normal lung tissue and 497 tumor tissue samples. The data were processed to find the differently expressed RBP (DERBP) genes by using R software. The RNA sequences of 1,542 RBPs (9) were ultimately included in the analysis; a total of 368 RBPs including 239 upregulated and 129 downregulated genes were deemed to be differentially expressed between normal tissues and those with tumor (P<0.05, |log2 fold change| >0.5, Table S1). A heat map and volcano plot of the results of the analysis are shown in .

Figure 1

Flowchart of lung adenocarcinoma (LUAD) RNA binding proteins (RBPs) analysis.

Figure 2

Differentially expressed RBPs (DERBPs) in LUAD. (A) Heat map of DERBPs between tumor and normal tissues; (B) Volcano plot of DERBPs: upregulated DERBPs are indicated by red dots, and downregulated DERBPs are indicated by green dots. LUAD, lung adenocarcinoma.

KEGG pathway and GO enrichment analysis

The enrichment of the identified RBP-encoding genes was evaluated by R software. The KEGG enrichment analysis indicated that the upregulated DERBP genes expression was dramatically enriched in three biological processes: ribosome, RNA transport, and ribosome biogenesis in eukaryotes (). Meanwhile, the downregulated DERBP genes expression was significantly enriched in mRNA surveillance pathway, RNA transport, and spliceosome (). GO analysis showed that of the upregulated DERBP genes enriched in BP, the majority included non-coding RNA (ncRNA) metabolic processing, ncRNA processing, ribonucleoprotein complex biogenesis, transfer RNA (tRNA) metabolic process, and ribosome biogenesis. Meanwhile, the downregulated DERBP genes were mainly involved in the regulation of translation, regulation of cellular amide metabolic process, regulation of mRNA metabolic process, RNA splicing, and RNA catabolic process. Under CC, we detected that upregulated DERBP genes were observably enriched in preribosome, ribosome, mitochondrial matrix, and ribosomal subunit. Meanwhile, the downregulated DERBP genes were significantly enriched in five biological processes: cytoplasmic ribonucleoprotein granule, ribonucleoprotein granule, nuclear speck spliceosomal complex, and P-body. In the MF analysis, the upregulated DERBP genes were notably enriched in catalytic activity, acting on RNA, catalytic activity, acting on a tRNA, ribonuclease activity, structural constituent of ribosome and nuclease activity, whereas the upregulated DERBP genes were significantly enriched in mRNA 3'-untranslated region (3'UTR) binding, translation regulator activity, translation regulator activity, nucleic acid binding, and single-stranded RNA binding ().

Figure 3

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology (GO) enrichment analysis of DERBPs. In KEGG analysis, the abundance of RBP is expressed in terms of dot size. Low q value is represented in red, and high q value is represented in blue. A darker red color displaying a smaller q value, in the same way, darker blue color indicating a higher q value. And the q value of GO terms are illustrated by color bar. BP, biological process; CC, cellular component; MF, molecular function. (A) KEGG pathway analysis of upregulated DERBPs. (B) KEGG pathway analysis of downregulated DERBPs. (C) GO terms of upregulated DERBPs. (D) GO terms of downregulated DERBPs.

PPI network structure and selection of key modules

To further investigate the mutual effect of DERBPs and identify the key RBPs in LUAD, we constructed a PPI network using Cytoscape software. This network, which is based on the data came from the STRING database, contained 349 nodes and 3,782 edges (). Next, the MODE tool in Cytoscape was used to manage the PPI networks and detect potential critical modules. From this, four modules were selected, the most important of which was composed of 44 nodes and 861 edges ().

Figure 4

Protein-protein interaction (PPI) network and modules analysis. (A) PPI network of DERBPs; (B) four critical modules from PPI network. Green circles: downregulation. Red circles: upregulation.

Selection of prognosis-related RBPs

A total of 349 DERBPs were recognized from the PPI network. In order to further explore the prognostic value of these RBP-related genes, we conducted a univariate Cox regression analysis and acquired 17 candidate hub DERBP genes associated with prognosis (). Subsequently, a multiple stepwise Cox regression analysis based on these 17 prognostic-associated candidate hub RBPs was performed to determine their influence on the overall survival (OS) of LUAD patients. This yielded nine hub DERBP genes that were independent predictors in LUAD patients (, ).

Figure 5

Forrest plot of univariate Cox regression analyses to identify hub RBPs in LUAD. RBPs, RNA binding proteins; LUAD, lung adenocarcinoma.

Figure 6

Forrest plot of multivariate Cox regression analysis for identification of prognosis related hub RBPs.

Table 1

Nine prognosis-related hub DERBPs identified by multivariate Cox regression analysis

RBP	Coef	HR	Lower 95% CI	Upper 95% CI	P value
CIRBP	0.239382	1.270463	0.931909	1.732012	0.130036
DARS2	0.245832	1.278685	0.948215	1.724329	0.10709
DDX24	−0.38377	0.681291	0.442956	1.047863	0.080617
GAPDH	0.301149	1.35141	1.080187	1.690734	0.008418
LARP6	0.212729	1.23705	0.956383	1.600082	0.105171
SNRPE	0.372009	1.450647	1.001597	2.101021	0.049021
WDR3	0.319847	1.376916	0.935578	2.026446	0.104755
ZC3H12C	0.416629	1.51684	1.111659	2.069702	0.0086
ZC3H12D	-0.41679	0.659161	0.401293	1.082733	0.099756

(DE)RBPs, differentially expressed RNA-binding proteins; Coef, coefficient; HR, hazard ratio; CI, confidence interval

Construction and analysis of an RBP-based prognosis

A prediction model was established based on nine key DERBP genes recognized by multivariate Cox regression analysis. The LUAD patient risk score was calculated by the following formula:

Next, we performed a survival analysis to evaluate the predictive power of the model. To account for the risk scores, the 497 LUAD patients were separated into high-risk and low-risk subgroups. The results confirmed that the survival curve of the high-risk subgroup was lower than that of the low-risk subgroup (). We further performed received operating curve (ROC) analysis to assess the predictive power of the OS-related prognostic models based on the nine RBPs and drew an ROC curve. The results showed that the area under the ROC curve (AUC) of this OS-related prognostic model was 0.707 for the training cohorts (), which demonstrated that it had certain predictive ability. The heat map of RBP expression, the status of patient survival, and the distribution of risk scores are summarized in . Next, we used the same formula to verify the predictive value of the OS-related model in the GSE13213 dataset and evaluated whether the model had similar predictive power in other datasets. We discovered that, in the GSE13213 cohorts, the test group and training cohorts had similar results ). This indicated that the nine DERBP genes-based model had good specificity and sensitivity in predicting the prognosis of patients with LUAD.

Figure 7

The nine RBP-related gene prediction model associated with overall survival (OS) of LUAD in TCGA dataset. (A) Survival curve for low-risk and high-risk groups in the prediction model; (B) receiver operating characteristic (ROC) analysis of the prediction model; (C) Heat map of 9 RBP gene expression in TCGA dataset; (D) risk score distribution; (E) survival status. The dotted line indicates the median risk score; most of the patients on the right died, showing an increased risk of death as the risk score increased. LUAD, lung adenocarcinoma.

Figure 8

The nine RBP-related gene prediction model associated with OS of LUAD in GSE13213 dataset. (A) Survival curve for low-risk and high-risk groups in the prediction model; (B) ROC analysis of the prediction model; (C) Heat map of 9 RBP gene expression in GSE13213 dataset; (D) risk score distribution; (E) survival status. The dotted line indicates the median risk score; most of the patients on the right died, showing an increased risk of death as the risk score increased. OS, overall survival; LUAD, lung adenocarcinoma.

Building a nine hub DERBP genes-based nomogram

To quantitatively evaluate the effect of the nine hub DERBP genes-based gene signature in the prognostic model, we constructed a nomogram to the predict survival time of LUAD patients (). From the results of multivariate Cox analysis, we built a nomogram that would be capable of predicting the 1-, 3-, and 5-year OS of LUAD patients using the R package “rms”.

Figure 9

A nomogram was built on the basis of the nine RBP genes for predicting OS of LUAD patients. RBPs, RNA binding proteins; OS, overall survival; LUAD, lung adenocarcinoma.

We used Cox regression analysis to evaluate the prognostic consequence of the different clinical characteristics of LUAD patients. The results of univariate Cox regression analysis revealed that in both the training cohort and test cohort, the risk score was a significant risk factor for OS and that tumor stage was correlated with the OS of LUAD patients (P<0.01). In contrast, sex and age were irrelevant to OS (). In short, multivariate Cox regression analysis revealed that risk score and tumor stage were independent risk factors for OS in LUAD patients (P<0.05; ).

Figure 10

Univariate analysis was performed to assess the clinicopathological prognostic value of the prediction model in LUAD. (A) Training cohort; (B) test cohort. LUAD, lung adenocarcinoma.

Figure 11

Multivariate analysis was performed to assess the clinicopathological prognostic value of the prediction model in LUAD. (A) Training cohort; (B) test cohort. LUAD, lung adenocarcinoma.

Validation the expression of key RBPs

To further investigate the expression of the nine risk-related genes in LUAD patients, we surveyed the protein expression of these DERBPs in the HPA database: DDX24, GAPDH, LARP6, and ZC3H12D were overexpressed in LUAD tumor tissue compared to normal tissue. By contrast, the expression of WDR3 in LUAD tumor tissue was relatively low (). Furthermore, the immunoreactivity of DARS2 and ZC3H12C was not significantly different between the LUAD tumor tissue and normal tissue. Information on CIRBP and SNRPE expression levels is not available in the HPA.

Figure 12

Validation of representative RBP expression of LUAD tumor and normal tissue in the Human Protein Atlas (HPA) database. RBP, RNA binding protein; LUAD, lung adenocarcinoma.

Discussion

An increasing amount of research has confirmed that the dysregulation of RBPs plays a critical role in cancer emergence, development, and progression (6). Nevertheless, merely a fraction of RBPs that contribute to tumorigenesis and progression have been extensively studied and identified (21). In the present study, we identified 368 DERBPs between normal and LUAD tissue of LUAD based on data from TCGA. Through KEGG and GO enrichment analysis, relevant biological pathways were systematically analyzed, and co-expression networks and PPI networks of DERBPs were constructed. In addition, univariate Cox regression and multivariate Cox regression were performed to build a risk prediction model. Survival analysis was implemented to further investigate the model’s clinical significance. We established a nomogram based on the nine hub DERBP genes to predict LUAD prognosis. These discoveries may help develop new biomarkers for the diagnosis and prognosis of LUAD patients.

The GO and KEGG pathway analysis demonstrated that the DERBPs were notably enriched in RNA transport, ribosome, ncRNA metabolic process, preribosome, and translation regulator activity. Furthermore, a DERBP protein interaction network was created, and we obtained modules containing 349 key DERBPs. Some of these hub DERBP genes have been revealed to serve a critical role in the progression of tumors (22,23). For instance, eukaryotic translation initiation factors (EIFs) have been used as a potential therapeutic target for patients with a variety of tumors. In previous studies, lung cancer patients with high EIF6 expression were predicted to have poor prognosis (24,25). Moreover, the absence of BOP1 was discovered to lead to the activation of MAPK pathway, leading to BRAF inhibitor resistance in melanoma (26). GNL3, as an oncogene, promotes the progression of osteosarcoma by regulating epithelial-mesenchymal transition (EMT) (27). BYSL was found to enhance glioblastoma cell invasion, cell migration, and EMT through the GSK-3β/β-catenin signaling pathway (28). Furthermore, Larp6-mediated mRNA localization was shown to be a key regulator of ribosomal biogenesis during cell migration, while EMT was demonstrated to upregulate the expression of LARP6 and promote the invasion and metastasis of tumor cells (29). The module analysis of PPI network in our study indicated that ribosome biogenesis in eukaryotes, ncRNA metabolic process, catalytic activity, acting on RNA, and preribosome were associated with LUAD.

Univariate and multivariate Cox regression analysis identified nine DERBP genes which showed prognostic value in LUAD, including CIRBP, DARS2, DDX24, GAPDH, LARP6, SNRPE, WDR3, ZC3H12C, and ZC3H12D. Previous studies have reported that the expression of SNRPE (30,31) and GAPDH (32) were relevant to tumor occurrence and development of lung cancer, which is consistent with our finding. Multivariate Cox regression analysis was then used to establish a risk model based on nine hub DERBPs to predict LUAD prognosis. This model used a TCGA cohort as a training set and a GEO cohort as a testing set. Subsequently, the ROC curves showed that the nine-gene signature had good diagnostic ability and could be used as a biomarker for screening LUAD patients with poor prognosis. However, the mechanisms by which these nine DERBP genes interact with tumorigenesis have not been clarified, and thus further studies are needed to fully elucidate the underlying mechanisms. Subsequently, we further built nomograms to estimate 1-, 3-, and 5-year OS. Our findings suggested that the nine DERBP genes-based gene model may be applied clinically for the diagnosis, prognostic management, and treatment of patients with LUAD.

However, there were some limitations in this study. First, the generation of prediction models depended on retrospective data from TCGA, which were not verified using clinical samples or prospective clinical study. Moreover, the clinical information acquired from the TCGA was insufficient, and may produce discrepancies in the results. Even so, our nine DERBP genes-related gene model shows immense potential for predicting prognosis in patients with LUAD. It may provide a useful reference for clinicians in determining customized clinical treatment options and providing personalized prognoses.

Conclusions

Bioinformatics analysis was performed to identify the potential prognostic key RPBs correlated with LUAD. The prognostic model we obtained was found to be an independent prognostic factor for LUAD. These nine DERBPs might be involved in LUAD tumor formation, development, and prognosis. Therefore, the nine-gene signature can be used as a marker to predict the postoperative prognosis and tumor progression of LUAD and may serve as novel a strategy for targeted intervention in LUAD.

The article’s supplementary files as