dataesr/scientific-paragraphs-categorization huggingface dataset:
While the HLCA is the most diverse lung study to date, it is still a reflection of the limited diversity in singlecell studies of the respiratory system. Thus, the HLCA core consists of data from 65% white individuals, highlighting the urgent need for diversification of the population sampled in lung studies. Similarly, the current HLCA core does not include BALF, autopsy or sputum samples. Expanding the atlas with these sample types will further diversify captured cell identities, and could improve the quality of the HLCA core as a reference for new datasets acquired by these sampling methods. Future extensions of the HLCA could further increase the power to detect the effects of demographic variables, such as age, BMI, and smoking status. Indeed, in the current HLCA the effect of smoking on cell types was relatively modest, also compared to other demographic covariates.
As metadata on pack years and smoking frequency was not available for most subjects, a coarser encoding of smoking status (current, ever, never) was necessary, increasing the within-group variance and thus decreasing statistical power. A larger number of subjects in the HLCA would counteract this limitation and could provide the power needed to detect more nuanced effects. Similarly, in the current HLCA core certain subsets of T cells (regulatory T cells, gamma-delta T cells) could not be identified as separate clusters. Mapping of additional, well-annotated immune cell data could provide the power to confidently detect these cell identities. Indeed, we showed that mapping datasets to the HLCA enabled refinement of HLCA annotations. The HLCA core is publicly available as a data portal and a model repository to explore, download, and use as a reference for new datasets.
As the atlassing community has multiple outlets for newly generated data, we made the atlas available in Sfaira 70 , Zenodo 71 , Azimuth 14 , CellTypist  72  , and FASTGenomics  73. Mapping new datasets to the HLCA core using scArches can be done via interactive portals, such as FASTGenomics and Azimuth, as well as locally using our Zenodo model  74  , which lends itself to integration into bioinformatics pipelines. As more lung datasets become available and the respiratory community progresses in charting the cellular landscape of the lung, for example via novel sampling techniques or multi-modal analysis such as combined spatial and singlecell analysis, knowledge of the lung and its cell types will also expand. This underlines the importance of viewing the HLCA as a 'live' resource that requires continuous updates.
Indeed, recent studies  1, 17, 22, 43  have proposed more, novel cell types that have not yet been included in the hierarchical, harmonized HLCA annotations as these do not yet have support from multiple studies or consensus within the field. Using recent advances in transfer learning, we can robustly and iteratively extend the current HLCA. We illustrated, by mapping a newly annotated dataset to the HLCA core, that incorporating novel data in the atlas enables the refinement of HLCA annotations. To consolidate this incorporation of new knowledge, the HLCA core will be continuously updated by adding new datasets as these become available. New builds of the HLCA will be established by fully reintegrating the original and new datasets into an updated HLCA core. In addition, versions of every build will include fast, transfer learning-based mappings of individual, newly published datasets to the latest available HLCA core build.
With the integration of new datasets and the open access of all HLCA infrastructure, the HLCA can serve as a community-and data-driven platform for open discussion on lung cell heterogeneity, thus providing the basis for reaching consensus on updated lung cell type annotations. Through this iterative process, the HLCA will function as a central reference for lung data, reflecting the latest knowledge of the lung and its cell identities in health and disease. In the future, multimodal data of the human respiratory system -including epigenomic, spatial, and imaging data -will become increasingly available. While first ideas for integrating such data have been proposed  75  , a systematic study of multimodal atlas building is not yet feasible due to limited data. Similarly, integrating singlecell atlases with large population cohort studies assaying human transcriptomic variation on the bulk level would greatly extend the represented diversity that is lacking in available single-cell studies.
Integration across bulk and single-cell data would however require methodological advances in bulk deconvolution and imputation. Overall, we expect future versions of the HLCA to also contain a multimodal and a spatial "track", which will enable deeper insights into cellular phenotypes within the organ. Taken together, the HLCA represents a major milestone for the respiratory community, providing a single-cell resource of unprecedented size and quality. It offers a model framework for building integrated, consensusbased, population-scale atlases for further organs within the Human Cell Atlas. The HLCA is publicly available, and combined with an open-access platform to map new datasets to the atlas, this resource paves the way towards a better and more complete understanding of both health and disease in the human lung.
Re-processing of published data The raw sequencing data of 4 previously published studies was re-aligned to GRCh38, ensembl 84 for the HLCA (Krasnow_2020, Barbry_Leroy_2020, Seibold_2020, Meyer_2021). For the Krasnow_2020 and Barbry_Leroy_2020 studies, the CellRanger 3. 1. 0-based HLCA pipeline was used for re-alignment. Seibold_2020 and Meyer_2021 data were processed as previously described  22, 43  , but with reference genome and genome annotation adapted to the HLCA (GRCh38, ensembl84). All other datasets in the HLCA core were originally already aligned to GRCh38, ensembl84, except data from the Lafyatis_Rojas_2019 study (GRCh38, ensembl93) (Supplementary Data Table  1 ).
Single-cell sequencing and preprocessing of unpublished data Unpublished data was generated as follows: Barbry_unpubl: Barbry_unpubl: Nasal and tracheobronchial samples were collected from IPF patients after obtention of their informed consent, following a protocol approved by CHU Nice, registered at clinicaltrials under reference NCT04529993. The libraries were prepared as described in Deprez et al  13  , and yielded an average of 61,000±11,000 cells per sample, with a viability above 95%. The single cell suspension was used to generate single cell libraries following the V3. 1 protocol for 3' chemistry from 10x Genomics (CG000204). Sequencing was performed on a NextSeq 500/550 sequencer (Illumina). Raw sequencing data were processed using the cellranger-6. 0. 0 pipeline, with reference genome and annotation GRCh38 and ensembl98.
Raw sequencing data was processed using the 10X Cell Ranger v3 pipeline and GRCh38 (hg38) reference genome. For downstream analysis, mitochondrial transcripts and doublets identified by DoubletDetection 79 v2. 4. 0 were removed. Samples were then combined and cell barcodes were filtered based on genes detected (>200, <7500) and gene UMI ratio (gene. vs. molecule. cell. filter) using Pagoda2 (github. com/hms-dbmi/pagoda2). Also using Pagoda2 for clustering, counts were normalized to total counts per nucleus. For batch correction, gene expression was scaled to dataset average expression. After variance normalization, all significantly variant genes 4519 were used for PCA. Clustering was done at different k values (50, 100, 200) using the top 50 principal components and the infomap community detection algorithm.
The single-cell suspensions used for 10X Genomics scRNA-seq analysis were cleared of red blood cells by using a red blood cell lysis bufer (eBioscience) followed by live cell counting and loading of 10,000 cells per lane. We used 10XGenomics Chromium Single-Cell 3' Reagent Kits v2 and v3 according to the manufacturers' instructions. Raw sequencing data were processed using the CellRanger 3. 1. 0-based HLCA pipeline, with reference genome and annotation GRCh38 and ensembl84. Ambient RNA correction was performed using FastCAR (https://github. com/LungCellAtlas/FastCAR), using an empty library cutoff of 100 UMI and a maximum allowed contamination chance of 0. 05, ignoring the mitochondrial RNA. Data were merged and processed using Seurat  28  , filtering to libraries with >500 UMIs and >200 genes and to the libraries containing the lowest 95% of mitochondrial RNA per sample and <25% mitochondrial RNA, normalized using sctransform  80  while regressing out % mitochondrial RNA.
In general, 15 principal components were used for the clustering, at a resolution of 0. 5 to facilitate manual annotation of the dataset. Clusters in the final object that were driven by single donors were removed. Raw counts and cell annotations were used as input for the HLCA. Jain_Misharin_2021: Nasal epithelial samples were collected from healthy volunteers who provided informed consent at Northwestern Medicine, Chicago, USA. Protocol was approved by Northwestern University IRB (STU00214826). Briefly, subjects were seated and asked to extend their neck. A nasal curette (Rhino-pro, VWR) was inserted into either nare and gently slid from posterior-to-anterior about one centimeter along the lateral inferior turbinate. Five curettes were obtained per participant. The curette tip was then cut and placed in 2 ml of hypothermosol and stored at 4C until processing.
Data and metadata collection To accommodate data protection legislation, single cell RNA sequencing datasets of healthy lung tissue were shared by dataset generators (Supplementary Data Table  2 ) as raw count matrices, thereby obviating the need to share genetic information. Unpublished raw count matrices and cell annotations were generated as described above. For the Barbry_2020 study, count matrices provided had ambient RNA removed. Count matrices were generated using varying software (Supplementary Data Table  2 ). All count matrices of the HLCA core except one study (see above, "re-processing of published data") were based on alignment to reference genome GRCh38, annotation ensembl 84. For all datasets from the HLCA core, a pre-formatted sample metadata form was filled out by the dataset providers for every sample, containing metadata such as subject ID of the donor from which the sample came, the donor's age, the type of sample, the sequencing platform, etc.
Cell type reference creation and metadata harmonization for the HLCA core To harmonize cell type labels from different datasets, a common reference was created to which original cell type labels were mapped (Supplementary Data Table  4 ). To accommodate labels at different levels of detail, the cell type reference was made hierarchical, with level 1 containing the coarsest possible labels (immune, epithelial, etc. ), and level 5 containing the finest possible labels (e. g. naive CD4 T cells). Levels were created in a data-driven fashion, recursively breaking up coarser level labels into finer ones where finer labels were available. Anatomical location of the sample was encoded into a continuous score, with 0 representing the most proximal samples (nose, inferior turbinate) and 1 representing the most distal possible sample (parenchyma, alveolar sacs) (Supplementary Data Table  3 ).
Preprocessing of data for integration benchmarking For computational efficiency, benchmarking was performed on a subset of the total atlas, including data from 10 studies split into 13 datasets (Lafyatis_Rojas_2020 was split into 10Xv1 and 10Xv2 data, Seibold_2020 was split into 10Xv2 and 10Xv3 data, and Banovich_Kropski_2020 was split into two based on processing site). The data came from 72 subjects, 124 samples and 372,111 cells. Preprocessing of the benchmarking data included the filtering of cells (minimum number of total UMI counts: 500) and genes (minimum number of cells expressing the gene: 5). For integration benchmarking, the scIB benchmarking framework was used  16. All benchmarked methods except scGen (i. e. BBKNN, ComBat, Conos, fastMNN, Harmony, Scanorama, scANVI, scVI and Seurat RPCA) were run at least twice: on the 2000 most highly variable genes, and on the 6000 most highly variable genes.
Of those methods, all methods that did not require raw counts as input were run twice on each gene set: once with gene counts scaled to have mean 0 and standard deviation 1, one with unscaled gene counts. scVI and scANVI, which require raw counts as input, were not run on scaled gene counts. scGen was only tested on 2000 unscaled highly variable genes. For highly variable gene selection, first, highly variable genes were calculated per dataset, using cellranger-based highly variable gene selection 87 (default parameter settings: min_disp=0. 5, min_mean=0. 0125, max_mean=3, span=0. 3, n_bins=20). 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 etc. until a predetermined number of genes was selected (2000 or 6000 genes).
Splitting of studies into datasets For integration of the data into the HLCA core, we first determined in 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 (e. g. 10Xv2 versus 10Xv3), or 2) processing of samples at different institutes within a study. To determine if these covariates caused batch effects within a study, we performed principal component regression 88. To that end, we preprocessed single studies separately (total counts normalization to median total counts across cells, and subsequent principal component analysis (PCA) with 50 principal components (PCs). For each study, we then calculated the fraction of the variance in the first 50 PCs that could be explained ("pcexpl") by the covariate of interest (i.
e. 10X version or processing institute): 𝑝𝑐𝑒𝑥𝑝𝑙 = ∑ !" #$% "#$(&' # |')") ∑ !" #$% "#$(&' # ) where var(pcI|cov) is the variance in scores for the ith PC 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 pcexpl was calculated for the randomized covariate. This was repeated over 10 random shufflings, and the mean and standard deviation of the pcexpl was then calculated for the covariate. If the non-randomized pcexpl was more than 1. 5 standard deviations above the randomized pcexpl, we considered the covariate a source of batch effect and split the study into separate datasets. Thus both Jain_Misharin_2021 and Lafyatis_Rojas_2019 were split into 10Xv1 and 10Xv2, Seibold_2020 was split into 10Xv2 and 10Xv3, while Banovich_Kropski_2020 was not split based on 10X version nor based on processing location.
Dataset integration For integration of the datasets into the HLCA core, coarse cell type labels were used as described for integration benchmarking, except cells with lacking annotations were set to "unlabeled" instead of being removed. scANVI was run on the raw counts of the 2000 most highly variable genes (calculated as described above), using datasets as the batch variable. 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, 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, lr_factor: 0.
1. For the semi-supervised 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, lr_factor: 0. 1. The latent embedding generated by scANVI was used for downstream analysis (clustering and visualization). For gene-level analyses (differential expression, covariate effect modeling) un-corrected counts were used.
Data embedding and clustering 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 resolution 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 resolution r=0. 2. The same was done for level 3, 4 and where needed 5, with k set to 15, 10, and 10 respectively, and resolution set to 0.
Variance between individuals explained by technical and biological covariates To quantify to what extent different technical and biological covariates correlated with inter-individual variation in the atlas, we calculated the "variance explained" by each covariate for each cell type. We first split the data in the HLCA core by cell type annotation, merging substates of a single cell type into one (Supplementary Data Table  10 ). For every cell type, we excluded samples that had fewer than 10 cells of the sample. We then summarized covariate values per sample for every cell type as follows: mean across cells from sample for scANVI latent components (integration results), UMI count per cell, and fraction of mitochondrial UMIs; for all other covariates (i. e. dataset, 3' or 5', BMI, cell ranger version (1, 2, 3 or 4), cell viability %, ethnicity, sampling method (e.
We furthermore filtered out datasets with fewer than 2 subjects, and refrained from modeling categories in covariates that had fewer than three subjects in their category, for that cell type. To determine if covariance between covariates was low enough for modeling, we calculated the variance inflation factor (VIF) between covariates at subject level. The VIF quantifies multi-collinearity 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 covariance was too high and excluded that cell type from the modeling. As many cell types lacked sufficient representation of other ethnicities than ''white", while including ethnicity in the analysis simultaneously reduced the samples that could be included in the analysis to only those with ethnicity annotations, we excluded ethnicity from the modeling.
Gene counts were summed across cells for every sample, within cell type. Sample-wise sums (i. e. pseudo-bulks) were normalized using edgeR's calcNormFactors function, using default parameter settings. We then used voom  91  , 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 or "Dream") that allows for mixed effect modeling, and modeled gene expression as: 𝑙𝑜𝑔(𝑛𝑜𝑟𝑚𝑐𝑜𝑢𝑛𝑡) ∼ 1 + 𝑎𝑔𝑒 + 𝑠𝑒𝑥 + 𝐵𝑀𝐼 + 𝑠𝑚𝑜𝑘𝑖𝑛𝑔 + 𝑛𝑜𝑠𝑒 + 𝑐𝑐𝑓𝑠𝑐𝑜𝑟𝑒 + (1 | 𝑑𝑎𝑡𝑎𝑠𝑒𝑡) where dataset is treated as a random effect, and all other effects are modeled as fixed effects. Resulting pvalues were corrected for multiple testing within every covariate using the Benjamini-Hochberg procedure.
Mapping of GWAS results to the lung cell atlas cell types To stratify GWAS results from several lung diseases by lung cell type, we applied stratified sc-LDSC, a method that can link GWAS results to cell types based on proximity of disease-associated variants to genes differentially expressed in the cell type  51. GWAS summary statistics of COPD  49  (GWAS catalog ID: GCST007692, dbGaP accession number: phs000179. v6. p2, n_cases=35,735, n_controls=222,076), and of lung cancer  48  (GWAS catalog ID: GCST004748, dbGaP accession number: phs001273. v3. p2, n_cases=29,266, n_controls=56,450) were made available on dbGap upon request. Summary statistics of lung function  50  (GWAS catalog ID: GCST007429, n=321,047 individuals), of asthma 47 (GWAS catalog ID: GCST010043, n_cases=88,486 n_controls=447,859), and of depression  98  (used as negative control, GWAS catalog ID: GCST005902, n_cases=113,769, n_controls=208,811) were publicly available.
A differential gene expression test was done for every grouped cell type (Supplementary Data Table  10 ) in the HLCA using a t-test, testing against the rest of the atlas. The top 1000 most significant genes with positive fold-changes were stored as genes characterizing that cell type ("cell type genes"), and used as input for sc-LDSC. Gene coordinates of cell type genes were obtained based on the GRCh37. 13 genome annotation. For single nucleotide-polymorphism (SNP) data (names, locations, and linkage-related information), the 1000 genomes European reference (GRCh37) was used, as previously described  51. Only SNPs from the HapMap3 project were included in the analysis. For identification of SNPs in the vicinity of cell type genes, we used a window size of 100,000 base pairs. Genes from X and Y chromosomes, as well as HLA genes, were excluded because of their unusual genetic architecture and linkage patterns.
Extension of the HLCA core by mapping of scRNA-seq and snRNA-seq lung data To map unseen scRNA-seq and snRNA-seq data to the HLCA, we used scArches, our transfer-learning based method that enables mapping of new data to an existing reference atlas  42. 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 re-analysis and de novo clustering of the joint data. The data to map was subsetted to the same 2000 highly variable genes ("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 input for scArches. Healthy lung data (Meyer_2021) was split into two datasets: 3' and 5'-based. Lung cancer data (Thienpont_2018) was 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 adapter weights that allowed for mapping of new data into the existing HLCA embedding, using the following parameters: freeze-dropout: true, surgery_epochs: 500, train base model: false, metrics to monitor: accuracy and elbo, weight-decay: 0, frequency: 1. The following earlystopping 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: 8, lr_factor: 0. 1.
We clustered at resolutions 0. 1, 0. 2, 0. 5, 1, 2, 3, 5, 10, 15, 20, 25, 30, 50, 80, and 100. Minimum recall (% of cells with the spatial cell type annotation captured in cluster) was set to 20%, and minimum precision (% of cells from Meyer_2021 study in the cluster that had the spatial cell type annotation) was set to 25%. The cluster with the highest recall was selected for every spatial cell type (unless this cluster decreased precision by more than 50% compared to the cluster with the second highest recall). If precision of the next best cluster was doubled compared to the cluster with highest recall, and recall did not decrease with more than 20%, this cluster was selected. 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  42.