then, genes that were highly variable in all datasets were considered overall highly variable, followed by genes highly variable in all datasets but one, in all datasets but two and so on until a predetermine d number of genes were selected ( 2, 000 or 6, 000 genes ). for scanvi and scvi, genes were subset to the hvg set and the resulting raw count matrix was used as input. for all other methods, scran - normalized ( as described above ) data were used. genes were then subset to the precalculated hvg sets. for integration of gene - scaled data, all genes were scaled to have mean of 0 and standard deviation of 1. two integration methods allowed for input of cell type labels to guide the integration: scgen and scanvi. as labels, level 3 harmonized cell type labels were used ( supplementary table 4 ), except for blood vessel endothelial, fibroblast lineage, mesothelial and smooth muscle cells, for which we used level 2 labels. since scgen does not accept unlabeled cells, cells that did not have annotations available at these levels ( that is, cells annotate d as cycling, epithelial, stromal or lymphoid cells with no further annotations; 4, 499 cells in total ) were left out of the benchmarking data. the dataset rather than the donor of the sample was used as the batch parameter. the maximum memory usage was set to 376 gb and all methods requiring more memory were excluded from the analysis. the quality of each of the integrations was scored using 12 metrics, with four metrics measuring the batch correction quality and eight metrics quantifying the conservation of biological signal after integration ( supplementary fig. 1; metrics previously described ). overall scores were computed by taking a 0. 4: 0. 6 weighted mean of batch effect removal to biological variation conservation ( bioconservation ), respectively. methods were ranked based on overall score ( supplementary fig. 1 ). for integration of the data into the hlca core, we first determined for which cases studies had to be split into separate datasets ( which were treated as batches during integration ). reasons for possible splitting were: ( 1 ) different 10x versions used within a study ( for example, 10xv2 versus 10xv3 ); or ( 2 ) the processing of samples at different institutes within a study. to determine whether these covariates caused batch effects within a study, we performed principal component regression. to this end, we preprocessed single studies separately ( total count normalization to median total counts across cells and subsequent pca with 50 principal components ). for each study, we then calculated the fraction of the variance in the first 50 principal components that could be explained ( pc expl ) by the covariate of interest ( that is, 10x version or processing institute ): [ 12pt ] $ $ } _ } } = _ ^ } ( } } ) } } _ ^ } ( } _ i } ) } } $ $ pc expl = ∑ i = 1 50 var cov ∑ i = 1 50 var pc i where var ( pc i | cov ) is the variance in scores for the i th principal component across cells that can be explained by the covariate under consideration, based on a linear regression. then, 10x version or processing institute assignments were randomly shuffled between samples and pc expl was calculated for the randomized covariate. this was repeated over ten random shufflings and the mean and standard deviation of pc expl were then calculated for the covariate. if the nonrandomized pc expl was more than 1. 5 standard deviations above the randomized pc expl, we considered the covariate a source of batch effect and split the study into separate datasets. thus, both jain _ misharin _ 2021 and ref. were split into 10xv1 and 10xv2; ref. was split into 10xv2 and 10xv3; and ref. and its pooled unpublished data were not split based on 10x version nor on processing location. for integration of the datasets into the hlca core, coarse cell type labels were used as described for integration benchmarking ( at1, at2, arterial endothelial cell, b cell lineage, basal, bronchial vessel 1, bronchial vessel 2, capillary, multiciliated, dendritic, fibroblast lineage, krt5 krt17 epithelial, lymphatic endothelial cell, macrophages, mast cells, megakaryocytes, mesothelium, monocytes, neutrophils, natural killer / natural killer t cells, proliferating cells, rare, secretory, smooth muscle, squamous, submucosal secretory, t cell lineage, venous and unlabeled ), except cells with lacking annotations were set to unlabeled instead of being removed. scanvi was run on the raw counts of the 2, 000 most hvgs ( calculated as described above ), using datasets as the batch variable to enable the conservation of interindividual variation. the following parameter settings were used: number of layers: 2; number of latent dimensions: 30; encode covariates: true; deeply inject covariates: false; use layer norm: both; use batch norm: none; gene likelihood: nb; n epochs unsupervised: 500; n epochs semi - supervised: 200; and frequency: 1. for the unsupervised training, the following early - stopping parameters were used: early stopping metric: elbo; save best state metric: elbo; patience: 10; threshold: 0; reduce lr on plateau: true; lr patience: 8; and lr _ factor: 0. 1. for the semisupervised training, the following early - stopping parameter settings were used: early stopping metric: accuracy; save best state metric: accuracy; on: full dataset; patience: 10; threshold: 0. 001; reduce lr on plateau: true; lr _ patience: 8; and lr _ factor: 0. 1. the integrated latent embedding generated by scanvi was used for downstream analysis ( clustering and visualization ). for gene - level analyses ( differential expression and covariate effect modeling ), uncorrected counts were used. to cluster the cells in the hlca core, a nearest neighbor graph was calculated based on the 30 latent dimensions that were obtained from the scanvi output, with the number of neighbors set to k = 30. this choice of k, while improving clustering robustness, could impair the detection of very rare cell types. coarse leiden clustering was done on the graph with a resolution of r = 0. 01. for each of the resulting level 1 clusters, a new neighbor graph was calculated using scanvis 30 latent dimensions, with the number of neighbors again set to k = 30. based on the new neighbor graph, each cluster was clustered into smaller level 2 clusters with leiden clustering at a resolution of r = 0. 2. the same was done for levels 3 and 4 and ( where needed ) 5, with k set to 15, 10 and 10, respectively, and the resolution set to 0. 2. clusters were named based on their parent clusters and sister clusters ( for example, cluster 1. 2 is the third biggest subcluster ( starting at 0 ) of cluster 1 ). for visualization, a 2d umap of the atlas was generated based on the 30 nearest neighbors graph. to quantify cluster cell type label disagreement for a specific level of annotation, the label shannon entropy was calculated on the cell type label distribution per cluster as [ 12pt ] $ $ - _ ^ k p ( ) } [ ) } ], $ $ − ∑ i = 1 k p x i log p x i, where x 1 … x k are the set of labels at that annotation level and p ( x i ) is the fraction of cells in the cluster that was labeled as x i. cells without a label at the level under consideration were not included in the entropy calculation. if < 20 % of cells were labeled at the level under consideration, the entropy was set to not available for the figures. the entropy of donors per cluster ( that is, diversity of donors in a cluster ) was calculated in the same way. to set a threshold for high label entropy, we calculated the label entropy of a hypothetical cluster with 75 % of cells given one label and 25 % of cells given another label, as a cluster with < 75 % of cells with the same label suggests substantial disagreement in terms of cell type labeling. clusters with a label entropy above that level ( 0. 56 ) were considered to have high label entropy.