Development and validation of a pyradiomics signature to predict initial treatment response and prognosis during transarterial chemoembolization in hepatocellular carcinoma

We aimed to develop and validate a pyradiomics model for preoperative prediction of initial treatment response to transarterial chemoembolization (TACE) in patients with hepatocellular carcinoma (HCC). To this end, computed tomography (CT) images were acquired from multi-centers. Numerous pyradiomics features were extracted and machine learning approach was used to build a model for predicting initial response of TACE treatment. The predictive accuracy, overall survival (OS), and progression-free survival (PFS) were analyzed. Gene Set Enrichment Analysis (GSEA) was further used to explore signaling pathways in The Cancer Genome Atlas (TCGA)-HCC cohort. Overall, 24 of the 1,209 pyradiomic features were selected using the least absolute shrinkage and selection operator (LASSO) algorithm. The pyradiomics signature showed high predictive accuracy across the discovery set (AUC: 0.917, 95% confidence interval [CI]: 86.93-96.39), validation set 1 (AUC: 0.902, 95% CI: 84.81-95.59), and validation set 2 (AUC: 0.911; 95% CI: 83.26-98.98). Based on the classification of pyradiomics model, we found that a group with high values base on pyramidomics score showed good PFS and OS (both P<0.001) and was negatively correlated with glycolysis pathway. The proposed pyradiomics signature could accurately predict initial treatment response and prognosis, which may be helpful for clinicians to better screen patients who are likely to benefit from TACE.

Introduction

Hepatocellular carcinoma (HCC) is a major cause of cancer-related death worldwide (1). Some patients with HCC are ineligible for liver transplantation and surgical resection because the curative surgery cannot be performed (2–4). For these patients with intermediate and advanced-stage HCC, transarterial chemoembolization (TACE) therapy is a promising treatment method following the National Comprehensive Cancer Network (NCCN) clinical practice guideline (5–7).

Initial treatment response has recently been reported to be a powerful indicator of favorable outcomes, such as longer progression-free survival (PFS) and overall survival (OS) (8–10). Studies on tumor burden were evaluated using magnetic resonance imaging (MRI) or computed tomography (CT) and found associations between imaging features (e.g., tumor size, tumor number) and treatment response to TACE in patients with HCC (11–14). However, imaging features have a limited accuracy of subjective judgment, and they do not reflect intra-tumor heterogeneity. Developing a robust and accurate algorithm to select patients who will show initial response to treatment remains challenging and important. Thus, an accurate model to identify patients with an initial response to TACE could be useful to optimize individualized treatment strategy.

Nowadays, radiomics have been a new and promising field that involves the extraction of large quantitative features from radiographic images (15, 16). The radiomics algorithm offers an unprecedented opportunity to improve cancer decision-making in a low-cost and non-invasive manner. Previous studies have shown that radiomics models of radiology images are significantly associated with clinical outcomes in cancer patients (17–21). We previously found that a radiomics model based on CT images could precisely predict microvascular invasion in HCC patients and the machine learning algorithm could be used to predict clinic outcome in cancer (22, 23). However, the standard method in this field was lacked and the potential mechanism of radiomics model was unclear in the HCC. Radiomics extracting from python package was named pyradiomics and provide the chance. The role of pyradiomics predictive models for initial response or prognosis in TACE treatment, and the association with signal transduction pathway remain unexplored.

In this study, based on the preoperatively CT images and machine learning algorithm, we aimed to develop and validate a robust and accurate pyradiomics signature for a noninvasive pretreatment prediction of initial response to TACE in HCC patients. The mutual relationships between pyradiomics score and clinical factors were further analyzed and validated in other patient cohorts. The subgroup analyses, clinical utility and prognosis of pyradiomics model were estimated. Gene Set Enrichment Analysis (GSEA) tool was used to reveal the association between pyradiomics model and Kyoto Encyclopedia of Genes and Genomes (KEGG), which contributes to interpretation of the potential mechanism but not “black box” of these machine learning models.

Materials and methods

Study design and patients

The flowchart of the machine learning model is presented in detail in Figure 1. This was a retrospective study of 313 patients with Barcelona Clinic Liver Cancer (BCLC) stage B HCC who underwent conventional TACE between February 2010 and December 2020. Patients were recruited from the Nanfang Hospital (n=141 patients, discovery set), Sun Yat-sen University Cancer Center (n=121, validation set 1), and the Second Affiliated Hospital of Gui Zhou Medical University (n=51) (validation set 2). The inclusion criteria were radiologically or pathologically confirmed HCC, initial TACE treatment, BCLC sage B, and arterial-phase CT images availability within 7 days before and 30 days after treatment. Patients who underwent loco-regional or whole-body therapies were excluded. According to the modified Response Evaluation Criteria in Solid Tumors (mRECIST), the initial response to TACE was classified as complete response (CR), partial response (PR), stable disease (SD), and progressive disease (PD) by an experienced radiologist as in the previous study (24). Initial treatment response and non-response were strictly defined as CR+PR and SD+PD, respectively. The OS was defined as the time from the start of TACE or hepatectomy treatment until death or the last contact; PFS was defined as the time between the beginning of TACE treatment and the progression or death of the tumor; disease-free survival (DFS) was defined as the time between the beginning of hepatectomy treatment and disease recurrence or death because of the tumor. This study was approved by the three institutional review boards of Nanfang Hospital, the Second Affiliated Hospital of Gui Zhou Medical University, and Sun Yat-sen University Cancer Center.

FIGURE 1

Figure 1 Flowchart for the development of the machine learning model. All patients have radiologically or pathologically proven HCC and undergo CT before TACE therapy. The radiologists manually segment the 3D ROIs. Thereafter, 1,167 features are extracted from the hepatic arterial CT images based on the “pyradiomics” package of python. Using LASSO method, 24 features are selected, and the pyradiomics model is built. Model validity is evaluated in two cohorts. The predictive performance of the pyradiomics model is compared with that of clinical factors. The association between the pyradiomics model and cancer signaling pathways is analyzed using the TCGA-HCC cohort. 3D, three-dimensional; HCC, hepatocellular carcinoma; CT, computed tomography; ROI, region of interest; LASSO, least absolute shrinkage and selection operator algorithm; TCGA, The Cancer Genome Atlas.

TACE procedure, CT acquisition and manual segmentation of the region of interest

TACE was performed under local anesthesia using the traditional femoral approach. TACE was performed under the guidance of digital subtraction angiography (Allura Xper FD 20, Philips) through the left and right hepatic arteries directly through the arteries supplying blood to the tumor when technically feasible. Hepatic arteriography, performed using a 5 Fr (RH or Yashiro) catheter, was first used to assess the location, number, size, and blood supply of the target tumor. The embolic emulsion agent, including epirubicin (30–60 mg), lobaplatin (30–50 mg), and lipiodol (10–30 mL), was injected into the artery supplying the tumor through a 2.7/2.8 Fr microcatheter. Thirty days after treatment, according to the modified Response Evaluation Criteria in Solid Tumors (ver. 1.1).

Contrast-enhanced computed tomography scans were performed as previously described (22). Contrast-enhanced computed tomography (CECT) was performed at hospital using an MDCT scanner and the detail information of CT image acquisition was described in the Supplementary Material. After the routine CT scanning, a contrast agent (Ultravist 370, Bayer SchL/s) was delivered via an injector (Ulrich CT Plus 150, Ulrich Medical, Ulm, Germany); and CECT was performed immediately after injection. Preoperative CT images were collected on the Picture Archiving and Communication System (PACS; Nanfang Hospital Network Center, China), with an optimal window setting adjustment (window width: 300, window level: 50). The CT images were downloaded through the Picture Archiving and Communication System. Two senior radiologists blinded to the treatment results manually segmented the three-dimensional (3D) regions of interest (ROI) in HCC using the ITK-SNAP (version 3.6, https://sourceforge.net/projects/itksnapx64/). Then, they saved and stored the main images and 3D segmented images for extraction of pyradiomics.

Extraction and reproducibility examination of pyradiomics features

MATLAB 2014b (https://ww2.mathworks.cn/) was used to standardize and reconstruct the segmented 3D ROI image. The thickness of the layer was 1 mm. Python 3.6 (https://www.python.org/downloads/release/python-360/) was used to install the package (https://github.com/Radiomics/pyradiomics) and extract the pyradiomics features from 3D images. These values included the texture, shape, size, and wavelet transform of the CT images. The intra- and inter-correlation coefficients (ICCs) of 50 hepatic artery CT images between two observers were used to evaluate the repeatability of pyradiomics feature extraction. To evaluate intra-observer reproducibility, two radiologists independently and manually segmented the ROIs of 50 patients. Meanwhile, to assess repeatability between observers, two readers extracted the high-dimensional pyradiomics features twice with at least 1-week interval. ICCs greater than 0.75 were set to indicate favorable consistency in pyradiomics extraction. These values, which described the texture, shape, size, and wavelet transform of the CT images, could be used to analyze the overall consistency and select the robust pyradiomics features with good reproducibility.

Development and validation of the pyradiomics signature for predicting therapy response

The least absolute shrinkage and selection operator (LASSO) is a powerful algorithm to choose the most important variables from high-dimensional features (22, 23). LASSO is a reduction method that shrinks the regression coefficient to a certain area. The main idea in using LASSO is constructing a first-order penalty function to obtain a refined model through the final determination of some variables coefficient 0 for feature screening. The penalty term of LASSO is:

$${\displaystyle \sum _{i=1}^n\left|wi\right|\le t}$$

This constraint uses the first-order penalty function of absolute value instead of the second-order function of square sum. Although the form is only slightly different, the results are very different. Some of the coefficients would generate to zero. In this study, LASSO based on 5-fold cross-validation was used to select 24 non-zero coefficients. Then, a pyradiomics score (PRS) was calculated based on a logistic method. A pyradiomics signature was consequently developed to predict TACE treatment response. The pyradiomics model’s performance was then evaluated in the discovery and two validation sets using receiver operating characteristic (ROC) analysis. The optimal cut-off value for predicting treatment response was calculated using the Youden’s index. According to the optimal cut-off PRS value, we divided the patients into two groups. The patients with high values (>-0.14) were defined as RS1, and those with low values (≤-0.14) as RS2. In our study, 112 patients from Nanfang Hospital and 29 from the Second Affiliated Hospital of Gui Zhou Medical University had the information of prognosis.

Gene Set Enrichment Analysis in the TCGA-HCC cohort

The purpose of Gene Set Enrichment Analysis (GSEA) is to explore the relationship between the level of imaging score and tumor-related signaling pathways and gene expression. The GSEA analysis data included four files: 1. The Gene Expression profiling data (Expression dataset) was derived from the mRNA Expression profiles of HCC in The Cancer Genome Atlas (TCGA), and contained 20,533 Gene Expression information; 2. The Gene sets contained 2,074 known Gene sets; 3. Chip annotations listed each probe on the DNA Chip and its matching Hugo gene symbols; 4. Phenotype labels were used to categorize samples into two classes for research purposes and to ensure that the order of the samples was consistent with that of the expression spectrum files. TCGA expression spectrum data were downloaded at: https://www.cancer.gov/about-nci/GDC. The data for a total of 46 HCC patients with preoperative CT images in TCGA database were downloaded from The Cancer Imaging Archive (TCIA) (http://www.Cancer.imaging.archive.net/). In this study, according to the value of PRS, 46 patients with HCC were divided into two groups: RS1 and RS2. The two groups were clustered to identify the distinct genes (fold change ≥2.0, P<0.05). Genes were identified using the “edgR” package. Based on the specific genes, pathway analysis was conducted to determine the potential mechanisms for the machine learning model (DAVID, https://david.ncifcrf.gov/).

Statistical analysis

The “pROC” package was used to plot the receiver operating characteristic (ROC) curves. A confidence interval (CI) of 95% for the area under the curve (AUC) was calculated in all cohorts. PRS was evaluated using the Mann-Whitney U test. The Akaike information criterion (AIC) was used to select the optimal model. The AIC is based on entropy and a measure of the goodness of a statistical model. The smaller the AIC, the better the model. The AIC can be expressed as: AIC = (2k-2L)/n. The Kaplan-Meier curves of DFS, PFS, and OS were analyzed using the “survminer” package. Decision curve analysis (DCA) was used to quantify the probabilities of net benefits at different threshold in patients with HCC, plotted by the “dca.R” package. All statistical analyses were performed using the R statistical software version 3.5.0 (R Core Team, 2018) and GraphPad prism 7.0. Two-sided P values<0.05 were considered significant.

Results

Patient characteristics

In total, 18 (12.76%), 13 (10.75%), and 10 (19.61%) patients in the training set, validation set 1, and validation set 2 were females, respectively. Furthermore, 93 (65.96%), 85 (70.25%), and 33 (64.70%) patients were aged less than 60 years, respectively. The baseline patient characteristics are shown in Table 1. Most patients (82.27, 86.78, and 82.35%) had Child-Pugh A disease and a low number of tumors (≤3) (85.82, 86.64, and 76.47%) in the training set, validation set 1, and validation set 2, respectively. Overall, 51.06, 52.06, and 49.01% of the patients in the training set, validation set 1, and validation set 2, respectively, had high alpha-fetoprotein (AFP) levels. Patients with small tumor size (≤5 cm) accounted for 13.47, 12.39, and 15.69% of the population in the discovery set, validation set 1, and validation set 2, respectively. There were 56 (39.72%), 51 (42.15%), and 19 (37.25%) patients in these three cohorts, respectively, who achieved CR/PR. There was no significant difference in the treatment response rate between the three sets.

TABLE 1

Table 1 Patient characteristics by study set.

Pyradiomics signature development and associated analysis of clinical factors

Pyradiomics features were extracted as described in previous studies (25–27). A total of 1,167 features were extracted from the hepatic arterial 3D-CT images. A total of 457 pyradiomics features were eliminated in the ICC analysis. The remaining 710 features were then subjected to feature selection and LASSO coefficient analysis. Based on 5-fold cross-validation via the maximum criteria, 24 coefficients were selected (Supplemental Figure S1). Twenty-four pyradiomics features were analyzed via multi-variable logistic regression and included to develop the pyradiomics signature (PRS) (Table 2).

TABLE 2

Table 2 Formula for calculating pyradiomics signature.

Treatment responders showed significantly higher PRS than non-responders in the three cohorts (all P<0.001) (Figures 2A-C). Tumor size was also associated with treatment response in the three cohorts (P<0.001, P<0.001, and P=0.017, respectively) (Supplemental Figure S2). However, we found there was no association between treatment response and other clinical factors, such as tumor number, AFP level. Meanwhile, PRS was also not correlated with age, sex, Child-Pugh classification, AFP level, and tumor number across all three cohorts (Figures 2D-F). In contrast, PRS was significantly associated with tumor size (r=-0.416, P<0.001; r=-0.514, P<0.001; r=-0.568, P<0.001, respectively) and treatment response (r=0.605, P<0.001; r=0.539, P<0.001; r=0.588, P<0.001, respectively).

FIGURE 2

Figure 2 Association between the pyradiomics score and clinical factors. (A-C) Correlations between the pyradiomics score and object response (CR+PR) in the discovery and two validation sets. (D-F) Correlation heatmaps of the pyradiomics score and clinical factors in the three cohorts. SD, stable disease; CR, complete response; PD, progressive disease; PR, partial response.

Evaluating classifiable accuracy of machine learning model by PRS

The area under the ROC curves of tumor size and the pyradiomics signature were analyzed. We used Youden’s index (defining as sum of sensitivity and specificity minus 1) to calculate the optimal cut-off value (-0.14) in the ROC analysis. The AUCs showed that the tumor size (AUC=0.752, 95% CI: 66.89-83.60, P<0.001) and pyradiomics signature (AUC=0.916, 95% CI: 86.93-96.39, P<0.001) could be the predictors of treatment response to TACE in patients with HCC (Figure 3A) (Supplemental Table S1). Furthermore, the model of the pyradiomics signature in the validation sets 1 and 2 also demonstrated a high AUC for predicting treatment response, with the AUC consistent with that in the discovery set (AUC=0.902, 95% CI: 84.81-95.59, P<0.001; AUC=0.911, 95% CI: 83.26-98.98, P<0.001, respectively). The AUCs of the tumor size were 0.778 (95% CI: 52.00-82.60, P<0.001) and 0.690 (95% CI: 54.66-83.33, P=0.024, Figures 3B, C). Based on boostrap (n=2000) analysis, the efficiency of the pyradiomics signature was significantly higher than that of the tumor size in the discovery (P=0.016) and the two validation sets (P=0.016 and P=0.008). The AIC of the model comprising a combination of the pyradiomics signature and tumor size was not superior to that of the model comprising only the pyradiomics signature in the discovery set (240.53 vs. 238.94). The examples of two patients with response (RS1) or no response (RS2) are shown in our study (Figure 3D). The patient 1 had higher pyradiomics score than patient 2 did (5.44 vs.-3.10).

FIGURE 3

Figure 3 Pyradiomics signature and tumor size predict treatment response to TACE. (A-C) ROC curves show the predictive performance of pyradiomics signature and tumor size for estimating CR and PR. The bootstrap (n=2000) test results of the two ROC curves indicate that the AUC of the sum of pyradiomics signature is significantly higher than those of the tumor size in the discovery and in validation 1 and 2 sets. (D) Examples of two patients with response (RS1) or no response (RS2). ROC, receiver operating characteristic; AUC, area under curve; CR, complete response; PR, partial response.

Clinical utility and subgroup analysis of PRS predictive accuracy

The DCA of the pyradiomics signature showed relatively good performance of the model regarding clinical application (Figure 4A). It was suggested from the DCA curve that when the threshold probability in a patient was 42%, more initial response could be achieved through a pyradiomics signature than either treat-all or treat-none strategies. The probability of acquiring treatment response ranged from 8 to 100%. Thus, a pyradiomics signature accurately identifies the patients who have the initial response and may benefit from TACE therapy. The patients were divided into subgroups based on six clinical variables to estimate the classification performance of the pyradiomics model further (Supplemental Table S2). The AUC was higher in the female patients than that in the male patients (0.966, 95% CI: 92.18-100.00 vs. 0.8985, 95% CI: 86.12-93.59, P=0.021) (Figure 4B). Meanwhile, the predictive accuracy was not affected by age, Child–Pugh classification, and AFP levels, compared with bootstrap=2000 (P=0.354, P=0.998, and P=0.424, respectively) (Figures 4C-E). Subgroup analysis by tumor number and tumor size also showed no significant difference in AUCs (P=0.443 and P=0.597, respectively) (Figures 4F-G).

FIGURE 4

Figure 4 DCA and subgroup analysis of the pyradiomics signature for predictive treatment response to TACE. (A) DCA of the pyradiomics signature. (B-G) Pyradiomics signature of clinical subgroup predicts therapy response in all patients (n=313) undergoing TACE treatment. The two ROC curves are compared using the bootstrap (n=2000) test. DCA, decision curve analysis; ROC, receiver operating characteristic; TACE, transarterial chemoembolization.

The prognostic value of PRS in patients undergoing TACE treatment

According to the optimal cut-off value of PRS, we divided the patients into two groups. The patients with high values (>-0.14) defined as RS1, and the patients with low values (≤-0.14) defined as RS2. In our study, 112 patients from Nanfang Hospital and 29 patients from the Second Affiliated Hospital of Gui Zhou Medical University had the information of prognosis. Therefore, based on stratification of RS1and RS2, PFS and OS were analyzed in a total of 141 patients. We found the patients with RS1 had a longer PFS and OS than those with RS2 (Median PFS: 25 vs. 9 months, hazard ratio [HR]=2.78, 95% CI: 1.93-4.00, P<0.001; Median OS: 49 vs. 22 months, HR=2.23, 95% CI: 1.47-3.38, P<0.001, respectively) (Figures 5A, B). Moreover, we used the biomarkers of RS1 and RS2 to predict the 1-, 3-, and 5-year-OS and PFS rates. In the time-dependent ROC curve estimation from censored survival, the higher accuracies were 1-year-OS/PFS predictions than 3- and 5-year-OS/PFS predictions (0.679 vs.0.676 vs.0.662; 0.737 vs. 0.602 vs.0.617, respectively) (Figures 5C, D).

FIGURE 5

Figure 5 Prognosis prediction of PRS in patients undergoing TACE treatment. (A, B) The OS and PFS of two classification (RS1 vs. RS2) are compared in the patients. (C, D) Time-dependent ROC curve analysis is performed in 1-year-OS/PFS, 3-year-OS/PFS and 5-year-OS/PFS predictions. PFS, progression-free survival; OS, overall survival; PRS, pyradiomics score; TACE, transarterial chemoembolization; ROC, receiver operating characteristic.

Association between PRS and glycolysis pathway genes stratify prognosis in TCGA-HCC patients (No-TACE treatment)

The same algorithm of machine learning in the training cohort was used in the CT images of 46 patients with HCC from the TCGA database (TCGA-HCC cohort) and the PRS was calculated by the pyradiomics formulation. According to the above cut-off value of pyradiomics score, all patients were then divided into two groups, which were defined as RS1 group (n=23) and RS2 group (n=23). In this experiment, high pyradiomics scores were mostly observed in the RS2 group and low pyradiomics scores were contrary. Through differential gene expression analysis (RS1 vs. RS2), between-group comparisons showed that 151 genes were significantly downregulated, while 167 genes were upregulated (FDR adjust P<0.05) (Figure 6A) (Supplemental Figure S3). The GSEA found several KEGG pathways were significantly associated with PRS, such as small molecule catabolic process, glycolysis, and recycling of bile acids and salts. In many pathways, we speculated that glycolysis resulting from tumor hypoxia status were mostly associated with TACE treatment. Therefore, we mainly focused on the glycolysis pathway genes and found that RS1 group was negatively associated with HK2 and PFKP presence (Figure 6B) (Supplemental Figure S4). We cound not determine the ability of the PRS in this study to evaluate clinical prognosis in patients with HCC who underwent hepatectomy but did not receive TACE therapy. We subsequently found that RS1 group had significantly longer DFS and OS time than the RS2 group did (P<0.001 and P=0.008, respectively) (Figures 6C, D).

FIGURE 6

Figure 6 Different genes, signaling pathways, and prognosis associated with predicted treatment response. (A) The heat map plot of mRNA expression from different genes. Patients with HCC are shown on the y-axis, and gene expression is shown on the x-axis. (B) Some genes are enriched in certain signaling pathways, such as in the glycolysis and gluconeogenesis pathways. (C, D) The responder group shows better DFS and OS than the non-responder group does. HCC, hepatocellular carcinoma; DFS, disease-free survival; OS, overall survival.

Discussion

Based on preoperatively CT images, this study developed a pyradiomics signature to accurately predict the initial treatment response and prognosis in HCC patients who underwent TACE therapy. The DCA and subgroup prediction analysis showed the good clinical performance of the model. We further analyzed the association between the pyradiomics model and KEGG pathway genes in TCGA-HCC database. Our findings provided a novel insight into the interpretability of the machine learning model for predicting therapy response or prognosis in different types of tumors.

TACE is a standard treatment modality for HCC patients with BCLC stage B disease (7). Although recent studies reported that there was a non-superiority of TACE with respect to bland embolization, and direct incremental costs of drug-eluting beads TACE (DEB-TACE) can be acceptable in hepatocarcinoma patients (28, 29), considering the extensiveness of method and health policy in our country, we chose the conventional TACE treatment in this study. Treatment response to the first TACE is a well-known predictor of clinical outcomes in patients with middle-stage HCC. This study used pyradiomics algorithms from the python package to extract features of tumor shape, texture, intensity, and wavelet transform characteristics from three-dimensional CT images in HCC patients who underwent TACE as reported by previous studies (30, 31). In our study, we used LASSO of 5-fold cross-validation to chose pyradiomics features for CT images, and 24 features were finally selected. We then used logistic regression analysis to calculate the pyradiomics score among HCC patients who underwent TACE. The results showed significantly higher pyradiomics scores in the response group. The response group also showed smaller tumor size, which was consistent with previous reports (32, 33). There was a significantly negative correlation between the pyradiomics score and tumor size. Compared with tumor size, the pyradiomics model showed higher accuracies for predicting treatment response. This could be because the pyradiomics model included larger amount of radiology information than the model based on tumor size did (34), and a machine learning method, such as LASSO, made a great contribution to enhance the predictive ability of diagnosis, therapy response or prognosis in solid tumors. However, the model combining the pyradiomics signature and tumor size did not show superior AIC to that using the pyradiomics signature alone. The machine learning model had a robust accuracy for predicting initial treatment response to TACE and the method could be used as an invasive stool in the other cancers.

The DCA revealed that the pyradiomics model could predict the treatment response when the probability of acquiring treatment response ranged from 8 to 100%. This result indicated that the pyradiomics signature could help determine clinical strategy and identify the patients who had initial treatment response according to the above the optimal cut-off value (>-0.14). Additionally, some clinical variables may affect the predictive accuracy of the pyradiomics model (6, 35), and the model’s predictive accuracy in patient subgroups was rarely reported and was still unclear. Interestingly, in this study, our subgroup analysis (Supplemental Table S2) showed that our machine learning model had a better predictive accuracy in females than in male patients. Of note, the sample size of female patients was significantly smaller than that of male patients (41 vs. 272). The different result of TACE treatment in male or female could be further confirmed by large patients. Meanwhile, subgroup analysis by age showed no significant difference in accuracies between patients aged ≤60 years and >60 years. There were also no significant differences in accuracy according to the Child–Pugh classification, AFP, tumor size, and number of tumors. Collectively, these results support that our pyradiomics model can accurately predict individual treatment response to TACE and may help identify patients who will benefit from the treatment. In addition, the correlation between predictive classification of pyradiomics model and prognosis was investigated. We found the RS1 group indicating well response of initial TACE treatment had better OS and PFS than RS2 group did. This result further demonstrated the initial response could improve the OS and PFS in the patients receiving TACE therapy. Moreover, our model based on the series CT images and analysis of time-dependent ROC curve could preoperatively predict the prognosis and to screening of patients who could benefit from TACE treatment among the BCLC stage B patients with HCC.

To explore the potential mechanism of the pyradiomics model, we used CT images from the TCGA-HCC cohort to screen the different genes and related cancer signaling pathways (Supplemental Figure S3). RS1 group was significantly associated with hypoxia. Previous studies have reported that the tumor hypoxia status is associated with resistance to chemotherapy, targeted therapy, and radiation therapy (36–38). This indicates that our machine learning model predicts treatment response to TACE by characterizing the hypoxia change in tumors from CT images. In the TCGA-HCC patients (No-TACE treatment), RS2 group with severe hypoxia also showed a significantly poorer prognosis than the RS1 group, suggesting the hypoxic status of HCC may be associated with resistant mechanism to TACE treatment, which is consistent with previous findings (39–41). Improving the hypoxic status in tumors is potentially one of the methods to improve the therapy effect and this mechanism should be analyzed in more detail.

Our study has some limitations. First, the number of patients with middle-stage HCC was relatively small. Second, the study was conducted retrospectively. However, we used CT images from four centers, and all the CT images were normalized before extracting pyradiomics features to develop a robust predictive model. Multi-center model analysis also showed robust predictive performance. However, the model still needs to be validated in larger prospective studies. Third, we trained and validated all the pyradiomics signatures in four medical centers. Extracting features from ROI images may have issues in reproducibility. In this study, we evaluated the reproducibility of pyradiomics and found a good agreement (ICC>0.75), markedly improving the predictive model’s robustness. Future studies should develop an automatic segmentation model for liver tumors and minimize the discrepancies among pyradiomics features.

In conclusion, the machine learning model based on pyradiomics features from 3D-CT images is a noninvasive yet highly accurate model for predicting the initial response and prognosis to TACE in patients with HCC. Thus, it may be a feasible tool for identifying patients who will benefit from TACE. The association between the pyradiomics model and cancer-related signaling pathways might help clinicians further understand the internal mechanism of machine learning. Finally, this radiology method could be used to improve the accuracy in clinical decision-making for other types of malignant tumors.

Data availability statement

For study involving human sample: The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author. For study involving TGCA data: The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

The studies involving human participants were reviewed and approved by Nanfang Hospital, the Second Affiliated Hospital of Gui Zhou Medical University, and Sun Yat-sen University Cancer Center. The patients/participants provided their written informed consent to participate in this study.

Author contributions

(I)Conception and design: JP. (II) Administrative support: None. (III) Provision of study materials or patients: JH and JP. (IV) Collection and assembly of data: JP, JZ, and FL. (V) Data analysis and interpretation: JP and JH. (VI) Manuscript writing: All authors and (VII) Final approval of manuscript: All authors.

Funding

This work was supported by the Science and Technology Fund Project of Guizhou Provincial Health Commission [gzwjkj2019-1-077], Qian Dong Nan Science and Technology Program [qdnkhJz2019-026], Open Funds of State Key Laboratory of Oncology in South China [HN2020-02], Qian Dong Nan Science and Technology Program [qdnkhJz2020-013], the Guizhou Medical University 2018 academic new talent cultivation and innovation exploration project [Grant No. 20185579-X], Science and Technology Foundation of Guizhou Province [Grant No. Qian ke he ji chu-ZK 2021, yi ban 454], and National Nature Science Foundation of China (Grant Nos. 82060327).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.853254/full#supplementary-material

Abbreviations

AFP, alpha-fetoprotein; AIC, Akaike information criterion; AUC, area under the curve; BCLC, Barcelona Clinic Liver Cancer; CI, confidence interval; CR, complete response; CT, computed tomography; DFS, disease-free survival; HCC, Hepatocellular carcinoma; HR, hazard ratio; ICCs, inter-correlation coefficients; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; OS, overall survival; PD, progressive disease; PFS, progression-free survival; PR, partial response; PRS, pyradiomics score; ROC, receiver operating characteristic; ROI, regions of interest; SD, stable disease; TACE, trans-arterial chemoembolization; TCGA, The Cancer Genome Atlas; NCCN, National Comprehensive Cancer Network; MRI, magnetic resonance imaging.

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2021) 71:209–41. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Shen JY, Li C, Wen TF, Yan LN, Li B, Wang WT, et al. Liver transplantation versus surgical resection for HCC meeting the Milan criteria: A propensity score analysis. Med (Baltimore) (2016) 95:e5756. doi: 10.1097/MD.0000000000005756

CrossRef Full Text | Google Scholar

3. Forner A, Reig M, Bruix J. Hepatocellular carcinoma. Lancet (2018) 391:1301–14. doi: 10.1016/S0140-6736(18)30010-2

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Pinna AD, Yang T, Mazzaferro V, De Carlis L, Zhou J, Roayaie S, et al. Liver transplantation and hepatic resection can achieve cure for hepatocellular carcinoma. Ann Surg (2018) 268:868–75. doi: 10.1097/SLA.0000000000002889

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Chai NX, Chapiro J. Therapy of intermediate-stage hepatocellular carcinoma: Current evidence and clinical practice. Semin Interv Radiol (2020) 37:456–65. doi: 10.1055/s-0040-1719186

CrossRef Full Text | Google Scholar

6. Chen M, Cao J, Hu J, Topatana W, Li S, Juengpanich S, et al. Clinical-radiomic analysis for pretreatment prediction of objective response to first transarterial chemoembolization in hepatocellular carcinoma. Liver Cancer (2021) 10:38–51. doi: 10.1159/000512028

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Prince DA-O, Liu K, Xu W, Chen M, Sun JY, Lu XJ, et al. Management of patients with intermediate stage hepatocellular carcinoma. Ther Adv Med Oncol (2020) 12:1758835920970840. doi: 10.1177/1758835920970840

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Kim BK, Kim SU, Kim KA, Chung YE, Kim MJ, Park MS, et al. Complete response at first chemoembolization is still the most robust predictor for favorable outcome in hepatocellular carcinoma. J Hepatol (2015) 62:1304–10. doi: 10.1016/j.jhep.2015.01.022

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Park C, Chu HH, Kim JH, Kim SY, Alrashidi I, Gwon DI, et al. Clinical significance of the initial and best responses after chemoembolization in the treatment of intermediate-stage hepatocellular carcinoma with preserved liver function. J Vasc Interv Radiol (2020) 31:1998–2006.e1. doi: 10.1016/j.jvir.2020.04.017

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Zhang Y, Zhang M, Chen M, Mei J, Xu L, Guo R, et al. Association of sustained response duration with survival after conventional transarterial chemoembolization in patients with hepatocellular carcinoma. JAMA Netw Open (2018) 1:e183213. doi: 10.1001/jamanetworkopen.2018.3213

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Bae JS, Kim JH, Lee DH, Kim JH, Han JK. Hepatobiliary phase of gadoxetic acid-enhanced MRI in patients with HCC: Prognostic features before resection, ablation, or TACE. Eur Radiol (2021) 31:3627–37. doi: 10.1007/s00330-020-07499-w

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Nam JY, Choe AR, Sinn DH, Lee JH, Kim HY, Yu SJ, et al. A differential risk assessment and decision model for transarterial chemoembolization in hepatocellular carcinoma based on hepatic function. BMC Cancer (2020) 20:504. doi: 10.1186/s12885-020-06975-2

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Vesselle G, Quirier-Leleu C, Velasco S, Charier F, Silvain C, Boucebci S, et al. Predictive factors for complete response of chemoembolization with drug-eluting beads (DEB-TACE) for hepatocellular carcinoma. Eur Radiol (2016) 26:1640–8. doi: 10.1007/s00330-015-3982-y

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Zhang WA-O, Xu AH, Wang W, Wu YH, Sun QL, Shu C. Radiological appearance of hepatocellular carcinoma predicts the response to trans-arterial chemoembolization in patients undergoing liver transplantation. BMC Cancer (2019) 19:1041. doi: 10.1186/s12885-019-6265-1

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Huang Y, Zeng H, Chen L, Luo Y, Ma X, Zhao Y. Exploration of an integrative prognostic model of radiogenomics features with underlying gene expression patterns in clear cell renal cell carcinoma. Front Oncol (2021) 11:640881. doi: 10.3389/fonc.2021.640881

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Shin Y, Nam Y, Shin T, Choi JW, Lee JH, Jung DE, et al. Brain MRI radiomics analysis may predict poor psychomotor outcome in preterm neonates. Eur Radiol (2021) 31:6147–55. doi: 10.1007/s00330-021-07836-7

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Song Z, Tang Z, Liu H, Guo D, Cai J, Zhou Z. A clinical-radiomics nomogram may provide a personalized 90-day functional outcome assessment for spontaneous intracerebral hemorrhage. Eur Radiol (2021) 31:4949–59. doi: 10.1007/s00330-021-07828-7

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Priya S, Agarwal A, Ward C, Locke T, Monga V, Bathla G. Survival prediction in glioblastoma on post-contrast magnetic resonance imaging using filtration based first-order texture analysis: Comparison of multiple machine learning models. Neuroradiol J (2021) 34:355–62. doi: 10.1177/1971400921990766

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Liu C, Gong J, Yu H, Liu Q, Wang S, Wang J. A CT-based radiomics approach to predict nivolumab response in advanced non-small-cell lung cancer. Front Oncol (2021) 11:544339. doi: 10.3389/fonc.2021.544339

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Chen BT, Jin T, Ye N, Mambetsariev I, Wang T, Wong CW, et al. Predicting survival duration with MRI radiomics of brain metastases from non-small cell lung cancer. Front Oncol (2021) 11:621088. doi: 10.3389/fonc.2021.621088

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Carles M, Fechter T, Radicioni G, Schimek-Jasch T, Adebahr S, Zamboglou C, et al. FDG-PET radiomics for response monitoring in non-small-cell lung cancer treated with radiation therapy. Cancers (Basel) (2021) 13:814. doi: 10.3390/cancers13040814

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Peng J, Zhang J, Zhang Q, Xu Y, Zhou J, Liu L. A radiomics nomogram for preoperative prediction of microvascular invasion risk in hepatitis b virus-related hepatocellular carcinoma. Diagn Interv Radiol (2018) 24:121–7. doi: 10.5152/dir.2018.17467

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Han L, Zhao K, Li Y, Han H, Zhou L, Ma P, et al. A gut microbiota score predicting acute graft-versus-host disease following myeloablative allogeneic hematopoietic stem cell transplantation. Am J Transplant (2020) 20:1014–27. doi: 10.1111/ajt.15654

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Tovoli F, Renzulli M, Negrini G, Brocchi S, Ferrarini A, Andreone A, et al. Inter-operator variability and source of errors in tumour response assessment for hepatocellular carcinoma treated with sorafenib. Eur Radiol (2018) 28:3611–20. doi: 10.1007/s00330-018-5393-3

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Cui S, Tang T, Su Q, Wang Y, Shu Z, Yang W, et al. Radiomic nomogram based on MRI to predict grade of branching type intraductal papillary mucinous neoplasms of the pancreas: A multicenter study. Cancer Imaging (2021) 21:26. doi: 10.1186/s40644-021-00395-6

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Wang T, Gao T, Guo H, Wang Y, Zhou X, Tian J, et al. Preoperative prediction of parametrial invasion in early-stage cervical cancer with MRI-based radiomics nomogram. Eur Radiol (2020) 30:3585–93. doi: 10.1007/s00330-019-06655-1

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Zhang W, Fang M, Dong D, Wang X, Ke X, Zhang L, et al. Development and validation of a CT-based radiomic nomogram for preoperative prediction of early recurrence in advanced gastric cancer. Radiother Oncol (2020) 145:13–20. doi: 10.1016/j.radonc.2019.11.023

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Facciorusso A, Bellanti F, Villani R, Salvatore V, Muscatiello N, Piscaglia F, et al. Transarterial chemoembolization vs bland embolization in hepatocellular carcinoma: A meta-analysis of randomized trials. United Eur Gastroenterol J (2017) 5:1196. doi: 10.3389/fonc.2020.01196

CrossRef Full Text | Google Scholar

29. Cucchetti A, Trevisani F, Cappelli A, Mosconi C, Renzulli M, Pinna AD, et al. Cost-effectiveness of doxorubicin-eluting beads versus conventional trans-arterial chemo-embolization for hepatocellular carcinoma. Dig Liver Dis (2016) 48:798–805. doi: 10.1016/j.dld.2016.03.031

PubMed Abstract | CrossRef Full Text | Google Scholar

30. He B, Song Y, Wang L, Wang T, She Y, Hou L, et al. A machine learning-based prediction of the micropapillary/solid growth pattern in invasive lung adenocarcinoma with radiomics. Transl Lung Cancer Res (2021) 10:955–64. doi: 10.21037/tlcr-21-44

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Van Griethuysen JJM, Fedorov A, Parmar C, Hosny A, Aucoin N, Narayan V, et al. Computational radiomics system to decode the radiographic phenotype. Cancer Res (2017) 77:e104–7. doi: 10.1158/0008-5472.CAN-17-0339

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Kim J, Choi SJ, Lee SH, Lee HY, Park H. Predicting survival using pretreatment CT for patients with hepatocellular carcinoma treated with transarterial chemoembolization: Comparison of models using radiomics. AJR Am J Roentgenol (2018) 211:1026–34. doi: 10.2214/AJR.18.19507

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Sheen H, Kim JS, Lee JK, Choi SY, Baek SY, Kim JY. A radiomics nomogram for predicting transcatheter arterial chemoembolization refractoriness of hepatocellular carcinoma without extrahepatic metastasis or macrovascular invasion. Abdom Radiol (NY) (2021) 46:2839–49. doi: 10.1007/s00261-020-02884-x

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Aerts HJ, Velazquez ER, Leijenaar RT, Parmar C, Grossmann P, Carvalho S, et al. Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach. Nat Commun (2014) 5:4006. doi: 10.1038/ncomms5006

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Jin Z, Chen L, Zhong B, Zhou H, Zhu H, Zhou H, et al. Machine-learning analysis of contrast-enhanced computed tomography radiomics predicts patients with hepatocellular carcinoma who are unsuitable for initial transarterial chemoembolization monotherapy: A multicenter study. Transl Oncol (2021) 14:101034. doi: 10.1016/j.tranon.2021.101034

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Chen J, Wang X, Yuan Y, Chen H, Zhang L, Xiao H, et al. Exploiting the acquired vulnerability of cisplatin-resistant tumors with a hypoxia-amplifying DNA repair-inhibiting (HYDRI) nanomedicine. Sci Adv (2021) 7:eabc5267. doi: 10.1126/sciadv.abc5267

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Codony VL, Tavassoli M. Hypoxia-induced therapy resistance: Available hypoxia-targeting strategies and current advances in head and neck cancer. Transl Oncol (2021) 14:101017. doi: 10.1016/j.tranon.2021.101017

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Nilsson MB, Robichaux J, Herynk MH, Cascone T, Le X, Elamin Y, et al. Altered regulation of HIF-1α in naive- and drug-resistant EGFR-mutant NSCLC: Implications for a vascular endothelial growth factor-dependent phenotype. J Thorac Oncol (2021) 16:439–51. doi: 10.1016/j.jtho.2020.11.022

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Shou Y, Yang L, Yang Y, Zhu X, Li F, Xu J. Determination of hypoxia signature to predict prognosis and the tumor immune microenvironment in melanoma. Mol Omics (2021) 17:307–16. doi: 10.1039/d0mo00159g

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zhang Q, Qiao L, Liao J, Liu Q, Liu P, Liu L. A novel hypoxia gene signature indicates prognosis and immune microenvironments characters in patients with hepatocellular carcinoma. J Cell Mol Med (2021) 24:3772–84. doi: 10.1111/jcmm.16249

CrossRef Full Text | Google Scholar

41. Zheng W, Zhang S, Guo H, Chen X, Huang Z, Jiang S, et al. Multi-omics analysis of tumor angiogenesis characteristics and potential epigenetic regulation mechanisms in renal clear cell carcinoma. Cell Commun Signal (2021) 19:39. doi: 10.1186/s12964-021-00728-9

PubMed Abstract | CrossRef Full Text | Google Scholar