Retrospective research was approved by the ethics committee of our hospital, with the informed consent requirement waived because of the retrospective nature. This study was conducted between Dec. 2017 and Dec. 2021, and included data from 110 patients diagnosed with spinal metastasis from primary BC. The patients were enrolled according to the following inclusion criteria: 1) pathological confirmation of spinal metastasis from primary BC, 2) T1-weighted contrast-enhanced (T1CE) MRI scans were performed before treatment, and 3) aged over 18 years. The exclusion criteria were: 1) lack of pathological data, 2) presence of other tumor diseases, 3) treated with phosphate drugs or chemoradiotherapy, 4) presence of vertebral compressed fractures, and 5) artifacts or diffuse spinal metastases in the MRI image. The included patients were divided into a training group and a validation group at a 2:1 ratio using stochastic stratified sampling. Figure 1 shows the process of recruiting patients, including the inclusion and exclusion criteria and the number of patients. Clinical characteristics were gathered for each patient from medical records, and included age, menopausal status, and family history. Pathological data included estrogen receptor (ER), progesterone receptor (PR), and lymph node metastatic (LNM) status. The expression status of ER, PR, HER-2, and Ki-67 was determined with standard IHC (Gnant et al., 2011). The staining of cells indicates the expression status of pathological indicators. The ER and PR expressions were deemed positive if the number of ER or PR positive-stained nuclei was greater than 1%, and the expression level of Ki-67 was considered high if its positive staining rate was greater than 14% (Goldhirsch et al., 2011). Cases with IHC staining intensity confirmed as 3+ were defined as HER-2 positive, and cases with staining intensity of 2+ required fluorescence in situ hybridization (Fehrenbacher et al., 2020) analysis to determine whether they were HER-2 positive.

The sagittal T1CE-MRI data were obtained using a Siemens 3.0T MRI device (Verio, Siemens, Germany) with a repetition time of 420 m, echo time of 9 m, flip angle of 150°, acquisition matrix with dimensions of 320 × 272, field of view of 100 × 100 mm, and thickness of 4 mm. The T1CE MRI data were acquired by intravenous injection of Gadolinium-DTPA contrast agent (0.1 mmol/kg, Omniscan, GE Healthcare). The MRI data were stored in DICOM format on the picture archiving and communication system. The ITK-Snap software (v.3.8, available for download at www.itk-snap.org) was used by a radiologist with 4 years’ working experience to segment the region of interest (ROI) along the tumor border on the MRI images. The delineated ROIs were stored in NII (Data Format Working Group, 2004) format for further analysis. Figure 2 shows examples of manually delineated ROIs, including different levels of Ki-67 (Figures 2A, 2B) and different HER-2 status (Figures 2C, 2D).

The radiomics features were all calculated using the PyRadiomics package (van Griethuysen et al., 2017). This is a comprehensive open-source platform that processes and extracts radiomics features from medical images using a large set of engineered hard-coded feature algorithms. The radiomics features are extracted in a four-step process: i) preprocessing of the images and ROIs; ii) application of enabled filters; iii) calculation of features; and iv) output of results. The two feature types are original features (first-order, shape and texture) and transformed features. The original features were calculated from the original MRI images, whereas the transformed features were calculated based on transformed MRI images obtained by applying various filters to the original images. In this study, the filter types used were the Exponential, Wavelet, Square, Squareroot, Local Binary Pattern, Logarithm, Gradient, and Laplacian of Gaussian filters. More information on the image feature extraction process can be found in the PyRadiomics documentation (https://pyradiomics.readthedocs.io/).

To assess the reliability of the radiomics features and to exclude unstable features, 30 patients’ data were randomly selected for intraclass correlation coefficient (ICC) analysis (Koo and Li, 2016). The features with an intraclass correlation coefficient of greater than 0.80 were further selected by least absolute shrinkage and selection operator (LASSO) regression with 10-fold cross-validation. The training set was used to fit the LASSO regression model, and the sparsity of features was controlled by adjusting the regularization parameter lambda during the fitting process. The coefficients of all features were obtained from the trained LASSO regression model. A larger lambda value will result in more features having a coefficient of zero, thereby reducing the complexity of the model and the risk of overfitting (Sauerbrei et al., 2007). The value of lambda was computed at the position of one standard error from the maximum AUC (area under the receiver operating characteristic (ROC) curve), then the regression coefficient was determined and the valuable features were screened.

The radiomics signature (RS) formula was calculated by integrating the final set of radiomics features and their corresponding coefficients using the glmnet package (Friedman et al., 2010) in R v.3.6. The performance of the RSs was assessed by ROC curve analysis, with the optimal cutoff values determined by the maximum Younden index (Ruopp et al., 2008) using the sklearn and matplotlib packages in Python v.3.6. The AUC values for the features were calculated based on logistic regression using the pROC package in R. Figure 3 depicts the workflow of this research, including ROI acquisition, feature extraction, feature selection, and model construction.

To identify statistically significant differences in the clinical and pathological characteristics between the high- and low-level Ki-67 and HER-2 positive/negative groups, the Chi-Square test and Fisher’s exact test were used to compare categorical variables. The normality of continuous variables was verified by the Shapiro–Wilk test. The Student’s t-test and Mann–Whitney U test were used to evaluate the continuous values. The hypothesis tests were two-sided with statistical significance set at 0.05.

Table 1 presents statistical results regarding the patients’ characteristics. Between the high-level Ki-67 and low-level Ki-67 groups, the age was found to be significantly different (p < 0.05). Between the HER-2 positive and negative groups, no significant differences were observed (p > 0.05), although the age produced p < 0.05 in the training cohort.

The most predictive radiomics features were selected from the T1CE MRI of the spinal metastasis. Figure 4 depicts the feature selection process with LASSO (Sauerbrei et al., 2007). LASSO regression determines the most valuable features by selecting the appropriate regularization parameter lambda. To predict the Ki-67 level and HER-2 status, five and two features were finally selected, respectively. Table 2 lists the prediction performance of each of these features. Two features have p-values of less than 0.05 in both the training and validation cohorts. Figure 5 shows boxplots of the selected features, describing the maximum, minimum, median, and upper/lower quartiles, as well as the outliers. A detailed explanation of each selected feature is shown in Supplementary Table S1.

The finally selected MRI features were used to build the RSs for predicting the Ki-67 level (RS-Ki-67) and HER-2 (RS-HER-2) status. The RSs were established based on the selected radiomics features weighted by the respective LASSO coefficients. The formulas for the RSs are as follows:

RS-Ki-67 = 0.6505 - wavelet-HHH_ngtdm_Contrast × 0.2005 + wavelet-LHL_firstorder_Skewness × 0.1802 + lbp-3D-m1_firstorder_InterquartileRange × 0.1750 + log-sigma-1-0-mm-3D_glcm_InverseVariance × 0.2552 + logarithm_glszm_SmallAreaEmphasis × 0.2373.

RS-HER-2 =−0.6762 + logarithm_gldm_LowGrayLevelEmphasis × 0.0960–lbp-3D-k_firstorder_Skewness × 0.0866.

Figure 6 depicts the ROC curves of the developed RSs, where the horizontal axis represents the false positive rate and the vertical axis represents the true positive rate. The AUC was used to evaluate the classification performance of the models. As listed in Table 3, RS-Ki-67 generated good prediction performance, with AUCs of 0.812 (sensitivity (SEN) = 0.667 and specificity (SPE) = 0.846) in the training group and 0.799 (SEN = 0.722 and SPE = 0.833) in the validation group. RS-HER-2 also generated good prediction performance, with AUCs of 0.796 (SEN = 0.720 and SPE = 0.776) in the training group and 0.705 (SEN = 0.733 and SPE = 0.762) in the validation group.