to check the extent to which covariates correlated with each other, thereby possibly acting as confounders in the principal component regression scores, we determined dependence between all covariate pairs for every cell type. if at least one covariate was continuous, we calculated the fraction of variance in the continuous covariate that could be explained by the other covariate ( dummying categorical covariates ) and took the square root ( equal to pearson ’ s r for two continuous covariates ). for two categorical covariates, if both covariates had more than two unique values, we calculated normalized mutual information between the covariates using scikit - learn, since a linear regression between these two covariates is not possible. to control for spurious correlations between interindividual cell type variation and covariates due to low sample numbers, we assessed the relationship between sample number and mean variance explained across all covariates for every cell type. we found that for cell types sampled in fewer than 40 samples the mean variance explained across all covariates showed a high negative correlation with the number of samples ( supplementary fig. 4a ). we reasoned that for these cell types correlations between interindividual variation and our covariates were inflate d due to undersampling. moreover, we note that at lower sample numbers technical and biological covariates often strongly correlate with each other across donors ( supplementary fig. 4c ). this might lead to the attribution of true biological variation to technical covariates, and vice versa, complicating the interpretation of observed interindividual cell type variation. consequently, we consider 40 a recommended minimum number of samples to avoid spurious correlations between observed interindividual variation and tested covariates, and excluded results from cell types with fewer samples. to select cell types for which covariate effects could be confidently modeled at the gene level, we followed the same procedure for every cell type: we filtered out all genes that were expressed in fewer than 50 cells and all samples that had fewer than ten cells of the cell type. we furthermore filtered out datasets with fewer than two donors and refrained from modeling categories in covariates that had fewer than three donors in their category for that cell type. we encoded smoking status as a continuous covariate, setting never to 0, former to 0. 5 and current to 1. anatomical region was encoded into anatomical region ccf scores as described earlier. as we noted that changes from the nose to the rest of the airways and lungs were often independent from continuous changes observed in the lungs only, we encoded nasal as a separate covariate, setting samples from the nose to 1 and all others to 0. bmi and age were rescaled, such that the 10th percentile of observed values across the atlas was set to 0 and the 90th percentile was set to 1 ( 25 and 64 years for age, respectively, and 21. 32 and 36, 86 for bmi ). to determine whether covariance between covariates was low enough for modeling, we calculated the variance inflation factor ( vif ) between covariates at the donor level. the vif quantifies multicollinearity among covariates of an ordinary least squares regression and a high vif indicates strong linear dependence between variables. if the vif was higher than 5 for any covariate for a specific cell type, we concluded that covariance was too high and excluded that cell type from the modeling. as many cell types lacked sufficient representation of harmonized ethnicities other than european, including harmonized ethnicity in the analysis simultaneously decreased the samples that could be included in the analysis to only those with ethnicity annotations; hence, we excluded harmonized ethnicity from the modeling. to model the effects of demographic and anatomical covariates ( sex, age, bmi, harmonized ethnicity, smoking status and anatomical location of the sample ) on gene expression, we employed a generalized linear mixed model. we used sample - level pseudo - bulks ( split by cell type ), since the covariates modeled also varied at the sample or donor level and not at the cell level. modeling these covariates at the cell level ( that is, treating single cells as independent samples even when they come from the same sample ) has been shown to inflate p values. first, we split the lung cell atlas by cell type annotation, pooling detailed annotations into one subtype ( for example, grouping all lymphatic endothelial cell subtypes into one ) ( supplementary table 5; includes the number of cells per cell type ). subsequent filtering, covariate encoding and exclusion of cell types due to covariate dependence are described above. gene counts were summed across cells for every sample, within cell type. sample - wise sums ( that is, pseudo - bulks ) were normalized using edger ’ s calcnormfactors function, using default parameter settings. we then used voom, a method designed for bulk rna - seq that estimates observation - specific gene variances and incorporates these into the modeling. specifically, we used a voom extension ( differential expression testing with linear mixed models ) that allows for mixed - effects modeling and modeled gene expression as: [ 12pt ] $ $ } [ } } ] 1 + } + } + } + } + } + } \, } \ \ + ( } } ) $ $ log normcount ~ 1 + age + sex + bmi + smoking + nose + ccf score + 1 ∣ dataset where dataset is treated as a random effect to correct for dataset - specific changes in expression and all other effects are modeled as fixed effects. resulting p values were corrected for multiple testing within every covariate using the benjamini – hochberg procedure. to identify more systematic patterns across genes and changes happening at the gene set level, a gene set enrichment analysis was performed using correlation - adjusted mean - rank gene set tests. the analysis was performed in r using the camerapr function in the limma package, with the differential expression test statistic. gene ontology biological process terms were tested separately for each comparison. these sets were obtained from msigdb ( version 7. 1 ), as provided by the walter and eliza hall institute ( https: / / bioinf. wehi. edu. au / msigdb / index. html ). to stratify gwas results from several lung diseases by lung cell type, we applied stratified linkage disequilibrium score regression in single cells ( sc - ldsc ). sc - ldsc can link gwas results to cell types by calculating, for each cell type, whether disease-associated variants are enriched in genomic regions of cell - type specific genes ( i. e. the region of each gene and its surrounding base pairs ), while taking into account the genetic signal of proximal linkage disequilibrium-associated regions. here cell - type specific genes are defined as genes differentially expressed in the cell type of interest. in contrast with simple enrichment testing of only significantly disease-associated genes from a gwas among genes differentially expressed in a cell type, this method takes into account all snps included in the gwas. thus, consistent enrichment of weakly disease-associated genes ( that would not individually pass significance tests ) in a cell type could still lead to a significant association between the disease and the cell type. in this way, sc - ldsc provides more statistical power to detect associations between cell types and heritable phenotypes such as lung diseases. to perform sc - ldsc on the hlca, first a differential gene expression test was performed for every grouped cell type ( supplementary table 5 ) in the hlca using a wilcoxon rank - sum test, testing against the rest of the atlas. the top 1, 000 most significant genes with positive fold changes were stored as genes characterizing that cell type ( cell type genes ) and used as input for ldsc.