今天所读的论文题目是《CGmapTools improves the precision of heterozygous SNV calls and supports allele specific methylation detection and visualization in bisulfite-sequencing data 》

论文地址点击蓝色链接：论文地址

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

动机：DNA甲基化对于植物和动物中的基因沉默和印迹非常重要。亚硫酸基测序的最新进步允许检测单核苷酸变异(SNVs)达到高灵敏度，但是准确识别来自部分C-to-T转换序列的杂合SNVs仍然具有挑战性。

结果：我们设计了两种方法，BayesWC和BinomWC，它们显著提高了从约80%到99%的杂合SNV调用的精度，同时保持了可比较的回忆。有了这些SNV调用，我们为等位基因特异性DNA甲基化（ASM）分析和可视化读取上的甲基化状态提供了功能。应用ASM分析到以前的数据集，我们发现平均1.5%的调查区域显示了等位基因甲基化，这些区域在转座子元素中显著富集，可能与相同的细胞类型共享。在低覆盖度数据中用于DMR分析的动态片段策略能够找到与肿瘤发生相关的关键基因的差异甲基化区域（DMRs）。最后，我们将40种应用程序集成到软件包CGmapTools中，以分析DNA甲基组。此包使用CGmap作为格式接口，并设计二进制格式以减小文件大小并支持快速数据检索，可以应用于上下文级别，基因明智，盒质明智，区域明智和样本明智的分析和可视化。

Introduction:

   1. DNA甲基化是一种已在哺乳动物和植物中被广泛研究的表观遗传标记，与基因和转座子元素(TEs)的活动有关。已经开发出许多技术用于描绘全基因组范围的DNA甲基化。DNA的亚硫酸化后进行测序（BS-seq）是构建DNA甲基组的当前黄金标准。全面的DNA甲基化分析为表观遗传调控提供了新的见解。

    2.    单核苷酸变异(SNVs)和DNA甲基化变异有助于基因组和表观基因组的多样性。SNVs可以通过全表观遗传组关联研究(EWAS)，和等位基因特异性DNA甲基化(ASM)计算来探讨。此外，从BS-seq数据中挖掘SNVs在许多应用中都很有用。然而，在BS-seq数据中，SNVs与C到T转化的胞嘧啶混淆，使得更难以区分SNVs与脱氨基化，突变和测序错误。已经开发出若干工具用于BS-seq数据中的SNV调用，已经能够达到良好的敏感度。Bis-SNP和 BS-SNPer都使用了贝叶斯策略，并通过最大后验概率预测显式基因型。由于观察到的读数可能根据C到T的转换而来自未甲基化胞嘧啶的不同基因型，因此无法总是精确地从BS-seq数据中预测出确切的基因型。从BS-seq数据中准确预测SNVs仍然具有挑战性。

     3.   DNA甲基化可以是等位基因特异的，并可以调控印记基因的表达。先前的研究发现，基因组中ASM是普遍存在的。识别ASM介导的基因表达需要先进的ASM调用工具。方等人设计了一个统计模型，无需基因型信息即可预测ASM。然而，这种方法无法区分样本细胞亚群之间的差异甲基化区域（DMR）和父源起源的ASM。 Gao等人开发了SMAP，它支持ASM分析，利用Bis-SNP预测的SNVs。

     4.   关于DNA甲基化的组织特异性调控，DMR分析已被广泛应用。在一些情况下，根据降低表征的亚硫酸测序（RRBS）文库构建了DNA甲基组，这是一种成本效益高的方法，用于测量样本中的部分CpG位点。在低覆盖率文库和间断覆盖区域中，特别是在RRBS数据中，始终很难识别真正的DMR。

      5  已发表的DNA甲基组是用于提取有趣特征和进一步探索的有用资源。由于单核苷酸分辨率的DNA甲基组通常非常大，存储，共享和处理大DNA甲基组数据集的挑战呼唤更好的方法和计算工具。已经开发了几个web服务器来显示收集的DNA甲基组数据，例如MethBase（，NGSmethDB和MethBank。然而，大多数数据库只报告CpG甲基化水平，并且只有限的访问覆盖信息和非CG甲基化。此外，由于带宽限制，从互联网上获取DNA甲基组数据总是一个瓶颈。用于对齐BS-seq数据的大多数工具可以产生SAM格式文件，但使用变异格式来存储甲基化信息这是共享数据的障碍。已经开发了几个用于BS-seq数据分析的下游流程，但大多数只能进行有限的一组分析。

        在这里，我们提出了一种新的方法，通过引入专门为BS-seq数据定义的通配符基因型来识别BS-seq数据中的SNVs。我们的方法显著提高了杂合SNV调用的精度。我们还根据杂合SNV调用开发了ASM分析管道，并提供了在读取中可视化等位基因特异性甲基化的功能。我们将SNV和ASM分析应用于以前的数据集，结果显示平均1.5%的区域是ASM，并且描述了ASM区域的分布和细胞类型特异性。我们还开发了一种动态片段DMR查找方法，该方法特别适用于RRBS和低覆盖WGBS图书馆。将这种方法应用于发布的癌症数据集，我们能够找到与肿瘤发生相关的关键基因相关的DMR。此外，我们将40个功能集成到CGmapTools包中，这为BS-seq数据的下游分析提供了进步。

2 Materials and methods

2.1 SNV calling strategies with wildcards

2.1使用通配符的SNV调用策略

    ATCGmap格式提供了每个碱基对的Watson链和Crick链上的读取计数，我们将其称为ATCGmap表（图1A）。使用ATCGmap表，CGmapTools可以计算基因型。由于在BS-seq数据中C可能被转换为U，因此读取中的T可能表明在未转换的基因组中为T或C。例如，如果ATCGmap表只在Watson链上有一个T，那么该位点可能来自如TT，CC或TC基因型等基因型（图1B）。因此，我们使用通配符来表示预测基因型中的模糊性（Y表示为T或C，R表示为A或G）（表1）。尽管通配符与模糊的基因型相关，但有时我们仍然可以估计它是否为SNV，甚至它是否为杂合SNV。当两条链都有很高的覆盖率时，我们可以解决这种模糊性并计算出确切的基因型。在这里，我们提供了两种SNV调用策略，（i）带有通配符的贝叶斯模型（BayesWC）策略，和（ii）带有通配符的二项式模型（BinomWC）策略。

图 1. 通过引入通配符基因型从BS-seq数据中进行SNV调用。 (A) 一个位置的ATCGmap表的定义。Aw是在Watson链上支持位置为A的Watson链的读取计数。 Ac是Crick链支持Watson链上位置为A（Crick链上的T）的读取计数。(B) 从ATCGmap表中进行基因型预测的示例。以RRBS为例，B的上图只覆盖一条链，读取计数可以来自基因型T，也可以来自基因型C，考虑到亚硫酸盐转化的影响，因此在基因型中引入一个通配符是必要的。B图下半，有两条链上的高覆盖率，来自反向链的信息有助于推断出一个明确的基因型。 (C) 当两条链都有足够的覆盖率时，BinomWC策略的模式。模糊的读取计数被添加到表中的相应位置，然后使用二项式测试从每条链中选择一组核苷酸；然后将两组的交集用作最终预测的基因型。 (D) 在模拟的WGBS数据集中进行杂合子SNV调用的精度分析四种策略。平均覆盖率水平为 10x，20x，30x，40x和50x。（E）在模拟的RRBS数据集中进行杂合子SNV调用的精度分析四种策略。

2.2 BayesWC strategy

在BayesWC中，我们假设基因组为二倍体。因此，基因组类型的后验概率为

其中O是在ATGCmap表中观测到的读数

并且基因型g是纯合的即或杂合的即，即:

其中

注：通配符字符定义为Y=T/C 和 R=A/G

先验定义为：

和

后验可以记作每个观察到的基因型的后验的乘积，即

其中

M和N是两个等位基因的核苷酸

假设将一个核苷酸误报为另一个核苷酸的概率为e，那么正确调用一个核苷酸的概率为p = 1 - 3e。因此我们可以为可能性绘制一个表格（补充表S1）。

我们引入通配符基因型为

其中，Y和R是基因型的通配符。对于通配符基因型，我们计算了后验结果，如下示例所示：

其中，系数的选择考虑了，以及。最后，从精确的基因型集和通配符基因型集中选择具有最高后验概率的基因型作为预测的基因型。

2.3 BinomWC strategy

        另外，我们提出了一种二项式通配符策略，这是前一种用于BS-seq数据的SNV调用方法的修改版本（Orozco等人，2015年）。基本思想是使用二项分布从ATCGmap表中预测基因型。(i) 当两条链都有足够的读数（达到10），会先在每条链上调用核苷酸，然后在两条链上交叉(Fig. 1C)。(ii) 当只有一条链有足够的读数时，将从高覆盖层链预测核苷酸，并可能引入模糊的基因型(Supplementary Fig. S1A)。(iii)当两条链都没有足够的读数时，两条链上的计数被合并为一个六元素向量进行基因型调用(Supplementary Fig. S1B)。

 2.4 Evaluation of SNV calling methods 

 2.4对SNV调用方法的评估

我们使用了一个50 Mbp的参考基因组，并从中生成了一个二倍体基因组，其中的同源和杂合SNVs的比率为0.1%。我们生成了100bp的读数，并按照Guo et al.（2014）的方法比对了全基因组亚硫酸测序（WGBS）和RRBS库。在不同的覆盖等级下，对BayesWC（命令：cgmaptools snv - m bayes - -dynamicP）、BinomWC（命令：cgmaptools snv - m binom）、Bis-SNP（v0.69，默认参数）和BS-SNPer（参数：- -minhetfreq 0.1 - -minhomfreq 0.85 - -minquali 15 - -mincover 0.4 - -maxcover 1000 - -minread2 2 - -errorate 0.02 - -mapvalue 20）的精度和召回率进行了评估。

2.5 Allele-specific methylation regions

2.5 等位基因特异性甲基化区域

        要找到ASM区域，我们首先使用BayesWC识别杂合SNVs。具有模糊基因型的SNVs被丢弃。然后，使用CGmapTools的“asr”模式计算ASM区域。所有通过映射读数连接的杂合SNV位点的区域都被调查。ASM区域应满足以下条件：(i) 至少覆盖两个CpG位点；(ii) 低甲基化等位基因上的mCG水平在0.2以下，并且在高甲基化等位基因上在0.8以上；(iii) 多重t检验的校正P值应为0.05以下。在不同基因组元素的富集研究中，至少覆盖一个碱基被认为是重叠的。为ASM区域和特定基因组区域的富集进行了超几何检验。所有杂合SNV链接的区域都被用作背景。在研究ASM的组织特异性时，进行了超几何检验以评估两个样本之间ASM区域的重叠，如果两个SNVs距离在500bp内，则认为两个样本中的ASM区域是重叠的。

3 Results

3.1 Improved precision of heterozygous SNV calls

        为了评估我们的SNV（单核苷酸变异）调用方法的性能，我们比较了BayesWC和BinomWC与Bis-SNP（Sun et al., 2013）和BS-SNPer（Gao et al., 2015a）在模拟的全基因组亚硫酸盐测序（WGBS）和缩减表示亚硫酸盐测序（RRBS）数据集上的表现。对于WGBS，结果显示BayesWC和BinomWC在杂合SNV调用的精度方面优于其他两种方法（图1D）。特别地，BayesWC达到约99%的精度，无论是对杂合还是纯合SNV调用（图1D和补充图S2）。对于纯合SNV调用，BayesWC与Bis-SNP的精度相当，但在平均覆盖率高于30时，其找回率优于Bis-SNP（补充图S2）。即使针对低覆盖率的数据，BayesWC还能达到高精度；而牺牲的是与其他方法相比较，其找回率低。对于RRBS，我们的结果显示BayesWC对于不同覆盖率范围的杂合SNV调用也具有高精度（图1E), 并且其找回率与Bis-SNP相似，但在覆盖率低时找回率较低。与BS-SNPer相比，BayesWC具有更高的精确度和类似的找回率（补充图S3）。

 3.2 广泛的等位基因特异DNA甲基化丰富在转座子中 

        基于精确预测的杂合SNVs，CGmapTools提供了一个识别ASM(等位基因特异甲基化)区域的功能。基于之前发表的DNA甲基化组（Guo等人，2016），我们计算了SNVs（补充表S2）并使用CGmapTools预测了ASM区域。大量的杂合SNV相关区域在人类和小鼠中显示了两个等位基因上的非对称甲基化水平（补充图S4）。因此，我们用严格的阈值定义了ASM区域。在人类中，平均有1.50%的区域被定义为ASM区域，而在小鼠中这个数字是1.45%（图2A）。此外，我们发现人类和小鼠的ASM区域都在L1重复序列中富集；人类的Alu元素也富集在ASM区域；小鼠的LTR元素在ASM区域中富集（图2B）。

        通过评估所有样本的ASM区域的成对重叠，我们发现ASM区域在相同的细胞类型中，如卵细胞、胚胎干细胞(ESC)和神经元，有显著的重叠。此外，我们发现ESC和神经元的ASM区域相似（图2C）。

图 2. 人类和老鼠的等位基因特异性 DNA 甲基化。(A)散点图显示了在人类和老鼠样本中定义的 ASM 事件的百分比和杂合 SNV 的数量。圆点代表人类样本，三角形代表老鼠样本。(B)在不同的基因组元素中进行 ASM 事件的富集分析，显示了种内基因组偏差的不同。颜色表示超几何检验富集显著性水平。(C)热图显示了人类（左侧面板）和老鼠（右侧面板）细胞类型间 ASM 的一致性。颜色表示使用超几何测试从低到高的一致性显著性水平。(D)已知 dbSNP130 中的 SNP 位点连接的 ASM 代表位点位于 FAM160A 的 5'UTR 和一个 CpG 岛中。根据两个杂合等位基因链接的读取在唐虎露图中为三个大脑样本和一个卵细胞样本代表性地选择。(E)由具有 C 到 T 变换的杂合 SNV 位点链接的代表 ASM 位点，该变换破坏了一个等位基因的 CpG 环境，位于一个已知的印记基因 H19 中。通过 T 链接的读取，通过一个灰色矩形背景标识，无法由于亚硫酸盐转变将模糊读取分配给等位基因 C 或 T。开放圈，未甲基化的 CpG 位点；实心圈，甲基化的 CpG 位点。

CGmapTools中的差异甲基化区域（DMR）分析。(A) 定义动态片段的示意表示。首先，只选择两个样本覆盖的位点；其次，通过定义片段进行基因组扫描，片段应具有最小的胞嘧啶（通常是CpG）计数量n，最大的片段大小S，以及两个相邻胞嘧啶之间的最大距离。灰色圆圈表示只被一个样本覆盖的胞嘧啶位点。然后，使用t检验比较每个片段中的胞嘧啶甲基化水平。实箭头表示片段的扩展，虚线箭头表示终止片段的扩展。(B) chr5上一个区域中的DMR和动态片段（背景）的图形表示。数据来自eRRBS数据集。(C) VCAN基因启动子区域的DMRs的棒棒糖图。箭头指示在(B)中的位置。位点特异性甲基化水平由条形的高度表示。从图中，两个动态片段（灰色盒子）被报告为DMRs，它们在AML中过度甲基化。

        CGmapTools还提供一个功能，用于可视化单个读取上的DNA甲基化状态。我们绘制了一个唐虎露图，该图是用于FAM160A1基因的5’UTR中的一个ASM区域，该区域包含位于CpG岛内的一个CpG-TpG SNV（图2D）。有趣的是，大脑神经元和卵母细胞的不同样本都显示CpG等位基因是低甲基化的，而TpG等位基因是高甲基化的。如上面讨论，映射到Watson链的读取中的T可能源自未甲基化的C的转换或者是T，因此当在唐虎露图中显示时，读取T被标记为模糊的（图2E）。

         之前的研究报告说ASM经常在发现在CpG二核苷酸的杂合SNVs处富集（Shoemaker等，2010）。我们分析了在一组高覆盖度的DNA甲基化组中，ASM区域在全基因组范围内的位置偏好。我们发现，ASM区域的百分比在细胞类型间变化，最小的率在人和老鼠的始原生殖细胞(PGCs)中（补充图S5）。

3.3 Identify DMRs with dynamic fragmentation strategy采用动态片段策略识别差异甲基化区域（DMRs）

        CGmapTools支持差异甲基化位点（DMS）分析和差异甲基化区域（DMR）分析。由于全基因组双硫酸盐测序（WGBS）资源密集，缩减表示双硫酸盐测序（RRBS）是一种更受欢迎的方法，它降低了每个样本的成本并被广泛应用于临床研究。因为RRBS库覆盖的区域是片段化的，我们提出了一种动态片段化策略来识别两个样本之间的DMRs。首先，只有在两个样本中都有足够的读取（≥5次）覆盖的CpG位点被选入DMR分析。利用它们包含最少数量的CpGs，最大长度的碱基，以及两个相邻CpGs之间的最大距离的标准，动态定义背景片段（图3A）。采用非配对的t检验来比较每个背景片段内共享CpG位点的甲基化水平。最后，通过p值和delta甲基化水平的阈值从背景片段中选取DMRs。

        我们将此DMR分析应用于急性髓系白血病（AML）样本和来自已发表的RRBS数据集的正常CD34+骨髓（NBM）对照样本（Akalin等，2012a）。总共选择了53061个动态区域，其中2961个区域是DMRs（P≤0.001和DmC≥0.5）（图3B），并且2161个DMRs在AML样本中过度甲基化。

        为了可视化局部区域的DNA甲基化水平，我们设计了棒棒糖图来区别未甲基化的位点和低覆盖率的位点。据报道，VCAN基因与肿瘤形成和癌症复发有关（Du等，2013）。以VCAN基因为例，棒棒糖图显示了两个动态确定的区域位于转录起始位点的下游，在AMLs中过度甲基化（图3C）。另一个例子显示了癌基因TRIM59的启动子中的一个DMR（补充图S6）。

3.4 Additional features of CGmapTools||CGmapTools的附加功能

3.4.1CGmap is used as the standard format

CGmapTools是BS-seq读取对齐后的下游分析包。因为CGmap和ATCGmap格式报道了与甲基化组相关的详尽信息，且适合进一步处理，我们在CGmapTools中将它们作为标准格式使用。与BS-Seeker2的输出格式类似（Guo等，2013），CGmapTools也提供了直接从BAMs生成两种格式的功能。ATCGmap格式为双链上的所有四种核苷酸提供可读信息。CGmap格式仅为CpG二核苷酸提供读取计数，可以用于大多数DNA甲基化分析（图4）。

受BAM格式（Li等，2009）的启发，我们设计了二进制格式ATCGbz和CGbz，存储排序的DNA甲基化组为omes; (iii)实现为命令行工具，方便并行处理和扩展;(iv) 支持基于二进制文件格式的即时检索以及（v）用户友好功能，可以在多个级别显示甲基化组，如设计Tanghulu图用于可视化原始读取上的甲基化状态，设计棒棒糖图以显示局部区域的低覆盖率和未甲基化的胞嘧啶。比起BAM文件，这些二进制格式占用的空间要少得多，与压缩的ATCGmap/CGmap格式相比，它们也节省了空间（补充图S7）。二进制格式的另一个优点是它们支持从磁盘快速检索。为了方便扩展，CGmapTools中的程序实现为命令行工具，这样用户可以轻松地将它们整合到自己的定制流程中。

CGmapTools生成的图表。 (A) DNA甲基化对样本hBrain_FcM55yNeun_Lis中不同环境的贡献的饼图。 (B) 各种情境中批量DNA甲基化的条形图。(C) 不同环境中DNA甲基化的分布图。(D) 单个样本hOocyte_Rrbs_GuoH的全基因组mCG分布图。 (E) 多样本中DNA甲基化的热图。图中的CpG的平均甲基化水平显示在右侧，基于斯皮尔曼相关系数构建了一个分级聚类树。(F) 基因体内碎片化区域中CpG甲基化水平的分布图。(G) AML样本中染色体范围的MEC（左），MEC的密度图以及MEC的累积分布图（右）为了深入探索等位基因水平上的DNA甲基化组，我们提供了使用精确的杂合子呼叫和dnovo探索SNV相关的ASM区域的流程。 为了深入探索等位基因水平的DNA甲基化组，我们提供了使用精确的杂合子呼叫和新探索SNV相关的ASM区域的管道。 我们的研究表明，ASM区域在转座子元素中显著富集。除了上述功能外，我们将CGmapTools设计成一个集成的DNA甲基化组分析包，具有以下优点：(i)实现动态片段发现的实现；

3.4.2 多功能的分析和可视化功能

关于文件处理，CGmapTools提供了转换文件格式和操作文件的功能，例如文件排序、合并、交集、拆分以及修补丢失的信息等。此外，CGmapTools提供了多样化的功能用于可视化DNA甲基化组。

对于单个样本，CGmapTools能生成展示不同序列上下文中甲基化贡献的饼图（见图5A），也能生成展示批量甲基化(图5B)及甲基化水平分布(图5D)和基因组的甲基化剖面图(图5D)的柱状图。正如Guo等人以前提出的，CW(W是A或T)在哺乳动物中是与CC有别的甲基化上下文（Guo等人，2016）。CGmapTools在植物研究中将CG、CHG和CHH上下文的DNA甲基化进行了报道，而在对哺乳动物的研究中，它报告CG、CW和CC上下文的甲基化。

对于多个样本，CGmapTools生成了展示染色体窗口中DNA甲基化水平的热图(图5E)。为了查看基因体中的mCG水平，CGmapTools还提供了功能，可以刻画基因或特定区域内的甲基化水平分布(图5F).

3.4.3 Evaluate the coverages in BS-seq library

        读取覆盖率是估计DNA甲基化水平和获取SNV调用的重要因素。SNV调用依赖于所有的核苷酸（A、T、C和G），而DNA甲基化水平仅取决于对齐到胞嘧啶的T和C读取计数。因此，我们定义了整体覆盖率（OAC）为所有核苷酸在两条链上的平均读取覆盖率，这可以从ATCGmap文件计算出来。我们还定义了甲基化有效覆盖率（MEC），即仅对胞嘧啶的平均读取覆盖率，这可以从CGmap文件计算出来。通常情况下，MEC稍高于OAC的一半（补充表S2）。CGmapTools还提供了多种工具用于可视化覆盖率的分布（图5G）。

4. 讨论

        在这项研究中，我们使用了新颖的BayesWC和BinomWC方法，这两种方法显著提高了杂合SNV调用的精确度，与以前的工具相比，并保持了可比较的召回率。为了深入探索等位基因水平的DNA甲基化组，我们提供了利用精确的杂合调用和对SNV相关ASM区域的新颖探究的流程。我们的研究表明，ASM区域在转座子元素中显著富集。

        除了上述功能外，我们设计了CGmapTools作为一个集成的DNA甲基化分析包，它的优点包括：（i）实施一种动态分段策略，用于探索低覆盖率数据中的DMR；（ii）使用标准的ATCGmap/CGmap格式，便于分享甲基化组；（iii）作为命令行工具实现,便于并行处理和扩展；（iv）支持基于二进制文件格式的即时检索；（v）用户友好的功能，可在多个层次上可视化甲基化组，例如设计一个Tanghulu图用于可视化原始读取上的甲基化状态，并设计一个Lollipop图来揭示局部区域内的低覆盖胞嘧啶和非甲基化胞嘧啶。最后，CGmapTools为分析BS-seq数据中的计算挑战提供了先进且丰富的解决方案。