Single-cell differential expression analysis between conditions within nested settings

Single-cell differential expression analysis between conditions within nested settings Leon Hafner1, Gregor Sturm2,3 and Markus List1,4 1 Data Science in Systems Biology, TUM School of Life Sciences, Technical University of Munich, Freising, Germany 2 Biocenter, Institute of Bioinformatics, Medical University of Innsbruck, Austria 3 Boehringer Ingelheim International Pharma GmbH & Co KG, Biberach, Germany 4 Munich Data Science Institute, Technical University of Munich, Garching, Germany Motivation & Background Benchmarking Workflow & Methods Results: Runtime Analysis Results: Method Performance Conclusion References: [1] Kip D. Zimmerman, Mark A. Espeland, and Carl D. Langefeld. A practical solution to pseudoreplication bias in single-cell studies. Nature Communications, 12(1), February 2021. [2] Jordan W. Squair et al. Confronting false discoveries in single-cell differential expression. Nature Communications, 12(1), September 2021. [3] Michael I Love, Wolfgang Huber, and Simon Anders. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), December 2014. [4] Gabriel E Hoffman and Panos Roussos. Dream: powerful differential expression analysis for repeated measures designs. Bioinformatics, 37(2):192–201, July 2020. [5] Rand R. Wilcox. Introduction to Robust Estimation and Hypothesis Testing. Elsevier, 2022. [6] Greg Finak, Andrew McDavid et al. MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biology, 16(1), December 2015. [7] Romain Lopez, Jeffrey Regier, Michael B. Cole, Michael I. Jordan, and Nir Yosef. Deep generative modeling for single-cell transcriptomics. Nature Methods, 15(12):1053–1058, November 2018. [8] Shiquan Ren, Hong Lai, Wenjing Tong, Mostafa Aminzadeh, Xuezhang Hou, and Shenghan Lai. Nonparametric bootstrapping for hierarchical data. Journal of Applied Statistics, 37(9):1487–1498, September 2010. [9] Simone Tiberi, Helena L Crowell, Pantelis Samartsidis, Lukas M Weber, and Mark D Robinson. Distinct: a novel approach to differential distribution analyses. bioRxiv, November 2020. [10] Hyun Min Kang, et al. Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nature Biotechnology, 36(1):89– 94, December 2017. [11] Christina B. Azodi, Luke Zappia, Alicia Oshlack, and Davis J. McCarthy. splatPop: simulating population scale single-cell RNA sequencing data. Genome Biology, 22(1), December 2021. Key Findings: - pseudobulking followed by classical bulk RNA-Seq differential expression analysis methods is often the most effective strategy - methods like MAST cannot cope with current atlas-level datasets sizes Recommendations: - use DESeq2 for simple non-hierarchical datasets und DREAM and the Permutation Test for more complex, nested datasets Future Research: - better ground truth datasets, as current datasets lack validated findings - include different simulation methods to avoid method effectiveness due to alignment with a single simulation Scenario A: - cell count fixed to 1000 - gene count variable (100 – 10,000) Scenario B: - cell count variable (100 – 10,000) - gene count fixed to 1000 Preprocessing and pseudobulking times included if used. System time within the nextflow pipeline from process startup to finish used. Fig 8: Decision tree supporting users in selecting an appropriate method for single-cell differential expression analysis. Fig 1: Withinand between-sample Spearman correlations. Shown are correlations of cells within the same sample (orange) and of cells from two different samples (blue). The middle line represents the median, and the lower and upper limits of the boxes specify the 25% quartile and the 75% quartile. The whiskers extend to the largest/smallest correlation point, within 1.5 × interquartile range. Method adapted from Ref. [1] Fig 3: Overview of the four main scenarios simulated for this study. Fig 7: Runtime of the different methods averaged over ten independent executions on (A) a variable number of genes and (B) on a variable number of cells. Fig 5: Negative control based on raw p-value cutoff and false positive rate (FPR). In a perfect case, the observed FDR corresponds to the p-value cutoff. Values below the diagonal indicates that a method is over-conservative, while values above the diagonal indicate that a method doesn’t appropriately control for false discoveries. Fig 6: Precision-recall plot of the performance on a real data set. Performance of all methods on a single-cell dataset of 1319 B-cells stimulated with interferon beta. The corresponding Reactome pathway serves as the ground truth for determining the precision-recall curve. Real data [10] -Dataset: publicly available scRNA-Seq dataset of 1,319 B-cells from 8 patients, each one stimulated with interferon-beta and the other as a control -Ground Truth: Interferon-beta pathway, obtained from the Reactome database as a silver standard for differential expression Simulated data with R package SplatPop [11] -Dataset: Two conditions with 5 samples per group and 250 cells per sample -Atlas: more complex scenario with 3 batches, each containing 5 samples, including batch effects -Dataset with varying cell numbers: variable number of cells per sample randomly drawn from a gamma distribution -Unbalanced Conditions in Atlas: extreme imbalances in batch composition with one batch having 9 samples of condition 1 and 1 sample of condition 2 Fig 4: Performance of the methods on the simulated data scenarios over 20 independent simulations and runs each for (A) the Dataset scenario, (B) the Atlas scenario, (C) the Dataset scenario with varying cell numbers, and (D) the Atlas scenario with unbalanced conditions. The random baseline is indicated by the dashed line. We excluded distinct from those benchmarks, as it computed the lowest possible p-value for the majority of the genes and therefore did not provide useful results. Fig 2: Overview of the implemented pipeline. DESeq2 & DREAM: minimal increase in runtime with more genes and no increase with more cells due to pseudobulking Permutation Test & Hierarchical Bootstrapping: high computational cost with increasing genes, but constant with increasing cells MAST, distinct, scVI: slow increase with variable genes and strong increase with variable cells (except scVI) Table 1: Seven methods that were benchmarked, differing in their implementation, input, and statistical approach. Pseudoreplication Bias: Cells from the same sample are treated as independent observations in differential expression analysis, violating the independence assumption of statistical tests [1, 2]. →Inflated type 1 error rates (false positives) →Reduced robustness and reproducibility Solutions: Pseudobulking: Aggregating data at the sample level to treat samples as independent units [1, 2]. Hierarchical Methods: Applying models that account for the nested data structure and do not assume independence, such as linear mixed effect models and hierarchical bootstrapping [1]. Many single-cell differential expression studies do not consider the pseudoreplication bias [2]. Method Implementation Input Statistical Approach Reference DESeq2 pseudobulk parametric [3] DREAM pseudobulk parametric [4] Permutation Test pseudobulk non-parametric [5] MAST single-cell parametric [6] scVI single-cell parametric [7] Hierarchical Bootstrapping single-cell non-parametric [8] distinct single-cell non-parametric [9] Aim: Comparing different strategies for differential expression analysis using simulation and real data. How is pseudoreplication bias in differential expression analysis of scRNA-Seq data addressed? How do we differ from previous efforts? Complex scenarios (atlas-level) and imbalances with consideration of pseudoreplication bias. Dataset: Parametric methods (DESeq2, scVI and MAST) significantly stronger than resampling-based methods (Hierarchical Bootstrapping & Permutation Test) Atlas: Permutation-Test performs best, but less overall difference between the methods Dataset with varying cell numbers: Varying cell numbers have little impact (even on pseudobulk methods) and parametric methods achieve higher AUPRCs Atlas with unbalanced conditions: Resampling-based methods excel, scVI with poor performance (due to missing second covariate parameter?) Permutation Test, DREAM and DESeq2 with good predictive power. scVI, Hierarchical Bootstrapping and MAST predict more false positives than expected. Methods show significantly lower performance in the real dataset compared to the simulated scenarios. This might partially be due to the incomplete Reactome ground truth.