High-dimensional profiling clusters asthma severity by lymphoid and non-lymphoid status

Abstract

Clinical definitions of asthma fail to capture the heterogeneity of immune dysfunction in severe, treatment-refractory disease. Applying mass cytometry and machine learning to bronchoalveolar lavage (BAL) cells, we find that corticosteroid-resistant asthma patients cluster largely into two groups: one enriched in interleukin (IL)-4⁺ innate immune cells and another dominated by interferon (IFN)-γ⁺ T cells, including tissue-resident memory cells. In contrast, BAL cells of a healthier population are enriched in IL-10⁺ macrophages. To better understand cellular mediators of severe asthma, we developed the Immune Cell Linkage through Exploratory Matrices (ICLite) algorithm to perform deconvolution of bulk RNA sequencing of mixed-cell populations. Signatures of mitosis and IL-7 signaling in CD206^-FcεRI⁺CD127⁺IL-4⁺ innate cells in one patient group, contrasting with adaptive immune response in T cells in the other, are preserved across technologies. Transcriptional signatures uncovered by ICLite identify T-cell-high and T-cell-poor severe asthma patients in an independent cohort, suggesting broad applicability of our findings.

Keywords: BAL; CyTOF; FceRI+; ICLite; IFN-g; RNA-seq; clusters; immune; multi-omics; severe asthma.

Figure 1.. Unsupervised clustering of BAL cells identifies immune populations

(A) Results of FlowSOM projected on t-stochastic neighbor embedding (t-SNE) space. (B) Relative staining intensity for indicated surface markers across t-SNE space. (C) Summary of cell types identified by FlowSOM. (D) Results of patient clustering projected on principal component analysis (PCA) plot. Patients are colored by group, with ellipses representing 90% of the confidence interval around group centroid. (E) Distribution of clinical disease severities across groups is represented by stacked bar chart of proportion with p value determined by Pearson’s chi-square testing of raw values. (F) Boxplot of FEV1 measured by spirometry across PGs with p value of variance calculated using Kruskal-Wallis testing. (G) Elastic net (EN) predicted lung function based on high-dimensional cell count versus measured values of FEV1% predicted. Gray area indicates the 95% confidence bounds around a linear regression model comparing the two. Spearman’s rho and p value are indicated in plot area. (H) Graphical representation of EN modeling of cellular determinants of lung function (FEV1). Coefficients are plotted in order of ascending value from left to right, with distance from the hashed line indicating magnitude of contribution to the model. Blue coloration of cluster ID denotes a negative coefficient, and red indicates positive.

Figure 2.. Patient groups are defined by divergent cell lineages

(A) Stacked bar plot of BAL immune cell composition of cohort participants from manual hierarchical gating, arranged by PG. (B) Boxplot of cell lineages differentially presented across PGs with p value of variance calculated using Kruskal-Wallis testing. Bars represent median values, with bounds of boxes representing interquartile range (IQR) and whiskers representing 1.5× the upper or lower IQR. (C) Network of correlated cell types across BAL of all patients in cohort. Nodes represent cells, and edges represent a Spearman correlation rho ≥0.45 with p < 0.05. Coloring of nodes for PG specificity is based on Dunn post hoc testing of cell lineages identified as variant by Kruskal-Wallis.

Figure 3.. Cellular immune phenotype relates to distinct cytokine expression patterns

(A) Graphical summary of cell lineages in t-SNE space. (B) Relative staining intensity of six cytokines across t-SNE-reduced space. (C) Pie charts demonstrating distribution of cytokine-positive cells in T-cell and non-T-cell compartments. (D) Density plots for IL-10 staining intensity of indicated cells demonstrate the relative distribution of events for a PG. Hashed line indicates the 85th quantile. Numeric values presented in plot area indicate the percentage of events in each PG falling above this cutoff. Cell types presented have passed Kolmogorov-Smirnov testing for variance in distribution with p value < 1e–10. Boxplots demonstrate cytokine-positive cells per million BAL cells per patient in respective groups with p value of variance calculated using Kruskal-Wallis testing. Bars represent median values, with bounds of boxes representing IQR and whiskers representing 1.5× the upper or lower IQR. (E) Plotting for IL-4 within indicated cell lineages as described in (D). (F) Plotting for IL-5 within indicated cell lineages as described in (D). (G) Plotting for IFN-γ within indicated cell lineages as described in (D).

Figure 4.. WGCNA of BAL links global transcriptional signatures to specific cell lineages

(A) GSEA using Hallmark and curated gene sets. The p values, normalized enrichment score (NES), and enrichment scores (ES) are listed in the figure. (B) Heatmap of correlation between WGCNA module eigengenes and immune cell log ratios from mass cytometry evaluation of BAL. The p value and Spearman’s rho are indicated in the figure. (C) Barplot of −log10(p values) for gene ontology (GO) term enrichment of WGCNA modules of interest based on correlation to cells from BAL. (D) Plotting pink module absolute gene significance (GS) for correlation to CD206⁺ CD11b⁺ CCR4⁺ macrophages versus module membership (MM). (E) Visualization of pink module network hub genes (diamonds) with next closest network members (circles). (F) Plotting magenta module GS for correlation to FcεRI⁺ CD127⁺ CCR4⁻ innate cells versus MM. (G) Visualization of magenta module network hub genes (diamonds) with next closest network members (circles). (H) Plotting green module GS for correlation to CD8 or CD4 T cells versus MM. (I) Visualization of green module network hub genes (diamonds) with next closest network members (circles).

Figure 5.. ICLite deconvolution of BAL sequencing improves GO term linkage

(A) Correlogram of gene module to cell cluster interactions using the ICLite algorithm. Spearman correlation of module scoring with immune cell log ratio was used construct a correlogram of BAL gene modules (x axis) versus cell lineages (y axis). Only associations with false discovery rate (FDR) corrected p < 0.05 are illustrated. (B) Barplot of −log10(p values) for GO term enrichment of ICLite modules 13 and 15, which correlate with CD206⁺ CD11b⁺ CCR4⁺ macrophages. Plotting of patient module score versus immune cell log ratio is adjacent, with Spearman’s rho and p value of comparison as indicated in the figure. Red hashed line indicates linear regression model for data. (C) Barplot of −log10(p values) for GO term enrichment of ICLite modules 10 and 11, which correlate with FcεRI⁺ CD127⁺ CCR4⁻ innate cells and CD206⁺ FcεRI⁺ macrophages, respectively. Plotting of respective patient module score versus immune cell log ratio is adjacent, with Spearman’s rho and p value of comparison as indicated in the figure. Red hashed line indicates linear regression model for data. (D) Barplot of −log10(p values) for GO term enrichment of ICLite modules 1 and 20, which correlate with CD4 and CD8 EMs, CMs, and TRMs. Plotting of respective patient module score versus CD4 TRM log ratio is adjacent, with Spearman’s rho and p value of comparison as indicated in the figure. Red hashed line indicates linear regression model for data.

Figure 6.. ICLite ascribes greater breadth of functional terms to BAL cells than WGCNA

(A) Circos plot demonstrating gene module membership across trait-association technologies. WGCNA modules (bottom) are connected to ICLite modules based off commonality in membership, with chords of diagram colored by WGCNA module. (B) Phylogram of GO term semantic similarity for ICLite modules. Distance is independent of cell associations and based only on functional enrichment from transcriptional data. Module coloring is based off hierarchical clustering of semantic similarity. (C) Barplots quantifying the total number of GO terms enriched in modules or the number of GO terms effectively linked to cells by respective technology. (D) Tree plot of GO term semantic clustering results for all modules effectively linked to BAL cells using either ICLite or WGCNA. Size of box and text indicates −log10(p value) of enrichment. Families are colored according to parent-child relationship in term clustering.

Figure 7.. Gene modules predict cellular phenotype in SARP cohort

(A) Schematic for machine learning model of patient classification in external cohort using training based on cell count and transcriptional profile (B) Receiver operating characteristics (ROC) curve of a sparse-partial least-squares discriminant analysis (sPLS-DA) model for cellular immune phenotype prediction using BAL gene module scoring within the IMSA (experimental) cohort. ROC curves were calculated as one class versus the others using 5-fold validation on the original training set. Reported area under the curve (AUC) values (right of plot) are based on comparison of predicted scores of one class versus the others using a two-component model. Wilcoxon test of predicted scores for one class versus the others met a significance threshold of p < 0.05 for all groups. (C) Distribution of clinical disease severities across predicted SARP PGs is represented as proportion in stacked bar chart with overall p value determined by Pearson’s chi-square testing of raw values. (D Boxplot of manual differential cell counts from BAL across predicted SARP PGs with p value of variance calculated using Kruskal-Wallis testing and represented on plot. Bars represent median values with bounds of boxes representing IQR and whiskers representing 1.5× the upper or lower IQR. (E) Plotting of patient module scores versus percent lymphocyte of participant BAL from the SARP cohort, with Spearman’s rho and p value of comparison as indicated in the figure. Red hashed line indicates linear regression model for data. (F) Boxplot of measured spirometry across predicted SARP PGs with p value of variance calculated using Kruskal-Wallis testing and represented on plot. (G) SARP BAL gene expression data was used for GSEA using Hallmark and curated gene sets. The p values, NES, and ES for gene sets enriched in predicted SARP PG2 are listed in the figure. (H) SARP BAL gene expression data were used for GSEA. The p values, NES, and ES for gene sets enriched in predicted PG3 are listed in the figure.