This study utilized X-Ray images from multicenter. Following thorough review and approval by the respective hospital ethics committees, the retrospective dataset was granted an exemption from the need for patient informed consent. In contrast, patients included in the prospective validation set were duly informed and provided written consent by signing informed consent forms.

To ensure the robustness and generalizability of our DLR model for classifying OVFs using X-ray images, we divided our dataset into four distinct subsets: training, internal validation, external validation, and prospective validation. The training set was used to develop the model, allowing it to learn to identify patterns and features indicative of OVFs.The internal validation set was employed to fine-tune the model parameters and mitigate overfitting, providing an initial assessment of the model’s performance. The external validation set was included to test the model’s generalizability to new, unseen data from different populations or settings, crucial for evaluating its applicability in diverse clinical environments. Finally, the prospective validation set was utilized to validate the model on prospectively collected data, offering insights into its real-world performance and ensuring its reliability and applicability over time and under varied conditions.

To form the training and validation sets (both internal and external), patient data from Center I and Center II, including X-ray, CT, and MRI related to OVFs, were gathered from December 2018 to December 2022. The selection criteria for the cases were defined as follows: Inclusion criteria included patients aged 50 years or older diagnosed with OVFs (22), encompassing those without any history of trauma or with only minor trauma incidents; availability of complete Dicom datasets for X-ray, CT, and MRI examinations, conducted within a maximum two-week interval; and comprehensive clinical data availability, including gender, age, and results from Dual-energy X-ray absorptiometry (DXA) exams. The exclusion criteria were suspected fracture cases due to infection or tumors; poor image quality or presence of artifacts; and patients with unclear health status or fracture classification. From January 2023 to June 2023, an independent prospective validation set was added from Center III, following the same inclusion and exclusion criteria. Figure 1 provides a detailed schematic of the case selection process, illustrating the random assignment of cases to the training set and internal validation set in an 8:2 ratio. For further information about the case collection process, grouping, image preprocessing, feature extraction and analysis, and model development, refer to the flowcharts and DLR workflow in Figures 2 and 3 .

All cases were classified and graded using the ASTLOF, which assigns scores to OVFs based on four aspects: vertebral morphology, MRI imaging, bone density values, and clinical symptoms (23). The descriptions are as follows: ① Morphological changes (CT scan): Normal = 0 points, compression = 1 point, burst = 2 points; ② MRI examination: Normal = 0 points, high signal changes on fat-suppressed imaging = 1 point, presence of vacuum phenomenon or fluid sign within the vertebral body = 2 points; ③ Bone density: T-score > -2.5 = 0 points, -3.5 ≤ T-score ≤ -2.5 = 1 point, T-score ≤ -3.5 = 2 points; ④ Clinical presentation: No significant pain = 0 points, position-induced lumbago = 1 point, persistent pain or presence of neurological symptoms = 2 points. Based on the total score, the classification is as follows: Class 0 (total score ≤ 3 points, conservative treatment), Class 1 (total score = 4, conservative or surgical treatment), Class 2 (total score ≥ 5 points, surgical treatment). Two radiologists, A and B, with 6 and 10 years of experience, respectively, evaluated the classification of OVFs. In instances of conflicting outcomes, a conclusive decision was reached through collaborative consultation.

Data on age, gender, and T- score from DXA were retrieved from the clinical health records system. Details about the imaging devices used for X-ray, CT, and MRI, along with the associated imaging parameters, are provided in Supplementary Table S1 . In the process of X-ray imaging for the thoracic or lumbar spine, it was essential to align the X-ray beam’s central ray perpendicular to either the 7th thoracic vertebra or the 3rd lumbar vertebra. The analyses and processing in this research were conducted using the lateral view images acquired from vertebral DR examinations.

Precision in segmenting vertebral bodies is crucial for the analysis that follows. In our research, physicians carried out the segmentation manually. Initially, physician A imported the X-ray images into the ITK-SNAP software (version 3.8.0, available at http://www.itksnap.org). During this process, the edges of the fractured vertebral bodies were carefully outlined and filled by hand. Care was taken to exclude adjacent intervertebral discs, pedicles, and any surrounding tissue

The outlining of the region encompassing the fractured vertebral body was carefully executed to ensure completeness and accuracy before being saved as a mask file in ‘nii’ format (refer to Figure 4 ). A month subsequent to this, 30 patients were randomly chosen from the training dataset, and their images were re-delineated by both physician A and physician B. To assess the consistency of vertebral body delineation, both intra-observer and inter-observer reliability were measured using the intraclass correlation coefficient (ICC) among the participants.

Each image underwent Z-score normalization to minimize variations across the images. The protocol for feature extraction adhered to the standards set by the Image Biomarker Standardization Initiative, as referenced in (24). Subsequently, the extracted radiomics features were harmonized across different centers utilizing the Combat method, detailed in (25), to mitigate discrepancies in the data. Utilizing the PyTorch deep learning framework within the Python 3.6 environment, the DTL method was executed, aligning with methodologies used in prior research (26). For this study, the ResNet50 model ( Figures 5 , 6 ) was selected as the foundational model for DTL, with a meticulously adjusted learning rate to enhance performance. Since the transfer features were selected from the second-to-last layer of the model (Average-Pooling layer), we divided the model parameters into two parts: the Backbone part and the Task-specific part. The initialization of the Backbone part used the pre-trained model parameters from RadImageNet (27) and ImageNet. The task-specific segment of the model received a random parameter initialization, consistent with the cosine annealing approach to learning rate decay (28). In our study, we employed the following hyperparameters for training our model: optimizer set to ‘Adam’, a batch size of 8, and training for 30 epochs. This learning rate was dynamically adjusted based on iteration count, as detailed in Supplementary S2 . For more details, please check the Supplementary S3 .

Initially, the selection of radiomic features was based on both their repeatability and minimal redundancy, focusing on those with an intraclass correlation coefficient (ICC) of 0.8 or higher (29). For features demonstrating high repeatability, Spearman’s rank correlation coefficient was employed to evaluate the inter-feature correlations. In instances where the correlation between any two features exceeded 0.9, only one of these features was preserved. To ensure maximum representation of the features, a strategy of greedy recursive elimination was applied. This involved systematically removing the feature with the greatest redundancy at each step. In the final step, the LASSO algorithm was employed. This method shrinks regression coefficients to zero by building a penalty function denoted as λ. Through this process, stable radiomics features were selected for inclusion in the LASSO-Cox analysis. For more details about the process of feature selection and fusion, please refer to Supplementary S4 .

To prevent data leakage, model training was exclusively conducted on the training dataset. Subsequent to the selection and integration of features, various machine learning classifiers were developed using the scikit-learn library. These classifiers encompassed a range of algorithms, including Logistic Regression (LR), Extremely Randomized Trees (Extratrees), Support Vector Machine (SVM), Light Gradient Boosting Machine (LightGBM), Multilayer Perceptron (MLP), Random Forests (RF), eXtreme Gradient Boosting (XGBoost) and k-Nearest Neighbor (KNN).

To optimize model training on the training dataset, the grid search algorithm was implemented, allowing for the adjustment of commonly used parameters within each model. A comparative evaluation of the performance across various classification models was conducted. To determine the most effective parameters for these models, 5-fold cross-validation was employed, leading to the identification of optimal fused feature labels. Furthermore, the significance of different features was assessed using the SHapley Additive exPlanations (SHAP) value (accessible at https://github.com/slundberg/shap). This method quantifies the contribution of each feature to the predictive outcomes, thereby elucidating their importance.

In the realm of machine learning, a fundamental method for handling multi-class classification tasks involves the concept of decomposition. This approach entails dividing the multi-class problem into multiple binary classification tasks. In our study, we employed the widely recognized ‘One-vs-Rest’ strategy for such multi-class scenarios. As an example, under this strategy, if Class 0 is designated as positive, then Class 1 and Class 2 are considered negative. Similarly, when Class 1 or Class 2 is assigned the positive label, the other two categories are automatically labeled as negative. Consequently, a trio of ‘One-vs-Rest’ (OvR) classification models were developed and trained. To assess the efficacy of these prediction models, Receiver Operating Characteristic (ROC) curves were constructed, and various metrics including the Area under the Curve (AUC), accuracy, sensitivity, and specificity were calculated. The models’ capability to generalize was then appraised using internal and external validation sets, along with a prospective validation set. For a comprehensive evaluation of the multi-class task, this study utilized both macro-averaged and micro-averaged AUC as the metrics of choice, which are different approaches to aggregate multiple confusion matrices. The macro-averaged AUC approach involves computing the AUC for each individual class and then averaging these AUCs across all classes. On the other hand, the micro-averaged AUC method amalgamates the predictive outcomes of all classes into a singular confusion matrix, from which the collective AUC is derived. Statistical analyses in this study were conducted using the R software (version 4.0.3). To conclude, the Delong test was applied to assess the AUC differences between prediction models, with a P-value less than 0.05 deemed indicative of a statistically significant variance.

In our study, 942 patients who satisfied the inclusion criteria were enrolled, ranging in age from 50 to 97 years, with an average age of 69.34 ± 10.19 years. The cohort consisted of 678 female and 264 male patients. Based on the DXA T-score classification, the study included 16 patients with normal bone mass, 205 with reduced bone mass, and 721 with osteoporosis. Within this group, 86 patients had experienced 2 osteoporotic vertebral fractures (OVFs), and 24 had 3 OVFs, bringing the total number of vertebral fractures analyzed to 1076. The demographic and clinical characteristics of these patients are detailed in Table 1 , while Table 2 outlines their treatment categorization across three different classes. In Class 0 (1-3 points), 357 patients (72.9%) underwent conservative treatment, and 113 patients (27.1%) underwent surgery. In Class 1 (4 points), conservative treatment was administered to 159 patients (34.5%), with the remaining 302 patients (65.5%) receiving surgical treatment. Lastly, in Class 2 (5-8 points), conservative treatment was given to 37 patients (29.6%), and surgical treatment was opted for in 85 patients (70.4%).

In the analysis, the LASSO-Cox regression model was applied to reduce the dimensionality of the combined features. The process of selecting the optimal penalty coefficient (λ set at 0.0222) and the feature selection methodology are illustrated in Supplementary Figure S5 . This figure depicts how the coefficients of the features vary with changes in λ. Post the final round of feature selection, a total of 14 radiomics features and 3 DTL features were selected to form the composite features. Utilizing these features and their respective regression coefficients, the DTL_Radscore was formulated, as shown in Supplementary Figure S6 . The detailed equation used to calculate the DTL_Radscore is available in Supplementary S7 .

In the application of the LASSO-Cox regression analysis, the selected penalty coefficient (λ) was determined to be 0.0126. The procedure used for selecting features, along with a curve graph that illustrates how the coefficients of the features change with λ, is presented in Supplementary Figure S8 . Upon completion of the final feature selection process, a combination of 6 radiomics features and 6 deep transfer learning (DTL) features were maintained as the fused features. These features, along with their associated regression coefficients, were utilized to develop the DTL_Radscore, as depicted in Supplementary Figure S9 . For a comprehensive understanding of the DTL_Radscore calculation, the formula is provided in Supplementary S10 .

According to the comparison of macro-averaged AUC, accuracy, and F1-score, the LR algorithm performed the best in the fused feature models trained on RadImageNet ( Supplementary Table S11 ) and ImageNet datasets ( Supplementary Table S12 ). The validation results of the two sets of fused feature prediction models for the three-class classification task can be found in Table 3 . Based on the RadImageNet dataset, the macro-averaged AUC for the training set, internal validation set, external validation set, and prospective validation set were all higher compared to the ImageNet dataset (0.913 vs 0.831, 0.926 vs 0.826, 0.940 vs 0.844, 0.913 vs 0.872). The Delong test revealed statistically significant differences (P<0.05) in all pairwise comparisons. Figure 7 displays the ROC curves of the two sets of fused feature models in predicting OVF classifications in the prospective validation set. Based on the binary “OvR” strategy, the RadImageNet dataset model showed the most effective prediction for Class 0, with an AUC of 0.969 and accuracy of 0.863. The AUC and accuracy for predicting Class 1 were 0.945 and 0.875, respectively, while for predicting Class 2, they were 0.809 and 0.692, respectively. Figure 8 shows cases where the prediction model based on the ImageNet dataset made classification errors, while the model based on the RadImageNet dataset made correct classifications.

To evaluate the contribution of each feature to the model predictions, the SHAP values for each feature were calculated. Supplementary Figure S13 displays the features ranked by their global SHAPley additive explanation values for the three-class classification and each individual class. This helps assess their contribution to the model predictions. The SHAP decision plots in Supplementary Figure S14 provide an intuitive visualization of the workflow of the prediction model in predicting Class 0, Class 1, and Class 2 in the prospective validation set.