A novel molecular-clinicopathologic nomogram to improve prognosis prediction of hepatocellular carcinoma

Data processing

The workflow of this article is shown in Figure 1. In the expression profiles of HCC tumors compared with the samples from normal tissues, we identified 669 differentially expressed lncRNAs (DElncNRAs) of |log Fold Change| ≥ 2 and p < 0.05(Supplementary Table 1 and Figure 2A). Of which, 595 lncRNAs were down-regulated and 74 lncRNAs were up-regulated. As shown in Figure 2B, significant differential expression was detected between the tumor and the adjacent normal groups. Subsequently, the DElncRNAs with P <0.05 were selected by univariate COX regression analysis. Therefore, a total of 191 OS-related lncRNAs and 86 recurrence-related lncRNAs were reserved for further study (Figure 2C). After taking the intersection with GSE76427 and GSE116174, 21 recurrence-related lncRNAs and 85 OS-related lncRNAs were finally obtained for classifier development. A total of 312 patients with OS data were randomized 1:1 into two groups, the training cohort (n=156) and the test cohort (n=156). The GSE116174 (n=64) was reserved as a validation cohort for predicting OS. Meanwhile, a total of 269 patients with recurrence data were randomly divided equally into two groups, the training cohort (n= 130) and the test cohort (n=139). GSE76427 (n= 81) was used as a validation cohort to validate recurrence-related models. LASSO COX selection method was applied to training cohort to develop a prediction model (OS: Figure 3A, 3B; recurrence: 3F, 3G). As shown in Supplementary Tables 2, 8 OS-related DElncRNAs and 14 recurrence-related DElncRNAs were identified by the LASSO COX selected model.

Figure 1. The workflow of this work.

Figure 2. Prognostic DElncRNAs identification process. (A) Volcano plot of differentially expressed lncRNAs in TCGA-LICH dataset; (B) Hierarchical clustering of HCC with or without tumor using 669 differentially expressed lncRNAs using Euclidean distance and average linkage clustering; (C) Venn diagram of prognostic DElncRNAs in prognostic lncRNAs (OS/recurrence multivariate cox p < 0.05) and DElncRNAs(|logFC| >2 and padj < 0.05); (D) Venn diagram of lncRNAs related to OS/recurrence. TCGA, The Cancer Genome Atlas; LICH, Liver hepatocellular carcinoma; HCC, hepatocellular carcinoma; DElncRNA, differentially expressed long non-coding RNA; OS, overall survival; LASSO, least absolute shrinkage and selection operator method.

Figure 3. Development and verification of 8-lncRNAs-based and 14-lncRNAs-based classifiers. (A) LASSO coefficient profiles of the 86 Significant difference lncRNAs in OS set. A vertical line is drawn at the value chosen by 10-fold cross-validation; (B) Ten-time cross-validation for tuning parameter selection in the LASSO model; (C, D) Time-dependent ROC curves and calibration curves of 8-lncRNAs-based classifier; (E) Time-dependent ROC curves of Liao’s biomarker for overall survival; (F) LASSO coefficient profiles of the 21 Significant difference lncRNAs in recurrence set, A vertical line is drawn at the value chosen by 10-fold cross-validation; (G) Ten-time cross-validation for tuning parameter selection in the LASSO model; (H) Time-dependent ROC curves and calibration curves of 14-lncRNAs-based classifier. LASSO, least absolute shrinkage and selection operator method; lncRNA, long non-coding RNA; OS, overall survival; ROC, receiver operating characteristic.

Multi-lncRNAs-based classifier

In order to contrive multi-lncRNAs-based classifiers for predicting OS and recurrence in HCC, LASSO COX selection method was performed with the 85 OS related lncRNAs and 21 recurrence related lncRNAs expression data. An 8-lncRNA-based classifier for OS (Figure 3A and 3B) and a 14-lncRNAs-based classifier for recurrence were constructed by training cohort (Figure 3E, 3F). All those lncRNAs are listed in Supplementary Table 2. 8-lncRNAs-based classifier = 0.0299*EXP (AC090921.1) + 0.0125*EXP (AC096637.2) + 0.1838*EXP(AP002478.1) + 0.2221*EXP (C10orf91)+ 0.0437*EXP (LINC01116) + 0.0251*EXP (LINC01224) + 0.0137*EXP (MAFG-DT) -0.1168*EXP(SERTAD4-AS1); 14-lncRNAs-based classifier = -0.0255*EXP (AC004477.1)+0.1647*EXP(AC010307.4)+0.0416*EXP(AC034229.4)+0.1580*EXP (AC209154.1)+0.3958*EXP(C10orf91)+0.0233*EXP(CDKN2A-DT)+0.0037*EXP(CDKN2B-AS1)+0.00057*EXP(FIRRE)-0.1140*EXP(LINC01549)+0.1813*EXP (LINC01572)+0.0958*EXP(MAFA-AS1)+ 0.1348*EXP(MAFG-DT)-0.365*EXP (MIR9-3HG)-0.0761*EXP (SNHG25). All patients were divided into low and high risk groups according to the optimal cut-off value calculated by X-TILE. The optimal cutoff value for the OS-related classifier was 0.2, and for the recurrence-related classifier was 0.1. The Kaplan-Meier log rank test illustrated that there were significant differences in OS and recurrence in the training cohort (Supplementary Figure 1A, 1E), the test cohort (Supplementary Figure 1B, 1F), the TCGA cohort (Supplementary Figure 1C, 1G), and the GEO cohort (Supplementary Figure 1D, 1H).

Patient characteristics

Since the training cohort and test cohort were equally randomly grouped, there was no significant difference or deviation between them. (Supplementary Tables 9–11, Table 1).

Overall survival

In the training cohort, 156 patients were enrolled. As shown in Supplementary Table 9, there were no significant differences were detected in the distribution of age (P = 0.598), neoplasm histologic grade (P = 0.179), vascular invasion (P = 0.872), performance status (P = 0.155), TNM T stage (P=0.523), TNM M stage (P = 0.298), adjacent hepatic tissue inflammation (P = 0.656), liver fibrosis Ishak score category (P = 0.923), family history (P = 0.922), race category (P = 0.968), HBV infection (P = 0.080), HCV infection (P = 0.139), alcohol consumption (P = 0.287), Child-Pugh classification (P = 0.068), AJCC pathological stage (P = 0.521) and gender (P = 0.849).

In the test cohort, the distribution of data was similar to the training cohort. The proportion of patients with performance status (2 + 3) (P = 0.018), TNM N stage (N1) (P = 0.031) and neoplasm histologic Grade (G3+G4) (P = 0.027) in the high-risk group was higher than the low risk group.

In the TCGA cohort, 312 patients were included for further study. The proportion of patients with performance status (2 + 3) (P = 0.018), tumor grade (G3 + G4) (P = 0.012), and HBV infection (P = 0.043) in the high-risk group was higher than the low risk group.

In the GSE116174 cohort, 64 patients were enrolled. As shown in Supplementary Table 9, there were no significant differences were detected in the distribution of age (P = 0.516), gender (P = 0.418), vascular invasion (P = 0.612), HBV infection (P = 0.849), HCV infection (P = 0.139), Alcohol consumption (P = 0.167), and AJCC pathological stage (P = 0.754).

As shown in Supplementary Figure 2A–2D), the AJCC pathological stage, performance status, HBV infection and neoplasm histologic grade are significantly correlated with the 8-lncRNAs-based classifier. The 8-lncRNAs-based classifier scores for the performance status (2 & 3 & 4), stage (III & IV), HBV positive, and tumor grade (G3 & G4) groups were higher than those of the performance status (0 & 1), stage (I & II), HBV negative, and tumor grade (G1 & G2) groups.

Recurrence

In the training cohort, 130 patients were enrolled. As shown in Supplementary Table 10, the proportion of patients with performance status (2 + 3) (P=0.026), Child-Pugh classification (C&D) (P=0.030), TNM T stage (T3 + T4) (P=0.006), and AJCC pathological stage (III & IV) (P=0.023) in the high-risk group was higher than the low risk group.

In the test cohort, the proportion of patients with tumor grade (G3 + G4) (P=0.006) in the high-risk group was higher than the low risk group. The remaining risk factors were not significantly different in distribution compared to the training cohort (Supplementary Table 10).

In the TCGA cohort, a total 269 patients were enrolled. As shown in Supplementary Table 10, The proportion of patients with tumor grade (G3 + G4) (P = 0.005), HBV infection (P = 0.001), TNM T stage (T3 + T4) (P = 0.029), TNM N stage (N1) (P = 0.023), and AJCC pathological stage (III & IV) (P = 0.007) in the high-risk group was higher than the low risk group.

In the GSE76427 cohort, 81 patients were enrolled. As shown in Supplementary Table 10, there were no significant differences were detected in the distribution of age (P = 0.960), AJCC pathological stage (P = 0.303) and gender (P = 0.117).

As shown in Supplementary Figure 2A–2D), the AJCC pathological stage, performance status, HBV infection and neoplasm histologic grade were significantly correlated with the 14-lncRNAs-based classifier. The 14-lncRNAs-based classifier scores for the performance status (2 & 3), stage (III & IV), HBV positive, and tumor grade (G3 & G4) groups were higher than those of the performance status (0 & 1), stage (I & II), HBV negative, and tumor grade (G1 & G2) groups. In addition, we also investigated the relationship between a total of 20 lncRNAs. The results are shown in Supplementary Table 9.

Prognosis value of the lncRNA-based classifiers

Additionally, we assessed the prognostic value of lncRNA-based classifiers.

Overall survival

Cox univariate analysis showed that the Performance Status, the tumor stage, TNM T classification, HBV infection, and the 8-lncRNA-based classifier were correlated with OS, whether in the training cohort, test cohort, or the TCGA cohort. After multivariable adjustment by these variables, Performance Status (HR: 2.589, 95% CI: 1.355-4.947; P = 0.004), TNM M stage (HR: 7.703, 95% CI: 1.603-37.021; P = 0.011), and the lncRNA-based classifier (HR: 15.483, 95% CI: 6.149-38.989; P < 0.001) remained to be powerful and independent factors for OS in the TCGA Cohort (Supplementary Table 11). In addition, multiple lncRNAs-based classifier was still an independent risk factor in the validation cohort (GSE116174).

In the time-dependent ROC curve, the 8-lncRNAs-based classifiers can effectively predict the 1-year, 3-year, and 5-year survival rates, and their AUC is 0.801, 0.789 and 0.755, respectively (Figure 3C). The average predicted probability (predicted survival rate) and Kaplan-Meier estimated (observed survival rate) were plotted, and the dotted line indicated the ideal reference line corresponding to the predicted survival rate and the actual survival rate. The calibration curve of 1-, 3- and 5-year survival probability based on 8-lncRNAs-classifier were in good agreement with the actual observed values. The C-index of 1-year, 3-year, and 5-year were 0.797, 0.751 and 0.746 respectively, indicating that the prediction model had high accuracy (Figure 3D). Compared with tdROC of liao et al., [13] a larger AUC indicated that our model had a good prediction ability (Supplementary Figure 4A–4C).

Recurrence

Cox univariate analysis showed that Performance Status, the tumor stage, TNM T stage, TNM M stage, HBV infection, and the 14-lncRNA-based classifier were correlated with recurrence, whether in the training cohort, test cohort, or the TCGA cohort. After multivariable adjustment by these variables, Performance Status (HR: 1.608, 95% CI: 1.213-2.131; P < 0.001), TNM M stage (HR: 5.782, 95% CI: 1.631-20.501; P = 0.007) and the lncRNA-based classifier (HR: 2.076, 95% CI: 1.457-2.957; P < 0.001) remained to be powerful and independent factors for recurrence in the TCGA cohort (Table 1).

In the time-dependent ROC curve, the 14-lncRNAs-based classifiers can effectively predict the 1-year, 3-year, and 5-year survival rates with AUC of 0.800, 0.686 and 0.789, respectively (Figure 3G). The average predicted probability (predicted survival rate) and Kaplan-Meier estimated (observed survival rate) were plotted, and the dotted line indicated the ideal reference line corresponding to the predicted survival rate and the actual survival rate. The calibration curve of 1-, 3- and 5-year survival probability based on 14-lncRNAs-classifier are in good agreement with the actual observed values. The C-index of 1-year, 3-year, and 5-year were 0.759, 0.717 and 0.719 respectively, indicating that the prediction model had good performance (Figure 3H). Compared with tdROC of liao et al., [13] a larger AUC indicated that our model had a good prediction ability. (Supplementary Figure 4D–4F)

Construction and evaluation of the nomogram

Subsequently, we constructed a gene-clinical nomogram (Figure 4A, 4B) by multivariate cox regression analysis (Supplementary Table 11, Table 1), combined with clinical characterization and lncRNAs-based classifier. TNM M stage, Performance Status and an 8-lncRNAs-based classifier were included in the gene-clinical nomogram of OS, while TNM M stage, Performance Status and 14-lncRNAs-based classifier were included in the gene-clinical nomogram of recurrence. Nomograms can visually predict the prognosis of patients according to their genes and clinical information, and accurately predict the survival and recurrence of patients at 1-, 3-, and 5 years. Moreover, the score of the nomogram was retained for development and validation of the performance of the nomogram risk score (Supplementary Tables 3 and 4). The risk score of OS-nomogram = 1/(-45.629*Classifier - 47*M - 16*Performance Status + 134.689). The risk score of recurrence-nomogram = 1/(-20*Classifier - 43*M -17.3*Performance Status + 94.7). And the accuracy of the nomogram in 1, 3 and 5 years was analyzed by tdROC, and the corresponding calibration curve was drawn (Figure 4).

Figure 4. Development and verification of OS-nomogram and recurrence-nomogram. (A) OS-nomogram based on 8-lncRNAs-based classifier, TNM M classifier and Performance Status; (B) recurrence-nomogram based on 14-lncRNAs-based classifier, TNM M classifier and Performance Status; (C–E) The 1, 3, and 5-year Time-dependent ROC curves compare the prognostic accuracy of the OS-nomogram; (F) 1, 3, and 5 year calibration curve and C-index of the OS-nomogram; (G–I) The 1, 3, and 5-year Time-dependent ROC curves compare the prognostic accuracy of the recurrence-nomogram; (J) 1, 3, and 5 year calibration curve and C-index of the recurrence-nomogram. OS, overall survival; lncRNA, long non-coding RNA; ROC, receiver operating characteristic; C-index, concordance index.

Overall survival

The OS-nomogram based on an 8-lncRNAs-based classifier combined with the TNM M stage and Performance Status has an AUC of 0.898, 0.834, and 0.814 for predicting 1, 3, and 5 years of OS, respectively. The C-index of 1, 3, and 5 years was 0.878, 0.838, and 0.828, respectively (Figure 4C–4F). The results indicated that the combination of the lncRNA-based classifier models, TNM M stage and Performance Status could enhance the capability to predict the prognosis of survival. Kaplan-Meier curve analysis showed that the two groups divided by cutoff value (0.006953) calculated by X-tile were still statistically significant in OS (Supplementary Figure 3A–3C).

Kaplan-Meier curve showed that patients in the training cohort, the test cohort and the TCGA cohort distributed by lncRNA-based classifiers with Performance Status had significantly different prognosis (p < 0.0001, Figure 5A–5C). As shown in Supplementary Table 5, the tests were performed by log rank test between groups.

Figure 5. Kaplan-Meier analysis in the training, validation and whole cohorts according to the lncRNA-based classifiers stratified by clinicopathological risk factors. (A–C) Kaplan-Meier survival curves of LIHC patients with combinations of lncRNA-classifier and TNM T classifier in the training, test and TCGA cohorts for OS; (D–F) Kaplan-Meier survival curves of LIHC patients with combinations of lncRNA-classifier and TNM stage in the training, test and TCGA cohorts for OS. lncRNA, long non-coding RNA; OS, overall survival; LIHC, Liver hepatocellular carcinoma.

Recurrence

The td-ROC showed that the recurrence-nomogram based on classifier, TNM M stage and Performance Status has an AUC of 0.786, 0.711, and 0.752 for predicting 1, 3, and 5 years of recurrence, respectively. The C-index of 1, 3, and 5 years was 0.719, 0.692, and 0.692, respectively (Figure 4G–4J). The Kaplan-Meier curve analysis also indicated that the prognosis of patients stratified by cutoff value (0.01470) calculated by X-tile was significantly different (Supplementary Figure 3D–3F). The lncRNA-based classifiers with TNM stage could distinguish patients in the training cohort, the test cohort and the TCGA cohort into the different risk of recurrence (p< 0.0001, Figure 5D–5F). As shown in Supplementary Table 6, the tests were performed by log rank test between groups.

GSEA identifies KEGG signaling pathway

In order to investigate different activated KEGG signaling pathways in HCC, GSEA was performed on 8 OS-related lncRNAs and 14 recurrence-related lncRNAs expression datasets. We considered the difference as statistically significant when | NES | ≥ 1, NOM p-value <0.01 and FDR q-val <0.25. All significant enrichment pathways were listed in Supplementary Tables 7 and 8. Figure 6 showed the most significant KEGG pathway function enrichment of 8-OS-related lncRNAs. Figure 7 showed the most significant KEGG pathway function enrichment of 14-recurrence-related lncRNAs. The results showed that Aminoacyl TRNA biosynthesis, Arginine and proline metabolism, Basal transcription factor, Base excision repair, Bladder cancer, Cytoplasmic DNA sensing pathway, DNA replication, Epithelial Signaling in Helicobacter Pylori Infection, Focal adhesion, Gap junction, homologous recombination, leukocyte transendothelial migration, Mapk signaling pathway, Nod like receptor signaling pathway, p53 signaling pathway, pathways in cancer, nucleotide excision repair, RNA degradation, cell cycle, spliceosome, VEGF signaling pathway, and WNT signaling pathways showed consistent enrichment in the up-regulated phenotypes of 8-OS-related lncRNAs (Figure 6). The results showed that Basal transcription factor, Base excision repair, DNA replication, mismatch repair, Lysosome, proteasome, homologous recombination, leukocyte transendothelial migration, P53 signaling pathway, pathways in cancer, Nucleotide excision repair, RNA degradation, cell cycle, and spliceosome signaling pathways showed consistent enrichment in the up-regulated phenotypes of 14-recurrence-related lncRNAs (Figure 7).

Figure 6. Gene set enrichment analysis of lncRNAs of overall survival related classifier. GSEA results showed in (A) AC090921.1, (B) AC096637.2, (C) AP002478.1, (D) C10orf91, (E) MAFG−DT, (F) SERTAD4−AS1, (G) LINC01116, and (H) LINC01224. GSEA, Gene set enrichment analysis.

Figure 7. Gene set enrichment analysis of lncRNAs of recurrence related classifier. GSEA results showed in (A) AC004477.1, (B) AC010307.4, (C) AC034229.4, (D) FIRRE, (E) LINC01572, (F) MAFG−DT, (G) CDKN2A−DT, (H) CDKN2B−AS1, and (I) MIR9−3HG. GSEA, Gene set enrichment analysis.

Patient data

Liver Hepatocellular Carcinoma (LIHC) read counts data was downloaded from TCGA, a publicly available portal (up to May 10, 2019, https://tcga-data.nci.nih.gov/tcga/). Forty-nine adjacent non-tumor samples and 368 HCC samples were obtained after the removal of non-HCC patients and patients who lost critical data. Clinical characteristics of patients were obtained from the cBioportal platform (http://cbioportal.org/) [31], A web resource for visual exploration and analysis multidimensional cancer genome data. The exclusion criteria were as follows: 1) not HCC samples; 2) samples with clinical data but without lncRNA sequence data; 3) samples missing important clinical or biological data; and 4) Patients were followed up for less than three months. After the removal of non-HCC patients and patients lacking critical information, 312 patients were finally reserved for further study. The expression matrices for GSE76427 and GSE116174 were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The expression matrices for GSE76427 and GSE116174 were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). GSE116174 contains 64 patients with OS > 90 days. In GSE76427, 81 patients with recurrence follow-up longer than 90 days were used for further analysis.

Data processing

LncRNAs of LIHC set was re-annotated by the gene annotation file " gencode.v30.long_noncoding_RNAs ", which is Downloaded from the gencode website. 15113 lncRNAs were identified from LIHC set. The expression value of each lncRNA was normalized with the TMM function of the ‘limma’ and ‘edgeR’ package for further analysis [32]. We used the ‘edgeR’ package to test all the data to identify lncRNAs that were differentially expressed by |logFC| >2 and padj < 0.05 in the tumor compared with normal samples. Then, the differentially expressed lncRNAs (DElncRNAs) were subjected to univariate Cox regression, and lncRNAs with p < 0.05 were identified as prognostic DElncRNAs for further research. GEO's lncRNAs were identified by the gpl platform annotation file and the Fasta file "gencode.v32.lncRNA_transcripts.fa".3494 lncRNAs were identified from GSE76427 and 5690 lncRNAs were identified from GSE116174. The OS-related prognostic DElncRNAs that intersect with GSE116174 are used for the development of OS-classifiers. Prognostic DElncRNAs associated with recurrences that intersect with GSE76427 are used to construct recurrence-classifiers. Because the detection methods used by the TCGA database and the GEO database are different, the background noise is also different. So we log2 (x + 1) the TCGA data set and normalize it using the zscore method. zscore normalization is also performed in the GEO dataset. The "sva" package is used to remove batch effects between TCGA and GEO datasets.

Construction of lncRNAs classifier

LASSO is a commonly used high-dimensional predictive regression method. The method is a compression estimation. By constructing a penalty function, it can get a more refined model, and make it compress some regression coefficients, that is, the sum of the absolute values of the forcing coefficients is less than a certain fixed value. Also, set some regression coefficients to be zero. Therefore, it retains the advantage of subset contraction and is a biased estimation for processing data with complex collinearity [33]. The lncRNAs related to OS and recurrence were identified by LASSO COX regression model [34]. The regression coefficients (β) of each related lncRNAs are reserved for the development of lncRNAs-based classifier. The lncRNAs based classifier = ∑ EXP(lncRNA) * β. Based on the optimal cut-off value calculated by x-tile software version 3.6.1 (Yale University School of Medicine, New Haven, CT, USA), the LIHC set was divided into high and low risk groups [35]. The time-dependent receiver operating characteristic (tdROC) curve analysis, calibration curve analysis, Kaplan-Meier survival analysis were used to evaluate predictive ability of the models in training cohort, test cohort, TCGA cohort, and the GEO cohort. After that, we began to construct genomic-clinical nomograms to predict the prognosis and mortality of each HCC patient individually [36].

Data analysis

The Chi-square test, COX survival analysis, and other data processing were completed by SPSS 19.0. Kaplan-Meier log rank test was calculated by medcalc (Version 19.0). Time-dependent ROC (tdROC) was used to assess the performance of lncRNA-based classifier with “time ROC” package in R software(Version 3.6.1). And area under ROC (AUC) was used to assess the accuracy of the prediction. Calibration curve and C-index are performed by ‘rms’ package. The larger C-index indicates that the prediction model has better accuracy [37]. ‘Hmisc’, ‘rms’, and ‘survival’ were used to develop a nomogram. When all the hypotheses are P < 0.05, the difference is statistically significant.

Gene set enrichment analysis

To identify the activation of different KEGG signaling pathways in HCC, we conducted GSEA between down-regulated and up-regulated phenotypes. The lncRNAs of the classifier were divided into up- and down-regulated groups by median. Gene Set Enrichment Analysis (GSEA) was performed by JAVA program GSEA 4.0.2 with the MSigDB Collection (c2.cp.kegg.v7.0.symbol). Normalized enrichment score (NES), nominal p-value and false discovery rate (FDR) were used to quantify enrichment magnitude and statistical significance, respectively [38]. When | NES | ≥ 1, FDR q-val <0.25 and NOM p-value <0.01 were considered significant.