一个R包搞定22种单细胞基因差异分析

参考文献：Confronting false discoveries in single-cell differential expression

杂志：Nature Communication

多细胞组织内 mRNA 丰度的变化被细胞类型的不同反应和这些细胞类型的相对丰度的变化所混淆。 因此，批量 RNA-seq 的分辨率不足以表征对生物扰动的多方面响应。单细胞 RNA-seq (scRNA-seq) 能够在单个细胞的分辨率下量化 RNA 丰度。单细胞技术的成熟现在能够对复杂组织内的细胞状态进行大规模比较，从而提供适当的分辨率来剖析细胞类型对扰动的特异性反应。单细胞数据的稀疏性和异质性最初鼓励开发专门的统计方法来识别差异表达的 mRNA。用于差异表达分析的统计方法的激增促使研究人员探究哪种方法产生的生物学结果最准确。

结果

A ground-truth resource to benchmark single-cell differential expression.

我们旨在比较可用的基于差异表达 (DE) 分析的统计方法生成生物学准确结果的能力。 我们推断，在已知实验基本事实的真实数据集中执行这种比较将忠实地反映这些方法的性能差异，同时避免模拟数据的缺点。 我们假设，可以从匹配的批量和 scRNA-seq 中获得与该基本事实最接近的近似值，该序列对相同的纯化细胞群进行，暴露于相同的扰动，并在相同的实验室中进行测序。 对文献的广泛调查确定了总共 18 个符合这些标准的“黄金标准”数据集。 这使我们能够在已知基本事实的实验环境中对 DE 方法进行大规模比较。

img

Pseudobulk methods outperform generic and specialized single-cell DE methods.

我们总共选择了14 种 DE 方法进行比较，这些方法代表了单细胞转录组学中使用最广泛的统计方法（方法，“差异表达分析方法”）。 近 90% 的近期研究都使用了这些方法。 我们根据批量与 scRNA-seq 数据集的 DE 结果之间的一致性评估了每种方法的相对性能。 为了量化这种一致性，我们计算了批量与 scRNA-seq 数据集结果之间的一致性曲线下面积 (AUCC)。

我们在 18 个黄金标准数据集的整个概要中比较了 14 种方法的性能。 该分析立即表明，所有六种表现最佳的方法都具有共同的分析特性。 在应用统计检验之前，这些方法在生物复制品中聚集细胞，形成所谓的“假体”。 相比之下，比较单个细胞的方法表现不佳。 pseudobulk方法和单细胞方法之间的差异非常显着，并且对用于量化一致性的方法具有鲁棒性。 此外，与匹配的蛋白质组学数据的比较表明，pseudobulk 方法还可以更准确地预测蛋白质丰度的变化

img img

我们探究DE 方法之间的差异是否也会影响转录组实验的功能解释。 为此，我们比较了批量基因本体论 (GO) terms富集分析与 scRNA-seq DE。 我们发现pseudobulk方法再次更忠实地反映了真实情况，如批量 RNA-seq 中捕获的那样。 例如，在比较用 poly(I:C)（一种合成双链 RNA）刺激的小鼠吞噬细胞时，单细胞方法未能确定相关的 GO terms。

Single-cell DE methods are biased towards highly expressed genes（单细胞 DE 方法偏向于高度表达的基因）

pseudobulk方法的意想不到的优越性迫使我们研究负责它们重述生物基本事实的能力的机制。 为了研究这些机制，我们制定并测试了几个可能解释这些性能差异的假设。

先前的研究表明，对于高度表达的基因，关于 DE 的推断通常更准确。 单细胞中基因表达的测量本质上是稀疏的。 通过在每个重复中聚集细胞，pseudobulk 方法显着减少了数据中的零数量，特别是对于低表达的基因。 因此，我们最初假设，pseudobulk 方法和单细胞方法之间的准确性差异可以通过pseudobulk 方法在低表达基因中的优越性能来解释。

img

• Schematic illustration of the creation of ‘pseudobulks’ from single-cell data. Top, biological replicate from which each cell was obtained. Bottom, simulated gene expression matrix. Read counts for each gene are aggregated across cells of a given type within each biological replicate。

为了验证这一假设，我们将基因分配到三个大小相同的bins中，包括低表达、中度和高度表达的基因。 然后，我们重新计算了每个 bin 内批量和 scRNA-seq DE 之间的一致性。 与我们的预测相反，我们观察到低表达基因的pseudobulk方法和单细胞方法之间的差异很小。 相反，pseudobulk方法和单细胞方法之间最显着的差异出现在高度表达的基因中。

img

• Mean AUCCs across eighteen ground-truth datasets after dividing the transcriptome into terciles of lowly, moderately, or highly expressed genes

这一意想不到的结果让我们质疑单细胞 DE 方法是否会对高表达基因产生系统性错误。 为了探索这种可能性，我们仔细检查了大量数据集，以识别 scRNAseq 数据中的每种方法都错误地称为 DE 的基因。 我们发现由单细胞 DE 方法识别的假阳性比由pseudobulk方法识别的假阳性表达更高。 相反，单细胞 DE 方法忽略的假阴性往往表达较低。 总之，这些发现暗示了单细胞方法将高度表达的基因识别为 DE 的系统趋势，即使它们的表达保持不变。

img

为了通过实验验证这一结论，我们分析了一个数据集，其中将合成 mRNA populations掺入每个含有单个细胞的孔中。 因此，这些单细胞文库中的每一个都包含相同浓度的每种合成 mRNA。 我们发现单细胞方法错误地将许多丰富的spike-ins 识别为DE。 相比之下，pseudobulk方法避免了这种偏差。

img

• d Spearman correlation between the mean expression of 80 ERCC spike-ins expressed in at least three cells and the –log10 p-value of differential expression assigned by each DE method. e Scatterplots of mean ERCC expression vs. –log10 p-value for exemplary single-cell and pseudobulk DE methods. Trend lines and shaded areas show local polynomial regression and the 95% confidence interval, respectively.

然后我们query这种偏差在单细胞转录组学中是否普遍存在。 我们汇编了 46 个 scRNA-seq 数据集的compendium，其中包含不同的物种、细胞类型、技术和生物扰动。 我们发现单细胞 DE 方法在整个compendium中显示出对高表达基因的系统偏好。

img

• 注：Rank and statistical significance of the GO term GO:0043330 (“response to exogenous dsRNA”) in GSEA analyses of mouse bone marrow mononuclear cells stimulated with poly-I:C, a type of synthetic dsRNA, for four h, using the output of fourteen DE methods. Source data are provided as a Source Data file.

Together, these experiments suggest that the inferior performance of single-cell methods can be attributed to their bias towards highly expressed genes.

img img

DE analysis of single-cell data must account for biological replicates.

这些发现意味着pseudobulk方法拥有一个共同的分析特性，可以避免这种偏差。 我们进行了一系列实验来确定这种机制。

经过多年的发展，用于识别pseudobulk数据（即 edgeR、DESeq2 和 limma）中的 DE 基因的统计工具已经得到改进。 因此，我们query这些方法是否结合了独立于跨细胞聚合基因表达程序的固有优势。 为了测试这种可能性，我们禁用了聚合过程并将pseudobulk方法应用于单个细胞。 引人注目的是，这个过程取消了pseudobulk方法的优越性。 对高表达基因的偏见的出现与这种表现的下降并行。

img

• 注：a Schematic illustration of the experiment shown in b, in which the aggregation procedure was disabled and pseudobulk DE methods were applied to individual cells。b Left, AUCC of the original fourteen DE methods, plus six pseudobulk methods applied to individual cells, in the eighteen ground-truth datasets. Right, Spearman correlation between ERCC mean expression and –log10 p-value assigned by six pseudobulk DE methods, before and after disabling the aggregation procedure.

这一结果提出了聚合过程本身直接导致pseudobulk方法优越性的可能性。 为了评估这个概念，我们将聚合程序应用于随机的细胞组，这产生了一个由“pseudo-replicates”组成的pseudobulk矩阵。这个实验导致pseudobulk方法的性能出现类似的下降，并结合了 偏向于高表达基因.

img img img

我们试图了解可以解释这两个实验中pseudobulk方法性能下降的共同因素。 我们认识到这两个实验都会丢失有关生物复制的信息。 聚集随机的细胞组形成pseudo-replicates，或在单个细胞的比较中完全忽略复制，都引入了对高表达基因的偏见和相应的性能损失。

在相同的实验条件下，重复表现出基因表达的内在差异，这反映了生物学和技术因素。 我们推断，未能考虑到这些差异可能会导致方法将重复之间的固有可变性错误地归因于扰动的影响。 为了研究这种潜在的机制，我们比较了假体和假重复中每个基因表达的差异。 最初，我们在用 poly-I:C12 刺激的骨髓单核细胞数据集中进行了这种比较。 我们发现重排复制品会导致基因表达方差的系统性降低，影响 98.2% 的基因。 我们接下来测试了这种方差的减少是否在我们的 46 个数据集概要中系统地发生。 每次比较都显示出基因表达方差的相同减少

img

• 注：e Variance of gene expression in pseudobulks formed from biological replicates and pseudo-replicates in mouse bone marrow mononuclear cells stimulated with poly-I:C. Shuffling the replicate associated with each cell produced a systematic decrease in the variance of gene expression. Right, pie chart shows the proportion of genes with increased or decreased variance in pseudo-replicates, as compared to biological replicates. f Decreases in the variance of gene expression in pseudo-replicates as compared to biological replicates across 46 scRNA-seq datasets. Points show the mean variance in biological replicates; arrowheads show the mean variance in pseudo-replicates。

基因表达方差的减少导致统计测试将基因表达的微小变化归因于扰动的影响。 例如，在 poly-I:C 数据集中，未能考虑重复中 Txnrd3 的可变表达导致将该基因误认为差异表达。 此外，我们发现高表达基因在假重复中表现出最大的方差下降，从而解释了单细胞方法对高表达基因的偏见

img

• g Left, expression of the gene c in biological replicates (points) and pseudo-replicates (arrowheads) from unstimulated cells and cells stimulated with poly-I:C, with the range of possible pseudo-replicate expression values shown as a density. Right, mean (horizontal line) and variance (shaded area) of Txnrd3 expression in biological replicates (left) and pseudo-replicates (right). P-values were calculated by edgeR-LRT.

这一系列实验共同揭示了pseudobulk方法出乎意料的优越性的原理。 差异表达的统计方法必须考虑生物复制的内在变异性，以在单细胞数据中产生生物学上准确的结果。 考虑到这种可变性，pseudobulk 方法可以正确识别由生物扰动引起的基因表达变化。 相比之下，未能解释生物复制会导致单细胞方法系统地低估基因表达的差异。 这种对方差的低估使单细胞方法偏向于高度表达的基因，从而损害了它们产生生物学准确结果的能力。

False discoveries in single-cell DE.

我们意识到，如果不考虑生物学重复之间的差异可能会在存在真实生物学扰动的情况下产生错误发现，那么在没有任何生物学差异的情况下也可能出现错误发现。为了测试这种可能性，我们在没有任何组间差异的情况下模拟了重复之间具有不同程度异质性的单细胞数据。我们将每个重复随机分配到人工“对照”或“处理”组，并在两种条件之间测试 DE。引人注目的是，单细胞方法在没有任何扰动的情况下鉴定了数百个 DE 基因。此外，根据我们对单细胞 DE 方法失败背后机制的理解，被错误称为 DE 的基因是那些在重复之间表达变化最大的基因。 Pseudobulk 方法消除了 DE 基因的错误检测。然而，创建pseudo-replicates导致虚假 DE 基因的重新出现，进一步证实了准确 DE 分析的要求。当额外的重复被引入数据集时，错误发现的数量减少了。相比之下，在模拟数据中引入额外的细胞只会加剧潜在的问题

这些发现迫使我们调查在真实的单细胞数据中是否会出现类似的错误发现。 为了探索这种可能性，我们最初分析了暴露于干扰素 5 的人外周血单核细胞 (PBMC) 的数据集。 我们提取了未接触过干扰素的对照样品，并将它们随机分成两组。 然后我们进行了 DE 分析。 未能解释生物复制的内在可变性在随机分配的复制之间产生了数百个 DE 基因

对这些错误发现感到不安，我们好奇这一观察是否反映了一个普遍的陷阱。 为了全面解决这个问题，我们确定了总共 14 个数据集，其中在控制条件下至少包含 6 个重复。 与之前的实验一样，我们将这些未受干扰的样本随机分为合成对照组和处理组，然后在这两组之间进行 DE 分析。 该系统分析证实，与pseudobulk方法相比，单细胞方法产生了系统性过多的假阳性。 尽管完全没有生物扰动，但产生的 DE 基因富集了数百个基因本体论 (GO) terms。 此外，我们再次证实，被错误识别为 DE 的基因对应于重复之间具有最高变异性的基因 。

总之，这些实验揭示了单细胞转录组学中 DE 分析的一个基本缺陷。 然而，我们的直觉是，这个陷阱可能会影响任何从每个生物复制中获得许多观察结果的技术。 例如，我们预计空间转录组学数据中也会出现错误发现。 为了测试这一预测，我们分析了一个空间转录组学数据集，该数据集对来自肌萎缩侧索硬化 (ALS) 模型的脊髓进行了分析。 我们将来自对照小鼠的数据随机分为两组，并在脊髓的每个区域内进行 DE。 未能解释生物复制之间变异性的统计方法在每个区域内鉴定了数千个 DE 基因。 相比之下，pseudobulk方法消除了这些错误发现。

这些实验表明，生物复制之间的变异性会混淆受生物扰动影响的基因的识别。在动物模型中，可以最大限度地减少在重复之间产生这种变异性的许多因素，包括遗传背景、环境、生物扰动的强度和时间以及样品处理。相比之下，在涉及人类受试者的实验中，这些变异来源本质上更难以控制。这种区别提出了一种可能性，即人体组织的单细胞研究在生物复制之间表现出更大的变异性，因此更容易受到错误发现的影响。为了系统地评估这种可能性，我们计算了 41 个人和小鼠 scRNA-seq 数据集中重复之间的变异性。与我们的假设一致，我们检测到人类数据集中重复之间的变异性明显更大。虽然我们表明考虑生物重复对于任何 DE 分析都是至关重要的，但这一结果强调了在人体组织的单细胞研究中解决这个问题的最重要的意义。

img img img

True and false discoveries in the injured mouse spinal cord

我们最终试图证明 DE 分析可以在以前未探索的生物组织中产生真假发现的程度。 为此，我们描述了脊髓损伤 (SCI) 对损伤下方细胞中基因表达的影响。 我们特别关注腰椎脊髓，因为该区域经历了多方面的变化，导致神经元功能不可逆转的退化。

我们在遭受中胸脊髓严重挫伤的小鼠身上进行了实验。 全身运动学的多因素量化揭示了小鼠产生运动能力的严重损害。 我们发现损伤引发了整个腰段新突触的异常生长，并伴有异常节段反射的出现。 SCI 下方回路的这种混乱重组与痉挛和神经元功能障碍有关

然后，我们收获了慢性 SCI 和未受伤对照小鼠的腰脊髓，并对这些组织进行了单核 RNA-seq (snRNA-seq)。我们总共对 19,237 个细胞进行了测序，这些细胞涵盖了腰脊髓的所有主要细胞类型 .

我们最初的目标是确定转录最受损伤干扰的细胞类型。 为了回答这个问题，我们使用 Augur 进行了细胞类型优先排序。 这种无偏见的分析表明，内皮细胞在损伤下方的脊髓中经历了最深刻的转录变化

这一意外发现促使我们研究这种优先排序背后的特定转录变化，以及不同统计方法揭示这些变化的能力。 为此，我们使用代表性的单细胞和假体方法对受伤和未受伤的内皮细胞进行了 DE 分析。 我们选择 Wilcoxon 秩和检验作为单细胞方法，因为该检验已成为单细胞转录组学领域中使用最广泛的方法，而 edgeR-LRT 因其高性能水平而作为pseudobulk方法。 这些方法确定了大部分不同的 DE 基因组，两种方法之间只有四个基因重叠。 相反，Wilcoxon 秩和检验和 edgeR-LRT 分别指定另外 44 个和 12 个基因作为 DE

迄今为止，我们的结果表明，未能解释重复之间的差异会导致单细胞 DE 方法产生错误的发现。因此，我们怀疑该数据集中通过 Wilcoxon 秩和检验确定的一些其他基因可能代表假阳性。为了阐明这些基因在受伤脊髓中的真实表达，我们进行了系统的体内筛选。我们获得了仅通过两种方法之一鉴定的 19 个假定 DE 基因的 RNAscope 探针，并量化了这些基因在受伤和未受伤小鼠的内皮细胞中的表达。 RNAscope 验证了被 edgeRLRT 称为 DE 的六个基因中的五个。与此形成鲜明对比的是，Wilcoxon 秩和检验称为 DE 的十三个基因中只有三个可以得到证实（p < 0.05，χ2 检验；图 5f-h）。一些经过验证的 edgeR-LRT 基因，包括 Slc7a11 和 Igfbp6，参与内皮细胞对缺氧的反应，支持在腰脊髓中建立慢性缺氧状态。与慢性缺氧的预期后果一致，我们检测到损伤水平以下存在大量萎缩性血管

img

总之，这些观察说明了单细胞DE方法产生错误发现的可能性。 相反，有效的单细胞DE分析，可以解释生物复制之间的差异，产生可重复的结论，可以在体内验证。

DE analysis with mixed models.

我们的实验表明，生物复制之间的差异决定了单细胞 DE 方法的性能。因此，我们对线性混合模型令人不满意的性能感到困惑。通过对生物复制内部和之间的变化进行明确建模，与pseudobulk方法相比，混合模型应该受益于更高的统计能力。为了澄清这种差异，我们评估了八个额外的泊松或负二项式广义线性混合模型 (GLMMs)。在 25-50 个细胞的数据集中，GLMMs 可以在非常特定的参数组合下产生准确的结果。然而，在包含 500 个或更多细胞的数据集中，它们的性能收敛于pseudobulk DE方法。此外，拟合最佳GLMMs所需的计算资源是巨大的。即使在下采样数据集中，单个细胞类型的DE分析平均需要13.5小时。相比之下，pseudobulk方法只需要在我们的 46 个数据集的概要中，每种细胞类型的分钟数。这些观察表明，在实践中，pseudobulk 方法为单细胞 DE 分析提供了速度和准确性之间的极好权衡。

img img img

Discussion

需要在单细胞转录组学中进行准确的 DE 分析，以剖析对疾病、创伤和实验操作的多方面反应背后的转录程序。 尽管统计方法对 DE 分析很重要，但决定其性能的原则仍然难以捉摸。 在这里，我们证明了有效 DE 分析的核心原则是统计方法能够解释生物复制的内在变异性。 考虑到这种可变性，决定了统计方法的生物学准确性。 相反，在没有任何生物学差异的情况下，未能解释生物复制变异性的方法可能会产生数百个错误发现。

研究人员研究单细胞以了解对生物扰动反应的更一般原理。澄清这些原则需要统计推断，这些推断超出构成任何特定数据集的单个细胞。我们的结果表明，通过在单个细胞水平上进行统计推断，单细胞 DE 方法将生物复制之间的变异性与生物扰动的影响混为一谈。重复之间存在差异是不可避免的，可以归因于技术因素和内在生物学差异。之前已经认识到，将重复之间的变异性与感兴趣的生物学效应混为一谈的可能性可能导致虚假的发现18,34。然而，这些研究几乎完全依赖于合成数据，并辅以一些说明性案例研究。因此，已发表的单细胞数据分析中错误发现的普遍性以及这些错误发现影响研究生物学结论的倾向仍不清楚。

在这里，我们表明错误发现的出现是一种普遍现象。利用具有实验基础事实的 18 个单细胞数据集的集合，我们证明使用不适当的统计方法会产生错误的发现，从而损害单细胞实验的生物学解释。这些错误的发现有可能浪费时间、精力和财力来追求误导性假设。例如，我们通过对受伤小鼠脊髓的系统体内筛选表明，通过最常用的统计方法鉴定的大多数 DE 基因都是错误的发现。此外，我们阐明了未能解释生物学和技术变异性使某些基因不成比例地被错误地识别为 DE 的机制的进展。我们在额外的 46 个单细胞 RNA-seq 研究的多方面数据集中证明了这些机制的普遍性。了解这些机制使我们发现，同样的基本问题会影响其他高维分析，包括空间转录组学，并且最有可能出现在人体组织的研究中，这表明生物复制水平的推断对于理解细胞和人类疾病的分子基础。

我们的结果表明，单细胞 DE 方法可能会产生错误的发现。这种理解揭示了该领域的巨大风险。我们的研究结果表明，许多已发表的研究结果可能是错误的。此外，如果不解决，可能会分配大量的研究资金来跟进这些错误的发现，从而损害科学。然而，使用解释重复之间变异性的 DE 方法可以直接纠正这种相关的可能性。其中，我们发现pseudobulk方法在转录组、蛋白质组和功能解释的水平上实现了对实验基本事实的最高保真度。因此，我们认为迫切需要对用于单细胞数据 DE 分析的统计方法进行范式转变。我们的观察强调了这种转变的必要性，即过去两年发表的大多数研究都使用了不适当的统计方法进行 DE 分析。此外，该领域使用最广泛的分析包目前默认采用易于错误发现的 DE 方法。多条件数据集的日益流行强调了采用适当的统计方法来防止错误发现扩散的重要性。为了促进这种转变，我们在 R 包中实现了这里测试的所有方法。

结论

pseudobulk找差异基因的方法优于单细胞直接找差异的方法，很多单细胞找差异基因的方法对高表达基因具有偏好性，主要是因为把每个细胞当成生物学复制的错误判断，从而造成了很多假阳性，错误的差异分析结果。 单细胞转录组学中的差异表达分析能够剖析细胞类型对疾病、创伤或实验操作等扰动的特异性反应。 虽然许多统计方法可用于识别差异表达的基因，但区分这些方法及其性能的原则仍不清楚。在这里，研究表明这些方法的相对性能取决于它们解释生物重复之间变异的能力。 忽略这种不可避免的变化的方法是有偏见的，并且容易出现错误的发现。 事实上，最广泛使用的方法可以在没有生物学差异的情况下发现数百个差异表达的基因（误差很大）。

作者github：neurorestore/Libra (github.com)（也是写出Augur包的实验室）

作者把这22种算法整合成了一个R包，可以直接使用

Libra is an R package to perform differential expression on single-cell data. Libra implements a total of 22 unique differential expression methods that can all be accessed from one function. These methods encompass traditional single-cell methods as well as methods accounting for biological replicate including pseudobulk and mixed model methods. The code for this package has been largely inspired by the Seurat and Muscat packages. Please see the documentation of these packages for further information.

System requirements

Libra relies on functions from the following R packages:

 dplyr (>= 0.8.0),
 purrr (>= 0.3.2),
 tibble (>= 2.1.3),
 magrittr (>= 1.5),
 tester (>= 0.1.7),
 Matrix (>= 1.2-14),
 pbmcapply (>= 1.5.0),
 lmtest (>= 0.9-37),
 tidyselect (>= 0.2.5),
 DESeq2 (>= 0.4.0),
 Seurat (>= 3.1.5),
 blme (>= 1.0-4),
 edgeR (>= 3.28.1),
 glmmTMB (>= 1.0.2.1),
 limma (>= 3.1-3),
 lme4 (>= 1.1-25),
 lmerTest (>= 3.1-3),
 matrixStats (>= 0.57.0),
 methods,
 stats,
 Rdpack (>= 0.7)

In addition, the Seurat, monocle3, or SingleCellExperiment packages must be installed for Libra to take Seurat, monocle, SingleCellExperiment objects as input, respectively. Methods that require additional packages may also require additional installs (e.g., MAST).

Libra has been tested with R version 3.6.0 and higher.

Installation

To install Libra, first install the devtools package, if it is not already installed:

install.packages("devtools")

Libra additionally requires the following installations to perform differential expression testing:

if (!requireNamespace("BiocManager", quietly = TRUE))
 install.packages("BiocManager")

 BiocManager::install(c("edgeR", "DESeq2", "limma"))

If Seurat is not installed this will be needed for single-cell methods.

install.packages("Seurat")

Finally, install Libra from GitHub:

devtools::install_github("neurorestore/Libra")

This should take no more than a few minutes.

Usage

The main function of Libra, run_de, takes as input a preprocessed features-by-cells (e.g., genes-by-cells for scRNA-seq) matrix, and a data frame containing metadata associated with each cell, minimally including the cell type annotations, replicates, and sample labels to be predicted. This means that in order to use Libra, you should have pre-processed your data (e.g., by read alignment and cell type assignment for scRNA-seq) across all experimental conditions.

Libra provides a universal interface to perform differential expression using 22 discrete methods. These methods are summarized as follows:

Single cell methods

• Wilcoxon Rank-Sum test

• Likelihood ratio test

• Student's t-test

• Negative binomial linear model

• Logistic regression

• MAST

Pseudobulk methods

• edgeR-LRT

• edgeR-QLF

• DESeq2-LRT

• DESeq2-Wald

• limma-trend

• limma-voom

Mixed model methods

• Linear mixed model

• Linear mixed model-LRT

• Negative binomial generalized linear mixed model

• Negative binomial generalized linear mixed model-LRT

• Negative binomial generalized linear mixed model with offset

• Negative binomial generalized linear mixed model with offset-LRT

• Poisson generalized linear mixed model

• Poisson generalized linear mixed model-LRT

• Poisson generalized linear mixed model with offset

• Poisson generalized linear mixed model with offset-LRT

By default Libra will use a pseudobulk approach, implementing the edgeR package with a likelihood ratio test (LRT) null hypothesis testing framework. Each of the 22 tests can be accessed through three key variables of the run_de function: de_family, de_method, and de_type. Their precise access arguments are summarized in the below table.

Method	de_family	de_method	de_type
Wilcoxon Rank-Sum test	singlecell	wilcox
Likelihood ratio test	singecell	bimod
Student's t-test	singlecell	t
Negative binomial linear model	singlecell	negbinom
Logistic regression	singlecell	LR
MAST	singlecell	MAST
edgeR-LRT	pseudobulk	edgeR	LRT
edgeR-QLF	pseudobulk	edgeR	QLF
DESeq2-LRT	pseudobulk	DESeq2	LRT
DESeq2-Wald	pseudobulk	DESeq2	Wald
limma-trend	pseudobulk	limma	trend
limma-voom	pseudobulk	limma	voom
Linear mixed model	mixedmodel	linear	Wald
Linear mixed model-LRT	mixedmodel	linear	LRT
Negative binomial generalized linear mixed model	mixedmodel	negbinom	Wald
Negative binomial generalized linear mixed model-LRT	mixedmodel	negbinom	LRT
Negative binomial generalized linear mixed model with offset	mixedmodel	negbinom_offset	Wald
Negative binomial generalized linear mixed model with offset-LRT	mixedmodel	negbinom_offset	LRT
Poisson generalized linear mixed model	mixedmodel	poisson	Wald
Poisson generalized linear mixed model-LRT	mixedmodel	poisson	LRT
Poisson generalized linear mixed model with offset	mixedmodel	poisson_offset	Wald
Poisson generalized linear mixed model with offset-LRT	mixedmodel	poisson_offset	LRT

If batch effects are present in the data, these should be accounted for, e.g., using Seurat or Harmony, to avoid biasing differential expression by technical differences or batch effects.

To run Libra with default parameters on a genes-by-cells scRNA-seq matrix expr, and an accompanying data frame meta, with cell_type, replicate, and label columns containing cell types, replicates, and experimental conditions, respectively, use the run_de function:

library(Libra)
DE = run_de(expr, meta = meta)

If your columns have different names, you can specify these using the cell_type_col, replicate_col, and label_col arguments:

DE = run_de(expr, meta = meta, cell_type_col = "cell.type", label_col = "condition")

If you would like to store the pseudobulk matrices in a variable, before running differential expression, you can do the following:

matrices = to_pseudobulk(expr, meta = meta)

Libra can also run directly on a Seurat object. For a Seurat object sc, with the sc@meta.data data frame containing cell_type and label columns, simply do:

DE = run_de(sc)

The same code can be used if sc is a monocle3 or SingleCellExperiment object instead.

Demonstration

To see Libra in action, load the toy single-cell RNA-seq dataset that is bundled with the Libra package:

data("hagai_toy")

This dataset consists of 600 cells, distributed evenly between six replicates and two conditions.

The hagai_toy object is a Seurat object with columns named cell_type and label in the meta.data slot, meaning we can provide it directly as input to Libra:

head(hagai_toy@meta.data)

 orig.ident nCount_RNA nFeature_RNA                                       sample replicate label                                  cell_type
2-CGAACATGTATAATGG SeuratProject         18           11 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes
2-ATTCTACAGTGGTAGC SeuratProject         23           12 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes
2-GCATACACAAACTGTC SeuratProject          3            3 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes
2-GAAGCAGAGATGCCAG SeuratProject         14            8 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes
2-ATCACGAGTCTAGCGC SeuratProject          3            3 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes
2-CTGCCTATCTGTCCGT SeuratProject         28           13 mouse1_lps4_filtered_by_cell_cluster0.txt.gz    mouse1  lps4 bone marrow derived mononuclear phagocytes

We run run_de, and inspect the differential expression:

DE = run_de(hagai_toy)
head(DE)

We can also use a different statistical framework, for example DESeq2:

DE = run_de(hagai_toy, de_family = 'pseudobulk', de_method = 'DESeq2', de_type = 'LRT')
head(DE)

Alternatively, we can use a mixed-model approach, which by default will use a negative binomial model structure:

DE = run_de(hagai_toy, de_family = 'mixedmodel')
head(DE)

However, this can be adjusted using the de_method argument:

DE = run_de(hagai_toy, de_family = 'mixedmodel', de_method = 'linear', de_type = 'LRT')
head(DE)

Running this example on a MacBook should be instantaneous. However, analyzing >20 real single-cell RNA-seq datasets, we found Libra takes a median of ~5 minutes. In general, runtime scales close to linearly with the number of cell types and cells. If using mixed models, by default, Libra runs on four cores, with each gene analyzed on a different core. To change the number of cores, use the n_threads argument. For example, running Libra on eight threads:

DE = run_de(hagai_toy, de_family = 'mixedmodel', de_method = 'linear', de_type = 'LRT', n_threads = 8)

... will run about twice as fast.

Calculating delta variance

We recently showed that statistical methods for differential expression must account for the intrinsic variability of biological replicates to generate biologically accurate results in single-cell data (Squair et al., 2021, Biorxiv; https://www.biorxiv.org/content/10.1101/2021.03.12.435024v1). Within the same experimental condition, replicates exhibit inherent differences in gene expression, which reflect both biological and technical factors. We reasoned that failing to account for these differences could lead methods to misattribute the inherent variability between replicates to the effect of the perturbation. To study this possibility, we compare the variance in the expression of each gene in pseudobulks and pseudo-replicates. We call this measure 'delta variance'. Users can use the calculation of delta variance to inform their differential expression results. For example, genes identified as differentially expressed by methods that do not account for biological replicate (i.e., 'singlecell' methods) that have a high delta variance should be treated with caution as they are likely to be false positives. Delta variance can be calculated as follows:

DV = calculate_delta_variance(hagai_toy)

This function will return a list of vectors, one for each cell type, each of which contains the delta variance for the genes present in the input expression matrix.

有关转换成pseudobulk的知识可以看这个链接：Single-cell RNA-seq: Pseudobulk differential expression analysis | Introduction to Single-cell RNA-seq - ARCHIVED