基于PacBio的甲基化检测方法
写在前面的
上一篇推送讲了cfDNA领域的明星团队卢煜明院士课题组建立了基于PacBio三代测序平台的长链cfDNA分子检测方法。该团队使用对应方法不仅发现了人体内确实存在>200bp甚至是>3kb的长链cfDNA分子，同时还使用了相同的方法比较揭示了HCC（肝癌）患者的长链cfDNA分子比例低于非肿瘤人群的长链cfDNA分子。我在上一篇的推送中只是讲了一半内容。该团队除了建立了PacBio-long-cfDNA检测技术，同时还建立了基于PacBio平台的5mC检测技术（胞嘧啶5号位碳原子甲基化）。本次的推送将详细分享该团队发表的如何利用PacBio测序平台检测5mC。
基于PacBio的5甲基腺嘌呤5mA检测技术
在介绍该团队的PacBio-5mC检测技术之前，首先需要介绍PacBio-6mA（腺嘌呤6号位氮原子甲基化修饰）检测技术。尽管5mC的检测技术并没有依赖6mA检测技术为基础，但是先介绍6mA检测技术主要目的是说明利用PacBio平台检测核苷酸甲基化修饰是可行的。介绍6mA检测技术之前，我们先了解一下PacBio测序中碱基识别的基本原理。
PacBio基于边合成边测序的原理（这种原理和illumina NGS测序原理一样），每一个碱基会修饰荧光基团，在合成时，游离的dNTP会被固定在底板上的酶捕获，激发光从玻璃底板发出，由于PacBio测序仪中测序小孔的直径很小，激发光的穿透力逐渐衰减，只能在小孔中传输很短的距离。所以只有dNTP足够靠近底部，荧光基团才会被激发光照射到发出荧光，其他游离的dNTP只有极小的可能性飘到小孔底部被照射到。在一个碱基合成结束后，带有荧光基团的磷酸基团从dNTP上掉落，发生淬灭。其他dNTP继续合成检测。因此，在PacBio测序中每一个碱基边合成边测序过程中均会产生由弱至强再由强至弱的荧光信号，那么PacBio仪器中存在一个接收实时接收荧光信号的传感器。这个传感器会记录不同时间点采集到的荧光信号类别及强弱，传感器进一步把这些连续光信号转化为数字信号，进而生成以时间-荧光信号强弱构成的脉冲曲线。进而可以根据脉冲曲线的特征识别碱基，因此PacBio测序还有额外的功能监测聚合酶合成碱基速度。
在前人的研究中，科学家发现当腺嘌呤的6号位氮原子发生甲基化修饰时，PacBio测序时荧光信号的脉冲曲线可以检出这种异常。具体的检测方法是，如果合成的模板链上存在一个6mA，那么合成的dTTP掺入合成的互补链时会显著地被延迟。这是因为6mA显著地影响了DNA合成中DNA聚合酶的动力学性能。下图中A和B分别展示了当模板链的相同位点存在甲基化或非甲基化的腺嘌呤时，对应的荧光脉冲曲线的曲线特征。因此PacBio检测6mA时只需要判断A的临近下文碱基掺入时间是否显著慢于理论平均时间即可。
基于PacBio的5mC检测技术
与6mA检测不同，5mC检测显得非常困难。因为5mC对于DNA合成时DNA聚合酶的动力学性质影响较轻微。难以产生如此显著的脉冲曲线差异。如果要设法找到5mC的检测方法，首先需要找到可以作为潜在区分5mC和非5mC的特征变量。但是目前已知的信息只有：
1.PacBio可以实时准确识别碱基类型（A/T/C/G），也就是PacBio可以知晓碱基及其合成时间。
2.PacBio在边合成边测序过程中会产生随着时间变化的不同碱基的荧光脉冲曲线。
如何根据如上两条信息挖掘出更多的特征变量？作者所在的团队通过思考找到了多个维度的可以描述PacBio胞嘧啶碱基及其临近环境的变量。
1.因为PacBio可以准确识别碱基类型，所以任意一个C，我们总能知晓其上下文碱基。
2.因为PacBio可以产生随着时间变化的不同碱基的荧光脉冲曲线，根据曲线的物理学定义，我们可以用任意两个相邻曲线峰的间隔描述相邻碱基参入合成的DNA链的时间差（Interpulse duration, IPD），可以用任意一个曲线的起峰和收峰的差异表征一个碱基参入合成的DNA链的时间差（Pulse width, PW）。
根据作者如上两点思考，对于任意一个碱基C，我们都能将以C碱基为中心，上下文碱基序列作为一个特征变量（Sequence Context），也能将PacBio上下文碱基掺入DNA链的时间差（IPD）作为第二个特征变量，还能将上下文碱基自身掺入DNA链的时间（PW）作为第三个特征变量。紧接着作者将上下文碱基序列及ATCG碱基联立为一个表，表格的列代表每一个C碱基及其上下文序列，行代表ATCG，表格中上下文序列中任意一个碱基在其对应的行中填入IPD及PW值。这样一个4x(2N)的表格被生成，由于双链DNA分为Watson链和Crick链，因此对于一段序列会有两个互补的表格。通过将Watson链及Crick链所对应的上下文序列的表格合并，最终产生了(8x2N)的表格。其中N的大小定义为窗口大小，作者选择窗口大小为21，即将待检测的C碱基及其上下游10bp序列纳入一个窗口，生成8x42的表格。表格中除了已填入的IPD和PW值以外，其余坐标均填入0。那么一个稀疏矩阵便生成。作者将8x42表格作为输入层，判断是否为甲基化或未甲基化的胞嘧啶作为输出层，作者建立了基于卷积神经网络（Convolutional neural network, CNN）的机器学习模型。由于输出层判断是否为甲基化碱基，整个模型属于二分类模型。整个CNN模型的卷积层共由两个1维卷积层组成，每个卷积层有64个滤波器，内核大小为4，使用ReLU作为连接上层神经元和下层神经元之间的激活函数。我推测，作者在这里选择ReLU激活函数的原因可能是该函数会使得一部分神经元的输出为0，这样造成网络具有稀疏性，并且还能减少参数的相互依存关系，缓解过拟合产生。如果使用sigmoid函数作为激活函数，计算量相对较大。因为sigmoid函数要进行浮点数四则运算，ReLU函数的表达式为f(x) = max(0, x)只需要一个if-else语句即可完成。不过由于本人目前暂无能力从全局理解该团队当时在这个层面的考量，因此这部分推测仅仅作为本人暂时的理解。
作者同时还对这种利用了上下文序列及IPD和PW等来自于PacBio测序过程中整体的动力学特征（Holistic Kinetic）的变量构建的模型取了一个非常有趣的名称，叫HK model，双关了作者团队所在的地点。
根据机器学习的一般步骤，策略模型选择好后，需要利用训练数据对模型进行训练。作者首先准备了使用未甲基化的dNTP扩增的全基因组（WGA）测序数据集作为阴性测试数据集，及使用了M.SssI酶处理的含有高度CpG甲基化的测序数据集作为阳性测试数据集。其中阴性测试数据集中大部分位点都没有甲基化，唯一可能存在甲基化信号的位点全部来自于本底基因组中的甲基化位点。M.SssI是来自于大肠杆菌的一个菌株，该菌株含有来自Sprioplasma sp. MQ1的甲基转移酶基因。该酶可以使得双链DNA中所有的CpG位点甲基化。作者选择了M.SssI酶处理的阳性数据集的一般数据作为阳性训练样本，并选择了阴性测试数据中等量的CpG位点的数据作为训练样本和阳性训练样本一同用于HK模型训练。阴阳性数据集剩余的一般样本用于模型的评估测试。作者还利用了PacBio的Sequel I、II、III三种测序试剂盒生成相应的测序数据进行模型训练。经过分析，作者发现HK模型能够很好地区分由不同试剂盒生成的测试数据的甲基化及未甲基化的胞嘧啶。通过AUC曲线分析，作者确定了一个用于区分甲基化或未甲基化胞嘧啶的临界cut-off值，这个值选择为0.5。在这个过程中作者还使用了另外一种隐马尔可夫模型（HMM）作为评估其中一个高测序深度的BC01样本的5mC检测性能。作者将这部分工作放于了补充材料中，由于本人暂未理解为什么对BC01样本要进行额外的分析研究，不过通过作者的论文描述看，隐马尔可夫模型对BC01样本的甲基化检测性能（83%灵敏度+84%特异性）也低于基于CNN的HK模型（87%灵敏度+92%特异性）。除此以外，作者还通过评估窗口大小及上下文序列长度甚至是测序深度变化对HK模型检测甲基化的影响。
HK模型对HCC肿瘤患者及非肿瘤群体甲基化差异的对比
作者建立好HK模型后在他们的第二项研究中进一步对肝癌患者及非肿瘤群体的甲基化程度进行了对比。从下图A提示，肿瘤患者的甲基化水平低于非肿瘤群体。紧接着，作者根据长链cfDNA分子的单分子CpG位点甲基化水平构建了HCC甲基化评估分数，用于评估HCC的风险。由于长链cfDNA分子长度分布差异会导致单个长链cfDNA分子中CpG位点存在数量上差异。上一篇推文中已经介绍了作者发现长度在1kb以内的长链cfDNA分子大部分含有6个以内的CpG位点。而长度在1kb以上的长链cfDNA分子中大部分含有7个以上的CpG位点。因此在构建肝癌甲基化分数前，作者将cfDNA分子分为了两组，一组是只含有1~6个CpG位点的cfDNA分子，一组是含有7个以上的CpG位点的cfDNA分子。通过作者构建的HCC甲基化分数计算，它们都能区分健康人和HCC患者。但是AUC分析揭示，使用至少7个CpG位点的长链cfDNA分子具有更好的HCC风险评估能力。
总结
作者不仅开发了基于PacBio平台的甲基化检测技术，还将其应用于HCC肿瘤患者和非肿瘤群体的cfDNA甲基化差异比较研究中。并揭示了使用至少7个CpG位点的长链cfDNA分子的CpG甲基化水平构建的HCC甲基化分数具有较好评估HCC风险的能力。为肿瘤早筛提供了更新的思路。
参考
Tse OYO, Jiang P, Cheng SH, Peng W, Shang H, Wong J, Chan SL, Poon LCY, Leung TY, Chan KCA, Chiu RWK, Lo YMD. Genome-wide detection of cytosine methylation by single molecule real-time sequencing. Proc Natl Acad Sci U S A. 2021 Feb 2;118(5):e2019768118. doi: 10.1073/pnas.2019768118. PMID: 33495335; PMCID: PMC7865158.
Choy LYL, Peng W, Jiang P, Cheng SH, Yu SCY, Shang H, Olivia Tse OY, Wong J, Wong VWS, Wong GLH, Lam WKJ, Chan SL, Chiu RWK, Chan KCA, Lo YMD. Single-Molecule Sequencing Enables Long Cell-Free DNA Detection and Direct Methylation Analysis for Cancer Patients. Clin Chem. 2022 Sep 1;68(9):1151-1163. doi: 10.1093/clinchem/hvac086. PMID: 35587130.
Flusberg BA, Webster DR, Lee JH, Travers KJ, Olivares EC, Clark TA, Korlach J, Turner SW. Direct detection of DNA methylation during single-molecule, real-time sequencing. Nat Methods. 2010 Jun;7(6):461-5. doi: 10.1038/nmeth.1459. Epub 2010 May 9. PMID: 20453866; PMCID: PMC2879396.
写在后面的
这一次推送的文献相对于以前的文献更加不好理解，主要是因为一方面利用机器学习构建PacBio平台的5mC甲基化检测技术本身具有不可解释性。什么是不可解释性？这是机器学习人工智能所有模型的一个特征。尽管机器学习人工智能有相关的算法，但是算法本身只是起到了局部拆分分析数据的能力。而全局每一个机器学习模型是如何影响整体模型的演化的这是不可解释的，就好比人类自己也不知道自己的大脑是如何理解传播信息的。我们只知道我们的大脑每一个脑区，每一个神经元在接触信息时产生了反应。但是这些反应是怎么整合加工的，我们无法理解。机器学习的本质就是在模仿人脑的学习模式。所以不可解释性就是机器学习模型的一大特征。就像这次推送中一样，HK模型就是基于CNN的机器学习模型。我们知道HK模型利用了什么特征变量，可以完成对三代数据中5mC的判断，但是HK模型是如何根据这三个特征做出的5mC判断的具有实际意义的过程是不可理解的。打一个比方，三体小说中，面壁者之所以不好被对付，就是因为人与人之间是无法知道对方内心真正的想法的，我们只能通过一个人的输入输出来判断这个的性格特点及好坏。