A retrospective study was conducted on a multicenter database, comprising data from fifteen Chinese tertiary hospitals from January 2012 to December 2021. The patients were categorized into TACE-PVTT0, TACE-PVTT1, and sorafenib-PVTT1 cohorts: (1) TACE-PVTT0, consisting of HCC patients without any portal vein involvement who underwent transarterial chemoembolization (TACE); (2) TACE-PVTT1, comprising patients with Type I PVTT treated with TACE; and (3) sorafenib-PVTT1, including patients with Type I PVTT who received sorafenib monotherapy. The objective was to identify the best group of patients for TACE treatment for HCC within the TACE-PVTT1 group. An RSF model was consequently created for the TACE-PVTT1 group, with 85% of the data set allocated to the training set for model construction, while the remaining 15% was reserved for validation. The mortality risk scores predicted by the RSF model were used to stratify patients into three risk tiers (low, intermediate, and high). Using the tertile partitioning method, we defined the cutoff values at the 33.3rd and 66.7th percentiles of the risk score distribution in the training cohort.

The primary objective was to identify a low-risk subgroup among patients with PVTT1 who optimal candidates for TACE would be. Validation cohorts were used to validate and assess the model’s reliability. TACE-PVTT0 and sorafenib-PVTT1 cohorts were utilized to compare with low-risk TACE-PVTT1, attempting to demonstrate that our selected patients achieved results close to those of the PVTT0 patients who underwent TACE and significantly outperformed sorafenib-PVTT1.

Finally, to further verify the positive survival outcome based on the reduction of selection bias, the low-risk TACE-PVTT1 cohort was compared against the TACE-PVTT0 and sorafenib-PVTT1 cohorts through the integration of PSM and Cox regression analysis (Figure 1).

A review of medical records was conducted to ascertain the clinical characteristics of patients undergone different treatments, and a total of 17 clinical parameters were collected. Tumor burden is an important factor in a range of variables, including the size and number of tumors. When collecting tumor size data, we recorded the maximum diameter of individual tumors.

We employed Cheng’s classification for the stratification of PVTT in HCC patients for its comprehensive consideration of the clinical characteristics and treatment practices specific to HCC patients. Chinese Cheng’s classification defines Type I PVTT based on the medical imaging results: the tumor thrombus invades the segmental portal vein or above. If the postoperative pathological result shows that the tumor thrombus is confined to microvascular, it is classified as type I (6, 18).

HCC Patients were diagnosed according to the American Association for the Study of Liver Diseases/European Association for the Study of the Liver guidelines.

The patients were incorporated based on the following inclusion criteria: (1) diagnosed by HBV-related HCC; (2) patients with Type I PVTT (tumor thrombus limited to segmental or sectoral branches of the portal vein) who were treated by TACE or sorafenib; (3) patients without PVTT treated by TACE; (4) Child-Pugh class A or B7 (score ≤ 7); (5) Eastern Cooperative Oncology Group (ECOG) performance status of 0, 1 or 2.

The exclusion criteria were: (1) irreversible coagulopathy; (2) tumor thrombus extending into the first-order branches (Type II), the main portal vein trunk (Type III), or the superior mesenteric vein (Type IV); (3) severe infection or concurrent active hepatitis that cannot be effectively controlled; (4) diffuse tumor infiltration or extensive distant metastasis; (5) the proportion of tumor occupying the total liver size is ≥70% (not an absolute contraindication); (6) significant reduction in peripheral blood white blood cells and platelets, with white blood cells <3×10⁹/L and platelets <50×10⁹/L (not an absolute contraindication); (7) severe renal dysfunction: serum creatinine level of 176.8 μmol/L or creatinine clearance rate <30 mL/min.

The study was approved by the Ethics Committee of the fifteen hospitals. Informed consent was obtained from all patients before collection of data.

During the TACE procedure, a vascular catheter was inserted selectively into the tumor-feeding artery with an injection containing a mixture of doxorubicin (10–50mg) and lipiodol (2–20 ml), cisplatin (10–110 mg), epirubicin (10–50 mg), and oxaliplatin (100–200 mg), which were selected according to the practice of each center. Subsequently, either gelatin sponge or polyvinyl alcohol foam (PVA) particles were introduced, and embolization was monitored until angiography observed reduced blood flow in the tumor arteries. Repeat TACE was allowed when a residual viable tumor or a new lesion was identified in a patient with normal liver function.

Sorafenib was administered orally at a dosage of 400mg twice daily. The treatment should be continued without interruption, ceasing only in the event of disease progression or intolerable side effects. The dose modifications of sorafenib would be individualized based on each patient’s tolerability and toxicity profile. Patients are required to visit the outpatient clinic every 3 or 4 weeks for ongoing assessments of safety and therapeutic response. In instances of drug-related grade 3 or 4 toxicities, sorafenib dosage would be appropriately reduced to mitigate adverse effects while maintaining therapeutic efficacy. This approach ensures patient safety without compromising the potential benefits of treatment.

All patients were followed up one month after treatment with TACE or sorafenib therapy, then at 3-month intervals for the first year and every 6–12 months after that. In clinical practice, the frequency of follow-up appointments is tailored based on the patient’s baseline characteristics and the response to the most recent treatment.

OS was defined as the interval from initial TACE or sorafenib therapy to all-cause death. Patients who survived or lost to follow-up at the last follow-up date were censored.

RSF is an ensemble machine learning method derived from the random forest algorithm, adapted to handle survival data, which may be censored. The method is applied to measure the minimal depth of a variable and rank the variable importance (VIMP) according to building multiple survival trees from bootstrap samples of the data. Under the assumption of minimal depth, the variates exerting a significant influence on predictive outcomes are predominantly those that consistently partition nodes close to the root. This occurs because these variables are responsible for segmenting the most expansive sample of the population. For VIMP, it quantifies the influence each predictor variable has on the outcome in a predictive model, thereby identifying the most impactful features within the dataset. We randomly selected a subset comprising 85% of the patients from the TACE-PVTT1 group to establish our training cohort. Within the training cohort, we employed RSF to ascertain the most influential variables contributing to the construction of decision trees which enables our subsequent work of model optimization. Linear predictive scores obtained by aggregating the predicted values of multiple decision trees inside the model were the indicators we used to quantify the results.

The PSM method was used to reduce the effect of selection bias and adjust for potential confounding factors, such as clinical stage, tumor size, and tumor number. Propensity scores were derived by fitting a logistic regression model based on age, gender, etiology, tumor size, therapy, AFP, ALB, TBIL, AST, ALT, PLT, INR, BUN, Cr, WBC, and HGB levels. The two groups were matched with a caliper width of 0.02, and 1:1 nearest neighbors matching without replacement was performed by the ‘MatchIt’ package in R software. This study used the PSM method to match patients with TACE-PVTT1 patients to those receiving sorafenib-PVTT1. Likewise, PSM was also used to match TACE-PVTT1 patients to TACE-PVTT0 patients.

The categorical variables in patients’ baseline demographics and clinical characteristics were described using frequencies and proportions. Additionally, continuous variables were described using measures of dispersion (median and interquartile range). The receiver operating characteristic Curve (ROC) and area under curve (AUC) were used to evaluate the predictive effectiveness of the model. K-M survival curves were applied to evaluate the probability of survival for different treatment groups, and hazard ratios (HR) and their 95% confidence intervals (CI) were calculated. The multivariate Cox regression analyses, applied to adjust the prognostic factors, were summarized by delineating the coefficient, HR, 95% CI, and corresponding P value for each covariable.

The RSF model was fitted using the “randomForestSRC” R package to predict the OS of patients based on the selected factors and the optimal hyperparameter combination. An “rms” R package was also used to develop a Cox regression analysis model. Two-tailed P-values < 0.05 were considered as statistically significant.

A total of 3948 patients with HCC were enrolled, including 3073 TACE-PVTT0, 763 patients with TACE-PVTT1, and 167 patients with sorafenib-PVTT1.

The baseline demographics and clinical characteristics of patients are shown in Table 1. Compared to the other two cohorts, patients in the TACE-PVTT0 group exhibited a significantly lower tumor burden, characterized by fewer tumor numbers and smaller tumor sizes. Additionally, the TACE-PVTT0 group demonstrated a lower incidence of elevated AFP levels (36.9% > 400 ng/mL) compared to the TACE-PVTT1 and sorafenib-PVTT1 groups (61.3% and 65.3% > 400 ng/mL, respectively). Furthermore, the age structure of the TACE-PVTT0 group is predominantly younger, hinting at a more advantageous demographic profile. In contrast, the sorafenib-PVTT1 group was identified as having the most substantial tumor burden, with an increased tumor number and expansive tumor sizes, signaling a more progressed stage of the disease.

Among the TACE-PVTT0, TACE-PVTT1 and sorafenib-PVTT1 groups, the median OS was 27.1 months, 8.03 months, and 6.9 months, and the median follow-up period was 36.5 months, 27.5 months and 45.1 months respectively. Patients in the TACE-PVTT0 subclassification had a significantly better OS than the others (Supplementary Figure 1).

The model’s predictive capability typically improves when it is constructed using an expanded variable space. The training cohort consisted of 85% of patients with TACE-PVTT1 (n=649) in this study. A comprehensive collection of 17 covariates, encompassing demographic and clinical parameters, was collected and analyzed to identify potential predictive factors. We used the RSF model to calculate the minimal depth of a variable and rank the VIMP. It compares the minimum depth, and a lower minimal depth value indicates that this variable separates a large set of observations and therefore has a larger impact on forest prediction (Figure 2A). The VIMP of predicting patient survival probability is demonstrated in Figure 2B.

We randomly selected features for each sample and trained an RSF model with 1000 corresponding survival trees. The training cohort was established to build the model, and the validation cohort was established for model testing. The out-of-bag (OBB) error rate tended to be stable when the survival trees increased to a certain number (Supplementary Figure 2A). The minimal depth and VIMP methods were integrated to filter the variables better. It reveals that tumor size, tumor number, AST, INR, and age are the predominant factors significantly associated with predictive accuracy (Supplementary Figure 2B).

The validation cohort consisted of 15% of those patients (n=114). Specifically, the AUC for the prediction of survival outcomes over the initial three-year period was 0.89, 0.92, and 0.89 in the training cohort, whereas that in the validation cohort reached a respective 0.88, 0.93, and 0.96, these results demonstrated an outstanding predictive capacity (Figures 3A, B). Furthermore, the calibration curves for these cohorts reveal a strong concordance between the model-predicted survival rates and the actual observed outcomes both the training cohort and the validation cohort showed excellent prediction ability, as presented in Figures 3C, D. The K-M survival curves showed significant differences when TACE PVTT1 group was stratified into low-risk, intermediate-risk and high-risk groups using trichotomy respectively (Supplementary Figure 3A). In addition to the survival probability, the age, ECOG, tumor size, AFP, ALB levels, and other clinical indicators of the low-risk group exhibited statistically significant disparities with medium- and high-risk groups (Supplementary Table 1, P<0.05). Patients were categorized into low, intermediate, and high risk using linear predictive scores, with ROC curves indicating survival probabilities at 1-, 2-, and 3- years, achieving AUCs of 0.84, 0.85, and 0.84, respectively. (Supplementary Figure 3B).

To enhance the interpretability and clinical applicability of our RSF model, we derived a simplified Nomogram based on the four most influential predictors: age, tumor number, INR, and tumor size (Supplementary Figure 4A). Each predictor was assigned a weighted score, and the total score was used to estimate the 1-, 2-, and 3-year survival rates. The calibration curves for the nomogram showed a high degree of consistency between the predicted and observed survival probabilities, indicating robust predictive performance and clinical reliability (Supplementary Figure 4B).

The total cohort was partitioned into training and internal testing sets using an 85:15 ratio. To ensure that this specific allocation did not bias the results, a sensitivity analysis was performed using a standard 80:20 split, which yielded consistent predictive performance (Supplementary Figure 5, AUC range: 0.884–0.945), thereby confirming the robustness of the model across different data-partitioning strategies.

In addition, we compared the OS in TACE-PVTT0 and low-risk TACE-PVTT1 groups which showed no significance (P = 0.19, Figure 4A). However, the effect of selection bias and potential confounding factors still existed, and some significant differences in some covariables were not excluded. Significant differences in clinical and pathological covariates are observed between the two cohorts (Supplementary Table 2).

After performing PSM to minimize selection bias, a total of 252 matched pairs of patients were identified with no significant difference in any covariates between the TACE-PVTT0 and low-risk TACE-PVTT1 groups, and the consistency of OS of the two groups was further confirmed (P = 0.54, Figure 4B). Supplementary Table 3 presents the baseline demographics and clinical characteristics between low-risk TACE-PVTT1 and TACE-PVTT0 groups following PSM. Consequently, patients categorized as low-risk group were deemed optimal candidates for TACE treatment in PVTT1.

To further evaluate the comparative efficacy of TACE-PVTT1 and low-risk TACE-PVTT0, we conducted a comprehensive multivariate Cox analysis. The results presented in Table 2 showed that the presence or absence of PVTT did not exert a significant effect on prognosis (P = 0.08). Conversely, a cohort of covariates, including ECOG, AFP levels, tumor number, and tumor size, demonstrated a substantial correlation with survival risks post-TACE treatment, achieving statistical significance (P < 0.05). However, these differences were eradicated after PSM (Supplementary Table 6).

A comparative analysis between sorafenib-PVTT1 and low-risk TACE-PVTT1 groups was also conducted, revealing a markedly superior OS rate in the low-risk TACE-PVTT1 group as opposed to sorafenib therapy, both pre-and post-PSM, as illustrated in Figure 5 (P < 0.05). These findings demonstrate that the model aids in the precise selection of PVTT1 patients who are more likely to benefit from TACE over sorafenib therapy. We provide evidence that the PSM method effectively adjusts for potential confounding factors, thereby minimizing irrelevant impacts of certain covariates on the comparison (Supplementary Tables 4, 5).

Moreover, we provided insights into the prognostic implications of TACE and sorafenib treatment approaches in the context of PVTT1 through a comprehensive multivariate Cox regression analysis, which showed that the choice of therapy had a significant influence on prognosis (P < 0.05), as detailed in Table 3. Notably, the HR for sorafenib therapy, with TACE as the reference category, is elevated to 3.4, which further escalates to 7.6 after PSM, as documented in Supplementary Table 7. This substantial increase underscores a significant difference in the efficacy of the two treatment modalities. Additionally, the remaining covariates did not exhibit statistically significant differences in their impact on prognosis. These results indicate the significant clinical relevance of our model, demonstrating its utility in guiding therapeutic decision-making and patient management.