Clinical parameters associated with childhood ALL patients up to June 10, 2019 were downloaded from the NCI TARGET database¹. A total of 1,842 patients were included in the first and second phases of the data set. Of these subjects, patients with cALL who had information of main observation indicators were target subjects of this study, which include survival time, survival status, age at diagnosis, gender, ETV6 RUNX1 fusion status, TRISOMY 4 10 status, MLL status, TCF3 PBX1 fusion status, karyotype, BCR-ABL1 fusion status, central nervous system (CNS) status at diagnosis, BMA blasts day 8, BMA blasts day 29, cell of origin, Down syndrome status, MRD at day 29, and DNA index. The specific screening process is presented in Figure 1. The clinicopathological characteristics of patients in the cohorts are listed in Table 1. All source data are presented in Supplementary Table S1.

Age refers to the age at diagnosis. The WBC count at diagnosis was the absolute peripheral WBC count (in ×10³/mcL). CNS status at diagnosis was determined according to the status of CNS leukemia (CNSL) at the time of diagnosis. Diagnosis and typing were according to Pediatric Acute Lymphoblastic Leukemia, Version 2.2020, NCCN Clinical Practice Guidelines in Oncology (23). MRD status was determined by flow cytometry in bone marrow on day 29. BMA blasts day 29 and day 8 represent the percent of blasts in bone marrow aspirate at day 29 or day 8 of induction therapy. The presence or absence of ETV6 RUNX1 fusion status, MLL status, TCF3 PBX1 status, and BCR-ABL1 fusion status were detected by fluorescence in situ hybridization (FISH), PCR, or cytogenetics. The type of cell origin is determined by bone marrow examination. Karyotype analysis was determined by chromosomal banding technique to find abnormal chromosome number of leukemia cells and structural changes such as translocation, inversion, and deletion (23).

Age was divided into two groups based on the cutoff value of 10 years old (24, 25), and WBC count at diagnosis was divided into three groups: “<50,” “50 to 100,” and “≥100” (3, 26, 27) according to the WBC count (in ×10³/mcL). Next, CNS status at diagnosis was divided into three groups: CNS1, CNS2, and CNS3. MRD day 29, which means minimum residual disease status at day 29 of induction therapy, was divided into four groups: <0.01, 0.01–0.1%, 0.1–1%, and >1% (28, 29). BMA blasts day 29, which represents the percent of blasts in bone marrow aspirate at day 29 of induction therapy, was divided into two groups with a 5% cutoff. Besides, BMA blasts day 8 was divided into two groups according to whether it is >20%. The karyotypes are divided into four groups according to the trisomy of chromosomes 4, 10, 17, and 18, which include “no trisomies in 4, 10, 17, 18,” “4, 10, 17, 18 have only one trisomy,” “double trisomies (DT),” and “triple trisomies (TT).” Cell of origin contained two subtypes, which included B cell ALL and T cell ALL. Finally, the DNA index, whose number represents the ratio of the DNA content or chromosome number in a tumor sample compared to a normal diploid sample, was divided into two groups: ≤0.8 and >0.8. ETV6 RUNX1 fusion status, MLL status, TCF3 PBX1 status, and BCR-ABL1 fusion status were divided into positive and negative according to whether its corresponding fusion gene were positive or negative.

All data including demographic and disease characteristics were expressed as count (%). Statistical analysis was performed using the R software (Version 3.6.1)².

The prognostic value of the 16 clinical characteristics was first calculated in the univariate Cox analysis; clinical features with a P < 0.1 in the multivariate Cox regression analysis were used to construct the nomogram via the “rms,” “survival,” and “foreign” packages of R (R version 3.6.1). Hazard ratio (HR) and 95% confidence interval (CI) were calculated. The performance of the nomogram was measured by the C-index and assessed by comparing the nomogram-predicted estimates versus the observed Kaplan–Meier estimates of survival probability (R package “rms”) (30). Based on the regression coefficients of the multivariate Cox regression analysis, a risk score composed of six clinical features in the nomogram was calculated, and the patients were divided into two groups by taking the corresponding median risk score as the cutoff point.

Kaplan–Meier curves and the log-rank test were used to compare the survival outcomes of the two groups with the R packages “survminer” and “survival” (R version 3.6.1). Receiver operating characteristic (ROC) curve analysis was employed to compare prediction concerning the accuracy and precision with the R package “survivalROC” (R version 3.6.1). Bootstrapping validation (1,000 bootstrap resamples) was used to calculate a relatively corrected C-index of the nomogram (31, 32). A P < 0.05 was considered significant.

The cohort included 673 cALL patients with complete clinical information of main observation indicators. The median of age at diagnosis was 5.5 years (range, 1.0–15.0 years). Three-hundred fifty-nine (53.3%) patients in the cohort were male, while 314 patients were female. The median of survival time was 3,016 days (range, 28–5,598 days). Survival status showed 120 deaths and 553 survivals.

In the primary cohort, we performed a univariate Cox regression analysis for each clinical factor (Table 2) and screened factors with P < 0.1, which included age, WBC count at diagnosis, CNS status at diagnosis, ETV6 RUNX1 fusion status, TRISOMY 4 10 status, karyotype, TCF3 PBX1 fusion status, BCR-ABL1 fusion status, MRD day 29, BMA blasts day 8, Down syndrome, and DNA index. Then, these factors were analyzed in a multivariate Cox regression model (Table 3). The results show that the patient and disease variables significantly associated with survival in multivariable modeling included positive ETV6 RUNX1 fusion (HR = 0.293, 95% CI = 0.124–0.695, P = 0.005), triple trisomies (TT) (HR = 0.211, 95% CI = 0.049–0.907, P = 0.036), MRD day 29 with 0.01–0.1% (HR = 1.942, 95% CI = 1.113–3.330, P = 0.016), MRD day 29 with 0.1–1% (HR = 2.140, 95% CI = 1.262–3.628, P = 0.005), MRD day 29 > 1% (HR = 2.508, 95% CI = 1.345–4.675, P = 0.004), and DNA index ≤0.8 (HR = 3.617, 95% CI = 1.233–10.612, P = 0.019).

The prognostic nomogram that integrated all significant independent factors for overall survival (OS) in the primary cohort is shown in Figure 2. The nomogram demonstrated that the karyotype had the largest contribution to prognosis, followed by ETV6 RUNX1 fusion status, DNA index, and CNS status. BMA blast day 8 and MRD day 29 showed a moderate effect on survival rate. Each category within these variables was assigned a point on the top scale based on the coefficients from multivariate Cox regression. By summing all of the assigned points for the eight variables and drawing a vertical line between total points and survival probability axis, we were easily able to obtain the estimated probability of 3-, 5-, and 7-year survival. The C-index for OS prediction were 0.754 (95% CI, 0.715–0.793) and 0.731 for the primary cohorts and bootstrapping validation, respectively. The calibration plot for the probability of survival at 3, 5, and 7 years showed an optimal agreement between the predictions by the nomogram and the actual observations (Figure 3).

We calculated the risk score, which was composed of six clinical features in the nomogram, based on the regression coefficients from the multivariate Cox regression analysis. We divided all patients in the training set into two groups, the high-risk score group (n = 336) and the low-risk score group (n = 337), by taking the corresponding median risk score as the cutoff point. Next, we compared survival predictions with regard to specificity and sensitivity according to the risk scores and clinical characteristics in the nomogram by ROC curve analysis. The results show that the areas under the curve (AUCs) of the risk score for 3-, 5-, and 7-year survival were 0.775, 0.776, and 0.772, respectively, which were higher than those of all of the clinical factors in the nomogram (Figures 4A–C). Patients with a high-risk score had a markedly worse OS than patients with a low-risk score (log-rank test P < 0.0001, Figure 4D).

In general, subtype patients with positive ETV6 RUNX1 fusion gene showed better OS than patients with negative ETV6 RUNX1 fusion gene, and patients in the high-risk group had worse OS than those in the low-risk group in ETV6 RUNX1 fusion-negative segment, while in the ETV6 RUNX1 fusion gene-positive segment, all patients were in the low-risk group (Figure 5A). Besides, the Kaplan–Meier curves of OS suggested that hypodiploid patients had significantly worse OS than patients not. What is more, patients in the high-risk group had worse OS than those in the low-risk group in non-hypodiploid segment, and all patients of the hypodiploid segment were in the high-risk group (Figure 5B). These above results indicated that the risk scores we have obtained are quite accurate and predictive.

To validate the robustness of our nomogram, we reviewed the data from the TARGET database again retrospectively, and the filtering process is presented in Supplementary Figure S1. We ended up with 299 valid cases of data as independent validation queues. Through a similar analysis process, we get supportive results. As Supplementary Figure S2 shows, in the validation cohort, the AUCs of the risk score for 3-, 5-, and 7-year survival were 0.683, 0.723, and 0.737, respectively, which were higher than those of all of the clinical factors in the nomogram (Supplementary Figure S2D). Patients with a high-risk score had a markedly worse OS than patients with a low-risk score (log-rank test P < 0.0001, Supplementary Figure S2A). The C-index of the nomogram for predicting survival was 0.703 (95% CI, 0.640–0.766). Besides, to further confirm the robustness of the model, we also validate the nomogram in the total cohort (N = 299 + 673 = 972). Patients with a high-risk score had a markedly worse OS than patients with a low-risk score (log-rank test P < 0.0001, Supplementary Figure S2C), and the AUCs of the risk score for 3-, 5-, and 7-year survival were 0.753, 0.783, and 0.773 (Supplementary Figure S2B), respectively. The C-index of the nomogram for predicting survival was 0.723 (95% CI, 0.684–0.762).