Explainable machine learning integrating CT-derived intratumoral heterogeneity for histological subtype prediction in pure ground-glass lung adenocarcinoma: multicenter validation and transcriptomic characterization

Introduction

The increasing detection of pure ground-glass nodules (pGGNs) has created an important need to distinguish relatively indolent from biologically higher-risk lung adenocarcinomas before surgery. pGGNs may represent adenocarcinoma in situ, minimally invasive adenocarcinoma, or invasive adenocarcinoma with different predominant growth patterns (Fu et al., 2021; Lee et al., 2019; Zuo et al., 2026c). Among invasive tumors, lepidic-predominant adenocarcinoma (LPA) generally has more favorable outcomes than acinar-, papillary-, micropapillary-, or solid-predominant disease (Kadota et al., 2014; Song et al., 2013; Hung et al., 2014). Imaging-based predictions may therefore provide useful preoperative context, but they should complement rather than replace definitive pathological assessment.

Routine CT characteristics, including tumor size, mean CT attenuation, shape, vacuole sign, and pleural indentation, provide useful clues to invasiveness and histological composition (Lee et al., 2018; Zhang et al., 2021; Aherne et al., 2018). Radiomics extends visual assessment by quantifying intensity, texture, and shape across the whole tumor, and multicenter studies have demonstrated its value for histological grading and invasive phenotypes in lung adenocarcinoma (Fan et al., 2026; Qi et al., 2026). However, whole-tumor summaries can average spatially distinct phenotypes. Habitat-based approaches address this limitation by partitioning tumors into radiologically distinct subregions and quantifying within-tumor variation, but their interpretation remains centered on subregional phenotypes rather than a unified description of spatial organization across the entire lesion (Zuo et al., 2025b). The ITH score was developed to bridge these scales by combining local radiomic variation with the global organization of connected imaging phenotypes. A higher score reflects greater fragmentation of connected imaging regions and should be interpreted as an imaging phenotype rather than a direct measurement of genomic clones or cellular diversity (Li et al., 2023).

Previous studies have associated ITH-derived imaging phenotypes with invasiveness, high-grade patterns, and histological subtype in lung adenocarcinoma, while more recent multicenter studies have extended the framework to lymphovascular invasion and prognosis (Zhang et al., 2024; Zheng et al., 2024; Zuo et al., 2025a; Zuo et al., 2025b; Zhao Y. et al., 2026; Zuo et al., 2026a; Zuo et al., 2026e). Nevertheless, three questions remain important. First, model selection must be separated rigorously from held-out evaluation to avoid optimistic bias. Second, the incremental value and segmentation reproducibility of the ITH score require direct evaluation. Third, the biological interpretation of an imaging-derived heterogeneity score remains incomplete without independent molecular characterization.

We therefore developed an explainable machine-learning model for preoperative histological subtype prediction and used SHAP to interpret the contribution and direction of the included predictors. We then performed an independent transcriptomic characterization to examine whether higher- and lower-ITH-score phenotypes were associated with coordinated molecular programs.

Materials and methods

Study participants and ethics approval

This retrospective multicenter study included patients from three institutions: The Affiliated Hospital of Southwest Medical University, Luzhou, Sichuan, China; Liuyang Traditional Chinese Medicine Hospital, Changsha, Hunan, China; and The Central Hospital of Xiangtan (The Affiliated Hospital of Hunan University), Xiangtan, Hunan, China. The study was conducted in accordance with the Declaration of Helsinki and was approved by the respective institutional ethics committees. Approval was obtained from the Medical Ethics Committee of Liuyang Traditional Chinese Medicine Hospital (Opinion No. 2022-016), the Medical Ethics Committee of Xiangtan Central Hospital (Approval No. 2023-08-18), and the Medical Ethics Committee of Southwest Medical University (Approval No. SWMU-IRB-2024-068). The Southwest Medical University ethics documentation explicitly included an application for waiver of informed consent for the retrospective analysis. All patient data used for analysis were anonymized before analysis.

The study cohort comprised 533 patients with stage IA lung adenocarcinoma presenting as pGGNs. Patients were eligible if they had a pathologically confirmed diagnosis and complete preoperative thin-section chest CT images with a reconstructed section thickness of 0.625–1.25 mm. Exclusion criteria were multifocal lung adenocarcinoma, rare histological subtypes such as colloid, fetal, or invasive mucinous adenocarcinoma, receipt of chemotherapy or radiotherapy before CT examination or surgery, and inadequate CT image quality. After eligibility assessment, the pooled multicenter cohort was stratified according to histological outcome and randomly divided in a 7:3 ratio into a training cohort (n = 373) and a held-out internal validation cohort (n = 160). The held-out cohort was not used for model selection or hyperparameter optimization and was evaluated only after the final model had been locked. The overall analytical workflow is shown in Figure 1.

CT acquisition and radiological assessment

CT scanners and acquisition parameters are summarized in Supplementary Table S1. Nonenhanced thin-section chest CT images were reviewed by two thoracic radiologists with more than 10 years of experience. Prespecified radiological variables were age, sex, tumor size (mm), mean CT attenuation (HU), margin, shape, lobulation, spiculation, vascular convergence sign, vacuole sign, and pleural indentation. Disagreements in qualitative CT assessment were resolved by consensus.

Lesion segmentation and reproducibility

Tumors were manually segmented using ITK-SNAP (version 3.9.0). The contour was placed along the tumor boundary while avoiding adjacent normal tissue. For repeat-segmentation analysis, a representative subset of 100 lesions was segmented twice by one reader and independently by a second reader. Reproducibility was assessed for both the derived ITH score and the final XGBoost-predicted probability using intraclass correlation coefficients (ICCs) with 95% confidence intervals.

ITH score calculation

The ITH score was calculated using the established ITHscore framework, which combines local radiomic phenotyping with global connected-region analysis (Li et al., 2023). For each pixel within the segmented tumor region, radiomic descriptors were extracted from a 2 × 2-pixel local neighborhood. A total of 104 local radiomics features were used, spanning first-order statistics, shape descriptors, GLCM, GLDM, GLRLM, GLSZM, and NGTDM classes; the complete feature list is provided in Supplementary Table S2. Local feature vectors were not z-standardized before clustering. K-means clustering with five clusters was applied to group pixels with similar local imaging phenotypes. Cluster labels were mapped back to tumor space, and spatially connected regions were identified using an 8-connectivity rule. The resulting connected-region pattern was summarized as the ITH score, with higher values indicating greater spatial fragmentation of clustered imaging phenotypes. Identical extraction parameters were used across centers. The implementation was based on the publicly available ITHscore Python package (https://pypi.org/project/ITHscore/).

Machine-learning model development

Twelve supervised algorithms were evaluated: logistic regression, linear discriminant analysis, classification and regression tree, random forest, gradient boosting decision tree, adaptive boosting, extreme gradient boosting (XGBoost), LightGBM, CatBoost, support vector machine, k-nearest neighbor, and multilayer perceptron. The analysis used 12 prespecified predictors: age, sex, tumor size, mean CT attenuation, margin, shape, lobulation, spiculation, vascular convergence sign, vacuole sign, pleural indentation, and the ITH score.

All model selection and hyperparameter optimization were confined to the training cohort. Stratified five-fold cross-validation was used, and GridSearchCV optimized candidate hyperparameters with ROC-AUC as the scoring metric. Continuous variables were standardized using parameters fitted within the corresponding training folds for scale-sensitive algorithms; tree-based methods used unscaled input features. Binary CT signs were encoded as 0/1 and multicategory variables were one-hot encoded within the training pipeline. No class weighting was applied. Random state was used for data splitting, cross-validation, and model initialization. XGBoost achieved the highest training-set cross-validation AUC and was selected as the primary classifier before the held-out cohort was evaluated. The predictor set, preprocessing rules, hyperparameters, and probability threshold of 0.50 were then locked, and the final model was refitted on the entire training cohort. Complete parameter specifications are provided in Supplementary Table S3.

Performance, incremental value, and center-level validation

The held-out validation cohort was evaluated once after model lock. Performance was summarized using ROC-AUC, accuracy, precision, sensitivity (recall), F1 score, specificity, negative predictive value (NPV), balanced accuracy, PR-AUC, and Brier score. Ninety-five percent confidence intervals were reported for the principal performance measures. The confusion matrix was also examined because non-LPA represented the clinically relevant positive class.

Incremental value of the ITH score was assessed using three prespecified models within the same XGBoost framework: M1, ITH score only; M2, clinical and conventional CT features; and M3, clinical and conventional CT features plus ITH score. M3 and M2 ROC-AUCs were compared using the paired DeLong test. To assess robustness to institution-level shifts, leave-one-center-out internal–external validation was performed by holding out each center in turn while using the remaining centers for model development.

Model interpretation

SHAP values were not used for feature selection, hyperparameter optimization, or model selection. SHAP analysis was performed only after the final XGBoost model had been locked and was used exclusively for post hoc model interpretation. Global SHAP summaries quantified cohort-level model attribution, whereas local dependence plots characterized the direction of contributions for the leading predictors. SHAP values were interpreted as model attributions and not as causal biological effects.

Independent transcriptomic characterization

To complement the stage I imaging analysis with molecular characterization, an independent radiogenomic analysis was performed using the publicly available NSCLC-Radiogenomics resource and its RNA-sequencing dataset GSE103584 (Zhou et al., 2018; Bakr et al., 2018). The resource contains paired thoracic imaging and molecular data from Stanford University and the VA Palo Alto Healthcare System. After matching available pathology and ITH information and restricting the analysis to pathological Tis/T1 adenocarcinoma, 43 cases were retained. This cohort was independent of the 533-patient machine-learning development and validation cohort. Cases were divided into high- and low-ITH-score groups using the cohort median ITH score; the cutoff was determined solely from the ITH distribution and without reference to transcriptomic outcomes.

Global expression structure was visualized using principal component analysis. Differential-expression analysis was performed using limma with high versus low ITH score as the coefficient of interest while adjusting for pathological T stage and patient affiliation (Ritchie et al., 2015). Genes were ranked by the moderated t statistic and analyzed by preranked gene set enrichment analysis (GSEA) using the MSigDB Hallmark collection (Subramanian et al., 2005; Liberzon et al., 2015). Positive normalized enrichment scores (NESs) indicate enrichment toward the High-ITH end of the ranked list. False-discovery rate (FDR) correction was applied for multiple testing.

Statistical analysis

Statistical analyses were performed in Python 3.6.5 and R 4.4.1. Group comparisons used appropriate parametric or categorical tests. Statistical significance for conventional hypothesis testing was defined as two-sided P < 0.05. ROC analyses used pROC, and GSEA results were interpreted using FDR-adjusted significance. Brier scores were used as an overall measure of probabilistic prediction error.

Results

Patient characteristics

A total of 533 patients were included; 392 (73.5%) had LPA and 141 (26.5%) had non-LPA histology. The training cohort included 373 patients (LPA, 276; non-LPA, 97), and the held-out validation cohort included 160 patients (LPA, 116; non-LPA, 44). Baseline characteristics were similar between the training and validation cohorts (Table 1). Compared with LPA, non-LPA tumors more frequently showed irregular shape, vacuole sign, and pleural indentation and were associated with older age, larger tumor size, higher mean CT attenuation, and higher ITH score (Table 2).

TABLE 1. Clinical, tumor-heterogeneity, and radiological characteristics in the training and validation cohorts.

Abbreviations: CT, computed tomography; ITH, intratumoral heterogeneity; LPA, lepidic-predominant adenocarcinoma.

TABLE 2. Clinical, ITH score, and radiological characteristics of LPA and non-LPA histology.

Abbreviations: CT, computed tomography; ITH, intratumoral heterogeneity; LPA, lepidic-predominant adenocarcinoma.

Training-only model selection and held-out validation

Within the training cohort, XGBoost achieved the highest mean five-fold cross-validation AUC (0.908 ± 0.021), followed by gradient boosting decision tree (0.901 ± 0.024), AdaBoost (0.893 ± 0.026), random forest (0.885 ± 0.023), and LightGBM (0.882 ± 0.025) (Supplementary Table S4). XGBoost was therefore selected as the primary classifier before the held-out validation cohort was evaluated.

In the held-out validation cohort, XGBoost achieved an AUC of 0.919 (95% CI, 0.870–0.968), accuracy of 0.850 (95% CI, 0.789–0.901), precision of 0.857 (95% CI, 0.641–0.968), sensitivity of 0.545 (95% CI, 0.337–0.753), and F1 score of 0.667 (95% CI, 0.482–0.812) (Table 3; Figures 2, 3). The confusion matrix comprised 24 true positives, 20 false negatives, 4 false positives, and 112 true negatives, corresponding to specificity 0.966, NPV 0.848, balanced accuracy 0.756, and PR-AUC 0.824. The Brier score was 0.098. GBDT showed a numerically similar held-out AUC of 0.917 (95% CI, 0.868–0.966); classifier selection was based on training-set cross-validation rather than this small held-out difference.

TABLE 3. Held-out performance of the five highest-ranking machine-learning algorithms.

Abbreviations: AUC, area under the receiver operating characteristic curve; CI, confidence interval. Brier score is an overall measure of probabilistic prediction error; lower values indicate lower error.

Incremental value and center-level robustness

The ITH score-only model achieved an AUC of 0.834 (95% CI, 0.765–0.903). The clinical/CT model achieved an AUC of 0.867 (95% CI, 0.806–0.928), whereas addition of the ITH score increased the AUC to 0.919 (95% CI, 0.870–0.968), an absolute increase of 0.052 (paired DeLong P = 0.012). PR-AUC increased from 0.731 to 0.824 and Brier score decreased from 0.126 to 0.098 (Supplementary Table S5).

In leave-one-center-out internal–external validation, AUCs were 0.882 for Center A, 0.851 for Center B, and 0.824 for Center C, with a mean of 0.852 ± 0.029. Corresponding specificities remained high (0.907–0.942), whereas sensitivities ranged from 0.417 to 0.491, indicating that discrimination generalized better than the fixed operating threshold (Supplementary Table S6).

Segmentation reproducibility

In 100 repeat-segmented lesions, the ITH score showed an intra-observer ICC of 0.91 (95% CI, 0.87–0.94) and an inter-observer ICC of 0.86 (95% CI, 0.80–0.90). XGBoost-predicted probabilities were similarly stable, with intra-observer and inter-observer ICCs of 0.94 (95% CI, 0.91–0.96) and 0.90 (95% CI, 0.85–0.93), respectively (Supplementary Table S7).

Post-hoc model interpretation

Post-hoc SHAP analysis of the locked XGBoost model showed that the ITH score was the dominant predictor, followed by tumor size and mean CT attenuation (Figure 4). Tumor heterogeneity accounted for 53.7% of total model attribution, CT radiological features for 41.4%, and clinical variables for 5.0%. Local SHAP dependence plots showed descriptive changes in attribution around ITH scores >0.55, tumor sizes >13 mm, and mean CT attenuation > −610 HU (Figure 5). These were not prespecified clinical thresholds and should not be interpreted as causal cutoffs.

Transcriptomic characterization of the ITH phenotype

The independent radiogenomic analysis included 43 pathological Tis/T1 adenocarcinomas, divided at the cohort median into high- and low-ITH-score groups. PCA showed partial overlap rather than complete separation of the two groups, indicating that ITH did not define a globally distinct transcriptomic subtype (Figures 6A,B). The top-ranked differential-expression heatmap nevertheless showed coordinated expression patterns across the high- and low-ITH-score groups (Figure 6C).

At the pathway level, Hallmark GSEA identified coordinated positive enrichment toward a high ITH score, including G2M checkpoint, E2F targets, mitotic spindle, epithelial–mesenchymal transition, interferon alpha response, and interferon gamma response (Figure 6D). Representative enrichment profiles showed G2M checkpoint (NES = 1.86, FDR <0.001), E2F targets (NES = 1.86, FDR <0.001), and epithelial–mesenchymal transition (NES = 1.55, FDR = 0.001) as prominent high-ITH-score programs (Figure 6E; Supplementary Table S8). These findings support an association between the imaging-derived ITH phenotype and coordinated proliferative, cell-cycle, phenotypic-plasticity, and interferon-related transcriptional programs.

Discussion

This multicenter study provides convergent evidence that CT-derived intratumoral heterogeneity adds information to pGGN histological subtype prediction beyond routine clinical and CT features. Training-set cross-validation identified XGBoost before held-out evaluation, and the locked model achieved an AUC of 0.919 in the validation cohort. More importantly, the contribution of the ITH score was demonstrated directly rather than inferred solely from feature attribution: adding ITH to clinical and CT features increased AUC by 0.052, increased PR-AUC by 0.093, and reduced the Brier score by 0.028.

The operating-point results also require a balanced interpretation. The model showed high specificity (0.966) and PPV (0.857) for non-LPA disease, but sensitivity was only 0.545 at the prespecified 0.50 threshold, meaning that a substantial proportion of non-LPA lesions were missed. Thus, the model should not be considered a stand-alone rule-out test or a substitute for pathological diagnosis. The leave-one-center-out analysis reinforces this point: discrimination remained reasonably stable across centers, whereas sensitivity decreased at the fixed threshold. Prospective work should therefore focus not only on AUC but also on threshold transportability and the clinical consequences of false-negative predictions.

The repeat-segmentation analysis addresses another important concern for a spatial heterogeneity metric. Both the ITH score and final predicted probability showed good-to-excellent intra- and inter-observer ICCs. These findings suggest that the principal signal was not dominated by small contour perturbations, although fully automated segmentation and prospective reproducibility studies remain desirable.

The broader relevance of spatially explicit heterogeneity quantification is supported by related work across several tumor types. Imaging-derived heterogeneity approaches have been explored in locally advanced nasopharyngeal carcinoma, colorectal cancer, breast cancer, and hepatocellular carcinoma, with reported associations involving prognosis, invasive phenotypes, treatment-related outcomes, or tumor microenvironmental features (Li et al., 2026; Zhou et al., 2026; Zuo et al., 2026d; Zuo et al., 2026b; Zhao J. et al., 2026). These cross-tumor applications support spatial heterogeneity as a general imaging concept while emphasizing that the biological meaning of an ITH score remains disease- and endpoint-specific.

The independent transcriptomic analysis provides a separate layer of biological characterization. A high ITH score was associated with coordinated enrichment of G2M checkpoint, E2F targets, mitotic spindle, EMT, and interferon-related programs. The convergence of G2M, E2F, and mitotic-spindle signals is consistent with a more proliferative and cell-cycle-active tumor state, in keeping with the central role of sustained proliferative signaling in malignant progression (Hanahan, 2022). EMT enrichment further suggests increased phenotypic plasticity and architectural remodeling, processes that permit transitions across epithelial–mesenchymal states and are closely linked to tumor progression, invasion, metastasis, and treatment resistance (Haerinck et al., 2023; Gu et al., 2023). The enrichment of interferon-alpha and interferon-gamma response programs may additionally reflect differences in the tumor–immune context, as both type I and type II interferon signaling exert context-dependent effects on antitumor immunity, immune escape, and treatment response within the tumor microenvironment (Gocher et al., 2022; Holicek et al., 2024). These findings should be interpreted as associations rather than evidence that the ITH score directly measures molecular heterogeneity or that imaging heterogeneity causally activates these pathways.

SHAP was restricted to post hoc interpretation and was not used for feature selection or model optimization. This separation allows feature attribution to improve transparency without introducing outcome-informed selection outside the training process. The observed SHAP dependence patterns remain descriptive and should not be treated as universal clinical thresholds.

Limitations

Several limitations remain. First, this was a retrospective study. Although leave-one-center-out internal–external validation was added, a fully independent prospective external cohort was not available. Scanner and reconstruction differences may still affect local radiomic measurements.

Second, the non-LPA group was heterogeneous, and the retrospective pathology reports did not permit reliable subdivision of every non-LPA tumor into acinar-, papillary-, micropapillary-, and solid-predominant categories. The binary comparison therefore does not resolve differences among individual non-LPA patterns.

Third, class imbalance and the moderate sensitivity for non-LPA limit interpretation at the prespecified 0.50 threshold. The Brier score quantified overall probabilistic prediction error, but formal calibration and decision-curve analyses were beyond the scope of this retrospective study. Prospective work should therefore evaluate probability calibration, threshold transportability, and net benefit before implementation in clinical workflows.

Finally, the transcriptomic analysis was exploratory and based on an independent public radiogenomic cohort with a modest sample size. It supports biological plausibility at the pathway level but does not establish causality or validate the ITH score as a direct molecular surrogate.

Conclusion

An explainable XGBoost model integrating CT-derived ITH with routine clinical and CT features showed favorable discrimination for differentiating LPA from non-LPA histology in pGGNs. The ITH score provided statistically significant incremental predictive information beyond clinical and CT features, demonstrated good repeat-segmentation reproducibility, and retained useful discrimination in center-level validation. Independent transcriptomic analysis further associated high ITH with cell-cycle/proliferative, EMT, and interferon-related programs. The moderate sensitivity and retrospective design support continued prospective external validation rather than immediate clinical use.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. OpenAI ChatGPT was used solely for English-language editing and document-formatting support. It was not used to generate study data, perform statistical analyses, or determine scientific conclusions. All AI-assisted text was reviewed and approved by the authors, who take full responsibility for the final manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2026.1958127/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Zhijian Huang, Sun Yat-sen University, China

Reviewed by: Xingwei Sun, Second Affiliated Hospital of Soochow University, China

Tianyang Zhong, Northwestern Polytechnical University, China

Footnotes

AUC, area under the receiver operating characteristic curve; CT, computed tomography; CV, cross-validation; DEG, differentially expressed gene; FDR, false-discovery rate; GSEA, gene set enrichment analysis; HU, Hounsfield unit; ICC, intraclass correlation coefficient; ITH, intratumoral heterogeneity; LPA, lepidic-predominant adenocarcinoma; ML, machine learning; NES, normalized enrichment score; NPV, negative predictive value; pGGN, pure ground-glass nodule; PR-AUC, area under the precision–recall curve; ROC, receiver operating characteristic; SHAP, Shapley additive explanations; VOI, volume of interest.