摘要

动机： 

     单细胞亚硫酸盐测序（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数据由一个集合表示，记为

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

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

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

2.1制定模型和测试统计量

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

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

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

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

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

以及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）。