Computed tomography-based radiomic features combined with clinical parameters for predicting post-infectious bronchiolitis obliterans in children with adenovirus pneumonia: a retrospective study

Patients

The research received approval from the Ethics Committee of Children’s Hospital of Chongqing Medical University (No.2023–151). Informed consent was waived due to the retrospective nature of the study. Data collection for this study began on April 3, 2023.

We enrolled children with adenovirus pneumonia who were hospitalized between October 2013 and February 2020. Inclusion criteria were: (i) aged 1 month to 18 years, (ii) acute lower respiratory symptoms, lung infltration and tested positive for adenovirus DNA in nasopharyngeal secretions or by direct immunofluorescence antigen test (Zhong, Lin & Dai, 2020), and (iii) had CT images acquired during acute adenovirus pneumonia prior to treatment. Exclusion criteria were: (i) incomplete clinical information, (ii) absence of CT images during acute adenovirus pneumonia or poor-quality CT images, (iii) concurrent infection with other pathogens during adenovirus pneumonia before the diagnosis of BO, and (iv) a previous history of organ transplantation.

Two pediatric respiratory specialists (5 and 8 years of experience) independently diagnosed the enrolled patients either as BO or non-BO patients according to the diagnostic criteria. In cases of disagreement, a consensus was reached through in-depth discussions. The children were diagnosed with BO after adenovirus pneumonia on the basis of the following criteria (Colom & Teper, 2009; Aguerre et al., 2010): (i) a history of adenoviral pneumonia, (ii) persistent or recurrent wheezing, coughing, shortness of breath, dyspnea, and exercise intolerance. Extensive wheezing and moist rales were auscultated bilaterally for more than 6 weeks with little response to bronchodilators. (iii) High-resolution computed tomography (HRCT) revealed a mosaic pattern, bronchiectasis, and thickening of the bronchial wall. (iv) Pulmonary function tests indicated obstructive or mixed ventilatory dysfunction in small airways with predominantly negative results in bronchodilation tests. (v) Other chronic lung diseases that could induce wheezing and coughing, such as bronchial asthma, congenital bronchopulmonary dysplasia, tuberculosis, and cystic fibrosis, were ruled out.

The detailed patient selection process is illustrated in Fig. 1. Patient data were randomly divided into training group and testing group at a ratio of 7:3.

CT image acquisition

Chest CT data were acquired prior to treatment via LightSpeed VCT 64-slice spiral CT (GE Healthcare, Chicago, IL, USA), Brilliance iCT 256-slice spiral CT (Philips Medical System, Eindhoven, Netherlands), Optima CT660 128-slice spiral CT (GE Healthcare, Chicago, IL, USA) scanners. The scanning parameters included tube voltage ranging from 80 to 100 kV, automatically regulated tube current, pitch of 0.984:1, collimation of 0.6 mm, slice thickness and slice interval of 5.0 mm, reconstructed slice thickness of either 1.25 or 1.0 mm, and a reconstructed slice interval of either 1.25 or 1.0 mm. The data were reconstructed via a conventional algorithm.

Clinical data

Clinical data were extracted from medical records by trained research staff via a standard collection form. The collected data included demographic characteristics, duration of fever (in days), length of hospitalization (in days), peak temperature (in °C), laboratory findings and clinical treatment.

Two experienced radiologists, with 6 and 12 years of experience, independently evaluated imaging features including the number of pneumonia lobes, pulmonary consolidation lobes, and the presence or absence of pleural effusion. In case of disagreement, consensus was reached through discussion.

Image segmentation

Digital Imaging and Communications in Medicine (DICOM) images were exported from the Picture Archiving and Communication System (PACS) and regions of interest for pulmonary lesions were manually delineated using ITK-SNAP (Version 3.8) software. The lesion margins on thin-slice CT images were manually delineated by a radiologist with 6 years of experience on the axial lung window (Fig. 2), and verified by a senior radiologist with 12 years of experience.

The lesion margins were manually delineated of thin-slice CT images by a radiologist with 6 years of experience. The segmentation results were then verified by a senior radiologist with 12 years of experience.

To ensure the availability of radiomic features due to the subjective nature of manual delineation (Fornacon-Wood et al., 2020), ROIs were redelineated for 30 randomly selected patients (BO: non-BO = 12:18). The intraclass correlation coefficient (ICC) was calculated for all extracted radiomic features from the two delineated RIOs in these 30 patients, and the features with an ICC > 0.85 were included in the subsequent study.

Radiomics feature extraction and selection

We uploaded the images and corresponding ROIs to the uAI research portal (uRP (Wu et al., 2023), Version: 20220915, United Imaging Intelligence, China), applied a series of preprocessing steps to ensure consistency and quality in radiomic feature extraction: we used the formula $f\left(x\right)={\displaystyle \frac{s\left(x-{\mu }_x\right)}{{\delta }_x}}$, $x$ and $f\left(x\right)$ are the original intensity and normalized intensity, ${\mu }_x$ and ${\delta }_x$ are the mean and standard deviation of the image intensity values, and $s$ is the optimal scaling defined by scale. Images were normalized centered on the standard deviation of the mean, resampled to voxel size 1 * 1 * 1 mm³ using b-spline interpolation, intensity normalization with a fixed bin width 25 Hounsfield units (HUs), and Z score normalization to obtain a standardized normal distribution of image intensities. Radiomic features such as first-order features, shape features, and texture features were extracted via PyRadiomics embedded within the uRP. The texture features included the gray-level co-occurrence matrix (GLCM), gray-level run-length matrix (GLRLM), gray-level dependence matrix (GLDM), neighborhood gray-tone difference matrix (NGTDM), and gray-level size-zone matrix (GLSZM).

To improve model performance and reduce overfitting in subsequent data analysis processes, it is essential to select the best feature set. We employed multivariate information (MI) and minimum redundancy maximum relevance (mRMR) methods to reduce redundant features. Subsequently, we selected non-zero coefficient features using the least absolute shrinkage and selection operator (LASSO), which employs an L1 regularizer as the cost function along with 10-fold cross validations. Ultimately, we selected the most relevant predictor features associated with prognosis. The radiomics score (Radscore) for each patient used in model construction was computed by linearly multiplying the selected features (F) with their corresponding coefficients (β) (Radscore = β₀ + β₁ * F₁ + β₂ * F₂ + … + β_n * F_n).

Models construction and evaluation

We constructed radiomic models using the Radscore, clinical models based on selected clinical parameters via multivariate logistic regression analysis in the training group, and combined models by integrating radiomic and clinical models. Logistic regression (LR), support vector machine (SVM), and random forest (RF) methods were utilized for model construction.

Model performance was assessed via the area under the receiver operating characteristic curve (AUC), and corresponding metrics, such as sensitivity, accuracy, specificity, positive predictive value (PPV), negative predictive value (NPV) and F1 score, were calculated. Calibration curves were plotted to assess model calibration. Decision curve analysis (DCA) was used to evaluate net benefit at different threshold probabilities and assess clinical utility. Model comparison was conducted via the Delong test, net reclassification improvement (NRI), and integrated discrimination improvement (IDI). Finally, a nomogram based on the best-performing model was created.

Statistical analysis

Statistical analysis was conducted via R statistical software (version: 4.1.1). All the statistical tests were two-sided, and a significance level of p < 0.05 was considered as statistically significant. Quantitative characteristics were evaluated via the Mann‒Whitney U test and are presented as medians with their corresponding quartiles. Qualitative characteristics were examined via the chi-square test and reported as counts (%). Influence characteristics that exhibited statistical significance at p < 0.05 in the univariate logistic analysis were included in the multivariate analysis using a stepwise selection approach.

Patient characteristics and clinical model construction

During the study period, data on 2,116 patients diagnosed with adenovirus pneumonia were collected. Ultimately, 165 patients who satisfied the inclusion and exclusion criteria were included in the study. Among them, 70 patients belonged to the BO group while 95 patients were part of the non-BO group. The participants were divided into a training group (n = 115, BO: non-BO = 49:66) and a testing group (n = 50, BO: non-BO = 21:29) at a ratio of 7:3. The detailed characteristics of the patients can be found in Table 1.

Several clinical characteristics, including age, gender, duration of fever, length of hospitalization, hemoglobin (Hb), albumin (ALB), lactate dehydrogenase (LDH), the number of pneumonia lobes, pulmonary consolidation lobes, the proportion of invasive mechanical ventilation (IMV), oxygen therapy, bronchoalveolar lavage, use of intravenous immunoglobulin (IVIG), and use of glucocorticoids were found to have statistically significant (p < 0.05) in the univariate analysis. Ultimately, length of hospitalization (OR = 1.203, p < 0.001) and number of pneumonia lobes (OR = 3.311, p < 0.001) were selected as clinical parameters via stepwise selection in multivariate logistic regression analysis (Table 2). These two variables were used to develop the clinical models. The AUCs of the LR, RF and SVM machine learning models were 0.878, 0.938, and 0.887, respectively, in the training cohort and 0.833, 0.849, and 0.804, respectively, in the testing cohort (Fig. 3).

Radiomics features extraction, selection and model construction

A total of 2,264 radiomic features were extracted from the ROIs on the CT images of 115 patients in the training group. Among the 30 patients with redelineated ROIs, 2,179 features (96.25%) extracted from the two delineated ROIs had an ICC > 0.85 and were included in the subsequent analysis.

Through the MI, mRMR, and LASSO methods, 10 optimal features (Table S1) were ultimately selected (Fig. 4). These 10 features include two first-order features, seven texture features, and one shape feature. The texture features included three GLSZM, one GLDM, and three GLCM. The radiomic score (Radscore) was calculated for each patient, forming the basis for constructing radiomic models.

The AUCs of the LR, RF and SVM machine learning models were 0.891, 0.960, and 0.919, respectively, in the training group and 0.734, 0.765, and 0.739, respectively, in the testing group (Fig. 3).

Combined model construction and comparison

We constructed combined models integrating the radiomic and clinical models via the LR, RF and SVM machine learning models. The detailed evaluation indicators for these models in the training and testing groups are presented in Table 3. In the testing cohort, the LR machine learning model performed exceptionally well with AUC (95% CI) of 0.890 [0.792–0.988], sensitivity of 0.810, specificity of 0.897, accuracy of 0.860, PPV of 0.850, NPV of 0.867, and F1 score of 0.829, outperforming the other two models. The ROC curves for all three machine learning models in both the training group and testing group are shown in Fig. 3. The AUCs of the LR, RF and SVM machine learning models were 0.946, 0.977, and 0.971, respectively, in the training group and 0.890, 0.859, and 0.885, respectively, in the testing group. The LR machine learning model achieved an AUC of 0.946 in the training group and 0.890 in the testing group, outperforming those from both clinical and radiomics models.

The Delong test was utilized to compare the LR machine learning models’ performance (Table 4). The combined model demonstrated significantly better predictive performance than the Radscore model in both the training and testing cohorts (p = 0.031 in the training cohort and p = 0.023 in the testing cohort). In the training cohort, the combined model’s predictive ability was notably superior to that of the clinical model (p = 0.002), whereas no significant difference between the two models in the testing cohort (p = 0.154).

We further assessed the LR machine learning model’s performance via net reclassification improvement (NRI) and integrated discrimination improvement (IDI) (Table 4). The IDI showed that the combined model outperformed both the clinical model and the radiomic model in both the training group and the testing group. In the testing group, the combined model had 0.338 (p < 0.001) better prediction accuracy than the clinical model and 0.219 (p = 0.011) better prediction accuracy than the radiomic model.

The Fig. 5 illustrates the calibration curve and decision curve analysis (DCA) of the LR machine learning model. The calibration curve indicates that the model’s predicted probabilities align well with actual probabilities, demonstrating good calibration. The DCA demonstrates that the combined model outperformed both the clinical and radiomic models in terms of predictive performance when considering a risk threshold between 0.2 and 1.0.

Given that the LR machine learning model showed excellent prediction performance, we used two independent risk factors (the number of pneumonia lobes and the length of hospital stay) and the Radscore to construct a nomogram (Fig. 6) based on LR machine learning model to facilitate clinical decision-making.