基于岭回归矫正MB-seq甲基化水平的方法及系统的制作方法
基于岭回归矫正MB-seq甲基化水平的方法及系统的制作方法
【技术领域】
[0001] 本发明属于基因工程技术领域，尤其涉及一种基于数学模型-岭回归矫正MB-seq 甲基化水平的方法及系统。
【背景技术】
[0002] DNA甲基化（DNAmethylation)是最早发现的修饰途径之一，大量研宄表明，DNA 甲基化能引起染色质结构、DNA构象、DNA稳定性及DNA与蛋白质相互作用方式的改变，从 而控制基因表达。早在1942年，C.H.Waddinton就提出了表观遗传学的概念，他指出，表观 遗传与遗传相对，主要研宄基因型和表型的关系。而现在，对于表观遗传学，比较统一的认 识是，其研宄在没有细胞核DNA序列改变的情况时，基因功能的可逆的可遗传的改变。也就 是说，在不改变基因组序列的前提下，通过DNA和组蛋白的修饰等来调控基因表达，其中又 以DNA甲基化（DNAmethylation)最为常见，DNA甲基化成为表观遗传学的重要组成部分。 随着人类基因组计划的开展，科学家们开始在基因组水平来研宄表观遗传学，逐步形成表 观基因组学（epigenomics)。表观基因组学就是要在整个基因组水平来研宄表观遗传过程 以及与这些过程密切相关的特定基因组区域的识别与鉴定。2000年10月，人类表观基因 组协会（HumanEpigenomeConsortium)由欧盟赞助，启动了旨在于人类6号染色体MHC区 域首先做出DNA的甲基化图谱的先导计划（PilotProject)。该计划顺利完成，引导启动 了 2003年的人类表观基因组计划（HumanEpigenomeProject，HEP)。2005年，美国国家卫 生院（NIH)下属的国立癌症研宄所启动了癌症基因组先导计划。2006年，该所与国立人类 基因组研宄所一起共同启动癌症基因组计划（CancerGenomeProject)。表观基因组学和 DNA甲基化与癌症的研宄成为新的热点。
[0003] 目前，人们认识到DNA甲基化对基因组正常功能维持是必要的，表观遗传水平的 改变可看作是复杂疾病（如癌症，精神疾病）发病机理的第一步（病因学上），以表观遗传 学为基础的药物对治疗复杂疾病将有巨大的潜力。但是，只有当我们的检测技术能够勾勒 出全基因组的DNA甲基化图谱时，对于DNA甲基化这一表观遗传修饰改变有关的疾病才有 可能得到全面的认识。因此，全基因组DNA甲基化检测技术的发展显得尤其重要，它的发展 是表观遗传学、表观基因组学研宄的重要基础，也将会给当今分子遗传学研宄带来新的变 革。
[0004] 而目前现有的检测技术，还是受到成本、分析周期、基因组覆盖度、分辨率及技术 可操作性等因素的影响：
[0005] (1)对于DNA甲基化检测技术领域，早期的检测手段是KristenH等结合454测序 仪对DNA甲基化进行了检测。该方法首先将基因组DNA进行亚硫酸盐处理，用含有公共接 头的引物对目的片段逐个扩增。利用454新一代测序仪对超过40例样本的25个相关基因 的CpG富集区靶点片段同时测序，共产生了 294631个序列。该方法首次利用高通量测序仪 精确定量检测了目的片段单个CpG位点的甲基化状态，初步展示了新一代测序仪的特点。
[0006] (2)ShawnJ等利用全基因组亚硫酸盐处理测序（WholeGenomeBisulfite Sequencing)，简称WGBS或MethylC-seq，对拟南芥全基因组的DNA甲基化谱进行高通量测 序。他们通过对结果进行比较分析后发现，无论甲基化位点的密度和序列所含碱基如何， 芯片技术发现甲基化位点的能力远不及亚硫酸氢盐测序法，后者甚至还在结构相对简单、 富含转座子的区域发现了甲基化位点，而芯片技术由于交联反应的限制，很难发现这类甲 基化位点。因此，MethylC-seq被视为DNA甲基化组检测的金标准，它能够实现对被测物种 或样本的全基因组DNA甲基化谱的全面、深度、单碱基分辨率的检测。但是，MethylC-seq 需要获得被测物种至少30倍覆盖度的测序数据量。以人类为例，需要至少90Gb的测序数 据，用目前Illumina最新的TrueseqV4的测序试剂价格来衡量，需要约4万人民币；并且 90Gb的亚硫酸盐处理后测序数据需耗费较长的运算时间，目前分析亚硫酸盐处理后测序数 据的常用软件BSMAP，在8核24GB内存的情况下，需要约22天的时间方能得到最终甲基化 图谱。另外，在大部分哺乳动物和植物中，甲基化的胞嘧啶（5mC)主要发生在胞嘧啶一鸟嘌 呤二核苷酸（CpG)上，约只占了全基因组所有碱基数量的1-6%，这使得MethylC-seq所得 到的数据中，仅有20 - 30%数据是有效地提供了DNA甲基化的信息。因此，高额的成本和 费时的运算，大大限制了MethylC-seq的大规模推广应用，尤其是在进行大型基因组的物 种或多样本的DNA甲基化图谱比较研宄中。
[0007] 然而，不同的甲基化DNA富集方法为基于新一代测序技术的DNA甲基化图谱检测 的成本控制奠定了基础。这些富集方法主要包括免疫共沉淀和限制性内切酶等。免疫共沉 淀方法特异性高，而结合限制性酶的新一代测序方法将具有较高的灵敏度。虽然每种方法 都有一定的局限性，但是它们为在基因组范围选择功能区来研宄DNA甲基化图谱提供了更 多的选择。
[0008] (3)基于限制性酶切的测序方法
[0009] 限制性内切酶也多应用于甲基化相关基因的鉴定，是表观遗传学研宄的重要工具 之一。已有3种限制性酶在甲基化研宄中得到应用：第一，甲基化敏感性内切酶，如BstUI、 HpaII、HhaI、SmaI以及Notl。这些酶能够识别CG富集区中的非甲基化位点，而甲基化位点 因为受到甲基的保护而不被识别。第二，甲基化依赖型限制性酶，它能识别并酶切甲基化的 CG位点。McrBC是最有代表性且最常使用的甲基化依赖型限制性酶。具体而言，McrBC可 以识别两个甲基化的半位点（RmC，R=A或G)且这两个位点存在于40~3000bp之内，酶 切发生在两个位点之间。第三，CG甲基化不敏感的同裂酶。例如MspI是Hpall的同裂酶， 他们的识别位点相同但前者不受识别位点甲基化状态的影响。与MspI相似，Xmal是Smal 的同裂酶，它们都识别'CCCGGG'位点，但是Smal酶切后生成钝端而Xmal生成5'粘性末 端。基于此，通过以上几种甲基化敏感或不敏感的限制性内切酶，对被测基因组进行酶切 后，能够特异性的富集高CpG密度的区域，再对其进行二代测序文库的制备以及高通量测 序，从而能够得到单碱基分辨率的DNA甲基化图谱。此类方法中具有代表性的是RRBS和 MRE-seq。针对人类物种而言，RRBS利用MspI对基因组DNA酶切后，然后选取40_220bp的 区域的DNA片段进行亚硫酸盐处理以及高通量测序，仅需大约3Gb的测序数据，便能够覆盖 CpG岛（CGI)的40%左右和启动子区域20%左右的CpG位点，并且由于CGI和启动子区域 在基因表达调控中的重要性，RRBS得到大规模的推广应用。同样，MRE-seq通过甲基化敏感 性内切酶（Bstn、Hpall、Smal)对基因组DNA进行酶切后，进行文库制备和高通量测序，大 约能涵盖人类基因组6%的CpG位点，只需3Gb的总数据量，便能够达到饱和。但目前大量 有关于肿瘤发生过程中DNA甲基化组改变的研宄中发现，差异甲基化区域（Differential methylationregion，DMR)主要发生在CGI的侧翼区域（CGIshore)，并且这些区域对于基 因表达的调控更为明显；同时重复序列元件（Repeatselement，RE) -般呈现高度甲基化， 这些元件上的甲基化状态与基因组的稳定性密切相关。而目前基于限制性内切酶的DNA甲 基化组测序技术，均不能很好的揭示CGIshore和RE区域上的DNA甲基化状态，因此它们得 到的DNA甲基化组不能够代表真正意义的全基因组DNA甲基化图谱。
[0010] (5)基于免疫共沉淀甲基化DNA片段的测序方法
[0011] 哺乳动物MBD家族由五个成员组成，包括MeCP2、MBD1、MBD2、MBD3和MBD4,甲基 化的CpG二核苷酸可被MBD特异性识别并结合。最近，一种基于重组抗体样蛋白MBD的结 合免疫沉淀和高通量测序技术的方法被应用于基因组DNA甲基化图谱的研宄，这种方法被 称为MBD-seq。另外，通过5-甲基胞嘧啶抗体也可用来进行富集甲基化DNA片段，而后结 合高通量测序技术，被称为甲基化DNA免疫共沉淀测序（MeDIP-seq)。基于免疫共沉淀原 理的DNA甲基化组检测技术，通过对甲基化修饰的DNA片段进行特异性富集，撇弃了非甲基 化修饰的DNA片段，对后续的数据产出量提供了良好的成本控制基础。以人类基因组为例， MeDIP-seq只需要大约25million的reads数量，便可涵盖80%的CpGs，大幅度的降低了测 序成本和分析周期。该方法具有较高的特异性，但是也有一定局限性一一它 们没有结合亚 硫酸盐处理或者甲基化敏感的限制性内切酶酶切，无法得到单碱基分辨率的DNA甲基化 图谱，其分辨率大约为l〇〇bp。
[0012] (6)针对免疫共沉淀甲基化DNA测序数据的分析方法
[0013] 由于MethylC-seq的高成本、RRBS和MRE-seq的基因组低覆盖度、MeDIP-seq 和MBD-seq的低分辨率，促使研宄人员利用生物信息算法去更好维持成本、基因组覆盖度 和分辨率之间的平衡，从而使得高通量测序更好服务于DNA甲基化组的研宄。基于免疫 共沉淀甲基化DNA的高通量测序数据，人们开发了一系列生物信息学算法，具体如下。针 对MBD-seq和MeDIP-seq的测序数据，MEDME和BayMeth分别被开发了出来，它们能够将 测序数据中得到的reads数量转换成窗口大小为100bp的区域甲基化水平（Riebleret al. 2014) ;MEDIPS可以将MelDP-seq数据进行计算，从而实现单碱基分辨率的DNA甲基化 图谱，但其得到的并非C的甲基化水平，而是介于1 一 1000的MEDIPS值，导致其得到的结 果无法和其他单碱基分辨率的DNA甲基化检测技术得到的结果相互比较；Batman是另一 款针对MeDIP-seq数据的生物信息学算法，其可实现单个CpG的DNA甲基化水平预测，但 所耗费的计算周期较长（Bock2012)，且算法较为复杂，开发人员也未提供完整的代码安 装文件，使得其他研宄者无法很好的重复该算法，另外batman所得到CpG甲基化水平往往 较真实甲基化水平偏低（Riebleretal. 2014);近来，一种基于条件随机场算法的机器自 动学习工具被开发出来，用于MeDIP-seq和MRE-seq整合后数据的单碱基分辨率DNA甲基 化图谱预测，但其应用于人胚胎干细胞系H1所得到的结果，通过与MethylC-seq测序所得 到的结果进行相关性分析发现，pearson系数仅达到0.77(Stevensetal. 2013)，并且由 于MethylCRF未考虑拷贝数变异对于DNA甲基化水平预测结果的影响，使得该方法在应用 于肿瘤发生时，可能会得到更失真的DNA甲基化组（Laird2010;Robinsonetal.2012; Riebleretal. 2014)〇
[0014] (7)做为本领域所公知的现有技术，MB-seq-甲基化DNA富集结合亚硫酸盐翻转 的甲基化检测技术（MeDIPbisulfitesequencing,MB-seq)拥有许多优点：它是一种高通 量的、单碱基分辨率的、低成本的、可适用于多种已知序列物种的DNA甲基化检测技术，但 是MB-seq存在甲基化水平的偏差，MB-seq甲基化水平被线性放大，所以MB-seq得到的单 个CpG位点的甲基化水平是相对甲基化水平。
【发明内容】
[0015] 针对现有技术中MB-seq-甲基化DNA富集结合亚硫酸盐翻转的甲基化检测技术 存在的甲基化水平的偏差的问题。发明人研发了一种基于岭回归矫正MB-seq甲基化水平 的方法及系统，考虑多种与DNA甲基化水平相关的因素，可以将MB-seq的相对甲基化水平 矫正到全基因组的胞嘧啶位点的绝对甲基化水平。
[0016] 本发明提供的基于岭回归矫正MB-seq甲基化水平的方法，包括以下步骤：
[0017] (1)提取信息
[0018] ⑵建模
[0019] (3)岭回归计算；
[0020] 其中，所述的步骤（1)需要提取的信息有：从参考基因组序列中提取基因组CpG密 度、GC含量和CpG-〇E值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每 个胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个 胞嘧啶的绝对甲基化信息；
[0021] 所述的步骤（2)建模如下：
[0023]其中：
[0024]y:目标函数；为以RRBS高通量测序数据唯一比对结果中提取到的覆盖到的每个 胞嘧啶的绝对甲基化信息；
[0025] x:回归变量矩阵；包括行、列；每行代表每个CpG变量；每列分别为每个变量的 CpG密度、GC含量、CpG-〇E值以及相对甲基化信息；
[0026] 所述的步骤（3)岭回归计算具体是
[0028] 对求导，结果为
[0029] 2XT (Y-XW) -2 入W
[0030] 令其为0,求得的值：
[0032] 输入新的回归变量矩阵X即可获得新Y值，即而获得全基因组的胞嘧啶位点的绝 对甲基化水平。
[0033] 优选的，在所述的步骤（1)中提取信息后，将提取到的信息进行阈值过滤，过滤低 质量碱基和序列，并过滤adapter污染序列。以得到更为精确合理的甲基化水平。
[0034]优选的，在所述的步骤（3)计算之前，采用交叉验证评估模型进行数据训练和测 试：
[0035]a).将预测特征变量和真实的甲基化水平分成训练和测试数据集；随机抽取50% 的CpG位点作为训练数据，剩下的50%作为测试数据；
[0036]b).先使用训练数据训练模型；再计算预测甲基化水平值和RRBS测量的甲基化水 平值之间的相关性系数；这个过程重复N次，N次的平均相关性系数用来表示模型的预测精 度；优选的，N多1000。可以获得更为精确的数据。
[0037] 对于每个基因组元件，单独进行训练和岭回归测试；而对同时位于多个基因组元 件的CpG位点，；取多个预测值的平均值；
[0038] c).甲基化水平的预测是全基因组范围的，并且对于RRBS原本就覆盖的位点，采 取RRBS的观测值作为最终的甲基化水平；所有未被RRBS覆盖的CpG位点，一律认为其未被 甲基化，并且不用于岭回归，甲基化水平预测值小于0或者大于分别基于岭回归的原则规 整到0和1。
[0039] 优选的，在在模型数据训练时：
[0040]a).当变量间存在共线性的时候，通过引入lambda表达式以解决最小二乘回归得 到的系数不稳定，方差很大的问题；
[0041]b).当模型包含常数项时，岭回归函数对y进行中心化，以y的均值作为因子；对 X进行中心化和归一化，以x中各个变量的均值和标准差作为因子；这样对x和y处理后，x 和y的均值为〇,这使得回归平面经过原点，即常数项为〇 ;
[0042]c).当模型不包含常数项时，因为要强制通过原点，该模型假设各个变量的均值为 0,因此不对X和y进行中心化，但是对x进行归一化，而且归一化因子也是假设变量均值为 〇计算出来的该变量的标准差。
[0043] 优选的，在使用该模型进行测试的时候，需要首先对x和y进行中心化和归一化， 此时因子是使用训练模型时候进行中心化和归一化的因子，然后再与系数相乘得到预测结 果。
[0044] 优选的，在步骤（3)岭回归计算之后，进行如下对异常点处理：
[0045] 1)将MB-seq检测深度为0的位点定义为甲基化水平为0;
[0046]2)结合MB-seq甲基化水平的观测值（MB level)，甲基化CpG个数（MB mCG)，MB-seq测序深度（MB depth)，当前CpG侧翼+/_100bp区域的MB-seq检测到的平均 甲基化水平（MB back),这200bp范围的甲基化CpG位点总数（MB mCG)，以及每一个CpG位 点上下游l〇〇bp的基因组CpG密度、GC含量，CpG-〇E值等对甲基化水平检测的影响，利用 岭回归导入到模型中，并且机器学习得到某一胞嘧啶位点甲基化水平；
[0047] 3)将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回归的甲 基化水平值小于〇的位点自动归为甲基化水平为0。
[0048] 优选的，所述的相对甲基化信息包括：MB-seq甲基化水平的观测值MBlevel，甲 基化CpG个数MBmCG,MB-seq测序深度MBdepth,当前CpG侧翼+/_100bp区域的MB-seq 检测到的平均甲基化水平MBback,这200bp范围的甲基化CpG位点总数MBmCG。
[0049] 岭回归矫正MB-seq甲基化水平的系统，包括以下模块：
[0050] 提取模块：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-〇E值；从 MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信 息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信 息；所述的相对甲基化信息包括 ：MB-seq甲基化水平的观测值MBlevel，甲基化CpG个数 MBmCG,MB-seq测序深度MBdepth,当前CpG侧翼+/_100bp区域的MB-seq检测到的平均 甲基化水平MBback,这200bp范围的甲基化CpG位点总数MBmCG ;
[0051] 建模模块：根据基因组信息和甲基化信 息，采用岭回归模型对真实甲基化水平 RRBSlevel和回归参数建立回归模型；
[0052] 回归模块：利用岭回归理论，并依据提取出来的基因组信息和甲基化信息，对基因 组上的胞嘧啶位点进行回归以得到甲基化水平的模块。
[0053] 与现有技术相比较，本发明具有如下有益效果：
[0054] 利用数学模型，对基因组上的RRBS未覆盖的胞嘧啶位点进行岭回归预测，使得甲 基化水平检测的准确度大于95%，从而消除MB-seq的偏差并得到全基因组甲基化图谱。利 用本发明可以从高通量测序MB-seq数据中，精确计算全基因组每一个CpG的甲基化水平。
[0055]
[0056]
【具体实施方式】
[0057] 基于岭回归矫正MB-seq甲基化水平的方法，包括以下步骤：
[0058] (1)提取信息⑵建模⑶岭回归计算；
[0059] 其中，所述的步骤（1)需要提取的信息有：从参考基因组序列中提取基因组CpG密 度、GC含量和CpG-〇E值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每 个胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个 胞嘧啶的绝对甲基化信息；
[0060] 所述的步骤（2)建模如下：
[0062]其中：
[0063]y:目标函数；为以RRBS高通量测序数据唯一比对结果中提取到的覆盖到的每个 胞嘧啶的绝对甲基化信息；
[0064] X:回归变量矩阵；包括行、列；每行代表每个CpG变量；每列分别为每个变量的 CpG密度、GC含量、CpG-〇E值以及相对甲基化信息；
[0065] 所述的步骤（3)岭回归计算具体是
[0067] 对求导，结果为
[0068] 2XT(Y-XW) -2 入W
[0069] 令其为0,求得的值：
[0071] 输入新的回归变量矩阵X即可获得新Y值，即而获得全基因组的胞嘧啶位点的绝 对甲基化水平。在所述的步骤（1)中提取信息后，将提取到的信息进行阈值过滤，过滤低 质量碱基和序列，并过滤adapter污染序列。以得到更为精确合理的甲基化水平。在所述 的步骤（3)计算之前，采用交叉验证评估模型进行数据训练和测试：a).将预测特征变量和 真实的甲基化水平分成训练和测试数据集；随机抽取50%的CpG位点作为训练数据，剩下 的50%作为测试数据；b).先使用训练数据训练模型；再计算预测甲基化水平值和RRBS测 量的甲基化水平值之间的相关性系数；这个过程重复N次，N次的平均相关性系数用来表 示模型的预测精度；优选的，N多1000。可以获得更为精确的数据。对于每个基因组元件， 单独进行训练和岭回归测试；而对同时位于多个基因组元件的CpG位点，；取多个预测值的 平均值；c).甲基化水平的预测是全基因组范围的，并且对于RRBS原本就覆盖的位点，采取 RRBS的观测值作为最终的甲基化水平；所有未被RRBS覆盖的CpG位点，一律认为其未被甲 基化，并且不用于岭回归，甲基化水平预测值小于0或者大于分别基于岭回归的原则规整 到0和1。在在模型数据训练时：a).当变量间存在共线性的时候，通过引入lambda表达 式以解决最小二乘回归得到的系数不稳定，方差很大的问题；b).当模型包含常数项时，岭 回归函数对y进行中心化，以y的均值作为因子；对x进行中心化和归一化，以x中各个变 量的均值和标准差作为因子；这样对x和y处理后，x和y的均值为0,这使得回归平面经 过原点，即常数项为0 ;c).当模型不包含常数项时，因为要强制通过原点，该模型假设各个 变量的均值为〇,因此不对x和y进行中心化，但是对x进行归一化，而且归一化因子也是 假设变量均值为〇计算出来的该变量的标准差。在使用该模型进行测试的时候，需要首先 对x和y进行中心化和归一化，此时因子是使用训练模型时候进行中心化和归一化的因子， 然后再与系数相乘得到预测结果。在步骤（3)岭回归计算之后，进行如下对异常点处理： 1)将MB-seq检测深度为0的位点定义为甲基化水平为0 ;2)结合MB-seq甲基化水平的观 测值（MBlevel),甲基化CpG个数（MBmCG)，MB-seq测序深度（MBd印th)，当前CpG侧翼 +/-100bp区域的MB-seq检测到的平均甲基化水平（MBback),这200bp范围的甲基化CpG 位点总数（MBmCG)，以及每一个CpG位点上下游100bp的基因组CpG密度、GC含量，CpG-〇E 值等对甲基化水平检测的影响，利用岭回归导入到模型中，并且机器学习得到某一胞嘧啶 位点甲基化水平；3)将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回 归的甲基化水平值小于0的位点自动归为甲基化水平为0。所述的相对甲基化信息包括： MB-seq甲基化水平的观测值MBlevel,甲基化CpG个数MBmCG,MB-seq测序深度MBdepth, 当前CpG侧翼+/_100bp区域的MB-seq检测到的平均甲基化水平MBback,这200bp范围 的甲基化CpG位点总数MBmCG。
[0072] 岭回归模型是采用岭回归实现的一种正则化线性回归。当多个预测因子含有非0 系数并且呈现正态分布的时候，岭回归是理想的回归方法，岭回归对每一个预测因子影响 小的模型尤其适用，并且它防止线性回归模型系数由于共线性而导致无法模拟和高变异 性。岭回归对共线预测因子的系数收缩并趋于零，例如，给出k相同的预测因子，都将获 得相同的系数等于1/k的单个因子的回归系数值。因此，岭回归不会使某个因子消失，不 能将某些因子摒弃来获得最优预测数据集。岭回归（2)估计解决回归问题（1)使用￡2惩罚 最小二乘法：
[0073]y=yln+X|3 +e: (1)
[0074]y=(yi，…，yn)T其中是观察表型的向量，ln是一个n维列向量的，y是一种常 见的截距，是nXp矩阵的表示，|3表示回归系数的向量，61是残差的向量和
是残差的误差。
[0076]其中
[0078] 是U-normi(二次方程式）损失函数（即残差平方和），X,7是X向量的第i个 行。
[0080] AH:2-norm坫于|3罚分，入彡0是调优（罚分，正规化，或复杂化）参数，这 些参数通过相对重要性决定经验误差和惩罚调节罚分的强度（即线性收缩）。A值越大， 收缩量越大。A的值依赖于数据，通过数据驱动的方法（交叉验证）进行确定使用。
[0081] CpG密度、GC含量以及CpG-〇E值三者计算方法分别为：
[0082] CpG密度：某一个CpG上下游各100bp范围内CpG个数除于201bp长度得到此CpG 位点的CpG密度；
[0083] GC含量：某一个CpG上下游各100bp范围内C和G总数除于201bp长度得到此CpG 位点的GC含量；
[0084] CpG-〇E值：CpG上下游各100bp范围内CpG个数乘于210bp，然后除于C和G个数 的乘积。
[0085] 岭回归矫正MB-seq甲基化水平的系统，包括以下模块：
[0086] 提取模块：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-〇E值；从 MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信 息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信 息；所述的相对甲基化信息包括 ：MB-seq甲基化水平的观测值MBlevel，甲基化CpG个数 MBmCG,MB-seq测序深度MBdepth,当前CpG侧翼+/_100bp区域的MB-seq检测到的平均 甲基化水平MBback,这200bp范围的甲基化CpG位点总数MBmCG ;
[0087] 建模模块：根据基因组信息和甲基化信息，采用岭回归模型对真实甲基化水平 RRBSlevel和回归参数建立回归模型；
[0088] 回归模块：利用岭回归理论，并依据提取出来的基因组信息和甲基化信息，对基因 组上的胞嘧啶位点进行回归以得到甲基化水平的模块。
[0089] 优选的，提取模块提取信息后，设置阈值过滤模块，还将提取到的信息进行阈值过 滤，过滤低质量碱基和序列，并过滤adapter污染序列。
[0090] 优选的，在回归模块计算之前，还设置数据训练和测试模块，采用交叉验证评估模 型进行数据训练和测试：
[0091] a).将预测特征变量和真实的甲基化水平分成训练和测试数据集；随机抽取50% 的CpG位点作为训练数据，剩下的50%作为测试数据；
[0092] b).先使用训练数据训练模型；再计算预测甲基化水平值和RRBS测量的甲基化水 平值之间的相关性系数；这个过程重复N次，N次的平均相关性系数用来表示模型的预测精 度；
[0093] 对于每个基因组元件，单独进行训练和岭回归测试；而对同时位于多个基因组元 件的CpG位点 ，；取多个预测值的平均值；
[0094] c).甲基化水平的预测是全基因组范围的，并且对于RRBS原本就覆盖的位点，采 取RRBS的观测值作为最终的甲基化水平；所有未被RRBS覆盖的CpG位点，一律认为其未被 甲基化，并且不用于岭回归，甲基化水平预测值小于0或者大于分别基于岭回归的原则规 整到0和1。
[0095] 优选的，N彡 1000。
[0096] 数据训练和测试模块对模型数据训练时：
[0097] a).当变量间存在共线性的时候，通过引入lambda表达式以解决最小二乘回归得 到的系数不稳定，方差很大的问题；
[0098] b).当模型包含常数项时，岭回归函数对y进行中心化，以y的均值作为因子；对 X进行中心化和归一化，以x中各个变量的均值和标准差作为因子；这样对x和y处理后，x 和y的均值为〇,这使得回归平面经过原点，即常数项为〇 ;
[0099]c).当模型不包含常数项时，因为要强制通过原点，该模型假设各个变量的均值为 0,因此不对X和y进行中心化，但是对x进行归一化，而且归一化因子也是假设变量均值为 〇计算出来的该变量的标准差。
[0100] 数据训练和测试模块对模型进行测试的时候，需要首先对x和y进行中心化和归 一化，此时因子是使用训练模型时候进行中心化和归一化的因子，然后再与系数相乘得到 预测结果。
[0101] 在回归模块计算之后，设置异常点处理模块，用于对异常点处理：
[0102] 1)将MB-seq检测深度为0的位点定义为甲基化水平为0 ;
[0103] 2)结合MB-seq甲基化水平的观测值（MBlevel)，甲基化CpG个数（MB mCG)，MB-seq测序深度（MBdepth)，当前CpG侧翼+/_100bp区域的MB-seq检测到的平均 甲基化水平（MBback),这200bp范围的甲基化CpG位点总数（MBmCG)，以及每一个CpG位 点上下游l〇〇bp的基因组CpG密度、GC含量，CpG-〇E值等对甲基化水平检测的影响，利用 岭回归导入到模型中，并且机器学习得到某一胞嘧啶位点甲基化水平；
[0104] 3)将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回归的甲 基化水平值小于〇的位点自动归为甲基化水平为0。
【主权项】
1. 基于岭回归矫正MB-Seq甲基化水平的方法，其特征在于，包括以下步骤： (1) 提取信息 (2) 建模 (3) 岭回归计算； 其中，所述的步骤（1)需要提取的信息有：从参考基因组序列中提取基因组CpG密度、 GC含量和CpG-OE值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个 胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞 嘧啶的绝对甲基化信息； 所述的步骤（2)建模如下：其中： y:目标函数；为以RRBS高通量测序数据唯一比对结果中提取到的覆盖到的每个胞嘧 啶的绝对甲基化信息； X:回归变量矩阵；包括行、列；每行代表每个CpG变量；每列分别为每个变量的CpG密 度、GC含量、CpG-OE值以及相对甲基化信息； 所述的相对甲基化信息包括：MB-seq甲基化水平的观测值MB level，甲基化CpG个数 MB mCG, MB-seq测序深度MB depth,当前CpG侧翼+/-IOObp区域的MB-seq检测到的平均 甲基化水平MB back,这200bp范围的甲基化CpG位点总数MB mCG ; 所述的步骤（3)岭回归计算具体是对求导，结果为 2XT (Y-XW) -2 λ W 令其为〇,求得的值：输入新的回归变量矩阵X即可获得新Y值，即而获得全基因组的胞嘧啶位点的绝对甲 基化水平。2. 如权利要求1所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：在所 述的步骤（1)中提取信息后，还将提取到的信息进行阈值过滤，过滤低质量碱基和序列，并 过滤adapter污染序列。3. 如权利要求1所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：所述 的步骤（3)计算之前，采用交叉验证评估模型进行数据训练和测试： a) .将预测特征变量和真实的甲基化水平分成训练和测试数据集；随机抽取50%的 CpG位点作为训练数据，剩下的50%作为测试数据； b) .先使用训练数据训练模型；再计算预测甲基化水平值和RRBS测量的甲基化水平值 之间的相关性系数；这个过程重复N次，N次的平均相关性系数用来表示模型的预测精度； 对于每个基因组元件，单独进行训练和岭回归测试；而对同时位于多个基因组元件的 CpG位点，；取多个预测值的平均值； c).甲基化水平的预测是全基因组范围的，并且对于RRBS原本就覆盖的位点，采取 RRBS的观测值作为最终的甲基化水平；所有未被RRBS覆盖的CpG位点，一律认为其未被甲 基化，并且不用于岭回归，甲基化水平预测值小于O或者大于分别基于岭回归的原则规整 到O和1。4. 如权利要求3所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于： N 彡 1000。5. 如权利要求3所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：在模 型数据训练时： a) .当变量间存在共线性的时候，通过引入lambda表达式以解决最小二乘回归得到的 系数不稳定，方差很大的问题； b) .当模型包含常数项时，岭回归函数对y进行中心化，以y的均值作为因子；对X进 行中心化和归一化，以X中各个变量的均值和标准差作为因子；这样对X和y处理后，X和 y的均值为〇,这使得回归平面经过原点，即常数项为〇 ; c) .当模型不包含常数项时，因为要强制通过原点，该模型假设各个变量的均值为0， 因此不对X和y进行中心化，但是对X进行归一化，而且归一化因子也是假设变量均值为0 计算出来的该变量的标准差。6. 如权利要求3所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：在使 用该模型进行测试的时候，需要首先对X和y进行中心化和归一化，此时因子是使用训练模 型时候进行中心化和归一化的因子，然后再与系数相乘得到预测结果。7. 如权利要求1所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：在步 骤（3)岭回归计算之后，进行如下对异常点处理： 1) 将MB-seq检测深度为0的位点定义为甲基化水平为0 ; 2) 结合MB-seq甲基化水平的观测值（MB level),甲基化CpG个数（MB mCG)，MB-seq测 序深度（MB depth),当前CpG侧翼+/-IOObp区域的MB-seq检测到的平均甲基化水平（MB back)，这200bp范围的甲基化CpG位点总数（MB mCG)，以及每一个CpG位点上下游IOObp 的基因组CpG密度、GC含量，CpG-OE值等对甲基化水平检测的影响，利用岭回归导入到模 型中，并且机器学习得到某一胞嘧啶位点甲基化水平； 3) 将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回归的甲基化 水平值小于〇的位点自动归为甲基化水平为0。8. 岭回归矫正MB-seq甲基化水平的系统，其特征在于：包括以下模块： 提取模块：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-OE值；从MB-seq 高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信息；从 RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信息；所 述的相对甲基化信息包括：MB-seq甲基化水平的观测值MB level,甲基化CpG个数MB mCG, MB-seq测序深度MB depth,当前CpG侧翼+/-IOObp区域的MB-seq检测到的平均甲基 化水平MB back,这200bp范围的甲基化CpG位点总数MB mCG ; 建模模块：根据基因组信息和甲基化信息，采用岭回归模型对真实甲基化水平RRBS level和回归参数建立回归模型； 回归模块：利用岭回归理论，并依据提取出来的基因组信息和甲基化信息，对基因组上 的胞嘧啶位点进行回归以得到甲基化水平的模块。
【专利摘要】基于岭回归矫正MB-seq甲基化水平的方法，属于基因工程技术领域，利用机器学习岭回归理论，并依据RRBS检测出来的绝对MB-seq甲基化水平进行数据训练并建立预测模型，对基因组上的RRBS未覆盖的胞嘧啶位点进行岭回归预测，使得甲基化水平检测的准确度大于95％，从而消除MB-seq的偏差并得到全基因组甲基化图谱。本发明还公开了一种基于岭回归的甲基化水平计算系统。利用本发明可以从高通量测序MB-seq数据中，精确计算全基因组每一个CpG的甲基化水平。
【IPC分类】G06F19/12
【公开号】CN104899474
【申请号】CN201510313520
【发明人】张保荣, 王晓东, 张久文
【申请人】大连三生科技发展有限公司
【公开日】2015年9月9日
【申请日】2015年6月9日
【技术领域】
[0001] 本发明属于基因工程技术领域，尤其涉及一种基于数学模型-岭回归矫正MB-seq 甲基化水平的方法及系统。
【背景技术】
[0002] DNA甲基化（DNAmethylation)是最早发现的修饰途径之一，大量研宄表明，DNA 甲基化能引起染色质结构、DNA构象、DNA稳定性及DNA与蛋白质相互作用方式的改变，从 而控制基因表达。早在1942年，C.H.Waddinton就提出了表观遗传学的概念，他指出，表观 遗传与遗传相对，主要研宄基因型和表型的关系。而现在，对于表观遗传学，比较统一的认 识是，其研宄在没有细胞核DNA序列改变的情况时，基因功能的可逆的可遗传的改变。也就 是说，在不改变基因组序列的前提下，通过DNA和组蛋白的修饰等来调控基因表达，其中又 以DNA甲基化（DNAmethylation)最为常见，DNA甲基化成为表观遗传学的重要组成部分。 随着人类基因组计划的开展，科学家们开始在基因组水平来研宄表观遗传学，逐步形成表 观基因组学（epigenomics)。表观基因组学就是要在整个基因组水平来研宄表观遗传过程 以及与这些过程密切相关的特定基因组区域的识别与鉴定。2000年10月，人类表观基因 组协会（HumanEpigenomeConsortium)由欧盟赞助，启动了旨在于人类6号染色体MHC区 域首先做出DNA的甲基化图谱的先导计划（PilotProject)。该计划顺利完成，引导启动 了 2003年的人类表观基因组计划（HumanEpigenomeProject，HEP)。2005年，美国国家卫 生院（NIH)下属的国立癌症研宄所启动了癌症基因组先导计划。2006年，该所与国立人类 基因组研宄所一起共同启动癌症基因组计划（CancerGenomeProject)。表观基因组学和 DNA甲基化与癌症的研宄成为新的热点。
[0003] 目前，人们认识到DNA甲基化对基因组正常功能维持是必要的，表观遗传水平的 改变可看作是复杂疾病（如癌症，精神疾病）发病机理的第一步（病因学上），以表观遗传 学为基础的药物对治疗复杂疾病将有巨大的潜力。但是，只有当我们的检测技术能够勾勒 出全基因组的DNA甲基化图谱时，对于DNA甲基化这一表观遗传修饰改变有关的疾病才有 可能得到全面的认识。因此，全基因组DNA甲基化检测技术的发展显得尤其重要，它的发展 是表观遗传学、表观基因组学研宄的重要基础，也将会给当今分子遗传学研宄带来新的变 革。
[0004] 而目前现有的检测技术，还是受到成本、分析周期、基因组覆盖度、分辨率及技术 可操作性等因素的影响：
[0005] (1)对于DNA甲基化检测技术领域，早期的检测手段是KristenH等结合454测序 仪对DNA甲基化进行了检测。该方法首先将基因组DNA进行亚硫酸盐处理，用含有公共接 头的引物对目的片段逐个扩增。利用454新一代测序仪对超过40例样本的25个相关基因 的CpG富集区靶点片段同时测序，共产生了 294631个序列。该方法首次利用高通量测序仪 精确定量检测了目的片段单个CpG位点的甲基化状态，初步展示了新一代测序仪的特点。
[0006] (2)ShawnJ等利用全基因组亚硫酸盐处理测序（WholeGenomeBisulfite Sequencing)，简称WGBS或MethylC-seq，对拟南芥全基因组的DNA甲基化谱进行高通量测 序。他们通过对结果进行比较分析后发现，无论甲基化位点的密度和序列所含碱基如何， 芯片技术发现甲基化位点的能力远不及亚硫酸氢盐测序法，后者甚至还在结构相对简单、 富含转座子的区域发现了甲基化位点，而芯片技术由于交联反应的限制，很难发现这类甲 基化位点。因此，MethylC-seq被视为DNA甲基化组检测的金标准，它能够实现对被测物种 或样本的全基因组DNA甲基化谱的全面、深度、单碱基分辨率的检测。但是，MethylC-seq 需要获得被测物种至少30倍覆盖度的测序数据量。以人类为例，需要至少90Gb的测序数 据，用目前Illumina最新的TrueseqV4的测序试剂价格来衡量，需要约4万人民币；并且 90Gb的亚硫酸盐处理后测序数据需耗费较长的运算时间，目前分析亚硫酸盐处理后测序数 据的常用软件BSMAP，在8核24GB内存的情况下，需要约22天的时间方能得到最终甲基化 图谱。另外，在大部分哺乳动物和植物中，甲基化的胞嘧啶（5mC)主要发生在胞嘧啶一鸟嘌 呤二核苷酸（CpG)上，约只占了全基因组所有碱基数量的1-6%，这使得MethylC-seq所得 到的数据中，仅有20 - 30%数据是有效地提供了DNA甲基化的信息。因此，高额的成本和 费时的运算，大大限制了MethylC-seq的大规模推广应用，尤其是在进行大型基因组的物 种或多样本的DNA甲基化图谱比较研宄中。
[0007] 然而，不同的甲基化DNA富集方法为基于新一代测序技术的DNA甲基化图谱检测 的成本控制奠定了基础。这些富集方法主要包括免疫共沉淀和限制性内切酶等。免疫共沉 淀方法特异性高，而结合限制性酶的新一代测序方法将具有较高的灵敏度。虽然每种方法 都有一定的局限性，但是它们为在基因组范围选择功能区来研宄DNA甲基化图谱提供了更 多的选择。
[0008] (3)基于限制性酶切的测序方法
[0009] 限制性内切酶也多应用于甲基化相关基因的鉴定，是表观遗传学研宄的重要工具 之一。已有3种限制性酶在甲基化研宄中得到应用：第一，甲基化敏感性内切酶，如BstUI、 HpaII、HhaI、SmaI以及Notl。这些酶能够识别CG富集区中的非甲基化位点，而甲基化位点 因为受到甲基的保护而不被识别。第二，甲基化依赖型限制性酶，它能识别并酶切甲基化的 CG位点。McrBC是最有代表性且最常使用的甲基化依赖型限制性酶。具体而言，McrBC可 以识别两个甲基化的半位点（RmC，R=A或G)且这两个位点存在于40~3000bp之内，酶 切发生在两个位点之间。第三，CG甲基化不敏感的同裂酶。例如MspI是Hpall的同裂酶， 他们的识别位点相同但前者不受识别位点甲基化状态的影响。与MspI相似，Xmal是Smal 的同裂酶，它们都识别'CCCGGG'位点，但是Smal酶切后生成钝端而Xmal生成5'粘性末 端。基于此，通过以上几种甲基化敏感或不敏感的限制性内切酶，对被测基因组进行酶切 后，能够特异性的富集高CpG密度的区域，再对其进行二代测序文库的制备以及高通量测 序，从而能够得到单碱基分辨率的DNA甲基化图谱。此类方法中具有代表性的是RRBS和 MRE-seq。针对人类物种而言，RRBS利用MspI对基因组DNA酶切后，然后选取40_220bp的 区域的DNA片段进行亚硫酸盐处理以及高通量测序，仅需大约3Gb的测序数据，便能够覆盖 CpG岛（CGI)的40%左右和启动子区域20%左右的CpG位点，并且由于CGI和启动子区域 在基因表达调控中的重要性，RRBS得到大规模的推广应用。同样，MRE-seq通过甲基化敏感 性内切酶（Bstn、Hpall、Smal)对基因组DNA进行酶切后，进行文库制备和高通量测序，大 约能涵盖人类基因组6%的CpG位点，只需3Gb的总数据量，便能够达到饱和。但目前大量 有关于肿瘤发生过程中DNA甲基化组改变的研宄中发现，差异甲基化区域（Differential methylationregion，DMR)主要发生在CGI的侧翼区域（CGIshore)，并且这些区域对于基 因表达的调控更为明显；同时重复序列元件（Repeatselement，RE) -般呈现高度甲基化， 这些元件上的甲基化状态与基因组的稳定性密切相关。而目前基于限制性内切酶的DNA甲 基化组测序技术，均不能很好的揭示CGIshore和RE区域上的DNA甲基化状态，因此它们得 到的DNA甲基化组不能够代表真正意义的全基因组DNA甲基化图谱。
[0010] (5)基于免疫共沉淀甲基化DNA片段的测序方法
[0011] 哺乳动物MBD家族由五个成员组成，包括MeCP2、MBD1、MBD2、MBD3和MBD4,甲基 化的CpG二核苷酸可被MBD特异性识别并结合。最近，一种基于重组抗体样蛋白MBD的结 合免疫沉淀和高通量测序技术的方法被应用于基因组DNA甲基化图谱的研宄，这种方法被 称为MBD-seq。另外，通过5-甲基胞嘧啶抗体也可用来进行富集甲基化DNA片段，而后结 合高通量测序技术，被称为甲基化DNA免疫共沉淀测序（MeDIP-seq)。基于免疫共沉淀原 理的DNA甲基化组检测技术，通过对甲基化修饰的DNA片段进行特异性富集，撇弃了非甲基 化修饰的DNA片段，对后续的数据产出量提供了良好的成本控制基础。以人类基因组为例， MeDIP-seq只需要大约25million的reads数量，便可涵盖80%的CpGs，大幅度的降低了测 序成本和分析周期。该方法具有较高的特异性，但是也有一定局限性一一它 们没有结合亚 硫酸盐处理或者甲基化敏感的限制性内切酶酶切，无法得到单碱基分辨率的DNA甲基化 图谱，其分辨率大约为l〇〇bp。
[0012] (6)针对免疫共沉淀甲基化DNA测序数据的分析方法
[0013] 由于MethylC-seq的高成本、RRBS和MRE-seq的基因组低覆盖度、MeDIP-seq 和MBD-seq的低分辨率，促使研宄人员利用生物信息算法去更好维持成本、基因组覆盖度 和分辨率之间的平衡，从而使得高通量测序更好服务于DNA甲基化组的研宄。基于免疫 共沉淀甲基化DNA的高通量测序数据，人们开发了一系列生物信息学算法，具体如下。针 对MBD-seq和MeDIP-seq的测序数据，MEDME和BayMeth分别被开发了出来，它们能够将 测序数据中得到的reads数量转换成窗口大小为100bp的区域甲基化水平（Riebleret al. 2014) ;MEDIPS可以将MelDP-seq数据进行计算，从而实现单碱基分辨率的DNA甲基化 图谱，但其得到的并非C的甲基化水平，而是介于1 一 1000的MEDIPS值，导致其得到的结 果无法和其他单碱基分辨率的DNA甲基化检测技术得到的结果相互比较；Batman是另一 款针对MeDIP-seq数据的生物信息学算法，其可实现单个CpG的DNA甲基化水平预测，但 所耗费的计算周期较长（Bock2012)，且算法较为复杂，开发人员也未提供完整的代码安 装文件，使得其他研宄者无法很好的重复该算法，另外batman所得到CpG甲基化水平往往 较真实甲基化水平偏低（Riebleretal. 2014);近来，一种基于条件随机场算法的机器自 动学习工具被开发出来，用于MeDIP-seq和MRE-seq整合后数据的单碱基分辨率DNA甲基 化图谱预测，但其应用于人胚胎干细胞系H1所得到的结果，通过与MethylC-seq测序所得 到的结果进行相关性分析发现，pearson系数仅达到0.77(Stevensetal. 2013)，并且由 于MethylCRF未考虑拷贝数变异对于DNA甲基化水平预测结果的影响，使得该方法在应用 于肿瘤发生时，可能会得到更失真的DNA甲基化组（Laird2010;Robinsonetal.2012; Riebleretal. 2014)〇
[0014] (7)做为本领域所公知的现有技术，MB-seq-甲基化DNA富集结合亚硫酸盐翻转 的甲基化检测技术（MeDIPbisulfitesequencing,MB-seq)拥有许多优点：它是一种高通 量的、单碱基分辨率的、低成本的、可适用于多种已知序列物种的DNA甲基化检测技术，但 是MB-seq存在甲基化水平的偏差，MB-seq甲基化水平被线性放大，所以MB-seq得到的单 个CpG位点的甲基化水平是相对甲基化水平。
【发明内容】
[0015] 针对现有技术中MB-seq-甲基化DNA富集结合亚硫酸盐翻转的甲基化检测技术 存在的甲基化水平的偏差的问题。发明人研发了一种基于岭回归矫正MB-seq甲基化水平 的方法及系统，考虑多种与DNA甲基化水平相关的因素，可以将MB-seq的相对甲基化水平 矫正到全基因组的胞嘧啶位点的绝对甲基化水平。
[0016] 本发明提供的基于岭回归矫正MB-seq甲基化水平的方法，包括以下步骤：
[0017] (1)提取信息
[0018] ⑵建模
[0019] (3)岭回归计算；
[0020] 其中，所述的步骤（1)需要提取的信息有：从参考基因组序列中提取基因组CpG密 度、GC含量和CpG-〇E值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每 个胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个 胞嘧啶的绝对甲基化信息；
[0021] 所述的步骤（2)建模如下：
[0023]其中：
[0024]y:目标函数；为以RRBS高通量测序数据唯一比对结果中提取到的覆盖到的每个 胞嘧啶的绝对甲基化信息；
[0025] x:回归变量矩阵；包括行、列；每行代表每个CpG变量；每列分别为每个变量的 CpG密度、GC含量、CpG-〇E值以及相对甲基化信息；
[0026] 所述的步骤（3)岭回归计算具体是
[0028] 对求导，结果为
[0029] 2XT (Y-XW) -2 入W
[0030] 令其为0,求得的值：
[0032] 输入新的回归变量矩阵X即可获得新Y值，即而获得全基因组的胞嘧啶位点的绝 对甲基化水平。
[0033] 优选的，在所述的步骤（1)中提取信息后，将提取到的信息进行阈值过滤，过滤低 质量碱基和序列，并过滤adapter污染序列。以得到更为精确合理的甲基化水平。
[0034]优选的，在所述的步骤（3)计算之前，采用交叉验证评估模型进行数据训练和测 试：
[0035]a).将预测特征变量和真实的甲基化水平分成训练和测试数据集；随机抽取50% 的CpG位点作为训练数据，剩下的50%作为测试数据；
[0036]b).先使用训练数据训练模型；再计算预测甲基化水平值和RRBS测量的甲基化水 平值之间的相关性系数；这个过程重复N次，N次的平均相关性系数用来表示模型的预测精 度；优选的，N多1000。可以获得更为精确的数据。
[0037] 对于每个基因组元件，单独进行训练和岭回归测试；而对同时位于多个基因组元 件的CpG位点，；取多个预测值的平均值；
[0038] c).甲基化水平的预测是全基因组范围的，并且对于RRBS原本就覆盖的位点，采 取RRBS的观测值作为最终的甲基化水平；所有未被RRBS覆盖的CpG位点，一律认为其未被 甲基化，并且不用于岭回归，甲基化水平预测值小于0或者大于分别基于岭回归的原则规 整到0和1。
[0039] 优选的，在在模型数据训练时：
[0040]a).当变量间存在共线性的时候，通过引入lambda表达式以解决最小二乘回归得 到的系数不稳定，方差很大的问题；
[0041]b).当模型包含常数项时，岭回归函数对y进行中心化，以y的均值作为因子；对 X进行中心化和归一化，以x中各个变量的均值和标准差作为因子；这样对x和y处理后，x 和y的均值为〇,这使得回归平面经过原点，即常数项为〇 ;
[0042]c).当模型不包含常数项时，因为要强制通过原点，该模型假设各个变量的均值为 0,因此不对X和y进行中心化，但是对x进行归一化，而且归一化因子也是假设变量均值为 〇计算出来的该变量的标准差。
[0043] 优选的，在使用该模型进行测试的时候，需要首先对x和y进行中心化和归一化， 此时因子是使用训练模型时候进行中心化和归一化的因子，然后再与系数相乘得到预测结 果。
[0044] 优选的，在步骤（3)岭回归计算之后，进行如下对异常点处理：
[0045] 1)将MB-seq检测深度为0的位点定义为甲基化水平为0;
[0046]2)结合MB-seq甲基化水平的观测值（MB level)，甲基化CpG个数（MB mCG)，MB-seq测序深度（MB depth)，当前CpG侧翼+/_100bp区域的MB-seq检测到的平均 甲基化水平（MB back),这200bp范围的甲基化CpG位点总数（MB mCG)，以及每一个CpG位 点上下游l〇〇bp的基因组CpG密度、GC含量，CpG-〇E值等对甲基化水平检测的影响，利用 岭回归导入到模型中，并且机器学习得到某一胞嘧啶位点甲基化水平；
[0047] 3)将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回归的甲 基化水平值小于〇的位点自动归为甲基化水平为0。
[0048] 优选的，所述的相对甲基化信息包括：MB-seq甲基化水平的观测值MBlevel，甲 基化CpG个数MBmCG,MB-seq测序深度MBdepth,当前CpG侧翼+/_100bp区域的MB-seq 检测到的平均甲基化水平MBback,这200bp范围的甲基化CpG位点总数MBmCG。
[0049] 岭回归矫正MB-seq甲基化水平的系统，包括以下模块：
[0050] 提取模块：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-〇E值；从 MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信 息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信 息；所述的相对甲基化信息包括 ：MB-seq甲基化水平的观测值MBlevel，甲基化CpG个数 MBmCG,MB-seq测序深度MBdepth,当前CpG侧翼+/_100bp区域的MB-seq检测到的平均 甲基化水平MBback,这200bp范围的甲基化CpG位点总数MBmCG ;
[0051] 建模模块：根据基因组信息和甲基化信 息，采用岭回归模型对真实甲基化水平 RRBSlevel和回归参数建立回归模型；
[0052] 回归模块：利用岭回归理论，并依据提取出来的基因组信息和甲基化信息，对基因 组上的胞嘧啶位点进行回归以得到甲基化水平的模块。
[0053] 与现有技术相比较，本发明具有如下有益效果：
[0054] 利用数学模型，对基因组上的RRBS未覆盖的胞嘧啶位点进行岭回归预测，使得甲 基化水平检测的准确度大于95%，从而消除MB-seq的偏差并得到全基因组甲基化图谱。利 用本发明可以从高通量测序MB-seq数据中，精确计算全基因组每一个CpG的甲基化水平。
[0055]
[0056]
【具体实施方式】
[0057] 基于岭回归矫正MB-seq甲基化水平的方法，包括以下步骤：
[0058] (1)提取信息⑵建模⑶岭回归计算；
[0059] 其中，所述的步骤（1)需要提取的信息有：从参考基因组序列中提取基因组CpG密 度、GC含量和CpG-〇E值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每 个胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个 胞嘧啶的绝对甲基化信息；
[0060] 所述的步骤（2)建模如下：
[0062]其中：
[0063]y:目标函数；为以RRBS高通量测序数据唯一比对结果中提取到的覆盖到的每个 胞嘧啶的绝对甲基化信息；
[0064] X:回归变量矩阵；包括行、列；每行代表每个CpG变量；每列分别为每个变量的 CpG密度、GC含量、CpG-〇E值以及相对甲基化信息；
[0065] 所述的步骤（3)岭回归计算具体是
[0067] 对求导，结果为
[0068] 2XT(Y-XW) -2 入W
[0069] 令其为0,求得的值：
[0071] 输入新的回归变量矩阵X即可获得新Y值，即而获得全基因组的胞嘧啶位点的绝 对甲基化水平。在所述的步骤（1)中提取信息后，将提取到的信息进行阈值过滤，过滤低 质量碱基和序列，并过滤adapter污染序列。以得到更为精确合理的甲基化水平。在所述 的步骤（3)计算之前，采用交叉验证评估模型进行数据训练和测试：a).将预测特征变量和 真实的甲基化水平分成训练和测试数据集；随机抽取50%的CpG位点作为训练数据，剩下 的50%作为测试数据；b).先使用训练数据训练模型；再计算预测甲基化水平值和RRBS测 量的甲基化水平值之间的相关性系数；这个过程重复N次，N次的平均相关性系数用来表 示模型的预测精度；优选的，N多1000。可以获得更为精确的数据。对于每个基因组元件， 单独进行训练和岭回归测试；而对同时位于多个基因组元件的CpG位点，；取多个预测值的 平均值；c).甲基化水平的预测是全基因组范围的，并且对于RRBS原本就覆盖的位点，采取 RRBS的观测值作为最终的甲基化水平；所有未被RRBS覆盖的CpG位点，一律认为其未被甲 基化，并且不用于岭回归，甲基化水平预测值小于0或者大于分别基于岭回归的原则规整 到0和1。在在模型数据训练时：a).当变量间存在共线性的时候，通过引入lambda表达 式以解决最小二乘回归得到的系数不稳定，方差很大的问题；b).当模型包含常数项时，岭 回归函数对y进行中心化，以y的均值作为因子；对x进行中心化和归一化，以x中各个变 量的均值和标准差作为因子；这样对x和y处理后，x和y的均值为0,这使得回归平面经 过原点，即常数项为0 ;c).当模型不包含常数项时，因为要强制通过原点，该模型假设各个 变量的均值为〇,因此不对x和y进行中心化，但是对x进行归一化，而且归一化因子也是 假设变量均值为〇计算出来的该变量的标准差。在使用该模型进行测试的时候，需要首先 对x和y进行中心化和归一化，此时因子是使用训练模型时候进行中心化和归一化的因子， 然后再与系数相乘得到预测结果。在步骤（3)岭回归计算之后，进行如下对异常点处理： 1)将MB-seq检测深度为0的位点定义为甲基化水平为0 ;2)结合MB-seq甲基化水平的观 测值（MBlevel),甲基化CpG个数（MBmCG)，MB-seq测序深度（MBd印th)，当前CpG侧翼 +/-100bp区域的MB-seq检测到的平均甲基化水平（MBback),这200bp范围的甲基化CpG 位点总数（MBmCG)，以及每一个CpG位点上下游100bp的基因组CpG密度、GC含量，CpG-〇E 值等对甲基化水平检测的影响，利用岭回归导入到模型中，并且机器学习得到某一胞嘧啶 位点甲基化水平；3)将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回 归的甲基化水平值小于0的位点自动归为甲基化水平为0。所述的相对甲基化信息包括： MB-seq甲基化水平的观测值MBlevel,甲基化CpG个数MBmCG,MB-seq测序深度MBdepth, 当前CpG侧翼+/_100bp区域的MB-seq检测到的平均甲基化水平MBback,这200bp范围 的甲基化CpG位点总数MBmCG。
[0072] 岭回归模型是采用岭回归实现的一种正则化线性回归。当多个预测因子含有非0 系数并且呈现正态分布的时候，岭回归是理想的回归方法，岭回归对每一个预测因子影响 小的模型尤其适用，并且它防止线性回归模型系数由于共线性而导致无法模拟和高变异 性。岭回归对共线预测因子的系数收缩并趋于零，例如，给出k相同的预测因子，都将获 得相同的系数等于1/k的单个因子的回归系数值。因此，岭回归不会使某个因子消失，不 能将某些因子摒弃来获得最优预测数据集。岭回归（2)估计解决回归问题（1)使用￡2惩罚 最小二乘法：
[0073]y=yln+X|3 +e: (1)
[0074]y=(yi，…，yn)T其中是观察表型的向量，ln是一个n维列向量的，y是一种常 见的截距，是nXp矩阵的表示，|3表示回归系数的向量，61是残差的向量和
是残差的误差。
[0076]其中
[0078] 是U-normi(二次方程式）损失函数（即残差平方和），X,7是X向量的第i个 行。
[0080] AH:2-norm坫于|3罚分，入彡0是调优（罚分，正规化，或复杂化）参数，这 些参数通过相对重要性决定经验误差和惩罚调节罚分的强度（即线性收缩）。A值越大， 收缩量越大。A的值依赖于数据，通过数据驱动的方法（交叉验证）进行确定使用。
[0081] CpG密度、GC含量以及CpG-〇E值三者计算方法分别为：
[0082] CpG密度：某一个CpG上下游各100bp范围内CpG个数除于201bp长度得到此CpG 位点的CpG密度；
[0083] GC含量：某一个CpG上下游各100bp范围内C和G总数除于201bp长度得到此CpG 位点的GC含量；
[0084] CpG-〇E值：CpG上下游各100bp范围内CpG个数乘于210bp，然后除于C和G个数 的乘积。
[0085] 岭回归矫正MB-seq甲基化水平的系统，包括以下模块：
[0086] 提取模块：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-〇E值；从 MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信 息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信 息；所述的相对甲基化信息包括 ：MB-seq甲基化水平的观测值MBlevel，甲基化CpG个数 MBmCG,MB-seq测序深度MBdepth,当前CpG侧翼+/_100bp区域的MB-seq检测到的平均 甲基化水平MBback,这200bp范围的甲基化CpG位点总数MBmCG ;
[0087] 建模模块：根据基因组信息和甲基化信息，采用岭回归模型对真实甲基化水平 RRBSlevel和回归参数建立回归模型；
[0088] 回归模块：利用岭回归理论，并依据提取出来的基因组信息和甲基化信息，对基因 组上的胞嘧啶位点进行回归以得到甲基化水平的模块。
[0089] 优选的，提取模块提取信息后，设置阈值过滤模块，还将提取到的信息进行阈值过 滤，过滤低质量碱基和序列，并过滤adapter污染序列。
[0090] 优选的，在回归模块计算之前，还设置数据训练和测试模块，采用交叉验证评估模 型进行数据训练和测试：
[0091] a).将预测特征变量和真实的甲基化水平分成训练和测试数据集；随机抽取50% 的CpG位点作为训练数据，剩下的50%作为测试数据；
[0092] b).先使用训练数据训练模型；再计算预测甲基化水平值和RRBS测量的甲基化水 平值之间的相关性系数；这个过程重复N次，N次的平均相关性系数用来表示模型的预测精 度；
[0093] 对于每个基因组元件，单独进行训练和岭回归测试；而对同时位于多个基因组元 件的CpG位点 ，；取多个预测值的平均值；
[0094] c).甲基化水平的预测是全基因组范围的，并且对于RRBS原本就覆盖的位点，采 取RRBS的观测值作为最终的甲基化水平；所有未被RRBS覆盖的CpG位点，一律认为其未被 甲基化，并且不用于岭回归，甲基化水平预测值小于0或者大于分别基于岭回归的原则规 整到0和1。
[0095] 优选的，N彡 1000。
[0096] 数据训练和测试模块对模型数据训练时：
[0097] a).当变量间存在共线性的时候，通过引入lambda表达式以解决最小二乘回归得 到的系数不稳定，方差很大的问题；
[0098] b).当模型包含常数项时，岭回归函数对y进行中心化，以y的均值作为因子；对 X进行中心化和归一化，以x中各个变量的均值和标准差作为因子；这样对x和y处理后，x 和y的均值为〇,这使得回归平面经过原点，即常数项为〇 ;
[0099]c).当模型不包含常数项时，因为要强制通过原点，该模型假设各个变量的均值为 0,因此不对X和y进行中心化，但是对x进行归一化，而且归一化因子也是假设变量均值为 〇计算出来的该变量的标准差。
[0100] 数据训练和测试模块对模型进行测试的时候，需要首先对x和y进行中心化和归 一化，此时因子是使用训练模型时候进行中心化和归一化的因子，然后再与系数相乘得到 预测结果。
[0101] 在回归模块计算之后，设置异常点处理模块，用于对异常点处理：
[0102] 1)将MB-seq检测深度为0的位点定义为甲基化水平为0 ;
[0103] 2)结合MB-seq甲基化水平的观测值（MBlevel)，甲基化CpG个数（MB mCG)，MB-seq测序深度（MBdepth)，当前CpG侧翼+/_100bp区域的MB-seq检测到的平均 甲基化水平（MBback),这200bp范围的甲基化CpG位点总数（MBmCG)，以及每一个CpG位 点上下游l〇〇bp的基因组CpG密度、GC含量，CpG-〇E值等对甲基化水平检测的影响，利用 岭回归导入到模型中，并且机器学习得到某一胞嘧啶位点甲基化水平；
[0104] 3)将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回归的甲 基化水平值小于〇的位点自动归为甲基化水平为0。
【主权项】
1. 基于岭回归矫正MB-Seq甲基化水平的方法，其特征在于，包括以下步骤： (1) 提取信息 (2) 建模 (3) 岭回归计算； 其中，所述的步骤（1)需要提取的信息有：从参考基因组序列中提取基因组CpG密度、 GC含量和CpG-OE值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个 胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞 嘧啶的绝对甲基化信息； 所述的步骤（2)建模如下：其中： y:目标函数；为以RRBS高通量测序数据唯一比对结果中提取到的覆盖到的每个胞嘧 啶的绝对甲基化信息； X:回归变量矩阵；包括行、列；每行代表每个CpG变量；每列分别为每个变量的CpG密 度、GC含量、CpG-OE值以及相对甲基化信息； 所述的相对甲基化信息包括：MB-seq甲基化水平的观测值MB level，甲基化CpG个数 MB mCG, MB-seq测序深度MB depth,当前CpG侧翼+/-IOObp区域的MB-seq检测到的平均 甲基化水平MB back,这200bp范围的甲基化CpG位点总数MB mCG ; 所述的步骤（3)岭回归计算具体是对求导，结果为 2XT (Y-XW) -2 λ W 令其为〇,求得的值：输入新的回归变量矩阵X即可获得新Y值，即而获得全基因组的胞嘧啶位点的绝对甲 基化水平。2. 如权利要求1所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：在所 述的步骤（1)中提取信息后，还将提取到的信息进行阈值过滤，过滤低质量碱基和序列，并 过滤adapter污染序列。3. 如权利要求1所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：所述 的步骤（3)计算之前，采用交叉验证评估模型进行数据训练和测试： a) .将预测特征变量和真实的甲基化水平分成训练和测试数据集；随机抽取50%的 CpG位点作为训练数据，剩下的50%作为测试数据； b) .先使用训练数据训练模型；再计算预测甲基化水平值和RRBS测量的甲基化水平值 之间的相关性系数；这个过程重复N次，N次的平均相关性系数用来表示模型的预测精度； 对于每个基因组元件，单独进行训练和岭回归测试；而对同时位于多个基因组元件的 CpG位点，；取多个预测值的平均值； c).甲基化水平的预测是全基因组范围的，并且对于RRBS原本就覆盖的位点，采取 RRBS的观测值作为最终的甲基化水平；所有未被RRBS覆盖的CpG位点，一律认为其未被甲 基化，并且不用于岭回归，甲基化水平预测值小于O或者大于分别基于岭回归的原则规整 到O和1。4. 如权利要求3所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于： N 彡 1000。5. 如权利要求3所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：在模 型数据训练时： a) .当变量间存在共线性的时候，通过引入lambda表达式以解决最小二乘回归得到的 系数不稳定，方差很大的问题； b) .当模型包含常数项时，岭回归函数对y进行中心化，以y的均值作为因子；对X进 行中心化和归一化，以X中各个变量的均值和标准差作为因子；这样对X和y处理后，X和 y的均值为〇,这使得回归平面经过原点，即常数项为〇 ; c) .当模型不包含常数项时，因为要强制通过原点，该模型假设各个变量的均值为0， 因此不对X和y进行中心化，但是对X进行归一化，而且归一化因子也是假设变量均值为0 计算出来的该变量的标准差。6. 如权利要求3所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：在使 用该模型进行测试的时候，需要首先对X和y进行中心化和归一化，此时因子是使用训练模 型时候进行中心化和归一化的因子，然后再与系数相乘得到预测结果。7. 如权利要求1所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：在步 骤（3)岭回归计算之后，进行如下对异常点处理： 1) 将MB-seq检测深度为0的位点定义为甲基化水平为0 ; 2) 结合MB-seq甲基化水平的观测值（MB level),甲基化CpG个数（MB mCG)，MB-seq测 序深度（MB depth),当前CpG侧翼+/-IOObp区域的MB-seq检测到的平均甲基化水平（MB back)，这200bp范围的甲基化CpG位点总数（MB mCG)，以及每一个CpG位点上下游IOObp 的基因组CpG密度、GC含量，CpG-OE值等对甲基化水平检测的影响，利用岭回归导入到模 型中，并且机器学习得到某一胞嘧啶位点甲基化水平； 3) 将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回归的甲基化 水平值小于〇的位点自动归为甲基化水平为0。8. 岭回归矫正MB-seq甲基化水平的系统，其特征在于：包括以下模块： 提取模块：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-OE值；从MB-seq 高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信息；从 RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信息；所 述的相对甲基化信息包括：MB-seq甲基化水平的观测值MB level,甲基化CpG个数MB mCG, MB-seq测序深度MB depth,当前CpG侧翼+/-IOObp区域的MB-seq检测到的平均甲基 化水平MB back,这200bp范围的甲基化CpG位点总数MB mCG ; 建模模块：根据基因组信息和甲基化信息，采用岭回归模型对真实甲基化水平RRBS level和回归参数建立回归模型； 回归模块：利用岭回归理论，并依据提取出来的基因组信息和甲基化信息，对基因组上 的胞嘧啶位点进行回归以得到甲基化水平的模块。
【专利摘要】基于岭回归矫正MB-seq甲基化水平的方法，属于基因工程技术领域，利用机器学习岭回归理论，并依据RRBS检测出来的绝对MB-seq甲基化水平进行数据训练并建立预测模型，对基因组上的RRBS未覆盖的胞嘧啶位点进行岭回归预测，使得甲基化水平检测的准确度大于95％，从而消除MB-seq的偏差并得到全基因组甲基化图谱。本发明还公开了一种基于岭回归的甲基化水平计算系统。利用本发明可以从高通量测序MB-seq数据中，精确计算全基因组每一个CpG的甲基化水平。
【IPC分类】G06F19/12
【公开号】CN104899474
【申请号】CN201510313520
【发明人】张保荣, 王晓东, 张久文
【申请人】大连三生科技发展有限公司
【公开日】2015年9月9日
【申请日】2015年6月9日
最新回复(0)