Predict Vaccine Response with Immune Landscapes Mapping

Comprehensive immunoprofiling of LAIV responses reveals distinct immunophenotypic groups

To define responder status to LAIV in 244 Gambian children (7), we focused on adaptive immune markers with paired baseline (day 0) and postvaccination (day 21) measurements, expressed as fold-change values (V21/V0; see Methods) (Figure 1A). Using fold-change accounted for interindividual variability in baseline immunity, capturing genuine vaccine-induced changes. We evaluated a comprehensive panel of antibody-mediated responses, including hemagglutination inhibition (HAI) titers, an indicator of antibodies that block the binding of the influenza virus to host cells (21). We used an influenza virus protein microarray to assess the breadth of antibody responses (22). This high-throughput platform profiles binding antibody responses across multiple influenza strains, including hemagglutinin (HA) proteins from various influenza A and B viruses. This allowed quantitative evaluation of serum antibody binding profiles before and after LAIV administration, providing insights into the specificity, magnitude, and breadth of the antibody responses, including cross-reactive responses. We also examined stalk-specific responses targeting conserved regions of the HA protein, including antibody-dependent cellular cytotoxicity activity measured against chimeric HA stalk constructs (e.g., cH6/1 and cH7/3) to assess cross-reactive immunity (23). Neuraminidase (NA) titers in blood and nasal mucosa offered insights into cross-protective responses (24). Complementing antibody profiles, we assessed T cell IFN-γ and IL-2 production upon stimulation with vaccine strain components (HA, NA, and matrix/nucleoprotein) to capture systemic cellular responses. Collectively, this panel of immunological assays provided a highly granular view of the magnitude and quality of immune responses elicited by LAIV administration, allowing us to capture a detailed immunophenotypic landscape.

Immune response landscape mapping of LAIV reveals distinct immunophenotypic groups. (A) Cohort overview depicting all features used for unsupervised machine learning analysis: 244 children (24–59 months of age) vaccinated with LAIV; mucosal and blood samples collected on day 0 (prevaccination) and day 21 (postvaccination). Vaccine-induced immune responses calculated as fold-change relative to prevaccination levels. (B) Workflow schematic for automated clustering pipeline applying t-SNE dimensionality reduction, KNN graph construction, and Louvain community detection to identify distinct immunophenotypic clusters. (C and D) Louvain resolution sweep results used to assess cluster stability and select optimal number of clusters. (C) Modularity score plotted against Louvain resolution parameter, colored by number of clusters identified (3–6). High modularity indicates well-separated clusters. Red diamond indicates selected clustering parameters. (D) Number of clusters identified plotted against Louvain resolution parameter, colored by modularity score. Stability of 3-cluster solution (red diamond) is observed across range where modularity is maximal (Q ≈ 0.717). (E) Clustered t-SNE plot of fold-change data (post/pre-LAIV) revealing 3 distinct LAIV response phenotypes: group 1 (green, n = 82), group 2 (orange, n = 88), and group 3 (purple, n = 74). (t-SNE parameters: perplexity: 30; exaggeration factor: 4; max iterations: 10,000; theta: 0; eta: 500; K: 60 for KNN graph; final silhouette score: 0.40). (F and G) Clustering patterns overlaid with demographic factors on t-SNE map. (F) Clustering by sex (female, green; male, orange). (G) Clustering by study year (2017, green; 2018, orange). (H) Heatmap and hierarchical clustering display fold-change data for key immune features across 3 clusters (columns: groups 2, 1, and 3 from left to right). Rows represent immune features, clustered using Euclidean distance and Ward’s D2 method. Heatmap cells are colored based on scaled FC values from –1 (blue, low FC) to 1 (red, high FC). The top color bar indicates responder groups (group 1, green; group 2, orange; group 3, purple). Side color bars indicate qualitative response classifications derived from assays: HAI (purple: high, dark; low, light), IgA (orange: high, dark; low, light), CD4+ T cell (blue: high, dark; low, light), and CD8+ T cell (green: high, dark; low, light). Column cluster ordering optimized for visual clarity.

This integrated, multimodal dataset served as input for the Immunaut machine learning framework (see Methods). To visualize patterns, we projected the high-dimensional data into a 2D space using t-distributed stochastic neighbor embedding (t-SNE) (Figure 1B). We then constructed a K-nearest neighbors (KNN) graph based on Euclidean distances in this reduced space. We applied the Louvain community detection algorithm to identify distinct immunophenotypic groups, which partitions the graph to maximize the modularity score (Q), a measure of clustering quality where a higher modularity score indicates more distinct and well-separated clusters. We systematically evaluated clustering stability by applying the algorithm across a range of resolution values (r), where lower resolutions yield fewer, larger clusters and higher resolutions produce more, smaller ones. This assessment revealed a resolution range where modularity reached a high and stable plateau (Q ≈ 0.717, Figure 1C), signifying a robust and well-defined community structure, and the number of clusters consistently converged to 3 (Figure 1D). This provides quantitative evidence that this partitioning reflects distinct biological subtypes rather than arbitrary divisions sensitive to parameter tuning.

The final 3-cluster partition is visualized on the t-SNE projection (Figure 1E), comprising group 1 (green, n = 82), group 2 (orange, n = 88), and group 3 (purple, n = 74). The average silhouette score of 0.4 indicates moderately distinct clusters. We observed no substantial enrichment of specific sexes (Figure 1F) or study years (Figure 1G) within any cluster, suggesting the clustering captures genuine immunophenotypic differences independent of these external biases known to affect immune responses to vaccines (25–27).

Individuals in group 1 (n = 82) displayed a distinct profile characterized by CD8+ T cell–mediated responses and notably low CD4+ T cell IFN-γ activity (Figure 1H). This group showed elevated IFN-γ and IL-2 production by CD8+ T cells upon stimulation, with the most pronounced responses against influenza B virus HA and matrix/nucleoprotein antigens (Figure 1H). Influenza B virus–specific CD8+ T cell IFN-γ responses were statistically significant compared with group 3 (Figure 2A). Conversely, humoral and IgA responses in group 1 were minimal or absent (Figure 1H and Figure 2, A and D). Although some N1-specific IgA responses were detected (Figure 1H), these responses were not statistically significant and were comparable to those observed in group 2, indicating that the N1 IgA responses were not a distinguishing feature of group 1 (Supplemental Figure 1A; supplemental material available online with this article; https://doi.org/10.1172/JCI189300DS1). Based on these results, we termed individuals in group 1 as “CD8+ T cell responders.”

Vaccine response immune signatures defining LAIV responder types. (A) Polar plot summarizing scaled median expression of key immune features in CD8+ T cell responders (group 1, green). CD8+ T cell responders are characterized by robust influenza B virus HA–specific CD8+ IFN-γ responses and limited humoral immunity, with median feature values represented in the polar plot and fold-change comparisons shown in the adjacent box plot. (B) Polar plot for mucosal responders (group 2, orange) illustrating strong mucosal IgA responses, particularly stalk-specific (cH7/3 IgA) and H3N2 virus HA–specific IgA antibodies and influenza B virus–specific responses. Box plots detail fold changes (shown as log10) for various immune features, highlighting systemic (influenza B virus HAI) and mucosal immune activation (IgA). (C) Polar plot depicting systemic, broad influenza A virus responders (group 3, purple), showing elevated systemic antibody responses to multiple influenza A virus strains (e.g., H1, H3), as well as cross-reactive IgG and antibody-dependent cellular cytotoxicity activity. Box plots show fold-change values (log10) for each immune marker across responder groups. (D) Integrated radar plot comparing scaled median immune expression profiles across all responder groups (CD8+ T cell responders in green, mucosal responders in orange, systemic broad influenza A virus responders in purple), emphasizing distinct immune feature distributions. This integrative visualization highlights the unique baseline and postvaccination immune landscapes that define each responder profile. Box plots denote minimum to maximum values, and points are all individuals within the group. *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001, by 1-way ANOVA Kruskal-Wallis test with Dunn’s multiple-comparison test to adjust for multiple testing.

In contrast, group 2 (n = 88) individuals exhibited a profile dominated by mucosal IgA responses (Figure 1H). This group showed statistically significant induction of mucosal IgA antibodies across all antigens and strains tested, including against chimeric HA stalk constructs, indicating the induction of antibodies targeting conserved HA regions (Figure 2B). The consistent IgA increase was unique to group 2 and not observed in groups 1 or 3 (Figure 2, B and D), validating their classification as mucosal responders. Group 2 also exhibited significant seroconversion to influenza B viruses (Figure 1H), evidenced by substantial increases in HAI titers (Figure 2B). Although antibody binding responses to influenza B viruses measured by influenza virus protein microarray were elevated (Figure 1H), they were not statistically significant (Supplemental Figure 1B). This suggests humoral immunity in group 2 included both mucosal and systemic antibody responses against influenza B virus.

Group 3 (n = 74) individuals showed robust systemic antibody responses to influenza A viruses (Figure 1H). This was evidenced by significant increases in HAI titers for both H1N1 and H3N2 strains (Figure 2C). The antibody responses demonstrated breadth, with substantial increases in binding to HA subtypes from multiple contemporary and historical H1N1 and H3N2 strains not present in the vaccine (Figure 2C). Elevated responses were also observed against the cH6/1 chimeric HA construct, including increased antibody-dependent cellular cytotoxicity (Figure 2C). Significantly higher N1 titers were also detected (Supplemental Figure 1C), supporting a coordinated response targeting conserved epitopes. Although some CD4+ T cell responses were elevated (Supplemental Figure 1C), significant IgA responses were absent, indicating predominantly systemic immunity (Figure 2, B and D). This finding supports classifying them as systemic broad influenza A virus responders.

Predictive modeling of LAIV response phenotypes based on baseline immune profiles

We next sought to determine whether prevaccination immune profiles could predict an individual’s response type. To achieve this, we used comprehensive baseline immunological measurements before vaccination (Figure 3A), including antibody profiles, T cell responses, Streptococcus pneumoniae load (pneumococcal carriage density), asymptomatic respiratory viral presence, RNA pathway scores from nasal samples, and frequencies of various immune cell subsets such as monocytes, plasmacytoid and myeloid DCs (pDCs and mDCs), and T follicular helper (Tfh) cells (28) relevant to LAIV, which relies on both innate and adaptive immune pathways to induce protection (5, 6, 29, 30).

Automated machine learning framework for mapping and predicting LAIV immunogenicity response phenotypes. (A) Overview of the automated machine learning framework developed to predict LAIV response phenotypes using baseline immune data from mucosal and blood samples, capturing multidimensional immune parameters such as transcriptomics, antibody titers, bacterial load, flu-specific T cell responses, and comprehensive immunophenotyping. (B) Step 1, balanced data partitioning: the dataset is split into training (80%) and testing (20%) sets, ensuring proportional representation of each immunophenotypic group (CD8+ T cell; mucosal; and systemic, broad influenza A responders) to maintain predictive accuracy across classes. Step 2, model optimization cycle: 10-fold cross-validation and hyperparameter tuning are applied across 141 machine learning models, each iteratively trained and validated to identify the best predictors of vaccine response. Step 3, model evaluation and scoring: predictive performance metrics, including specificity, sensitivity, and AUC, are calculated on the test set (20%) for model validation. Feature importance scores are computed for each baseline variable, providing a ranked analysis of each immune parameter’s contribution to LAIV response prediction. (C) Multiclass ROC plot of the gradient boosting machine model evaluated on the test set (20%), displaying predictive accuracy across all 3 classes: CD8+ T cell responders (green); mucosal responders (orange); and systemic, broad influenza A responders (purple) in a one-versus-all comparison. (D) Variable importance score table for the gradient boosting machine model, showcasing the cumulative importance of the selected baseline features across the 3 predicted classes, highlighting the most influential parameters in LAIV immunogenicity prediction.

To model the mapped vaccine responses, we applied the Sequential Iterative Modeling OverNight (SIMON) platform, which is designed for high-dimensional datasets with substantial interindividual variability (17, 31) (Figure 3B). We systematically tested 141 machine learning algorithms to ensure the selection of the most accurate and biologically meaningful model (17, 31–33). We employed 10-fold cross-validation during model training to enhance robustness and mitigate overfitting, and we assessed performance on a held-out test set to ensure generalizability. Out of the 141 models tested, 26 achieved a test set area under the receiver operating characteristic curve (AUC) above 0.7, underscoring the predictive strength of our baseline immune profiles (Supplemental Table 1).

Among all models, the gradient boosting machine model was the top performer (Supplemental Table 1). It achieved an accuracy of 59.57% (exceeding the null accuracy of 36.17%; P = 0.0009), a balanced accuracy of 71.67%, an F1-score of 0.6286, a precision of 0.6902, and a recall of 0.6471, highlighting its capacity to balance false-positives and false-negatives and an overall AUC of 0.8. One-versus-all AUCs confirmed robust performance across individual classes: 0.80 for CD8+ T cell responders, 0.77 for mucosal responders, and 0.73 for systemic broad influenza A virus responders (Figure 3C). Training gradient boosting machine models on individual or pairwise data types showed that removing any primary data modality reduced performance, demonstrating that integration of diverse features was essential for high accuracy (Supplemental Table 2). The gradient boosting machine model’s capacity for feature importance estimation, its ability to manage high-dimensional data, and its robustness to missing data (Supplemental Figure 2) make it a powerful tool for this classification task.

Next, we identified the baseline features that were most critical for classification (Figure 3D). The top predictor was the baseline HAI geometric mean titer against H3N2 (score 100), indicating that preexisting systemic immunity drives the response type. However, high baseline mucosal IgA against various influenza antigens, including influenza B/Victoria/2/87-like lineage HA and NA (61), pH1N1 HA (48), N1 (20), H3N2 NA (18), and cH7/3 IgA (44) was also pivotal, underscoring the complementary roles of systemic and mucosal immunity. Key cellular parameters, such as IFN-γ–producing T cells (e.g., influenza A virus matrix/nucleoprotein CD4 IFN-γ, score 64; H1N1 and H3N2 HA CD4 IFN-γ, 33 and 36; influenza B HA CD8+ IFN-γ, 21) and Tfh cell frequencies (34), also surfaced as prominent predictors. Additionally, innate immune cells, pneumococcal carriage density (36), and asymptomatic respiratory viral infection (6) emerged as critical modulating factors. Notably, baseline nasal RNA-derived Gene Ontology (GO) pathways, encompassing metabolism (GO:0072521, variable importance score 46), morphogenesis (GO:0060562, score 40), and Hedgehog signaling (GO:0008589, GO:0007224, scores 22 and 6), contributed substantially, pointing to a context-dependent model of vaccine responsiveness where tissue-level processes shape the immune response.

Collectively, these observations suggest that LAIV response phenotypes arise not from a single dominant factor but emerge from a finely tuned network of systemic and local immunity, innate and adaptive cellular components, and underlying tissue-level processes.

Identifying preexisting immune landscapes that shape LAIV responses

To delineate the preexisting immune landscapes that define each group, we hypothesized that specific baseline conditions characterize each responder class. To test this, we combined machine learning–derived insights with exploratory analyses of baseline seropositivity, viral shedding, and detailed immunological profiling. The resulting patterns suggest that historical exposure to influenza strains plays a pivotal role in shaping the immune response to LAIV (Figure 4).

Baseline immune landscape and viral shedding profiles predictive of LAIV response groups. (A) Heatmap of baseline immune features predictive of LAIV response groups, organized by hierarchical clustering to show feature relationships and variations across groups (Euclidean distance, Ward’s D2 clustering method). Each cell reflects a scaled expression level, with red representing high expression and blue indicating low expression, revealing the distribution of immune features at baseline across the identified immunophenotypic clusters. (B) The proportion of seropositive children (HAI titer ≥10) at baseline (before vaccination) within each responder group and across all 3 LAIV-strains, pH1N1, H3N2, and influenza B virus. (C) The proportion of children that shed LAIV strains (pH1N1, H3N2, and B) on day 2 and day 7 after vaccination across all 3 responder groups. (D) Box plots showing baseline features, including H3N2 HAI geometric mean titer (gmt), titer of antibodies binding H3 HA from A/Switzerland/9715293/2013 analyzed by influenza virus protein microarray (H3 HA SWISS IVPM), titer of antibodies binding NA from group 2 (N2) and frequency of influenza B virus HA–specific CD8+ T cells producing IFN-γ across all 3 responder groups, CD8+ T cell responders (green); mucosal responders (orange); and systemic, broad influenza A virus responders (purple). Box plots denote minimum to maximum values, and points are all individuals within the group. *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001, by 1-way ANOVA Kruskal-Wallis test with Dunn’s multiple-comparison test to adjust for multiple testing. (E) Forest plots showing log-odds estimates from a logistic regression model. The plots illustrate the association of the mucosal responder (orange) and systemic, broad influenza A virus responder (purple) groups with the outcomes of viral shedding (day 2 and 7) and HAI seropositivity, relative to the CD8 T-cell responder group which serves as the reference category. The analysis is stratified by LAIV strain (H3N2, pH1N1, and B), and the error bars represent the confidence intervals for the log-odds estimates.

Children who became CD8+ T cell responders had a distinctive baseline signature (Figure 4A). Before vaccination, this group had significantly higher seropositivity for H1N1 (48%) and H3N2 (72%) (P = 0.049 and P < 0.0001) and elevated baseline HAI responses (Figure 4, A and B). They also had elevated baseline levels of influenza virus–specific IgA in nasal secretions, targeting multiple LAIV strains and cH7/3 chimeric stalk construct, features identified by the machine learning model as important (Figure 4A). After vaccination, this group showed significantly reduced shedding of H3N2 by day 7 (17% shedding rate), an association confirmed by logistic regression (β = 1.21, P = 0.0078) (Figure 4E). Baseline nasal transcriptional analysis identified enrichment of pathways, including purine metabolism (GO:0072521; score 46) and regulation of defense response (GO:0031347; score 16) (Figure 4A). These children also had elevated baseline S. pneumoniae loads (score 36), Hedgehog signaling pathway (GO:0007224, score 6), and asymptomatic respiratory viruses (predominantly rhinovirus) detected before vaccination (χ2 test P = 0.57; Supplemental Figure 3). The frequency of circulating classical monocytes and mDCs was also elevated at baseline (Figure 4A).