期刊： Proceedings of the National Academy of Sciences

DOI： 10.1073/pnas.2019768118

作者： O Y Olivia Tse; Peiyong Jiang; Suk Hang Cheng; Wenlei Peng; Huimin Shang; et al

出版日期 ：2021-01-25

网址： https://doi.org/10.1073/pnas.2019768118

摘要

5-甲 基胞嘧啶 (5mC)是 一种重要的表观遗传修饰。亚硫酸氢盐测序 (BS-seq)有局限性，如严重的 DNA 降解。利用单分子实时测序，我们开发了一种直接检测 5mC 的方法。该方法在测量窗口内全面检查 DNA 聚合酶的动力学信号 (包括脉冲间隔时间和脉冲宽度 )和每个核苷酸的序列背景，称为整体动力学 (HK)模型。每个分析的双链 DNA 分子的测量窗口由 21个核苷酸组成，在中心的 CpG位点上有一个胞嘧啶。我们使用扩增的 DNA(未甲基化 )和 M.SssI 处理的 DNA(甲基化)(M.SssI 是一种 CpG 甲基转移酶 )来训练卷积神经网络。我们使用这些样本来区分甲基化状态的曲线下面积高达 0.97。单碱基分辨率全基因组 5mC 检测的灵敏度和特异性分别达到 90%和 94%。然后在人鼠杂交片段上测试 HK模型，其中杂交的每个成员具有不同的甲基化状态。该模型还在从各种生物样本中提取的人类基因组 DNA 分子上进行了测试，如褐毛、胎盘和肿瘤组织。 HK 模型推断的总体甲基化

引言

该方法使用较温和的条件将 5mC 转化为胸腺嘧啶，试图克服 BS- seq 中存在的局限性。然而，TAPS涉及多个酶和化学反应步骤，包括 TET 氧化，吡啶硼烷还原和 PCR 扩增。在任何 DNA 转化步骤中发生的不期望的转化效率都会对 5mC 分析的准确性产生不利影响。

我们设想一种理想的测量碱基修饰的方法是一种可以直接应用于天然 DNA 的方法，在测序之前不需要任何 DNA 的化学/酶转化和 PCR 扩增。第三代测序技术，如纳米孔测序(如牛津纳米孔技术公司)和单分子实时(SMRT)测序(如太平洋生物科学公司，PacBio)，使单分子实时测序成为可能，为探索这些检测碱基修饰的方法提供了机会。

水平与 BS-seq 模型推断的甲基化水平具有良好的相关性(r = 0.99;P<0.0001)， 并允许测量印迹基因中等位基因特异性的甲基化模式。综上所述，这种方法提供了一个同时进行全基因组遗传和表观遗传分析的系统。

表观遗传学|表观遗传学碱基

第三代测序修饰

意义

单分子实时(SMRT)测序理论上提供了直接评估天然 DNA 分子的某些碱基修饰的机会， 而无需任何事先的化学 /酶转化和 PCR扩增，使用 DNA 聚合酶的动力学信号。然而， 5mC 修饰引起的动力学信号变化是非常微妙的。因此，还没有实现 5mC 修饰的全基因组测量。我们通过整体分析 DNA聚合酶的动力学信号和测量窗口内每个碱基的序列背景，使用 SMRT 测序增强了5mC 检测。我们使用卷积神经网络来训练甲基化分类模型，从而实现全基因组 5mC 检测。该方法的灵敏度和特异度分别达到 90%和 94%， 总体甲基化水平与亚硫酸氢盐测序的相关性为99%。

DNA 甲基化是一种将甲基共价添加到 DNA 分子上的生物过程。

该过程最常见的形式发生在胞嘧啶嘧啶环的第 5 位:即 5-甲基胞嘧啶(5mC)。DNA 甲基化在细胞的表观遗传调控中起着重要作用，包括基因组印迹、x 染色体失活和致癌作用(1,2)。最广泛使用的检测 5mC 的方法包括亚硫酸氢盐处理，然后是 PCR 或大规模平行 DNA 测序(即亚硫酸氢盐测序[BS-seq])(3,4)。然而，这种基于亚硫酸氢盐的技术存在显著缺陷。例如，亚硫酸氢盐处理的恶劣反应条件可以降解大部分输入 DNA(5)。这种 DNA降解使得长 DNA 分子测序具有挑战性。另一个缺点是亚硫酸氢盐诱导的 DNA 降解优先作用于富含未甲基化 cy-的基因组区域不同亚硫酸氢盐处理方案导 致整体甲基化的高估 和特定基因组区域的实 质性差 异 (4 ) 。最近， 一种无 亚硫酸 氢盐 的方法 ( 称为纳米孔测序检测碱基修饰的可行性(7 )。然而，测序结果往往伴随着高测序错误，如插入和缺失(8 )。这些错误会导致引入许多位点，缺失了甲基化分析所需的信号。这样的局限性目前的纳米孔测序技术可能会阻碍在单分子水平上解码甲基化模式，尤其是对于人类基因组这样的大基因组。相比之下，纳米孔测序最多读取两次 DNA 模板(即，包括沃森和克里克)

图 1 所示。 利用单分子测序和 HK 模型检测 5mC 的示意图。双链 DNA 分子与发夹接头连接，形成圆形 DNA 模板。ZMW 中的 DN A 聚合酶会将不同荧光团 标记的核苷酸结合到 DNA 模板的互补链中， 从而发出不同的荧光颜色，表示核苷酸信息:例如，红色、黄色、绿色和蓝色分别代表 G、C、T 和 a。光脉冲 信号反映了 DNA 聚合酶动力学，取决于碱基修饰。脉冲信号包括 IPD 和 PW。对于甲基化分析的胞嘧啶， ipd、PWs 和围绕胞嘧啶的序列上下文被组织成一 个数据矩阵，称为测量窗口。为了说明目的，在 CpG 位点内的胞嘧啶上游和下游的 10 nt 表示为沃森链的 5'-G[CCA TGC] A T A C G T T [ GA TG CA ] a -3 '。为了 简单起见，括号内的碱基被省略(用“…” 表示)。在这种情况下，测量窗口大小，包括中间的询问胞嘧啶，为 21 nt。对于-3 位置对应腺嘌呤(“ a”)的碱基， 与“ a” 相关的 IPD(1.8 )和 PW(0 .7)填充在“-3” 列和“ a” 行之间的相应细胞中。同列的其他细胞用“ 0” 填充。其余与 21-nt 序列相关的 ipd 和 PWs 根据相 同的规则填充在该测量窗口中。来自克里克链(' 5- T[ TTG CA T] CAA CG TA [ TG CA TG] G-3 ')的动力学信号和序列背景也被类似地处理。结合两个互补的 CpG位点(即沃森链和克里克链)的测量窗口用于下游分析。 使用来自甲基化和非甲基化胞嘧啶的多个组合测量窗口来训练 CNN， 以区分测试样本中的甲基化和非甲基化胞嘧啶。 CN N 包括输入层、卷积层和输出层。将测量窗口输入输入层，然后进行卷积层处理;然后， 通过基于 sig moid 函数的输出层生成 CpG 的甲基化概率(范围:0 到 1)。这种方法被称为“ 整体动力学(HK)模型”(HK 模型)。

SMRT 测序依赖于环形 DNA 模板的创建，该模板允许对分子进行多次测序，从而大大提高碱基调用的准确性(9)。理论上， 在此， 我们将 以沃森链的 动力学信号 和序列上 下文的数据 处理碱基修饰会影响 DNA 合成过程中 DNA 聚合酶的动力学。例如，当胸腺嘧啶(T)与模板上相反的 n6 -甲基腺嘌呤(6mA)结合时，DNA 聚合酶的反应速度会减慢，导致当前碱基与下一个碱基结合的时间间隔增加(10)。染料标记的核苷酸发出的脉冲信号可用于监测聚合速度的这些变化，从而能够检测碱基修饰 (10)。例如，脉冲间持续时间(IPD)(即两个连续荧光脉冲之间的时间间隔)可用于识别腺嘌呤(6mA)的 n6 甲基化。与 6mA 检测不同，据我们所知，目前还没有报道过使用 SMRT 测序的方法来实现天然 DNA 分子 5mC 全基因组检测的实际有意义的准确性。5mC 检测的挑战是由 DNA 聚合酶动力学的微妙变化引起的，通过这种变化，鸟嘌呤与 5mC 相反。例如，Clark 等人报道使用 IPD 对 CpG 位点内的胞嘧啶进行 5mC 的检出率很低，从 1.9 到 4.3%不等(11)。在本研究中，我们试图开发一种方法，通过整体利用序列上下文和与 DNA 聚合酶动力学相关的脉冲信号(称为整体动力学 (HK)模型)，利用 SM RT 测序实现 5mC 的准确检测。基于 HK模型，我们使用甲基化和未甲基化数据集训练卷积神经网络 (CNN)来检测 5mC 修饰。

结果

HK 模 型检测 5mC 的原理

如图 1 所示，用发夹适配器连接双链天然 DNA 分子，形成拓扑圆形的 DNA 模板。测序引物通过发夹接头上的互补位点退火成环状 DNA 模板。环状 DNA 模板与 DNA 聚合酶结合，形成复合物，每个复合物随后固定在零模波导(zmw)的底部。ZMW 中的 DNA 聚合酶分子催化了用不同荧光团标记的核苷酸结合到 DNA 模板的互补链中。DNA聚合酶在聚合过程中的动力学变化可以在单分子基础上进行监测。

使用不同的荧光染料测定碱基含量。例如，红色、黄色、绿色和蓝色分别代表 G、C、T 和 A(图 1)。荧光标记的核苷酸发出的光脉冲信号反映了 DNA 聚合酶动力学，这取决于碱基修饰。因此，适当使用脉冲信号将有可能确定一个胞嘧啶是否被甲基化。脉冲信号包括 IPD 和脉冲宽度(PW)，前者表示连续两次碱基结合之间的时间间隔，后者表示与碱基结合相关的荧光信号发射的时间持续时间。脉冲信号与聚合反应发生的序列环境相关。在此，我们开发了一种通过使用脉冲信号来确定 DNA 甲基化的方法，包括 ipd、PWs 和序列上下文。序列上下文是指一段 DNA 中的碱基组成(A、C、G 或 T)和碱基顺序。对于 CpG 位点内的胞嘧啶，围绕胞嘧啶的 ipd、PWs 和序列上下文被组织成一个数据矩阵，称为 a

为例。在模板 D NA 的 C p G 位点内被询问的胞嘧啶的位置表示为位置 0 。为便于说明， 沃森和克里克链模板包含 10 个核苷酸 (nt ) 的 上 游 和 下 游 的 胞 嘧 啶 分 别 表 示 为 5 '- G[C C AT GC ]AT ACGT T [G AT GC A] A -3  '  和  5  '  - T [ T T GCAT ]C A AC GT A [T GC A T G] G -3 ' 。为 了 简单 起见 ， 图 1

中省略了括号内的碱基。在这种情况下，包括胞嘧啶本身 ( 在中心) 在内的测量窗口大小为 2 1 n t 。对于-3 对应腺嘌呤( “ A ”)碱 基的位置，与“ A ” 相关的 IP D(1 .8 ) 和 P W(0 .7 ) 填充在“-3 ”列 和“A ” 行之间的交叉处( 称为细胞) 。其他位于-3 列和胞嘧啶(C )、鸟嘌呤(G) 、胸腺嘧啶( T )之间的细胞用“0 ” 填充。其他与 21 -nt 序列相关的 i p d 和 P Ws 在该测量窗口的相应细胞中填充。对源 自克里克链的动力学信号和序列上下文进行了类似的处理。

由于几乎所有甲基化的 CpG 位点都对称地出现在两条链上 (12)，我们将沃森链上的 CpG 位点侧翼的测量窗口与克里克链上配对的 CpG 位点侧翼的测量窗口结合起来，形成了一个用于下游分析的组合测量窗口。我们使用了一些来自甲基化和未甲基化胞嘧啶的组合测量窗口来训练 CNN。训练后的 CNN 模型将用于区分测试样本中的甲基化和未甲基化胞嘧啶。这种 5mC 检测的分析框架整体地利用了测量窗口内单个核苷酸之间 DNA 聚合酶的动力学信号，以及序列上下文(即核苷酸信息和顺序)，因此被称为“整体动力学(HK) 模型”(HK 模型)。 HK 模型包括一个输入层、卷积层和一个输出层。来自每个测量窗口的 HK 模型所需的数据(即序列上下文、IPD 和 PW)被输入到输入层，然后由卷积层进行处理(图 1)。基于 s 型函数的输出表示甲基化的概率，称为 CpG 位点胞嘧啶的甲基化评分，范围从 0 到 1。由于是二元分类，CpG 位点内胞嘧啶未甲基化的概率为 1 -甲基化评分。甲基化得分越大，CpG 位点越有可能被甲基化。基于接收算子特征(ROC)曲线，定义了甲基化评分阈值，用于对分析的 DNA 分子中每个 CpG 位点的甲基化状态进行分类。有关培训和测试程序的详细信息见材料和方法。

使用扩增和m.s sssi处理的 DNA 训练 5mC 检测的 HK 模型。

为了证明 使用 HK 模型以全基因组方式确定甲基化状态的可行性和性 能，该模型使用 SMRT 测序数据集进行训练和验证，包括未 甲基化数据集(即阴性数据集)和甲基化数据集( 即阳性数据集)。未甲基化数据集包含通过全基因组扩增(WGA) 制备的扩增 DNA 的测序结果(表示为 WGA 数据集)。在 WGA 中使用未 修饰的核苷酸导致扩增的 DNA 几乎不含碱基修饰(除了少量 输入基因组 DNA)。甲基化数据集包含测序结果测量窗口(图 1) 。由于圆形分子可以多次测序， 因此使用测量窗口内每个核苷酸的平均 IP D 和 PW 值进行下游分析。m.s sssi(一种 CpG 甲基转移酶，从含有 Sprioplasma sp.菌株 MQ1 甲基转移酶基因的大肠埃希菌菌株中分离出来)处理的 DNA 结果表明，在此之前，甲基化双链 DNA 中的所有 CpG位点

测序(用 m.s sssi 处理数据集表示)。m.s sssi 甲基转移酶使 CpG位点甲基化(13)。在 m.s sssi 处理样本数据集中测序的 CpG 位点中，有一半用于训练 HK 模型。在 WGA 数据集中，随机抽取相同数量的 CpG 站点用于训练 HK 模型。使用 m.s sssi 处理样本数据集中剩余的一半 CpG 位点和 WGA 数据集中相同的数量来验证模型。在本研究中，我们在 PacBio Sequel I 测序仪上使用 1)Sequel I 测序试剂盒 3.0(正式名称为 Sequel 测序试剂盒 3.0)，在 PacBio Sequel II 测序仪上使用 Sequel II 测序试剂盒 1.0 和 Sequel II 测序试剂盒 2.0，获得 WGA 和 m.s sssi 处理 DNA 数据集，用于本研究中不同试剂和测序仪间的 HK 模型评估。

对于 Sequel I 测序试剂盒 3.0，我们使用了来自 m.s sssi 处理的 DNA 样本(完全甲基化)的 328,233 个 CpG 位点和来自 WGA样本(完全未甲基化)的 328,233 个 CpG 位点来训练 HK 模型。

m.s sssi 处理数据集的甲基化评分(中位数:0.99;四分位数间距 [IQR]: 0.93 至 1.0)与 WGA 数据集的结果分离(中位数:0.04;IQR: 0.02 ~ 0.1) (P 值< 0.0001,Mann-Whitney U 检验)(图 2A)。ROC曲线下面积(AUC)为 0.97(图 2B)。

进一步分析不同试剂盒制备的 S M RT 测序数据集。在 S eq u elII ki t 1 .0 和 S eq u el II k i t 2 .0 制备的两个训练数据集中， 在甲基化评分方面， W GA 和 m .s s s si 处理数据集之间的分离也可以清楚地观察到( 图 2 A) ， 两个数据集分别使用了 1 1 ,2 7 2 ,5 5 2 和 3 2 5 ,7 8 0个 C p G 位点。S eq u el II 测序试剂盒 1 .0 和 2 .0 制备的数据集的 AUC 值分别为 0 .9 6 和 0 .9 4 (图 2 B )。

HK 模型使用扩增和m.ssi处理的 DNA检测 5mC的性能。

图 2 C 和 D 显 示了 HK 模型在测试数据集中的性能。Sequel I 测序试剂盒 3.0 和 Sequel II 测序试剂盒 1.0 和 2.0 的 AUC 值分别为 0.97、0.96 和 0.93。这些结果表明，HK 模型可以准确地确定甲基化状态。只要训练和测试过程基于相同的实验条件，HK 模型适用于由 不同测序试剂盒和测序仪产生的数据。

在所有三个测试数据集(SI 附录， 图 S 1 a -c)中，基于 HK 模型的 AUC 值(范围: 0 .9 3 至 0 .9 7 )远远大于基于 IP D 或 PW 值的 C p G位点的 AU C 值(0 .5 3 至 0 .6 7 )，表明 HK 模型在查询基础上大大优于使用动力学信号的传统方法。

我们定义了一个甲基化评分截止值来对 CpG 位点的甲基化 状态进行分类。我们选择 0.5 作为甲基化评分的截断值，这是 在训练数据集中接近每个 ROC 曲线左上角的点。甲基化评分 大于 0.5 的 CpG 位点被归类为甲基化;否则归为非甲基化。使 用 Sequel I 试剂盒 3.0 产生的数据集，我们可以达到 94%的特 异性和 90%的敏感性。对于 Sequel II 试剂盒 1.0 生成的数据集，特异度和灵敏度分别为 92% 和 87%。对于 Sequel II 测序试剂盒 3.0 生成的数据集，特异度和灵敏度分别为 89%和 83%。

除了 CNN 模型外，我们还尝试使用隐马尔可夫模型(HMM) 对 SMRT-seq 具有高深度测序平均值的 BC01 样本进行 5mC 检测的性能评估(SI 附录，表 S2)。结果，我们发现 HMM 模 型(83%的敏感性和 84%的特异性)的表现似乎比 HK 模型(87% 的敏感性和 92% 的特异性)差。关于 HMM 实现的细节在 SI 附 录，方法和材料(SI 附录，图 1)中描述。si 附录，图 S2 和 S3)。

窗口大小和子读深度对 5mC 检测性能的影响。

为了研究测量窗口的窗口大小和子读深度如何影响 HK  模型的性能，我们改变了测量窗口大小， 覆盖了 1 ,3 ,5 ,7 ,9 ,1 1 ,2 1 ,3 1 ,4 1 ,5 1 和 6 1 n t 。对于特定的测 量 窗 口 大 小 ， 我 们 进 一 步 改 变 了 子 读 深 度 ， 覆 盖 了 1 ,2 ,3 ,4 ,5 ,1 0 ,1 5 ,2 0 ,2 5 和 30 倍。HK 模型首先使用比较来自 W GA和 m .s ss si 处理数据集的 1 0 0 ,0 0 0 个测量窗口的训练数据集进行训练。对于窗 口大小和子 读深度的每 种组合， 我们从与训 练数据集不重叠的完整数据集中随机抽取 2 0 0 0 个 C p G 位点， 从而形成测试数据集。我们分析了子读深度为 10×的 Sequel II 测序试剂盒 1.0 生成的数据集。测量窗口大小为 1 时，AUC 为 0.70(SI 附录，图 S4A)。当我们将测量窗口大小增加到 3,7,21 和 31 nt 时，AUC C   D    值分别增加到 0.84,0.90,0.93 和 0.93(SI 附录，图 S4A)。此外，

当应用 21 nt 的测量窗口大小时，使用 1 的亚读深度发现 AUC 为 0.75(SI 附录，图 S4B)。当我们将亚读深度增加到 5 和 10 时，观察到 AUC 值分别增加到 0.85 和 0.93(SI 附录，图 S4B)。这 些数据表明，HK 模型的性能可以通过调整测量窗口大小和子 读深度要求来提高。SI 附录，图 S4C 显示了不同的性能

图 2 所 示。 使用扩增 DNA 和 m.s sssi 处理 DNA 生成的数据集对 HK 模型进行训练和验证。(A) 来自全基因组扩增 DN A ( WGA D NA 数据集)和 m.s sssi处理 DNA (m.s sssi 处理 DNA 数据集)的训练数据集中甲基化分数的箱形图不同测序试剂盒的基础包括 Sequel I 测序试剂盒 3.0 和 Sequel I I 测序试剂盒1.0 和 2.0。( B)基于不同测序试剂盒的训练数据集 ROC 曲线。( C)测试数据集中甲基化得分的箱形图。(D)测试数据集的 ROC 曲线。

从未甲基化的胞嘧啶中分离甲基化的胞嘧啶达到一个平台， AUC为 0 .9 6 ，窗口大小为 2 1 nt ，亚读深度为 30 倍。在子读段深度为 10 ×的情况下，使用 2 1nt 的测量窗口大小也使我们能够达到一个好的结果。

方法

通过单分子实时测序对胞嘧啶甲基化进行全基因组检测

AUC 为 0.93。为了平衡适合下游分析的分子数量和准确性，我们在本研究中采用了 21  nt  的窗口大小和至少 10×的子读深作为默认设置。人类单倍体基因组中有 2820 万个 CpG 位点，产生 2820 万个测量窗口。其中，69.4%的 21  nt  测量窗口包含一个 CpG 位点。分别有 21.2%和 6.5%的测量窗口包含两个和三个 CpG 位点。只有 3%的测量窗口包含三个以上的 CpG 位点。因此，我们认为大多数测量不会受到位于同一测量窗口内的两个附近 CpG 位点引起的动力学信号的潜在相互作用的影响。结果与 S eq u el I 测序试剂盒 3 .0 和 S eq u el II 测序试剂盒 2 .0 生在 SMR T 测序过程中。因此， 在这样的训练中，基因组序列背景可能存在一定程度的冗余。

使用人 -鼠杂交片段分析不同的甲基化状态。

由于上述验证过程依赖于 WGA 和 m.s sssi 处理的 DNA 样本，这些样本在理论上是片段均质甲基化或未甲基化的，因此我们进一步测试了 HK 模型是否可以推广到携带异质甲基化状态的片段(即同时包含甲基化和未甲基化 CpG 位点的片段)。为此，我们基于限制性酶切 (HindIII 和 NcoI，均为 6 碱基切割)和 DNA 结扎(材料和方法)

成的数据集相关(S I 附录，图 2) 。S5 -S 7 )得出了一个一致的结论， 生成了包含人鼠杂交片段的两个数据集，如 SI 附录，图 S8 所即 HK 模型的性能将取决于窗口大小和子读深度。亚读深度的增加通常会增加区分甲基化和非甲基化胞嘧啶的 A UC 值。2 1 n t 的测量窗口大小是甲 基化分析的稳 健参数，因为 这样的窗口大小似乎在亚读深度为 30 倍时达到平台值(S I 附录， 图 s 7 a 和 B) 。有趣的 是 ， 在 窗 口大 小 和亚 读 深度 范 围内 ， 相 对 较早 的 试剂 盒示。第一个数据集包含人类部分被 m.s sssi 甲基化而小鼠部分被 WGA 甲基化的杂交 DNA 分子，命名为人类(甲基)-小鼠 (unmeth)数据集。第二个数据集包含具有相反甲基化模式的杂交 DNA 分子:即人类部分未甲基化而小鼠部分甲基化，命名为人类(unmeth) -小鼠(meth)数据集。我们使用 Sequel II 测序仪和S eq u el I 测序试剂盒 3 .0 似乎优于其他两个更新的试剂盒。例如， Sequel II 测序试剂盒 1.0 对样本 H01 和 H02 进行测序，得到S eq u el I 测序试剂盒 3 .0 (SI 附录，图 S 7 A) 和 S eq u el II 测序试剂盒 1 .0 (S I 附录， 图 S 4 C ) 和 2 .0 (S I 附录， 图 S 7 B )的测量窗口大小为 21 nt，亚层深度为 30 时，A UC 值分别为 0 .9 8 ,0 .9 6 和 0 .9 4 。

序列上下文的数量对 5mC 检测性能的影响

人类参考基因组共有 2820

570 万(中位数:1.3 kb;子读取深度中位数:10×)和 330 万(大小中位数:1.2 kb;中位子读深度:10.5×)分子，分别用于人类(meth) -小鼠(unmeth)和人类(unmeth) -小鼠(meth)数据集。我们应用从具有均匀甲基化状态的数据集训练的 HK 模型来个 CpG 位点(University of California Santa Cruz hg19)。其中， 确定人(冰毒)-小鼠(冰毒)数据集中每个人-小鼠杂交 DNA 分子

以 CpG 位点为中心的 21-nt 序列上下文共 2070 万个。如 SI 附录表 S1 所示，在 2070 万个上下文中，Sequel I kit v3、Sequel II kit v1 和 Sequel II kit v2 制备的数据集，用于 HK 模型训练的序列上下文的百分比分别为 2.7%、32.7%和 1.3%。由于我们对使用 Sequel II 试剂盒 v1 制备的训练样本获得了更高的测序通量，因此该样本中经经验覆盖的序列上下文更多。每个测试数据集包含与其对应的训练数据集相似数量的序列上下文。值得注意的是，即使在训练和测试过程中涵盖了不同数据集的不同数量的上下文，所得到的 HK 模型的性能似乎变化不大，接收器工作特征曲线下的面积(AUC)值在 0.93 至 0.97 之间。

为了进一步研究序列上下文的数量如何影响 HK 模型的性能，的 CpG 位点的甲基化状态。根据与限制性内切酶识别位点 (HindIII 或 NcoI)最近碱基的相对位置(即距离)，我们共收集了限制性内切位点上下游 50 个碱基对(bp)内的 104,896 个 CpG 位点。来自分子人类部分的位置被定为上游(负值)，来自小鼠部分的位置被定为下游(正值)。被确定为甲基化的 CpG 位点的百分比被视为甲基化水平。图 3A 显示，在这个人(甲基)-小鼠 (unmeth)数据集中，人类部分显示为甲基化，甲基化水平范围为 85.9 ~ 93.0%，而小鼠部分显示为非甲基化，甲基化水平范围为 6.7 ~ 9.6%。在人类(unmeth) -小鼠(meth)数据集中发现了相反的模式(图 3B)。

我们进一步分析了限制性内切酶位点两侧最近的两个 C p G 位我们对序列上下文进行了降采样分析，随机抽取 1,000、5,000、 点， 评估了邻近 C p G 位点的动力学信号的潜在相互作用对 HK50,000 、 100,000 、 200,000 、 300,000 、 400,000 、 500,000 、1,000,000、5,000,000 和 10,000,000 个序列上下文。SI 附录，图 S4D 显示，随着序列上下文数量的增加，HK 模型的性能逐渐提高。例如，使用 1,000 个序列上下文时，AUC 为 0.73，而使用 300,000 个序列上下文时， AUC 增加到 0.95 。当使用 300,000 个序列上下文时，性能达到平台期。换句话说，基因组中 1.45%的 21-nt 序列上下文(即 300000 /20.7

模型性能的影响。由于限制性内切酶识别位点长度为 6 b p ，且不包含 C p G 位点， 因此切割位点周围最近的两个 C p G 位点之间的最小核苷酸数限制为 6 个碱基( 不包括 C p G 位点内的 4 个碱

基)(S I 附录， 图 s 9 a 和 B)。两个最近的 C p G 位点之间的最大核 苷酸数为 1 7 个( 即 21 − 4) ， 因为该评估考虑了 2 1 nt 的窗口大小。对于人类( 甲基)- 小鼠(u n m et h )数据集，这两个最近的 C p G 位点

中有 8 2 .4 % 具有“ M -U ” 模式( 图 3 C ) ，这表明来自人类部分的 C p G 位点内的第一个胞嘧啶被甲基化( M ) ，而 C p G 位点内的第二个胞嘧啶未甲基化(U) 。这些结果表明 HK 模型甚至可以稳健地解码 DNA 分子中每个 C p G 位点的甲基化

百万× 100%)足以训练 HK 模型很好地区分 CpG 位点的甲基化和非甲基化胞嘧啶。我们推测，许多序列上下文可能对 DNA 聚合酶的动力学特征有类似的影响

具有不同的甲基化状态。在人类(unmeth) -小鼠(meth)数据集中，这两个最近的 CpG 位点中有 82.0%具有“U-M”模式，这一事实进一步证明了这一结论(图 3D)。C

除了全基因组的甲基化水平，我们进一步分析了 1 兆碱基(Mb)分辨率下的甲基化水平。Circos 图(19)显示了对褐皮 DNA、胎 盘 DNA 和 HepG2 细胞系 DNA 样本的分析(图 5 A-C)，通过 HK 模型(图 5 A-C，内环)推断出的 1 mb 基因组箱的甲基化水 平谱与 BS-seq(图 5 A-C，外环)确定的甲基化水平谱高度一致。 HK 模型与 BS-seq 之间的一致性在散点图( 图 5 D-F)中得到进一 步证明，黄皮毛 DNA、胎盘 DNA 和 HepG2 细胞系 DNA 样本 的相关系数分别为 0.85、0.94 和 0.98。HCC 肿瘤样本及其配对

的邻近非肿瘤组织样本的结果见 SI 附录，图。图 S10 和图 S11。

众所周知，在转录起始位点(tss)附近的区域会观察到较低的甲基化密度(12)。值得注意的是，在 HK 模型确定的结果中确实看到了 TSS 区域周围甲基化水平的“谷型”，这在 BS-seq结果中得到了证实(图 5G)。

图 3所示。 人-鼠杂交片段甲基化模式分析。( A)人(冰毒)-鼠( 冰毒) 数据集中存在的人-鼠杂交片段中 CpG 位点的甲基化水平。 CpG 位点根据与限制性内切位点( HindIII 或 N coI)最近碱基的相对距离汇集在一起。(B)人(u n meth )-鼠( met h)数据集中存在的人-鼠杂交片段中 CpG 位点的甲基化水平。( C)在人(冰毒)- 鼠( 冰毒) 数据集中存在的人- 鼠杂交片段的两个最近的 Cp G 位点的甲基化模式，紧邻一个限 制性切切位点 (Hind III 或 N coI)。(D ) 人类 (un meth ) -小鼠( met h)数据集中存在的人类- 小鼠杂交片段的两个 CpG 位点的甲基化模式， 紧邻限制性切切位点(Hin dIII/N coI)。“ M- M” 表示人类和小鼠部位的第一和第二 CpG 位点都被甲基化。“ M-U ” 表示人类部分的第一个 Cp G 位点甲基化， 而小鼠部分的第二个 Cp G 位点未甲基化。“ U- M” 表示人类部分的第一个 CpG 位点未甲基化， 而小鼠部分的第二个 Cp G 位点已甲基化。“ U -U ” 表示人类和小鼠部位的第一和第二 Cp G 位点都未甲基化。

使用 HK 模型测定生物样品的甲基化。

为了进一步验证训练后的 HK 模型是否可以用于分析真实的生物样本，我们使用 S eq u el II 测序仪和 S eq u el II 测序试剂盒 1 .0 (P acBi o )对 11 个组织 DN A 样本进行了测序(S I 附录，表 S 2 )。我们获得了 600 万个已测序分子的中位数， 大小中位数为 5 .9 千碱基(k b )。子读长中位数为 4 .3 × (IQR :3 .6 ~ 6 .7 ×)。每个样本也通过 BS -s eq 测序， 中位数为 5000 万对 reads 。通过 M et h y -Pi p e 软件确定 C p G 位点的甲基化状态(1 4 )。

我们比较了 HK 模型和 BS -s eq 两种测量结果之间的总体甲 基化水平。总体甲基化水平定义为在所有测序的 C p G 位点中确定被甲基化的 C p G 位点的百分比。图 4 显示，通过 HK 模型分析的样品的总体甲基化水平与 BS -s eq 量化的结果相关性良好(r =0 .9 9;P 值< 0 .0 0 0 1 )。胎盘 DN A( 样本 P L 0 1 ) 、肝细胞癌( HC C )肿瘤组织 D N A (HC C 0 1 和 HC C 0 2 )和 Hep G2 细胞系 DNA 的甲基化水平( 范围: 4 8 .4 % 至 5 8 .4 %) 低于邻 近非 肿瘤 DN A ( NT 0 1 和 NT 0 2 )和黄皮毛 DN A (B C 0 1 至 B C 0 5 )(范围: 6 9 .0 % 至 7 5 .7 %) 。在胎盘 DN A 、H C C 肿瘤组织 D NA 和 H ep G2 中观察到低甲基化细胞系 DN A 与先前的研究一致(1 5 -18 )，进一步表明 HK 模型在区分来自各种生物样品的天然 DNA 分子中的甲基化和非甲基化胞嘧啶方面具有稳健性。

HK 模型与 BS-Seq 单碱基分辨率甲基化相关性研究

。为了比较单碱基分辨率下的相关性，我们计算了样本 BC01 的 SMRT- seq 和 BS- seq 结果中至少 20 个测序分子覆盖的每个 CpG 位点的甲基化水平。由于存在大量的 CpG 位点，因此使用平滑散点图来可视化 HK 模型和 BS-seq 推断的甲基化水平的相关性(SI 附录，图 S12)。HK 模型与 BS-seq 的 Pearson 相关系数为 0.8 (P 值< 0.0001)。选取具有代表性的区域(ch r1 : 1 4 5 ,0 7 1 ,3 6 9 ~ 1 4 5 ,0 7 5 ,7 00 ) ， 进行单碱基分辨率下 HK 模型与 BS -s eq 的比较。 如图 6A 所示， 16 个测序分子来自该区域，通过 HK 模型进行分析，中位读长为 3103 nt(范围为 14 8 4 ~8490 nt)。分子与 Cp G 岛重叠的部分(C G I)

图 4 所示。通过 BS-seq 和 HK 模型量化的总体甲基化水平的相关性。每个点代表一个样本。

图 5 所示 。 BS-seq 和 HK 模型在 1- Mb 分辨率下定量的甲基化水平。Circo s 图显示了由 HK 模型(内环) 和 BS-s eq(外环)确定的人类基因组中不同 1- Mb 区域的甲基化水平，包括白毛(A)、胎盘(B)和 Hep G2 H CC 细胞系( C)。散点图显示了由 HK 模型和 BS-seq 确定的每个 1- Mb 基因组区域的甲基化水平的相关性，包括白毛(D)、胎盘(E)和 HepG2 HCC 细胞系(F)。(G)围绕 tss 的甲基化模式。该区域主要被确定为非甲基化区域 (中位数:163 nt;范围:30 - 599 nt)。图 6B 说明了这一点

CGI 区域以外的部分分子(即 CGI 岸)倾向于甲基化(图 6A)。这种明显的模式在 BS-seq 结果中得到证实，有 102 个结果序列HK 模型可以提供完整的基因型信息，包括 A、C、G 和 T(即四字母信息)和 CpG 二核苷酸的甲基化状态。然而，对于 BS- seq，基因型而与“T”等位基因相连的片段则未甲基化。等位基因之间的差异甲基化模式通常在非印记区域( 例如从基因组中随机选择的区域[chr12: 21,729,541 ~ 21,739,542]) 中观察不到(图 7B)。与非印记区域相比，在两个等位基因之间，所有 4 个印记基因均有差异甲基化区域(图 7C)。SI 附录，图 S13 显示了覆盖其他三个印迹基因(NAP1L5、ZIM2 和 PLAGL1)的每个 DNA 分子的甲基化模式，它们都表现出等位基因特异性甲基化模式。

讨论

我们开发了一种通过 SMRT 测序全面利用动力学信号和序列背景来实现全基因组胞嘧啶甲基化检测的方法。的

图 6 所示。 单碱基分辨率下的甲基化 模式。 (A ) CG I 重叠区 域 chr1: 145,07 1,3 69 ~ 145 ,075 ,70 0 的甲基化模式。 CGI 的基因组坐标以蓝色突出显示。“ (I )” 和“( II )” 表示两个序列读数， 用于突出 HK 模型和 BS-s eq之间读数的差异。(B)使用 HK 模型(表示“ I”) 和 BS-s eq (表示“ II” )生成的遗传和表观遗传信息。 为了便于可视化， A、 C、 T 和 G 用不同的颜色表示。对于 HK 模型， 直接从结果中同时读出原始 基因组序列和甲基化信息。对于 BS -seq，“ TG ” 读数的解释( 即 T 是否意味着未甲基化的胞嘧啶， 或者 T 是否存在于基因组的该位置 )只能在与参考基因组序列进行比较后才能进行。填充的棒棒糖， 甲基化的 C;未 填充棒棒糖， 未甲基化 C。

信息主要局限于三个字母的信息(即 A、G 和 T)。

代表性印迹基因的甲基化测定。DNA 甲基化对于在父本或母本等位基因上建立印记标记很重要(20)，通常显示等位基因特异性甲基化模式。因此，我们期望 SMRT 测序能够使用 HK 模型在单分子分辨率下分析等位基因特异性甲基化模式。我们选择了四个具有代表性的印迹基因，SNURF、PLAGL1 、NAP1L5 和 ZIM  2，这些基因在一项研究中被报道在各种组织中普遍印迹(21)。由于样品 BC01 具有较高的测序深度(SI 附录，表 S2)，我们使用 HK 模型来确定样品 BC01 中与这四个印迹基因重叠的分子的甲基化状态。例如，印迹基因 SNURF，显示了等位基因特异性甲基化模式，覆盖了 15 号染色体上25,200,004 ~ 25,201,976 的已知印迹控制区(22)(图 7A)。在印记控制区，与“C”等位基因相连的片段被甲基化

图 7 所示。 每个单分子的甲基化模式来源于印迹区域。 (A )显示与 SNU RF 基因印迹区相关的每个 DNA 分子甲基化模式的一个。 x 轴表示 CpG 位点 的坐标。以蓝色高亮显示的坐标表示 cgi。 红点表示甲基化的 Cp G 位点。 绿点表示未甲基化的 Cp G 位点。 嵌入在每个水平序列红色和绿色点之间 的字母(即 Cp G 位点)表示 SNP 位点的等位基因。每个水平点序列右侧括 号中的数字表示片段的大小。虚线矩形表示与已知印迹控制区重叠的区域。 (B)一个前任样品显示了源自非印迹区域的每个 DNA 分子的甲基化模式。 虚线矩形表示突出显示用于比较的 SNP 位点周围的区域。(C)印迹区和非印迹区之间的甲基化水平。