Magnetic resonance imaging and molecular features associated with tumor-infiltrating lymphocytes in breast cancer

Study design

We carried out this institutional review board-approved, Health Insurance Portability and Accountability Act (HIPAA)-compliant retrospective study in three steps (Fig. 1). First, we characterized both tumor and parenchymal enhancement patterns at dynamic contrast-enhanced (DCE) MRI and evaluated their association with TILs. Second, we built a composite model to predict TILs by integrating imaging with molecular and clinicopathological data. Third, we tested the prognostic significance of the TIL model in an independent cohort.

Patient cohorts

We analyzed two breast cancer cohorts from The Cancer Genome Atlas (TCGA) project [29] and the I-SPY 1 (Investigation of Serial Studies to Predict Your Therapeutic Response with Imaging And moLecular Analysis) trial [30]. For this study, the inclusion criteria for TCGA cohort were (1) pathologically proven invasive carcinomas, (2) pretreatment DCE-MRI data available, (3) H&E-stained whole-tumor tissue sections available, and (4) tumor gene expression data from RNA-sequencing (RNA-seq) and mutational data from whole-exome sequencing available. We applied similar inclusion criteria to select patients from the I-SPY 1 cohort, except that outcomes were available, but H&E-stained slides and mutational data were not required. After selection, 126 patients from TCGA and 105 patients from I-SPY 1 were eligible for the proposed study. The detailed selection procedures are shown in Additional file 1: Figure S1. Clinical and imaging data are publicly available for both cohorts from The Cancer Imaging Archive (TCIA) (www.cancerimagingarchive.net).

Evaluation of tumor-infiltrating lymphocytes

TILs were evaluated for TCGA cohort, for which detailed biospecimen collection and processing protocols have been described elsewhere [29]. In brief, the tumor sections were collected from surgical specimens and reviewed by a board-certified pathologist to confirm the presence of invasive carcinoma. The H&E-stained whole-slide tumor sections were digitally scanned and are available from the Cancer Digital Slide Archive (http://cancer.digitalslidearchive.net/).

Two pathologists (XT and XL, with 30 and 5 years of experience, respectively, in reading breast cancer tissue slides) evaluated TILs in consensus based on the recommendations from the International Immuno-Oncology Biomarker Working Group on Breast Cancer [16]. Two pathologists simultaneously reviewed the digital slides of each patient from the Cancer Digital Slide Archive, and the TILs were measured as the percentage of lymphocytes and macrophages in the area of total intratumoral stromal compartments. In addition, three discrete categories are defined, with ≤ 10%, > 10% to ≤40%, and > 40% to ≤ 90% TILs indicating tumors with no/minimal, intermediate, and high lymphocyte infiltration, respectively [16]. To assess interrater variability, we calculated the intraclass correlation coefficient (ICC) between our TIL percentage and those reported in a previous study focused on TNBC [31] for 15 overlapped cases in TCGA cohort.

Imaging protocols

The detailed imaging protocol for TCGA cohort has been reported elsewhere [27]. In brief, the scans were performed between September 1999 and June 2006 at six participating centers with a 1.5-T or 3-T GE Healthcare (Milwaukee, MI, USA), Siemens (Erlangen, Germany), or Philips (Amsterdam, The Netherlands) whole-body MRI system with a standard double-breast coil. The dynamic protocol of DCE-MRI was in accordance with the American College of Radiology guidelines, which included one precontrast and two to seven postcontrast scans (with a gadolinium-based contrast agent), in either the sagittal or axial view.

The detailed imaging protocol for the I-SPY cohort was reported elsewhere [32, 33]. To match the MRI from TCGA cohort, we focused on the scans acquired before neoadjuvant chemotherapy (i.e., baseline scans). MRI was performed through a 1.5-T GE Healthcare, Siemens, or Philips system, with a dedicated breast radiofrequency coil. The DCE-MRI protocols include one precontrast scan and two postcontrast phases with one ~ 2.5 minutes and another one ~ 7.5 minutes.

Image processing and harmonization

Given the diverse imaging protocols within the multicenter TCGA data and I-SPY 1 cohorts, we developed a pipeline to normalize the imaging data before extracting quantitative feature. First, we applied the N4 bias correction to correct for shading artifacts. Next, we standardized the temporal resolution of DCE-MRI scans in TCGA and I-SPY cohorts. In particular, for each patient, we included DCE-MRI before contrast agent administration and two postcontrast scans, with one having a 2–3-minute delay and the other having an ~ 7.5-minute delay. Third, to explicitly account for heterogeneous imaging protocols, for each individual, the voxel values of DCE-MRI were normalized by the parenchyma without contrast (i.e., the average value of interquartile voxel from parenchyma before administrating contrast). Finally, the MRI scans were resized to have an isotropic voxel resolution of 1 mm to assure consistent and meaningful computation of 3D textural features.

Tumor and background parenchyma segmentation

The detailed process used for segmentation was reported elsewhere [28, 34]. Briefly, two radiologists with 14 and 11 years of experience, respectively, in breast imaging manually delineated the 3D tumor slice-by-slice and reached consensus regarding 3D tumor contours. The ipsilateral parenchyma was segmented automatically through Fuzzy C-means clustering. The 3D parenchymal segmentation was inspected by two radiologists, and they manually revised it when necessary.

MRI feature exaction

Molecular features related to tumor-infiltrating lymphocytes

Tumor mutation burden is an important genetic factor in mediating antitumor immunity. Tumors with a higher mutation load are associated with a higher neoantigen level and thus are more immunogenic and likely to have higher immune infiltration and more TILs [39]. The cytolytic activity reflects local immune effector function and can indicate the presence of TILs. Indeed, cytolytic activity computed from the gene transcript levels of two critical immune cytolytic effectors [40], perforin (PRF1) and granzyme A (GZMA), has been shown to be closely related to immune infiltration and CD8+ T-cell activation [41, 42]. For TCGA breast cancer cohort, gene expression data from RNA-seq and mutational data from whole-exome sequencing are available in the Genomic Data Commons (https://gdc.cancer.gov/). On the basis of these data, we computed the nonsynonymous somatic mutational burden and cytolytic activity score, defined as the geometric mean of the expression of two genes: GZMA and PRF1 [40]. Similarly for the I-SPY 1 cohort, we computed the cytolytic activity score on the basis of microarray gene expression data available from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/; [GEO:GSE22226]) [43]. The ComBat algorithm [44] was implemented to harmonize the gene expression data from TCGA and I-SPY.

Association with tumor-infiltrating lymphocytes and predictive modeling

We first evaluated the Pearson linear correlation between individual imaging features and percentage of TILs in TCGA cohort. Next, we built a predictive model for TILs by combining multiple imaging features into an imaging signature. For this purpose, we used linear regression with feature selection via LASSO (least absolute shrinkage and selection operator) [45] to avoid overfitting. In addition, tenfold cross-validation was applied and repeated 100 times to minimize the selection bias. The most frequently selected imaging features (> 90%) were used to fit the final model. Further, we investigated whether combining the imaging signature with immune-related molecular features (cytolytic score and somatic mutation burden) would improve prediction accuracy for TILs by fitting a composite model via multivariate linear regression.

Performance evaluation

To evaluate the prediction models, we calculated the Pearson linear correlation between pathologist-rated and estimated percentage of TILs. In addition, patients were divided into three recognized TIL categories (low, intermediate, and high immune infiltration) [16], and pairwise classification among the three categories was evaluated. We compared the performance of the composite model with molecular features based on cytolytic score and imaging signature. In particular, the ROC analysis and AUC were used to assess the binary prediction accuracy of the models. The threshold used to separate different prediction models was defined on the basis of Youden’s J statistics [46], and the corresponding sensitivity, specificity, and accuracy were reported. Finally, we tested prognostic significance of the imaging signature as well as the composite TIL model by assessing their association with recurrence-free survival (RFS) in the entire I-SPY 1 cohort as well as in clinically relevant subgroups according to the receptor status. Because the prognostic value of TILs seems to be strongest in TNBC [11, 13], we expect that the composite model would also be prognostic within the TNBC subgroup in the I-SPY 1 cohort.

Statistical analysis

In univariate analysis, to adjust for multiple statistical testing, the Benjamini-Hochberg method was used to control the false discovery rate (FDR). The Mann-Whitney U statistic was used to assess the statistical significance of binary classification of TIL categories by comparing the prediction models with a random guess with an AUC of 0.5. The DeLong test was used to determine the 95% CIs and compute P values for the comparison of ROC curves. The Cox proportional hazards model was used to build survival models. Kaplan-Meier analysis was used to estimate survival probability. The log-rank test and concordance index were used to assess prognostic performance. All statistical tests were two-sided. P value < 0.05 and FDR < 0.2 were considered to be statistically significant. Statistical analysis was performed in R (R Foundation for Statistical Computing, Vienna, Austria).

Patient characteristics and tumor-infiltrating lymphocyte evaluation

Imaging features associated with tumor-infiltrating lymphocytes

Each of the 17 imaging features independently characterizes the cancer phenotypes, and their pairwise correlation map is shown in Additional file 1: Figure S2. Figure 2 shows the heat map of 17 imaging features for 126 patients in TCGA cohort ranked on the basis of their TILs, monotonically increasing from left to right. In the univariate analysis, 4 of 17 imaging features were significantly associated with the percentage of TILs (P < 0.05 and FDR < 0.2), as shown in Fig. 3. Among these four features, the tumor volume was positively correlated with TILs, whereas cluster shade of signal enhancement ratio (SER) map, mean SER of tumor surrounding BPE, and proportion of BPE were negatively correlated with TILs (Additional file 2: Table S2). Next, we built an imaging signature for TILs by fitting a linear model, which consisted of five imaging features: 4.4 × M1 − 3.14 × TEX2 − 2.0 × TS − BPE2 − 2.62 × BPE1 − 0.72 × BPE3 + 13.02, where M1 = tumor volume, TEX2 = cluster shade of SER map, TS-BPE2 = mean SER of tumor surrounding BPE (2 cm), BPE1 = BPE volume (percentage enhancement or PE > 20%), and BPE3 = BPE proportion (PE, > 20%). The mean and SD values of the five selected imaging features are shown in Additional file 2: Table S3. This imaging signature had a moderate linear correlation with TILs (ρ = 0.40; 95% CI, 0.24–0.54; P = 4.2E-6). Moreover, the imaging signature is able to separate three TILs categories in pairwise fashion (Fig. 4a), with prediction accuracy of 0.73, 0.71, and 0.71, respectively (Table 3). Figure 5 showed the details of three representative patients where there is good agreement between the predicted TILs from proposed imaging signatures and TIL readings by two pathologists.

	Specificity	Sensitivity	Accuracy
Low vs intermediate TIL groups
Imaging signature	0.93	0.36	0.73
Cytolytic activity	0.86	0.68	0.83
Imaging + cytolytic activity	0.82	0.68	0.85
Low vs high TIL groups
Imaging signature	0.70	0.89	0.71
Cytolytic activity	0.68	1	0.71
Imaging + cytolytic activity	0.91	0.88	0.91
Intermediate vs high TIL groups
Imaging signature	0.64	0.89	0.71
Cytolytic activity	0.84	0.44	0.74
Imaging + cytolytic activity	0.68	0.78	0.75

TIL Tumor-infiltrating lymphocyte

Relationships between imaging, molecular signatures, and tumor-infiltrating lymphocytes

Composite model for tumor-infiltrating lymphocytes

On multivariate analysis, the imaging signature and cytolytic score remained as independent predictors of TILs (P = 0.004 and P < 0.0001, respectively) after adjusting for stage, estrogen receptor/progesterone receptor/HER2 status, and mutation burden (Table 4). We retained both significant variables and refitted a composite model for predicting TILs: 5.86 × Imaging Signature + 7.78 × Cytolytic Score + 13.0. The linear correlation between the composite model and TILs was improved (ρ = 0.62; 95% CI, 0.50–0.72; P = 9.7E-15). Detailed box plots of inferred TILs from the composite model vs the original pathologists’ readings are presented in Fig. 4c.

We tested the composite model for predicting three predefined TILs categories. As shown in Fig. 6a, cytolytic score alone could not differentiate between intermediate and high TILs (AUC, 0.63; P = 0.14). By integrating imaging signature and cytolytic score, the composite model successfully separated these two TILs groups (AUC, 0.76; P = 0.01), and the improvement was statistically significant (DeLong test P = 0.039). Similar results were observed for differentiating low vs high TILs groups (AUC, 0.88 vs 0.94) (Fig. 6b). For distinguishing low and intermediate groups, there was no significant improvement using the composite model over cytolytic score (AUC, 0.77 vs 0.79) (Additional file 1: Figure S3). In addition, we performed a detailed evaluation of the proposed composite model, as in Table 3, where the composite model’s accuracy is 0.85, 0.91, and 0.75, respectively.

Clinical validation of the composite model

The previously developed composite model was used to infer TILs based on imaging and molecular data in an independent cohort from the I-SPY 1 trial. We found that hormone receptor-negative (HR−)/HER2− or TNBC had significantly higher predicted TILs than HR+/HER2− breast cancer (P = 0.049) (Additional file 1: Figure S4). With the threshold values obtained from the training cohort, we divided the patients into three groups based on the predicted TIL values. Then, we investigated the relationship between predicted TIL groups and outcomes. In TNBC, patients without recurrence had significantly higher predicted TILs than those who developed recurrence (P = 0.024) (Fig. 7a). Within the TNBC group, distinct RFS exists between the predicted no/minimal TIL group and the predicted high/intermediate TILs group (log-rank P = 0.0008) (Fig. 8a), where the group with lower TILs had significantly worse prognosis. However, predicted TIL groups were not associated with RFS in HR+/HER2− or HER2+ breast cancer (Fig. 7b and c and Fig. 8b and c, respectively).

Additionally, to validate the imaging signature for TILs, we applied it to 44 patients with TNBC who had images publicly available in the I-SPY 1 cohort. Similarly, with the threshold value obtained from the training cohort, we classified the patients into three TIL categories. A trend similar to that for the composite mode was observed, where no/minimal TILs had a significantly worse prognosis than high/intermediate TILs regarding their RFS (log-rank P = 0.042) (Fig. 8d).