Development and Validation of an RNA Binding Protein-associated Prognostic Model for Hepatocellular Carcinoma

Clinical data and RNA sequencing data of the patients

We downloaded the RNA high-throughput sequencing data from The Cancer Genome Atlas (cancergenome.nih.gov/; TCGA) database, including 374 tumor tissue samples and 50 normal liver tissue samples (www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga ). We also downloaded the clinical data of 377 HCC patients from the TCGA. After excluding six patients with incomplete data, a total of three hundred and seventy-one patients were finally included in the study, with each having data on follow-up time, survival status, disease stage, etc. The RNA high-throughput sequencing data included the expression data of 60,483 genes, and we extracted the expression levels of 1,473 RBPs from such. As this research did not involve human participants, no research ethics review was necessary.

Identification of the differentially-expressed RBPs (DERs) in HCC patients

We used R software to analyze expression of the extracted 1,473 RBPs. A total of 325 RBPs showed expression differences between the tumor tissues and normal tissues, including 203 up-regulated and 122 down-regulated DERs. The threshold for the DERs was set as |log fold-change (FC)|>1 and false discovery rate (FDR) <0.05.

Enrichment analysis and protein-protein interaction (PPI) network of the DERs

We divided the DERs into an up-regulated group and a down-regulated group, and then performed gene ontology (GO) function analysis on the two groups of RBPs for three GO domains: molecular functions (MF), biological processes (BP) and cellular components (CC). Next, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis on the two groups of patients. We used the clusterProfiler package in R software for the GO function enrichment analysis and the KEGG pathway enrichment analysis. Then, we uploaded these DERs to the online tool STRING (https://string-db.org ) to construct the PPI network. We deleted the disconnected nodes, and the remaining RBPs were used for the next analysis. To further study the role of the DERs in HCC, we used Cytoscape software to create a PPI network that incorporated the nodes from the STRING database. At the same time, we also used the MCODE tool in the Cytoscape software to create a PPI subnet. Then, enrichment analysis and coexpression analysis were performed on the sub-network.

Screening of prognostic RBPs, construction and testing of the prognostic models

We used univariate Cox regression analysis to screen out 29 RBPs related to the patients’ outcomes, from among all of the DERs. The threshold of the Cox regression analysis was set to 0.001. All HCC patients were divided into a training group and a test group. Then, in the training group, a multivariate Cox regression analysis was performed to screen out nine RBPs. A prognostic model was constructed based on the relationship between the expression levels of these nine RBPs and the patients’ outcomes in the training group. The risk score of the HCC patients in the test group was calculated according to the prognostic model, and the patients in the test group were divided into high-risk and low-risk groups according to the risk scores. Then, we used survival curves, receiver operating characteristic (ROC) curves, risk curves and independent prognostic analysis to test the predictive power of the prognostic model. In addition, we still used the data GSE76427 in the Gene Expression Omnibus (GEO) database to verify our model externally.

Construction of the nomogram

The coefficients obtained by the Cox regression model were used to construct the nomogram of overall survival. To construct the nomogram, we first determined a scale axis of 0–100 to represent the score (in order to make the expression of the RBPs correspond to the scores on the scale axis) and we calculated the score of each RBP. We added the scores of each RBP to obtain the total score. Based on the correspondence between the total score axis and the survival rate, the 1-, 2-, 3- and 5-year survival rates of the patients could be predicted.

Cell culture, RNA isolation and quantitative real-time PCR (qRT-PCR) analysis

Because the risk ratio of enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2) is greater and the coefficient in the risk score calculation formula is greater, we choose EZH2 for further verification. A total of 20 HCC patient samples and paired non-tumor liver tissue samples were collected from the patients hospitalized at the Zhongnan Hospital, Wuhan University (Hubei Province, China). All patients provided written informed consent to the use of tissues for scientific research in the Department of Hepatobiliary and Pancreatic Surgery. All cell lines were purchased from the Chinese Type Culture Collection. All the cell lines were maintained in Dulbecco’s modified Eagle’s medium/high glucose (GE, USA) containing 10% fetal bovine serum (Gibco, USA). RNA was extracted using TRIzol reagent (TaKaRa, Japan) according to the manufacturer instructions. cDNAs were generated by the reverse transcription synthesis kit (TaKaRa) and the SYBR Green PCR Kit (TaKaRa) was used for qRT-PCR analysis. The primers used for EZH2 were 5′-GAGTTGGTGAATGCCCTTGGT-3′ and 5′-CATCTCGGTGATCCTCCAGATC-3′.

Gene knock-down

For generation of EZH2 knock-down cells, small interfering (si)RNA transfections were performed using the siRNA transfection reagent. The following siRNA sequences were used: EZH2-siRNA, CGGCUUCCCAAUAACAGUATT, UACUGUUAUUGGGAAGCCGTT.

Cell proliferation, migration, invasion, apoptosis, and cell cycle analyses

A CCK8 Kit (Dojindo, China) was used to measure cell viability. Transwell assay was used to assess cell migration. Cell migration was assessed by in vitro scratch wound assay. A BioCoat Matrigel Invasion Chamber (BD Biosciences, USA) was used to assess cell invasion; the number of cells migrating and invading was counted in three random areas. Apoptosis was assayed using the Annexin V-FITC Apoptosis Detection Kit (Invitrogen, USA), according to the manufacturer’s instructions, and the percentage of apoptotic cells was verified by flow cytometry (Beckman-Coulter, USA). To detect the cell cycle, 48 hours after transfection, the cells were stained with PI (propidium iodide) and assessed.

Statistical analysis

The datasets generated and analyzed during the current study are available in TCGA (https://cancergenome.nih.gov/ ) and GEO (www.ncbi.nlm.nih.gov/geo ). Difference analysis and regression analysis of the data were processed in the R3.6.3 software (https://www.r-project.org ). Statistical analysis of clinical data was performed using SPSS v.23.0 (IBM Corp., USA). Variables were compared by t-test or chi-square test.

DERs in HCC patients

A flowchart showing the steps of data processing and integrative data analysis is provided in Figure 1A. The original data downloaded from TCGA included the RNA sequencing data of 50 normal tissues and 374 liver cancer tissue samples. The sequencing data itself comprised 60,483 sets of RNA expression data, including that of 1,473 RBPs. After data processing, it was found that 325 RBPs had expression differences (Fig. 1B), including 203 up-regulated and 122 down-regulated RBPs (Fig. 1C). The results of the DERs are presented in the form of heat maps and volcano maps.

GO functional and KEGG pathway enrichment analysis

The GO function enrichment analysis was performed on the obtained up-regulated RBPs and down-regulated RBPs. The GO terms included BP, MF, and CC. For BP, the up-regulated RBPs were mainly enriched in non-coding RNA metabolic process, RNA splicing, and non-coding RNA processing. For CC, the up-regulated RBPs were mainly enriched in spliceosomal complex, cytoplasmic ribonucleoprotein granules, and ribonucleoprotein granules. For MF, the RBPs were mainly enriched in catalytic activity and acting on RNA (Fig. 1D). Down-regulated RBPs were mainly enriched in regulation of translation, regulation of cellular amide metabolic processes, cytoplasmic ribonucleoprotein granules, ribonucleoprotein granules, and catalytic activity (Fig. 1E). KEGG pathway enrichment analysis of up-regulated RBPs showed that they were mainly enriched in spliceosome, RNA transport, and mRNA surveillance pathways (Fig. 1F). Down-regulated RBPs were enriched in herpes simplex virus 1 infection, influenza A, and RNA degradation (Fig. 1G). The results of the enrichment analysis are also displayed in the form of tables (Table 1 and Table 2). The results in the chart are sorted by p-values, and the top 15 results of the GO enrichment analysis are shown.

PPI network and the coexpression network of RBPs

To study the interactions among the 325 RBPs, we used STRING to construct a PPI network. The connections between the molecules represent the possible interactions between two protein molecules, and the different colors of the lines represent different levels of evidence. After deleting 21 disconnected RBPs, there were a total of 304 nodes and 2,794 connections in the PPI diagram (Fig. 2A). Then, we used Cytoscape software to create a coexpression network of interactions among all nodes in the PPI. Among them, the block of proliferation 1 (i.e. BOP1) had the most interactive RBPs, and it had confirmed or potential interactions with 50 RBPs (Fig. 2B). We applied the MCODE tool in Cytoscape to construct a sub-network of the coexpression network (Fig. 2D–H), and then we selected the first five sub-networks according to their association score to identify the first important module, including 115 nodes and 1,295 edges (Fig. 2C).

Screening of prognosis-related RBPs and construction of a risk scoring model

A total of 305 RBPs were included in the PPI network. To screen-out the RBPs related to prognosis, we used univariate Cox regression analysis to screen-out a total of 29 RBPs, and we calculated their hazard ratios (Fig. 3A). Among them, 19 RBPs had a hazard ratio greater than 1, which means they had a negative impact on patient prognosis, and 10 RBPs had a positive impact on the prognosis. We then divided the patients into a training group and a test group at a ratio of 7:3, and we performed multivariate Cox regression analysis on the expression of the selected RBPs of patients in the training group (Fig. 3B). Then, we built a prediction model based on the relationship between the expression of RBPs in the training group of patients and the patient’s survival and survival status. The formula for calculating the patient risk score was as follows:

To test the effectiveness of the predictive model, we used the survival curve method to evaluate the prognostic model. In the training group, patients with high-risk scores had a worse prognosis than patients with low-risk scores (p<0.01; Fig. 3C). In the test group, patients were divided into a high-scoring group and a low-scoring group according to the risk score model, and the survival rates of the two groups were also significantly different (p<0.01; Fig. 3D). In addition, we adopted the ROC test method and calculated the area under the ROC curve (AUC). The AUC of the training group was 0.735 and the AUC of the test group was 0.740, indicating that the risk scoring model we constructed has good diagnostic performance (Fig. 3F, G). Patients in the GEO cohort were also divided into high-risk groups and low-risk groups based on this model. There were also significant differences in survival rates between the two groups (p<0.05), and the AUC of the GEO cohort was 0.740 (Fig. 3E, H). This shows that the risk scoring formula we established can accurately divide patients into a high-risk group with a poor prognosis and a low-risk group with a good prognosis.

Independent prognostic analysis of the risk score

To evaluate whether the risk score is an independent prognosis-related factor, we conducted an independent prognostic analysis of risk scores for patients in the training group and test group. The prognostic analysis used univariate and multivariate Cox regressions. The univariate prognostic analysis of the training group showed that the clinical stage and risk score of the disease were independent factors affecting the prognosis (p<0.001; Fig. 3I). The multivariate prognostic analysis showed that the risk score (p<0.001) and clinical stage (p<0.05) were independent factors affecting the prognosis, and the risk score had a higher hazard ratio (Fig. 3J). In the test group of patients, both univariate and multivariate prognostic analysis showed that the risk score was an independent factor affecting the prognosis (p<0.01). This showed that the risk scoring model we built has good predictive ability.

Validation of the predictive performance of the risk scoring model

The expression levels of the nine RBPs in the training group and the test group were significantly different between the high-risk group and the low-risk group (Fig. 4A, B). Fig. 4C and D show the distribution of the patients’ risk scores. Fig. 4E and F show the relationship between the patient’s survival status, survival time, and risk score. The red dots represent high-risk patients, and the green dots represent low-risk patients. It can be seen that the higher the risk score, the higher the proportion of patients whose follow-up outcome is death and the shorter the follow-up time.

Clinical features of the high-risk group and low-risk group

We obtained the clinical characteristics of the two groups of patients, and performed a statistical analysis of the surgical methods, alpha-fetoprotein values, degree of liver cirrhosis, and other adjuvant treatments. The above-mentioned and other features of the two groups of patients found no significant difference, indicating risk score is an independent predictor of prognosis (Table 3).

Nomogram construction

To better establish the relationship among RBPs’ expression, risk score and patient survival, we developed a nomogram. According to the nomogram, the expression levels of 9 RBPs can be converted into corresponding scores, and then the scores can be added to obtain the total risk score of the patient. The risk score corresponds to the estimated survival rate, including 1-, 2-, 3- and 5-year survival rates. According to the nomogram, the prognostic model can be applied in the clinic, and the long-term survival rate of a single patient can be predicted based on the expression of RBPs (Fig. 4G).

Expression of RBPs in HCC tissue

To further determine the expression of RBPs in liver cancer tissues for constructing this prognostic model, we used the immunohistochemical staining results in the database to show that BRCA1, CSTF2, EZH2 and XPO5 are highly expressed in liver cancer tissues. EIF2AK4 and MRPL54 are expressed at low levels in liver cancer tissues (Fig. 4H).

Verification of expression of EZH2 in tissues

The primers for EZH2 mRNA were designed, and the expression of 20 pairs of HCC and adjacent tissues was further verified by qRT-PCR. EZH2 showed a significant increase in liver cancer (Fig. 5A). In addition to the mRNA level, at the protein level, we also verified the high expression of EZH2 in tumors by western blotting (Fig. 5D). In addition, we also tested the data of EZH2 expression in normal liver cell lines and various HCC cell lines (Fig. 5B). Among them, the Hep3B cell line showed the highest expression of EZH2. We chose the Hep3B cell line for subsequent in vitro experiments.

Effect of EZH2 on the malignant behaviors of liver cancer cells

After transfection with siRNA, the mRNA and protein levels of EZH2 decreased significantly (Fig. 5C, E). We used the scratch test and the Transwell assay to determine whether reducing EZH2 affects the invasion and migration of HCC cells. Compared with the control group (siRNA-control), the migration ability of the HCC cells at the edge of the scratch in the siRNA-EZH2 group was significantly reduced (p<0.05; Fig. 5F). In addition, the number of HCC cells in the si-EZH2 group that passed through the Transwell chamber was decreased significantly (p<0.05; Fig. 5G). CCK8 experiment demonstrated that knock-down of EZH2 reduced the proliferation ability of cells (Fig. 5H). However, flow cytometry did not find a significant effect of EZH2 on cell apoptosis and cell cycle (Fig. 5I, J).

Acknowledgement

The authors would like to thank the medical research center, Zhongnan Hospital of Wuhan University, for providing equipment.

Data sharing statement

All data are available upon request.

Funding

Our work was supported by the Research Fund of the Health Commission of Hubei Province (WJ2021M255); the Cancer Research and Translational Platform Project of Zhongnan Hospital of Wuhan University (ZLYNXM202004); the Translational Medicine and Interdisciplinary Research Joint Fund Project of Zhongnan Hospital of Wuhan University (ZNJC201918); and a grant from the National Key Research and Development Program of China (SQ2019YFC200078/02).