Patients

The institutional review boards of three centers approved this retrospective study and waived patient informed consent.

Subsequently, a total of 230 patients with histopathologically confirmed HCC who underwent preoperative gadobenate-enhanced (Gd-BOPTA) MRI before hepatectomy were retrospectively included. Among them, 144 patients from center 1 were identified as the training set, 54 patients from center 2 were identified as the test set, and the validation set comprised 32 patients from center 3. The patient selection process is illustrated in Fig. 1.

Fig. 1

Flowchart shows the patient selection process

Patients who underwent Gd-BOPTA MRI prior to hepatectomy from three centers between August 2015 and October 2022 were retrospectively recruited. The inclusion criteria were as follows: (a) Gd-BOPTA MRI performed within two months before surgery; (b) presence of histologically confirmed HCC with VETC results.; (c) without previous history of HCC-related treatments (chemotherapy, radiotherapy, transarterial chemoembolization, radiofrequency ablation, or liver transplantation). Exclusion criteria included: (a) Patients with distant metastasis, vascular invasion, tumor thrombus, or other malignant tumors; (b) Incomplete clinical and pathological information; (c) Missing image data or severe artifacts affecting imaging quality.

Laboratory examinations and histopathology

The baseline clinical characteristics were retrospectively collected from medical records of three hospitals. Two pathologists with over 10 years of experience in liver pathological analysis independently evaluated all available histological specimens. They were blinded to other information of the patients. The VETC pattern is defined as a continuous lining of CD34 positive endothelium that encapsulate individual tumor clusters, forming a cobweb-like pattern [6,7,8]. The VETC region was assessed semi-quantitatively on a scale in 5% increments [21]. VETC pattern appearing on 55% or more of the tumor surface (ranging from 0%−100%) was defined as VETC positive (VETC +) [22]. In case of disagreement, consensus was reached through discussion.

MR image acquisition

Before the scan, patients fasted for 4–6 h and then received an intravenous bolus injection of Gd-BOPTA at a dose of 0.1 mmol/kg, administered at a rate of 2.0 mL/s, followed by a 20-mL saline flush. T1-weighted image (T1WI), T2-weighted image(T2WI) and diffusion-weighted image (DWI) were collected. Arterial, portal, delayed, and hepatobiliary phases (AP, PVP, DP, HBP) were obtained at 20–30 s, 50–60 s, 90–120 s, and 60 min after the injection of GD-BOPTA, respectively. Detailed scanner and scan parameters for the three centers are summarized in Supplementary Table S1-3.

Radiology feature analysis and clinicoradiological risk factor

MR images were reviewed independently by two abdominal radiologists with 8 and 10 years of MR imaging experience, respectively. If there was any disagreement, a consensus was reached after discussion. The interobserver agreement for qualitative features was evaluated using the kappa test and variables with K coefficient of ≤ 0.85 were removed. The two radiologists independently assessed the following features: (a) tumor diameter, defined as the maximum diameter measured across the largest plane of the tumor, including the capsule [23]; (b) tumor number, classified as solitary (n = 1) or multiple (1 < n ≤ 3) [24]; (c) tumor margin, categorized as smooth or non-smooth. A smooth margin describes tumors that exhibit a round or oval shape with a continuous and even contour. In contrast, a non-smooth margin signifies tumors with an irregular shape, characterized by uneven or budding portions at the tumor's periphery [25]; (d) shape, round- or oval-like were defined as regular, while others are defined as irregular, such as lobulated, star awn, and needle-like; (e) radiological capsule enhancement, which is a peripheral rim of uniform and smooth hyperenhancement in the PVP or DP, classified into three groups (absent, incomplete, complete) [26]; (f) restricted diffusion, identified by hyperintensity on DWI with a b-value of 600 and hypointensity on ADC [27]; (g) rim arterial phase hyperenhancement (APHE), defined as the irregular ring-like enhancement seen around hypovascular central areas of a tumor during the arterial phase of imaging [28]; (h) Nonrim APHE [29]; (i) Arterial peritumoral enhancement, characterized by a crescent-shaped or polygonal area located outside the tumor margin, which displays hyperintensity during the arterial phase (AP) of imaging and isointensity during the portal venous phase (PVP) [30]; (j) nonperipheral"washout", refers to hypointensity of the lesions compared to the surrounding liver parenchyma in any late phase other than the AP [31]; (k) enhancement pattern, classified as typical dynamic enhancement with arterial hyperintensity and PVP/DP washout, or atypical enhancement [32, 33]. (l) intratumoral necrosis, the area in the tumor without enhancement, high signal on T2WI, low signal on T1WI; (m) tumor hypointensity on HBP, presenting as hypointense tumor on HBP images, when compared with the surrounding liver parenchyma [34]; (n) peritumoral hypointensity on HBP, presenting as hypointense areas of liver parenchyma located outside of the tumor margin in crescent shape, wedge shape, or flame-like shape on HBP images [30, 35].

Radiomics analysis of MR images

Radiomics workflow comprised image segmentation, feature extraction, feature selection, and model construction and evaluation (Fig. 2).

Fig. 2

Flowchart of radiomics analysis

Image segmentation

The entire tumor segmentation exercise was accomplished by a radiologist (reader A, 8 years of abdominal imaging experiences) with ITK-SNAP software. Regions of interests were manually drawn on the T1-weighted images (T1WI), T2-weighted images(T2WI), diffusion-weighted images (DWI), enhanced arterial phase (AP), portal venous phase (PVP), delayed phase (DP), and HBP (Hepatobiliary phase) images, covering the whole tumor. The intratumoral region was defined as the area within radiologist-annotated tumor boundaries. The peritumoral area ≤ 2 mm was generated by python algorithm based on the intratumoral ROI (Tumor), a portion of the ROI beyond the liver parenchyma was removed manually. Feature post-fusion is the combination of intratumoral and peritumoral area ≤ 2 mm as two separate ROIs (Tumor + Peri2mm). Feature pre-fusion refers to the delineation of intratumoral and peritumoral area ≤ 2 mm as a whole ROI (TumorPeri2mm). A total of 30 patients were randomly selected, then the same procedure repeated by another radiologist (reader B, 10 years of abdominal imaging experiences). Reader A repeated the same procedure one month later. The inter- and intra-class correlation coefficient (ICC) was calculated.

Radiomics feature extraction and selection

Radiomic features of four ROI modes were extracted using PyRadiomics. Radiomics features were extracted including first-order features, shape features, texture features, wavelet-transformed and laplacian of gaussian features within four ROIs, which were ROI(Tumor), ROI(Peri2mm), ROI(TumorPeri2mm) and ROI(Tumor + Peri2mm). A set of 1046 radiomic features for each MRI phase (T1WI, T2WI, DWI, AP, PVP, DP, HBP) were extracted for each ROI, which means a total of 7322 radiomics features for ROI(Tumor), ROI(Peri2mm), ROI(TumorPeri2mm) and 14,644 radiomics features for ROI(Tumor + Peri2mm). For each ROI, the radiomics features with ICC less than 0.85 were removed. Then KNN (K-NearestNeighbor) imputation was performed for the remaining radiomics features and then they were standardized into a normal distribution with z-scores to eliminate index dimension differences of the data. The selections of radiomics features under four ROI modes were all as follows. First, the Mann–Whitney U or Student t test was employed to ascertain the difference between VETC positive(VETC +) and VETC negative(VETC-). The significant variables with p < 0.05 were selected. Second, we used univariate logistic regression analysis to select the useful radiomics features. The significant variables with p < 0.05 were retained. Third, Spearman’s rank correlation was implemented to eliminate features with a correlation coefficient greater than 0.75. Finally, we applied the least absolute shrinkage and selection operator (LASSO) method, which can reduce the dimension to select the features with the largest amount of information. The Alpha was optimized through the fivefold cross-validation method (Fig. 3).

Fig. 3

Radiomics feature selection using the LASSO regression algorithm. The fivefold cross-validation process was performed to obtain the optimal tuning parameter Alpha. A The optimal alpha value of 10–1 for the subsequent analysis. B Three features from MRI images with non-zero coefficients were chosen for the radiomics model construction according to the optimal alpha

Model construction and evaluationClinicoradiological model construction

First, univariate logistic regression analysis was performed on all clinical variables and MR imaging indicators. Then, features with P values less than 0.05 were included in the multivariate logistic regression analysis to identify independent risk factors associated with VETC and construct a clinicoradiological model.

Radiomics model construction

Machine learning (ML) models were developed by seven different ML classifiers: logistic regression, support vector machine (SVM), decision tree (DT), random forest (RF), Light Gradient Boosting Machine (LightGBM), eXtreme Gradient Boosting (XGBoost) and multilayer perceptron (MLP), which respectively combined the radiomics signatures selected from the four modes of ROIs, and compared the performance between them in three cohorts, respectively. When most of the ML models perform well in three cohorts in that ROI, the ML model with the best performance in corresponding ROI was the ultimate radiomics model.

Fusion model construction

Subsequently, a radiomics score (Rad_score) was calculated, which combining clinicoradiological features constituted the fusion model using multivariate logistic regression.

Model evaluations and comparisons

The discrimination performance of the models was measured by area under the receiver operating characteristics (ROC) curve. Calibration curves were employed to evaluate predictive accuracy. And, decision curve analysis (DCA) was performed to assess the clinical utility of the nomograms by quantifying the net benefits across various threshold probabilities. The performance of three models in three sets was compared by AUC, sensitivity, specificity, PPV (positive predictive value), NPV (negative predictive value) and accuracy. The Net Reclassification Improvement (NRI) and Integrated Discrimination Improvement (IDI) were used to evaluate the improvement in model performance. Finally, we select the model with the best performance in predicting VETC as the optimal predictive model, and a nomogram was built as a graphical presentation.

Statistical analysis

Statistical analyses were performed with IBM SPSS Statistics (version 25.0) and R software (version 4.3.2).

Continuous variables following a normal distribution with homogeneous variance are expressed as mean ± standard deviation (SD) and compared using Student’s t test. While for variables that do not meeting these criteria, which are expressed as median (interquartile range) and compared using the Mann–Whitney U test. Categorical variables were expressed as case numbers (percentages) and compared using the chi-square test. The DeLong test was used to compare the AUCs of these models. All analyses were considered signifcant at p values of less than 0.05 (two-tailed).

SHapley Additive Explanations (SHAP) is a method for interpreting the results of predictive models based on cooperative game theory [36], which involved both global and local strategies. At the global level, two primary visualizations are leveraged: a summary bee swarm plot and a bar graph depicting the average SHAP values for each feature. The bee swarm plot offers a panoramic view of feature significance, illustrating how each feature influences predictions across the entire dataset. Similarly, the bar chart underscores the importance of key features that significantly drive the model's outcomes. For local interpretation, waterfall plots and force plots were used to explain predictions on an individual basis. Waterfall plots are utilized to visually represent how each individual feature contributes to a specific prediction. Meanwhile, force plots offer a more intricate and granular view, displaying the exact impact of each feature on the predicted outcome for a single patient.