XGBoost Classifier Based on Computed Tomography Radiomics for Prediction of Tumor-Infiltrating CD8+ T-Cells in Patients With Pancreatic Ductal Adenocarcinoma

Objectives This study constructed and validated a machine learning model to predict CD8+ tumor-infiltrating lymphocyte expression levels in patients with pancreatic ductal adenocarcinoma (PDAC) using computed tomography (CT) radiomic features. Materials and Methods In this retrospective study, 184 PDAC patients were randomly assigned to a training dataset (n =137) and validation dataset (n =47). All patients were divided into CD8+ T-high and -low groups using X-tile plots. A total of 1409 radiomics features were extracted from the segmentation of regions of interest, based on preoperative CT images of each patient. The LASSO algorithm was applied to reduce the dimensionality of the data and select features. The extreme gradient boosting classifier (XGBoost) was developed using a training set consisting of 137 consecutive patients admitted between January 2017 and December 2017. The model was validated in 47 consecutive patients admitted between January 2018 and April 2018. The performance of the XGBoost classifier was determined by its discriminative ability, calibration, and clinical usefulness. Results The cut-off value of the CD8+ T-cell level was 18.69%, as determined by the X-tile program. A Kaplan−Meier analysis indicated a correlation between higher CD8+ T-cell levels and better overall survival (p = 0.001). The XGBoost classifier showed good discrimination in the training set (area under curve [AUC], 0.75; 95% confidence interval [CI]: 0.67–0.83) and validation set (AUC, 0.67; 95% CI: 0.51–0.83). Moreover, it showed a good calibration. The sensitivity, specificity, accuracy, positive and negative predictive values were 80.65%, 60.00%, 0.69, 0.63, and 0.79, respectively, for the training set, and 80.95%, 57.69%, 0.68, 0.61, and 0.79, respectively, for the validation set. Conclusions We developed a CT-based XGBoost classifier to extrapolate the infiltration levels of CD8+ T-cells in patients with PDAC. This method could be useful in identifying potential patients who can benefit from immunotherapies.


INTRODUCTION
The microenvironment of pancreatic ductal adenocarcinoma (PDAC) is highly immunosuppressive and heterogeneous, characterized by an abundant desmoplastic stroma, inflammatory response, and neovascularization (1). Even with surgical resection, the radical resection rate is only approximately 18% (2), and the prognosis remains poor (3). Traditional chemotherapy is minimally effective, despite some recent success (4).
Tumors are a proliferation of abnormal cells that can escape immune eradication (5). The occurrence of immune escape is a key process in cancer progression. Immunotherapy, which aims to stimulate the body's immune system against tumor cells, can overcome this problem. The recent success of immunotherapy targeting immune checkpoint inhibitors (ICI), such as the programmed cell death protein 1 (PD1) and PD1 ligand (PD-L1) pathways, has shed new light on the treatment of patients with tumors (6,7). Nonetheless, treatment with these drugs has failed to show significant clinical benefit in unselected patients with PDAC, whose objective response rate to ICI therapy has been approximately 5% in previous clinical trials (8,9). Therefore, there is a clear need to develop related predictive biomarkers to identify subsets of patients who may benefit from ICI therapy. An effective ICI therapy prerequisite is a high level of CD8 + tumor-infiltrating lymphocytes (TILs) in the tumor tissues, suggesting the importance of investigating CD8 + TILs (10). Immunohistochemistry is the gold standard for evaluating CD8 + TILs. However, the clinical application of immunohistochemistry is limited by its invasiveness, time consumption, tumor heterogeneity, and unrepeatability. In recent years, liquid biopsy is a hot spot of research. As a rapid and noninvasive alternative to tissue biopsy, liquid biopsy can capture circulating leukocytes to reflect cancer immunity (11). In general, cancer immunity consists of the local immunity in the tumor microenvironment and the systemic immunity in circulating peripheral blood (12). However, it is unclear whether systemic immune response always correlates with local immune response (12,13). Takahiro Tsujikawa et al. have emphasized the utility of local immune monitoring for patient stratification, which could improve immunotherapy's success rate (14).
Computed tomography (CT) is widely used for tumor detection, staging, and treatment response monitoring in clinical practice. Recently, radiomic biomarkers have been of great interest. They may extract spatial and temporal features from images that are useful in predicting the underlying molecular mechanisms, the tumor-immune microenvironment, and clinical outcome. Studies dealing with glioma, esophagus, lung, and liver cancers have shown that several imaging features extracted by radiomics were closely related to CD8 + TIL density (15)(16)(17)(18)(19). While reports are predicting clinicopathological results from tissue sections in PDAC (20)(21)(22), so far, there are no radiomic studies revealing the immune environment in PDAC. Subtyping of the immune microenvironment in PDAC will help design personalized immunotherapy for patients with PDAC.
Thus, we aimed to develop and validate a radiomic signature of immune infiltration in PDAC using radiomic data extracted from contrast-enhanced CT images in this study, which might help us identify the novel predictors of immunotherapy efficacy.

Patients
This retrospective single-center cross-sectional study was reviewed and approved by the Biomedical Research Ethics Committee of our institution. The requirement for informed consent was waived by the Institutional Review Board. Data were obtained from consecutive patients treated for pancreatic cancer at our institution between January 2017 and April 2018 ( Figure 1).

CT Scanning
Multiphasic CT was performed with a pancreas-specific protocol using 320-slice multidetector-row CT scanners (Aquilion ONE, Canon Medical Systems, Tokyo, Japan). The details are shown in Appendix 1.

Pathological Image Analysis
All specimens were analyzed by two pathologists, one with 30 and the other with 20 years of experience in pancreatic pathology. Pathological examination and analysis were standardized as described previously (23). A CD8 antibody (DakoCytomation, Glostrup, Denmark) was used in pathological examinations. Each CD8-stained section was converted to digital pathological images by the scanner (NanoZoomer S60, Hamamatsu Healthcare, Japanese). The tumor boundaries were manually delineated, after which a customizable digital microscopy analysis platform (Visiopharm, Hørsholm, Denmark) was used to quantify CD8 in the tumor. The two pathologists examined the results, and the outcomes were determined by consensus. Subsequently, the proportion of the area of CD8 was calculated in the tumor.

Radiological Imaging Analysis
The details are shown in Appendix 1.

Radiomics Workflow
The radiomics workflow included: (1) image segmentation, (2) feature extraction, and (3) feature reduction and selection. The detailed method is shown in Figure 2.  We used the draw tool, which is available in the Editor module of 3D Slicer version 4.8.1 (open source software; https://www.slicer.org/), to delineate the tumors in multiple slices. The details are shown in Appendix 1.
To assess interobserver reliability, ROI segmentation was performed in a blinded fashion by two radiologists (readers 1 and 2, respectively). To evaluate intraobserver reliability, reader 1 repeated the feature extraction twice during a week period. This reader completed the remaining image segmentations, and the readout sessions were conducted over 2 weeks period. Assessments of interobserver and intraobserver reliability were performed by obtaining the intraclass correlation coefficient (ICC). ICC values >0.75 were selected for subsequent investigation.

Statistical Analyses
Normal distribution and variance homogeneity tests were performed on all continuous variables. Those with normal distribution were expressed as mean and standard deviation, while those with non-normal distributions were expressed as medians and ranges. We evaluated the overall survival (OS). Deaths were set as events, and deaths attributed to other causes were set as censored observations. Survival times were calculated from surgery date to the time of death or the end of follow-up (August 1, 2020). First, the optimal cut-off CD8 level was determined with the help of X-tile (25). The X-tile program divided the patients into CD8-low and CD8-high groups, according to the optimal cut-off value. Kaplan−Meier estimates were applied to graph the survival curves, and the log-rank test was performed to analyze the differences between the curves. Second, we examined the differences in all variables between the CD8-low and CD8-high groups. Student's t-test (normal distribution), Kruskal−Wallis H test (skewed distribution), and the chi-square test (categorical variables) were used to determine the intergroup statistical differences. Third, univariate regression analysis was applied to estimate the effect size between all variables and the CD8 groups. Fourth, the prediction model was constructed using an extreme gradient boosting classifier (XGBoost). XGBoost was performed using R software supplemented with the XGBoost package. The discrimination of the model was evaluated using a receiver operating characteristic (ROC) curve. The area under the curve (AUC) was calculated concurrently. The calibration of the model was assessed using the calibration curves and Hosmer−Lemeshow test. Finally, the model's clinical usefulness was tested with a decision-curve analysis (DCA) by quantifying the net benefit at different threshold probabilities. A two-tailed p-value <0.05 was considered statistically significant. All analyses were performed using R software (version 3.3.3, The R Foundation for Statistical Computing, Vienna, Austria).

Clinical Characteristics
Based on the optimal CD8 level cut-off determined by X-tile (18.69%; Figures 3A, B), all patients were divided into CD8-high (CD8 >18.69%, n = 101; 54.89%) and CD8-low (CD8 ≤18.69%, n = 83; 45.10%) groups ( Figure 3C). CD8 expression was 28.07 ± 9.12% and 14.17 ± 2.93% in the CD8high and CD8-low groups, respectively. Forty-six patients in the CD8-high group and 48 patients in the CD8-low group died. The Kaplan−Meier curves of the two groups were significantly distinct (p = 0.001). A log-rank test showed that the survival duration in the CD8-high group (22.63 months, 95% CI: 20.20-36.20) was significantly longer than that in the CD8-low group (14.67 months, 95% CI: 12.13-22.37) ( Figure  3D). Among the clinical, pathological, and imaging characteristics that we investigated, T and N stage in the training set differed significantly between the two groups. The patient characteristics are shown in Table 1.

Radiomics Analysis
A total of 1409 radiomics features were extracted from arterial and portal venous phases, respectively. The ICC interobserver and intraobserver were good, with 0.70-0.93 and 0.85-0.90, respectively.  The radiomics features were reduced and selected in the arterial and portal venous phase images. The radiomics features that did not significantly differ between the groups or did not show significant correlations with CD8 expression were excluded. The remaining 67 radiomics features were further reduced using a LASSO logistic regression model. Finally, the radiomics characteristics were reduced to 10 features (Supplemental Figures 1A, B), and the LASSO logistic regression formula was used to obtain the rad-score ( Table 2). The rad-score was significantly lower (p < 0.001) in the CD8-high group (median: -0.43; range: -1.61−1.42) than in the CD8-low group (median: -0.16; range: -1.16−2.35) (Supplemental Figure 1C).

Univariate Analysis
The results of the univariate analysis ( Table 3) demonstrated that the rad-score and T stage were significantly associated with CD8 expression.

Development, Performance, and Validation of the Prediction Model
The performance of the prediction model combining radiomics features and tumor size is shown in

Clinical Utility of the Prediction Model
The decision curve of the rad-score is shown in Figure 6. The decision curves show that with a threshold probability >0.16, using the XGBoost classifier to predict CD8 + T-cell added more benefit than the "treat all patients as high CD8 + T-cell" scheme or the "treat none as low CD8 + T-cell" scheme.

DISCUSSION
Immunotherapy has emerged as a promising treatment in cancer; assessing patients' different immune statuses with PDAC can better help physicians identify those who can benefit from immune therapies. Although relevant genetic subtypes have been identified (26,27), clinicians still lack reproducible and biologically meaningful biomarkers to identify patients with favorable prognoses at initial diagnosis. We focused on the radiomic features extracted from the pancreatic protocol CT scan, which is widely used in practice, to identify such a biomarker. Compared with histopathologic and molecular biomarkers, radiomics has the potential to predict the molecular profiles of tumors from image phenotypes inexpensively, non-invasively, and easily. In this study, we observed that the infiltration of CD8 + TILs is associated with the prognosis of patients with PDAC. Further, we established a CT-based radiomic score to extrapolate the tumor immune infiltration levels in patients with PDAC.
Cellular immunity is important for the immune system and plays a critical role in eliminating cancer and preventing inflammation. CD8 + T-cells can lyse tumor cells directly that expose tumor-specific antigens in various cancers, including PDAC (28). Quantification of CD8 + TILs, known as the immunoscore, was developed to evaluate the association between the infiltration level of CD8 + TILs and patients with PDAC survival (27,(29)(30)(31)(32), with results consistent with those in our study. Our study used X-tile plots (25), a new bioinformatics tool for biomarker assessment and outcomebased cut-point optimization, to provide a global assessment of every possible way of dividing the patients with PDAC into lowand high-level CD8 expression. All patients were divided into either CD8-high (CD8 >18.69%, n = 101; 54.89%) or CD8-low (CD8 ≤18.69%, n=83; 45.10%) groups, based on the optimal cut-off of CD8 level, as determined by x-tile (18.69%). Furthermore, a log-rank test showed that the survival duration in the CD8-high group (22.63 months, 95% CI: 20.20-36.20) was significantly longer than that in the CD8low group (14.67 months, 95% CI: 12. 13-22.37).
Compared to the immunoscore of surgical tissue samples, measuring the level of CD8 + TILs by radiomics is more convenient, which is especially important in patients with unresectable PDAC. Sun et al. built a CT-based radiomic signature to assess CD8 + TIL infiltration determined by RNAseq data (33). More than fifteen types of tumors were included in this study, but not PDAC. We are the first to have investigated the possibility of extrapolating the infiltration levels of CD8 + TILs in PDAC using radiomics based on CT in both a training and validation cohort. The rad score can reflect the infiltration level of CD8 + TIL, and the association between lower rad scores and higher CD8 + TIL infiltration can be observed, suggesting that the rad score may be an important prognosis biomarker for patients with PDAC. Furthermore, the DCA test showed that the rad-score could effectively facilitate clinical decision-making. OR, odds ratio; CI, confidence interval; Rad-score radiomics score; BMI, body mass index; LVSI, lymphvascular space invasion; PD, pancreatic duct; CBD, common bile duct; Rad-score, radiomics score.
The intra-tumor heterogeneity assessed by radiomics may reflect genomic heterogeneity, and tumors with more genomic heterogeneity are more likely to resist therapy and develop distant metastasis; thus, they tend to predict a worse prognosis (14,(34)(35)(36).
Texture analysis is an objective mathematical method based on their gray levels and spatial relationships (37). The most widely used texture analysis methods are the gray-level co-occurrence matrix (GLCM) and gray-level run-length matrix (GLRLM) (38). GLCM can describe the pixel distribution within a region and indicate the frequency of various combinations of grey values observed (38). GLRLM describes the relationships in linear onedimensional terms (39). Chen et al. observed that highly immune infiltrated HCCs were more homogenous, explaining the high value of GLCM (17). Sun et al. observed that GLRLM could be representative of inflammatory infiltrate, which could reflect homogeneity or heterogeneity of an image (33). In our study, the radiomic signature comprised textural features from the gray-level size-zone matrix (GLSZM). GLSZM is an extended version of GLRLM that describes the size and intensity of voxels clusters in a region of interest (40), which has proven useful when the main characteristic is heterogeneity (40). There are several limitations to this study. First, our validation cohort was from the same center as the training cohort, which restricts our findings' generalizability to other centers. Second, as a retrospective single-center study, the relatively small sample size may weaken our conclusion. The sample size should be increased to help draw a more reliable result. Third, a few studies have found the importance of joint analysis of PD-L1 expression with CD8 expression, which may explain the mechanism of the immunosuppressive microenvironment of PDAC (20,41,42). However, several studies (20,41) have observed that the PD-L1 -/CD8 high subtype had the best survival, whereas patients with low CD8 expression had similar survival regardless of PD-L1 status, which means the endogenous CD8 + TIL-mediated antitumor immune response may play a key role in the prognosis of patients with PDAC. Therefore, evaluating CD8 infiltration levels should be prioritized in a limited timeframe. Fourth, a few recent studies have suggested that the combination of intratumoral and peritumoral radiomics is more effective in predicting therapeutic outcomes (17,43). Therefore, in the future, further studies involving peritumoral radiomics in larger populations are needed. In addition, the prediction performance of XGBoost in this study is not fully satisfactory, so we will continue to explore other deep learning models to improve the diagnostic efficiency in the future.

CONCLUSION
In conclusion, our study established and validated an enhanced CT-based rad-score for predicting the infiltration level of CD8 + TILs in patients with PDAC. This rad-score may be useful in the pretreatment prediction of individual patient immunoscores to guide accurate prognosis prediction and precision immunotherapy for patients with PDAC.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Biomedical Research Ethics Committee of Changhai hospital. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
Guarantors of the integrity of entire study: YB and CW. Study design or data acquisition or data analysis/interpretation: all authors. Manuscript drafting or manuscript revision: JL, ZS, and   The black line shows the hypothesis that all patients had low CD8 + T-cell infiltration. The x-axis shows the threshold probability, which is where the expected benefit of treatment is equal to the expected benefit of avoiding treatment. The decision curves show that with a threshold probability greater than 0.16, using the prediction model to predict CD8 + T-cell infiltration adds more benefit than the treat-all-patients as high CD8 + T-cell infiltration scheme or the treat-none as low CD8 + T-cell infiltration scheme in the training set.
Research Plan of SHDC (SH DC20 20CR40 73), 234 Platform Discipline Consolidation Foundation Project (2019YPT001), and Shanghai Science and Technology Innovation Action Plan Medical Innovation Research Project (20Y11912500).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.671333/ full#supplementary-material Supplementary Figure 1 | Radiomic feature selection using a parametric method, the least absolute shrinkage and selection operator (LASSO). (A) Selection of the tuning parameter (l) in the LASSO model via 10-fold cross validation based on minimum criteria. Binomial deviances from the LASSO regression crossvalidation procedure have been plotted as a function of log (l). The y-axis indicates binomial deviances. The lower x-axis indicates the log (l). Numbers along the upper x-axis represent the average number of predictors. Red dots indicate the average deviance values for each model with a given l, and vertical bars through the red dots show the upper and lower values of the deviances. The vertical black lines define the optimal values of l, where the model provides its best fit to the data. The optimal l value of 0.05 with log (l) = -3.05 has been selected. (B) LASSO coefficient profiles of the 67 texture features. The dotted vertical line is plotted at the value selected using 10-fold cross-validation in A. The 10 resulting features with non-zero coefficients are indicated in the plot. (C) The error-bar chart of the 10 radiomics features and radiomics score.