Identification of an eight-gene signature for survival prediction for patients with hepatocellular carcinoma based on integrated bioinformatics analysis

HCC patients and mRNA expression profiles

One dataset including gene expression profiles and associated corresponding clinical information of HCC patients analysed in this study was downloaded from The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov/). The other dataset from GSE14520 with mRNA expression profiling was obtained from the Affymetrix HT Human Genome U133A Array (Affymetrix, Santa Clara, CA, USA).

Gene exclusion criteria were as follows: genes with missing expression values in more than 30% of samples or patients and genes whose expression values were 0 in all samples. The k-nearest neighbour method was used to calculate the remaining missing gene expression values. As the data were from a public database, further approval by an ethics committee was not required. This study met the publication guidelines provided by TCGA (http://cancergenome.nih.gov/publications/publicationguidelines).

Construction of a prognostic PCG signature in the training dataset

Univariate Cox proportional hazards regression analysis was performed to screen out those genes with a significant relationship with patients’ OS in the training dataset (Guo et al., 2018). Then, we used the random survival forests-variable hunting (RSFVH) algorithm to filter prognostic genes until twelve PCGs were screened out (Guo et al., 2016).

Subsequently, based on the above prognostic genes, we performed a multivariate Cox regression analysis and constructed the risk score model as follows: ${\displaystyle \mathrm{Risk Score}\left(RS\right)={\sum }_{i=1}^N\left({\mathrm{Exp}}_i\ast \mathrm{Co}{\mathrm{e}}_i\right)}$ where N is the number of prognostic PCGs, Exp_i is the expression value of PCGs, and Coe_i is the estimated regression coefficient of PCGs in the multivariate Cox regression analysis. Twelve PCGs could form 2¹²−1 = 4,095 combinations or signatures. The corresponding risk scores for the patients from both training and validation dataset were calculated using the risk score system. Receiver operating characteristic (ROC) curves were plotted based on the risk score and survival status of each patient. Then, we removed the signature with the maximum area under the curve (AUC) in the training dataset.

Validation experiments in cell lines

We isolated total RNA from the normal (293T, L02) or cancer (Hep-G2, Hep-3B, PLC/PRF/5, Huh-7) cell lines with a QIAGEN RNeasy Mini Kit. Total RNA was then reverse-transcribed by 5 ×PrimeScript RT Master Mix according to the manufacturer’s recommendation. PCR was performed in the presence of 1 µl complementary DNAs (cDNAs), 0.5 µl forward and reverse primers and 2 ×Taq PCR Master Mix in a total volume of 20 µl. We used GAPDH (Forward: GGAGCGAGATCCCTCCAAAAT, Reverse: GGCTGTTGTCATACTTCTCATGG) as an internal reference. The reaction conditions were carried out according to the manufacturer’s instructions. PCR products were analysed on a 2% agarose gel.

Statistical analysis

Kaplan–Meier (KM) curves were plotted when the median risk score in each dataset was used as the cutoff value to compare survival risk between high-risk and low-risk groups. Multivariate Cox regression analysis was performed to test whether the PCG signature was an independent prognostic factor. Significance was defined as P < 0.05. All were performed in R (R Core Team, 2018) with R packages pROC, survival, TimeROC and randomForestSRC, which were downloaded from Bioconductor.

Furthermore, the co-expressed relationships between the prognostic PCGs in the signature and all other protein-coding genes were computed using the Pearson test, and those genes with P value < 0.05 and absolute value of the Pearson coefficient > 0.5 were selected for Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analyses, which were performed with the clusterProfiler package (Yu et al., 2012).

Patient characteristics and expression profiles

From the TCGA database, 423 HCC samples with mRNA expression profiles and clinical data were downloaded, including 50 normal tissues and 373 cancer tissues. Simultaneously, we obtained the mRNA expression profiles of 445 HCC samples from GSE14520 and their clinical characteristics, including 220 normal and 225 cancers. Then, 332 TCGA and 221 GEO cancer samples with corresponding overall survival data (OS) were selected as training and test sets to explore a prognostic PCG signature and validate the power of the signature in predicting the survival of HCC patients, respectively. After the initial analysis as described in the methods, a total of 16,101 PCG expression values of HCC patients were obtained. All these gene expression values were log₂ transformed. Of the enrolled 553 HCC patients, the median age was 62 years (17–90 years), and 528 were stage I, II, III, and IV, while the stages of 25 patients were unknown. The clinical information of the two datasets is summarized in Table 1.

Identification of prognostic mRNAs from the training dataset

We conducted a univariate Cox proportional hazards regression analysis of the 16,101 PCGs in the training dataset and revealed a subset of 2,231 PCGs that was significantly correlated with patients’ OS (P value < 0.05). To display these selected genes, a volcano plot using the univariate Cox coefficient as the X-axis and −log₁₀ (P value) as the Y-axis was constructed (Fig. 1A). As shown in Fig. 1A, we identified 2,231 genes with significant differences (P < 0.05), which are represented by blue dots; black dots represent the remaining genes with no statistically significant differences. Then, we further performed random forest supervised classification algorithm using the 2,231 genes and screened out 12 PCGs (PVRIG, CSF1, MAFG, S100A9, LRP10, TNNI3K, XRN2, GPR18, DCAF13, EID3, SGCB, and FAM163A) strongly related to patient survival according to the permutation important score by random survival forests-variable hunting (RSFVH) algorithm (Fig. 1B).

Construction of the prognostic PCG signature in the training dataset

These 12 screened PCGs could construct 4,095 risk score models and form 4,095 signatures. To select a signature with the largest prediction power in the training dataset, we performed ROC analyses using patients’ survival status and signature risk scores as variables. Through comparing their areas under the respective ROC curves (AUC), a PCG signature comprising eight genes (PVRIG, S100A9, LRP10, TNNI3K, GPR18, DCAF13, SGCB, and FAM163A) with the max AUC was screened out (Fig. 1C, Table 2).

To validate the expression of these genes, we extracted total RNA from six cell lines (see ‘Methods’). PCR and agarose gel results using the eight gene primers shown in Table 2 revealed that these eight genes were highly abundant in normal and liver cancer-related cell lines (Fig. 1D).

Performance evaluation of the PCG signature for survival prediction

In the training dataset, with the median value of risk score as cutoff, patients were divided into a high-risk group (n = 166) and a low-risk group (n = 166). Kaplan–Meier survival analyses were performed to compare the overall survival of two groups of patients. The low-risk group had significantly better clinical outcomes than the high-risk group (2.20 vs. 8.93 years, log-rank test P < 0.001; Fig. 2A) in the training dataset. The OS rate of patients in the high-risk group was 32.25% and that of patients in the low-risk group was 75.60%.

To validate the prognostic prediction power of the signature, we applied the PCG signature-based risk model to the test set. The median risk scores of 221 HCCs were calculated as the cutoff point in the test set using the Affymetrix gene chip rather than RNA-seq. The test dataset was divided into high-risk and low-risk groups. Kaplan–Meier curves for the high- and low-risk groups in the test dataset are shown in Fig. 2B (median 2.68 vs. 4.24 years, log-rank test P = 0.004, n = 221).

The risk score, gene expression heat-map, and survival status of each HCC patient in the training or test cohort are plotted in Figs. 2C and 2D. For patients with low-risk scores, the expression of PR18, PVRIG, and TNNI3K was upregulated, and DCAF13, FAM163A, LRP10, SGCB, and S100A9 were expressed at low levels.

Comparing the survival prediction power of the PCG signature with other models

In clinical practice, TNM stage is an available prognostic biomarker routinely used as a non-invasive method for predicting the survival of HCC patients. To compare the survival prediction power of TNM stage and the signature, we performed ROC analysis. In the training dataset, the AUCs of the PCG signature were larger than the TNM stage (Signature-AUC Training = 0.73, 95% CI [0.68–0.78] vs. TNM-AUC Training = 0.61, 95% CI [0.54–0.67], Fig. 3A). This result was further verified in both TCGA and GEO datasets (Signature-AUC Test = 0.66 vs. TNM-AUC Test = 0.64, n = 553, Fig. 3B). We also compared the survival predictive power at 3, 5 and 9 years of the PCG signature with TNM staging by TimeROC analysis in both the TCGA dataset and in the entire datasets. In TCGA, the respective AUCs of the PCG signature and TNM staging were 0.78 (0.71–0.84) and 0.65 (0.57–0.73) at 3 years, 0.77 (0.69–0.85) and 0.61 (0.51–0.71) at 5 years and 0.84 (0.76–0.92) and 0.51 (0.22–0.81) at 9 years (Figs. 3C and 3D). In the entire TCGA and GEO dataset, the respective AUCs of the PCG signature and TNM staging were 0.71 (0.67–0.76) and 0.67 (0.62–0.72) at 3 years, 0.69 (0.63–0.75) and 0.67 (0.61–0.74) at 5 years and 0.80 (0.72–0.88) and 0.52 (0.22–0.83) at 9 years (Figs. 3E and 3F).

Moreover, we also compared our eight-gene model with the three-gene prognostic signature (0.384 × RTN3 − 0.561 × SOCS2 − 0.434 × UPB1) (Li et al., 2017) in the entire TCGA and GEO dataset (n = 553), as they were all composed of PCGs. KM analysis showed the two models had good ability in stratification of HCCs (log-rank P < 0.0001); however, our eight-gene signature performed better, as it could divide the 553 HCCs into high-risk and low-risk groups (Fig. 4A), unlike the three-gene model, at more than 7.5 years (Fig. 4B). ROC analysis suggested that our signature outperformed the three-gene model in survival prediction (0.66, 95% CI: 0.62 ∼0.7 vs. 0.6, 95% CI: 0.56 ∼0.64, Fig. 4C). In addition, by TimeROC analysis, we compared the survival predictive power of the three- and eight-gene signatures at 3, 5 and 9 years, showing similar results to those above (Figs. 3E and 3F).

The PCG signature is an independent prognostic factor

To assess whether the PCG signature was an independent risk factor for survival prediction, we performed a Chi-squared test and found no association between the PCG signature with a series of clinical parameters in the training or test groups (Table 3). Then, univariate Cox and multivariate Cox regression analyses were performed in the training dataset and showed that the PCG signature was an independent prognostic factor after adjusting for other clinical features, including sex, age and pTNM (high-risk group vs. low-risk group, HR = 7.23, 95% CI [4.19–12.47], P < 0.001, n = 332, Table 3). The same result was observed in the GEO dataset when adjusting for other clinical features, including sex, age ALT, AFP, cirrhosis, BCLC staging, CLIP staging and TNM staging (high-risk group vs. low-risk group, HR = 1.69, 95% CI [1.06–2.70], P = 0.03, n = 221, Table 3).

Stratification analysis

To obtain a better understanding of the clinical significance of the signature in HCC patients, we correlated the signature with a series of parameters in the two groups (n = 332). As seen in Tables 3 and 4, there was an association between the signature and TNM stage or pathologic T stage (tumour size) in the TCGA and GEO datasets (Chi-square test P < 0.05, Tables 4 and 5). Then, we integrated the TCGA and GEO datasets together and stratified the TNM stage and T stage by the PCG signature risk score. Kaplan–Meier curves showed that patients of low TNM stage (TNM I + II, n = 411) or high TNM stage (TNM III + IV, n = 117) and tumour size ≤ 5 cm (T stage I, II) or tumour size >5 cm (T stage III, IV) were further stratified into two different risk subgroups. The log-rank test showed that the high-risk patients of TNM or T low stage subdivided by the signature had shorter survival than the low-risk patients (log-rank test P < 0.001, Figs. 5A and 5C). Similarly, the TNM or T high stage patients were also divided into a high-risk group with lower OS and a low-risk group with higher OS (log-rank test P = 0.0019/< 0.001, Figs. 5B and 5D). In addition, the GEO set also showed an association between the signature and BCLC staging or AFP; KM analysis indicated that the HCCs could be subgrouped into four clusters by the signature and BCLC staging or alpha fetoprotein (AFP) with different OS: BCLC staging 0/A and low risk, BCLC staging B/C and low risk, BCLC staging 0/A and high risk, and BCLC staging B/C and high risk (log-rank test P < 0.001, Fig. 5E); or AFP ≤ 300 ng/ml and low risk, AFP > 300 ng/ml and low risk, AFP ≤ 300 ng/ml and high risk, and AFP > 300 ng/ml and high risk (log-rank test P = 0.0092, Fig. 5F).

Functional characterization of the selected prognostic PCGs

To explore the functional implications of these selected eight PCGs, we performed Pearson correlation analyses between the eight PCGs and protein-coding genes based on their expression levels in the TCGA and GEO datasets. A total of 776 protein-coding genes were highly correlated with at least one of the selected PCGs (Pearson correlation coefficient >0.5/≤0.5, P< 0.05). GO and KEGG of these co-expressed protein-coding genes were performed and revealed that co-expressed protein-coding genes in the two datasets were significantly enriched in 878 different terms, including 42 KEGG and 836 Go terms (P < 0.05). Most were immunity related and showed the eight genes might be involved in the immune system process through interacting with those co-expressed protein-coding genes that affect important biological processes such as T cell receptor signalling pathway, T cell activity and chemokine signalling pathway (Fig. 6).