To prepare for the next COVID-19 pandemic, we established a clinical severity stratification model using machine learning with immune signatures. Three COVID-19 cohorts (one discovery cohort and two validation cohorts) and 16 age- and sex-matched healthy volunteers (negative for SARS-CoV-2 and virus-specific Immunoglobulin M [IgM] and Immunoglobulin G [IgG], as indicated by the reverse transcription-polymerase chain reaction [RT-PCR] test) were included in this study. According to the clinical severity classification criteria ( Supplementary Table S1 ), which was modified from World Health Organization guidelines (2), patients in the discovery cohort were classified into mild and critical cases. We screened potential variables by longitudinally comparing the levels of anti-SARS-CoV-2 antibodies, inflammatory cytokines, plasma complement components, and cellular immune signatures between critical and mild cases. A self-designed 42-parameter panel, including nine energy metabolism enzymes, was applied to phenotypic immune signatures using CyTOF technology. The most clinically relevant immune signatures and plasma parameters were introduced into machine learning.

Indices of interest, including the levels of inflammatory cytokines, complement components in plasma, and anti-SARS-CoV-2 antibodies, were extracted from electronic medical records ( Table 2 ). Specifically, they were the systemic LDH, lactate, complement component 3 (C3), complement component 4 (C4), 50% hemolytic unit of complement (CH50), IgG, immunoglobulin A (IgA), IgM, immunoglobulin E (IgE), interleukin-1 (IL-1), interleukin-2 (IL-2), interleukin-4 (IL-4), interleukin-5 (IL-5), IL-6, interleukin-8 (IL-8), interleukin-9 (IL-9), interleukin-10 (IL-10), interleukin-12 (IL-12), interleukin-13 (IL-13), interleukin-17 (IL-17), interferon-α (IFN-α), interferon-γ (IFN-γ), tumor necrosis factor-α (TNF-α), granulocyte colony-stimulating factor, granulocyte macrophage colony-stimulating factor, vascular endothelial growth factor, macrophage inflammatory protein-1-α (MIP1-α), and monocyte chemotactic protein-1. All data were collected and verified by two experienced doctors.

PBMCs were isolated from peripheral blood using Ficoll density gradient centrifugation. To sort cell precipitates, they were combined with 5 mL of fluorescence-activated cell sorting (FACS) buffer (1×phosphate buffered saline supplemented with 0.5% bovine serum albumin) and centrifuged at 400×g for 5 min at 4°C. The supernatant was discarded and the cell precipitates were resuspended in FACS buffer. To examine the samples, the viability rate must be greater than 85%, and the number of cells must not be less than 3×10⁶.

To ensure homogeneous staining, approximately 2×10⁶–3×10⁶ PBMCs were used for each patient. PBMCs were stained with cisplatin (Fluidigm) (0.1 uL, 2 min, room temperature) for live/dead, washed with cell staining buffer (CSB) (Fluidigm), and spun down (300×g, 5 min, room temperature). PBMCs were then incubated with human TruStain FcX (BioLegend) for 10 min at room temperature. After incubation, PBMCs were stained with 50 uL surface receptor staining mix (30 min, room temperature) and washed twice with CSB (300×g, 5 min, room temperature). Next, the PBMCs were incubated with FixL buffer (Fluidigm) for 15 min at room temperature and washed twice with Perm-S buffer (Fluidigm) (800×g, 5 min, room temperature). PBMCs were stained with 50 uL intracellular mix (30 min, room temperature) and washed twice with CSB (800×g, 5 min, room temperature). PBMCs were fixed in 1 mL 1.6% paraformaldehyde. Samples were fixed and permeabilized by incubating 1 mL Fix and Perm buffer (Fluidigm) with 1 uL nucleic acid Ir-Intercalator (Fluidigm) overnight at 4°C. Metal-conjugated antibodies and other reagents are listed in Supplementary Table S2 .

Before acquisition, PBMCs were washed twice with CSB and resuspended at a concentration of 1.1×10⁶ cells/mL in the Cell Acquisition Solution (Fluidigm) containing 10% EQ Four Element Calibration Beads (Fluidigm). PBMCs were acquired using a Helios CyTOF Mass Cytometer (Fluidigm) equipped with a SuperSampler fluidics system (Victorian Airships), and data were collected as previously described. fcs files.

After acquisition, data were concatenated using the fcs concatenation tool from Cytobank and manually gated to retain live, singlet, and valid immune cells. CytoNorm was used in two steps according to the instructions provided in the R library CytoNorm to normalize the data (16). For the downstream analysis, the fcs files were loaded into R. The signal intensities for each channel were arcsinh-transformed with a cofactor of 5 (x_transf = asinh(x/5)). To visualize high-dimensional data, t-distributed stochastic neighbor embedding analysis (t-SNE) (17) and flow self-organizing map (FlowSOM) (18) algorithms were performed on all samples. Approximately 10,000 cell events in each sample were pooled and included in the t-SNE analysis, with a perplexity of 30 and theta of 0.5. The R t-SNE package for the Barnes Hut implementation of the t-SNE was used in this study. To study the developmental trajectory of natural killer (NK) cells and classical monocytes, dynamic immunometabolic states and cell transitions were analyzed using the Monocle algorithm (19). Data are displayed using the ggplot2 R package.

Since the target variable (clinical severity) for model training was labelled data, provided by clinical experts. The supervised learning methods are more appropriate than unsupervised-, semi-supervised-, and reinforcement learning methods. By comparing the advantages of different supervised methods (20–30)( Supplementary Table S3 ), we finally employed AdaBoost, Back Propagation, Gradient Boosting Decision Tree, Random Forest, and Support Vector Machine algorithms to construct classifiers for discriminating patients with critical COVID-19 from mild ones. The important immune and systemic features (cDCs and LDH) were introduced to the model as inputs. Five-fold cross-validation (with four folds for training and one-fold for validation) and external validation were performed. For five-fold cross-validation, all the training data were randomly split into five parts. Each part was considered as the training part and the others were used for validation. Here, we performed the five-fold cross-validation five times and the averaged values of AUC were adopted. For the external validation, Back Propagation algorithm was performed.

Statistical analyses were performed using the R software (v.4.0.4). The normality of patient data was tested using the Shapiro–Wilk normality test. Statistically significant differences between phenotypes were calculated using two-sided multiple Student’s t-tests for variables with a normal distribution and Wilcoxon rank-sum tests for other variables. Spearman’s correlation analysis was performed on significantly different clusters, cytokines, and clinical indicators to assess their correlations using the R package stats (4.1.0). Receiver operating characteristic (ROC) analysis was performed with the R package pROC (1.16.2), and a heatmap was generated with the R package ggplot2 (4.0.5).

A total of 39 individuals diagnosed with COVID-19 (15 mild and 24 critical cases) admitted to our ICU were included in cohort 1 as the discovery cohort to determine potential predictive parameters. As shown in Table 1 , the basic information of the critical and mild cases was comparable. The Murry lung injury score, length of mechanical ventilation, and length of ICU stay were significantly high in critical cases ( Table 1 ). Longitudinal comparisons of inflammatory cytokines, antibodies, and complement components revealed that systemic IL-6, IL-12, and LDH levels were important in distinguishing mild cases from critical cases. The variation trends in these parameters were consistent across all sampling points ( Table 2 ).

To acquire a full landscape of the immune signatures of PBMCs and identify the potentially important clusters for the stratification of COVID-19, we performed CyTOF analysis with a 42-parameter panel (consisting of 33 surface markers and 9 intracellular metabolic markers) ( Supplementary information, Figure S1 ). The obtained data were subjected to a FlowSOM clustering algorithm and t-distributed stochastic neighbour embedding (t-SNE) analysis, which enabled the identification of distinct clusters representing different immune cell types. According to the dimensional reduction results of the marker expression level, 34 clusters were obtained ( Figure 1A ). Then, to provide reference for other similar studies, which may apply different panels, we further classified these 34 clusters into “eleven major immune cell populations” (CD4⁺ T, CD8⁺ T, γδT, DPT, DNT, pDC, cDC, NK, NKT, B, and Monocytes), which were often studied ( Supplementary information, Table S4; Figure S2 ).

We found that the composition of PBMCs in patients with COVID-19 varied significantly from that in healthy volunteers. The total counts of PBMCs (in per millilitre of peripheral blood) and the counts of the main immune cell types, such as T, B, and NK cells, of patients with COVID-19 decreased significantly. However, the number of monocytes increased ( Supplementary information, Figure S2 ). Comparison of the percentages of all defined 34 clusters further confirmed that, fifteen immune cell subsets were significantly differed between COVID-19 patients and healthy volunteers ( Figures 1B, C ). Most of these subsets were acquired immune cell subsets and were significantly decreased in COVID-19. In addition, variations in two major innate immune cell subsets (NK cells (C03) and classical monocytes (C12), with the average percentages more than 5% in healthy volunteers) were also found ( Figures 1B, C ). As the host innate immunity is the first line of defense, we further investigated these two subsets’ metabolic status. As shown in Figure 2 , the metabolic markers participating in the process of glucose (such as CS, GLS, PFKFB3, and PDk1_pS241) and tryptophan metabolism (IDO1 and KAT1) were significantly altered in both NK cells (C03) and classical monocytes (C12). The developmental trajectories further demonstrated that under COVID-19, NK cells gradually transformed from C01 to C03, namely, from a relative metabolic steady state to a disturbed state with decreased oxidative phosphorylation but boosted glycolysis and tryptophan catabolism ( Figures 2A–C ). For classical monocytes, C12 gradually transformed to C09, namely, to tryptophan exhaustion ( Figures 2D–F ).

As described in the Methods section, to identify the important clusters distinguishing critical cases from mild cases, we compared the cell counts and percentages of each cluster within PBMCs at all sampling points. In total, five candidate clusters were found, and the differences in cDCs (C07), Tregs (C20), CD4⁺ Tcm (C24), pDCs (C05), and DPT (C29) were shared by the results from all sampling points and the first day samples ( Figure 3A ). As the percentages of pDCs and DPT were below 0.5%, they were not considered in subsequent analyses. Next, we investigated whether these clusters were associated with clinical parameters and prognosis. The results demonstrated that the counts of cDCs, Tregs, and CD4⁺ Tcm were significantly decreased in the critical cases and patients who ultimately died ( Figures 3B, C ). Their levels were positively or negatively correlated with systemic parameters, lung injuries, and the length of mechanical ventilation ( Figures 3D–F , and Supplementary information, Table S5 ). Within each severity group, the longitudinal analysis showed that the counts of these three clusters were not significantly different among different sampling points ( Supplementary information, Figure S3 ). These findings indicated that altered cDCs, Tregs, and CD4⁺ Tcm were stable/sensitive predictive biomarkers because their level wouldn’t be significantly influenced by sampling timing and/or transient condition relief. Specifically, cDCs was the most important cluster, negatively correlated with LDH and positively correlated with IL-2, IL-12, TNF-α, and MIP1-α ( Figure 3E ). Receiver operating characteristic analysis further revealed that the single variable cDCs was effective in predicting critical COVID-19 ( Figure 3G ). And the level of LDH was the most important systemic parameter because of its strong negative correlation with cDCs, Tregs, and CD4⁺ Tcm ( Figure 3E ).

Considering the potential of machine learning for disease severity stratification, we developed clinical severity stratification models based on important key clusters (cDCs, Tregs, and CD4⁺ Tcm) and systemic parameters (LDH, IL-6, IL-12). As we expected, machine learning models with six parameters as inputs showed good effects in predicting clinical severity ( Figure 4A ). Among these parameters, cDCs and LDH were the most important immune signature and systemic signature, respectively ( Figure 4B ). The model using cDCs and LDH as individual input also performed well, with an average AUC of approximately 0.8 in the discovery cohort ( Figures 4C, D ). The validation of machine learning models with single input (with Back Propagation algorithm) further demonstrated that the clinical severity stratification model based on single cDCs had an AUC of 0.77 ( Figure 4E ). And the model based on systemic LDH had an AUC of 0.88 ( Figure 4F ). Notably, patients in validation cohort 1 were recruited in 2020 and infected with a different strain compared with the patients in the discovery cohort. These results indicate that our models, based on single biomarker (cDCs or LDH), performed well in COVID-19 severity stratification, with good robustness and generalization.

To provide detail reference for clinicians in quick decision-making for the next pandemic, we analyzed the effect of cDCs and LDH in severity prediction in validation cohorts and tried to find out the optimal reference limits. In validation cohort 1 (adopted from Stephenson et al.’s published work (15)), the proportion of cDCs decreased in critically ill participants across the three UK centers ( Supplementary information, Figures S4A–C ). The percentage of cDCs showed good effects in predicting clinical severity (AUC = 0.74, Figure 4G ). The optimal cutoff point was 1.2%, and the sensitivity was 0.93 (95% CI 0.70-0.99). In the validation cohort 2 (adopted from Peking University Third Hospital), similar with the findings in the discovery cohort, significant increase of LDH ( Supplementary information, Figure S4D ) and its predictive effect was found (AUC = 0.89, Figure 4H ). The cutoff point was 270.5 U/L and the sensitivity was 0.92 (95% CI 0.86-0.95). Accordingly, the reference limits of cDCs and LDH for critical COVID-19 were less than 1.2% and more than 270.5 U/L.