请用加权甲基化水平来计算区域 DNA 甲基化

链接: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3523709/

随着测序技术的发展，测序越来越便宜，高通量测序的文章数不胜数，然而我自己以后也不知道猴年马月才用 WGBS 技术了，刚好疫情在家，看点代表性综述（翻译），记录一下下，巩固一下基础知识。如有不当之处，还恳请轻喷。。

本文为大牛 Joseph R. Ecker 实验室于 2012 年在 Trends in Genetics 上发表的一篇关于 DNA 甲基化值定义的综述（估计是当时看到各种各样的计算方法，实在看不下去了）‘Leveling’ the playing field for analyses of single-base resolution DNA methylomes。本文主要介绍了几种计算位点、区域 DNA 甲基化值的方法，以及作者自己的看法。

大佬实验室直达：https://ecker.salk.edu/publications/

更好阅读体验移至：https://www.yuque.com/docs/share/b084b7b3-32ef-4af6-8453-2120293308b2?#

太长不看系列：

加权甲基化水平应该是研究 DNA 甲基化的默认指标

摘要：

在过去的几年里，DNA 测序的成本大幅下降，而测序 reads 的数量和长度却增加了。这一测序革命导致了广泛采用的方法来研究全基因组的 DNA 甲基化模式，统称为全基因组重亚硫酸氢盐测序( WGBS: whole-genome bisulfite sequencing )。目前，通过高通量测序和重亚硫酸氢钠转化（用于检测胞嘧啶 DNA 甲基化的黄金标准方法）相结合的方法，可以对单碱基分辨率的 DNA 甲基化进行常规解码 [1,2]。随着对 DNA甲基化数据的获取和分析的增加，越来越需要对特定胞嘧啶或区域的甲基化量的定义达成共识。“甲基化水平”的定义通常很模糊，而且根据实验者和所处理的问题的不同会有很大的不同。因此，这里提出了一套从 WGBS 数据分析“甲基化水平”时需要考虑的指导原则。

Single-site methylation level：单位点甲基化水平

WGBS 允许检测单个胞嘧啶的甲基化状态。该过程使用重亚硫酸氢钠将未甲基化的胞嘧啶转化为尿嘧啶（甲基化的胞嘧啶不受影响，仍然为胞嘧啶），最终通过 PCR 将其转化为胸腺嘧啶 [3]。然后，通过对转换后的产物进行排序并将数据 mapping 到参考基因组，就可以检测到这些信息。假设参考基因组中该位点为胞嘧啶碱基 C ，若所得到的 reads 中该位点为胸腺嘧啶 T ，那么说明该位点胞嘧啶 C 未甲基化，若所得到的 reads 中该位点为胞嘧啶 C ，那么说明该位点胞嘧啶 C 甲基化。尽管在单个细胞中，胞嘧啶要么甲基化，要么未甲基化，但为了达到测序或检测包含杂合区域的细胞所需的 DNA 的最小数量，通常得到的是多个细胞的混合样 DNA。因此，这些实验产生了一个不同测序 reads 的集合，其中一些表明一个特定的胞嘧啶被甲基化，另一些表明同一位点未被甲基化。在 CG 位点的背景下，DNA 两条链链之间的甲基化是高度一致的。所以合并两条链的 reads 数通常来说是没问题的（把对称的位置算作一个单位）。通常，使用二项式检验来确定观察到的甲基化频率是否高于重亚硫酸氢盐转化率无效和测序错误所预期的背景值。使用该检验时，给定位点的甲基化量通常表示为甲基化 reads（例如：reads 中该位点为 C 的 reads 数目） 与覆盖该位置的 reads 总数比值（例如：reads 中该位点为 C 或者 T 的 reads 数目）（ 图 1 ）。我们将这种位点特异性指标称为位点的甲基化水平。

图 1：一个甲基化区域的示例场景和计算“甲基化水平”的多种方法。
(a) 纯合区域 a 等位基因和 (b) b 等位基因和 (c) 杂合子的不同 profiles的 DNA 甲基化区域的例子。
绿色的“C”表示甲基化的胞嘧啶，蓝色的“C”表示二项式检验确定的未甲基化胞嘧啶。黑色的“C”方框表示没有测序 reads 覆盖特定的胞嘧啶。
每个胞嘧啶上面的分数表示一个胞嘧啶的 reads 除以该位置的胞嘧啶和胸腺嘧啶的 reads 总数。
(d) 使用 a-c 区域，使用每种讨论的方法计算“甲基化水平”。
C = 支持甲基化胞嘧啶的 reads 数，T = 支持未甲基化胞嘧啶的 reads 数，i =胞嘧啶的位置，n =基因组该区域胞嘧啶的数目，M = 二项式检验确定位置为甲基化时的指示变量 1。
从左侧第 2 位开始计算的结果被用来计算“位点甲基化水平”。
此外，对于任何含有蓝色‘C’的位置，将含有胞嘧啶的 reads 进行平方，并不包括在文中所述的甲基化水平的计算中。
缩写：mC, methylcytosine。

虽然单位点的数量可能与生物学有关，但研究人员通常对某一区域而不是特定位点的甲基化水平感兴趣。有越来越多的可以鉴定差异的甲基化区域的计算机程序（它本身是变化的来源，但由于空间限制，这里不考虑它），但是一旦这些区域被检测到，就有多种方法来计算甲基化水平，即使使用相同的 WGBS 数据，通常也会导致非常不同的值。

Fraction of methylated cytosines：甲基化胞嘧啶的比例

最简单的方法是结合一个区域的位点特异性甲基化信息，计算在统计上显示显著甲基化量的胞嘧啶的比例（例如：甲基化胞嘧啶的比例）。如果一个人对不同区域的甲基化位点的潜力感兴趣，这个指标是有用的，因为一个重要的二项式检验表明，群体中至少有一个细胞/等位基因在该区域有一个甲基化胞嘧啶。如前所述，每个位点的 reads 数据代表了对整个人群中甲基化状态的调查，而甲基化胞嘧啶的比例指标并没有捕捉到关于每个位点甲基化水平差异的信息。这可能很重要，因为在特定位点甲基化的细胞比例的改变可能表明基本的表型改变（比如：在癌症 [4] 和发育 [5] 中）。此外，将这种方法广泛应用于杂合区域可能会有问题，因为如果其中一个等位基因甲基化，大量检测到的甲基化胞核嘧啶将通过二项式检验（图 1 ）。这一缺陷可以通过筛选那些包含可以将甲基化信息与特定等位基因连接起来的杂合遗传变异的测序 reads 来克服。

# A allele 中因为没有通过二项检验的位点甲基化看成 0，通过的则为 1，求甲基化的比例
> (1 + 1 + 1 + 0 + 1 + 1 + 0 + 1 + 1 + 1 + 1 + 1) / 12
[1] 0.8333333

约等于 0.833

Mean methylation level：平均甲基化水平

在一个区域的甲基化水平度量中包含额外信息的一种方法是取该区域内各位点的甲基化水平的算术平均值（图一）。虽然甲基化胞嘧啶的比例有所提高，但这种方法并没有考虑到区域内不同位点的可变测序覆盖率。由于每个位点对总体平均值的贡献相等，因此该方法假设每个位点的信息内容也相等；然而，更高测序深度的的测序位点将为该区域的平均甲基化水平提供更准确的估计。

# 因为没有通过二项检验的位点甲基化 reads 应该看成 0， A allele 中有两个蓝色 C 
> (2/2 + 3/3 + 20/21 + 0/18 + 4/30 + 4/25 + 0/25 + 3/17 + 3/20 + 6/7 + 4/6 + 5/5) / 12
[1] 0.5079995

约等于 0.508

Weighted methylation level：加权甲基化水平

因此，可能希望通过每个位点的测序深度来加权每个位点对一个区域中水平的贡献量（图 1 ）。例如，假设在一个区域中有两个位点，其中一个位点 100 条 reads 中甲基化的 reads 数为 90，而另一个位点总共才 2 条 reads，但是甲基化 reads 数为 1。该区域的甲基化水平将有很大的不同，这取决于是对第一个位点施加更多的权重还是对两个位点赋予相同的权重（平均甲基化水平）（图 1）。

Total_C = (2 + 3 + 20 + 2 + 4 + 4 + 2 + 3 + 3 + 6 + 4 + 5) - (2 + 2) = 54
# 因为没有通过二项检验的位点甲基化 reads 应该看成 0， A allele 中有两个蓝色 C ，故减去对应的 C 的数目

Total_C_T = 2 + 3 + 21 + 18 + 30 + 25 + 25 + 17 + 20 + 7 + 6 + 5 = 179

mC = 54 / 179 = 0.301676 ≈ 0.302

其他注意事项和标准化方法

计算平均值和加权甲基化水平的一个关键点是，尽管二项式检验认为未甲基化的位点仍然应该包括在这些计算中，但它们不应该为计算提供任何甲基化 reads（ Box 1 ）。换句话说，未通过二项式检验的位点的甲基化 reads 应该设置为零，而不管在该位点检测到的甲基化 reads 是多少。同样重要的是，所描述的所有指标都只考虑了许多样本中特定区域内的胞嘧啶。人们可能会试图通过一个区域的大小（例如，一个区域的甲基化水平）来标准化该区域的甲基化水平。但这并没有考虑到不同基因组片段间碱基组成的差异。因此，甲基化的显著差异可能仅仅是由碱基组成的差异造成的。(例如，一个胞核嘧啶含量较少的区域)，这通常不能提供信息。

Box 1:
计算甲基化水平的其他注意事项
在某些情况下，在计算平均值或加权甲基化水平时，将来自二项式检验未确定具有统计学意义的位置的胞嘧啶计数包括在内。我们不包括这些计数，因为二项式检验的目的是确定在高于实验噪音的地方是否存在甲基化。因此，如果检验失败，在该位点发现的胞嘧啶应该被认为是未甲基化的、未转化的碱基。如果将这些位点包括在内，将会夸大甲基化水平的估计。这一变化的幅度将取决于该区域的大小和覆盖范围以及未甲基化位点的数量，但为了保守起见，我们省略了它们。诚然，一些通过二项式检验的位点也会包含这些未转化的胞嘧啶，但一旦它们通过了这项测试，该位点至少会有一些甲基化。如果对未转化的胞嘧啶有顾虑，可以从平均值或加权甲基化水平减去估计的未转化率。此外，由于 WGBS 是链特异的，我们倾向于忽略支持在胞嘧啶位置进行 A 或 G 碱基调用的 reads 操作，因为这些操作代表了测序错误，不能准确地用于计算。基于这些原因，我们认为 C 和 T 在统计上有显著意义的位点的 reads 代表了我们对真实甲基化水平的最佳估计。

结束语

作者建议，加权甲基化水平应该是研究 DNA 甲基化的默认指标，因为它是最广泛适用的。这并不是说这里描述的替代指标没有价值（例如：如上所述的甲基化水平比例）。因此，对科学家来说，精确地定义他们试图用 DNA 甲基化分析来回答的问题是至关重要的，这样他们就可以正确地计算出捕获生物学兴趣的甲基化水平，从而对结果进行准确的评估。

个人的一点看法

尽管 DNA 甲基化分析开发不计其数的软件，但是软件之间的计算 DMR 吻合性几乎很差，仍然有很大的进步空间，不论在实验技术还是分析方法上。尽管现在已经发表了各式各样的文章，但是。。。 草，算了。