An Apoptotic Gene Signature for the Prognosis of Hepatocellular Carcinoma

Introduction

Liver cancer is the sixth most frequently diagnosed malignant tumor and the fourth leading cause of cancer-related deaths worldwide, with about 841,000 new cases and 782,000 deaths per annum.¹ Hepatocellular carcinoma (HCC) accounts for about 90% of cases of primary liver cancer; the remaining pathologic types are cholangiocarcinoma and mixed hepatocellular-cholangiocarcinoma.² Surgical resection of the tumor completely and safely is the most effective treatment method for HCC, however, most patients are diagnosed at an advanced stage and become unresectable.³ And thus, the 5-year survival rate for HCC patients was reported to be less than 12%.⁴ Although drugs such as sorafenib and immune checkpoint inhibitors could improve patients’ prognoses, the effect is still unsatisfactory.⁵ Therefore, more sensitive biomarkers are earnestly required for the diagnosis and prognosis of HCC patients.

The initial events of HCC include uncontrolled proliferation and reduction in cell death.⁶ Apoptosis, also known as programmed cell death type I, is precisely regulated by numerous checkpoints to maintain equilibrium.⁷ Suppression of apoptosis was reported to cause cancers, autoimmune diseases, propagation of intracellular pathogens, and neurodegenerative diseases.⁸ An elevated level of apoptosis could be observed in acute liver injury, but apoptosis generally occurs to a lesser extent in non-alcoholic fatty liver disease and liver cancer.^9–11 An early study demonstrated that overexpression of SMAD3 down-regulated Bcl-2 and sensitized hepatocytes to apoptosis through the TGF-β pathway in HCC.¹² Remarkably, a decrease in apoptosis is not only responsible for HCC progression but also for tumor resistance to therapies.¹³ And the apoptosis-inducing therapy showed a potential treatment value for HCC, such as the TRAIL-mediated therapy and inhibition of the negative endogenous apoptosis regulators.^14–16

In our present study, the apoptosis-related differentially expressed genes (AR-DEGs) were defined as the intersection of the DEGs and GSEA. A Lasso-COX regression model was then conducted with the AR-DEGs and survival data. When the partial likelihood deviance reached its minimum in the model, five genes (ie, SOX4, UNG, NRP1, IL2RB, and IVNS1ABP) were selected to calculate the risk score. Patients were divided into the high-risk group and low-risk group based on the median risk score. Also, univariate and multivariate COX regression models were conducted using risk scores and clinical characteristics. ROC curve, Calibration curve, and DCA certified the prognostic significance of models. Moreover, the signature including 5 genes was verified with 76 HCC patients in our center. Finally, protein-protein interaction (PPI) and Gene Ontology (GO) enrichment analysis were performed to investigate the potential mechanism. Briefly, we identified a novel apoptotic gene signature and it might serve as a tool for survival prediction, treatment adjustment, and even clinical classification.

Materials and Methods

Data Collection and Process

The GSE14520 dataset was obtained from the Gene Expression Omnibus (GEO) database and used as a training set. The training set consisted of 218 HCC and 221 adjacent normal tissues. All the microarray raw data in the training set were normalized and transformed with log2(x+1). The mRNA sequencing data of 364 patients with HCC and 50 normal samples were accessed from the UCSC Xena (https://tcga.xenahubs.net) and selected as a validation set. Meanwhile, the 1982 apoptotic genes list in “GO_APOPTOTIC_PROCESS” was downloaded from GSEA. Hepatocellular Carcinoma Expression Atlas Database (HCCDB) contained over 4000 clinical samples.¹⁷ Normalized microarray and sequencing data of GSE10143 (HCCDB7), GSE76427 (HCCDB17) and ICGC_LIRI_JP (HCCDB18) were downloaded to increase the validity of data.

Clinical Specimens

Total 76 HCC tissues were collected from the Department of Hepatobiliary and Pancreatic Surgery in The First Affiliated Hospital of Zhengzhou University, including 7 tumor tissues and paired adjacent normal tissues. The 76 HCC patients with follow-up information were defined as the test set. The clinical specimens were immediately stored in the liquid nitrogen for further qRT-PCR and Western blotting assays.

Identification of AR-DEGs

Based on the “Limma” R package, DEGs were obtained from the training and validation sets.¹⁸ And the cutoff values were set as follows: the adjusted P-value < 0.05 and |log2 fold change (FC)| > 1.5. Then, genes were defined as apoptosis-related differentially expressed genes (AR-DEGs) with the “Venn Diagram” R package and enrolled in the subsequent analysis, if they were the intersection of DEGs and apoptosis-related genes.

Construction of the Apoptotic Signature and Risk Score Calculation

A Lasso-penalized Cox regression model was applied in Rstudio software with AR-DEGs to shrink and choose variables further.¹⁹ The “glmnet” R package²⁰ helped determine the optimal penalty parameter λ and partial likelihood deviance. And finally, a five-apoptotic-gene signature was conducted and the risk score was calculated. The calculation formula was as follows:

The Expression_i was the expression data of the five genes and Coef_i represented the corresponding coefficient.

Evaluation of the Prognostic Value of the Apoptotic Gene Signature

The training set was equally divided into the high-risk group (n=109) and low-risk group (n=109) based on the median risk score. And the survival analysis between the high and low-risk groups was performed with the Log rank test. The receiver operating characteristic (ROC) curves were plotted to obtain the area under the curve (AUC) for 3-year and 5-year overall survival with “pROC” R package. Also, the validation set was divided into the high-risk group (n=182) and low-risk group (n=182), and then survival and ROC analyses were repeated for accuracy.

Assessment of the Apoptotic Signature in Subgroup Analysis

To identify the prognostic value of the signature in different populations, subgroup analysis was performed in the training set. The clinical characteristics were stratified into the following groups: male versus female, age younger than 50 years versus 50 years and older, AFP less than 300 ng/mL (low) versus higher than 300 ng/mL (high), ALT less than 50U/L (low) versus higher than 50U/L (high), early BCLC_stage (0+A) versus advanced BCLC_stage (B+C), main tumor size less than 5 cm versus main tumor size bigger than 5 cm, multinodular negative versus multinodular positive, cirrhosis negative versus cirrhosis positive, and early TNM_stage (I+II) versus advanced TNM_stage (ш+IV). Then, survival curves were drawn to compare the median survival time between the high-risk and low-risk subgroups.

Construction and Calibration of Nomogram

To determine the appropriate variables for nomogram, univariate and multivariate Cox regression analyses were conducted in the training and validation sets separately. Nomograms were widely used for cancer prognosis because of their ability to reduce complicated models into a single numerical estimate of the probability. Thus, the independent prognostic variables were included in the nomogram by the “rms” and “ggplot2” R packages.²¹ Then, calibration of the nomogram was assessed visually using 3-year and 5-year calibration plots. The 3-year and 5-year decision curves analysis (DCA) were used to evaluate the net benefits of risk score and traditional indicators.²² The x-axis in DCA was threshold probability, while the y-axis was net benefits.

Quantitative Real-Time PCR (qRT-PCR)

Trizol RNA isolation reagent (Invitrogen, Carlsbad, CA, USA) was used to isolate the total RNA. The samples were reverse transcribed into cDNA with cDNA Synthesis kit (Takara, Shiga, Japan). Then, the qRT-PCR was conducted to detected the gene expressions with SYBR Premix Ex Taq (Takara) and β-actin was chosen as internal reference. And the primers were shown as follows: β-actin, Forward (F): 5ʹ- GGAGCGAGATCCCTCCAAAAT-3ʹ, Reverse (R): 5ʹ-GGCTGTTGTCATACTTCTCATGG-3ʹ; SOX4, Forward (F): 5ʹ- AGCGACAAGATCCCTTTCATTC-3ʹ, Reverse (R): 5ʹ- CGTTGCCGGACTTCACCTT-3ʹ; NRP1, Forward (F): 5ʹ- GGCGCTTTTCGCAACGATAAA-3ʹ, Reverse (R): 5ʹ- TCGCATTTTTCACTTGGGTGAT-3ʹ; UNG, Forward (F): 5ʹ- CCCCACACCAAGTCTTCACC-3ʹ, Reverse (R): 5ʹ- TTGAACACTAAAGCAGAGCCC-3ʹ; IVNS1ABP, Forward (F): 5ʹ- GATGTTCGACTTCAGGTCTGTG-3ʹ, Reverse (R): 5ʹ- CGTGAGAAATTCCATGAGGATCA-3ʹ; IL2RB, Forward (F): 5ʹ- CTGCTTCACCAACCAGGGTTA-3ʹ, Reverse (R): 5ʹ- GGGGTCGTAAGTAAAGTACACCT-3ʹ.

Western Blotting

The samples were lysed in RIPA lysis buffer with protease inhibitor cocktail (Beyotime Biotechnology, Beijing, China) for 30 min on ice. The proteins were separated by SDS-PAGE and then transferred onto 0.45 μm NC membranes (Millipore) at 300 mA for 2 h. The membranes were blocked with 5% non-fat milk and then incubated with primary antibodies. Antibodies used were: GAPDH (60004-1-Ig, Proteintech), SOX4 (ab236557, Abcam), NRP1 (PA1-29513, Invitrogen™), UNG (PA5-20265, Invitrogen™), IVNS1ABP (ab127566, Abcam), and IL2RB (ab61195, Abcam). The membranes were then incubated with appropriate HRP-conjugated secondary anti-rabbit IgG antibodies for 1 h at room temperature after washed four times with TBST. The results were visualized by the ChemiDoc Imaging System (Bio-Rad, American).

Extra Validation of the Apoptotic Signature in HCCDB and Test Set

The signature was checked against the datasets containing follow-up information in the HCCDB (ICGC_LIRI_JP, GSE10143, and GSE76427). Based on the results of qRT-PCR, the risk score in the test set was calculated and the survival probabilities were computed by the Kaplan–Meier method. Furthermore, the genetic alterations and protein levels of the five apoptotic genes were investigated in the cBioPortal for Cancer Genomics and the Human Protein Atlas respectively.^23,24

Protein-Protein Interaction, GO and KEGG Analysis

To uncover underlying mechanisms in the apoptotic gene signature, protein-protein interaction (PPI) prediction was performed and the analytical conditions were as follows: high confidence (interaction score>0.400) and no more than 20 interactors per shell.²⁵ Then, the results of PPI were applied in Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG).²⁶ The outcome was visualized by Cytoscape (Version 3.72) and Rstudio (Version 1.2.5033).

Statistical Analysis

All the statistical analyses were processed with R software and SPSS (IBM, SPSS, version 22), and P-value < 0.05 was considered statistically significant unless otherwise specified. The least absolute shrinkage and selection operator (LASSO) regression was used to construct the apoptotic signature. The Log rank test was applied to analyze overall survival. Univariate and multivariate Cox proportional hazard models were constructed to obtain the hazard ratios and to select independent prognostic variables. The AUC value of greater than 0.70 was considered acceptable The nomogram was conducted with a multivariate Cox proportional hazard model and then visualized. When evaluating the nomogram, we have utilized the calibration curve and the DCA curve to justified the validity of the model.

Results

Identification of DEGs and Establishment of the Apoptotic Signature

According to the criterion, there were 415 DEGs in the training set and 666 DEGs in the validation set (Supplementary Figure 1A and B). The apoptosis-related differentially expressed genes (AR-DEGs) were defined as the intersection of DEGs and apoptotic genes (Figure 1A, Supplementary Table 1). The 36 AR-DEGs were considered as the key genes that played crucial roles in both HCC and apoptosis. Then, AR-DEGs, survival status, and survival time were used to construct Lasso-penalized Cox regression model (Figure 1B). When the partial likelihood deviance in model reached its minimum, the optimal model consisted of five genes, namely SOX4, UNG, NRP1, IL2RB and IVNS1ABP (Figure 1C). Based on the apoptotic signature, the risk score was calculated as follows: risk score=0.364×SOX4+ 0.2101×IVNS1ABP+0.0890×NRP1+0.1165×UNG+(−0.1386×IL2RB). The detailed description and coefficient of the five genes could be found in Supplementary Table 2. Also, unpaired Students t-test suggested that SOX4, UNG, NRP1 and IVNS1ABP were up-regulated and IL2RB was down-regulated in training and validation sets (Figure 1D and E).

Figure 1 Construction and validation of the apoptotic signature. (A) Venn diagram of the apoptotic genes and differentially expressed genes (DEGs) in GSE14520 and TCGA_LIHC. Black bold fonts indicate the apoptosis-related differentially expressed genes (AR-DEGs). (B) Lasso-penalized Cox regression model to shrink and choose variables. (C) Five genes are selected from 36 AR-DEGs in the model when the partial likelihood deviance reached its minimum. The bold dash-dotted line represents the optimal parameter. (D) The mRNA expression levels of SOX4, UNG, NRP1, IL2RB and IVNS1ABP in the training set. (E) The mRNA expression levels of SOX4, UNG, NRP1, IL2RB and IVNS1ABP in the validation set.

The Apoptotic Signature Predicted Survival in HCC

Patients were divided into two groups based on median risk score.

In the training set, HCC patients with high risk scores had shorter median survival time than those with low risk scores (Figure 2A, p<0.0001). And the 3-year AUC and 5-year AUC values were 0.805 and 0.836 respectively, both of which were above 0.7 (Figure 2B). As for the validation set, the low risk score group also tend to survive than the high risk score group (Figure 2C, p=0.0016). The time-dependent ROC revealed the prognostic accuracy of the signature was 0.751 at 3 years, and 0.821 at 5 years (Figure 2D). Moreover, the prognostic value of signature was checked in HCCDB again. The patients with high risk scores had poor prognosis in GSE10143, GSE76427 and ICGC_LIRI_JP (Supplementary Figure 2A–C, p=0.003, 0.00071 and 0.014, respectively). Taken together, the results demonstrated the prognostic value of the apoptotic gene signature in multiple sets.

Figure 2 The survival analysis and time-dependent ROC curves based on the signature-derived risk score. (A) Kaplan–Meier plot reveals that patients at high risk have a worse prognosis in training set. (B) The 3-year and 5-year ROC curves of the risk score in training set. (C) Kaplan–Meier plot reveals that patients at high risk have a worse prognosis in validation set. (D) The 3-year and 5-year ROC curves of the risk score in validation set.

Survival Analyses in Subgroup Using the Signature

The purpose of subgroup survival analysis was to detect whether the apoptotic signature could benefit specific populations. Based on the risk score, patients were stratified according to gender, age, tumor size, cirrhosis, multinodular, AFP, ALT, TNM_stage, and BCLC_stage. The prognostic value of the signature was identified in all the subgroup analyses (all p<0.01, Figure 3A–D, Supplementary Figure 3), except for female groups (Supplementary Figure 3D).

Figure 3 Subgroup survival analysis in specific populations. (A) Survival analysis in patients with high AFP level using risk score. (B) Survival analysis in patients with low AFP level using risk score. (C) Survival analysis in patients with high ALT level using risk score. (D) Survival analysis in patients with low ALT level using risk score. Abbreviations: AFP, alpha-fetoprotein; ALT, alanine aminotransferase.

Construction and Validation of the Nomogram

In the training set, univariate Cox regression analysis revealed that tumor size, cirrhosis, pathological TNM_stage, BCLC_stage, AFP, and signature-derived risk score were the prognostic factor for survival (Table 1). The six significant variables were then included in multivariate Cox regression analysis but only BCLC_stage, cirrhosis, and risk score were still independent prognostic factors (Table 1). And the nomogram was constructed based on BCLC_stage, cirrhosis, and risk score (Figure 4A). When patients met specific conditions, they received corresponding points in the nomogram. The higher total points got, the worse outcome was. Compared with the BCLC_stage or cirrhosis, risk score showed higher net benefit at 3 years and 5 years whatever the threshold probability was (Figure 4B and C). Moreover, calibration curves containing three internal validation groups (70 samples per group) demonstrated that the nomogram performed well. The 45° dotted lines represented an ideal model and the nomogram almost approached the best prediction (Figure 4D and E).

Table 1 Univariate and Multivariate Cox Regression Analyses of OS in GSE14520

Figure 4 Nomogram, calibration curve and decision curve analyses (DCA) in the training cohort. (A) Construction of nomogram to predict overall survival at 3 years and 5 years. (B and C) DCA curves suggest that risk score has higher net benefit than cirrhosis status and BCLC_stage at 3 years and 5 years. (D and E) Calibration curves containing three internal validation groups demonstrate that the nomogram performs well. The 45° dotted lines represent an ideal model and the nomogram almost approaches the best prediction.

Similarly, after the univariate and multivariate Cox regression analyses (Table 2), pathological N_stage, M_stage, and risk score were selected to build nomogram in validation set (Figure 5A). DCA curves at 3 years and 5 years also proved the theoretical accuracy and clinical utility of the apoptotic signature (Figure 5B and C). And we also generated three internal validation groups (120 samples per group) randomly to testify the deviations from the ideal model (Figure 5D and E).

Table 2 Univariate and Multivariate Cox Regression Analyses of OS in TCGA

Figure 5 Nomogram, calibration and DCA curves in the validation cohort. (A) Construction of nomogram to predict overall survival in the validation cohort. (B and C) DCA curves suggest that risk score has a higher net benefit than pathological N_stage and M_stage for survival prediction. (D and E) Calibration curves reveal that the nomogram approaches the ideal model.

Potential Mechanisms of the Apoptotic Signature

Protein-protein interaction prediction was performed using the STRING public database and the cut-off of the combined score was set at 0.4 (high). There were 42 nodes and 269 interactions in the interactome of SOX4, UNG, NRP1, IL2RB and IVNS1ABP (Supplementary Figure 4). The interacting proteins were then analyzed in Gene Ontology (GO) and KEGG databases. The GO provided structured and computable annotation regarding known genes. The potential interaction proteins were annotated in three aspects: cellular component (CC), molecular function (MF), and biological process (BP). The enrichment results suggested that the apoptotic signature mainly mapped to the chromosome and took part in growth factor, death receptor, tumor necrosis factor receptor binding, and so on (Figure 6A). KEGG was originally designed as an integrated database for the biological interpretation of completely sequenced genes. The KEGG enrichment revealed that the signature might play crucial roles in “PI3K−Akt signaling pathway”, “Apoptosis”, “Cell cycle”, “MAPK signaling pathway”, “TNF signaling pathway” and “Ras signaling pathway” (Figure 6B). These results suggested that the signature might participate in apoptosis and related pathways indirectly or directly.

Figure 6 GO and KEGG analysis. (A) Molecular function, biological process, cellular component enrichment analysis. (B) KEGG enrichment analysis.

Extra Validation of the Apoptotic Signature in the Test Set

In the test set, the mRNA expressions of SOX4, UNG, NRP1, IL2RB and IVNS1ABP were validated using qRT-PCR, and the protein expressions level with Western blotting (Figure 7A and B). And the expressions of SOX4, UNG, NRP1 and IVNS1ABP were significantly up-regulated in tumor tissues, while the IL2RB was down-regulated (Figure 7C). The Log rank test suggested that the low risk group also showed a better prognosis in the test set (Figure 7D, p=0.002). Briefly, the signature-derived risk score had a robust and reliable value for prognostic prediction in the training, validation, test and HCCDB sets. In the Human Protein Atlas (HPA), the protein expressions of IVNS1ABP, NRP1, SOX4, and UNG were significantly increased using immunohistochemistry (Supplementary Figure 5A). Unfortunately, the IL2RB was not recorded in HPA. Meanwhile, IVNS1ABP possessed the most frequent genetic alterations (8%), followed by SOX4 (2.4%), NRP1 (2.2%), IL2RB (1.9%) and UNG (0.5%). The most frequent alteration type was amplification and missense (Supplementary Figure 5B).

Figure 7 Extra validation in the test set. (A) Relative mRNA expression of SOX4, UNG, NRP1, IL2RB and IVNS1ABP in 7 HCC and paired normal tissues. (B) Representative Western blot images of the five proteins. (C) The Western blot results are normalized to GAPDH. (D) Survival analysis in the test set. *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001.

Discussion

According to the annual projections of the World Health Organization, over 1 million people will die from liver cancer in 2030.⁴ HCC is the second most lethal tumor after pancreatic cancer, with a 5-year survival rate less than 12%.²⁷ Although the molecular mechanisms of HCC remain unclear, apoptosis plays a key role in tumorigenesis and cellular immortality.²⁸ With the development of high-throughput sequencing, numerous biomarkers were reported,²⁹ but only few studies highlighted the predictive value of multiple gene signatures, especially the apoptotic signature in HCC.

In our present study, AR-DEGs were identified and used to construct Lasso-penalized Cox regression model. The Lasso-COX model avoided the possible collinearity impact and selected five genes (SOX4, UNG, NRP1, IL2RB and IVNS1ABP) to build the apoptotic signature. In training and validation sets, the signature-derived risk score showed robust prognostic values and patients at low risk had longer median survival time than those at high risk (Figure 2). All the AUC values at 3 years and 5 years were over 0.7 in training and validation sets (0.805, 0.836; 0.751, 0.821, respectively), which suggested that the signature shown quite high prognostic accuracy. The signature was also effective in subgroup analyses (age, TNM_stage, BCLC_stage, etc.) and this would benefit specific populations. It is noteworthy that the signature does not perform well in the female population. And one possible reason may be the deficiency of samples to detect the statistical difference.

Our previous study already identified four genes for the prediction of HCC patient survival.³⁰ However, the novel and effective nomogram was plotted and the signature was testified in a larger sample size (a total of 6 datasets) in the present study. The nomogram provides a potential method for surgeons to improve communication with patients. The numerical and visualized survival rate means a more intuitive interpretation of HCC prognosis. The patients who obtain high scores in the nomogram could receive more aggressive treatment approaches. With a nomogram, the outcome and treatment for each patient would be tailored regarding the personal score. In the training and validation sets, the nomogram was plotted respectively using risk score and traditional clinical characteristics (Figures 4A and 5A). And the complicated statistical predictive models were reduced into user-friendly graphical interfaces. The DCA curves at 3 years and 5 years revealed that the risk score had higher net benefit than traditional clinical characteristics at a certain threshold probability. The internal validation also suggested that the nomogram was very near to an ideal model for prognosis prediction (Figure 4D and E; Figure 5D and E). Therefore, the apoptotic signature might have clinical application in the prognosis evaluation when it was applied in nomogram.

In the subsequent research, GO and KEGG enrichment analyses were performed using the five apoptotic genes and their potential interactomes. The results suggested that the genes were enriched in several biological processes and pathways, including cell cycle, apoptosis, Hepatitis B, PI3K−Akt, MAPK, TNF and Ras signaling pathway. Previous studies demonstrated the above processes and pathways were involved in liver cancer.^31–34 Also, the five apoptotic genes were reported to be correlated with HBV infection or liver cancer. Over-expression of SOX4 might promote HCC intrahepatic metastasis by targeting semaphorin 3C.³⁵ And the upregulation of Neuropilin-1 (NRP-1) was related to the boost of arterial blood supply in primary carcinoma and it converted the hepatic vasculature to the arterial phenotype.³⁶ The previous study demonstrated that the decrease of IL2RB might lead to the escape from immune surveillance in HCC.³⁷ To increase the validity of data, the mRNA and protein expressions of the five genes were detected with traditional experimental methods. The results of qRT-PCR and Western blot analyses demonstrated that SOX4, UNG, NRP1 and IVNS1ABP were up-regulated in tumor tissues, while IL2RB was down-regulated (Figure 7A and B). The experimental verification of the five apoptotic genes was consistent with previous bioinformatic analyses.

However, there are still some weaknesses that deserve comment. First, given the complicated processes of HCC, some crucial apoptotic genes might be omitted and thus reduces the reliability of the signature. Then, the apoptotic signature needs to be verified in more independent validation sets with data integrity. This process might address the lack of samples in the subgroup (female, non-cirrhosis). Moreover, the basic experimental studies about the signature were warranted to investigate the elusive mechanisms of HCC.

Briefly, a novel signature was identified for the prognosis prediction of HCC patients (Figure 8). This five-gene signature might reflect the deregulation of apoptosis and provide prognostic biomarkers in HCC. Further basic research and clinical trials would be conducted to confirm our findings.

Figure 8 Flow chart for the experimental design.