This single-center retrospective cohort study enrolled 537 patients with Barcelona Clinic Liver Cancer (BCLC) C stage HCC who received PD-1/PD-L1 inhibitors (atezolizumab, camrelizumab, etc.) as monotherapy or combined with local therapies (TACE, HAIC, radiofrequency ablation), targeted agents (anti-angiogenics, TKIs), or radiotherapy at Sun Yat-sen University First Hospital (January 2019–May 2023). Inclusion criteria: (1) age ≥ 18 years; (2) histologically/radiologically confirmed unresectable or recurrent metastatic advanced HCC (BCLC C); (3) ≥ 2 cycles of PD-1/PD-L1 therapy; (4) ECOG performance score 0-1; (5) Child-Pugh A/B; (6) complete baseline clinical, imaging, and laboratory data. Exclusion criteria: (1) mixed/sarcomatoid HCC; (2) concurrent malignancies; (3) first immunotherapy administered outside our hospital; (4) missing key data or incomplete follow-up; (5) active autoimmune disease requiring systemic therapy. Patients were randomized 6:4 into training (n=322) and validation (n=215) cohorts (Figure 1).

This study was approved by the Institutional Research Ethics Committee of the First Affiliated Hospital of Sun Yat-sen University. Due to the retrospective nature, informed consent was waived. The study was conducted in accordance with the Helsinki Declaration.

Demographic characteristics, treatment-related details, and clinical baseline characteristics from the week prior to the first immunotherapy were collected through the electronic medical record system. Demographic data included age, gender, weight, height; etiology included smoking history, alcohol consumption history, hepatitis B virus (HBV) and hepatitis C virus (HCV) infections; imaging features included liver cirrhosis, lymph node metastasis, extrahepatic metastasis, portal vein tumor thrombus (PVTT), hepatic vein tumor thrombus (HVTT), and ascites; laboratory indicators included liver function [total bilirubin (TBIL), aspartate aminotransferase (AST), alanine aminotransferase (ALT), albumin (ALB), globulin (GLO)], blood routine [hemoglobin (Hb), platelet distribution width (PDW), neutrophil count, lymphocyte count, platelet count, monocyte count], AFP, LDH, and FIB. Based on these indicators, composite scores were calculated, including body mass index (BMI), NLR, platelet-to-lymphocyte ratio (PLR), lymphocyte-to-monocyte ratio (LMR), albumin-to-globulin ratio (AGR), PNI, SII, aminotransferase-to-platelet ratio index (APRI), albumin-bilirubin index (ALBI), aminotransferase-to-neutrophil ratio index (ANRI), Hb, albumin, lymphocyte, and platelet (HALP) score, and the modified Gustave Roussy Immune (GRIm). Standardized Formulas for Composite Scores were as follows:

Patients were assessed every 2–3 months during the first year, and every 6 months thereafter. Follow-up evaluations included routine blood tests, liver function tests, AFP levels, and CT/MRI scans. Two radiologists independently analyzed the imaging results, and in case of any discrepancies, the final decision was made through discussion or by consulting a third-party expert. OS was defined as the time from the start of immunotherapy to death from any cause or the last follow-up, while progression-free survival (PFS) was defined as the time from the start of immunotherapy to radiological progression (RECIST 1.1), death, or the last follow-up. The follow-up period ended on May 30, 2023.

Normality of continuous variables was assessed using the Shapiro-Wilk test. Normally distributed variables were expressed as mean ± standard deviation and compared via Student’s t-test, while non-normally distributed variables were reported as median (interquartile range [IQR]) and analyzed using the Mann-Whitney U test. Categorical variables were presented as frequencies (percentages) and compared via Pearson’s χ² test or Fisher’s exact test, as appropriate. Optimal cutoff values for continuous biomarkers were determined using R packages based on OS. Univariate Cox proportional hazards regression was first performed to identify potential prognostic factors (P < 0.05 threshold for inclusion). Variables meeting this criterion were subsequently incorporated into multivariate Cox regression. To address multicollinearity and enhance model robustness, LASSO regression with 10-fold cross-validation was used to optimize variable selection. The optimal penalty parameter (λ) was determined using the λ.1se criterion, which selects the most parsimonious model whose performance is within one standard error of the minimum cross-validated error. This approach prioritizes model simplicity and generalizability over maximal fitting of the training data. Multicollinearity among the LASSO-selected variables was explicitly assessed by calculating variance inflation factors (VIFs). All VIF values remained below the conservative threshold of 5, confirming that significant multicollinearity was not present in the final model. A prognostic nomogram was then constructed using selected predictors. Discrimination was evaluated using time-dependent ROC curves. Calibration curves assessed agreement between predicted and observed survival probabilities. Clinical utility was quantified via DCA. Risk stratification was performed using the median value of calculated risk scores in the training cohort as the cutoff. This threshold was then applied to the validation cohort to categorize patients into high- and low-risk groups. Survival differences across risk strata were visualized using Kaplan-Meier curves and compared via log-rank tests. All analyses were conducted using SPSS 26.0 and R 4.3.1 (packages: survival, ggplot2, time-ROC.etc). A p-value < 0.05 was considered statistically significant.

The study ultimately included 537 patients, with 322 in the training group and 215 in the validation group (Figure 1). The median PFS was 7.7 months (95% CI: 6.5-9.7), and the median OS was 27.3 months (95% CI: 21.9-35.0) (Supplementary Figure S1).

Baseline demographic and clinical characteristics are summarized in Table 1. The cohort comprised predominantly males (91.1%), with a median age of 53.31 ± 11.17 years and median BMI of 22.87 ± 3.31 kg/m². The majority (85.1%) had HBV infection, while HCV infection was rare (3.4%). Imaging studies demonstrated cirrhosis in 40.6% of cases, PVTT in 50.1%, and extrahepatic metastases in 54.0%. Hepatic function assessments classified 77.7% as Child-Pugh grade A and 76.2% as ALBI grade 2. Elevated AFP levels exceeding 400 ug/L were observed in 50.8% of patients. The most commonly used immunotherapy agents are camrelizumab (36.1%) and tislelizumab (29.4%). The training and validation cohorts demonstrated balanced baseline characteristics across all key variables, with no statistically significant intergroup differences (all p > 0.05).

The optimal cutoff values for baseline prognostic biomarkers associated with OS were determined based on the highest log-rank statistic (Supplementary Figure S2), including the following: age (38 years), BMI (19.78 kg/m²), AGR (0.94), LDH (191.00 U/L), FIB (3.13 g/L), Hb (128.00 g/L), PDW (10.50 fl), NLR (3.34), SII (1356.42), PLR (228.96), LMR (3.77), PNI (43), ANRI (27.01), ALRI (42.36), APRI (1.26), and HALP (17.52). Patients were classified into high-risk and low-risk groups based on these cutoff values for OS (Supplementary Figure S3).

Through univariate Cox regression analysis in the training cohort, we identified 21 clinical factors significantly associated with OS (P < 0.05) (Table 2).

In the multivariate Cox regression model, LDH ≤ 191 U/L (HR = 0.51, 95% CI: 0.28-0.93, P = 0.028) and APRI ≤1.26 (HR = 0.44, 95% CI: 0.23-0.85, P = 0.015) demonstrated significant survival-protective effects, while Liver cirrhosis (HR = 1.72, 95% CI: 1.09-2.70, P = 0.020), Age ≤38 years (HR = 1.92, 95% CI: 1.04-3.56, P = 0.038) and Hb ≤ 128 g/L (HR = 1.71, 95% CI: 1.11-2.63, P = 0.015) was independently associated with a shortened OS (Table 2). Kaplan-Meier survival analysis for OS further validated the prognostic value of these biomarkers (Supplementary Figure S4).

Univariate Cox regression analysis (with a screening threshold of P < 0.2) was first performed to identify prognostic factors associated with OS. Subsequently, LASSO regression with 10-fold cross-validation and the λ.1SE criterion was used to optimize the variable selection. Thirteen independent prognostic factors with non-zero coefficients were identified: age, BMI, Child-Pugh grade, liver cirrhosis, AFP, FIB, Hb, LDH, PDW, PNI, SII, HALP and APRI (Figure 2; Supplementary Figure S5, Supplementary Table S1). The nomogram constructed based on these variables allows for the intuitive quantification of the contribution of each clinical factor to predicting the 3/6/12/24-month survival probability of advanced HCC patients after immunotherapy (Figure 3A).

The nomogram demonstrated good discriminatory ability in the training cohort, with a C-index of 0.72 (95% CI: 0.67–0.78). Time-dependent ROC curve analysis showed that the model predicted the 3/6/12/24-month OS with AUC values of 0.90 (95% CI: 0.84–0.95), 0.81 (95% CI: 0.73–0.89), 0.75 (95% CI: 0.67–0.82), and 0.77 (95% CI: 0.69–0.85) in the training cohort, respectively. The corresponding AUC values in the validation cohort remained stable (0.71 [95% CI:0.61–0.82], 0.74 [95% CI:0.63–0.85], 0.69 [95% CI:0.60–0.78], 0.66 [95% CI:0.54–0.77]) (Figure 3B). Calibration curves demonstrated a high concordance between the predicted survival probabilities and the Kaplan-Meier observed values in both the training and validation cohorts (Figure 3C). DCA further confirmed its clinical utility: when the threshold probability ranged from 20% to 60%, the nomogram showed a significantly higher net clinical benefit compared to traditional single biomarker models in both the training and validation cohorts (Figure 4).

Using the nomogram developed in this study, patients were categorized into low-risk and high-risk groups based on calculated risk factors. In the training cohort, the high-risk group had an OS hazard ratio (HR) of 3.22 (95% CI: 2.11–4.91; P < 0.0001) (Figure 5). In the validation cohort, the HR for OS in the high-risk group was 2.47 (95% CI: 1.48–4.13; P < 0.001). While the nomogram was developed for OS, it also predicts PFS with clinical relevance. In the training cohort, the HR for PFS in the high-risk group was 1.69 (95% CI: 1.27–2.26; P < 0.001) (Figure 5), and in the validation cohort, the HR for PFS was 1.59 (95% CI: 1.11–2.28; P = 0.01).

The model exhibited consistent discriminative ability across different treatment strategies. For the ICI + anti-angiogenic therapy group (n=75), the AUCs at 3, 6, 12, and 24 months were 0.78 (95% CI: 0.67–0.88), 0.95 (0.86–1.00), 0.65 (0.45–0.86), and 0.66 (0.40–0.93), respectively. In the ICI + TKIs group (n=342), the corresponding AUCs were 0.85 (0.77–0.92), 0.78 (0.70–0.86), 0.73 (0.65–0.80), and 0.76 (0.67–0.85) (Supplementary Figure S6). There were no significant differences in OS (P = 0.91) or PFS (P = 0.76) between the treatment groups, confirming that the model’s risk stratification remained consistent regardless of the therapeutic regimen (Supplementary Figure S7).

In the HBV-positive cohort (n=457, 85.1%), the model maintained stable performance, with AUCs of 0.78 (0.70–0.86), 0.79 (0.72–0.86), 0.75 (0.69–0.81), and 0.74 (0.66–0.81) at 3, 6, 12, and 24 months. Calibration curves showed excellent alignment with observed outcomes, and decision curve analysis demonstrated significant net benefit within clinical decision thresholds, further validating the model’s clinical utility in HBV-endemic populations (Supplementary Figure S8).