This retrospective study was conducted in accordance with ethical standards and was approved by the institutional review boards of [Jiangjin Central Hospital of Chongqing (Center 1)] (Approval No: KY20241204-001) and [Chongqing General Hospital (Center 2)] (Approval No: KY S2024-058-01), with a waiver of informed consent. It analyzed patient data from Center 1 and Center 2 between June 2022 and September 2024.

Inclusion criteria (1): Pathologically confirmed NSCLC by surgical resection or biopsy, with documented Ki-67 IHC results (2); Chest CT scans obtained before biopsy or surgery. Exclusion criteria (1): Patients with poor CT image quality or incomplete clinical/imaging data (2); Patients with unclear delineation between lesions and adjacent obstructive pneumonia or atelectatic tissue (3); Patients who received any treatment or any invasive examination before CT examination.

A total of 259 patients were enrolled from Center 1 (163 males, 96 females, mean age 65.4 ± 9.7 years), randomly divided into a training set (N = 181) and a validation set (N = 78) at a ratio of 7:3. At Center 2, an independent test set was established with 112 patients (60 males, 52 females, mean age 65.5 ± 9.6 years). The overall study design and analytical pipeline are summarized in Figure 1.

IHC analysis was used to assess Ki-67 expression, defined as the percentage of positively stained cells. Under high magnification (*400), the number of Ki-67 positive tumor cells was quantified; Ki-67 expression was calculated as: (number of positive tumor cells in each zone/total number of tumor cells in each zone) ×100%. Values from five different areas were recorded and averaged. Previous studies have reported significant prognostic differences in NSCLC patients when a 40% cutoff value for Ki-67 expression is used (23–26). Furthermore, several radiomics studies aiming to predict Ki-67 levels have also adopted this 40% threshold (27, 28). Thus, this study defined high Ki-67 expression as ≥40% and low expression as <40%.

CT examination was performed using two scanners: a UCT530 40-row helical CT scanner (United Imaging Healthcare, Shanghai, China) at Center 1, and a Philips IQon Spectral 64-row helical CT scanner (Philips Medical Systems, Best, Netherlands) at Center 2. The scanning parameters were configured as follows: Center 1: Tube voltage of 120 kV, automatic tube current modulation (reference: 80 mAs), pitch factor of 1.075, rotation time of 0.6 s/rev, collimation width of 22 mm, acquisition matrix of 512×512, reconstruction slice thickness/interval of 1 mm, and lung window settings of window width 1200 HU/window center -600 HU; Center 2: Tube voltage of 120 kV, automatic tube current modulation (reference: 129 mAs), pitch factor of 1.000, rotation time of 0.5 s/rev, maintaining identical collimation width (22 mm), acquisition matrix (512×512), reconstruction parameters (1 mm thickness/interval), and lung window settings (window width 1200 HU/window center -600 HU).

Baseline clinical data of patients were systematically documented, encompassing demographic characteristics (gender, age, alcohol consumption, smoking history), tumor profiles (clinical stages, histological type), clinical symptoms (cough, sputum production, fever, breathing difficulty, chest tightness, headache, lymphadenopathy), and laboratory data (white blood cell count, lymphocyte count, platelet count, hemoglobin level, red blood cell count, creatinine level).

Two radiologists independently assessed the CT radiological features in a double-blind method, which were defined as expert-based descriptions of lesion morphology. The assessed features included tumor location, shape, maximum cross-sectional diameter, bronchial cutoff sign, lobulation, air bronchogram, spiculation, pleural tag, vacuole (or cystic component), and tumor density. Any disagreements were resolved by consensus.

The uAI Research Portal (29) (version: 20240730, https://urp.united-imaging.com/; Shanghai United Imaging Intelligence Co., Ltd., Shanghai, China), a licensed commercial platform, was employed for image segmentation and all subsequent classification modeling. This one-stop analysis system for multimodal clinical research images was developed using Python (version 3.7.3) and incorporated the widely used PyRadiomics package. Access to the platform was granted under our institutional license agreement with the vendor. Using the platform’s integrated VB-Net model, regions of interest (ROIs) were automatically segmented and reconstructed into 3D volumes of interest (VOIs). This model, previously validated for pulmonary nodules [average Dice similarity coefficient (DSC): 0.915] (30). Additionally, a three-tier arbitration mechanism was established: two radiologists with 5 and 8 years of experience, respectively, independently performed slice-by-slice corrections of the ROIs. The inter-observer DSC was 0.887. For regions with high inter-observer agreement (DSC ≥ 0.700), the contours from the more experienced radiologist (8 years) were used. In cases of inter-observer discrepancies (DSC < 0.700), a senior radiologist with over 15 years of experience arbitrated and finalized the contours.

All modeling was conducted within the uAI Research Portal using the open-source Scikit-learn library (Scikit-learn 0.23.2). To ensure a fair comparison of predictive performance across feature modalities, clinical-radiological (based on independent clinical-radiological features), radiomics (using rad-score), and deep learning (using deep-score) models were developed with a uniform support vector machine (SVM) classifier. SVM is well-suited for handling high-dimensional data and potential non-linear relationships. For the combined model integrating clinical-radiological features, rad-score, and deep-score, the SVM classifier was compared against four other algorithms: Decision Tree, Random Forest, XGBoost, and Logistic Regression. The SVM demonstrated better performance on the validation set (Supplementary Table S1) and was consequently selected. All models were optimized for hyperparameters via grid search, with the specific parameters detailed in Supplementary Table S2.

Receiver operating characteristic (ROC) curves were plotted to calculate the area under the curve (AUC), specificity, sensitivity, and accuracy. Calibration curves, along with the Brier score, assessed the agreement between model predictions and pathologic outcomes. Decision curve analysis (DCA) was applied to quantify the clinical net benefit across a range of threshold probabilities from 0% to 100%.

SHAP is an open-source python library (RPID: SCR_021362) for machine learning model interpretation. It quantitatively assesses the contribution of each feature to model outputs through Shapley value calculation. In this study, feature summary plot, beeswarm plot, decision plot and force plot were generated to provide global and local interpretations.

The Kolmogorov–Smirnov test was applied to assess data distribution normality. Normally distributed data were presented as mean ± standard deviation ( x¯ ± s), whereas non-normally distributed data were reported as median with interquartile range [M (Q25, Q75)].

A univariate analysis was performed in the train set to assess each clinical and radiological feature’s association with the outcome. Continuous variables were compared using the student’s t-test or Mann-Whitney U test, while categorical variables were compared using the Chi-square test or Fisher’s exact test. Features with a p-value < 0.05 in the univariate analysis were subsequently included in the multivariate analysis. Multivariate analysis was conducted to identify clinical and radiological features independently associated with Ki-67 expression level. Model comparison between the combined model and the individual models was performed using the integrated discrimination improvement (IDI) index to quantify the magnitude of the net improvement in predicted probabilities. Statistical analyses were executed using the uAI Research Portal and SPSS statistical software (version 26.0, https://www.ibm.com). A p-value<0.05 was considered statistically significant.

The clinical-radiological features of the three datasets are summarized in Table 1. IHC stratification showed high Ki-67 expression (≥40%) was present in 42.0% (76/181), 48.7% (38/78), and 29.5% (33/112) of cases in the training, validation, and independent test sets, respectively. Multivariate analyses identified histologic type (p<0.001), diameter (p=0.026), density (p=0.009) and gender (p=0.041) as independent clinical-radiological features predicting Ki-67 expression level (Table 2). The clinical-radiological model achieved the AUC values of 0.763 (training set), 0.803 (validation set), and 0.820 (independent test set), as detailed in Table 3.

Following feature selection, four robust radiomics features were retained to compute the rad-score (Supplementary Table S3). The high-expression group exhibited a significantly greater rad-score compared to the low-expression group across all three cohorts (training, validation, and independent test; all p<0.001), as shown in Figures 2A–C. The radiomics model achieved AUC values of 0.817 in the training set, 0.784 in the validation set, and 0.750 in the independent test set (Table 3).

The CAM module was used to visualize the attention regions that the deep learning network focused on when making classification decisions (Figure 3; a detailed explanation is provided in Supplementary Text 2). Four discriminative deep learning features were retained to compute the deep-score (Supplementary Table S4). The high-expression group exhibited a significantly greater deep-score compared to the low-expression group across all three cohorts (all p<0.001), as shown in Figures 2D–F. The AUCs of the training, validation, and independent test sets of the deep learning model were 0.912, 0.800, and 0.817, respectively (Table 3).

The combined model performance metrics revealed AUC values of 0.929 in the training set, 0.825 in the validation set, and 0.892 in the independent test set (Table 3, Figures 4A–C). In the independent test set, the combined model demonstrated significant improvements in overall predictive performance, with IDI values of 0.115, 0.288, and 0.095 over the clinical-radiological, radiomics, and deep learning models, respectively (Table 4). The DCA demonstrated that the combined model had a superior clinical net benefit than all the other models (Figures 4D–F). The calibration curves revealed that the combined model displayed strong agreement between predictions and pathology results (Figures 4G–I), with Brier scores of 0.101 (training set), 0187 (validation set) and 0.125 (independent test set), indicating excellent probabilistic calibration. Among all the features, deep-score and rad-score contributed the most to the prediction (Figures 5A–C).

For global feature analysis: the importance plot (Figure 6A) showed the relative significance ranking of features, with Deep-score being the most influential; the bee-swarm plot (Figure 6B) displayed SHAP values for each sample across features (positive values increased the probability of high Ki-67 expression, while negative values reduced it); the decision plot (Figure 6C) illustrated cumulative feature impacts on model predictions, with gray lines denoting baseline values, red trajectories mapping positive samples (high-expression group), and blue trajectories tracking negative samples (low-expression group).

For individual visualization, force plots (Figures 7A, B) illustrated the model’s decision-making process for both high- and low-expression cases. Score calculation started from the baseline value, with each feature represented as a directional force via SHAP values. Feature contributions were quantified by colored bars: bar length corresponded proportionally to the feature’s influence magnitude on the final prediction value f(x). Specifically, red bars indicated increased probability of high expression, while blue bars denoted decreased probability. These directional forces were then cumulatively aggregated to yield the overall effect.