The inclusion criteria were as follows: i) patients who underwent gadoxetic acid-enhanced MRI within 2 weeks before surgery and ii) patients who had a postsurgery diagnosis of HCC with definite CK19 status (+/−). The exclusion criteria were as follows: i) concurrent malignancies or distal metastasis; ii) clinical data or MR images were missing or of low quality; and iii) small tumors (maximum diameter <1 cm). The patient enrollment process is illustrated in Supplementary Figure S1 .

The study included patients from the following three independent medical centers: 1) First Affiliated Hospital of Zhejiang University School of Medicine (FAH-ZJU; patients treated between January 2016 and December 2017); 2) Shulan Health (Hangzhou) Hospital (SLH; patients treated between January 2017 and December 2018); and 3) Lishui Central Hospital (LSCH; patients treated between January 2018 and December 2018). FAH-ZJU, from which the majority of patients were enrolled, was set as the training cohort, and SLH and LSCH were the independent external validation cohorts.

The following baseline (presurgery) patient data were collected: age, sex, tumor number, tumor size, hepatitis B virus infection, serum tumor markers, and serum liver and kidney function indices. Histopathologic data included liver cirrhosis, tumor differentiation grade, microvascular invasion, macrovascular invasion, and lymph node metastasis.

A flow diagram of the study is shown in Figure 1 . The study was divided into four parts: 1) patient enrollment; 2) image preprocessing and segmentation; 3) multisequence MRI radiomics feature extraction and model construction; and 4) model validation. The delineation of the region of interest (ROI) is shown in Figure 2 .

Presurgery gadoxetic acid-enhanced MR images were retrieved from the Picture Archiving and Communication Systems of the participating centers. Patients fasted for 4–6 h before MR scanning. The scan was performed after injection of the contrast agent (0.025 mmol/kg gadoxetic acid, 15 ml, 7.04 g; Guangzhou Consun Pharmaceutical Co., Guangzhou, China) at a rate of 2 ml/s via an injector. The contrast agent was flushed with 20 ml saline injected at the same rate. The scanners used at the three medical centers were as follows: FAH-ZJU, 3.0-T MR scanner (Signa HDxt; GE Healthcare, Milwaukee, WI, USA); SLH, 1.5-T or 3.0-T MR scanner (Signa HDxt); and LSCH, 1.5-T MR scanner (Aera; Siemens Healthcare, Erlangen, Germany) or 3.0-T MR scanner (Ingenia; Philips Healthcare, Best, Netherlands).

Tumor delineation was performed using ITK-SNAP software (http://www.itksnap.orge) (25) by two senior radiologists specializing in abdominal diagnoses with 6 and 8 years of work experience. ROIs were manually contoured on each slice in both T2-weighted and diffusion-weighted imaging (DWI) phase images. The two radiologists checked each other’s delineation results for concordance, and a final check was performed by a chief radiologist.

Radiomics features in ROIs were extracted from the MR images using MATLAB 2016a software (MathWorks, Natick, MA, USA). Image preprocessing was performed to minimize gray range variation between different MR scanners, including voxel size resampling and gray level normalization. For each patient, 968 image features were extracted including 484 T2- and 484 DWI-derived features. The features belonged to four categories: intensity, texture, wavelet, and shape. The feature extraction task was performed as described in our previous work (26, 27).

Histopathologic examination and diagnosis were carried out according to World Health Organization criteria. Tumor specimens were obtained along with surrounding liver tissue at a volume ratio of 1:1, fixed with formalin, and embedded in paraffin (28). Positive CK19 expression was defined as membranous or cytoplasmic immunoreactivity in ≥5% of tumor cells (29). Pathology data were obtained from the pathology department of each medical center.

We established three radiomics-based models, i.e., T2, DWI, and combined (both T2 and DWI radiomics features), using an artificial neural network (ANN) algorithm (hidden size = 2, initial random weight = −0.1 to 0.1, weight decay = 5e−4, and maximum number of iterations = 200). We used a two-step process to select optimal radiomics features. We first performed a univariate analysis to identify features showing a significant difference between CK19+ and CK19− groups in the training cohort. The maximum relevance minimum redundancy (MRMR) algorithm (“mRMR” package in R CRAN) was used to rank feature importance. Briefly, input features were ranked by maximizing the mutual information (MI) to class labels and minimizing the MI with other features. To obtain the best prediction model, four different machine learning algorithms, namely, multiple logistic regression (MLR), support vector machine (SVM), random forest (RF), and ANN, were applied to establish radiomics classifiers with optimal features. For each patient, a radiomics score was calculated to determine the probability of positive CK19 status. The diagnostic ability of the classifiers was evaluated based on the area under the receiver operating characteristic (ROC) curve (AUROC). Given the utility of clinical features, we also developed a clinical classifier to determine CK19 status. A model incorporating clinical factors was developed in the training cohort by multivariate logistic regression analysis. Backward stepwise selection was applied using the likelihood ratio test with Akaike’s information criterion (AIC) as the stopping rule.

This retrospective study was approved by the institutional review board of all participating centers (FAH-ZJU, SLH, and LSCH). Written informed consent was waived for the retrospective use of patients’ clinical and medical imaging data.

Differences between groups were evaluated with the t-test or Mann–Whitney U test for continuous variables, with the chi-square test for qualitative variables, and with the Wilcoxon test for nonparametric variables. Data were analyzed using SPSS v20 software (IBM, Armonk, NY, USA) and R v3.4.1 software (R Core Team, Vienna, Austria). All statistical tests were two-sided, and a p-value <0.05 was considered statistically significant.

The study ultimately enrolled 257 HCC patients including 143 from FAH-ZJU (CK19+:CK19− = 64:79), 75 from SLH (CK19+:CK19− = 34:41), and 39 from LSCH (CK19+:CK19− = 8:31). Baseline clinical and histopathologic characteristics of the participants are shown in Table 1 .

We identified 64 T2 and 64 DWI radiomics features that differed significantly between the CK19+ and CK19− groups in the training cohort according to the Wilcoxon test. These features were selected with the MRMR algorithm, and the top 15 were retained for further analysis ( Supplementary Figure S2 ). In the T2 model, eight radiomics features were identified as optimal features; in the DWI model, 11 features were selected to construct the classifier; and in the combined model, the top 12 features were used to develop the predictive model. The selected features included skewness_DWI, energy_LLL_T2, uniformity_DWI, denth_T2, SZE_LHH_DWI, maxpr_LLH_T2, SZE_HLH_DWI, SZE_HLH_DWI, and idmnc_T2, Busyness_LLL_DWI, and RP_HLH_DWI.

Of the three radiomics models, the combined model showed the best predictive performance for CK19+ status with AUROCs of 0.857 [95% confidence interval (CI): 0.792–0.922] in the FAH-ZJU training cohort, 0.726 (95% CI: 0.610–0.842) in the SLH cohort, and 0.790 (95% CI: 0.639–0.941) in the LSCH cohort. We also found that the performance of the DWI model [FAH-ZJU: 0.854 (95% CI: 0.782–0.927); SLH: 0.635 (95% CI: 0.507–0.764); LSCH: 0.734 (95% CI: 0.555–0.913)] was superior to that of the T2 model [FAH-ZJU: 0.802 (95% CI: 0.729–0.875); SLH: 0.623 (95% CI: 0.492–0.754); LSCH: 0.605 (95% CI: 0.418–0.792)]. The ROC curves of the three radiomics models are shown in Figure 3 .

We also developed classifiers using three other machine learning methods, namely, the MLR, SVM, and RF algorithms. A comparison of the predictive performance of the corresponding models in independent cohorts showed that these three machine learning method-derived classifiers were inferior to the ANN classifier. The ROC curves of the four classifiers are shown in Figure 4 .

An optimal clinical model was constructed based on the multivariate analysis that used the minimum AIC. Six clinical variables including sex, lymph node status, total bilirubin, direct bilirubin, α-fetoprotein, and log α-fetoprotein were used for model construction. The predictive ability of the clinical model was worse than that of the radiomics-based model in the training and independent validation cohorts [FAH-ZJU: 0.737 (95% CI: 0.654–0.819); SLH: 0.560 (95% CI: 0.425–0.694); LSCH: 0.585 (95% CI: 0.313–0.857)]. The ROC curves for the models are shown in Figure 3 .