dataesr/scientific-paragraphs-categorization huggingface dataset:
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  28  ). 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 ).
Splitting of studies into datasets 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  101. 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): PC expl = ∑ 50 i=1 var (cov) ∑ 50 i=1 var (PC i ) where var(PC i |cov) is the variance in scores for the ith 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. 22 were split into 10xv1 and 10xv2; ref. 25 was split into 10xv2 and 10xv3; and ref. 21 and its pooled unpublished data were not split based on 10x version nor on processing location.
Integration of HLCA core datasets with scANVI 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.
Thresholds for high label/donor entropy and doublet clusters 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. Six small clusters with high label entropy even at the coarsest level of annotation highlighted doublet populations (identified via simultaneous expression of lineage-specific marker genes; for example, expression of both epithelial (AT2) and stromal (smooth muscle) marker genes) not labeled as such in the original datasets. These clusters were removed from the HLCA core, bringing the total number of clusters to 94.
The number of cells per cell type is shown for all levels in Supplementary Table  5. Mislabeling of original cells was identified by comparing final annotations with original harmonized labels and checking whether these were contradictory (and not only done at different levels of detail). Out of 61 final cell types, 18 included mostly mislabeled cells, 12 of which were previously known cell types. Despite consisting of mostly mislabeled cells from the original datasets, individual experts agreed on the annotation of these cell types: for all previously known cell types with a high proportion of mislabeled cells, the majority of experts agreed on the final annotation for the atlas, or only differed in the granularity of annotation.
Next, we performed principal component regression on every covariate, as described previously (see the section 'Splitting of studies into datasets'), but now using scANVI latent component scores instead of principal component scores for the regression, yielding a fraction of latent component variance explained per covariate. Samples that did not have a value for a given covariate (for example, where the BMI was not recorded for the donor) were excluded from the regression. Categorical covariates were converted to dummy variables. Cell typecovariate pairs for which only one value was observed for the covariate were excluded from the analysis. Quantification of the correlation or dependence between variables within a cell type and identification of the minimum number of samples needed to control for spurious correlation are described below.
Cell type filtering covariate encoding for gene-level modeling 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.
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  106  , 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: 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  107. The analysis was performed in R using the cameraPR function in the limma package  108  , with the differential expression test statistic. Gene Ontology biological process terms  109, 110  were tested separately for each comparison. These sets were obtained from MSigDB (version 7. 1)  111  , as provided by the Walter and Eliza Hall Institute (https://bioinf. wehi. edu. au/MSigDB/index. html).
Generating cell type signature matrices for deconvolution To enable deconvolution of bulk expression data on the basis of the HLCA, HLCA cell type signature matrices were generated. One generic-purpose signature matrix was created per sublocation of the respiratory system (that is, one parenchyma, one airway and one nose tissue matrix; Supplementary Table  10 ). Additionally, a script to generate custom reference sets from the HLCA data is provided together with the HLCA code on GitHub (https://github. com/LungCellAtlas/ HLCA) to tailor the deconvolution signature matrix to any specific research question. Cell types were included in the bulk deconvolution signature matrix on the basis of cell proportions (constituting >2% of cells within samples of the corresponding tissue in the HLCA core).
Deconvolution of bulk expression data using the HLCA core The parenchymal signature matrix was used to deconvolve RNA expression data of samples from the Lung Tissue Database 52 (GEO accession number GSE23546) using only lung tissue samples from patients with COPD GOLD stages 3 and 4 (n = 27 and 56, respectively) and matched controls (n = 281). The Lung Tissue Database dataset was run on the Rosetta/Merck Human RSTA Custom Affymetrix 2. 0 microarray platform (HuRSTA-2a520709; GPL10379). As this platform has multiple probe sets for each gene, we focused on the probe sets that were derived from curated RefSeq records (with NM_ accession prefixes) when present to maximize the accuracy of the deconvolution. Where genes did not have probe sets based on curated RefSeq records or had multiple probe sets mapping to curated RefSeq records, the probe set with the highest average microarray intensity across samples was selected.
Quantile normalization of the data and subsequent deconvolution were performed using CIBERSORTx. A Wilcoxon rank-sum test between control and GOLD stage 3/4 samples was performed to identify statistically significant compositional changes in advanced-stage COPD compared with control tissue. GOLD 3/4 and control samples were matched for age and smoking history. Cell types with >60% of samples predicted to have a proportion of zero were excluded from the Wilcoxon test, as the high number of tied ranks (zeros in both groups) would result in inflated P values. P values were multiple testing corrected using the Benjamini-Hochberg procedure. The same procedure was followed for a dataset of nasal brush bulk RNA-seq samples from asthmatic donors pre-and postinhalation of corticosteroids (n = 54 and 26, respectively)  50  and a dataset of airway biopsy bulk RNA-seq samples from asthmatic donors and controls (n = 95 and 38, respectively)  51.
Extension of the HLCA core by mapping of new data To map unseen scRNA-seq and single-nucleus RNA-seq data to the HLCA, we used scArches, our transfer learning-based method that enables mapping of new data to an existing reference atlas  71. scArches trains an adaptor added to a reference embedding model, thereby enabling it to generate a common embedding of the new data and the reference, allowing reanalysis and de novo clustering of the joint data. The data to map were subsetted to the same 2,000 HVGs that were used for HLCA integration and embedding, and HVGs that were absent in the new data were set to 0 counts for all cells. Raw counts were used as Resource https://doi. org/10. 1038/s41591-023-02327-2 input for scArches, except for the ref. 40 dataset, for which ambient RNA removal was run previously on the raw counts.
Healthy lung data  40  were split into two datasets: 3′ and 5′ based. Lung cancer data  41  were also split into two datasets: 10xv1 and 10xv2. The model that was learned previously for HLCA integration using scANVI was used as the basis for the scArches mapping. scArches was then run to train adaptor weights that allowed for mapping of new data into the existing HLCA embedding, using the following parameter settings: freeze-dropout: true; surgery_epochs: 500; train base model: false; metrics to monitor: accuracy and elbo; weight-decay: 0; and frequency: 1. The following early-stopping criteria were used: early stopping metric: elbo; save best state metric: elbo; on: full dataset; patience: 10; threshold: 0. 001; reduce lr on plateau: True; lr patience: 8l and lr_factor: 0. 1.
Gene name harmonization To enable cross-dataset gene-level analysis, harmonization of gene names from different datasets (using different reference genome builds and genome annotations; Supplementary Table  1 ) was necessary. Both annotation sources (for example, Ensembl or RefSeq) and annotation versions (for example, Ensembl release 84 or Ensembl release 91) contribute to the variation between different gene naming schemes. Therefore, both annotation sources and versions, including outdated ones, need to taken account to enable the mapping of all gene names to a single naming scheme. For the harmonization of gene names, we aimed to map all original gene names to the target scheme HUGO Gene Nomenclature Committee gene name, corresponding to the Ensembl release 107 publication. To find the most likely match between each original gene name and a target gene name, we first downloaded Ensembl releases 79 to 107, which included for each release: (1) all Ensembl gene IDs from the downloaded release (for example, ENSG00000081237.
21); (2) corresponding Ensembl transcript and protein IDs (for example, ENST00000442510. 8 and ENSP00000411355. 3); (3) matching Ensembl IDs from the previous release; (4) matching gene IDs from other genome annotation sources (for example, RefSeq); and (5) matching gene, transcript and protein identifiers from various external resources, such as UniProt, the HUGO Gene Nomenclature Committee and the Consensus Coding Sequence Project. We then constructed a graph, with each Ensembl ID, other genome annotation ID and external resource identifier represented by a single node per release. Nodes were then connected based on the matching (points 2-5) provided by Ensembl, weighing edges based on Ensembl similarity scores where available. For each original gene name from the HLCA datasets, the path with the lowest mean edge weight from that gene name to a gene name from the target names (Ensembl release 107) was selected to find the most likely matching gene name from the target (Supplementary Table  17 ).
Identification of genes associated with common batch effects To identify the genes most commonly exhibiting batch-specific expression, the HLCA was split by cell type and a differential expression analysis was performed (based on a Wilcoxon rank-sum test) in each cell type, comparing cells from one dataset (batch) with those from all other datasets and repeating this for all datasets. Datasets that had fewer than ten cells of the cell type or fewer than three samples with cells of the cell type were excluded from the test. For each test, genes were filtered such that only genes that were significantly upregulated were retained. Next, the fraction of included datasets in which a gene was significantly upregulated in the cell type (affected dataset fraction) was calculated for all genes. To find genes that were often batch effect associated across many cell types, the mean of the affected dataset fractions was calculated across cell types for each gene.
Cell type label transfer from the HLCA core to new datasets To perform label transfer from the HLCA core to the mapped datasets from the extended HLCA, we used the scArches k nearest neighbor-based label transfer algorithm  71. Briefly, a k nearest neighbor graph was generated from the joint embedding of the HLCA core and the new, mapped dataset, setting the number of neighbors to k = 50. Based on the abundance and proximity in a cell's neighborhood of reference cells of different types, the most likely cell type label for that cell was selected. Furthermore, a matching uncertainty score was calculated based on the consistency of reference annotations among the k nearest neighbors of the cell of interest u c,y,N c = 1 -p (Y = y|X = c, N c ) where u c,y,Nc is the uncertainty score for a query cell c with transferred label y; N c is its set of k nearest neighbors; and p(Y = y|X = c, N c ) is the weighted (by edge weights) proportion of N c with label y, as previously described  113.
For example, γδ T cells, ILCs, megakaryocytes, natural killer T cells and regulatory T cells were not annotated in the HLCA core, as these cell types could not be distinguished with confidence in the integrated object and were often lacking in the constituent datasets. Indeed, cell types from the T cell/ILC/natural killer lineages are known to be particularly difficult to annotate using transcriptomic data only  16. Therefore, cells annotated with these labels in the projected dataset were largely incorrectly annotated as CD4 + T cells, CD8 + T cells and natural killer cells through label transfer (Fig. 5b  and Extended Data Fig. 6e ).
Calibration of uncertainty cutoff for classifying as unknown For the extended atlas, we calibrated the uncertainty score cutoff by determining which uncertainty levels indicate possible failure of label transfer. To determine the uncertainty score at which technical variability from residual batch effects impairs correct label transfer, we evaluated how label transfer performed at the level of datasets, as these predominantly differ in experimental design. To determine an uncertainty threshold indicative of possible failure of label transfer, we harmonized original labels for 12 projected datasets  54, 58, 59, 64, 66  (one unpublished: Duong_lungMAP_unpubl) and assessed the correspondence between original labels with the transferred annotations. Only cells with level 3 or 4 original annotations were considered, as these levels represent informative annotations while not representing the finest detail.
Level 5 annotations will often display high uncertainty levels due to high annotation granularity rather than remaining batch effects. To assess the optimal uncertainty cutoff for labeling a new cell as unknown, we used these results to generate a receiver operating characteristic curve. We chose a cutoff around the elbow point, keeping the false positive rate below 0. 5 (uncertainty cutoff = 0. 2; true positive rate = 0. 879; false positive rate = 0. 495) to best distinguish correct Resource https://doi. org/10. 1038/s41591-023-02327-2 from incorrect label transfers (Supplementary Fig. 10a ). False positives were either due to incorrect label transfer or incorrect annotations in the original datasets. Cells with an uncertainty higher than 0. 2 were set to unknown.
Cross-dataset analysis of IPF-associated cell states To uncover the cell identities affected in IPF, label transfer uncertainty was analyzed for three mapped datasets from the extended HLCA  58, 62, 64  that included both IPF and healthy samples. Cell types of interest were determined based on the largest increase in mean label transfer uncertainty in IPF compared with healthy samples, while checking for consistency in increments across the three datasets. This highlighted alveolar fibroblasts as the main cell type of interest. To find IPF-specific alveolar fibroblast states, all alveolar fibroblasts from the abovementioned datasets and two more IPF datasets  21, 24  (for which no healthy data were mapped, as these were already in the core) were clustered together with the alveolar fibroblasts from the HLCA core.
For clustering, a k nearest neighbor graph was calculated on the joint scArches-derived 30-dimensional embedding space setting k = 30, after which the cells were clustered using the Leiden algorithm with a resolution of 0. 3. The resolution was chosen such that datasets were not isolated in single clusters, thus avoiding clustering driven solely by dataset-specific batch effects. One cluster (cluster 5) was small (n = 460 cells) and displayed low donor entropy (0. 17), indicating that it almost exclusively came from a single donor (96% of cells from HLCA core donor 390C). It was therefore excluded from further analysis. To perform differential gene expression analysis, gene counts were normalized to a total of 7,666 counts (the median number of counts across the HLCA) and then log transformed with a pseudocount of 1.
Multidisease analysis To investigate whether the HLCA can be used to identify disease-associated cell states shared across multiple diseases, MDMs from the HLCA core, together with all cells from the mapped datasets labeled as MDMs based on label transfer, were jointly analyzed. Datasets and diseases with fewer than 50 MDMs were excluded from the analysis. The cells were subsequently clustered as described above for the cross-dataset IPF analysis. Finally, a Wilcoxon rank-sum test was used on the normalized data to detect differentially expressed genes per cluster (number of cells per cluster: n = 64,915 (cluster 0), 47,539 (cluster 1), 32,027 (cluster 2), 31,097 (cluster 3), 25,267 (cluster 4), 1,998 (cluster 5) and 307 (cluster 6)). The results were filtered as described above (Supplementary Table  15 ).
Version information The following tools and versions were used: R (version 4. Extended Data Fig. 7  | Mapping of unseen lung cancer data to the HLCA. a, UMAP of the jointly embedded HLCA (dark blue, plotted on top) and lung cancer data (gray). b, Same as but now plotting cells from the HLCA core in gray. Cells from the mapped data are plotted on top, and colored by the cancer type of the patient. c, Same as a, but now coloring cells from the HLCA core by their final final annotation, and coloring cells from the mapped cancer data in black. Cells from each of the compartments are outlined to ease visual identification of cell types by colors. d, Uncertainty of label transfer, shown for all cells from the mapped data. Regions dominated by high-uncertainty cells are labeled by the original cell type label.
Cells from the HLCA core are colored in gray. e, Uncertainty of label transfer (ranging from 0 to 1) for the mapped cells, subdivided by original cell type label. Number of cells per label is shown between brackets. Boxes of cell type labels not present in the HLCA core are colored red. Cell types are ordered by mean uncertainty. Boxes show median and interquartile range of uncertainty. Cells with uncertainties more than 1. 5 times the interquartile range away from the high and low quartile are considered outliers and plotted as points. Whiskers extend to the furthest non-outlier point. f, Sankey plot of original labels of the mapped data versus predicted annotations based on label transfer. Cells with uncertainty >0. 3 are labeled 'unknown'. g, Percentage of cells from newly mapped healthy lung dataset that are either annotated correctly or incorrectly by label transfer annotation (matched at the level of the original labels), or annotated as unknown, subdivided by original cell type label.