The image datasets were gathered from The Cancer Imaging Archive (TCIA) open-access dataset, and the respective gene-expression profiles were acquired through The Cancer Genome Atlas (TCGA). Motivated by prior studies (16–18) that indicated that TIME is correlated with the prediction of breast cancer, we created and verified the association of imaging phenotypes with TIME by using three datasets. There was zero patient overlap across the three datasets, and descriptive and clinical statistics of all three cohorts are shown below in Table 1 .

The first dataset, called the TCGA cohort, consists of the data from 335 individual cases gathered from the TCGA database, together with RNA sequencing data from cancer specimens as well as Recurrence-free survival (RFS) and overall survival (OS); however, there are no data with respect to imaging. This cohort was then randomly partitioned into a training set (~70%) and a validation set (~30%).

The second dataset, referred to as the Radiogenomic Cohort, originally consisted of 137 cases of patients who had usable DCE-MRI images of TCGA-BRCA, along with the appropriate gene expression information in the TCGA dataset. One patient with no usable gene expression data, seven without usable clinical details, and nine whose imaging was not complete were deleted from the study. The finalized dataset consisted of 120 patients, and these were allocated to training and validation sets in a ratio of 8:2.

We enrolled a validation cohort consisting of 222 breast cancer cases (from the I-SPY 1 TRIAL in the TCIA database), together with usable DCE-MRI and appropriate RFS and OS information. We eliminated 26 cases of patients who had incomplete image sequences, 10 patients without measurable neoplasms, and 31 manifesting poor image quality. The resulting dataset thus encompassed a panel of 155 breast cancer patients.

Data in the TCGA and TCIA databases are open access, and our study adhered to the data- accessibility policies and release guidance with respect to both databases, and therefore did not require approval from the local ethics committee.

As Figure 1 illustrates, the framework of our research consisted of two blocks: (i) calculation of an immunoscore based on RNA sequencing information, and (ii) development of a radiomic feature (radiomics ImmunoScore [RIS]) for non-invasive assessment of the cancer immunoscore, and evaluation of the capability of the RIS to predict survival.

We used the CIBERSORT algorithm to calculate the patient’s immunoscore as previously described, the procedure of which can be found in Supplementary Material (19–21). This method is designed to work effectively and is already proven on gene expression spectra as measured with microarrays. In the present study, scores of 22 immune cell categories were assessed using CIBERSORT on a series of breast cancer gene expression profiles according to clinical notes. The imputed scores of immune cell groups generated from CIBERSORT were deemed to be exact under a threshold of P < 0.050. A patient was therefore deemed qualified to undergo further analysis only if CIBERSORT P <0.050. The immunoscore was then created using a Least Absolute Shrinkage and Selection Operator(LASSO) Cox regression analysis.

MR images were downloaded from the TCIA database (for MR acquisition parameters, see Supplementary Material ). Two diagnostic imaging physicians (XR and WT, with 5 and 9 years of diagnostic breast MRI experience, respectively) manually displayed the primary neoplasm on MRI images with ITK-SNAP software. In order to obtain infiltrative-margin data, a ring was formed around the primary neoplasm, with the tumor margin automatically expanded outward by 2 mm and the tumor margin contracted inward by 1 mm to form a ring of 3 mm thickness. The macrovasculature, neighboring organs, and air spaces were excluded ( Figure S1 ).

Radiomics features were defined based on the PyRadiomics Python package, version 1.2.0 (22), and analyzed using the recommended set of defaults. We extracted 479 quantified features of a patient’s MR images separately in each of the regions of interest, i.e., peritumoral and intratumoral regions, and evaluated them by calculating the AUC ( Figure S2 ). The extracted features included 90 first-order features, 14 shaped features, as well as 375 second- and higher-order textural features. The four textured features we studied were based on wavelet decomposition of the grayscale co-occurrence matrix, grayscale run-length matrix, grayscale size-region matrix—as well as the neighborhood grayscale-difference matrix. Image features with various spatial scales were derived by rotating the filter parameters (2.0, 2.5, 3.0, and 3.5) between 2.0 and 3.5 using a Gaussian spatial bandpass filter (∇2G).

Inter-observer and intra-observer consistency was performed by analyzing all radiomics features extracted based on intra and interclass correlation coefficients (ICCs). Thirty patients were randomly selected, features of which were extracted by radiologists XR and WT. The same steps were also repeated for two weeks by radiologist XR. ICC > 0.8 suggested good agreement.

Recursive feature elimination was employed for selecting the most helpful prediction features among the primary dataset.

Using the training set of Radiogenomic Cohort, we built a random forest model to predict the RNA-Seq-based immunoscore [radiomics ImmunoScore (RIS)], which was selected as the optimal model by a five-fold cross-validation. This model was executed on the validation set, with an optimal threshold for the RIS using the Youden index, which optimized the total sensitivity and specificity.

Comparisons between the two groups were made by Student’s t-test for continuous variables and either Chi-squared or Fisher exact-probability tests for categorical variables. Kaplan-Meier method-based survival curves were produced and compared using log-rank tests. We used Cox proportional risk models for univariate and multivariate analyses. LASSO-Cox regression analysis was performed for constructing the immunoscore for breast cancer. A random forest classifier model was used to classify the immunoscore. Model accuracy was evaluated with the AUC. Inter-observer and intra-observer consistency was performed by ICCs. Error detection rates were computed to obtain corrected P-values in multiple comparisons. We employed R 3.4.0 and SPSS 22.0 for statistics, and bilateral P-values <0.05 were regarded as significant.

The selected protocols for the TCGA cohort patients are presented in Figure S3 . Following application of the data-screening criterion, overall survival data from 335 clinically annotated breast cancer specimens were accessible for additional analyses. Details of the patient demographics are shown in Table 1 .

Detailed clinicopathologic features for individuals in the Radiogenomic Cohort (n=120) and Validation Cohort (n=155) are presented in Table 1 . The median age (interquartile range) of the 275 patients enrolled in this study was 51.0 (44.0–59.0) years.

We employed the survminer software package on the TCGA cohort training set (235 patients) for generating the best cutoff values per immune cell fraction. A forest plot showing the correlation between every immune cell sub-population and overall survival is shown in Figure 2A . The immunoscore was modeled on the training set using LASSO-Cox regression analyses ( Figures 2B, C ) (see Supplementary Material for the formula used to calculate the immunoscore). Time-dependent ROC analyses were performed at the 2-, 3-, and 5-year time-points to study the accuracy of the prognosis of the immunoscore as a continuous variable in the training set ( Figure 2D ), and the corresponding AUC values and calibration curve are shown in Figures S4A, B . The cut-off (-0.115) derived by the survminer package was then utilized to classify the patients in the training set into high and low immunoscore groups. The results of the five-year survival analysis of different immunoscore groups, different age groups, and different pathologic stages are shown in Figure S5 . The results of our multivariate Cox regression analysis regarding the correlation between immunoscore and overall survival are depicted in Table S1 .

The identical equation was used in the validation data set of the TCGA cohort in order to verify a similar prognostic value for the constructed immunoscore model across populations. In the validation set, the prognostic precision regarding the immunoscore used as a continuous variable was also evaluated via time-dependent ROC analyses ( Figure 2E ).

Both intra- and inter-observer ICCs were greater than 0.8, indicating good reproducibility of feature extraction.

A random forest was used to construct the classification for the IS_BC in the training set of the Radiogenomic Cohort. Select the top 10 features in terms of feature importance and plot the feature relative importance histogram ( Figure S6 ). The resulting radiomics signature (RIS) consisted of 21 predictors with six marginal features and 15 intratumoral features ( Table S2 ). In the training set, the capability of the RIS to classify high IS_BC and low IS_BC showed an AUC of 0.899 (95% confidence interval [CI], 0.832–0.966) ( Figure 3A ). The radiomics signature revealed a similar accuracy in predicting the IS_BC in the validation set with an AUC of 0.815 (95% CI, 0.607–1.000) ( Figure 3B ). The RIS, however, exhibited a higher AUC value than any single radiomics feature ( Figure S7 ). In addition, the AUC of the RIS was compared with the volume and diameter of the core and infiltration zones of the validation set ( Figure S8 ). The optimal cut-off for the RIS in the training set was 0.686 as defined by the ROC curve ( Figure 3A ). Therefore, patients were classified into a low-RIS group when their RIS was <0.686, and a high-RIS group when their RIS was ≥0.686. The association between the RIS and clinicopathologic characteristics is shown in Table S3 .

The prognostic value of the RIS was then evaluated in our validation cohort. The five-year RFS and OS in the low-RIS group were 91.49% and 91.32%, respectively, and these survival indices in the high-RIS group were 84.09% and 82.94%, respectively ( Figure 4 ), indicating that the prognoses for patients who were stratified on the basis of the RIS were significantly different.

We conducted multivariate Cox regression analyses and adjusted for clinicopathologic variables. In the validation cohort, the RIS continued to be a strong and independent prognostic predictor of RFS and OS ( Table 2 ). Finally, we noted that the combined model-integrating radiomics and clinicopathologic features enhanced the prognostic accuracies of the validation cohort ( Figure S9 and Tables S4, S5 ).