前言

差异表达分析是单细胞RNA测序（scRNA-seq）数据分析的一个关键步骤，是后续很多高级分析的基石。在早期的研究中，差异表达分析主要用于识别不同亚群之间的表达差异，以确定各亚群的标记。然而，随着多样本、多组别的scRNA-seq数据集越来越多。因此，传统的单细胞差异分析工具可能不太适合检出此情况中更有价值的差异基因，而样本级别的差异分析可能更适合这样的情形。

我们今天要介绍的主角：muscat，就是适用于对多样本、多组别的scRNA-seq数据进行差异分析工具，其更加关注如何解释样本间和细胞间的变异性，并从中得出可以推广到样本层面的结论，即差异状态（differential state，DS）分析。为了探索目前什么样的统计方法是最适合这类分析的，这项研究模拟多样本多状态的scRNA-seq数据集，并对16种不同的DS分析方法进行了评估。这些方法包括基于细胞级别的混合模型和基于聚合伪批量数据的方法。

那么接下来，让我们先通过对应的文献一起学习一下这项研究（ps：后续还会有代码实操部分对它的功能进行介绍）。

主要内容

一、问题起源

在学习这项研究的具体结论之前，我们先搞清楚以下两个问题：什么是差异状态（DS）分析？以及为什么基于细胞级别的混合模型和基于聚合伪批量数据的方法可以检测细胞的状态变化？

1. 什么是差异状态分析？

DS分析的三个关键方面：

• 能够检测特定亚群的变化。

• 与聚类和亚群分配的分析正交。即DS分析与聚类或细胞亚群分配分析之间具有独立性和互补性。独立性：DS分析的结果不依赖于聚类或亚群分配分析的结果；互补性：DS分析的结果与聚类或亚群分配分析的结果相互补充。与聚类或亚群分配不同，DS分析关注的是细胞状态的动态变化，而不是静态的细胞类型，DS分析可以揭示这些细胞在不同条件下的变化。

• 差异状态（DS）分析与跨条件下亚群的差异丰度（DA）分析是两种不同的分析。DS分析更关注细胞的功能状态和行为的变化而DA分析更关注细胞亚群的数量或丰度的变化。

2.为什么基于细胞级别的混合模型和基于聚合伪批量数据的方法可以检测细胞的状态变化？

相比于bulk RNA-seq数据，scRNA-seq数据包含着更为丰富的信息，每个细胞的基因表达模式往往可以反应出细胞的状态（静息或活化），这在bulk RNA-seq数据中常常被掩盖。基于细胞级别的混合模型和基于聚合伪批量数据的方法可以帮助我们挖掘类似结论。

• 基于细胞级别的混合模型推断细胞状态变换

  我们首先对scRNA-seq数据进行降维聚类分群，对于一个我们感兴趣的特定亚群，我们可以使用混合模型对每个细胞的基因表达建模。具体来说，我们可以假设静息状态的细胞的基因表达服从一个分布，而活化状态的细胞的基因表达服从一个不同的分布，从而推断每个细胞处于活化状态的概率。值得注意的是，相比于伪批量的方法，基于细胞级别的混合模型通常需要更复杂的算法，并且可能需要更多的计算资源。

• 基于聚合伪批量数据的方法推断细胞状态变换

  同理，我们首先对scRNA-seq数据进行降维聚类分群并聚合为伪批量数据，也就是将来自同一亚群和同一条件的所有细胞的基因表达数据进行合并，对于一个我们感兴趣的特定亚群， 这个亚群在对照组中的细胞大多数是静息状态，而在疾病组中的细胞大多数是活化状态。在这种情况下，我们可能会观察到，这个亚群在疾病组的伪批量数据中，与细胞活化相关的基因的表达水平比在对照组中的伪批量数据中要高。

  那么，在这个例子中，我们为什么不直接在单细胞层面比较疾病组和对照组与细胞活化相关的基因表达的水平差异呢？主要原因可能有两个：1、直接在细胞层面进行比较会受到一些挑战。由于单细胞数据固有的噪音、drop out的存在或样本中细胞状态存在较大的异质性等原因阻碍着我们对细胞状态转变的识别。2、伪批量处理还可以使我们更容易地使用一些针对bulk数据的分析工具，如edgeR或DESeq2等。

二、一个灵活可靠的模拟模型框架

在本研究中，作者基于参考数据集开发了一个模拟框架，模拟了scRNA-seq数据的各种特性，并用它评估了16种DS分析方法，涵盖了广泛的模拟情景，包括样本数量的变化、每个亚群的细胞数量的变化，以及引入的差异表达模式的大小和类型。首先，作者使用负二项分布（NB）模拟scRNA-seq数据集，参考数据集有两套：（i）来自8名红斑狼疮患者的外周血单个核细胞（PBMCs）的scRNA-seq数据，这些数据是在IFN-β治疗前后6小时测量的（总共16个样本），其中的细胞已经被标注为各种免疫亚群；以及（ii）来自8只被分为载体和LPS治疗组的小鼠的大脑皮层组织的单核RNA-seq数据。接下来，作者通过增加平均表达的变化（DE）、低和高表达状态组分的比例变化（DP）、差异模态（DM）或比例和模态的变化（DB）模拟细胞状态的变化。不受状态变化影响的基因要么等效表达（EE），要么在两种条件下由相等比例的细胞以低和高表达状态表达（EP）。通过调整参数（亚群特异基因比例、DS基因比例、logFC值等）可以控制模拟数据集中的亚群特异性和状态变化。总的来说，作者构建了一个灵活可靠的模拟单细胞RNA测序数据的框架，为后续测试16种不同的DS分析方法提供基础。

image.png

三、评估16种DS分析方法的性能

对于基于聚合的方法，作者考虑了各种输入数据（对数转换的表达值，残差，计数）、汇总统计（均值，总和）和差异测试方法（limma-voom，limma-trend，edgeR）的组合。对于非聚合的方法，作者考虑了四种方法：混合模型MM、MAST，scDD、AD。MAST用于对数化的count矩阵；AD测试和scDD用于对数化的count矩阵和标准化残差（vstresiduals）。对于AD测试，作者考虑了两种不同的方法来测试分布，替代假设样本在样本方面或组方面有所不同。接下来，作者评估了不同方法在检测单细胞RNA测序数据中的差异状态方面的性能。总体而言，所有方法对DE类别的基因表现最好，其次是DM、DP和DB。对于DE、DM和DP，大多数伪批量方法和细胞水平MM模型表现良好。为了研究亚群大小对DS检测的影响，作者使用每个亚群样本20-400个细胞的子集，在包含10% DE基因的模拟上运行了方法。对于大多数方法，FDR控制随细胞数量的变化而剧烈变化，而所有方法的TPR随着细胞数量的增加而改善。对于基于聚合的方法，大约100个细胞足以达到不错的性能；特别是从20个细胞增加到100个细胞（每个亚群每个样本）会有相当大的性能提升，但从200增加到400个细胞只有适度的增益。

image.png

为了研究整体方法一致性，作者在每个DS类别的五个模拟重复实验中，将每种方法返回的排名最高的DS检测（FDR < 0.05）相交。结果表明，方法之间的整体一致性很高，共同识别的DS基因确实存在真实差异，而被某个方法单独识别的DS基因更有可能为假阳性发现。在运行时间上， MMs是最慢的，其次是AD测试、MAST，然后是scDD。基于聚合的DS方法是最快的。

image.png

四、LPS处理小鼠皮质的DS分析

为了探究LPS如何影响大脑皮质，作者使用snRNA-seq研究外周LPS给药对小鼠前额叶皮质中所有主要细胞类型的影响，识别在神经细胞和非神经细胞中受LPS影响的基因和通路。在这里，作者使用伪批量和edgeR将DS分析框架应用于四个对照（载体）和四个LPS处理的小鼠的snRNA-seq数据。首先，作者将snRNA-seq注释为8个亚群：星形胶质细胞，内皮细胞，微胶质细胞，髓鞘形成细胞前体细胞（OPC），脉络丛室管膜（CPE）细胞，髓鞘形成细胞，兴奋性神经元和抑制性神经元。作者在至少一个亚群中鉴定出915个具有DSs（FDR < 0.05，∣logFC∣ > 1）的基因，其中751个只在一个亚群中检测到（补充图13）。由于仅依赖阈值容易产生偏见，接下来作者将所有DEGs并集的FC值进行聚类。结果表明，一组明显的基因集（共识聚类ID 3）在所有亚群中都被上调，并且富集响应（外部）生物刺激、防御和免疫反应相关功能。接下来，作者计算了效应系数，总结了每个细胞在多大程度上反映了整体的倍数变化。对于内皮细胞和胶质细胞，效应系数分布在载体和LPS样本之间明显分离，表明大多数细胞受到影响。相反，神经元分布的大量重叠表明只有少数细胞受到影响。

image.png

五、DS方法性能总结

最后，作者对这16种DS检测的敏感性和特异性，p值分布的均匀性，模拟和估计的logFCs之间的一致性，适应复杂实验设计的能力，以及运行时间等方面进行了评估。

image.png

小结

目前，我们对scRNA-seq数据进行差异分析时，更多关注的是细胞水平的差异，而较少关注样本水平，这导致我们得到的差异结果主要是围绕DA（差异丰度）的相关结果。这项研究的作者引入了DS（差异状态）的概念，即关注细胞的功能状态或行为的变化；并且基于这个概念，开发了一个灵活的DS模拟框架。这个框架评估了16种DS分析方法，使研究人员可以使用这个框架识别scRNA-seq中更可靠的DS结论。为了让大家可以更全面地了解到这个工具，我们将在后续的推文中介绍muscat的代码实操部分。

今天的分享就到这里啦，让我们下期再会~

[参考文献]

Crowell HL, Soneson C, Germain PL, Calini D, Collin L, Raposo C, Malhotra D, Robinson MD. muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data. Nat Commun. 2020 Nov 30;11(1):6077. doi: 10.1038/s41467-020-19894-4. PMID: 33257685; PMCID: PMC7705760.