Bronchial brushings were reconstructed in silico from the single-cell data by taking the sum of all transcript counts for each gene across all cells procured from each donor. Negative binomial generalized linear models were built using the MASS R package (v7.3-45), modeling transcript counts as a function of smoking status (FDR q < 0.05: n = 593 genes). In parallel, using never and current smoker bulk bronchial brushing microarray data (GEO series GSE7895), linear models were built using the stats R package (R v3.2.0), modeling gene-level expression values as a function of smoking status (FDR q < 0.05: n = 689 genes). The correlation between test statistics generated from both models was then measured to compare differential expression results (fig. S4A). Using the overlap among smoking-associated genes identified in both models (n = 155 genes), correlations (Spearman) among in silico bronchial brushings and bulk bronchial brushings were examined (fig. S4B).