In this study, a total of 506 pediatric patients diagnosed with KD between February 2013 and November 2023 were enrolled at the Children’s Hospital of Fudan University. Venous blood samples were collected from KD patients at the time of initial evaluation in hospital. The blood analysis of the samples was conducted in the laboratory department of our hospital. The demographic and laboratory data were extracted from the medical record. All patients met the diagnostic criteria outlined in the Expert Consensus on the Diagnosis and Acute-Phase Treatment of Kawasaki Disease (Subspecialty Group of Rheumatology, 2022), which include persistent fever for ≥5 days and at least four of the five principal clinical features: polymorphous rash, bilateral nonexudative conjunctival injection, changes in lips and oral cavity, changes in the extremities, and cervical lymphadenopathy. The diagnosis of CALs was established based on echocardiographic findings, defined as a Z score ≥2 mm (Kuo, 2023).

Patients were excluded if they had received immunosuppressive therapy within the previous 3 months or had evidence of cardiac, hepatic, or renal insufficiency; active infections; or immunodeficiency disorders. All diagnoses and treatment decisions were made by one of two experienced pediatric clinicians specializing in KD. Detailed demographic and clinical characteristics of the enrolled patients are presented in Table 1.

For a fair comparison of performance across different input feature sets, rigorous data preprocessing procedures were implemented using the scikit-learn library in Python (version 3.9.13). Missing values in continuous variables were imputed using median substitution, a robust univariate method that reduces sensitivity to outliers while preserving the central tendency of the data. Continuous features were standardized using z-score normalization (mean = 0, standard deviation = 1) to ensure comparability across variables and to enhance algorithmic convergence. Where applicable, categorical variables were transformed using one-hot encoding to enable compatibility with machine learning models.

Following preprocessing, the dataset was randomly partitioned into training and testing subsets in an 80:20 ratio. Stratified sampling was employed to maintain consistent class distributions across subsets, thereby minimizing potential sampling bias due to class imbalance. All preprocessing steps were conducted prior to model training and cross-validation to prevent data leakage and ensure methodological rigor.

We employed five distinct algorithms to identify the most informative predictors from a large pool of candidate variables: Corrected as Gradient Boosting Decision Tree (GBDT), Extra Tree (ET), Random Forest (RF), Logistic Regression (LR), and Least Absolute Shrinkage and Selection Operator (LASSO) regression. These algorithms were selected for their capacity to rank feature importance based on different theoretical foundations—tree-based ensemble methods, linear coefficients, and regularization penalties. The aim was to adopt a data-driven approach that retains features with high predictive value while eliminating redundant or irrelevant variables, thereby enhancing model stability and performance on unseen data.

To develop predictive models for the diagnosis of CALs in KD patients, we implemented eight distinct ML algorithms: Support Vector Classifier (SVC), LR, RF, ET, Decision Tree Classifier (DT), Gaussian Naïve Bayes (Gaussian NB), K-Nearest Neighbors (KNN), and Multi-Layer Perceptron Classifier (MLPC). These algorithms were selected to represent a diverse range of classification paradigms, encompassing linear models, ensemble methods, probabilistic models, distance-based learning, and neural networks. All selected models are widely used in biomedical research and offer complementary strengths. Each algorithm was chosen for its balance between interpretability and capacity to capture complex linear or nonlinear relationships among input features.

To assess the models’ discriminative power, receiver operating characteristic (ROC) curves were plotted, and the area under the curve (AUC) was calculated. Mul-tiple evaluation metrics were computed to comprehensively assess classification performance, including precision, recall, accuracy, and F1-score, all derived from the confusion matrix.

To rigorously evaluate model performance and minimize the risk of overfitting, stratified 10-fold cross-validation was conducted on the training dataset ensuring class balance across folds. Average performance across all folds was reported to ensure robustness and generalizability of the models. Model calibration was assessed by generating calibration (reliability) curves, and the Brier score (BS) was computed as a quantitative measure of the accuracy of probabilistic predictions. Lower BS indicate better calibrated models. In addition, the Kolmogorov–Smirnov (KS) test was applied to evaluate the separation between predicted probability distributions of the positive and negative classes.

Finally, to evaluate the temporal robustness and real-world applicability of the developed models, external validation was performed using a temporally independent test cohort collected after the model development period. This prospective validation strategy provided further evidence of the model’s generalizability to future clinical data.

Summary data on outcomes were collected from published GWAS meta-analyses and publicly available data. These summary data were analyzed by MR to determine if there was a causal association between selected features and the risk of coronary artery disease. In order to increase the reliability of the study results, the causal relationship between selected features and coronary artery disease risk was investigated using five Mendelian randomization methods. They are MR Egger, weighted median, inverse variance weighted (IVW), simple mode and weighted mode. IVW, which assumes that each genetic variant exists independently and can influence outcome only through the exposure of interest and combines the Wald ratios of individual SNPs, was employed as the principal method of analysis in this study. However, causality may be biased in the presence of pleiotropy (Bowden et al., 2015; Grover et al., 2017). The remaining four methods were used as complementary methods to IVW, although they are less powerful (Chen et al., 2020). A statistically significant association between exposure and outcome was deemed to be present when the p-value was found to be less than 0.05.

Heterogeneity tests were carried out for statistically significant results using Cochran’s Q-test (p < 0.05 was considered heterogeneity). Meanwhile, we used MR-Egger intercept tests to detect pleiotropy (p < 0.05 was considered pleiotropy) (Burgess and Thompson, 2017; Chen et al., 2023). Finally, the leave-one-out sensitivity analysis was performed to examine if one single SNP drove the causal association. In this study, R software and the “Two Sample MR” package were used for all MR analyses.

SPSS 25.0 was used for data analysis. For measurement data, the D'Agostino-Pearson omnibus test was first used to assess normality. The measurement data conforming to normal distribution were expressed as mean ± standard deviation ( X ¯ ± S), and non-normally distributed measurement data were expressed as median (Interquartile range) (M (Q25, Q75), %). The differences between two groups were compared using an independent samples t-test or the non-parametric Mann-Whitney U test. P < 0.05 was considered statistically significant.

A comparative analysis delineated multiple laboratory biomarkers between KD patients with and without CALs. The two groups showed significant differences in a variety of indicators (Figure 1A). KD patients with CALs predominantly fell into the high-risk of abnormalities in coagulation system, including activated partial thromboplastin time (APTT), thrombin time (TT) and prothrombin time (PT). We further performed a correlation analysis of laboratory biomarkers in KD patients to explore potential physiological or pathological associations. The results, shown in Figure 1B, illustrate the correlations among several key laboratory biomarkers in KD patients.

To reduce the risk of overfitting and improve the generalizability of the ML models, a comprehensive feature selection strategy was implemented. The results of feature importance ranking and selection across different classifiers are presented in Figure 2 and Supplementary Tables S1, S2. By integrating the results of ROC analysis with differential expression analysis, we identified a core set of seven key laboratory biomarkers as the most informative for model construction (Table 2). Cholinesterase (CHE) was excluded because it did not reach statistical significance in ROC curve analysis. These selected features included creatine kinase-MB (CK-MB), fibrinogen (FIB), international normalized ratio (INR), APTT, PT, red cell distribution width-standard deviation (RDW-SD), and C-reactive protein (CRP).

Based on selected features, different methods were used in order to get the best CALs risk prediction model. Figure 3 and Table 3 shows the results of the 8 ML models testing. Concerning the whole dataset, the RF model achieved the best performance with AUC of 0.90, whereas most other models gave an AUC above 0.8 (Figure 3A). However, LR and MLPC stand out with the high prediction accuracy 0.84. It means these models have excellent discriminating power in predicting CALs.

These parameters are defined as follows: Precision = True Positive/(True Positive + False Positive), Recall = True Positive/(True Positive + False Negative), Accuracy = True Positive + True Negative /(True Positive + True Negative + False Positive + False Negative), F1-score = 2×Precision × Recall/(Precision + Recall).

We further evaluated the models using the BS to assess the accuracy of probabilistic predictions. Among all models, RF, MLPC, and LR achieved the lowest BS, with values of 0.097, 0.110, and 0.113, respectively (Figure 4A), indicating superior probability calibration. Additionally, we applied K-fold cross-validation to compare model performance based on average test error. The top-performing models achieved a training accuracy converging around 0.8. RF demonstrated the highest overall accuracy on the test set (84.5%), followed by LR (83.9%) and ET (83.9%) (Figure 4B). The KS test results are shown in Figures 4C–J. A larger KS statistic indicates stronger discrimination between positive and negative classes. Interestingly, while RF performed well in overall accuracy, it showed relatively poor discrimination capability based on the KS statistic. In contrast, DT and MLPC models yielded KS curves closest to the true distribution, suggesting better class separation performance.

Based on the results of the above model evaluation, we selected three models with better performance for data validation. They were MLPC, ET and RF. We proposed a time-cross validation in total of 74 KD patients collected from future (2022–2023), including 36 patients with CALs and 38 patients without CALs. The external validation results were shown in Figure 5 and Table 4. The accuracy of MLPC is 0.78, ET is 0.76, and RF is 0.72. Therefore, we chose MLPC as our final prediction model. To facilitate the use of our prediction models, we developed this model (http://127.0.0.1:5000).

To further confirm the reliability of the model, Mendelian Randomization analysis was performed to confirm the relationship between selected features and CALs. Among all selected features, three features have a clear causal relationship with coronary heart disease. Figure 6A displays the other four additional MR analysis techniques results (MR-Egger, weighted median, simple model, and weighted model). IVW method results revealed evidence of a significant connection between APTT and CALs risk (OR = 0.809, 95% CI = 0.70–0.94, p = 0.004), as well as RDW (OR = 0.935, 95% CI = 0.88–1.00, p = 0.04). This association was further supported by the scatter plot (Figure 6B). Finally, we performed the leave-one-out analysis by removing each instrumental SNP to ensure that no single SNP heavily influenced the causal estimate. Forest plots were generated for each variable (Supplementary Figure S1). No SNP had a substantial effect size on the study’s estimation, indicating the robustness of the findings. The above findings demonstrate a consistent, genetically-based causal relationship between APTT or RDW and CALs.