Modelling group heteroscedasticity in single-cell RNA-seq pseudo-bulk data

Abstract

Group heteroscedasticity is commonly observed in pseudo-bulk single-cell RNA-seq datasets and when not modelled appropriately, its presence can hamper the detection of differentially expressed genes. Most bulk RNA-seq methods assume equal group variances which will under- and/or over-estimate the true variability in such datasets. We present two methods that account for heteroscedastic groups, namely voomByGroup and voomWithQualityWeights using a blocked design (voomQWB). Compared to current gold-standard methods that do not account for heteroscedasticity, we show results from simulation studies and various experiments that demonstrate the superior performance of both voomByGroup and voomQWB in error control and power when group variances in pseudo-bulk scRNA-seq data are unequal. We recommend the use of either of these methods over established approaches, with voomByGroup having the advantage of accurate variance estimation since group variance trends can take on different "shapes", whilst voomQWB has the advantage of catering to complex study designs.

Group Heteroscedasticity - A Silent Saboteur of Power and False Discovery in RNA-Seq Differential Expression

“…BCV, a key indicator of biological variability in gene expression among biological replicates, is widely used for estimating gene expression variance in RNA-Seq data 45 . Following previous studies 72 , here we demonstrate group heteroscedasticity by comparing the variability in group-specific BCV values and visualize this effect in the form of distinctive positions and shapes of group-specific mean-variance trend curves across these twelve diverse datasets ( Fig. S1 ).…”

“…1 ). Our approach mirrors that of a previous study 72 which identified group heteroscedasticity through higher biological coefficient of variation (BCV) values in one group compared to another. BCV, a key indicator of biological variability in gene expression among biological replicates, is widely used for estimating gene expression variance in RNA-Seq data 45 .…”

“…Two of these methods were based on negative binomial (NB) assumptions, which include the widely employed DESeq2 40 and edgeR.LRT (edgeR + likelihood-ratio test) 53 , while two operate under a nonparametric setup, including dearseq 21 and the classical Wilcoxon rank-sum test. The remaining two methods involve two variants of the limma model for differential expression analysis of RNA-Seq data, referred to as voom 33 (limma + voom transformation) and voomByGroup 72 (limma + group-specific trends).…”

“…Another study by Chen et al 10 highlighted unequal group variance in single-cell data and emphasized the under-utilization of tests accounting for varying variances. More recently, a study by You et al 72 highlighted that group heteroscedasticity is present in pseudo-bulk data and subsequently proposed a statistical model to account for it. However, the extent to which heteroscedasticity impacts false discovery and power, especially in large samples, remains to be investigated.…”

Population-level comparisons of gene regulatory networks modeled on high-throughput single-cell transcriptomics data

“…Pseudo-bulk profiles are frequently calculated in differential gene expression analysis to take into account biological variation between samples (16, 17). However, the biological variability between transcription factors and their target gene interactions must be accurately modeled across multiple samples to identify consistent mechanistic patterns causing phenotypic changes across samples within a population (15, 18, 19).…”