摘要

动机：

单细胞亚硫酸盐测序（scBS-seq）方法的使用，允许对DNA甲基化模式进行精确的单细胞级别分析，能够识别稀有的种群，揭示细胞特异的表观遗传变化，并提高差异甲基化分析的精确度。然而，由于限制的测序深度和覆盖范围，数据稀疏和零和一的过多，使得使用scBS-seq进行差异甲基化检测时，精确度准确率经常降低。因此，急需一种创新的差异甲基化分析方法，有效处理这些数据特性并提高识别准确率。

结果：

结果：我们提出了一种名为scDMV的新型贝塔混合方法，用于分析单细胞亚硫酸盐测序数据中的甲基化差异，该方法有效地处理了过多的零和一，并能适应低输入测序。我们的大量模拟研究表明，scDMV方法在敏感性、准确性和控制假阳性率方面都优于其他几种方法。此外，在实际数据应用中，我们观察到即使在低输入样本中，scDMV在识别差异甲基化区域方面也表现出更高的精度和敏感性。另外，scDMV揭示了使用单细胞全基因组测序数据进行GO富集分析时，其他方法常常忽视的重要信息。

介绍

表观遗传学研究的是与DNA序列变化无关的可遗传的基因表达变化。关键的表观遗传修饰，如DNA甲基化、组蛋白修饰、启动子-增强子相互作用和非编码RNA调控，都起着至关重要的作用，可能导致疾病的发生。在这些修饰中，由于DNA甲基化的可逆性和其作为药物靶点的潜力，它已经受到了大量的关注。在哺乳动物中，DNA甲基化主要发生在CpG位点，那里的胞嘧啶的第五个碳原子被DNA甲基转移酶甲基化，形成5-甲基胞嘧啶。CpG位点可以分布在整个DNA序列中，也可以集中在位于调控区域的CpG岛上。理解DNA甲基化对阐明其对细胞发育、疾病进展和基因调控的影响至关重要。

分析样本间的DNA甲基化差异对于理解疾病发病机制、预防疾病和诊断疾病至关重要。常用的两种甲基化差异分析方法分别是差异甲基化位点（DMS）分析和差异甲基化区域（DMR）分析。DMS分析侧重于单一样本内的个别甲基化位点，而与基因表达的关联较小。相比之下，DMR分析考虑的是由一个或多个DMS组成的连续区域，并允许在多个样本组之间进行比较，从而提供对基因表达更多的理解。

近年来，用于识别差异甲基化的基于测序的方法有所增加。这些方法包括多种方法，如逻辑回归、贝塔-二项分布（beta-binomial distribution）、隐马尔可夫模型、香农熵和二元分割平滑（binary segmentation）。现有的算法包括“eDMR”，“RADMeth”，“BSmooth”和"CGmapTools"。

当依赖来自多个细胞的平均数据时，使用传统策略研究DNA甲基化多样性的能力是有限的。单细胞全基因组亚硫酸盐测序（Single-cell whole-genome bisulfite sequencing）(scWGBS 和 scRRBS)已经成为评估单个细胞和稀有细胞类型的DNA甲基化多样性的有前景的方法。然而，单细胞DNA甲基化测序数据的稀疏性和独特特性（包括低覆盖率和过多的零和一（例如在2.4节的真实例子中，0和1的甲基化率之和超过0.9)），使得传统的统计方法不充分。因此，需要新的方法来进行使用scBS-seq数据的差异甲基化分析。

在这篇文章中，我们提出了一种名为scDMV（(zero–one inflated beta mixture model）的策略，用于分析单细胞亚硫酸盐数据。我们假定，根据表示甲基化率的细胞和区域特异效应的条件，scBS-seq数据遵循二项分布。此外，我们使用零一通胀的贝塔分布来模拟效应分布，以解释零和一的过度存在，以及在scBS-seq数据中观察到的过度离散化。我们采用EM算法来估计模型参数，并利用Wald检验来进行差异甲基化分析。我们通过包括模拟实验和真实数据应用在内的数值研究，将scDMV的性能与两种现有的方法，methylpy和CGmapTools进行了比较。结果表明，scDMV的性能优越，特别是在捕获单细胞全基因组测序数据进行GO富集分析的重要信息方面。

2.材料与方法

观察到的scBS-seq数据由一个集合表示，记为

$n_{gij}$ 表示从第j个区域的g类型的第i个细胞获得的总读取， $x_{gij}$ 表示从第j个区域的g类型的第i个细胞获得的甲基化读取，g对应于细胞类型； $N_{g}$ 表示属于g类型的细胞数，i表示样本或细胞，j对应于不同的CG区域，M表示考虑的CG段的总数。

设 $p_{gj}$ 表示在第j区域属于g类型的细胞的甲基化率。我们将Pg定义为g类型细胞的甲基化率向量，由 $P_{g} = (p_{g1},...,p_{gM})^T$ 给出。差异甲基化分析(differential methylation analysis)的主要目标是检验两个甲基化率向量之间的平均甲基化水平是否相等的零假设( null hypothesis)是否成立。替代假设(alternative hypothesis)，记为H1，暗示存在特定区域，记为 $m \in (1,...,M )$ ，其中两种细胞类型的平均甲基化率不同。

为了处理这个假设检验问题，我们首先构造一个检验统计量。随后，我们开发了一种识别DMR集的程序。

2.1制定模型和测试统计量

假定总读数 $n_{gij}$ ，有理由假设甲基化读数 $x_{gij}$ 遵循二项分布，可以表示为：

其中 $p_{gj}$ 代表区域j中g类型细胞的甲基化率。值得注意的是，细胞的甲基化率显示出显著的异质性( heterogeneity)，通常以过多的零和一为特征。为了捕捉这种变异性，我们将 $p_{gj}$ 建模为随机效应( random effect)，并定义其均值为 $E_{p_{gj}} = \mu_{gj}$ 。对于每个区域j，我们测试假设

观察到单个细胞中的DNA甲基化率pgj倾向于在0和1之间集中，我们假设pgj遵循0-1通胀的混合贝塔分布。这个分布可以如下特征化：

在这个模型中，评估两种样本类型之间的差异甲基化表达归结为检查：

重要的是要强调，与两种样本类型相关的四个参数可以在CG区域j中变化(vary across CG)。现在，让我们继续推导出估计这些参数的算法。对于给定的g和j值，我们计算似然函数，并使用EM算法估计参数 $\theta_g = (\pi_{gj0},\pi_{gj1},\alpha_{gj},\beta_{gj} ).$ 。信息矩阵 $I(\pi_{gj0},\pi_{gj1},\alpha_{gj},\beta_{gj})$ 的形式为

以及g类细胞在区域j中的甲基化率的期望值

在零假设下，得到的Wald检验统计量为

渐近的遵循标准正态分布。总的来说，基于等式（1）导出了检验统计量，图1中的流程图描绘了scDMV方法。

2.2 DMR 识别

接下来，我们评估所有M个区域以得到M个P值。P值低于指定截止值的区域被认为显著不同，并被分类为DMR。随后，常用的方法利用样本间的DNA甲基化差异精化最初获得的DMR，以提高准确性。在我们的研究中，使用加权DNA甲基化来计算甲基化程度差异，记为D，该差异在每个区域内的两种样本之间存在。然后为D设置一个截止值，最后的DMR必须满足D超过截止值的条件。通常，模型识别的DMR是那些P值低于P值截止值和D超过D截止值的区域。

2.3 仿真研究

现有的方法，如CGmapTools和methylpy，是专门设计用于分析由大量细胞形成的细胞群之间的差异的，这些细胞使用的是传统的亚硫酸盐测序技术。因此，它们可能不能直接应用于单细胞数据。相比之下，我们在这项研究中提出的新的统计方法，称为scDMV，考虑到了scBS-seq测序数据中深度和覆盖面低的特定特性。为了评估和比较CGmapTools、methylpy和新的scDMV方法在识别单细胞DNA甲基化数据中的DMR方面的性能，我们基于模拟的scBS-seq数据进行了几次模拟实验。

我们生成了一个由73个样本组成的模拟数据集，其中两种明显的细胞类型的样本大小分别为48和25。每个样本包括10000个位点，其中每个位点的甲基化读取由x表示，总读取由n表示。数据集分为1000个区域，每10个连续位点形成一个区域。数据模拟过程涉及三个主要步骤。首先，我们为两种样本类型分配四个参数的值（p10; p11; a1; b1）和（p20; p21; a2; b2）。其次，我们根据实际的总阅读数获得模拟的总阅读数。最后，根据基础理论模型生成模拟的甲基化阅读数。

为了说明生成模拟数据的过程，我们以一个特定的区域为例。首先，我们通过从真实总读取数据的每一列中随机抽取10个非零值，生成模拟的总读取数据，记为n。然后，根据前面提到的先验分布生成该区域的甲基化读取x。按照这个程序，我们获得了特定区域的模拟数据。

上述过程重复1000次，产生1000个区域的模拟数据。每个区域的甲基化读取x和总读取n数据分别存储在单独的列表中，结果得到1000个列表的集合。然后这1000个列表被合并成一个综合列表，代表最后的模拟数据。模拟数据中0和1处的甲基化率之和是0.66。

进行了两种类型的模拟实验：差异实验（difference experiments）和无差异实验（indifference experiments）。在差异实验中，参数值（p10; p11; a1; b1）和（p20; p21; a2; b2）特意设置为不同，产生了模拟数据，展现了组间的差异。另一方面，在无差异实验中，（p10; p11; a1; b1）和（p20; p21; a2; b2）被设定为完全相同，产生了没有组差异的模拟数据。对于每种类型的实验，都模拟了五组实验数据。

为了评估scDMV方法在识别DMR方面的比较准确性，我们进行了模拟实验，并将其性能与methylpy和CGmapTools进行了比较。样本组的区域平均甲基化水平被定义为甲基化读取总和与该区域内所有样本中所有位点的总读取总和的比率。甲基化水平的差异，记为D，表示了同一区域内两个样本组之间的区域平均甲基化水平的差异。对于每个方法，分别为每个区域计算P值和D值。我们采用了不同的P值（0.001，0.005，0.01和0.05）和D（0，0.1，0.15和0.2）截止值。满足P值不超过指定截止值和D超过定义阈值的区域被认为是识别出的DMR。

2.4 真实数据应用

为了探索早期胚胎发育期间的甲基化模式，我们在一个公开可用的数据集（Benjamini和Hochberg，1995）上使用了scDMV方法，并将其结果与两种替代方法进行了比较。

2.4.1 真实数据实验设计

这个数据集（GEO ID: GSE81233）包括来自两个连续发育阶段的73个样本，包括25个4-细胞样本和48个8-细胞样本。进行了组内和组间实验。在组间实验中，涉及到有差异的样本，我们使用所有的73个样本来识别25个4-细胞胚胎样本和48个8-细胞胚胎样本之间的DMR。通过应用不同的P值截止值和不同的DNA甲基化水平差异（D）阈值，获得了三种方法的实验结果。在组内实验中，样本之间没有差异，我们从48个8-细胞胚胎样本中选择了40个，将它们平均分为两组来识别DMR。由于同一发育阶段的样本的整体甲基化模式倾向于保持稳定（Benjamini和Hochberg，1995），所以这两组之间不应存在甲基化的个体差异。因此，在这种情况下识别出的DMR可以被认为是假阳性DMR。

为了提高识别DMR的准确性，我们进行了数据预处理和过滤程序。首先，我们评估每个位点的显著性，以确定其对区域的影响。被排除的位点包括丢失值超过总样本数50%的位点。位点过滤后，将数据划分为最大长度为300bp的区域，确保每个区域包含至少3个位点。按照上述数据组织方法，每个区域由两个列表表示：一个用于总阅读数据，另一个用于甲基化阅读数据。每条染色体的所有区域都存储为一个名为“testRegion”的集合列表，从而为每条染色体形成一个列表。随后，使用“testRegion”中的数据进行使用scDMV方法的实验。在实验过程中，我们采用加权方法计算区域甲基化水平。我们对每条染色体执行实验，为每个区域生成P值和D值。最后，我们应用之前定义的截止值来相应地过滤区域。

在最后一步，我们将两个实验的结果合并，进行全面的分析和比较三种方法，从而评估scDMV方法在识别DMR方面的性能。

2.4.2 DMR的注释

通过ChIPSeeker R包（版本1.24.0）（Yu et al. 2015），基于它们在人类基因组（hg19）中对应的区域，对8细胞阶段和4细胞阶段之间确定的DMR进行了注释。注释过程包括基于它们相对于基因转录的位置对DMR进行分类，包括启动子区域（距离转录起始位点或TSS 2kb范围内）、内含子、外显子以及基因间区域。此外，DMR还基于它们与CpG岛、CpG岸（距离岛2kb范围内）、CpG搁浅区（距离岸2kb范围内）以及开放海域（在之前三个区域之外）的关联进行注释。CpG岛的注释信息是从UCSC Genome Browser网站（http://genome.ucsc.edu）（Karemaker和Vermeulen，2018）获得的。

2.4.3功能富集分析

使用Metascape软件（http://metascape.org）（Karolchik等人，2003）识别富集的基因本体（GO）术语。功能富集分析的基因列表包括在它们的启动子（距离TSS 2kb范围内）和/或基因体内含有DMR的基因。对于GO富集分析，只选择与“生物过程”相关的术语。使用Benjamini-Hochberg方法调整P值，以控制假发现率（FDR）（Zhou et al.2019）。

在差异实验中，所有1000个模拟数据区域都被划定为差异甲基化区域（DMRs），而在无差异实验中，所有1000个区域表现出无差异甲基化。为了确保模拟数据的有效性，我们将其分布与真实数据（具体来说，是单细胞甲基化测序数据中染色体1的数据）进行了比较，如补充附录中的补充图表S1所示。图表呈现了以下观察结果：(i) 模拟数据和真实数据的范围都被扩大到0-1；(ii) 模拟数据的分布与真实数据的分布紧密对应。

我们通过分析五个模拟实验的平均结果，来评估scDMV和对比算法（Methylpy, CGmapTools）的集合性能。P值截止值为0.01的结果在表1中展示，而所有五个实验的详细结果可以在补充附录的补充表格S1-S3中找到。

在模拟实验中，我们使用敏感性和精确性作为算法的性能指标。敏感性是以P值不超过定义的截止值的区域数量与在差异实验中总区域数量之比计算的。另一方面，精确性是以在差异实验中确认的DMRs数量与所有实验中识别的DMRs总数之比确定的。

为了可视化结果，我们在各个截止点上绘制了每种方法的五个实验结果，以假阳性发现率(FDR)表示在横坐标轴，以敏感性表示在纵坐标轴。补充附录的补充图表S2展示了这幅图，其中黑色垂直线对应FDR为0.005，表示当FDR低于此阈值时，统计显著性。

—————————————————

通过检查图表，可以明显看出，scDMV算法在控制I型错误时始终表现出更高的敏感性，尤其是当D（甲基化水平的差异）为0或更大时。此外，scDMV在比较其他两种算法时展示出更优的精确性，对于D的多个截断点（特别是0、0.1和0.15）都能保持高敏感性。

此外，我们还展示了三种方法的精确性盒图，如补充附录图表S3所示。如补充图S3A所示，无论筛查条件如何，scDMV始终能达到比其他两种方法更高的精确性。值得注意的是，在所有情景中，scDMV的精确性始终高于0.98。

除了敏感性和精确性，研究人员通常也对算法的假阳性率（FPR）进行评估。补充图表S3B展示了三种方法的FPR盒图，明显显示出scDMV始终能维持低于其他两种方法的假阳性率。值得注意的是，methylpy的假阳性率显著高，这可能由该方法的实验原理造成。

基于模拟结果，我们可以得出以下结论：scDMV在识别单细胞数据中的DMRs时，表现出了非凡的精确性和敏感性。总的来说，前述的模拟结果表明，scDMV方法在精准检测单细胞亚硫酸盐测序数据中的DNA甲基化差异方面，超过了其他两种方法。

————————

3.2

在评估scDMV的性能时，我们使用精确性作为评价标准。精确性被定义为在所有识别出的DMRs中正确检测到的DMRs的比例。在前述的实验中，我们假设在组间实验中检测到的DMRs都是正确识别的，而在组内实验中检测到的DMRs都是误识别的。因此，我们将组间实验中识别出的DMRs数量定义为真阳性（TP0），将组内实验中发现的DMRs数量定义为假阳性（FP0）。因此，精确性可以使用以下公式进行计算：

数据处理后，第一次实验产生了总共11083个区域，而第二次实验得到了13193个区域。我们设置P值截断点为0.001、0.005、0.01和0.05，以及D截断点为0.1、0.15和0.2。

因此，我们得到了在不同阈值下的三种方法的实验结果（完整的实验结果参见补充附录的补充表S4）。

结果表明，随着阈值变得更严格，CGmapTools和scDMV方法的精确性逐渐降低，而methylpy方法的精确性相对稳定。具体来说，当P值的截断点设置为0.001、0.005和0.01时，scDMV方法始终保持了0.71以上的精确性水平，而其他两种方法未能达到0.66的精确性。当我们将P值截断点设置为0.05时，scDMV的精确性下降，但仍高于0.65，而两种相应方法的精确性低于0.59。

在比较不同阈值的结果时，scDMV相比其他两种方法始终表现出较高的精确性。换句话说，scDMV在确保误识别区域减少的同时，可以识别出更多的DMR。

3.3 差异甲基化的表征

为了准确描述8-细胞阶段和4-细胞阶段之间的DMRs，我们设定了严格的阈值选择基于P值≤0.001和D≥0.2的DMRs。结果发现，全基因组共有1457个DMRs。在这些DMRs中，有1446个（99.25%）在8-细胞胚胎中显示出甲基化过度（图2A）。这个结果与先前的研究一致，后者报告了8-细胞胚胎与4-细胞胚胎相比全局DNA甲基化水平的显著上升（Benjamini和Hochberg 1995，Zhu等。2018）。

为了研究DMRs的基因组分布，我们使用了clusterProfiler包（Hanna等人。2016）基于人类hg19参考基因组注释区域。分析显示，连续发育阶段之间的大部分DMRs（49.69%）位于转录本的内含子区域中（图2B）。此外，11.94%的DMRs被发现在启动子区域，这与启动子甲基化与转录沉默的关联是一致的（Siegfried等人。1999，Yu等人。2012）。值得注意的是，当比较scDMV方法与其他两种方法时，也观察到了DMRs基因组分布的相似模式（图2D）。

此外，一部分明显的比例（4.95%）由scDMV工具检测到的DMRs被发现位于CpG岛内（图2C），这与由CGmapTools工具识别的DMRs的基因组分布相符（图2E）。相反，methylpy识别的DMRs对CpG密度高的区域有更大偏好，如CpG岛（9.02%）和CpG海滨（14.7%）（图2E）。methylpy在DMR识别过程中的这种偏见，包括首先识别差异甲基化位点然后将其合并为DMRs，可能有助于这种观察。可能像scDMV这样的方法，直接定义候选甲基化区域，通过避免这些偏见达到更高的准确性（Baylin 2005）。

3.4 scDMV有效地捕获了CGmapTools忽视的关键信息

由于methylpy方法的假阳性率较高，我们将比较重点放在scDMV和CGmapTools方法之间。其中，scDMV报告了总共1457个DMRs，而CGmapTools报告了535个DMRs。值得注意的是，scDMV捕获了CGmapTools识别的95.7%的DMRs，如图3A（左）所示。此外，在基因层面，两种方法识别的DMRs之间有308个基因的重叠，如图3A（右）所示。总体来说，scDMV提供了更大量的信息，几乎捕获了CGmapTools报告的所有DMRs。

考虑到scDMV识别的DMR基因（417）数量比CGmapTools多，对两个基因列表中不同的DMR基因进行功能注释就变得很有趣。通过GO富集分析，发现共有的DMR基因在与发育调控相关的功能，如解剖结构大小的调节和小GTP酶介导的信号转导等方面显著富集（图3B）。

相比之下，scDMV特殊识别的DMR基因在与蛋白磷酸化和神经系统发育调节等过程中高度富集（图3B）。值得注意的是，早期研究已经显示，磷酸化动态在早期发育过程中的调控蛋白质中起主导作用，8-细胞胚胎中的磷酸化蛋白质与翻译后机制相关（Bloom和McConnell 1990、Ju¨ hling等人 2016、Peuchen等人 2017）。DMR基因的功能富集表明，与蛋白质磷酸化相关的基因的DNA甲基化变化可能在胚胎发育中起关键作用，特别是在8-细胞阶段。这些发现为连续发育阶段之间的DNA甲基组动态提供了宝贵的见解。

总的来说，scDMV方法成功捕获了几乎所有CGmapTools识别的DMR。此外，与CGmapTools相比，scDMV揭示了更广泛的重要生物事件，显示出其提供更全面见解的能力。

scDMV: 一种用于处理单细胞二硫化物测序数据的DNA甲基化变异的零–一通胀β混合模型。原文：“scDMV: a zero–one inflated beta mixture model fo...