用甲基化模式估计全基因组DNA甲基化异质性

　　在异质细胞群中，单个细胞的行为不同，对环境的反应也不同。这种细胞多样性可以通过测量DNA甲基化模式来评估。具有可变甲基化模式的位点是细胞异质性的信息，可能作为疾病和发育进展的生物标志物。细胞间甲基化异质性可以通过单细胞甲基组或池细胞的计算技术进行评估。然而，这些方法精确估计甲基化异质性的可行性和性能需要进一步评估。

　　在这里，我们提出了基于模型的方法，采用了最初来自生物多样性的数学框架，以估计全基因组DNA甲基化异质性。通过特征对比对模型和现有方法的性能进行了评价，并在合成数据集和真实数据集上进行了测试。总的来说，我们的方法比其他方法有优势，因为它们与实际异质性有更好的相关性。我们还证明甲基化异质性提供了不同于传统甲基化水平的额外生物学信息层。在案例研究中，我们发现CG和非CG甲基化的不同甲基化异质性可以预测拟南芥基因组元件之间的调节作用。这为植物表观基因组学研究开辟了新的方向。最后，我们证明了我们的评分可能能够识别人类癌症样本中的基因座，作为早期癌症检测的假定生物标志物。

　　我们将生物多样性的数学框架引入三种基于模型的方法来分析全基因组DNA甲基化异质性，以监测细胞异质性。我们的方法，即MeH，已经被实施，用现有的方法进行评估，并对研究界开放。

　　DNA甲基化是发生在胞嘧啶上的一种可遗传的表观遗传修饰，在转录调控、发育规划[1]和疾病进展[2]等许多生物学过程中起着至关重要的作用。全基因组DNA甲基化及其与基因表达的关联已经通过最新的下一代测序方法[3]得到了广泛的研究，如亚硫酸盐测序(BS-seq)[4,5]和酶促甲基测序(EM-seq)[6]。当从甲基reads数据转换的reads(即从BS-seq或EM-seq的reads)与参考基因组对齐时，可以确定特定胞嘧啶的甲基化状态(甲基化或未甲基化)。在批量甲基化测序(如BS-seq和EM-seq)中，汇集了数百万个细胞，代表了可能与其可变DNA甲基化谱相关的异质细胞混合物(图1A)。在给定的基因组区域内对齐的甲基读取产生由多个胞嘧啶形成的甲基化模式，代表单个细胞(图1B)。基因组区域的甲基化模式可以从完全甲基化到完全未甲基化。中间模式可能表明细胞间DNA甲基化的变化。

　　图1

　　DNA甲基化模式的插图作为监测细胞发育的代理。DNA甲基化模式用于监测细胞异质性，可能由细胞发育、细胞类型混合、分化、遗传变化、疾病、压力引起。黑(白)点表示甲基化(未甲基化)胞嘧啶。不同的颜色对应不同的细胞亚型(使用BioRender.com创建)。B DNA甲基化模式由排列的甲基读段(通过BS-seq或EM-seq)与颜色编码的甲基化状态形成。一组甲基化单倍型被圈起来，以表示覆盖特定位点的特定分数所考虑的读数

　　DNA甲基化异质性在一个特定的位点被定义为DNA甲基化模式之间的变化在这个位点，在一个细胞池。甲基化异质性可能由遗传或表观遗传因素的多种表观遗传调控引起[7]。例如，启动子的可变甲基化通常与对环境刺激或细胞发育进程的转录反应有关[8,9](图1A)。由于群体中单个细胞的行为可能不相同，这可能是由于遗传变化通常伴随着表观基因组变化，或者在细胞分化的情况下，表观遗传变化。通过监测甲基化模式的变化，甲基化异质性可能能够在生物发育或疾病进展过程中捕获遗传或表观遗传因素的指纹。

　　已经开发了评估甲基化异质性的实验方法和计算方法。在实验策略中，单细胞BS-seq (scBS-seq)主要用于研究小鼠稀有细胞(如卵母细胞)的甲基化异质性[7,10]。scBS-seq可以直接测量甲基化异质性，通过单细胞亚硫酸盐测序，然后使用BPRMeth[11]、Melissa[12]和scMET[13]等工具进行计算分析，通过甲基化谱对单细胞进行估算和聚类。尽管如此，scBS-seq也有其自身的成本和挑战，例如在制备文库时需要分离单个细胞，读取作图比率低，由于必须测序的细胞数量而导致的成本高[7]，亚硫酸氢盐处理导致的DNA严重损失。此外，由于植物细胞的细胞大小和细胞壁较大，原有的哺乳动物单细胞测序方案不能直接适用于植物单细胞测序，阻碍了植物的研究。尽管使用第三代测序技术的DNA甲基化谱分析不需要亚硫酸盐转换[14,15]，但据报道碱基调用的高测序错误率超过15%[15]，甲基化调用的错误率高达40%[16]。因此，已经尝试使用甲基化序列汇集细胞基因组区域胞嘧啶的甲基化状态来计算量化甲基化异质性。

　　迄今为止，只有少数计算方法[3,17,18,19,20,21,22,23]和数据库[24]被开发出来，用于使用来自池细胞的全基因组甲基化数据来估计甲基化异质性(参见表1的方法及其特征列表)。其中大多数被设计用来比较基因组区域内一致或不一致甲基化的甲基化和未甲基化的胞嘧啶或读对的数量(methyl- concurrency;MC[22]，不一致读数比例;PDR[20]，不一致读对分数;FDRP[18])，并扩展到基于概率的设置(Epipolymorphism;EP[19])或纳入模式相似性(定量FDRP;qFDRP [18]， Methylation Haplotype Load [21];MHL)。MHL考虑模式相似性，因为它计算所有可能长度的子串在每个读段中完全甲基化的比例，这也使其具有与DNA甲基化水平相似的特征[18]。qFDRP比较读对的相似性，对于可能来自中间甲基化区域的不一致的读对，其权重更高。因此，qFDRP评分可能不完全独立于甲基化水平。此外，设计了基于香农熵的方法来估计类似于异质性的混沌程度(甲基化熵;我[23])。

　　表1计算结果全基因组DNA甲基化异质性评分的几种方法

　　虽然这些方法在各自的思想和优势上存在差异，但对方法的评估与DNA甲基化的性质和下一代测序的数据格式有关。因此，这里列出了在实现过程中需要考虑的建议特性。首先，分析非CG位点甲基化的能力，即CHG和CHH背景下的胞嘧啶(H=A, C或T)。事实上，在植物和动物的基因组中，非CG位点的数量超过CG位点。在植物甚至真菌中，非cg甲基化在许多重要的生物过程中起着至关重要的作用[1]。其次，评分线性，即得分与潜在甲基化异质性之间的线性相关。线性评分能够公平地评估所有基因组区域的异质性，以及样本之间不同异质性的位点。如附加文件1:图S1所示，非线性评分不太可能忠实地反映低异质性和高异质性位点之间的潜在异质性(即偏斜)。第三，考虑不同甲基化模式之间的相似性。具有相同甲基化模式的细胞可能来自相同的细胞亚群。两种高度相似的模式可能是由同一群体中细胞的一些胞嘧啶引起的甲基化逐渐变化引起的。通过将所有模式区分开来而忽略模式相似性可能会丢失关于细胞发育的微妙信息[25]。第四，在评分中混淆甲基化异质性和甲基化水平;这样的评分很容易与甲基化水平的估计相混淆，在甲基化水平的估计中，一些模式比其他模式更重要(例如，完全甲基化vs.未甲基化)。混淆甲基化水平的评分可以从估计所有模式都被平等考虑的甲基化异质性的原始想法中转移出来。最后，全基因组筛选的能力对用户群体尤为重要。该方法的实施应允许全基因组筛选和多个样本之间的比较，从而能够检测到具有可变异质性的位点。根据甲基化异构估计器的这些特性考虑，我们总结了一个涵盖几种流行实现的特性表(参见表1)。

　　在本研究中，我们引入了一系列基于Chao等人[26]的数学模型的多样性指数，这些指数在量化生物多样性方面已被证明是成功的。生物多样性可以解释为物种或类型的有效数量。我们采用该框架及其特定的变体模型来量化甲基化异质性。Chao等人[27]的模型(见式1)中，Hill数[28]即有效类型数是一个q阶多样性指数的参数族，表示相同丰富类型的数量。希尔数需要使类型的平均比例丰度等于在感兴趣的数据集中观察到的丰度。在Chao等人的框架中，set被认为是实体的集合。对于中的每个实体，其属性值为，其丰度为。中实体的总丰度由其对应的丰度加权的属性和给出。因此，实体的相对丰度为，中所有实体的属性与对应丰度的乘积之和等于1，。基于特定属性的set (order)的属性多样性由统一框架给出如下:

　　（1）

　　其中是决定模型对相对丰度敏感性的参数。当，属性的丰度对公式没有贡献，并给出所有实体时的香农熵指数[29]，它根据属性的丰度对属性进行加权。时，它是辛普森指数的倒数形式[30]，它被发现在不同情况下提供了对多样性的稳健估计。Chao等人通过改变集合和属性值，基于物种多样性[31]、系统发育多样性[28]和基于距离的功能多样性[27]等不同属性，提出了一个涵盖Hill数主要变量的统一框架。

　　当基础模型(Eq. 1)用于测量生物多样性时，set被视为不同物种的集合，该物种的属性值是该物种的函数，用于描述例如该物种的种群大小，或与另一个物种或与集合中所有物种的相对相似性。物种的丰度就是该物种的种群大小。如上所述的其他变量要么是归一化因子，要么是与物种不直接相关的模型参数。

　　同样，当考虑在特定位点(如基因组区域)观察到的甲基化模式时，我们假设估计该位点的甲基化异质性类似于测量特定领域内的生物多样性。在不考虑甲基化模式之间相似性的简单设置下，甲基化异质性基本模型(Eq. 1)中的变量可以翻译为:

　　1.

　　集合被认为是从特定位点的比对中观察到的不同甲基化模式的集合;和

　　2.

　　是pattern的属性值。它可以是模式的丰富度，或者在扩展设置中，模式之间的成对相似性，或者在对齐中所有观察到的模式之间的相似性(参见方法;方程式。3 - 6);

　　3.

　　是模式的丰度，可以通过具有该模式的读取次数来估计。

　　可以修改这种设置以适应模式相似性，使用模式对而不是模式作为实体以及相应的属性值。

　　Chao等人的数学基础模型(附加文件1:注S1)具有尺度不变性、弱单调性和倍增性等数学性质。其中，弱单调性是指当增加一个与现有物种有较大差异的新物种时，多样性会增加[32]。倍增性是指将具有相同多样性的两个相互排斥的群体以相同的权重混合在一起，导致基于距离的多样性增加一倍或增加四倍[27]。倍增特性被认为是分集估计所需的基本特性。请注意，虽然所有基于Hill数的变量都表现出加倍性质，但香农熵[29]和辛普森指数[30]都没有表现出这种性质。

　　我们发现Chao等人的希尔数框架(附加文件1:注释S1)中的数学性质是使用甲基化模式量化甲基化异质性的理想选择，原因有几个。首先，用甲基化模式测量表观遗传多样性类似于测量物种的生物多样性，在分析甲基化异质性方面，数学模型在概念上是新颖的。其次，数学性质及其意义使得异质性得分在样本间更具可解释性和可比性;因此，我们采用这些特性来估计DNA甲基化异质性。结果表明，基于甲基化异质性丰度(AB)、基于配对相似性(PWS)和基于系统发育树(PHY)的三种模型方法可以通过分析合成数据，并以scBS-seq数据为基准，帮助克服现有方法的缺点。最后，我们提供的例子显示了该方法在拟南芥非cg甲基化异质性分析和监测癌症样本疾病进展方面的优势。此外，我们的模型的计算程序是为了实现甲基化异质性(MeH)的基因组筛选估计[33]，并比较在不同条件下获得的样本。我们的MeH实现和教程可以在https://github.com/PaoyangLab/MeH上公开获得。

　　摘要。

　　背景

　　结果

　　讨论

　　结论

　　方法

　　数据和材料的可用性

　　缩写

　　参考文献。

　　致谢。

　　作者信息

　　道德声明

　　# # # # #

　　我们首先展示了现有方法和我们提出的方法的行为与玩具的合成排列的例子。随后，为了证明评分的主要特征，即线性，我们合并了多个单细胞甲基组，以测试异质性的估计是否随着细胞数量的增加而增加。此外，我们表明，比较甲基化异质性可以揭示样品之间的差异，这些差异可能仅通过甲基化水平无法检测到。最后，为了在实际数据上验证我们的方法，我们分析了拟南芥甲基组以分析非cg甲基化异质性和人类结直肠癌数据。

　　我们编制了一个评估表，对几种流行的现有方法和我们基于模型的方法进行了评估(表1)。该表包括甲基化异质性估计和实现中需要考虑的具体特征。总的来说，我们基于模型的方法相对于现有方法的主要优势在于可能扩展到非cg甲基化位点，评分线性和考虑甲基化模式之间的相似性，以进行无偏和有意义的评估。

　　为了评估我们提出的方法和其他现有方法精确检测甲基化异质性变化的能力，我们创建了具有可变甲基化模式的玩具示例。为了确保两种方法之间的公平比较，构建了完全对齐的reads来表示完整的甲基化模式，并模拟了类似基因组区域比对的甲基化模式组合(见图2A-C，顶部面板)。首先，我们检验了甲基化异质性会随着新模式的出现而增加的假设。因此，我们提出的三种模型以及FDRP、qFDRP、MHL、ME和EP都呈现出单调增加的趋势(见图2A虚线)，而其他方法表现出不同的趋势。接下来，我们检查了模式相似性在模型中的重要性，因为甲基化模式可能是由与甲基化维持相关的逐渐变化引起的。如图2B和C所示，我们期望当模式变得更加多样化(从左手到右手)时，理想分数会增加。我们发现只有PWS、MHL和qFDRP能够检测到甲基化模式的这种变化(见图2B和C)。因此，只有PWS和qFDRP两种方法符合这两个假设。qFDRP的一个特别关注的问题是，其评分的设计使其很容易受到前面描述的甲基化水平的影响(见表1)。我们基于模型的方法，特别是PWS方法，证明了有效平衡所有这些特征的能力。因此，使用PWS进行以下分析。

　　图2

　　甲基化异质性方法的评价。A-C用合成数据集估计甲基化异质性。上图列出了不同位点甲基化模式的组合。圆形是基于模型的方法，三角形是现有的方法。虚线表示有增加趋势的方法。四种类型的分数用于比较，基于模型的方法:AB, PWS和PHY;基于一致性的方法:MC、PDR、FDRP;基于熵的:ME，基于概率的:EP;现有的考虑模式相似度的方法有qFDRP和MHL。D . PWS和ME估计合并小鼠ESC和肌肉单细胞甲基化异质性。E针对不同数量的ESC单细胞甲基化组绘制了全基因组甲基化异质性比率。黑线表示假设合并的单元格都是异质的期望值，而红色表示线性

　　在检查真实数据时，我们预计甲基化异质性水平会随着新模式的引入而增加。我们处理了来自两种不同细胞类型(肌肉细胞和胚胎干细胞)的小鼠[10]的许多单细胞甲基组。我们使用PWS方法进行分析，因为它是图2A-C中唯一通过前面评价的方法。我们还将ME和EP纳入单细胞分析，因为它们共享类似的全基因组分析数据输入格式。

　　我们预计甲基化异质性评分会从一种细胞类型的单细胞甲基化组增加到混合两种细胞类型的甲基化组。为此，在估计异质性之前，将两种细胞类型的6个单细胞甲基化组合并(即，将多个单细胞甲基化组的原始reads合并成一个合并的甲基化组，然后将其与参考基因组对齐)。我们发现，正如预期的那样，PWS的甲基化异质性在混合细胞类型中增加了(见图2D)。相比之下，在相同的数据集上运行ME方法，混合细胞类型显示出较低的甲基化异质性评分。

　　接下来，我们评估了已知组成的不同数量的ESC单细胞甲基化组之间的甲基化异质性(见图2E)。在一个完美的环境中，添加更多相同类型的细胞不会增加异质性。然而，在ESC的实际数据中，每个单细胞甲基组可能无法覆盖ESC的所有预期模式;基于亚硫酸氢盐处理导致的大部分DNA损伤，观察到的模式可能在这些单细胞甲基化组之间非常不同。因此，随着新模式(来自新添加的单个细胞，甚至来自同一细胞类型)的加入，甲基化异质性评分可能会增加，我们还希望观察到甲基化模式逐渐饱和，异质性趋于稳定。首先，将不同数量的单细胞甲基组，即6、8、10、12、14、16和18组合为合并甲基组，以模拟批量测序数据。对于那些来自许多细胞的甲基化组，我们预计总体上甲基化异质性(表明细胞异质性)高于那些来自较少细胞的甲基化组。

　　我们计算了每种选择方法的全基因组甲基化异质性比率(有关程序和计算的详细信息，附加文件1:注释S3)。总的来说，我们观察到随着单细胞甲基组数量的增加，所有方法的比例都单调增加。为PWS绘制了从6到8个甲基组的外推线(红色)，以证明每添加2个甲基组，预期的线性增加。这条线还显示，ME和EP很可能很快达到明显偏离线性的平台，这表明这两个分数在检测新模式方面不太敏感。这种非线性所揭示的低灵敏度在实际数据应用中不太受欢迎(参见附加文件1:图S1的演示)，特别是在不同样本或区域进行比较时。此外，我们发现这些方法都没有完美地显示加倍特性(见图2E黑色实线)。这可能发生，因为在实际数据中，ESC的这些单细胞通常不是相互排斥的组。尽管如此，我们发现PWS的异质性与其他方法相比是相对线性的。当单细胞甲基化组的数量增加时，它也显示出较小的偏差，这使它成为一个合理的评分。

　　为了确定甲基化异质性和常用的甲基化水平度量之间的差异，我们绘制了通过PWS方法估计的甲基化异质性与人类结直肠癌(CRC)(图3A)和拟南芥野生型甲基组(附加文件1:图S2)的3个重复样本的甲基化水平的对比图。如图3A所示，散点图表明甲基化异质性和甲基化之间的关系在不同的胞嘧啶环境下(即CG, CHG和CHH, H=A, C或T)有所不同。我们观察到甲基化异质性和甲基化水平在CpG甲基化区域之间呈曲线状关系，甲基化异质性较高的区域具有中等甲基化水平，在人类和拟南芥中都发现了这种关系。这些区域可能反映了在基因区域中常见的表观基因组变化的动态过程(图3C)。在拟南芥中，我们还分析了非cg位点的甲基化异质性。我们发现，与CG位点相比，非ch位点(即CHG和CHH)与甲基化水平的关系非常不同(附加文件1:图S2)。虽然非cg位点的甲基化程度较低，但其中一小部分甲基化程度较高，且甲基化异质性较高。值得注意的是，中度甲基化区域在所有情况下都保持着多样化的异质性，这在单独使用甲基化水平进行评估时很容易被忽视。此外，我们的PWS评分能够在甲基化水平变化不明显的情况下检测甲基化异质性的变化。简而言之，甲基化异质性可以潜在地补充甲基化水平的使用，以识别使用甲基化水平无法检测到的微小变化，并提供与传统甲基化水平不同的生物信息层。

　　图3

　　全基因组甲基化异质性谱。平均甲基化异质性与相邻正常CRC样本的3个重复的平均甲基化水平绘制。B拟南芥基因组不同基因组特征中高(前10%)和低(后10%)甲基化异质性区域的比例。拟南芥高表达基因和低表达基因(上、下25%)间CG甲基化异质性的C - meta图。拟南芥高、低表达TEs及其邻近区域CHG甲基化异质性的D meta图

　　为了揭示植物全基因组甲基化异质性，我们使用PWS对拟南芥野生型甲基组进行了分析，其覆盖范围为58X[34]。我们发现高甲基化异质性区域优先靶向CG和非CG位点的转座元件(te)，这与低甲基化异质性区域不同(图3B)。此外，与非CG位点相比，CG位点的高甲基化异质性在基因体上大量富集;表明对CG和非CG的偏好有所不同。随后，我们比较了高表达和低表达的基因和TEs(顶部和底部25%)的甲基化异质性，见图3C和d。我们观察到CG甲基化异质性与转录起始位点(TSS)附近的基因表达呈负相关，其次是转录末端位点(TEs)的正相关;这表明DNA甲基化在启动子和基因体上的动态表观遗传调控。我们还发现，与高表达TEs相比，低表达TEs表现出更高的CHG甲基化异质性(见图3D)，这表明活性TEs的甲基化模式在植物细胞中是高度可变的。

　　我们的研究产生了植物中甲基化异质性的第一张地图。高甲基化异质性区域位于特定的基因组特征，在CG和非CG甲基化异质性之间存在差异。甲基化异质性被证明与转录调控有关。我们的研究结果阐明了CG和非CG甲基化异质性在拟南芥基因组中的独特功能。

　　接下来，我们想要证明具有差异甲基化异质性的基因组区域也可以被认为是感兴趣表型的生物标志物。我们下载并处理了来自CRC的人类减少亚硫酸氢盐测序(RRBS)甲基组数据[35]，其中包括不同阶段，包括III-IV期CRC冷冻肿瘤(肿瘤)，来自同一患者的病原体显示的正常粘膜(正常)，以及从切除肿瘤两侧边缘收集的组织学证实的匹配正常样本(邻近正常)。最初的研究分析了每个阶段的10个样本，发现特定癌症基因的启动子甲基化增加了40%，从而触发肿瘤的转录变化，而在邻近的正常基因中，启动子甲基化仅增加了20%，而表达没有变化，可能是由于启动子甲基化的较低变化不足以触发转录变化。

　　为了验证我们的方法，我们使用PWS方法对每个正常、邻近正常和肿瘤样本进行了3个重复的甲基化异质性分析。目的是看看我们的PWS方法是否能够识别假定的生物标志物，作为现有方法(如EWAS)的替代方法。许多DNA甲基化水平的研究已经表明，在正常组织和正常邻近组织之间存在甲基化差异[36]，或者在正常组织和有癌症发展风险的正常组织之间存在甲基化差异[37]。因此，我们并行分析DNA甲基化水平，以评估甲基化异质性和甲基化水平之间的可预测性(附加文件1:图S4为差异甲基化区域DMRs和差异异质区域DHRs的维恩图)。

　　在相邻的正常和正常样本之间(n=911)以及肿瘤和正常样本之间(n=1558)共鉴定出2319例dhr。这些DHR主要存在于基因体中(图4A，左面板)。在对RRBS基因组进行正常化后，我们发现DHR在启动子、外显子、5 ' -和3 ' UTR上富集，但在内含子上不富集;表明可能与转录有关。

　　图4

　　甲基化异质性与甲基化水平的比较及PWS异质性的评价。A不同基因组特征下dhr的比例(左)和富集图(右)。B dhg和dmg的维恩图。C CRC 3期dhg甲基化异质性热图。使用PWS估计CPXM2周围甲基化异质性的IGV图，DHR为橙色阴影。每个蓝色条表示400 bp箱中的平均甲基化异质性，并且标记了DHR条的确切值。E .在DHR内特定CG窗口4内所有样品中甲基化模式的组成，以确定潜在的疾病模式

　　我们将dhr与共定位基因关联后，鉴定出953个差异异质基因(dhg)(图4B)，而仅检测到14个差异甲基化基因(dmg)(甲基化变化15%，p值< 0.05)，其中包括原研究中报道的2个基因C9orf69和RAPGEFL(见dhg和dmg的鉴定方法)。只有一种DMG，即fk506结合蛋白10 (FKBP10)，被发现也是DHG。这可能表明基于甲基化水平和异质性的分析实际上针对的是不同的基因集。

　　为了跟踪不同阶段异质性的变化，我们绘制了一个甲基化异质性热图，使用肿瘤dhg和邻近的正常dhg(图4C)。热图显示，从正常、接近正常到肿瘤，甲基化异质性有明显的变化，其中大多数基因增加了它们对肿瘤的异质性。类似的dmg热图无法准确显示样本组之间的差异，因为一个正常样本被归类为癌症组(附加文件1:图S5)。通过独创性通路分析[38]，我们分析了邻近正常(附加文件1:图S6)或肿瘤(附加文件1:图S7)特异性非重叠dhg的丰富功能。丰富的疾病和功能清楚地表明，通过比较邻近正常样本和正常样本鉴定的dhg参与结直肠癌相关疾病;这表明这些基因甲基化异质性的变化与癌症进展高度相关，邻近的正常dhg可预测结直肠癌肿瘤。总之，DHG分析在选择与感兴趣表型相关的区域方面补充了传统的DMG方法。

　　我们总共鉴定出162个基因在肿瘤dhg和邻近的正常dhg之间显著重叠(重叠测试;p < 0.00001)(图4B;附加文件1:图8)。这些基因显示出DNA甲基化对肿瘤形成的强烈和持续的异质性变化。这些基因甲基化异质性的变化可能表明，由于癌症形成或细胞分化等原因，特定的甲基化模式随着细胞类型的变化而出现。

　　作为一个例子，我们从162个重叠的dhg中找到了CPXM2。CPXM2是一种蛋白质编码基因，据报道与多种人类疾病有关，如发育性疾病[39]、阿尔茨海默病和精神分裂症[40]，并且在激活时促进肿瘤侵袭性[41]。如CPXM2甲基化异质性的截图所示(图4D)，在相邻正常样本、正常样本、肿瘤样本和正常样本之间的比较中，在启动子处不断发现重叠的dhr。这个特定DHR的甲基化模式的组成(图4E)揭示了一个特定的甲基化模式，标记为橙色的“1111”(一排完全甲基化的胞嘧啶)，似乎是一种“疾病”模式。它不存在于正常样本中，但它开始出现在邻近的正常样本中，并在肿瘤样本中稳定下来。此外，当完全甲基化的“1111”模式增加或其他部分甲基化的与“1111”非常相似的模式(如“0111”或“1011”等)存在时，显示这种模式的读取比例增加;与在正常样本中观察到的大多数模式类似于未甲基化的“0000”模式不同，读取开始变得类似于模式“1111”。这验证了我们的模型检测甲基化模式变化的能力，甲基化模式可能作为疾病早期检测的生物标志物。

　　在本研究中，我们提出了一种采用Chao等人的生物多样性数学框架[26]的方法来估计甲基化异质性。我们基于模型的方法随后被应用于MeH程序，从甲基序列数据中估计全基因组甲基化异质性。我们的研究结果证明了MeH在多个甲基组亚群中突出不同甲基化模式的能力。与仅检测不同模式的现有方法不同，我们基于Hill数的方法具有实现无偏估计的数学特性，并且能够分析读取之间的模式相似性。此外，MeH可以用作评估CG和非CG甲基化异质性的工具，并询问在癌症发展过程中预先指定的细胞群中甲基化模式的变化。

　　我们的分析表明，现有的一些甲基化异质性方法给出了非线性评分，有些方法不适应甲基化模式之间的相似性;两者都会导致不太有利的结果。换句话说，当使用基于ME或EP的方法评估甲基化异质性时，在样本或基因组位置的比较中，当不同模式的数量较少时，变化往往被高估，而当不同模式的数量较多时，变化往往被低估。不同的量表之间会有差异，如果我们平等对待它们，可能会导致错误的发现(重大变化)。其他分数中对模式相似性的忽略，例如基于一致性的方法，原则上可以进行修改，以便在未来的开发中考虑到这一点。此外，由于测序错误已被证明会引入异质性评分的偏差[18]，因此使用模拟亚硫酸酯测序数据评估了包括PWS、ME和MP在内的评分的行为(见附加文件1:注释S4)。虽然所有方法都受到测序误差的影响，但当引入5%的测序误差时，PWS的变化较小(3.7%)，重复之间的变化较小。考虑到目前普遍的测序误差小于1-2%，PWS仍然是一种很好的选择。

　　另一方面，实现基于窗口筛选的方法(如基于模型的，ME和EP)的一个警告是，只考虑完整的甲基化模式(即，完全对齐的读取)。这可能会导致大量数据丢失。然而，在目前的NGS技术下，通常使用的150-200个碱基对的读长足以覆盖4个CpG二核苷酸，并且测序成本在过去几年中显著降低。还有一些插入工具可以通过添加模拟邻近甲基化模式的模拟reads来克服低覆盖率的问题，例如BSImp[42,43]，以最大限度地利用有限测序深度的数据。

　　我们基于一个框架开发了三种变体，该框架考虑了使用甲基化模式量化甲基化异质性的不同属性，并展示了每种变体的优势和劣势，同时使用玩具示例和单细胞甲基组对现有方法进行了评估。虽然每种方法都有自己的弱点，有些是高度相关的，但没有一个是可以替代的。[18]进行了全面的比较，并指出了每个评分应该使用的场景。特别是，我们发现我们的模型改善了非线性问题，因此有利于跨样本的系统全基因组比较，并为进一步分析确定dhr。

　　在我们使用单细胞甲基化组进行分析之后，我们还展示了基于散装甲基化组和合并甲基化组来估计细胞类型组成的方法的潜力，因为这种效果似乎是令人上瘾的。我们预计，添加额外的单细胞甲基组将导致甲基化异质性增加，这可以通过使用8个和6个单细胞甲基组估计的全基因组甲基化异质性之间的差异来估计(图2E)。然后，给定特定甲基组的全基因组甲基化异质性，我们可以潜在地估计细胞组成。这种对一种类型细胞的分析，但最终可能会详细阐述纳入不同类型细胞的策略(图2D)，这可能有利于诸如区分细胞群体组成的研究，以提高疾病或发育背景下细胞分化的准确性。

　　图4E所示的甲基化模式的组成证明了使用DHR识别可能存在的特定甲基化模式的潜力。这种特定DHR位置的甲基化模式可以被认为是特定表型的生物标志物。完全甲基化胞嘧啶的“1111”模式在邻近的正常和肿瘤样本中出现的比例非常低，并且只有在正常细胞中常见的甲基化模式发生很大比例变化后才会出现。这有两个含义。首先，甲基化模式的变化要么是一个渐进的过程，要么甲基化模式的出现依赖于中间模式的存在。在任何一种情况下，甲基化异质性可能比甲基化水平更敏感。这进一步表明在量化甲基化异质性时考虑甲基化模式之间的相似性的重要性，因为并非所有模式都具有同样的独特性。

　　通常很难确定甲基化与基因表达之间的关联，更不用说甲基化与基于基因表达的表型之间的关联了;然而，甲基化异质性被认为与表型直接相关，因此做出了假设，并对数据进行了分析。同样，这种方法是独立于甲基化水平，并被证实是如此。虽然我们没有研究表观遗传调控的机制，但拟南芥的甲基化谱表明甲基化异质性与基因和TE的表达呈负相关，这实际上可能为甲基化与转录组学的联合分析开辟了新的方向。最后，对结直肠癌病例模式的追踪确实表明，在研究疾病进展时，甲基化异质性可能是比甲基化水平更有效的指标。

　　最终，MeH(如图5所示的工作流程)可以使用提出的基于模型的方法来描述全基因组甲基化异质性。该方法为用户提供了根据胞嘧啶位点和甲基化环境(包括所有CG, CHG和CHH环境)的数量来指定窗口大小的自由，用于甲基化异质性的评估，并且是同类方法中的第一个。我们还嵌入了最近开发的一种甲基组插入方法[42]，以最大限度地覆盖评估甲基化异质性的有限偏差，如附加文件1所示:图S9。在估计甲基化水平时，该方法估算甲基化状态的准确率超过85%，偏差仅为~ 3%。尽管许多关于细胞异质性的研究都集中在哺乳动物数据上，但涉及不同DNA甲基转移酶的甲基化调控等重要主题的研究只能在植物上进行，而甲基化在其他情况下很常见，如CHG和CHH[44]。最后，我们提供了一个应用为生物多样性开发的数学模型来估计甲基化异质性的例子。

　　图5

　　MeH工作流程概述。MeH的函数在黄色区域显示。每个方框代表一个组件，对应于一系列任务。数据处理步骤用红色表示，步骤中使用的工具用斜体表示，虚线框表示该步骤是可选的。输入数据(蓝色框)和输出数据(橙色框)的方括号中表示文件类型。

　　ccDownload: /内容/ pdf / 10.1186 / s13072 - 023 - 00521 - 7. - pdf