DEMs and DEIRGs in patients with HCC
Fig. 1 displays the flow diagram of our study. The analysis of mRNA expression profiles between HCC tissues and adjacent normal liver tissues identified 1416 DEMs (Fig. 2a, b). A total of 1016 mRNAs were significantly upregulated and 400 mRNAs were significantly downregulated in HCC tissues compared to normal liver tissues. The IRGs among the 1416 DEMs were identified, and a total of 90 DEIRGs, including 49 upregulated mRNAs and 41 downregulated mRNAs, were filtered for further analysis (Fig. 2c, d).
Construction of the immune‐related risk signature
In the TCGA cohort, the training set was used to construct the immune-related risk signature. The testing set and entire set were used for validation. Using univariate Cox regression analysis, we identified 33 IRGs with prognostic prediction capacity from a total of 90 DEIRGs (Additional file 2: Figure S1). The 33 IRGs were enrolled in LASSO Cox regression analysis. After 100 rounds of 10-fold cross-validation, 8 optimal IRGs were identified when lambda took the minimum value of 0.059, and the 8 IRGs were used to establish an immune‐related risk signature for HCC patients (Fig. 3a, b). The results of the univariate and multivariate Cox regression analyses of the 8 optimal IRGs are shown in Table 1. The formula of the risk score is as follows:
Table 1 Univariate and multivariate analysis of the 8 optimal IRGs
Gene ID
|
Univariate Analysis
|
|
Univariate Analysis
|
|
HR (95%CI)
|
P-value
|
|
HR (95%CI)
|
P-value
|
APLN
|
1.298 (1.066, 1.582)
|
0.010
|
|
1.178 (0.941, 1.476)
|
0.154
|
CDK4
|
1.721 (1.360, 2.176)
|
0.000
|
|
1.289 (0.967, 1.717)
|
0.083
|
CXCL2
|
0.826 (0.724, 0.943)
|
0.005
|
|
0.924 (0.790, 1.079)
|
0.318
|
ESR1
|
0.514 (0.321, 0.821)
|
0.005
|
|
0.790 (0.479, 1.302)
|
0.355
|
IL1RN
|
0.739 (0.624, 0.875)
|
0.000
|
|
0.819 (0.680, 0.987)
|
0.036
|
PSMD2
|
2.334 (1.610, 3.383)
|
0.000
|
|
1.180 (0.739, 1.886)
|
0.488
|
SEMA3F
|
1.591 (1.101, 2.299)
|
0.013
|
|
1.279 (0.877, 1.865)
|
0.201
|
SPP1
|
1.131 (1.060, 1.208)
|
0.000
|
|
1.078 (0.993, 1.170)
|
0.073
|
Risk score= (0.164*APLN) + (0.254*CDK4) + (-0.079*CXCL2) + (-0.236*ESR1) + (-0.199*IL1RN) + (0.166*PSMD2) + (0.246*SEMA3F) + (0.075*SPP1)
The risk scores for patients were calculated, and patients were divided into high- and low-risk groups using the median risk score of the training set as the cut-off threshold in each set. As shown in Fig. 3c-j, among the 8 genes of the immune-related risk signature, significantly elevated APLN, CDK4, PSMD2, SEMA3F, and SPP1 expression was found in high-risk group compared with low-risk group, while CXCL2, ESR1, and IL1RN expression was significantly decreased in high-risk group compared with low-risk group (all P<0.001). The ROC curves displayed that the AUCs of the immune-related risk signature at 1, 2 and 3 years were 0.813, 0.717, and 0.732 for the training set; 0.631, 0.690, and 0.755 for the testing set; and 0.764, 0.710, and 0.728 for the entire set, respectively (Fig. 4a-c left panel). The distribution of the risk scores, the survival status of the patients ranked according to the risk scores and the expression values of the 8 genes are shown in Fig. 4a-c middle panel. The Kaplan-Meier survival curve combined with the log-rank test showed that patients in the high-risk group had a significantly shorter survival time compared to patients in the low-risk group (all P < 0.05) (Fig. 4a-c right panel). Our results indicated a good performance of the immune-related risk signature for the survival prediction of HCC patients.
External validation of the immune-related risk signature
To avoid bias arising from only investigating the TCGA cohort, we enrolled a total of 232 HCC samples in the ICGC database to validate the predictive value of the signature. We utilized the same formula to calculate the risk score for each patient and then classified the patients into high- and low-risk groups using the same cut-off threshold. As shown in Fig. 4d, the AUCs for 1, 2, and 3 years were 0.704, 0.758, and 0.787, respectively. Consistent with the results in the TCGA cohort, patients in the high-risk group had significantly poorer overall survival than those in the low-risk group. Overall, our signature was able to accurately predict the overall survival of HCC patients.
Independent prognostic role of the immune-related risk signature
A total of 238 HCC patients in the TCGA cohort with complete clinical information (including sex, age, BMI, AFP, histological grade, BCLC stage, and vascular invasion), were included for further analysis. Univariate and multivariate Cox regression analyses indicated that BCLC stage (HR=2.135, P=0.004), vascular invasion (HR=2.142, P=0.003) and risk score (HR=1.612, P<0.001) calculated from the immune-related risk signature were independent prognostic factors for the overall survival of HCC patients (Fig. 5a, b). In addition, the prediction accuracy of the signature was compared with all clinical factors. AS shown in Fig. 5c-e, the signature had the largest AUC at 1, 2, and 3 year, this suggested the predictive performance of the signature was superior to any other clinical factors. Subsequently, all patients were stratified into subgroups based on their clinicopathological factors, such as sex (female/male), age (<60/≥60), BMI (<25/≥25), AFP (<400/≥400), tumor grade (G1+2/G3+4), BCLC stage (I+II/III+IV), and vascular invasion (no/yes), and each subgroup was further divided into the low- and high-risk groups using the immune-related risk signature (Table 2). The results of stratification analysis showed that the high-risk patients in each stratum of all clinicopathological factors had significantly shorter survival times than the low-risk patients (all P < 0.05) (Fig. 5f-s). Our findings suggest that the immune-related risk signature is a promising independent risk factor for the prognosis assessment of HCC.
Table 2 Clinicopathological factors of HCC patients from TCGA cohort in the low- and high-risk groups
Factor
|
Risk score
|
P value
|
|
Low
|
High
|
|
Sex
|
|
|
|
Male
|
81(72.3%)
|
80(63.5%)
|
0.189
|
Female
|
31(27.7%)
|
46(36.5%)
|
|
Age
|
|
|
|
< 60
|
54(48.2%)
|
61(48.4%)
|
1.000
|
≥60
|
58(51.8%)
|
65(51.6%)
|
|
BMI
|
|
|
|
< 25
|
64(57.1%)
|
60(47.6%)
|
0.181
|
≥25
|
48(42.9%)
|
66(52.4%)
|
|
AFP
|
|
|
|
<400
|
95(84.8%)
|
87(69.0%)
|
0.007
|
≥400
|
17(15.2%)
|
39(31.0%)
|
|
Grade
|
|
|
|
G1+2
|
76(67.9%)
|
57(45.2%)
|
0.001
|
G3+4
|
36(32.1%)
|
69(54.8%)
|
|
BCLC stage
|
|
|
|
Ⅰ+Ⅱ
|
96(85.7%)
|
99(78.6%)
|
0.207
|
Ⅲ+Ⅳ
|
16(14.3%)
|
27(21.4%)
|
|
Vascular invasion
|
|
|
|
Yes
|
80(71.4%)
|
72(57.1%)
|
0.031
|
No
|
32(28.6%)
|
54(42.9%)
|
|
Establishment of a new predictive nomogram
We then established a new nomogram to predict 1-, 2-, and 3-year overall survival in 238 HCC patients from the TCGA cohort using three independent prognostic factors, including BCLC stage, vascular invasion and risk score (Fig. 6a). According to the total points of the nomogram, the patients were divided into high- and low-point groups using the median point value as the cut-off threshold. Survival curves for the high- and low-point groups were plotted by Kaplan-Meier analysis, and the results showed that those with higher points had significantly poorer overall survival (P =0.0014) (Fig. 6b). The AUCs of the 1-, 2-, and 3-year overall survival predictions for the nomogram were 0.750, 0.693, and 0.733, respectively (Fig. 6c). The C-index for the nomogram was 0.689 (95% confidence interval [CI] 0.629–0.749). The calibration curve showed that the nomogram performed well at predicting overall survival in HCC patients, but the nomogram might underestimate or overestimate mortality (Fig. 6d-f).
Clinical correlation of the immune-related risk signature
To investigate the ability of our signature to predict the progression of HCC, we analyzed the relationships between the risk factors from our signature (the risk score and 8 risk genes) and clinical parameters (sex, age, BMI, AFP, histological grade, BCLC stage, and vascular invasion) in the entire set. The results indicated that the expression of CDK4 and SEMA3F was higher in female patients, and the expression of ESR1 was higher in male patients (all P < 0.05) (Fig. 7a-c). In patients older than 60 years, the expression of CXCL2 and ESR1 was higher than in those younger than 60 years (both P < 0.05) (Fig. 7d, e). The expression of CXCL2 was higher in patients with BMI ≥ 25 than in those with BMI < 25 (P < 0.05) (Fig. 7f). Patients with AFP ≥ 400 had a higher expression of CDK4 and SEMA3F but a lower expression of CXCL2 and ESR1 than patients with AFP < 400 (Fig. 7g-j). As the values of certain factors increased (CDK4, PSMD2, and SPD1 levels and the risk score), the histological grade of HCC patients increased, and histological grade decreased in patients with increased CXCL2 and ESR1 expression (all P < 0.05) (Fig. 7k-p). And the higher the expression of SEMA3F, the higher the patient’s BCLC stage (P < 0.05) (Fig. 7q). Patients with vascular invasion had higher SPP1 levels and risk scores and lower ESR1 expression than patients without vascular invasion (all P < 0.05) (Fig. 7r-t). These results indicated that the dysregulation of immune-related risk gene expression is correlated with the development of HCC.
Mutation and immune analysis of the signature
As shown in Fig. 8a, the mutation count was greater in the high-risk group than in the low-risk group (P<0.05). Previous studies have suggested that the tumor mutation burden (TMB) is significantly related to the clinical effectiveness of immunotherapy[9]. The results demonstrated that our signature may stratify patients with different sensitivities to immunotherapy. To further identify more immune-related mechanisms of the signature, we analyzed the correlation between the risk score and immune checkpoints, including programmed cell death protein 1 (PD-1), programmed cell death ligand 1 (PD-L1), and cytotoxic T-lymphocyte-associated protein 4 (CTLA-4). We found that HCC patients in the high-risk group tended to have higher expression levels of PD-1 (P<0.01), PD-L1 (P<0.05), and CTLA-4 (P<0.05) (Fig. 8b-d). This finding suggests that high-risk HCC patients might better respond to immune checkpoint blockade (ICB) immunotherapy targeting PD-1 and CTLA-4.The analysis of immune cell infiltration showed that as the risk score increased, the abundances of immune cells (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells) in HCC tissues also increased (all P<0.05) (Fig. 8e-j).
Functional enrichment analysis
To investigate the mechanisms underlying the immune-related risk signature, we examined the expression correlation between the risk genes from the signature and all IRGs using two-sided Pearson correlation coefficients in the entire set. AS shown in Additional file 3: Table S2, a total of 37 IRGs were expressed as significantly correlated with at least one of the risk genes (|Pearson correlation coefficient| > 0.5 and P < 0.01). Then, we performed functional enrichment analysis with the clusterProfiler package to explore the underlying biological processes and pathways of the 37 IRGs, revealing 444 enriched biological processes (Additional file 4: Table S3) and 116 enriched KEGG pathways (Additional file 5: Table S4). The GO results of the top 10 immune‐related biological processes are shown in Fig. 9a,including Fc receptor signaling pathway, the immune response-regulating cell surface receptor signaling pathway, antigen receptor-mediated signaling pathway, regulation of innate immune response, response to tumor necrosis factor, cellular response to interleukin-1, neutrophil mediated immunity, regulation of leukocyte cell-cell adhesion, regulation of T cell activation, and regulation of leukocyte activation. The top 10 immune-related KEGG pathways are shown in Fig. 9b, including the T cell receptor signaling pathway, B cell receptor signaling pathway, viral carcinogenesis, chemokine signaling pathway, natural killer cell-mediated cytotoxicity, PD-L1 expression and PD-1 checkpoint pathway in cancer, IL-17 signaling pathway, NF-kappa B signaling pathway, Th17 cell differentiation, and Th1 and Th2 cell differentiation. Therefore, our results indicated that the immune‐related risk signature displays an intensive immune phenotype.