Pan-human Azimuth

This repository describes our pipeline for constructing a pan-tissue single-cell RNA-seq atlas that adopts an unified ontology for cell type annotations, and subsequently processing this data to train Panhuman Azimuth, a neural network-based cell type classifier. The Python implementation of our classifier can be accessed here, and the R interface here. The following sections describe how the training data is collected, reannotated, and processed, as well as other relevant preprocessing tasks such as feature selection and empty cell simulation for quality control.

1. Downloading Data download/

• cellxgene_python_download.ipynb

  • Downloads CZI objects in Python pipeline. Mostly recently ran with census release 2023-05-15.

• hubmap_datasets.ipynb

  • Downloads HuBMAP objects in python pipeline. UUID numbers come from the HuBMAP data portal.

• download_disco.sh, download/links.txt

  • Downloads DISCO dataset objects with links provided.

2. Preprocessing preprocessing/

• get_czi_reference_data.R

  • Preprocesses and merges downloaded CZI objects.

• get_azimuth_reference_data.R

  • Preprocesses available Azimuth references.

• get_gtex_reference_data.R

  • Preprocesses downloaded GTEX snRNA-seq data.

• get_hubmap_reference_data.R

  • Preprocesses and merges downloaded HuBMAP objects.

• disco_refs.R

  • Preprocesses DISCO tissue-specific references.
  • Part of the preprocessing is also done in duplicates/overlap_pipeline.R

3. Dealing with Duplicate Datasets duplicates/

• overlap_pipeline.R

  • Identifies datasets that potentially include duplicate data (i.e. meta analysies) using random projection and nearest neighbor search.
  • Performs some preprocessing of DISCO datasets, as well as metadata harmonization.

• data/

  • Table of HuBMAP IDs and associated publications
  • Table of dataset IDs and # of cells that are likely to be duplicate

4. Annotating annotating/

The following functions are called in annotating/run_annotation.R with annotating/config.csv.

Tissue-specific set up and reference mapping

• create_full_tissue_objects.R
  • Loads datasets from various sources (i.e. CZI, HuBMAP) and processes them by writing to BPCells matrices, creating directories, harmonizing tissue names, removing duplicate datasets, and ensuring metadata name consistency.
  • Reannotate all data using MapQuery with the reference dataset (i.e. DISCO or an alternative reference).
  • Sketch each source-specific assay, and join the sketched layers to return a tissue-specific Seurat object with one sketched assay.
  • The sketched assay is processed to use only "valid" genes as variable features for scaling and performing dimensionality reduction.
  • Performs per-dataset clustering on the full tissue-specific object and projects cluster IDs from the sketched to the full assay.
  • Identifies unique pairs of cell types that exist within each projected cluster in the full assay as a result of reannotation.

Perform pairwise comparisons to confirm/revise reannotations

• building_pairwise_models.R

  • For each unique pair of cell types, find markers using FindMarkers on the specified reference (usually DISCO).
  • Balance classes using SMOTE if necessary
  • Fit logistic regression model using cv.glmnet

• applying_pairwise_models.R

  • Load the full tissue-specific object and iterate through dataset-specific clusters
  • Identify a pair of cell types (majority and minority, based on number of cells) in each cluster
  • Retrieve corresponding pre-trained model from building_pairwise_models.R
  • Apply the model to the minority cells, and switch cell type annotations if probability of being the other cell type is above 0.5
  • If the minority cell type is the second most common in the cluster, also apply the model to the majority cells

Perform one-vs-all comparisons to assess reannotations

• building_one_vs_all_models.R

  • Retrieve 60 positive and 40 negative markers per cell type from the reference dataset.
  • Build models for each cell type and data source.
  • Address imbalance by duplicating minority class and adjusting majority class size, followed by applying SMOTE.
  • Train logistic regression models (cv.glmnet) on the balanced data.
  • Process the above in parallel for speed.

• applying_one_vs_all_models.R

  • For each cell type and dataset, find valid cells and select models using check_markers_by_celltype, which assesses whether cells from a specific source (i.e. HuBMAP, CZI) do not differentially express markers that are representative of its reannotated identity, and defaults to the DISCO-based model if not.
  • Predict and store probabilities as one_vs_all_score and one_vs_all_score_ogct metadata based on updated and original annotations. These scores are used to assess annotation quality for downstream filtering.

5. Training data prep train_prep/

The following scripts are run in this order to filter data, reformat labels, select features/genes, and create merged objects for training.

• process_tissue_objects.R

  • Processes each tissue-specific object based on filtering configs.
  • Filter cells based on their one_vs_all_scores, first by a 0.75 cutoff, then an optional cutoff to keep the top n cells when sorted by scores.
  • Conversion from the sparse matrix to the HDF5 file is done using the write_matrix_anndata_hdf5 function from the BPCells package, and metadata is saved separately

• select_features.R

  • Select best features for each cell type by performing Wilcoxon rank-sum tests using the presto package
  • Computes a power metric to evaluate the discriminatory strength of each feature and identifies top features for each cell type
  • Find the unique set of features across all cell type groups to use as the features to train on for that tissue object

• select_nsforest_features.py

  • Following this tutorial, run NSForest by tissue and extract the binary genes from the output to be used as training features

• process_labels.py

  • Based on a CSV file of all hierarchical cell type relationships , extracted from the DISCO cell type tree, create a mapping of parent-child relationships and build a tree
  • For each tissue, update each cell's annotation by retrieving the corresponding node in the tree and using its full path from the root node to format a hierarchical label

• merge_objects.py

  • Merges all tissue-specific objects with updated hierarchical labels into one AnnData object with metadata.

6. Quality-control qc/

To enable the classification model to simultaneously assess whether a queried cell should be filtered (i.e. low-quality/dying cells, empty droplets) for QC and cell selection alongside cell type annotation, we collected and simulated such low-quality cells to be added to our training data. These cells are annotated as "Empty" and added using the following scripts.

• collect_empty_cells.R

  • With a directory of downloaded 10X gene expression h5 files, identifies junk cells/debris by finding the inflection point in the barcode rank plot, and subsets out the junk cells/debris to be combined across all files into one Seurat object

• simulate_empty_cells.R

  • Simulates ambient RNA from an scRNA-seq experiment based on the notion that in ambient cells, the relative share of highly expressed genes increases.
  • Using Tabula Sapiens v1 tissue-specific objects, calculate a pseudobulk profile, choose a total UMI counts drawn between a lower and upper threshold, amplifying abundant gene expression by reweighting, and then simulate counts from a negative binomial distribution.
  • Includes visualization steps to help with adjusting parameters.

• merge_empty_cells.py

  • Merge both collected and simulated empty cells with the training data after processing in Step 5, while holding out some datasets for validation.