CoChat
Sign Up
Cited by 194Single-cell and spatial transcriptomicsGene expression and cancer classificationGene Regulatory Network Analysis

The DOI matches its Crossref record; 87 of the 87 references it cites were found in a registry or an index.

NEBULA is a fast negative binomial mixed model for differential or co-expression analysis of large-scale multi-subject single-cell data

Open access
He, Liang, Davila-Velderrain, Jose, Sumida, Tomokazu S., Hafler, David A., Kellis, Manolis, Kulminski, Alexander M.·2021·Communications Biology
DOI 10.1038/s42003-021-02146-6PMC PMC8155058

Published version of a preprint

Open this paper in CoChat — free accountOpen PDF

Abstract

The increasing availability of single-cell data revolutionizes the understanding of biological mechanisms at cellular resolution. For differential expression analysis in multi-subject single-cell data, negative binomial mixed models account for both subject-level and cell-level overdispersions, but are computationally demanding. Here, we propose an efficient NEgative Binomial mixed model Using a Large-sample Approximation (NEBULA). The speed gain is achieved by analytically solving high-dimensional integrals instead of using the Laplace approximation. We demonstrate that NEBULA is orders of magnitude faster than existing tools and controls false-positive errors in marker gene identification and co-expression analysis. Using NEBULA in Alzheimer’s disease cohort data sets, we found that the cell-level expression of APOE correlated with that of other genetic risk factors (including CLU, CST3, TREM2 , C1q, and ITM2B ) in a cell-type-specific pattern and an isoform-dependent manner in microglia. NEBULA opens up a new avenue for the broad application of mixed models to large-scale multi-subject single-cell data.

Overview

The claims and statistics of this paper haven’t been checked yet.

Run the check — free, no account needed. It reads the open-access text and checks the sources it cites.

A few free checks a day · papers already verified open instantly · results are public

Integrity signals

Registry

Crossref · 10.1038/s42003-021-02146-6

Cites

87 of 87 references verified0 citation errors · 0 not checkable

Statistics

Figures

Fig. 1 Overview of the input, output, and computational efficiency of NEBULA. a A flowchart of NEBULA including the input data, the major estimation steps in NEBULA, and the analyses that were conducted using NEBULA in the application to the real data5,27,78. b The computational time (measured in log10(seconds)) of fitting an NBMM for 10,000 genes with respect to CPS by NEBULA, glmer.nb18, and glmmTMB19. The number of subjects was set at 50. c The computational time (measured in log10(seconds)) of fitting an NBMM for 10,000 genes with respect to the number of subjects. The error bars represent one standard deviation of n = 10,000 genes. The CPS value was set at 200. Two fixed-effects predictors were included in the NBMM. The average benchmarks were summarized from scenarios of varying subject-level and cell-level overdispersions and the CPC value of a gene ranging from exp(−4) to 1. CPS: cells per subject. CPC: counts per cell.
Fig. 1 Overview of the input, output, and computational efficiency of NEBULA. a A flowchart of NEBULA including the input data, the major estimation steps in NEBULA, and the analyses that were conducted using NEBULA in the application to the real data5,27,78. b The computational time (measured in log10(seconds)) of fitting an NBMM for 10,000 genes with respect to CPS by NEBULA, glmer.nb18, and glmmTMB19. The number of subjects was set at 50. c The computational time (measured in log10(seconds)) of fitting an NBMM for 10,000 genes with respect to the number of subjects. The error bars represent one standard deviation of n = 10,000 genes. The CPS value was set at 200. Two fixed-effects predictors were included in the NBMM. The average benchmarks were summarized from scenarios of varying subject-level and cell-level overdispersions and the CPC value of a gene ranging from exp(−4) to 1. CPS: cells per subject. CPC: counts per cell.CC BY · PubMed Central
Fig. 2 Comparison of estimated cell-level and subject-level overdispersion parameters between NEBULA-LN and NEBULA-HL. a The cell-level overdispersion estimated by NEBULA-LN and NEBULA-HL under different combinations of CPS, β0 (a lower β0 corresponding to a lower CPC value), and ϕ. The number of subjects was set at 50. b The subject-level overdispersion estimated by NEBULA-LN and NEBULA-HL under different combinations of β0, σ², and ϕ. The number of subjects was set at 50. The CPS value was set at 400. c The subject-level overdispersion estimated by NEBULA-LN and NEBULA-HL under different combinations of the number of subjects, σ², and ϕ. The CPS value was set at 400, and β0 was set at 0.05. d The subject-level overdispersion estimated by NEBULA-LN and NEBULA-HL under different combinations of CPS, σ², and ϕ. The number of individuals was set at 50, and β0 was set at 1. The summary statistics were calculated from n = 500 simulated replicates in each of the scenarios. The estimates from NEBULA-HL presented in these panels were based on the first-order Laplace approximation, and the results comparing the first-order and higher-order LA methods in NEBULA-HL were presented in Supplementary Fig. S17. CPC: counts per cell. Center line: median. Box limits: upper and lower quartiles.
Fig. 2 Comparison of estimated cell-level and subject-level overdispersion parameters between NEBULA-LN and NEBULA-HL. a The cell-level overdispersion estimated by NEBULA-LN and NEBULA-HL under different combinations of CPS, β0 (a lower β0 corresponding to a lower CPC value), and ϕ. The number of subjects was set at 50. b The subject-level overdispersion estimated by NEBULA-LN and NEBULA-HL under different combinations of β0, σ², and ϕ. The number of subjects was set at 50. The CPS value was set at 400. c The subject-level overdispersion estimated by NEBULA-LN and NEBULA-HL under different combinations of the number of subjects, σ², and ϕ. The CPS value was set at 400, and β0 was set at 0.05. d The subject-level overdispersion estimated by NEBULA-LN and NEBULA-HL under different combinations of CPS, σ², and ϕ. The number of individuals was set at 50, and β0 was set at 1. The summary statistics were calculated from n = 500 simulated replicates in each of the scenarios. The estimates from NEBULA-HL presented in these panels were based on the first-order Laplace approximation, and the results comparing the first-order and higher-order LA methods in NEBULA-HL were presented in Supplementary Fig. S17. CPC: counts per cell. Center line: median. Box limits: upper and lower quartiles.CC BY · PubMed Central
Fig. 3 Empirical type I error rate of testing cell-level and subject-level variables using NEBULA-LN and NEBULA-HL. a The empirical type I error rate of testing a cell-level variable using NEBULA-LN (blue) and NEBULA-HL (red) under different combinations of CPS, σ², and ϕ. The number of subjects was set at 50. b The empirical type I error rate of testing a subject-level variable using NEBULA-LN (blue) and NEBULA-HL (red) under different combinations of the number of subjects, σ², and ϕ. The CPS value was set at 400. The type I error rate was calculated from n = 500 simulated replicates in each of the scenarios and was evaluated at the significance level of 0.05 (the dashed lines). We used the higher-order Laplace approximation in NEBULA-HL for scenarios of low CPC. CPS: cells per subject. CPC: counts per cell.
Fig. 3 Empirical type I error rate of testing cell-level and subject-level variables using NEBULA-LN and NEBULA-HL. a The empirical type I error rate of testing a cell-level variable using NEBULA-LN (blue) and NEBULA-HL (red) under different combinations of CPS, σ², and ϕ. The number of subjects was set at 50. b The empirical type I error rate of testing a subject-level variable using NEBULA-LN (blue) and NEBULA-HL (red) under different combinations of the number of subjects, σ², and ϕ. The CPS value was set at 400. The type I error rate was calculated from n = 500 simulated replicates in each of the scenarios and was evaluated at the significance level of 0.05 (the dashed lines). We used the higher-order Laplace approximation in NEBULA-HL for scenarios of low CPC. CPS: cells per subject. CPC: counts per cell.CC BY · PubMed Central
Fig. 4 Decomposed cell-level and subject-level overdispersions in major neural cells. a The cell-level overdispersion and b the subject-level overdispersion of 16,207 genes in the excitatory neurons estimated by NEBULA across different CPC values. No covariates other than the intercept were included in the model. The color indicates the two-dimensional kernel density computed using the kde2d R function with n = 100. c The cell-level overdispersion of 16,207 genes in the excitatory neurons estimated by NEBULA versus their subject-level overdispersion. d The estimated cell-level overdispersion and e the estimated subject-level overdispersion in each of the six brain cell types before and after adjusting for age, sex, race, AD status, the total number of features of the cell, and the percentage of ribosomal genes. Only genes with CPC > 0.1% in the relevant cell type were included in the above results. The cells in this analysis are from the 48-subject snRNA-seq data set5 in the frontal cortex.
Fig. 4 Decomposed cell-level and subject-level overdispersions in major neural cells. a The cell-level overdispersion and b the subject-level overdispersion of 16,207 genes in the excitatory neurons estimated by NEBULA across different CPC values. No covariates other than the intercept were included in the model. The color indicates the two-dimensional kernel density computed using the kde2d R function with n = 100. c The cell-level overdispersion of 16,207 genes in the excitatory neurons estimated by NEBULA versus their subject-level overdispersion. d The estimated cell-level overdispersion and e the estimated subject-level overdispersion in each of the six brain cell types before and after adjusting for age, sex, race, AD status, the total number of features of the cell, and the percentage of ribosomal genes. Only genes with CPC > 0.1% in the relevant cell type were included in the above results. The cells in this analysis are from the 48-subject snRNA-seq data set5 in the frontal cortex.CC BY · PubMed Central
Fig. 5 Identification of cell marker genes using NEBULA. a A UMAP plot of the cells in the PBMC scRNA-seq data27. The automatic cell-type annotation was performed using NEBULA and scCATCH29. The corrected or refined cell type by the manual annotation was shown after the arrow. Cluster 16 and 22 were not annotated because of either too few cells or marker genes. b Comparison between log(FC) estimated by NEBULA in the scRNA-seq data27 and the effect size estimated by glm.nb in the cell-sorting bulk RNA-seq data30 in the analysis of differentially expressed genes between naive CD8+ and CD4+ T cells. Genes having both log(FC) and an effect size >2.5 were highlighted in blue for naive CD8+ and red for naive CD4+ T cells. Comparison between the log(FC) and p-values from c NEBULA and glmer.nb and d NEBULA and glm.nb in the differential expression analysis of detecting marker genes for Cluster 1 (memory CD4+ T cells) in the scRNA-seq data27. e Q–Q plots of p-values for testing the association between simulated counts and cell clusters using NEBULA and a simple negative binomial regression (shown as NB in the plot). The counts were simulated under the null model (i.e., no association between the counts and cell clusters) based on an NBMM with the subject-level overdispersion (σ²) ranging from 0 to 1 and a cell-level overdispersion fixed at 1. f Comparison of the p-values of the top 30 genes identified by glm.nb associated with each of the four astrocyte subpopulations (AST0-AST3) in the snRNA-seq data5 with those reported by NEBULA. The labels of the four astrocyte subpopulations are corresponding to those presented in ref. 5. Subject OD: the subject-level overdispersion estimated by NEBULA. Genes highlighted in red showed substantially reduced significance in NEBULA compared to the p-values in the simple negative binomial regression.
Fig. 5 Identification of cell marker genes using NEBULA. a A UMAP plot of the cells in the PBMC scRNA-seq data27. The automatic cell-type annotation was performed using NEBULA and scCATCH29. The corrected or refined cell type by the manual annotation was shown after the arrow. Cluster 16 and 22 were not annotated because of either too few cells or marker genes. b Comparison between log(FC) estimated by NEBULA in the scRNA-seq data27 and the effect size estimated by glm.nb in the cell-sorting bulk RNA-seq data30 in the analysis of differentially expressed genes between naive CD8+ and CD4+ T cells. Genes having both log(FC) and an effect size >2.5 were highlighted in blue for naive CD8+ and red for naive CD4+ T cells. Comparison between the log(FC) and p-values from c NEBULA and glmer.nb and d NEBULA and glm.nb in the differential expression analysis of detecting marker genes for Cluster 1 (memory CD4+ T cells) in the scRNA-seq data27. e Q–Q plots of p-values for testing the association between simulated counts and cell clusters using NEBULA and a simple negative binomial regression (shown as NB in the plot). The counts were simulated under the null model (i.e., no association between the counts and cell clusters) based on an NBMM with the subject-level overdispersion (σ²) ranging from 0 to 1 and a cell-level overdispersion fixed at 1. f Comparison of the p-values of the top 30 genes identified by glm.nb associated with each of the four astrocyte subpopulations (AST0-AST3) in the snRNA-seq data5 with those reported by NEBULA. The labels of the four astrocyte subpopulations are corresponding to those presented in ref. 5. Subject OD: the subject-level overdispersion estimated by NEBULA. Genes highlighted in red showed substantially reduced significance in NEBULA compared to the p-values in the simple negative binomial regression.CC BY · PubMed Central
Fig. 6 Analysis of cell-level co-expression of APOE in astrocytes and microglia. a and b Top 20 genes whose expression in astrocytes (panel a) and microglia (panel b) were correlated with the expression of APOE at the single-cell level. The correlation measured by log(FC) was obtained by NEBULA in two snRNA-seq data sets5,78 in the frontal cortex from the ROSMAP. c and d Top 10 KEGG pathways in which the genes most correlated with APOE in astrocytes (panel c) and microglia (panel d) were enriched. The size of the circles indicates how many genes were found in that pathway. The color indicates the p-value of the enrichment.
Fig. 6 Analysis of cell-level co-expression of APOE in astrocytes and microglia. a and b Top 20 genes whose expression in astrocytes (panel a) and microglia (panel b) were correlated with the expression of APOE at the single-cell level. The correlation measured by log(FC) was obtained by NEBULA in two snRNA-seq data sets5,78 in the frontal cortex from the ROSMAP. c and d Top 10 KEGG pathways in which the genes most correlated with APOE in astrocytes (panel c) and microglia (panel d) were enriched. The size of the circles indicates how many genes were found in that pathway. The color indicates the p-value of the enrichment.CC BY · PubMed Central

Verification

Claims the paper makes

The claims and statistics of this paper haven’t been checked yet.

Run the check — free, no account needed. It reads the open-access text and checks the sources it cites.

A few free checks a day · papers already verified open instantly · results are public

Claims other papers make about this one

Each is checked against this paper’s own text: does it say what the citing paper says it does?

Self-check

Percentages, sums and reported differences re-computed from the paper’s own numbers.

The self-check needs the paper’s full text.

Citations

Cites

Show all 85 references

Research works printed but not linked: 1

  1. Huber, P. J. Robust Statistics. (John Wiley & Sons, 2004).

Cited by

Sorted by citations

See all 194

Related

Metadata

Authors
He, L., Davila-Velderrain, J., Sumida, T. S., Hafler, D. A., Kellis, M., & Kulminski, A. M.
Affiliations
  • Duke University (US)
  • Yale University (US)
  • Broad Institute (US)
Corresponding author
  • He, Liang
  • Kellis, Manolis
  • Kulminski, Alexander M.
ORCID
Venue
Communications Biology
Publisher
Springer Science and Business Media LLC
Published
2021-05-26
Volume
4
Issue
1
Pages
629-629
ISSN
2399-3642 (ISSN-L)
DOI
10.1038/s42003-021-02146-6
PMC
PMC8155058
OpenAlex
W3164861786
Cited by
181
References
85
Citations by year
2026: 42 · 2025: 40 · 2024: 48 · 2023: 27 · 2022: 19 · 2021: 4
FWCI
8.67
Type
journal article
Language
en
Open access
Licence
Locations
  • Communications Biology (publishedVersion, cc-by)
  • bioRxiv (Cold Spring Harbor Laboratory) (acceptedVersion, cc-by-nc-nd)
  • PubMed Central (publishedVersion, cc-by)
  • PubMed (publishedVersion)
  • DOAJ (DOAJ: Directory of Open Access Journals) (submittedVersion, cc-by-sa)
  • DSpace@MIT (Massachusetts Institute of Technology) (submittedVersion, cc-by-nc)
Funders
  • National Institutes of Health (5r01ag070488-05, UL1TR001863)
  • National Institute on Aging (R01-AG070488, R01AG017917, AG058002, R01AG030146, R01AG058002, R01AG062335, R01AG036836, RF1AG054012, R01AG065477, R01AG048015, U01AG046152, R01AG036042, R01AG015819, P30AG010161, RF1AG062377, U01AG032984, U01AG061356, RC2AG036547, U01AG046161, R01AG061853, RF1AG057473, R01AG067151)
  • National Institute of Mental Health (U01MH119509, R01MH109978)
  • National Cancer Institute (R01CA227473, P50CA121974)
  • National Institute of Allergy and Infectious Diseases (R01AI072222, P01AI073748, U24AI118672, P01AI039671, U19AI089992)
  • National Institute of Neurological Disorders and Stroke (R25NS079193, U01NS110453, UG3NS115064)
  • National Center for Advancing Translational Sciences
Topics
  • Life Sciences › Biochemistry, Genetics and Molecular Biology › Molecular Biology › Single-cell and spatial transcriptomics
  • Life Sciences › Biochemistry, Genetics and Molecular Biology › Molecular Biology › Gene expression and cancer classification
  • Life Sciences › Biochemistry, Genetics and Molecular Biology › Molecular Biology › Gene Regulatory Network Analysis
Keywords
nebulae, differential expression analysis, co-expression analysis, large sample approximation, computational efficiency, marker gene identification, Alzheimer's disease, apolipoprotein E, microglia, cell type-specific gene expression, isoform expression, genetic risk factors
MeSH terms
Alzheimer Disease, Apolipoproteins E, Humans, Models, Statistical, Gene Expression, Binomial Distribution, Microglia, Computational Biology, Gene Expression Profiling, Single-Cell Analysis

About these checks

Registry, retraction, Cites and statistics results are computed, not judged, with one exception: a printed reference that no index lists is read once by a model before it counts as a citation error. Claims are extracted by a model and kept only when their quote appears verbatim in the text; each cited source is then checked, with a model’s help, for whether it says what the claim relies on. Those verdicts can be wrong — read the claim’s quote and the cited source before relying on one.

Statistics and arithmetic are recomputed, not judged.

Confidence bands are shown only for fields calibrated on our annotated test set. None is calibrated yet, so no band is shown.

What counts as a claim ›How we measure confidence ›How we check this ›

Are you an author, or do you see a problem with this page? Request a review