Cell type prioritization in single-cell data

Abstract

We present Augur, a method to prioritize the cell types most responsive to biological perturbations in single-cell data. Augur employs a machine-learning framework to quantify the separability of perturbed and unperturbed cells within a high-dimensional space. We validate our method on single-cell RNA sequencing, chromatin accessibility and imaging transcriptomics datasets, and show that Augur outperforms existing methods based on differential gene expression. Augur identified the neural circuits restoring locomotion in mice following spinal cord neurostimulation.

Extended Data Fig. 1. Augur overcomes confounding factors to cell type prioritization in simulated cell populations

a-b, Area under the receiver operating characteristic curve (AUC) of a random forest classifier trained in three-fold cross-validation to distinguish two simulated populations of cells, with the total number of cells increasing from n = 100 to n = 1,000 and the proportion of differentially expressed genes between the two populations varying from 0% to 100%, a, or the location parameter of the differential expression factor log-normal distribution varying from 0.1 to 1.0, b. c-d, As in a-b, but with the naive random forest classifier replaced with the subsampling procedure employed by Augur. e-f, Relationship between Augur AUC and the proportion of differentially expressed genes, e, or the location parameter of the differential expression factor log-normal distribution, f, in distinguishing two simulated populations (n = 200 cells total). The mean and standard deviation of n = 10 independent simulations are shown. Inset, two-sided Pearson correlation. g, Cell type prioritizations (AUC or number of differentially expressed genes) for a naive random forest classifier, Augur, and an exemplary single-cell differential expression test, the Wilcoxon rank-sum test, for two simulated populations of cells with 50% of genes differentially expressed and a log-normal location parameter of 0.5, with the total number of cells increasing from n = 100 to n = 1,000 cells. Like a naive random forest strategy, the number of differentially expressed genes detected by the Wilcoxon rank-sum test scales linearly with the number of cells. The mean and standard deviation of n = 10 independent simulations are shown. Dotted lines show linear regression; shaded areas show 95% confidence intervals. h-i, Number of differentially expressed genes detected by six tests for single-cell differential gene expression between two simulated populations of cells, with the total number of cells increasing from 100 to 1,000 and the proportion of differentially expressed genes between the two populations varying from 0% to 100%, h, or the location parameter of the differential expression factor log-normal distribution varying from 0.1 to 1.0, i. j, Relationship between number of differentially expressed genes detected by five tests for single-cell differential gene expression and the proportion of differentially expressed genes simulated between the two populations, for simulated populations of between 100 and 1,000 cells (see also Fig. 1e). All single-cell differential expression tests detect a larger number of differentially expressed genes in a large population of cells with modest transcriptional perturbation (20% of genes differentially expressed) than in a smaller population of cells with more profound perturbation (70% of genes differentially expressed).

Extended Data Fig. 2. Augur overcomes confounding factors to cell type prioritization in a compendium of published single-cell RNA-seq datasets

a, Overview of n = 22 published scRNA-seq datasets comparing two or more experimental conditions, used to verify the relationship between cell type prioritizations from a random forest classifier, Augur, or single-cell differential expression tests. Left, heatmap indicating the species of origin, the sequencing protocol, and whether cells or nuclei were sequenced. Right, properties of each dataset, including the total number of cell types identified in the original studies; the total number of cells sequenced; the number of cells per type (red bars indicate mean); and the mean number of reads for cells of each type. b, Pearson correlations between the AUC of each cell type, and the number of cells of that type sequenced, across 22 datasets for Augur, bottom, and a naive random forest classifier without subsampling, top, as shown in Fig. 2c. c, Pearson correlations between the number of differentially expressed genes per cell type, at 5% FDR, and the number of cells of that type sequenced, across 22 datasets for six statistical tests for single-cell differential expression. d, Number of cells in the top-ranked cell type across 22 datasets for Augur, bottom, and a naive random forest classifier without subsampling, top. e, Number of cells in the top-ranked cell type across 22 datasets for six statistical tests for single-cell differential expression. f, Jaccard index between the top-ranked 1 to 5 cell types across 22 datasets, comparing Augur and six statistical tests for single-cell differential expression. g, Cell type prioritizations in the Grubman et al., 2019 dataset by Augur and a representative test for single-cell differential expression, the Wilcoxon rank-sum test (“DE”). h, Relationship between AUC and number of differentially expressed genes per cell type, at 5% FDR, in the Grubman et al., 2019 dataset. Dotted line shows linear regression. i, Relationship between AUC and number of cells sequenced in the Grubman et al., 2019 dataset. Augur cell type prioritizations are uncorrelated with the number of cells per type. Dotted line shows linear regression; inset shows two-sided Pearson correlation. j, Relationship between number of differentially expressed genes and number of cells sequenced in the Grubman et al., 2019 dataset. Cell type prioritizations based on the number of differentially expressed genes are strongly correlated with the number of cells per type. Dotted line shows linear regression; inset shows two-sided Pearson correlation.

Extended Data Fig. 3. Augur overcomes confounding factors to cell type prioritization in simulated tissues and across single-cell modalities

a, Number of cells within each of eight cell types in a simulated tissue with increasingly unequal cell type proportions, as quantified by the Gini coefficient. b, Cell type prioritization in simulated scRNA-seq data from a tissue with 5,000 cells distributed in eight cell types, with 10-80% of genes DE in response to perturbation, and increasingly unequal numbers of cells per type (as quantified by the Gini coefficient). The correlation to simulation ground truth (proportion of DE genes) is shown for five tests for single-cell differential gene expression. The mean and standard deviation of n = 10 independent simulations are shown. Dashed line shows mean Gini coefficient of cell type frequencies across 22 published scRNA-seq datasets. **, p < 0.01; ***, p < 0.001, two-sided paired t-test. c, Inequality of cell type proportions in published scRNA-seq data. Top, Gini coefficient of cell type proportions across 22 published scRNA-seq datasets. Horizontal line and shaded area show the mean and standard deviation of the Gini coefficient across all datasets. Bottom, number of cells of each type across 22 published scRNA-seq datasets. d, Comparison of cell type prioritization in independent scRNA-seq and single cell imaging transcriptomics (STARmap) studies of the mouse visual cortex after light exposure. Left, Augur cell type prioritization in the STARmap dataset. Bottom, Augur cell type prioritization in the scRNA-seq dataset. Center, correspondence between cell types defined in the scRNA-seq and STARmap datasets, quantified as the Spearman correlation coefficient between average profiles for each cell type across 139 genes present in both datasets.

Extended Data Fig. 4. Differential cell type prioritization in single-cell RNA-seq data

a, Schematic overview of the permutation-based test for differential prioritization with Augur. First, cell type prioritization is performed within each of two conditions separately, yielding condition-specific AUCs for each cell type. Next, sample labels are randomly permuted within each cell type, and cell type prioritization is performed on shuffled data, yielding a null distribution of AUCs for each cell type and condition. AUCs for matching cell types are compared across conditions to calculate a ‘ΔAUC score’ for each cell type, and a null distribution of ΔAUC scores is calculated using the permuted data. Permutation p-values can then be calculated for each cell type, enabling the identification of statistically significant differences in cell type prioritization between conditions, as well as the condition in which the cell type is more transcriptionally separable. b, Neuron subtypes with statistically significant differences in AUC between female and male mice during parenting, in a single-cell imaging transcriptomics experiment employing multiplexed error robust fluorescence in situ hybridization (MERFISH) (n = 79 subtypes). Eleven subtypes have significantly higher AUCs in female parents, whereas two have significantly higher AUCs in male parents. c, Relationship between differential prioritization ΔAUC for parenting between male and female mice, and AUC for sex in naive mice. Several neuronal subtypes preferentially activated during parenting in female mice are also transcriptionally distinct in naive mice, such as the I-32 cluster, which is enriched for aromatase expression, and expresses multiple sex steroid hormone receptors. d, Neuron subtypes with statistically significant differences in AUC in response to whisker lesioning in Cx3cr1^+/− as compared to Cx3cr1^−/− mice, in a single-cell RNA-seq experiment (n = 28 subtypes). Four subtypes are have significantly higher AUCs in homozygous mice, whereas one subtype has a significantly higher AUC in heterozygous mice.

Extended Data Fig. 5. Cell type prioritization from transcriptional dynamics in acute experimental perturbations

a, Left, schematic overview of the scSLAM-seq workflow. Cells are exposed to the nucleoside analogue 4-thiouridine (4sU), which is incorporated during transcription and converted to a cytosine analogue by iodoacetamide prior to RNA sequencing. This labeling permits in silico deconvolution of RNA molecules transcribed before and after 4sU exposure (‘old’ and ‘new’, respectively), and calculation of the ratio of new to total RNA (NTR), an experimental analogue to the computationally determined ‘RNA velocity’^,. Right, AUCs for mouse fibroblasts exposed to lytic mouse cytomegalovirus (CMV) at 2 h post-infection, calculated by applying Augur to either total RNA or the NTR. The greater separability for the NTR reflects additional information specifically captured by the temporal dynamics of RNA expression in the context of this acute perturbation. b-e, Cell type prioritization based on exonic reads, total RNA, or RNA velocity for cells of the mouse visual cortex after exposure to light for 1 h, b-c, or 4 h, d-e, in the Hrvatin et al., 2018 dataset⁵. The AUC is significantly higher for RNA velocity than for either exonic reads (1 h, n = 34 cell types, 4 h, n = 35 cell types; two-sided paired t-tests: b, 1 h, p = 6.9 × 10^-7; d, 4 h, p = 8.2 × 10^-7) or total RNA (c, 1 h, p = 2.8 × 10^-7; e, 4 h, p = 3.0 × 10^-6), reflecting additional information specifically captured by acute transcriptional dynamics. f-g, Cell type prioritization based on exonic reads, total RNA, or RNA velocity in an Act-seq dataset, which minimizes transcriptional changes induced by single-cell dissociation. Cell types of the medial amygdala in mice subjected to 45 min of immobilization stress and control mice were profiled by Drop-seq after treatment with the transcription inhibitor actinomycin D. The AUC is higher for RNA velocity than for either exonic reads (f, p = 0.026, n = 6 cell types) or total RNA (g, p = 0.053), reflecting the additional information specifically captured by acute transcriptional dynamics, and indicating this is not an artefact related to the transcriptional perturbations known to be induced by conventional dissociation procedures. h-i, Cell type prioritization based on exonic reads, total RNA, or RNA velocity in a chronic perturbation. Cell types of the lateral hypothalamic area were profiled by Drop-seq in mice after 9-16 weeks of maintenance on either high-fat diet or control diet. No significant difference in AUCs was observed for RNA velocity compared to either exonic reads (h, p = 0.22, n = 13 cell types) or total RNA (i, p = 0.98), consistent with the time scale of the experimental perturbation.

Extended Data Fig. 6. Subclustering of single-neuron transcriptomes identifies 39 neuron subtypes in the mouse lumbar spinal cord

See also Extended Data Fig. 7a. a, Dot plot showing expression of one marker gene per cell type for the 39 neuron subtypes of the mouse lumbar spinal cord. b, Neuron subtype detection across experimental conditions (n = 6,035 neurons). TESS, targeted electrical epidural stimulation of the lumbar spinal cord. c, Proportion of neurons of each subtype detected in each experimental condition. d, Neuron subtype detection across experimental replicates (n = 3 mice per condition). e, Proportion of neurons of each subtype detected in each experimental replicate.

Extended Data Fig. 7. Robustness of Augur cell type prioritizations for mouse lumbar spinal cord neurons

a, Clustering tree of mouse spinal cord neurons over seven clustering resolutions, revealing the hierarchical relationships between spinal cord neuron subtypes. Node color reflects AUCs for cell type prioritization in targeted electrical epidural stimulation. b, AUCs for each of 37 neuron subtypes represented by at least 20 cells in both control and TESS-treated mice. c-e, Robustness of cell type prioritization for neuron subtypes of the mouse lumbar spinal cord. c, Impact of systematically withholding cells from each of six replicates (n = 3 per group) on cell type prioritization. Left, cell type prioritization with all six replicates, as in Fig. 2f. Grey tiles indicate neuron subtypes that were not represented by at least 20 cells in each condition after removal of cells from an experimental replicate. d, Impact of varying Augur parameters, including the number of subsamples and the size of each subsample; random forest-specific hyperparameters (number of trees, minimum split size, number of features sampled per split); and the choice of classifier (random forest, RF; L1-penalized logistic regression, LR) on cell type prioritization. Grey tiles indicate sample sizes larger than the number of cells of that type in the dataset. e, Impact of varying RNA velocity parameters, including exonic and intronic expression filters, the number of cells in the k-nearest neighbors pooling, and the extreme quantiles used to fit γ coefficients, on cell type prioritization.

Extended Data Fig. 8. Absence of colocalization of canonical marker genes for cell types not prioritized by Augur and Fos by RNAscope in situ hybridization

Schematic indicates imaging location for each marker within the spinal cord. Bottom, proportion of cells expressing Fos from cell types prioritized by Augur (n = 3 cell types) or not prioritized by Augur (n = 6 cell types). Cell types prioritized by Augur are significantly more likely to express Fos after walking with TESS, compared to controls (p = 0.01, two-sided Fisher’s exact test), whereas cell types not prioritized by Augur do not display a statistically significant difference (p = 0.74). Error bars show standard deviation of the sample proportion.

Extended Data Fig. 9. Impact of mean gene expression level on cell type prioritization

Cell type prioritizations were performed using both Augur and a representative single-cell differential expression method, the Wilcoxon rank-sum test, using the entire transcriptome (left column) or genes divided into five quintiles based on mean expression (right columns). Insets show two-sided Pearson correlations throughout. a, Relationship between Augur cell type prioritizations (AUC) and the proportion of differentially expressed genes between two simulated populations of cells (n = 200 cells total), as shown in Supplementary Fig. 1e. The mean and standard deviation of n = 10 independent simulations are shown. b, As in a, but with Augur applied to each quintile of gene expression separately. The AUC remains strongly correlated with the ground-truth perturbation intensity, regardless of mean expression levels (r ≥ 0.92). c, Relationship between Augur cell type prioritizations (AUC) and the location parameter of the differential expression factor log-normal distribution between two simulated populations of cells (n = 200 cells total), as shown in Supplementary Fig. 1f. The mean and standard deviation of n = 10 independent simulations are shown. d, As in c, but with Augur applied to each quintile of gene expression separately. The AUC remains strongly correlated with the ground-truth perturbation intensity, regardless of mean expression levels (r ≥ 0.95). e-f, As in a-b, but showing the number of differentially expressed genes detected by a Wilcoxon rank-sum test at 5% FDR, either across the entire transcriptome, e, or within each expression quintile, f. No differentially expressed genes are detected at 5% FDR outside of the top expression quintile. g-h, As in c-d, but showing the number of differentially expressed genes detected by a Wilcoxon rank-sum test at 5% FDR, either across the entire transcriptome, g, or within each expression quintile, h. No differentially expressed genes are detected at 5% FDR outside of the top expression quintile. i, Cell type prioritization in simulated scRNA-seq data from a tissue with 5,000 cells, distributed in eight cell types, with increasingly unequal numbers of cells per type, as quantified by the Gini coefficient and shown in Fig. 1f. The correlation to simulation ground truth (proportion of DE genes) is shown for Augur and a representative test for single-cell DE (Wilcoxon rank-sum test). The mean and standard deviation of n = 10 independent simulations are shown. j, As in i, but with both Augur and the Wilcoxon rank-sum test applied to each quintile of gene expression separately. k, Pearson correlation between Augur cell type prioritizations (AUC) and simulation ground truth (proportion of DE genes) in simulated scRNA-seq data from tissue with eight cell types, subjected to perturbations of varying intensity, as quantified by the the location parameter of the differential expression factor log-normal distribution. The mean of n = 10 independent simulations is shown for each perturbation intensity.. l, As in k, but with Augur applied to each quintile of gene expression separately. Augur incorporates information from lowly expressed genes even in subtle perturbations. m, Number of differentially expressed genes detected by a Wilcoxon rank-sum test at 5% FDR for each cell type in the Kang et al. dataset, within each expression quintile, confirming the simulations in a-l reflect trends in real data.

Extended Data Fig. 10. Impact of batch effects on cell type prioritization

Two populations of cells (n = 200 cells total) were simulated, with each condition sequenced in two batches, and varying degrees of perturbation-dependent differential expression and/or technical batch effects were introduced according to five different batch effect scenarios. For each of the five scenarios, the following panels are shown from left to right: i, Principal component analysis (PCA) of a representative simulation. ii, Correlation between AUC and magnitude of simulated batch effect with 0% of genes differentially expressed in response to perturbation, reflecting the introduction of a spurious difference between conditions where none exists (inset, two-sided Pearson correlation). iii, Correlation between AUC and magnitude of simulated batch effect when the random forest classifier is tasked with predicting batch rather than condition (AUC_batch), confirming the batch effect introduces the expected separability. iv, Correlation between proportion of genes differentially expressed in response to perturbation and AUC for simulated populations of cells with no batch effect, and batch effects of three different magnitudes. v, Cell type prioritizations in simulated populations of cells with varying perturbation intensity (% DE genes) and batch effect magnitudes. vi, As in i, but after computational batch effect correction by alignment of mutual nearest neighbors. vii, As in v, but after computational batch effect correction by alignment of mutual nearest neighbors. a, Impact of batch effects on cell type prioritization when technical batch is unconfounded with either condition or differential expression. b, Impact of batch effects on cell type prioritization when batch #1 is twice as large as batch #2. c, Impact of batch effects on cell type prioritization when perturbation-dependent differential expression is stronger in one of the two batches. d, Impact of batch effects on cell type prioritization when technical batch is mildly confounded with condition (simulated cells are overrepresented in batch 1 by a factor of 20%). e, Impact of batch effects on cell type prioritization when technical batch is moderately confounded with condition (simulated cells are overrepresented in batch 1 by a factor of 50%). f, Impact of batch effects on cell type prioritization when technical batch is severely confounded with condition (simulated cells are overrepresented in batch 1 by a factor of 80%).

Fig. 1. Augur correctly prioritizes cell types in synthetic and experimental single-cell datasets.

a, Schematic overview of Augur. b, AUCs of Augur and a naive random forest classifier without subsampling in simulated scRNA-seq datasets containing increasing numbers of cells. Cell type prioritizations are confounded by training dataset size for the naive classifier, but Augur abolishes this confounding factor. The mean and standard deviation of n = 10 independent simulations are shown. Dotted lines show linear regression; shaded areas show 95% confidence intervals. c, Pearson correlations between the AUC of each cell type, and the number of cells of that type sequenced, across a compendium of 22 scRNA-seq datasets, for Augur and a naive random forest classifier without subsampling. d, Augur AUCs scale monotonically with both the proportion of DE genes and the magnitude of DE in simulated cell populations of n = 200 cells. e, Relationship between number of DE genes detected by a representative test for single-cell differential gene expression (Wilcoxon rank-sum test), and the proportion of DE genes simulated between the two populations, for simulated populations of between n = 100 and n = 1,000 cells. f, Cell type prioritization in simulated scRNA-seq data from a tissue with 5,000 cells, eight cell types and increasingly unequal numbers of cells per type, as quantified by the Gini coefficient. The Pearson correlation to the simulation ground truth (proportion of DE genes) is shown for Augur and a representative test for single-cell DE (Wilcoxon rank-sum test). The mean and standard deviation of n = 10 independent simulations are shown. Dashed line shows mean cell type Gini coefficient across n = 22 published scRNA-seq datasets (0.52). **, p < 0.01; ***, p < 0.001, two-sided paired t-test. g, Pearson correlation between cell type prioritizations (AUC/number of DE genes) and simulation ground truth for Augur and six tests for single-cell DE in simulated tissues containing eight cell types (n = 5,000 cells) with a cell type Gini coefficient of 0.50, approximately equal to the mean of 0.52 in 22 published scRNA-seq datasets. h, Cell type prioritizations of bone marrow-derived mononuclear phagocytes from four species stimulated with LPS for between two and six hours for Augur and a representative test for single-cell DE (Wilcoxon rank-sum test). i, Pearson correlation between cell type prioritizations and duration of LPS exposure for Augur and six tests for single-cell DE (in the same order as in g). Pearson correlations with a two-sided p-value less than 0.05 are shown in orange. j, Left, Augur cell type prioritizations mirror the number of DE genes in a microarray dataset of FACS-purified cells (two-sided Pearson correlation, n = 6 cell types matching between bulk and single-cell datasets). Right, the number of DE genes detected in the scRNA-seq dataset by a Wilcoxon rank-sum test is uncorrelated with the FACS gold standard. k, Pearson correlation between cell type prioritizations and the FACS gold standard for Augur and six tests for single-cell DE (n = 6 cell types matching between bulk and single-cell datasets). Pearson correlations with a two-sided p-value less than 0.05 are shown in orange. l, Reproducibility of cell type prioritization in two independent studies of Alzheimer’s disease (n = 6 cell types matching between single-cell datasets, two-sided Pearson correlation). m, Augur cell type prioritizations in a scATAC-seq dataset track with the number of DE genes in an RNA-seq dataset of FACS-purified cells (n = 3 cell types matching between bulk and single-cell data).

Fig. 2. Augur identifies neuron subtypes that enable walking after paralysis.

a, Top, single-nucleus RNA-sequencing experimental design to prioritize neuron subtypes recruited by TESS. Middle, chronophotography of mice in the presence or absence of TESS and monoaminergic agonists. Bottom, stick diagram decompositions of right leg movements; leg endpoint trajectory with acceleration at toe-off; activity of extensor and flexor muscles of the ankle. b, Principal component analysis (n = 3 mice) of gait parameters for each condition (small circles). Large circles show the average per group. c, Bar plot shows the average scores on principal component 1 (PC1), which quantify the locomotor performance of paralyzed mice (n = 3) and mice walking with TESS (n = 3). d, Uniform manifold approximation and projection (UMAP) visualization of 18,514 nuclei, revealing the six major cell types of the mouse lumbar spinal cord. e, UMAP visualization of 6,035 neurons subjected to an additional round of sub-clustering and the 39 identified neuron subtypes. f, UMAP visualization of 6,035 neurons, colored by Augur cell type prioritization (AUC). The seven prioritized neuron subtypes with the highest AUCs are highlighted. g, Monosynaptically restricted anterograde tracing in Vsx2-Cre mice reveals V2a interneurons densely innervating motor neurons (ChAT). Similar results were obtained from three independent experiments. h, Dot plot showing expression of the immediate early gene Fos in neuron subtypes prioritized by Augur. i, Confirmation of colocalization of V2a, V1/V2b, and Spp1 marker genes (Vsx2, Slc6a5, and Spp1 respectively) and Fos by RNAscope in situ hybridization. Schematic indicates location of imaging for each marker within the spinal cord to aid specificity. Similar results were obtained from two independent experiments.