Integrated gene expression analysis of PBMCs
To investigate the transcription profiles, PBMCs from 8 patients with HCC and 4 healthy controls, as well as PBMCs from 3 healthy controls co-cultured with HCC cancer cells (Huh7) were collected to perform RNA-Seq. Baseline characteristics of patients and healthy controls were presented in Supplementary Table S1. Our results showed that a total of 290 genes were identified as differentially expressed genes (DEGs) in PBMCs of patients with HCC compared with healthy controls, which included 213 up-regulated (P<0.05, log2FC>1.5) and 77 down-regulated genes (P<0.05, log2FC<1.5; Fig. 1A). In co-culture of PBMCs with Huh7 cells, we found a total of 367 DEGs with 222 up-regulated genes and 145 down-regulated genes in PBMC co-culture with HCC compared with PBMCs without HCC (Fig. 1B).
To identify whether the DEGs of RNA-Seq data from PBMCs of HCC represent as cancer-induced genes, DEGs from patients with HCC and co-culture model were compared with published microarray data sets (GSE 58208 and 49519)12 using Connection Up-and Down-Regulation Expression Analysis of Microarrays eXtension (CU-DREAMX)13. The intersection among three data sets were shown in Venn diagram (Fig. 1C). Our results showed that a total of 24 DEGs with 18 up-regulated and 6 down-regulated genes were overlapped among these data (Fig. 1C), and intersected genes across data were listed in Supplementary Table S2. To assess whether the above-mentioned 24 intersected genes differed significantly between patients with HCC and healthy controls, a hierarchical clustering using heatmap analysis was further performed. In this respect, the results demonstrated that these two groups could be clearly discriminated by these genes (Fig. 1D). In addition, heatmap analysis of 24 intersected genes in PBMCs from co-culture model and published microarray data were shown in Supplementary Figure S1. Together, these 24 DEGs in PBMCs were speculated as cancer-induced genes due to the alteration of genes when interacting with HCC in the co-culture model and presenting in PBMCs of patients with HCC.
Functional gene annotation and pathway enrichment analysis
Functional enrichment analyses have shown to play an important role in the identification of biological characteristics in transcriptome data. In this study, we performed Gene Ontology (GO) and gProfiler analysis to identify the functional and signaling pathway of DEGs. The majority of DEGs were significantly enriched in the molecular function such as cytokine activity and biological process such as positive regulation of leucocyte cell-cell adhesion, which involved in immune regulation (Fig. 2). Furthermore, enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were similarly involved in immunological response, which included tumor necrosis factor (TNF), pro-inflammatory interleukin (IL)-17 and toll-like receptor signaling pathways.
Selection and validation of candidate biomarkers in PBMCs of HCC patients
To investigate whether the identified DEGs in PBMCs could be used as biomarkers in clinical setting, 5 upregulated candidate genes that were consistently up-regulated across PBMCs samples in our RNA-seq analysis were selected. These genes including BHLHE40, AREG, SOCS1, CCL5 and DDIT4 (green labelled in Fig. 1A-B and 1D) were chosen for validation by qRT-PCR method. The validated cohort consisted of 100 patients with HCC, 100 patients CHB and 100 healthy controls. Baseline characteristics in each group were shown in Table 1. Among these individuals, healthy controls and the non-HCC group were matched for age and gender with patients with HCC. Patients with HCC had significantly higher in aspartate aminotransferase (AST) and AFP levels than patients without HCC. There was no significant difference between groups in terms of serum alanine aminotransferase (ALT) and total bilirubin levels.
Our results demonstrated that BHLHE40, AREG, SOCS1, CCL5 and DDIT4 expression in patients with HCC were higher than those detected in the non-HCC group and healthy controls (Fig. 3). The average relative expression level of BHLHE40 in the HCC group (6.86±4.89) were significantly higher than the non-HCC group (0.62±3.65, P<0.001) and healthy controls (0.00±3.44, P<0.001). Similarly, AREG levels in the HCC group were significantly higher than the non-HCC group (2.34±2.53 vs 0.41±3.92, P<0.001) and healthy controls (0.00±4.88, P<0.001). For SOCS1 expression level, there was significantly higher in patients with HCC than healthy controls (2.55±3.28 vs. 0.00±5.39, P<0.001), but did not reach statistical significance when compared with the non-HCC group (1.13±4.94, P=0.053). For CCL5, the levels of this gene in the HCC group (4.32±4.81) were significantly different from the non-HCC (1.25±4.46, P<0.001) and healthy controls (0.00±4.19, P<0.001). A similar trend was found for DDIT4 expression in the HCC group (6.62±4.48) when compared with the non-HCC group (0.16±3.42, P<0.001) and healthy controls (0.00±3.37, P<0.001). There was no statistically significant difference of these 5 genes between patients without HCC and healthy controls.
Selected candidate genes as diagnostic markers of HCC
To investigate a diagnostic performance of candidate genes in discriminating HCC from non-HCC, the ROC curves were calculated. The area under the ROC curve (AUROC) was 0.83 [95 % confidence interval (CI); 0.78-0.89, P<0.001] for BHLHE40, 0.69 (95 % CI; 0.62-0.77, P<0.001) for AREG, 0.54 (95 % CI; 0.46-0.62, P=0.363) for SOCS1, 0.69 (95 % CI; 0.61-0.76, P<0.001) for CCL5, 0.85 (95 % CI; 0.80-0.90, P<0.001) for DDIT4 and 0.81 (95 % CI; 0.75-0.87, P<0.001) for AFP (Fig. 4).
According to the ROC analysis, BHLHE40 and DDIT4 were considered the best biomarkers among the studied candidate genes, their diagnostic role was further assessed. In this respect, the optimal cut-off value of BHLHE40 in differentiating HCC from non-HCC was 1.80 with a sensitivity of 84.0% and specificity of 64.0%. Similarly, the cut-off value of DDIT4 was 2.10 with a sensitivity of 75.0% and specificity of 76.0%. The combination of BHLHE40 and DDIT4 without AFP slightly increased AUROC to 0.86 and improved accuracy to78.24% which was the highest accuracy in this study (Supplementary Table S3). Our results show that BHLHE had best sensitivity of 84% and AFP had highest specificity of 88%. Furthermore, CU-DREAMX analysis demonstrated that BHLHE40 and DDIT4 of PBMCs were not increased in patients with head and neck cancers and pancreatic cancer (Supplementary Figure S2). These results suggested that BHLHE40 and DDIT4 could be used as potential biomarkers for HCC.
The diagnostic role of BHLHE40 and DDIT4 in AFP-negative HCC and small HCC
Among patients with HCC, there was a strong correlation between BHLHE40 and DDIT4 (r=0.826; P<0.001). However, either BHLHE40 and DDIT4 was not correlated with AFP values (r=0.197; P =0.050 and r=0.154; P=0.126, respectively). Using the normal upper limit of AFP (20 ng/mL) as a reference, there were 51 (51%) and 49 (49%) patients with AFP-negative HCC and AFP-positive HCC, respectively (Fig. 5). Among the AFP-negative group, 78.4% (40/51) of patients had elevated BHLHE40 level (≥1.8) and 66.7% (34/51) of patients had high DDIT4 level (≥2.1). For the AFP-positive group, high levels of BHLHE40 and DDIT4 were detected in 89.8% (44/49) and 83.7% (41/49), respectively. Furthermore, 82.6% (42/51) of patients with AFP-negative HCC had elevated BHLHE40 and/or DDIT4 (Supplementary Figure S3). Among early HCC (stages 0 and A), we found that 27.8% (10/36) patients had elevated AFP concentration, while 86.1% (31/36) and 69.4% (25/35) patients had elevated levels of BHLHE40 and DDIT4, respectively. Together, these results suggested that BHLHE40 and DDIT4 could be promising biomarkers for detecting AFP-negative HCC and early HCC, as well as they could be complementary to AFP in diagnosis of HCC.
Prognostic performance of BHLHE40 and DDIT4 in patients with HCC
The potential prognostic values of BHLHE40 and DDIT4 in terms of overall survival were also analyzed. Based on Kaplan-Meier analysis, the median overall survival of patients with low BHLHE40 levels (<1.8) was significantly better than that of patients whose levels were ≥1.8 (33.8 vs. 15.2 months, P<0.001 by log rank test) (Fig. 6A). Similarly, the median overall survival of patients with low DDIT4 levels (<2.1) was significantly better than that of patients whose levels were elevated (26.7 vs. 15.0 months, P=0.001) (Fig 6B).
BHLHE40 and DDIT4 were entered into the multivariate analysis together with other variables that might influence overall survival of the patients. These factors included age, gender, platelet counts, serum TB, AST, ALT, albumin, AFP level, presence of cirrhosis, tumor size and BCLC stage. The multivariate analysis using the Cox proportional hazards model revealed that high BHLHE40 and BCLC stage were independent predictive factors of overall survival. However, DDIT4 was not selected as an independent factor associated with overall survival (Table 2).