The current study was performed as a retrospective study. The ethics committee of West China Hospital, Sichuan University approved this retrospective analysis, and waived the requirement for informed consent (2020 − 850). Our study was conducted in accordance with the Declaration of Helsinki, TRIPOD, and Guidelines of Luo et al [26, 27].
Patient demographic
All the patients who underwent radical surgical treatment for pCCA in West China Hospital, Sichuan University during January 2010 and December 2018 were identified. Patients with pCCA confirmed histologically were included in the study. Exclusion criterion: lack of preoperative imaging; only had CT images; examined at 3T MRI; movement artifacts and breathing artifacts; accepted preoperative transarterial chemoembolization; lost in follow-up since discharge; postoperative mortality within 90 days; incomplete pathologic reports or lack of preoperative serum test.
Standard patient demographic was collected. The admission notes, operation records, pathologic reports and radiologic findings were reviewed for each patient. Potential prognostic factors that might be associated with the prognosis of pCCA were collected: age; gender; maximum diameter of the tumor; preoperative blood glucose level; platelet count; white blood cell count; neutrophile count; lymphocyte count; monocyte count; albumin level; globulin level; ratio of albumin to globulin; total bilirubin (TBIL), direct bilirubin (DBIL), indirect bilirubin (IBIL), gamma-glutamyl transpeptidase (GGT), aspartate amino transferase (AST), alkaline phosphatase (ALP), alanine aminotransferase (ALT), carcinoembryonic antigen (CEA); carbohydrate antigen 19 − 9 (CA19-9); hepatitis B virus (HBV) infection; preoperative bile duct drainage; Bismuth type; resection margin status; histologic grade; T stage; operation details and survival status. The endpoint of this study was overall survival (OS), which was computed as the interval between the date of operation and the date of death or last follow-up. All laboratory indicators were examined within one week before surgery.
Treatment and follow-up
Surgical procedures were finally determined and conducted according to preoperative multidisciplinary team (MDT) discussion and intraoperative exploration. The operative technique and management were described in our previous study [28]. R1 was defined as microscopically positive. R0 was defined as no residual tumor. Routine histopathological workup was conducted for all resected pCCA by the Department of Pathology. The TNM stage of pCCA was mainly determined by operation and pathologic records.
Patients were followed-up until January 2020. Outpatient follow-up was executed every 2–3 months for the first year post operation, and every 3–6 months thereafter. At each visit, measurement of tumor markers and liver function, examination of CT and/or MRI were performed for assessment.
MRI Protocol and Conventional MRI Analysis
All enrolled patients underwent dynamic contrast-enhanced (DCE) MRI. All MR scans were performed at a 1.5 T MR scanner (Avanto, Siemens Healthcare, Erlangen, Germany) with an 8-channel body array coil. All patients were instructed to fasting for 6–8 h before the examination. Details regarding the MRI retrieval procedure and acquisition parameters are presented in Appendix E1 (online).
All data were transferred to a devoted workstation (Advantage Workstation 4.6, GE Healthcare) for analysis. Two radiologists with 16 (Wei Zhang, reader1) and 15 (Jun Zhang, reader2) years of experience in abdominal MRI diagnosis reviewed the imaging features of all patients independently. The radiologists were aware that the patients had pCCA, but they were blinded to the exact pathologic type and TNM stage. Differences in their findings were resolved by consensus. Conventional MRI findings were assessed focus on invasion of hepatic artery, invasion of portal vein, clinical lymph nodules status (cN stage) and growth patterns of the tumors. Invasion of hepatic artery and portal vein were considered present at following signs: (Ⅰ) overt tumor thrombus in the vessel. (Ⅱ) the vessel was deformed, stenosis or occluded. (Ⅲ) 50% or more of the circumference of vessel was surrounded by tumor. (Ⅳ) the border between the tumor and the relevant vessel was irregular. cN stage was considered according to the 8th edition of AJCC staging system. Stage cN0 represents free of metastatic regional LNs. Stage cN1 was defined as patients with 1–3 metastatic regional LNs. Stage cN2 was defined as > 3 metastatic regional LNs [3]. Metastatic regional LNs was defined as: (Ⅰ) the lymph node has a short diameter of more than 1 cm. (Ⅱ) fusion of the lymph nodes. (Ⅲ) the occurrence of a centrally necrosis.
Univariate and multivariate COX regression analysis of preoperative variables
After univariate analysis, variables with P < 0.10 were inputted to multivariate analysis. Independent variables were test for the proportional hazards (PH) assumption and used for construction of model.
Image acquisition, Region-of-Interest segmentation
All images were stored in a Digital Imaging and Communications in Medicine (DICOM) format and anonymized. Prior to image analysis, the MR images at contrast enhanced artery phase (CE-AP) and portal vein phase (CE-PVP) were preprocessed by two steps as follows: (1) Intensity standardization (Z-score normalization) was conducted for all the MR images to reduce the variability across different image acquisitions and improve the radiomics feature reproducibility to some extent [29]. (2) Every standardized MR image was resampled to a uniform voxel dimension of 1.0 × 1.0 × 1.0 mm3 before region-of-interest (ROI) segmentation. Image registration for preoperative MR images at artery and portal vein phases was firstly performed with in-house software (Artificial Intelligence Kit, AK, version 3.2.2, GE Healthcare).
ROI segmentation was performed using ITK-SNAP software (version 3.6.0, open source software; https://itk.org/). ROI was delineated for the entire tumor on each slice of both axial CE-AP and CE-PVP images in each patient. To ensure intra- and inter-observer reproducibility, layer-by-layer images segmentation was performed along the tumor contour separately on CE-AP and CE-PVP MR images by two radiologists (twice by reader1 and once by reader2) for 30 patients selected randomly. The two radiologists segmented the images with double-blinded manner. The previously mentioned radiologist (Wei Zhang, reader1) then finished all tumor segmentation.
Radiomics feature extraction and selection, and model development
Radiomics features were extracted from CE-AP and CE-PVP images separately based on segmented tumor ROI by using in-house software (Artificial Intelligence Kit, AK, version 3.2.2, GE Healthcare). For respective CE-AP and CE-PVP images, 396 radiomics features were extracted, including 42 histogram features, 9 morphological features, 144 gray-level co-occurrence matrix (GLCM) features, 180 gray- level run length matrix (GLRLM) features, 11 gray-level size zone matrix (GLSZM) features and 10 Haralick features. The description of the radiomics features which are involved in the radiomics OS predicting models were described in Appendix E2.
Radiomics feature with missing values were firstly replaced with median values. Prior to feature selection, we used the intra- and interclass correlation coefficient (ICC) with the 95% confidence interval (CI) to assess the agreement for each extracted feature. After ICC calculation, features (ICC > 0.7, indicated good agreement) of CE-AP and CE-PVP were kept for further selection. Then the dataset was randomly divided by the ratio of 3:2 into training (110 patients) and testing (74 patients) cohorts, and the feature selection was performed in the training cohort.
The independent significant radiomics features for OS prediction were selected separately based on CE-AP and CE-PVP features. Following the Z-score normalization of radiomics feature in the training and testing cohorts, the redundant features in the training set were excluded by intra-correlation analysis at a cut value of 0.9, which means only one feature was retained when the correlation coefficient between two features larger than 0.9. Next, the univariate COX regression (P < 0.05) and the least absolute shrinkage and selection operator (LASSO) COX regression with 10-fold cross validation was used to select independent features for OS. Two radiomics signatures (SignatureAP and SignaturePVP) were built through a linear combination of selected features weighted by their corresponding coefficients generated from LASSO Cox regression models.
Finally, SignatureAP, SignaturePVP and independent clinical variables were used independently or in combination for construction of different models including: The clinical model (Modelclinic), the AP radiomics model (ModelAP), the PVP radiomics model (ModelPVP) and four combination model (Modelclinic&AP, Modelclinic&PVP, ModelAP&PVP, Modelclinic&AP&PVP).
Assessment of models
Then, we compared the performance of the models. The C-index was employed to assess the discrimination power of the models. The C-index was commonly used to evaluate the discriminative ability of prognostic models in survival analysis. The value of the C-index can range from 0.5 to 1.0 (0.5 indicates no discriminative ability and 1 indicates perfect performance). The best predictive model was determined by performance of comparison and selected for subsequent use. Meanwhile, comparison with the AJCC staging system was carried out using concordance probabilities.
Development of nomogram and validation of the prediction model
To provide a more visualized and individual predictive model, we drafted nomograms for the best predictable model and the clinical model. The discriminant accuracy of the model was further evaluated using the area under curves (AUC) of time-dependent receiver operating characteristic curve (ROC). Calibration curves were assessed by plotting the nomogram-predicted probabilities against the observed rates graphically via a bootstrap method with 1000-iteration resampling. Decision curve analysis (DCA) was employed to evaluate the clinical usefulness of models, through calculation of the net benefit for a range of threshold probabilities [30]. In addition, risk stratification using the constructed model was conducted. Kaplan-Meier survival curves were plotted, and the log-rank tests were used to compare OS between the low-risk and high-risk groups. The performance of the models was then validated in an independent testing set by using the formula generated from the training cohort.
Sample size estimation
As limited by the available sample size and no generally accepted methods to estimate the sample size for developing radiomics-based pCCA risk prediction models, we ensured that the number of positive outcome events in our study was near or larger than 10 times that of covariates (predictors) referring to the TRIPOD statement [26]. In the current study, there are 73 positive outcome events in the training set (110 patients in total) which meet the required condition. Meanwhile, there are 8, 2 and 3 independent predictors retained in the arterial-phase, portal-vein-phase radiomics and preoperative clinical OS predicting model, respectively.
Statistical analysis
All statistical analyses were performed using Statistical Package for Social Sciences software (SPSS, version 25.0, IBM, Armonk, NY, USA), R software (Version: 3.5.3, https: www.r-project.org). Continuous variables were presented as the mean and standard deviation (SD). The Mann-Whitney U-test or the t-test was performed to compare continuous variables. Cut-off values were determined by using the X-tile software (Version: 3.6.1, Yale University, New Haven, CT) [31]. Categorical variables were used Chi-squared (χ2) test or Fisher’s exact test to compare between groups. The Schoenfeld residual test was used to test the proportional hazards (PH) assumption for the selected clinical features [32]. All C-indexes and HRs were reported with 95% CIs. A two-tailed P-value < 0.05 was considered statistically significant. The following R packages were mainly used: “glmnet” for LASSO COX regression algorithm; “survival” for survival analyses including PH assumption of COX regression, COX regression analysis, Kaplan-Meier analysis and decision curve analysis; “tdROC” for time-dependent ROC analyses.