dataesr/scientific-paragraphs-categorization huggingface dataset:
We show that CellTypist is able to integrate diverse datasets to alleviate their batch effects and maintain the cell type relationships across datasets. Having a consistent database and search architecture to quantitatively compare cell types is critical to the single-cell field  44, 45. The seamless integration of the two functionalities (harmonization and integration) in CellTypist allows us to semi-automate the atlas building effort. We harmonized ~250 high-quality cell types from more than 3. 6 million cells for 12 human tissues, providing, to our knowledge, the first multi-organ resource that allows information of cell type harmonization, integration, and classification to be simultaneously inspected and queried. We also anticipate that this highly structured CellTypist architecture, together with the multi-organ landmark resource, will aid in the assembly of diverse single-cell resources, facilitate the propagation of cell knowledge across communities, and improve the interpretability of organto-organ and health-to-disease cell relationships.
Data and code availability All datasets used in this study are publicly available. The five immune datasets were downloaded from the tissue stability cell atlas (https://www. tissuestabilitycellatlas. org), crosstissue immune cell atlas (https://www. tissueimmunecellatlas. org), gut cell atlas (https://www. gutcellatlas. org), HCA data portal (https://data. humancellatlas. org/explore/projects/cc95ff89-2e68-4a08-a234-480eca21ce79), and developmental cell atlas (https://developmental. cellatlas. io/fetal-liver). The four lung datasets were downloaded from Gene Expression Omnibus: GSE136831, GSE128033, GSE135893, and GSE122960. The 38 datasets used for building the multi-organ atlases are also publicly available, with their detailed information listed in Table  S1. CellTypist.
Induction of the PCT-based tree A direct calculation of the transcriptomic distance between two cells (or cell types) usually suffers from substantial batch effects when they originate from two different studies. To overcome this, CellTypist leverages the PCT algorithm to predict the transcriptional dissimilarities instead of calculating them. Since there is no ready-made PCT algorithm for use in the single-cell field, we adopt the commonly used tree structure, that is, parallel arrays from the Python scikit-learn package  46  (version 0. 24. 1), and implement the PCT framework adapted for single-cell transcriptomics. A distance matrix 𝑑 𝑐,𝑐𝑡 is first constructed for each dataset independently, where 𝑐 and 𝑐𝑡 represent the cells and cell types (defined as the cell centroids) in a given dataset, respectively.
The general assumption here is that cells within one dataset are less impacted by batch effects in comparison to those across datasets. The metric of correlation distance (1 -Pearson's 𝑟) or Euclidean distance is used to calculate this matrix based on the gene expression space or low-dimensional latent representations (e. g. , PCA), respectively. For each cell 𝑐 (each row), this matrix records its dissimilarity profile across cell types. We next define two key components of the clustering tree: distance measure and prototype  47. In CellTypist, the dissimilarity between two cells in the matrix 𝑑 is assessed using the Euclidean distance, and the prototype 𝑝(𝑛) is defined as the dissimilarity profile averaged across cells that are involved during the tree induction of a given node 𝑛. Accordingly, the intra-node variance 𝜎 2 (𝑛) is defined as the sum of squared distances between each cell 𝑐 and the prototype 𝑝(𝑛): where 𝑖 enumerates each cell contained in the node 𝑛, 𝑗 enumerates each cell type in the reference dataset of interest, and 𝑐(𝑖, 𝑗) is the dissimilarity value for cell 𝑖 and cell type 𝑗.
The proxy impurity improvement is then quantified by summing the impurities of the two daughters of the node 𝑛 and adding a negative sign. The clustering tree is split at each possible node by a "test" corresponding to an expression cutoff in a given gene that mostly separates cells within this node. By incorporating all above optimizations, a clustering tree is induced for each dataset, with a stringent requirement for node sizes (i. e. , the number of cells to be greater than 20 which is adjustable in CellTypist). However, the clustering tree is not exhaustively split to avoid overfitting. In the meantime, the resulting tree is not pruned by cost complexity pruning (CCP) as used by many tree-related algorithms, but by a more intuitive F-test as detailed below.
Impurity reduction test and tree pruning CellTypist uses an F-test to determine whether the clustering tree needs to be further induced to a lower hierarchy. Compared to CCP, this approach underscores the clustering tendency of similar cells which likely correspond to a homogeneous cell state but places less emphasis on the tree complexity (i. e. , number of leaves generated by possible splits). The level of significance for a possible split at the node 𝑛 in the reference dataset can be assessed by an Ftest: 𝐹 |𝑐𝑡|×(|𝑛|-1),|𝑐𝑡|×(|𝑛|-2) = 𝜎 2 (𝑛) |𝑐𝑡| × (|𝑛| -1) 𝜎 2 (𝑙) + 𝜎 2 (𝑟) |𝑐𝑡| × (|𝑛| -2) where 𝜎 2 (𝑛), 𝜎 2 (𝑙), and 𝜎 2 (𝑟) are the intra-cluster variances in the node, left daughter, and right daughter, respectively. Replacing all the three variances in this equation with the impurity values (𝐼(𝑛), 𝐼(𝑙), 𝐼(𝑟)) in CART allows us to infer the F-statistics in PCT as: 𝐹 |𝑐𝑡|×(|𝑛|-1),|𝑐𝑡|×(|𝑛|-2) = 𝐼(𝑛) × |𝑛| |𝑛| -1 𝐼(𝑙) × |𝑛𝑙| + 𝐼(𝑟) × |𝑛𝑟| |𝑛| -2 A p-value cutoff which defaults to 0.
Prediction of dissimilarity profiles across datasets After a PCT is built for each dataset separately, it will be used to predict the cell-by-cell-type distance matrix 𝑑 for cells in the remaining datasets. Each query cell will be, according to its gene expression profile, assigned to one of the leaves in the clustering tree and predicted to be the prototype. The resulting distance matrices are then concatenated across cells from all datasets to form an expanded matrix 𝐷 ̂𝑐,𝑐𝑡 where 𝑐 covers cells from all datasets and 𝑐𝑡 represents the cell types from the reference dataset of interest. Prediction accuracy is calculated as the coefficient of determination ( 𝑅 2 ) between the actual ( 𝐷 𝑐,𝑐𝑡 ) and predicted ( 𝐷 ̂𝑐,𝑐𝑡 ) dissimilarity profiles: 𝑅 2 = 1 - ∑ ∑ (𝐷 𝑖,𝑗 -𝐷 ̂𝑖,𝑗 ) 2 |𝑐𝑡| 𝑗=1 |𝑛| 𝑖=1 ∑ ∑ (𝐷 𝑖,𝑗 -𝐷 ̅ 𝑗 ) 2 |𝑐𝑡| 𝑗=1 |𝑛| 𝑖=1 𝐷 ̂𝑐,𝑐𝑡 is next turned into a rank matrix starting from 1 to the number of elements in this matrix, and is further maximum-normalized to lie between 0 and 1.
This is also the strategy adopted by MetaNeighbor  48. Finally, all 𝐷 ̂𝑐,𝑐𝑡 are concatenated across datasets, forming a global matrix 𝐺 𝑐,𝑐𝑡 of dissimilarities among all cells (𝑐) and all cell types (𝑐𝑡). Of note, because a cell type is usually defined at the cluster level and theoretically no cluster is pure, CellTypist can optionally filter out cells whose gene expression profiles do not correlate most with the cell type they belong to. This will speed up the run as only a subset of cells are used, but will render these filtered cells unharmonized in the downstream workflow (see the section below). Meta-analysis is performed by averaging cells (i. e. , rows) from the same cell type in the matrix 𝐺, leading to a cell-type-by-cell-type distance matrix among all cell types studied. This matrix is further transposed and averaged with the original matrix to form the final matrix that enables a consistent visualization given its symmetry.
Cell type harmonization Before undergoing the following stages, all cells are predicted into the most likely cell types in each dataset, on the basis of the global matrix 𝐺 mentioned above. Specifically, for each dataset, except for the cells which originate from this dataset, all other cells will be assigned one of the cell types that has the smallest distance with the query cell. We denote the result as the assignment matrix 𝐴 ∈ 𝑅 𝑐×𝑑 in CellTypist depicting how each cell 𝑐 is designated in each dataset 𝑑.
Stage I: initialization of the dataset core and definition of cell type relations Two datasets are first selected, named as D1 and D2. Using their cross-prediction information in 𝐴, CellTypist calculates two confusion matrices 𝑀12 and 𝑀21, representing how cell types in D1 are assigned across cell types of D2, and vice versa. CellTypist performs a pairwise alignment to form the dataset core for the two datasets with an iterative cell type removal procedure: 1. Initialize an empty harmonization table 𝐻𝑇 with three columns (D1, relation, D2) 2. Normalize each row of 𝑀12 and 𝑀21 to a sum of 1 From step 3 onwards, if 𝑀12 or 𝑀21 becomes an empty matrix in any step, the workflow will immediately end. A novel cell type C in D1 is appended as "C = NONE"; a novel cell type C in D2 is appended as "NONE = C"; an 1:1 alignment between C1 from D1 and C2 from D2 is appended as "C1 = C2"; an 1:N alignment between C1 from D1 and C2 and C3 from D2 is appended as "C1 ∋ C2" and "C1 ∋ C3"; a N:1 alignment between C2 and C3 from D1 and C1 from D2 is appended as "C2 ∈ C1" and "C3 ∈ C1"; an unharmonized cell type C from D1 is appended as "C = UNRESOLVED"; an unharmonized cell type C from D2 is appended as "UNRESOLVED = C".
1), and these matched cell types in D2 have a back-match proportion greater than minimum_unique_percent (the same parameter used in 3), this will be an 1:N match. 5) If based on 𝑀21, one cell type C in D2 has more than two cell types aligned in D1 with a match proportion greater than minimum_divide_percent, and these matched cell types in D1 have a back-match proportion greater than minimum_unique_percent, this will be a N:1 match. 6) Based on 𝑀12, if after all these categorizations, a cell type in D1 remains unharmonized (at this time, 𝑀12 should have no columns available), this will be an unharmonized cell type in D1. 7) Based on 𝑀21, if after all these categorizations, a cell type in D2 remains unharmonized, this will be an unharmonized cell type in D2. Of note, during multi-dataset (>=3) alignments, all three cutoffs used above are dynamically chosen to find the best alignments during iterations (see stage IV).
Stage III: expansion of the dataset core and harmonization table After the new D1 and associated annotations of cells in D1 are obtained, an additional dataset D2 (which is still named D2 because the earlier D1 and D2 have been combined and renamed to D1) is incorporated. The dataset core is expanded by aligning the D2 with D1 using a procedure below: 1. Form the 𝑀12 matrix (rows) as defined in the earlier 𝐻𝑇. This new 𝐻𝑇 will be expanded to have multiple columns with one dataset per column joined by a "relation" column, serving as input for step 5. When all datasets are iteratively incorporated, the final output will be a 𝐻𝑇 covering all cell types across all datasets.
Stage IV: tuning the optimal dataset order and cell reannotation Although all pairwise predictions between datasets are considered by the above stages (information from the assignment matrix 𝐴), an optimal order of datasets is still useful with respect to grouping of similar datasets and visualization. CellTypist uses a greedy strategy to determine the optimal order of datasets. CellTypist first calculates the similarity between any pair of datasets. Next, two datasets with the highest similarity are added in the first iteration of the harmonization pipeline, followed by addition of another dataset which has the highest total similarity with the datasets already added in the previous round. Here, the inter-dataset similarity is quantified by the numbers of alignments between two datasets (with adjustable weights of 2, 1, -1, and -2 for one-to-one, one-to-many, novel, and unharmonized matches, respectively) divided by the total number of cell types, denoting the degree of cell type share between the two given datasets.
During stage I to III, combinations of five minimum_unique_percent (0. 4, 0. 5, 0. 6, 0. 7, 0. 8) and three minimum_divide_percent (0. 1, 0. 15, 0. 2) are exhaustively tried to find the best combo, by choosing the parameters that result in the least alignments when a new dataset is added. This procedure is repeatedly conducted in each iteration. After the final 𝐻𝑇 is obtained, each cell will be reannotated to one of the cell types (rows) in 𝐻𝑇 based on the assignment matrix 𝐴.
Stage V: reorganization and visualization of 𝐻𝑇 and its cell type hierarchy An interconnected (i. e. , connections spanning multiple rows of 𝐻𝑇) group of cell types is further defined at the high hierarchy for 𝐻𝑇, with each group being independent of other groups. For each cell type group, if it contains novel ("NONE") or unharmonized ("UNRESOLVED") representations, these blank labels will be renamed in order to distinguish the many blank labels in the 𝐻𝑇. Specifically, a blank element in the first column (dataset) will be searched against cell types from the other columns across the given row until a non-blank one is found and renamed to it. A blank element in other columns (except the first column) will be searched in the opposite direction and renamed to the first non-blank element in the given row.
The newly formed 𝐻𝑇 forms the basis for the Sankey plot and the flat tree plot in CellTypist. Next, columns (datasets) of each cell type group will be reordered according to how many unique values each column has, including the renamed blank cell types. To disentangle the cross connections in 𝐻𝑇 for visualizations, rows of each group are further reordered to enable consistent alignments to lie at bottom of each group. This reordered 𝐻𝑇 is used to plot the cell type hierarchy by cross-linking unique cell types including the renamed blank cell types across datasets.
Annotation-aware data integration CellTypist employs the approach of restricted neighborhood search to integrate cells from different datasets (or batches), a concept similar to that used by BBKNN  16. After cell type harmonization (or other procedures that can result in consistent cell type labels across datasets), all cells are placed into a common naming schema. In CellTypist, for each cell belonging to a specific cell type, its nearest neighbors can only be found out of cells that belong to the same cell type. Since this strategy cannot produce a fully connected neighborhood graph and each cell type is strongly separated with other cell types, CellTypist uses the cell type hierarchy and meta-neighbors to further expand the search space and preserve cell type spectrum. Specifically, the cell type hierarchy, especially the cell types grouped at the high hierarchical level, contains information about what cell types are transcriptionally similar.
During the neighborhood search, each cell will search, in addition to the cell type it originates from, other cell types at the high hierarchy. Moreover, CellTypist leverages the distance matrix among cell types (see the section "Prediction of dissimilarity profiles across datasets") to find the nearest neighbors of cell types. For each cell type in each dataset, two nearest cell types will be located (the default parameter of n_meta_neighbors is three, meaning that two other cell types will be incorporated). These cell types will be next combined across datasets to form a final cell type group, against which neighbors of each cell will be searched. Due to the presence of dataset-specific cell types, some cell type groups may not cover all datasets. For these groups, the neighbors of a query cell will be contributed only by datasets (or batches) that qualify, through providing more cell neighbors per dataset.
Analysis of five immune datasets The five immune datasets were downloaded from the CellTypist database  8  and combined into one cell-by-gene count matrix. After normalization and log-transformation (with a pseudocount of 1), the highly variable genes (HVGs) were selected using the function scanpy. pp. highly_variable_genes (batch_key is set to the dataset of origin) from Scanpy  49  (version 1. 7. 2), followed by data scaling (scanpy. pp. scale) and principal component analysis (PCA) (scanpy. tl. pca). Cell type harmonization and cell reannotation were conducted by following the aforementioned pipeline in the Methods section "Cell type harmonization", collectively wrapped as the function celltypist. harmonize in CellTypist. The result was visualized as a tree plot or a cell type hierarchy using the function celltypist.
treeplot, details of which can be found in the stage V of the Methods section "Cell type harmonization". Metaanalyses of cell types across datasets were performed by averaging the cell-by-cell-type matrix 𝐺 across rows (cells) that belong to the same cell types. To confirm the existence of cell subtypes (e. g. , ABCs) harmonized and detected by CellTypist, associated cells (e. g. , memory B cells) were extracted and re-clustered using a canonical PCA-based clustering approach. Marker genes used for inspection and validation of cell subtypes were from the CellTypist immune cell encyclopedia (https://www. celltypist. org/encyclopedia/Immune/v2).
Analysis of four diseased lung datasets The four lung datasets were downloaded from respective studies  [27] [28] [29] [30]. For each dataset, cells with fewer than 200 expressed genes were discarded. We next combined cell type annotations with disease conditions to define the statuses of cells in the four datasets, and further removed the cell status with a size fewer than 20 cells. Normalization, log-transformation, HVGs selection, scaling, and PCA were performed as in the Methods section "Analysis of five immune datasets". The harmonization graph and cell type hierarchy were obtained using the function celltypist. harmonize in CellTypist, and visualized with celltypist. treeplot. Label transfer analysis (e. g. , between "Aberrant_Basaloid" and "KRT5-/KRT17+") was conducted using the logistic regression framework in CellTypist v1.
0 based on the model training (celltypist. train) followed by cell type prediction (celltypist. annotate) with the default parameters. Existence of cell subtypes harmonized and revealed by CellTypist was confirmed in a similar manner as that for the five immune datasets. The following gene signatures were used. Pulmonary fibrosis: CDKN2A, COL1A1, FN1, MMP7, MUC5B, SMAD3, ITGB6, GDF15, and EPHB2; ECM: from the Gene Ontology (GO) term of "extracellular matrix structural constituent"; profibrotic signature: SPP1, LIPA, SPARC, GPC4, PALLD, CTSK, MMP9, and CSF1. Scores of gene signatures were calculated using the Scanpy function scanpy. tl. score_genes to account for background gene expression levels. For the analysis of molecular alterations, differential expression analysis was performed using the Scanpy function scanpy.
Integration of datasets with different tools The five immune datasets or four diseased lung datasets were integrated in an annotation-aware manner using CellTypist, with the algorithmic details in the Methods section "Annotationaware data integration". The dataset of origin was considered as a batch confounder and the cell reannotation information from the harmonization pipeline was treated as a biological label key. The default parameters (e. g. , n_meta_neighbors = 3) were used. The Scanpy function scanpy. tl. umap was further utilized to generate UMAP embeddings. We also used the following eight methods to integrate the datasets. 1) FastMNN: normalized expression data (HVGs only) was used. The data was split by the dataset of origin, and served as input for the wrapper function RunFastMNN in the R package Seurat 9 (version 4.
3. 0). After this, neighbors were located using the function FindNeighbors (tag-value pairs: dims = 1:50, reduction = "mnn") and UMAP was generated with the function RunUMAP (tag-value pairs: dims = 1:50, reduction = "mnn"). 2) Harmony: normalized and scaled expression data (HVGs only) with pre-calculated PCA coordinates were used. The wrapper function RunHarmony in Seurat was used with the dataset of origin being a batch confounder. After this, neighbors were located using the function FindNeighbors (tag-value pairs: dims = 1:50, reduction = "harmony") and UMAP was generated with the function RunUMAP (tag-value pairs: dims = 1:50, reduction = "harmony"). 3) Seurat V4: normalized expression data (HVGs only) was used. The data was split by the dataset of origin, and for each dataset gene expression was scaled separately.
Next, the function FindIntegrationAnchors (tag-value pairs: normalization. method = 'LogNormalize', scale = FALSE, dims = 1:50) followed by IntegrateData (tag-value pairs: normalization. method = 'LogNormalize', dims = 1:50) were used to generate a batch-corrected expression matrix. This matrix underwent a canonical pipeline of data scaling (ScaleData), PCA (RunPCA), neighborhood search (FindNeighbors), and UMAP generation (RunUMAP). 4) BBKNN: pre-calculated PCA coordinates were used. The function bbknn. bbknn was applied to the data with the dataset of origin being a batch confounder. The resulting neighborhood graph was used to generate UMAP embeddings with the function scanpy. tl. umap in Scanpy. 5) Scanorama: scaled expression data (HVGs only) was used. The wrapper function scanpy. external.
pp. scanorama_integrate in Scanpy was used to generate the latent representations with the dataset of origin being a batch confounder, followed by neighborhood construction (scanpy. pp. neighbors) using the derived latent space and UMAP generation (scanpy. tl. umap). 6) scVI: raw count expression (HVGs only) was used with the batch covariate being the dataset of origin. A scVI model was built using the function scvi. model. SCVI (tag-value pairs: n_layers = 2, n_latent = 30, gene_likelihood = "nb") and further trained to generate the latent representations, followed by neighborhood construction (scanpy. pp. neighbors) using the derived latent space and UMAP generation (scanpy. tl. umap). 7) scANVI: raw count expression (HVGs only) was used with the batch covariate being the dataset of origin.
The function scvi. model. SCVI was first applied as in scVI, and a scANVI model was next built using the function scvi. model. SCANVI. from_scvi_model (the biological key is set as the cell annotations derived from the CellTypist harmonization procedure; both high-and low-hierarchy annotations were used). This model was trained with 20 epochs to generate the latent representations, followed by neighborhood construction (scanpy. pp. neighbors) using the derived latent space and UMAP generation (scanpy. tl. umap). 8) scGen: normalized expression (HVGs only) was used with the batch covariate being the dataset of origin and the biological covariate (high-or low-hierarchy annotations) being the cell annotations derived from the CellTypist harmonization procedure. The function scgen. SCGEN in scGen was used to build the model, which after training, resulted in the latent representations, followed by neighborhood construction (scanpy.
Benchmarking analysis of CellTypist We used the metrics summarized in scIB 50 (version 1. Simpson's index metrics were not selected, and other metrics requiring a cluster key were also not used here. Due to an elongated runtime, only 50,000 cells were randomly sampled for assessing the six metrics. Considering that not all tools produce batch-corrected lowdimensional latent representations, we used the ultimate UMAP embeddings (see the previous section), together with the neighborhood graph constructed from the UMAP, to directly assess biological conservation and batch removal at the level of UMAP. Cells from the five immune datasets were next downsampled to 10,000, 50,000, 100,000, 200,000, 300,000 or all cells (417,866). The integration procedures as detailed in the Methods section "Integration of datasets with different tools" were followed for each cell subset.