This study was conducted utilizing the multicenter prospective CSTAR cohort. An inception cohort of 2404 newly diagnosed lupus patients without thrombocytopenia at registration was formed. Among these patients, 41 cases were lost to follow-up, and 93 cases had missing core data. Consequently, 2270 patients were included in the final analysis. In our cohort, there existed missing data in 213 participants, visually presented in Figure S1 (see Supplementary Material 2). The missing proportions of predictive variables were all lower than 10%. The MICE method was employed for multiple imputations, generating five imputed datasets. Variable selection and model fitting were conducted separately on the five datasets, and the results were combined according to Rubin’s rules. Notably, 2057 patients without missing data were defined as the complete dataset and used for sensitivity analysis. Detailed information is available in the study flowchart (Fig. 1).

In the study, 2270 patients were included, of whom 80 developed TP during follow-up, resulting in an incidence of 3.52%. The baseline characteristics of this cohort are presented in Table 1. The median follow-up time was 629 (203–1305) days, approximately 1.72 years. The mean age at lupus onset was 33.79 ± 12.48 years, with 91.9% female patients. The median SLEDAI score at registration was 6 (2–10), with 287 cases (12.6%) showing baseline organ damage. Clinical manifestations included mucocutaneous (69.1%), musculoskeletal (54.8%), renal (40.6%), neuropsychiatric (8.6%), leukopenia (36.9%), and hemolytic anemia (7.4%). Approximately 76.7% of patients were discovered with anti-dsDNA antibody positivity, followed by anti-SSA (58.9%), anti-RNP (45.3%), and anti-Sm (43.4%). In our cohort, 87.5% of patients received glucocorticoids, with 110 cases (4.8%) receiving GC pulse therapy. HCQ (82.1%) was prescribed as a cornerstone treatment, and the five most common immunosuppressants were MMF (24.2%), CYC (17.2%), MTX (12.8%), tacrolimus (5.9%), and CsA (3.5%). Additionally, monoclonal antibodies and IVIG were used in 2.7% and 1.8% of patients, respectively.

Patients were divided into the SLE-TP group (n = 80) and the SLE-non-TP group (n = 2190). As shown in Table 1, the univariate analysis indicated that patients with baseline organ damage (SDI ≥ 1, 23.7% vs. 12.2%, p < 0.001), hemolytic anemia (16.2% vs. 7.1%, p = 0.004), hypocomplementemia (55.8% vs. 39.3%, p = 0.002), and serositis (26.2% vs. 16.2%, p = 0.026) were more likely to develop TP during follow-up. Higher positivity rates of anti-Histone (34.6% vs. 23.5%, p = 0.034), anti-β2GPI (29.9% vs. 18.6%, p = 0.020), and Coombs test (38.4% vs. 19.1%, p < 0.001) were also observed in the SLE-TP group. Additionally, more patients in the SLE-TP group received GC pulse, CYC, tacrolimus, and IVIG therapy. The characteristics of patients with different TP severity were displayed in Table S1 (see Supplementary Material 2). Higher positivity rates of anti-Sm and ACL were discovered in patients with moderate to severe TP.

Incorporating literature review and expert opinions, 28 variables available in the CSTAR cohort were selected as preliminary candidate variables (see Table S2 in Supplementary Material 2). To enhance accuracy and reduce overfitting, the Lasso regularization method was employed for further variable selection (Figure S2, Supplementary Material 2). The β coefficients from the five imputed datasets are displayed in Table S3 (see Supplementary Material 2).

After confirming the clinical significance of the predictors, six variables were selected: baseline SDI score ≥ 1 (HR 2.207, 95% CI 1.350–3.609, p = 0.002), hemolytic anemia (HR 1.953, 95% CI 1.017–3.750, p = 0.044), low complement (HR 2.351, 95% CI 1.004–5.505, p = 0.049), anti-β2GPI antibody positivity (HR 1.805, 95% CI 1.084–3.004, p = 0.024), Coombs test positivity (HR 1.878, 95% CI 1.123–3.141, p = 0.017), and anti-histone antibody positivity (HR 1.595, 95% CI 1.017–2.587, p = 0.059). Definitions for each variable are detailed in Table S2 (see Supplementary Material 2). The final model was constructed using a multivariable Cox proportional hazards regression model. The β coefficients, estimated HRs, and 95% CIs of the risk predictors are presented in Table S4 (see Supplementary Material 2) and Fig. 2. For examination of the proportional hazard assumption, Schoenfeld individual tests of each variable are shown in Figure S3 (see Supplementary Material 2). Figure S4 (see Supplementary Material 2) plotted the cumulative risk curves for new-onset TP in lupus patients stratified by each predictor.

Forest plot of thrombocytopenia onset prediction model in SLE patients

The cumulative thrombocytopenia risks of lupus individuals can be calculated using the following formulas: P (1-year) = 1 − 0.9955971exp(onset index), P (2-years) = 1 − 0.9940235exp(onset index), P (3-years) = 1 − 0.9905769exp(onset index). The onset index is defined as 0.7917 × SDI ≥ 1 + 0.6694 × hemolytic anemia + 0.8549 × low complement + 0.5905 × anti-β2GPI + 0.6304 × Coombs + 0.4666 × anti-Histone. The predictive model algorithm is illustrated in a nomogram (Fig. 3).

Nomogram for thrombocytopenia onset prediction model in SLE patients

Nomogram aids in predicting the probability of developing thrombocytopenia for a lupus patient. Points for baseline SDI score ≥ 1, hemolytic anemia, low complement, Coombs test positivity, anti-β2GPI antibody positivity, and anti-histone antibody positivity can be obtained using a point caliper and summed for total points. The total points can be matched to the 1-year, 2-year, and 3-year cumulative incidence scales.

Model performance was evaluated among 2270 patients with 80 events. The apparent total C-index of the prediction model was 0.714 (0.653–0.775). The apparent C-index values for risk prediction at one, two, and three years were 0.741 (0.660–0.823), 0.730 (0.655–0.805), and 0.710 (0.643–0.777); while the Brier scores were 0.018 (0.012–0.024), 0.025 (0.017–0.032), and 0.037 (0.027–0.046), respectively. After enhanced Bootstrap resampling, the optimism-corrected C-index values for total, 1-year, 2-year, and 3-year predictions were 0.702, 0.726, 0.715 and 0.695. The optimism of the Brier score was lower than 0.001, resulting in minimal differences between apparent and optimism-corrected values. Detailed information is provided in Table S5 (see Supplementary Material 2). The calibration curve (1000 Bootstraps) depicted the comparison between predicted and observed risk probabilities (Fig. 4). Sensitivity analysis was conducted using the complete dataset (2057 patients, 73 events), as shown in Table S6 (see Supplementary Material 2).

The receiver operating characteristic (ROC) curves of the model and calibration curves in internal validation. A. Time-dependent ROC curves of the predictive model. B. 1-year calibration curve. C. 2-year calibration curve. D. 3-year calibration curve

The nomogram scores for SDI ≥ 1 (baseline organ damage), hemolytic anemia, low complement, anti-β2GPI positivity, anti-histone positivity, and Coombs test positivity were 91, 80, 100, 60, 59, and 71 points, respectively. The total point was calculated for every patient. A total risk score of 211 was determined as the optimal cut-off point. According to the nomogram-derived risk scores, participants were classified into a high-risk group (> 211) and a low-risk group (0–211), with TP-onset risks of 7.59% and 2.52%, respectively. The risk stratifications for new-onset TP are shown in Table S7 (see Supplementary Material 2). The cumulative risk curves of two different risk groups were drawn in Fig. 5.

Cumulative risk curve of thrombocytopenia onset in SLE patients after risk stratification. A. Imputed data. B. Complete data