A unified metric of human immune health

Abstract

Immunological health has been challenging to characterize but could be defined as the absence of immune pathology. While shared features of some immune diseases and the concept of immunologic resilience based on age-independent adaptation to antigenic stimulation have been developed, general metrics of immune health and its utility for assessing clinically healthy individuals remain ill defined. Here we integrated transcriptomics, serum protein, peripheral immune cell frequency and clinical data from 228 patients with 22 monogenic conditions impacting key immunological pathways together with 42 age- and sex-matched healthy controls. Despite the high penetrance of monogenic lesions, differences between individuals in diverse immune parameters tended to dominate over those attributable to disease conditions or medication use. Unsupervised or supervised machine learning independently identified a score that distinguished healthy participants from patients with monogenic diseases, thus suggesting a quantitative immune health metric (IHM). In ten independent datasets, the IHM discriminated healthy from polygenic autoimmune and inflammatory disease states, marked aging in clinically healthy individuals, tracked disease activities and treatment responses in both immunological and nonimmunological diseases, and predicted age-dependent antibody responses to immunizations with different vaccines. This discriminatory power goes beyond that of the classical inflammatory biomarkers C-reactive protein and interleukin-6. Thus, deviations from health in diverse conditions, including aging, have shared systemic immune consequences, and we provide a web platform for calculating the IHM for other datasets, which could empower precision medicine.

Extended Data Fig. 1 |. Subject demographics and further characterization of the serum protein and transcriptomic modules.

a, Density plot of patient and healthy participants’ age distributions (Kolmogorov-Smirnov test assessing difference between the two distributions, p = 0.41). Extended Data Fig. 1a–c only show data for individuals in primary set of participants; data for set-aside participants not shown but included in Extended Data Table 1. b, Box plots of participant ages by disease condition. Box plot center lines correspond to the median value; lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles), and lower and upper whiskers extend from the box to the smallest or largest value correspondingly, but no further than 1.5X inter-quantile range. c, Proportion of females (F) and males (M) in each disease/condition group. d, Pearson correlation between the protein (PM) or transcriptomic (TM) modules (columns) and cellular CBC/TBNK parameters (rows), computed with 198 (197) subjects with both whole blood transcriptome (serum protein) and CBC/TBNK data. *Adjusted P < 0.05. Abbreviations are same as in Fig. 1g. e, Conceptual illustration of parameter temporal stability, defined by low intra-subject variation relative to inter-subject variation. f, Variance assigned to the subject term in the variance partition analysis fit using only a subject random intercept (see Methods), run across each cellular parameter, protein module (PM), and transcriptomic module (TM). g, Percent variation explained by the subject term in the variance partition model in the protein and transcriptomic features using the variance partition model with only a subject random intercept (see Methods) as in (f). Proteins (left) and genes (right) are ordered on the x-axis by the percent variation explained by the subject term. h, Percent variation explained by the patient and medication covariate (showing effect of each medication individually) for each protein (left) and gene (right) measured. Medications were included in the model if they were used by many patients and not highly confounded with one of the condition groups. PLT, platelets; RBC, red blood cell parameters; WB: whole blood.

Extended Data Fig. 2 |. Disease heterogeneity in circulating cell frequencies.

a, Heatmap of complete blood count (CBC) and lymphocyte (T, B, NK cell) phenotyping (TBNK) parameters (rows) across patients and healthy participants (columns); columns and rows are ordered by hierarchical clustering. Top annotation row shows the participant age, middle row shows the large condition groups (n > 10 participants), and third row shows all condition groups regardless of number of participants. b, Patients and healthy participants shown in PC1 and PC2 space of CBC and TBNK parameters. Each parameter was standardized to unit variance and mean of zero prior to computation of the principal components. The text denotes the subject’s condition, and the color denotes larger condition groups. Large dots and text denote the centroid of that disease group. Only conditions with greater than three participants have a centroid shown. AI, autoinflammatory diseases; Telo, telomere disorders; PID, primary immunodeficiencies.

Extended Data Fig. 3 |. Characteristics of the individual and joint PCs from the JIVE analysis.

a, Similar to Fig. 2d, but with subjects projected onto the transcriptomic individual PC (iPC) 1 vs. iPC2 space. b, Similar to (a) but showing the serum protein iPCs. c, Gene set enrichment of transcriptomic (left) and serum protein (right) features positively correlated with jPC1 (enrichment calculated using the CameraPR function in limma; genes/proteins ranked by the Spearman correlation with the jPCs). Gene sets from KEGG pathways, GO biological process gene sets, Reactome pathways, and the blood transcriptomic modules and Human Protein Atlas tissue gene sets. d, Similar to (c) but for negatively correlated features. e, Scatterplot of a hematopoietic composite score (see Methods) vs. jPC2. Left panel displays the trend across all conditions including healthy controls and the right set of panels focus on individual disease groups whose clinical presentation may include marrow failure or lymphopenia. Inset focuses on GATA2 patients, highlighting those with abnormal bone marrow biopsies. Spearman correlation and associated unadjusted P values are shown. f, Scatterplot of Median Absolute Deviation (MAD) of jPC1 and jPC2 scores for each condition in the study. A higher MAD corresponds to greater variation within a disease for that jPC. G2BMD, GATA2 deficiency-associated bone marrow disorder; MDS, myelodysplastic syndrome.

Extended Data Fig. 4 |. Supporting data for the development and characterization of the IHM.

a, ROC curves for RF classifiers from leave-one-out-cross-validation using temporally stable features of individual or the indicated combinations of data modalities. b, ROC curve for the RF classifier (trained on all data modalities in the primary dataset) applied to the set of unseen, independent set-aside patients and healthy participants. c, Negative log₁₀ adjusted P values of feature GVI in each RF classifier. P values were determined through permutation (see Methods). Labels are shown for parameters passing an FDR cutoff of 0.2 for each classifier. FDR adjustment was performed on p values for parameters within a classifier. Features used in each classifier are shown on x-axis. d, Gene set enrichment of transcriptional surrogate signatures for the PM2 protein module and serum amyloid A (SAA) measurement. e, Scatterplots with regression lines of 182 participants’ IHM scores vs. the first three JIVE PC scores from the joint (jPCs) along with transcriptomic and proteomic individual iPCs. Pearson correlation and associated unadjusted P values are shown. f, Overview of the strategy for assessing the differential expression of the IHM transcriptional surrogate between healthy and autoimmune diseases in each of 28 cell types from Ota et al. g, Comparison between autoimmune diseases and healthy for the IHM and jPC1 transcriptional surrogate signature scores in 28 cell types from Ota et al. A negative effect size represents a decrease in the signature score in individuals with autoimmune disease relative to healthy. Effect sizes and adjusted P values are derived from linear regression models (see Methods). h, Box plots of IHM transcriptional surrogate signature scores comparing healthy controls vs. disease subjects from Ota et al. highlighting selected cell types from (g): classical monocytes (n = 333 patients and 77 healthy controls), neutrophils (n = 328 patients and 78 healthy controls), and plasmacytoid dendritic cells (pDC; n = 330 patients and 76 healthy controls). Effect size (Δ) and adjusted P value from linear regression models are shown. Box plot center lines correspond to the median value; lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles), and lower and upper whiskers extend from the box to the smallest or largest value correspondingly, but no further than 1.5X inter-quantile range. BTM, blood transcriptomic modules; CBC, complete blood count; CD8⁺ TEMRA, CD8⁺ T effector memory CD45RA⁺ cells; NK, natural killer; RDW, red cell distribution width; TBNK, lymphocyte (T, B, NK cell) phenotyping.

Extended Data Fig. 5 |. Supporting data for assessing the IHM.

a, Scatterplot with trendline showing the age dependence of jPC1 in healthy individuals (n = 34) only. Pearson correlation and P values shown. b, Overlap between proteins in the IHM protein surrogate signature and those reported in the Baltimore Aging Study to be significantly associated with chronological age (odds ratio and P value from one-sided Fisher’s exact test). c, Relationship between the IHM protein surrogate score and serum relative IL-6 level (as measured by Somalogic) in the Baltimore Aging study (Spearman correlation and associated P value shown; n = 240 participants). d, Similar to (c) but in the 34 healthy participants of the monogenic cohort. e, Correlation between age and the IHM protein surrogate score after regressing out CRP in the Baltimore Aging Study. Pearson correlation and associated P value shown. f, Scatterplots showing association between the Somalogic relative serum level of CXCL9/monokine induced by gamma (MIG) and the IHM in the healthy participants (left; n = 148) and patients only (right; n = 34) in the monogenic cohort (with Spearman correlation and P values shown). g, The IHM was re-derived without including PM2 (which contains CXCL9/MIG and correlated proteins). Scatterplot shows the correlation between age and this alternative IHM (without PM2) in the 34 healthy participants (with Pearson correlation and P values shown). h, Scatterplots showing the correlation between the IHM (x-axis) and the IMM-AGE (y-axis) transcriptomic surrogate scores in our monogenic cohort (training set on left, holdout set on right). Pearson correlations and P values shown. I, Classification performance for distinguishing healthy individuals from patients in the holdout set of the monogenic cohort using the IHM and IMM-AGE transcriptomic surrogate scores. The area under the curve (AUC) is listed next to the corresponding signature. j, Using a cell type-specific gene expression dataset with 10 common autoimmune diseases (Ota et al.), box plots show the AUC achieved by the IHM and IMM-AGE transcriptomic surrogate scores as well as pan-disease signature scores from Tuller et al. in distinguishing patients from healthy controls. Each circle represents one of the 28 cell types/subsets, with lines connecting the same cell subset across signatures. Unadjusted P values from paired, two-sided Wilcoxon test are shown. Box plot center lines correspond to the median value; lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles), and lower and upper whiskers extend from the box to the smallest or largest value correspondingly, but no further than 1.5X inter-quantile range. k, Alternative representation of the classification performance shown in (j) for classical monocytes (left) and TEMRA CD8⁺ T-cells (right). The AUC is listed next to the corresponding signature.

Extended Data Fig. 6 |. Supporting data for assessing the IHM.

a, Scatterplot showing the Spearman correlation of serum proteins with the IHM transcriptional surrogate signature within 34 healthy individuals (x-axis) vs. 154 disease individuals (y-axis) from the monogenic cohort. The names of the 20 proteins with the highest absolute correlations on the x- or y-axes are shown. b, Left: graphical overview of the analyses that assess whether correlation between individual serum proteins and IHM transcriptional surrogate scores depends on age. Right: Scatterplot showing the Spearman correlation values of serum proteins with the IHM transcriptional surrogate signature, with raw correlation coefficients on the y-axis and partial correlations accounting for age on the x-axis. The names of the 20 proteins with the highest absolute correlations on the x- or y-axes are shown. Correlations were computed with 34 healthy participants only. c, Similar to (b) but showing the correlation and partial correlation computed in patients only (n = 154). d, Scatterplot of IHM transcriptional surrogate signature score vs. Neurotrophin-3 in healthy controls (n = 34) from this study. Spearman correlation and associated are P value shown.

Extended Data Fig. 7 |. Supporting data for validating the IHM in external datasets.

a, Similar to Fig. 4d, but for a selection of transcriptional surrogate signatures capturing the status of the indicated parameters and expanded to also show the effect sizes and associated standard errors in each of the 21 studies in the meta-analysis. Size of circles indicates the relative sample numbers of each study. Summary meta-effect sizes are shown at the bottom. b, Assessment of whether genes in a given transcriptional surrogate signature had significantly lower P values in the meta-analysis results compared with genes not in the signature. Significance is determined by two-sided Wilcoxon rank sum test. c, Similar to Fig. 4g (showing the IHM protein surrogate score), but here showing the IHM transcriptional surrogate score. d, Classification performance measured by area under the curve (AUC) of the IHM surrogate scores (orange), several published aging scores (blue; see Extended Data Table 2) and individual components of the IHM (red) in distinguishing RA patients from healthy controls in the RA-MAP cohort. Unadjusted P values are from bootstrap sampling comparing IHM protein surrogate score and respective metric. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001, NS = not significant. e, Effect sizes and associated unadjusted P values in distinguishing RA patients from healthy controls in the RA-MAP cohort for each of the indicated metrics (similar to (d)) in a linear model including CRP as a covariate (see Methods). f, Similar to Fig. 4i, but for studies that included participants less than 50 years old (n = 630 across 18 studies). g, Similar to Fig. 4h, but with age regressed out of the IHM transcriptional surrogate score prior to performing the meta-analysis. h, Box plots of IHM protein surrogate scores for healthy participants and heart failure patients. Cohort 1: subjects were recruited at baseline, sampled, and then followed longitudinally. Some developed heart failure in the future (‘New onset HF’; n = 185) while some did not (‘Non-HF’; n = 583). IHM scores were derived from initial, baseline samples. Cohort 2: Participants (n = 85) recruited from an active heart failure clinic. Unadjusted P values are from two-tailed Wilcoxon tests. Box plot center lines correspond to the median value; lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles), and lower and upper whiskers extend from the box to the smallest or largest value correspondingly, but no further than 1.5X inter-quantile range. i, Relationship between IHM protein surrogate score and BMI from a cohort of 745 sedentary adults. Spearman correlation and associated P value shown.

Fig. 1 |. Study and data overview.

a, Patient groups and data collected. Individual disease groups are shown in c. b, Conceptual overview of the study and analysis approaches. Both disease group centric (top-down) and individual participant-based (bottom-up) analyses were pursued. c, Breakdown of cohort by disease and sample type. Data are broken down into the number of ‘primary’ samples (equal to the number of participants analyzed in this study), participants reserved (‘set-aside’) for independent follow-up analyses (Methods) and samples from the primary participants (‘repeat’) collected at additional timepoints. d, Gene–gene correlation heatmap of whole blood transcriptomic data. Modules of correlated genes (or TMs; k = 12) are annotated by color at the top and left. Modules were created using all transcriptional features; however, only the temporally stable genes are shown in the heatmap (see f and g below). Only modules with significant enrichments are annotated. e, Similar to d but for serum protein data. Modules of correlated proteins (PMs; k = 10) are annotated by color at the top and left. The gray module contains a large, weakly correlated set of proteins. f, Distribution of the percentage of variance assigned to each covariate across all measured proteins (1,305) and transcripts (15,729). The transcriptomic (protein) data had 276 (271) samples with 62 (64) participants with repeated sampling. g, Percentage of variance assigned to each variable in the variance partition analysis, run across each TM (blue), serum PM (magenta) and CBC parameter (green). This analysis used participants with repeat samples collected at different timepoints. The CBC/TBNK data consisted of 271 samples with 63 participants with repeated sampling. AI, autoinflammatory diseases; Telo, telomere disorders; PID, primary immunodeficiencies; CBC/TBNK, CBC and lymphocyte (T, B, NK cell) phenotyping; TM, whole blood TMs; PM, serum PMs; WBC, white blood cell count; MCHC, mean corpuscular hemoglobin concentration; HGB, hemoglobin; PLT, platelet count; MCH, mean corpuscular hemoglobin; MCV, mean corpuscular volume; RBC, red blood cell count.

Fig. 2 |. Bottom-up integration of transcriptomic and serum protein personal immune profiles reveals an emergent axis of immune health.

a, Conceptual overview of JIVE analysis integrating whole blood transcriptome and serum protein data. JIVE was performed using the subject-level data (Methods; n = 188 participants who had both measurements). b, Variation explained by the joint (shared by both data types), individual data type and residual latent factors in JIVE analysis. c, Heatmaps showing Pearson correlation between jPCs (rows) and major peripheral immune parameters and module scores (columns). Correlation was computed using the subject-level data (n = 182 participants who had serum protein, whole blood transcriptomic and CBC/TBNK data). Abbreviations are the same as for Fig. 1g. *Adjusted P < 0.05, FDR adjustment performed across all comparisons. d, Projection of patients (n = 154) and healthy (n = 34) participants onto the jPC1 versus jPC2 space. Small text labels show the disease groups to which individual patients belong whereas large dots and text denote the centroid (mean jPC1 and jPC2 values) of the indicated disease group. Only conditions with greater than three participants have a centroid shown. Box plots show projections onto single PC dimensions with patients grouped by disease condition (jPC1 below the centroid plot; jPC2 to the right of the centroid plot). Each participant’s score is represented as a single point. Unadjusted P values from two-sided Wilcoxon test comparing each disease and the healthy group: *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. e, Box plot of jPC1 scores comparing 154 patients (all disease conditions combined) with 34 healthy participants (P value computed using two-sided Wilcoxon test). f, Scatterplot of jPC1 scores derived using all participants (x axis) versus jPC1 scores derived using patients only (y axis; left) or healthy participants alone (right). Pearson correlation and associated P value shown (same group of participants as in d and e). For all box plots, center lines correspond to the median value, lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles), and lower and upper whiskers extend from the box to the smallest or largest value correspondingly, but no further than 1.5 × inter-quantile range. WB, whole blood.

Fig. 3 |. Top-down supervised machine learning classification analysis independently reveals an IHM highly concordant with that from unsupervised analysis.

a, Conceptual overview of the supervised machine learning analysis of healthy participants versus patients using RF classifiers to obtain an immunological health measure known as the IHM. The number of temporally stable features used from each data modality is shown. Models were trained using the subject-level data (n = 182 participants with serum protein, whole blood transcriptomic and CBC/TBNK data). b, ROC curve for distinguishing healthy participants versus patients in the training samples using the approach shown in a. c, Bar plot of the −log₁₀ adjusted P values for features passing a 0.2 FDR significance cutoff (gray dashed line; P values estimated through permutation testing of GVI from the RF classifiers); these top features contributed to the classifier used to derive the IHM. Direction was determined as the sign of the average difference between heathy participants and patients from all disease groups. d, Correlation between IHM and jPC1 scores across participants. Least squares regression lines are included for healthy participants (n = 34) and patients (n = 148) separately with Pearson correlations and associated P values shown. e, IHM scores of individual participants grouped by condition. Unadjusted P values from two-sided Wilcoxon test comparing each disease to the healthy group: * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001. f, Similar to e, but here showing smoothed density of IHM scores for each of the groups with at least ten participants. g, Scatterplot with trendline showing the age dependence of the IHM in healthy individuals only (Pearson correlation and P values shown; n = 34 healthy participants in the training set with serum protein, whole blood transcriptomic and CBC/TBNK data). For all box plots, center lines correspond to the median value, lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles), and lower and upper whiskers extend from the box to the smallest or largest value correspondingly, but no further than 1.5 × inter-quantile range.

Fig. 4 |. Assessing the IHM in independent datasets.

a, Graphical depiction of the creation of a protein IHM surrogate signature (Methods) followed by evaluation of the IHM and healthy aging association in the independent Baltimore Aging Study. b, Scatterplot with trendline showing the negative correlation between chronological age and the IHM protein surrogate score in 240 healthy participants from the Baltimore Aging Study. Pearson correlation and associated P value are shown. The shaded area represents the 95% confidence interval (95% CI). c, Graphical depiction of the evaluation of surrogate IHM signatures in (from left to right): (1) meta-analysis of four common, nonmonogenic immunological diseases across 21 independent studies; (2) the DA in a pediatric SLE cohort; (3) meta-analysis comparing vaccination response across age groups in independent studies; and (4) the DA and response to treatment in an RA cohort. d, Meta-effect sizes of the average difference between autoimmune patients and healthy participants for the IHM and jPC1 transcriptional surrogate scores. **Unadjusted P < 0.01. e, IHM transcriptional surrogate signature scores stratified by DA in a pediatric SLE cohort. Lines connecting samples collected from different timepoints of the same individual (n = 43). Averaged scores are used for patients having multiple visits with the same DA classification. Unadjusted P values from paired, two-tailed Wilcoxon test are shown. f, Correlation between baseline (pre-treatment) IHM transcriptional surrogate score and DA (DAS28, DA score in 28 joints) of 221 RA patients in the RA-MAP study. Spearman correlation and associated P value are shown. The shaded area represents the 95% CI. g, IHM protein surrogate scores at baseline and 6 months after treatment for patients with RA classified by their treatment response status (n = 33 responders and 33 nonresponders). Significance of group differences is determined by a mixed-effect model correcting for age and gender. The horizontal gray dashed line denotes the median IHM score of healthy controls (n = 26). h, Forest plot (top) of effect sizes from the meta-analysis across 11 independent vaccination cohorts that included both younger (<50 years; n = 217) and older (≥50 years; n = 356) participants, testing whether the IHM transcriptional surrogate signature evaluated at baseline before vaccination was associated with antibody titer responses to vaccination. Position and size of individual squares denote the effect size and relative number of participants in the indicated study, respectively. The diamond at the bottom represents the overall meta-effect size combining evidence across the cohorts. Width of the diamond indicates the standard error of the meta-effect. The density plot (bottom) shows the age distribution of the younger and older groups across all studies listed in the forest plot. i, Similar to h but for studies that included participants at least 50 years old (n = 492 across 9 studies). For all forest plots, each point shows the estimated mean effect and error bars show the 95% CI (1.96 × standard error). The meta-analysis P value is derived using the MetaIntegrator R package (Methods). For all box plots, center lines correspond to the median value, lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles), and lower and upper whiskers extend from the box to the smallest or largest value correspondingly, but no further than 1.5 × inter-quantile range.