Development and validation of a RNA binding protein-associated prognostic model for head and neck squamous cell carcinoma

Flow chart of our study

Figure 1 shows the detailed workflow of the study. Firstly, DERBPs that are differentially expressed between normal samples and HNSCC are found. Then, we used the training group to construct a specific prognostic model. The prognostic model was further confirmed and optimized in the testing group. The prediction power of these models was checked using time-dependent ROC analysis. Then, the RBPs in the prognosis model were comprehensively verified. The prognostic value and key role of RBPs was analyzed by our study systematically in HNSCC.

Figure 1. The Flowchart for identification the survival-related RBPs in HNSCC. RBPs: RNA binding proteins. HNSCC: Head and Neck Squamous Cell Carcinoma.

Differential expression and functional annotation analysis of RBPs in HNSCC

The corresponding clinical data and mRNA expression data of 44 nontumor samples and 502 HNSCC tissue samples were downloaded from TCGA database (Table 1). After abstracting 1495 RBP expression values, we acquired differentially expressed RBPs and showed the expression pattern of differentially expressed RBPs in HNSCC and non-tumor tissues by volcano map and thermogram. In HNSCC, 190 differentially expressed genes were obtained in tumour tissues, of which 109 genes were upregulated and 81 genes were downregulated (Figure 2). Then, to understand the biological characteristics of these genes, we performed functional enrichment analysis on upregulated and downregulated RBPs. Figure 3 summarizes the GO terms and enriched KEGG pathways of these genes. In HNSCC, our results indicate that in the related biological processes, upregulated RBPs are significantly enriched in nucleic acid phosphodiester bond hydrolysis, RNA catabolic process, negative regulation of viral genome replication, regulation of translation and DNA methylation or demethylation (Figure 3A), while downregulated RBPs are significantly enriched in the regulation of RNA splicing, regulation of the mRNA metabolic process, regulation of mRNA splicing via the spliceosome, regulation of mRNA processing, and regulation of cellular amide metabolic process (Figure 3B). In terms of molecular function, upregulated RBPs were significantly enriched in acting on RNA, catalytic activity, double-stranded RNA binding, helicase activity, ribonuclease activity and nuclease activity (Figure 3A), while downregulated RBPs were significantly enriched in mRNA binding, mRNA 3' UTR binding, poly(A) binding, single-stranded RNA binding and mRNA 3' UTR AU-rich region binding (Figure 3B). Through cell composition (CC) analysis, we found that the upregulated RBPs were mainly enriched in ribonucleoprotein granules, cytoplasmic ribonucleoprotein granules, P-bodies, spliceosomal complexes and P granules (Figure 3A), while the downregulated RBPs were significantly enriched in cytoplasmic stress granules, cytoplasmic ribonucleoprotein granules, ribonucleoprotein granules, RNA cap binding complexes and mRNA cap binding complexes (Figure 3B). In addition, the results of KEGG pathway enrichment analysis of differentially expressed RBPs showed that the upregulated RBPs were enriched in the spliceosome, RNA transport, the mRNA surveillance pathway, the RIG-I-like receptor signaling pathway and microRNAs in cancer (Figure 3C), while the downregulated RBPs were enriched in RNA transport, the mRNA surveillance pathway, RNA degradation and progesterone-mediated oocyte maturation (Figure 3D).

Figure 2. Differentially expressed RBPs in HNSCC and non-tumour samples. The volcano plot (A) and Clustered heatmap (B) of differentially expressed RBPs in HNSCC and normal tissues.

Figure 3. Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) analysis of the differentially expressed RBPs. Results of GO functional annotation analysis of up (A) and down (B) regulated differentially expressed RBPs; Results of KEGG pathways enrichment analyses of up (C) and down (D) regulated differentially expressed RBPs. The dot size represents the enriched gene number.

Construction of the PPI network and screening of key modules

To better understand the potential molecular functions of these differentially expressed RBPs in HNSCC, a PPI network was constructed using the Cytoscape software and STRING database. The PPI network consists of 161 nodes and 581 edges (Figure 4A). Then, we further analyzed the co-expression network and used the plug-in pattern in Cytoscape to detect the potential key modules and determine the first three important modules (Figure 4B). Module 1 consists of 16 nodes and 113 edges, module 2 consists of 10 nodes and 45 edges, and module 3 consists of 12 nodes and 35 edges. Functional enrichment showed that the genes of module 1 were mainly enriched in the regulation of multi-organism processes, RNA binding and nucleic acid binding. The genes of module 2 were significantly enriched in mRNA splicing, the spliceosome, and the mRNA surveillance pathway. However, the genes of module 3 were significantly enriched in rRNA processing, RNA processing, and translation (Supplementary Table 1).

Figure 4. Protein-protein interaction (PPI) network and modules analysis. (A) PPI network of differentially expressed RBPs; (B) 3 critical modules from PPI network. The outer circle is critical module 1, the middle circle is critical module 2, and the inner circle is critical module 3. Red circles: up regulation RBPs; Green circles: down regulation RBPs.

Construction and verification of the HNSCC-specific predictive prognosis model

A total of 161 key RBPs with differential expression were identified from the PPI network. In order to study the relationship between the expression of RBPs and the prognosis of HNSCC patients, we established the prognosis model of HNSCC patients in the training group. Univariate Cox regression analysis was initially executed to acquire genes that were significantly related to prognosis, and LASSO regression and multivariate Cox regression analyses were then used to create the final prognosis model (Table 2, and Figure 5A, 5B). After creating the prognosis model, the patients were divided into a low-risk group and a high-risk group, and Kaplan-Meier survival analysis was then carried out on the testing set and training set. The results demonstrated that the overall survival time of patients with high risk score was obviously poorer than those with low risk score in the HNSCC dataset (Figure 5C, 5D). In HNSCC patients, a 9-gene model (CELF2, CPEB1, DDX39B, EIF3L, EZH2, KHDRBS3, RNASE10, RNASE3 and SIDT1) was successfully obtained. Using this model to anticipate the risk score of each patient, we found that CPEB1, EIF3L, KHDRBS3 and RNASE3 are positive risk-related genes, whereas CELF2, DDX39B, EZH2, RNASE10 and SIDT1 are negative risk-related genes. The AUC (area under the ROC curve) in the training and testing groups was remarkable (3-year AUC, 0.752 vs. 0.669; 5-year AUC, 0.781 vs. 0.687) (Figure 5E, 5F). The model can accurately predict the OS of HNSCC patients. Besides, we ranked all HNSCC patients according to risk score to analyze the survival rate distribution. From the scatter plot, we can see the survival status of patients with different risk score; the mortality rate of patients rises with increasing risk score. Heat maps show that the expression of RBPs is related to the increase of patient risk score (Figure 6A–6F).

Figure 5. Construction and validation of the prognostic risk model in HNSCC patients. (A) Screening of optimal parameter (lambda) at which the dotted vertical lines were drawn. (B) Lasso coefficient profiles of the candidate RBPs with non-zero coefficients determined by the optimal lambda. Kaplan-Meier plot of the high-risk (red) and low-risk (blue) HNSCC patients in the training group (C) and testing group (D). The 3-year (red) and 5-year (blue) ROC curves in the training group (E) and testing group (F) of HNSCC patients.

Figure 6. RBPs-related prognostic characteristics in patients with HNSCC. Risk score distribution of HNSCC patients with different risks in the training group (A) and testing group (B) (low, green; high, red). Dot plots showing the survival time and risk score in training group (C) and testing group (D). The heatmap of the 9 key genes expression profiles in the training group (E) and testing group (F) (low, blue; high, red).

The prognosis model of HNSCC patients is independently associated with OS

Cox regression analysis was used to analyze the correlation between OS and clinical parameters such as histological grade, age, risk score and pathological stage. Univariate Cox regression analysis demonstrated that the histological grade, age, pathological stage, and risk score of patients were associated with OS (P <0.05). However, by multivariate regression analysis, we found that only risk score, N stage and age were independent prognostic factors associated with OS (P <0.05) (Figure 7A, 7B). At the same time, to construct a quantitative model of the prognosis of HNSCC patients, we also combined 9 RBP markers to construct a nomogram (Figure 7C and Supplementary Figure 3). Based on multivariate Cox analysis, the point scale of the nomogram was used to assign points to each variable. We drew a horizontal line to determine the points of each variable, calculated the total points of each patient by summing the points of all variables, and then standardized the total points of each patient to a distribution of 0 to 100. By constructing a vertical line between each prognostic axis and total point axis, the estimated 1-, 3-, and 5-year survival rates of HNSCC patients could be calculated, which may help relevant practitioners make clinical decisions for HNSCC patients. Our consequences suggest that the established specific prognostic models and genes can be used to predict the OS of HNSCC patients.

Figure 7. The prognostic value of different clinical parameters. Univariate (A) and multivariate (B) Cox regression analyses of OS in HNSCC. (C) Nomogram for predicting 1-, 3-, and 5-year OS of HNSCC patients.

Comprehensive analysis of genes in the RBP prognosis model

From the prognosis model, we acquired 9 genes and then further assessed the prognosis value of these genes in other databases. Correlation analysis of picked genes in TIMER (Tumor Immune Estimation Resource) database showed that most genes were closely associated with mRNA expression (Supplementary Table 2 and Supplementary Figure 1). The genes were analyzed by the GEPIA (Gene Expression Profiling Interactive Analysis) database using Kaplan-Meier analysis. The results show that in HNSCC, CELF2, EZH2, RNASE10 and SIDT1 are positively correlated with OS, indicating that these highly expressed genes correlate with a good prognosis (Figure 8). All in all, the Kaplan-Meier analysis results agree with the univariate Cox analysis results, which means that most genes are infused into specific prognosis models and have a potent predictive capacity. Next, the protein expression pattern of genes in the prognosis model were analyzed by the HPA database (Figures 9, 10A–10C). The results demonstrated that the expression of the CELF2 protein was low in normal head and neck tissues, but this protein was not detected in HNSCC tissues (Figure 10A). The EZH2 protein was low expressed in normal tissues and moderately expressed in tumor tissues (Figure 10B). The SIDT1 protein was moderately expressed in normal head and neck tissues but was not detected in HNSCC tissues (Figure 10C). We then used the cBioPortal database to find out the CNV (Copy number variations) and mRNA expression changes in these genes (Figure 10D). The results demonstrated that CNVs were associated with the mRNA expression changes in these genes. It is noteworthy that SIDT1 and KHDRBS3 demonstrated the highest mRNA expression changes and CNV in the entire analysis sample, which may suggest that CNVs are the leading driving forces for the mRNA expression changes in these genes. To verify our analysis, the relative mRNA expression of CELF2, EZH2 and SIDT1 in HNSCC cell lines was evaluated by qRT-PCR. As shown in Figure 10E, compared with that in the human bronchial epithelial cell line HBE, the mRNA expression of the CELF2 gene was relatively lower in the laryngeal squamous carcinoma cell line Hep2 (P<0.05) but showed no difference in that in the nasopharyngeal carcinoma cell line HK1 (P>0.05). The EZH2 mRNA expression was higher in both Hep2 and HK1 cells (P<0.01), and the SIDT1 mRNA expression was lower in both Hep2 and HK1 cells (P<0.01); these results are consistent with the protein expression patterns detected in HNSCC tissues in HPA database. Besides, EZH2 protein expression was higher in HK1 and Hep2 than HBE cells which was verified by western blot assay as showed in Figure 10F. In addition, the clone forming ability of 2 tumor cell lines was significantly higher than that of the control cell lines (Figure 10G, 10H). Our results agreed with findings we have observed previously, which further verified the reliability of our experiment. The GSEA method was used to calculate the pathways and enriched features between low-risk and high-risk patients, as we found out that low-risk and high-risk patients have significant prognostic differences in OS. In the results of GSEA enrichment, we observed that the high-risk groups were enriched in response of EIF2AK4 GCN2 to amino acid deficiency, ribosome, SRP-dependent cotranslational protein targeting to membrane, rRNA modification in the nucleus and cytosol, Myc targets v1 and allograft rejection. (Supplementary Figure 2A–2F). The low-risk group was enriched in the transcriptional regulation by RUNX1, interferon signaling, signaling by Wnt, MAPK family signaling cascades, hallmark interferon gamma response and hallmark interferon alpha response (Supplementary Figure 2G–2L). Some studies have shown that these pathways are related to the development of HNSCC. In short, the GSEA results indicate that RBP-related signals are associated with the progression and development of HNSCC.

Figure 8. Kaplan-Meier analyses of ARGs in prognostic model. Kaplan-Meier analyses of (A) CELF2, (B) EZH2, (C) RNASE10 and (D) SIDT1. The statistical significance was determined by Log-rank test.

Figure 9. Survival associated RBPs expression in the human protein atlas database (HPA). Kaplan-Meier curves of survival associated (A) CELF2, (B) EZH2, (C) RNASE10 and (D) SIDT1 for HNSCC patients. Pink line indicates high expression group while blue line indicates low expression group.

Figure 10. Immunohistochemistry (IHC), mutations and verification results in the prognosis-related RBPs. The protein levels of (A) CELF2, (B) EZH2 and (C) SIDT1 were determined by immunohistochemistry using indicated antibodies in HPA database, the staining strengths were annotated as Not detected, Low, Medium and High. The bar plots indicating the number of samples with different staining strength in HPA database. (D) OncoPrint showing the copy number alterations and mRNA expression alterations of 9 RBPs in prognostic model. (E) The mRNA expression levels of CELF2, EZH2 and SIDT1 in HNSCC cell lines were determined by qRT-PCR, GAPDH was verified as a housekeeping gene. (F) Western blot was conducted to evaluate the protein expression of EZH2 in HNSCC cell. (G, H) Colony-forming assay was conducted to evaluate the growth of HNSCC cell lines.

Data preprocessing and screening of differentially expressed RBPs

RNA sequencing data and corresponding clinical data of 502 HNSCC and 44 normal head and neck tissue specimens were acquired from the TCGA database (https://portal.gdc.cancer.gov/) (Supplementary Table 3). A total of 1542 RNA-binding proteins that have been confirmed to exist in the human genome were extracted from the literature. See Supplementary Table 4 for details. To ensure a unified standard, the RNA sequencing results (FPKM) of HNSCC were transformed and standardized by the formula log2 (x+1). The Wilcox test in R (version 4.0.0, https://www.r-project.org/) was used to calculate the differential expression of RBPs between HNSCC and control samples. The genes with at least 1.5-fold alteration and corresponding P value less than 0.05 were chosen as the differentially expressed RBPs (DERBPs).

Gene enrichment analysis and PPI network construction and module selection

The main biological characteristics of these RBPs were detected by GO enrichment and KEGG pathway analysis. Go analysis terms include biological process (BP), molecular function (MF) and cell composition (CC). Then differently expressed RBPs are submitted to the STRING database (http://www.STRING-db.org/) to identify protein-protein interaction information. The PPI network is further constructed and visualized by using the software of Cytoscape 3.6.1. The important modules and genes were screened out in PPI network by using molecular complex detection (mcode) plug-in. The MCODE scores and the number of nodes are both greater than 5. Statistically significant differences P<0.05 are denoted as starred values.

Construction and verification of prognosis model in HNSCC

We first integrated the expression data of DERBPs with the corresponding clinical information, selected RBPs co-expressed in the PPI network, and then randomly divided the data into a testing group and a training group for following verification. Univariate Cox regression analysis was used to analyze the expression data of RBPs in the training group, and RBPs that were significantly related to survival were obtained (P <0.05). Finally, we used lasso regression and multivariate Cox regression analysis to constructed a prognosis model. In the training group and the testing group, the risk score of each patient was calculated according to the regression coefficient and expression value of each gene in the model. The calculation formula is: β1 * Exp1 + β2 * exp2 + βi * EXPi, in which β represents coefficient value and exp represents gene expression level. The risk score is an indicator to measure the prognostic risk of each HNSCC patient. We used the median risk score to draw a Kaplan-Meier survival curve and generated receiver operating characteristic (ROC) curves to determine the accuracy of the prediction model. Finally, the rms R package was used to analyze the nomogram plot and predict the possibility of OS. P <0.05 is a significant difference.

Comprehensive analysis of RBPs in HNSCC prognosis model

Correlation analysis was performed on the selected RBPs in the TIMER database, and the Pearson correlation coefficient between each gene pair was calculated. Then use the GEPIA2 database (http://gepia.cancer-pku.cn/) to perform Kaplen-Meier analysis on the RBPs in the risk-specific model, and use the log-rank test to determine the statistical significance. The expression of RBPs at translation level was analyzed by comparing the immunohistochemical staging images in HPA database (http://www.protein.atlas.org/). On the basis of the staining intensity, it is marked as high, medium, low and undetected. The cbioProtal database (http://www.cbioportal.org/) was used to further analyze the nine RBPs in the risk-specific model to evaluate alterations in mRNA expression and copy number. Last, GSEA (gene set enrichment analysis) was executed to examine the ways and characteristics of enrichment in the predicted low-risk and high-risk population. Using GSEA, this research studied whether the characteristics of activation/inhibition genes were abundant in low-risk and high-risk patients. The standardized enrichment score (NES) and standardized P value are used to calculate the enrichment of the Hallmarks and canonical pathways. Terms with | NES |> 1 and P <0.05 are considered significantly rich.

Verification of screened RBPs by qRT-PCR

The human bronchial epithelial cell line HBE, laryngeal squamous carcinoma cell Hep2 and nasopharyngeal carcinoma cell line HK1 were obtained from the Shanghai Zhong Qiao Xin Zhou Biotechnology. Hep2 and HK1 cells were maintained at 37° C in a 5% CO₂ incubator in Dulbecco’s modified Eagle’s medium (DMEM) with 10% fetal bovine serum (FBS) (Biological Industries, Kibbutz Beit Haemek, Israel) and 1% penicillin–streptomycin (Keygen Biotech, Nanjing, China). HBE cells were cultured with specific keratinocyte medium supplemented with 1% keratinocyte growth factor and 1% penicillin–streptomycin (Zhong Qiao Xin Zhou Biotech, Shanghai, China). For quantitative analyses of screened RBPs, total RNA was extracted using TrizoL (Life Technologies, Carlsbad), and cDNA was prepared using High Capacity cDNA Reverse Transcription Kit (Thermo Fisher, MA) and qPCR were performed with NovoStart® SYBR qPCR SuperMix Plus (Novoprotein, Shanghai, China). The qPCR primer sequences of EZH2 gene were 5’-CCCGCTGAGGATGTGGATAC-3’ and 5’-CATGGTTAGAGGAGCCGTCC-3’, primer sequences of SIDT1 were 5’-AGCCCCTCTCAACCTCAGTA-3’ and 5’-GCAGCTTTCTTGGTCATGGA-3’, primer sequences of CELF2 were 5’-CAGCACCAATGCAAACCCTC and TCCCGAGAGAGGTCAAGGAG-3’, which were designed by Primer Premier 5.

The quality of cDNA samples was verified using GAPDH as a housekeeping gene.

Verification of screened RBPs by western blot

Protein lysates were prepared using RIPA lysis buffer (0.1% SDS, 1% NP-40, 1 mM EDTA, 50 mM Tris PH 7.5, 150 mM NaCl, 0.25% deoxycholate) with protease and phosphatase inhibitors (Roche, Welwyn Garden City, UK). Protein concentration was determined using BCA Protein Assay Kit (Beyotime, Shanghai, China). Following 10% SDS gel electrophoresis and subsequent immunoblotting, bound anti-EZH2 antibody (1:1000), protein expression was detected by ECL Chemiluminescence substrate (Biosharp, Shanghai, China).

Verification of screened RBPs by colony-forming assay

HBE, HK1 and Hep2 cells were diluted and plated into six-well plates in triplicate to execute a colony-formation assay (CFA). Inoculate 200 cells into each well of a 6-well plate, change the medium every 3 days, and culture for 2 weeks. After that, the colonies containing more than 100 cells were fixed by methanol, stained using 0.3% crystal violet (Sigma, USA) and counted manually. And the colony formation rate was reckoned using formula as below: colony formation rate = (number of colonies / number of seeded cells) × 100%.