CN104899474A - 基于岭回归矫正MB-seq甲基化水平的方法及系统 - Google Patents

Abstract

基于岭回归矫正MB-seq甲基化水平的方法，属于基因工程技术领域，利用机器学习岭回归理论，并依据RRBS检测出来的绝对MB-seq甲基化水平进行数据训练并建立预测模型，对基因组上的RRBS未覆盖的胞嘧啶位点进行岭回归预测，使得甲基化水平检测的准确度大于95％，从而消除MB-seq的偏差并得到全基因组甲基化图谱。本发明还公开了一种基于岭回归的甲基化水平计算系统。利用本发明可以从高通量测序MB-seq数据中，精确计算全基因组每一个CpG的甲基化水平。

Description

基于岭回归矫正MB-seq甲基化水平的方法及系统

技术领域

本发明属于基因工程技术领域，尤其涉及一种基于数学模型-岭回归矫正MB-seq甲基化水平的方法及系统。

背景技术

DNA甲基化(DNA methylation)是最早发现的修饰途径之一，大量研究表明，DNA甲基化能引起染色质结构、DNA构象、DNA稳定性及DNA与蛋白质相互作用方式的改变，从而控制基因表达。早在1942年，C.H.Waddinton就提出了表观遗传学的概念，他指出，表观遗传与遗传相对，主要研究基因型和表型的关系。而现在，对于表观遗传学，比较统一的认识是，其研究在没有细胞核DNA序列改变的情况时，基因功能的可逆的可遗传的改变。也就是说，在不改变基因组序列的前提下，通过DNA和组蛋白的修饰等来调控基因表达，其中又以DNA甲基化(DNA methylation)最为常见，DNA甲基化成为表观遗传学的重要组成部分。随着人类基因组计划的开展，科学家们开始在基因组水平来研究表观遗传学，逐步形成表观基因组学(epigenomics)。表观基因组学就是要在整个基因组水平来研究表观遗传过程以及与这些过程密切相关的特定基因组区域的识别与鉴定。2000年10月，人类表观基因组协会(Human Epigenome Consortium)由欧盟赞助，启动了旨在于人类6号染色体MHC区域首先做出DNA的甲基化图谱的先导计划(Pilot Project)。该计划顺利完成，引导启动了2003年的人类表观基因组计划(Human Epigenome Project，HEP)。2005年，美国国家卫生院(NIH)下属的国立癌症研究所启动了癌症基因组先导计划。2006年，该所与国立人类基因组研究所一起共同启动癌症基因组计划(Cancer Genome Project)。表观基因组学和DNA甲基化与癌症的研究成为新的热点。

目前，人们认识到DNA甲基化对基因组正常功能维持是必要的，表观遗传水平的改变可看作是复杂疾病(如癌症，精神疾病)发病机理的第一步(病因学上)，以表观遗传学为基础的药物对治疗复杂疾病将有巨大的潜力。但是，只有当我们的检测技术能够勾勒出全基因组的DNA甲基化图谱时，对于DNA甲基化这一表观遗传修饰改变有关的疾病才有可能得到全面的认识。因此，全基因组DNA甲基化检测技术的发展显得尤其重要，它的发展是表观遗传学、表观基因组学研究的重要基础，也将会给当今分子遗传学研究带来新的变革。

而目前现有的检测技术，还是受到成本、分析周期、基因组覆盖度、分辨率及技术可操作性等因素的影响：

(1)对于DNA甲基化检测技术领域，早期的检测手段是Kristen H等结合454测序仪对DNA甲基化进行了检测。该方法首先将基因组DNA进行亚硫酸盐处理，用含有公共接头的引物对目的片段逐个扩增。利用454新一代测序仪对超过40例样本的25个相关基因的CpG富集区靶点片段同时测序，共产生了294631个序列。该方法首次利用高通量测序仪精确定量检测了目的片段单个CpG位点的甲基化状态，初步展示了新一代测序仪的特点。

(2)Shawn J等利用全基因组亚硫酸盐处理测序(Whole Genome Bisulfite Sequencing)，简称WGBS或MethylC-seq，对拟南芥全基因组的DNA甲基化谱进行高通量测序。他们通过对结果进行比较分析后发现，无论甲基化位点的密度和序列所含碱基如何，芯片技术发现甲基化位点的能力远不及亚硫酸氢盐测序法，后者甚至还在结构相对简单、富含转座子的区域发现了甲基化位点，而芯片技术由于交联反应的限制，很难发现这类甲基化位点。因此，MethylC-seq被视为DNA甲基化组检测的金标准，它能够实现对被测物种或样本的全基因组DNA甲基化谱的全面、深度、单碱基分辨率的检测。但是，MethylC-seq需要获得被测物种至少30倍覆盖度的测序数据量。以人类为例，需要至少90Gb的测序数据，用目前Illumina最新的Trueseq V4的测序试剂价格来衡量，需要约4万人民币；并且90Gb的亚硫酸盐处理后测序数据需耗费较长的运算时间，目前分析亚硫酸盐处理后测序数 据的常用软件BSMAP，在8核24GB内存的情况下，需要约22天的时间方能得到最终甲基化图谱。另外，在大部分哺乳动物和植物中，甲基化的胞嘧啶(5mC)主要发生在胞嘧啶－鸟嘌呤二核苷酸(CpG)上，约只占了全基因组所有碱基数量的1-6％，这使得MethylC-seq所得到的数据中，仅有20－30％数据是有效地提供了DNA甲基化的信息。因此，高额的成本和费时的运算，大大限制了MethylC-seq的大规模推广应用，尤其是在进行大型基因组的物种或多样本的DNA甲基化图谱比较研究中。

然而，不同的甲基化DNA富集方法为基于新一代测序技术的DNA甲基化图谱检测的成本控制奠定了基础。这些富集方法主要包括免疫共沉淀和限制性内切酶等。免疫共沉淀方法特异性高，而结合限制性酶的新一代测序方法将具有较高的灵敏度。虽然每种方法都有一定的局限性，但是它们为在基因组范围选择功能区来研究DNA甲基化图谱提供了更多的选择。

(3)基于限制性酶切的测序方法

限制性内切酶也多应用于甲基化相关基因的鉴定，是表观遗传学研究的重要工具之一。已有3种限制性酶在甲基化研究中得到应用：第一，甲基化敏感性内切酶，如BstUI、HpaII、HhaI、SmaI以及NotI。这些酶能够识别CG富集区中的非甲基化位点，而甲基化位点因为受到甲基的保护而不被识别。第二，甲基化依赖型限制性酶，它能识别并酶切甲基化的CG位点。McrBC是最有代表性且最常使用的甲基化依赖型限制性酶。具体而言，McrBC可以识别两个甲基化的半位点(RmC，R＝A或G)且这两个位点存在于40～3000bp之内，酶切发生在两个位点之间。第三，CG甲基化不敏感的同裂酶。例如MspI是HpaII的同裂酶，他们的识别位点相同但前者不受识别位点甲基化状态的影响。与MspI相似，XmaI是SmaI的同裂酶，它们都识别‘CCCGGG’位点，但是SmaI酶切后生成钝端而XmaI生成5′粘性末端。基于此，通过以上几种甲基化敏感或不敏感的限制性内切酶，对被测基因组进行酶切后，能够特异性的富集高CpG密度的区域，再对其进行二代测序文库的制备以及高通量测序，从而能够得到单碱基分辨率的DNA甲 基化图谱。此类方法中具有代表性的是RRBS和MRE-seq。针对人类物种而言，RRBS利用MspI对基因组DNA酶切后，然后选取40-220bp的区域的DNA片段进行亚硫酸盐处理以及高通量测序，仅需大约3Gb的测序数据，便能够覆盖CpG岛(CGI)的40％左右和启动子区域20％左右的CpG位点，并且由于CGI和启动子区域在基因表达调控中的重要性，RRBS得到大规模的推广应用。同样，MRE-seq通过甲基化敏感性内切酶(BstUI、HpaII、SmaI)对基因组DNA进行酶切后，进行文库制备和高通量测序，大约能涵盖人类基因组6％的CpG位点，只需3Gb的总数据量，便能够达到饱和。但目前大量有关于肿瘤发生过程中DNA甲基化组改变的研究中发现，差异甲基化区域(Differential methylation region，DMR)主要发生在CGI的侧翼区域(CGIshore)，并且这些区域对于基因表达的调控更为明显；同时重复序列元件(Repeats element，RE)一般呈现高度甲基化，这些元件上的甲基化状态与基因组的稳定性密切相关。而目前基于限制性内切酶的DNA甲基化组测序技术，均不能很好的揭示CGIshore和RE区域上的DNA甲基化状态，因此它们得到的DNA甲基化组不能够代表真正意义的全基因组DNA甲基化图谱。

(5)基于免疫共沉淀甲基化DNA片段的测序方法

哺乳动物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甲基化图谱，其分辨率大约为100bp。

(6)针对免疫共沉淀甲基化DNA测序数据的分析方法

由于MethylC-seq的高成本、RRBS和MRE-seq的基因组低覆盖度、MeDIP-seq和MBD-seq的低分辨率，促使研究人员利用生物信息算法去更好维持成本、基因组覆盖度和分辨率之间的平衡，从而使得高通量测序更好服务于DNA甲基化组的研究。基于免疫共沉淀甲基化DNA的高通量测序数据，人们开发了一系列生物信息学算法，具体如下。针对MBD-seq和MeDIP-seq的测序数据，MEDME和BayMeth分别被开发了出来，它们能够将测序数据中得到的reads数量转换成窗口大小为100bp的区域甲基化水平(Riebler et al.2014)；MEDIPS可以将MeIDP-seq数据进行计算，从而实现单碱基分辨率的DNA甲基化图谱，但其得到的并非C的甲基化水平，而是介于1－1000的MEDIPS值，导致其得到的结果无法和其他单碱基分辨率的DNA甲基化检测技术得到的结果相互比较；Batman是另一款针对MeDIP-seq数据的生物信息学算法，其可实现单个CpG的DNA甲基化水平预测，但所耗费的计算周期较长(Bock 2012)，且算法较为复杂，开发人员也未提供完整的代码安装文件，使得其他研究者无法很好的重复该算法，另外batman所得到CpG甲基化水平往往较真实甲基化水平偏低(Riebler et al.2014)；近来，一种基于条件随机场算法的机器自动学习工具被开发出来，用于MeDIP-seq和MRE-seq整合后数据的单碱基分辨率DNA甲基化图谱预测，但其应用于人胚胎干细胞系H1所得到的结果，通过与MethylC-seq测序所得到的结果进行相关性分析发现，pearson系数仅达到0.77(Stevens et al.2013)，并且由于MethylCRF未考虑拷贝数变异对于DNA甲基化水平预测结果的影响，使得该方法在应用于肿瘤发生时，可能会得到更失真的DNA甲基化组(Laird 2010；Robinson et al.2012；Riebler et al.2014)。

(7)做为本领域所公知的现有技术，MB-seq—甲基化DNA富集结合亚硫酸盐翻转的甲基化检测技术(MeDIP bisulfite sequencing,MB-seq)拥有许多优点：它是一种高通量的、单碱基分辨率的、低成 本的、可适用于多种已知序列物种的DNA甲基化检测技术，但是MB-seq存在甲基化水平的偏差，MB-seq甲基化水平被线性放大，所以MB-seq得到的单个CpG位点的甲基化水平是相对甲基化水平。

发明内容

针对现有技术中MB-seq—甲基化DNA富集结合亚硫酸盐翻转的甲基化检测技术存在的甲基化水平的偏差的问题。发明人研发了一种基于岭回归矫正MB-seq甲基化水平的方法及系统，考虑多种与DNA甲基化水平相关的因素，可以将MB-seq的相对甲基化水平矫正到全基因组的胞嘧啶位点的绝对甲基化水平。

本发明提供的基于岭回归矫正MB-seq甲基化水平的方法，包括以下步骤：

(1)提取信息

(2)建模

(3)岭回归计算； 

其中，所述的步骤(1)需要提取的信息有：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-OE值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信息；

所述的步骤(2)建模如下： 

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\underset{j=0}{\overset{P}{\mathrm{\Σ}}}{w}_j{x}_{\mathrm{ij}})}^2+\mathrm{\λ}\underset{j=0}{\overset{P}{\mathrm{\Σ}}}{\mathrm{\ω}}_j^2$

其中：

y:目标函数；为以RRBS高通量测序数据唯一比对结果中提取到的覆盖到的每个胞嘧啶的绝对甲基化信息；

x:回归变量矩阵；包括行、列；每行代表每个CpG变量；每列分别为每个变量的CpG密度、GC含量、CpG-OE值以及相对甲基化信息；

所述的步骤(3)岭回归计算具体是

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\underset{j=0}{\overset{P}{\mathrm{\Σ}}}{w}_j{x}_{\mathrm{ij}})}^2+\mathrm{\λ}\underset{j=0}{\overset{P}{\mathrm{\Σ}}}{\mathrm{\ω}}_j^2$

对求导，结果为

2X^T(Y-XW)-2λW

令其为0，求得的值： 

$\hat{w}={(X^{T}X+\mathrm{\λI})}^{-1}{X}^{T}Y$

输入新的回归变量矩阵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).当模型不包含常数项时，因为要强制通过原点，该模型假设各个变量的均值为0，因此不对x和y进行中心化，但是对x进行归一化，而且归一化因子也是假设变量均值为0计算出来的该变量的标准差。

优选的，在使用该模型进行测试的时候，需要首先对x和y进行中心化和归一化，此时因子是使用训练模型时候进行中心化和归一化的因子，然后再与系数相乘得到预测结果。

优选的，在步骤(3)岭回归计算之后，进行如下对异常点处理：

1)将MB-seq检测深度为0的位点定义为甲基化水平为0；

2)结合MB-seq甲基化水平的观测值(MB level),甲基化CpG个数(MB mCG),MB-seq测序深度(MB depth),当前CpG侧翼+/-100bp区域的MB-seq检测到的平均甲基化水平(MB back),这200bp范围的甲基化CpG位点总数(MB mCG)，以及每一个CpG位点上下游100bp的基因组CpG密度、GC含量,CpG-OE值等对甲基化水平检测的影响，利用岭回归导入到模型中，并且机器学习得到某一胞嘧啶位点甲基化水平；

3)将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回归的甲基化水平值小于0的位点自动归为甲基化水平为0。

优选的，所述的相对甲基化信息包括：MB-seq甲基化水平的观测值MB level,甲基化CpG个数MB mCG,MB-seq测序深度MB depth,当前CpG侧翼+/-100bp区域的MB-seq检测到的平均甲基化水平MB back,这200bp范围的甲基化CpG位点总数MB mCG。

岭回归矫正MB-seq甲基化水平的系统，包括以下模块：

提取模块：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-OE值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信息；所述的相 对甲基化信息包括：MB-seq甲基化水平的观测值MB level,甲基化CpG个数MB mCG,MB-seq测序深度MB depth,当前CpG侧翼+/-100bp区域的MB-seq检测到的平均甲基化水平MB back,这200bp范围的甲基化CpG位点总数MB mCG；

建模模块：根据基因组信息和甲基化信息，采用岭回归模型对真实甲基化水平RRBS level和回归参数建立回归模型；

回归模块：利用岭回归理论，并依据提取出来的基因组信息和甲基化信息，对基因组上的胞嘧啶位点进行回归以得到甲基化水平的模块。

与现有技术相比较，本发明具有如下有益效果：

利用数学模型，对基因组上的RRBS未覆盖的胞嘧啶位点进行岭回归预测，使得甲基化水平检测的准确度大于95％，从而消除MB-seq的偏差并得到全基因组甲基化图谱。利用本发明可以从高通量测序MB-seq数据中，精确计算全基因组每一个CpG的甲基化水平。

具体实施方式

基于岭回归矫正MB-seq甲基化水平的方法，包括以下步骤：

(1)提取信息(2)建模(3)岭回归计算； 

其中，所述的步骤(1)需要提取的信息有：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-OE值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信息；

所述的步骤(2)建模如下： 

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\underset{j=0}{\overset{P}{\mathrm{\Σ}}}{w}_j{x}_{\mathrm{ij}})}^2+\mathrm{\λ}\underset{j=0}{\overset{P}{\mathrm{\Σ}}}{\mathrm{\ω}}_j^2$

其中：

y:目标函数；为以RRBS高通量测序数据唯一比对结果中提取到的覆盖到的每个胞嘧啶的绝对甲基化信息；

x:回归变量矩阵；包括行、列；每行代表每个CpG变量；每列分别为每个变量的CpG密度、GC含量、CpG-OE值以及相对甲基化信息；

所述的步骤(3)岭回归计算具体是

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\underset{j=0}{\overset{P}{\mathrm{\Σ}}}{w}_j{x}_{\mathrm{ij}})}^2+\mathrm{\λ}\underset{j=0}{\overset{P}{\mathrm{\Σ}}}{\mathrm{\ω}}_j^2$

对求导，结果为

2XT(Y-XW)-2λW

令其为0，求得的值： 

$\hat{w}={(X^{T}X+\mathrm{\λI})}^{-1}{X}^{T}Y$

输入新的回归变量矩阵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).当模型不包含常数项时，因为要强制通过原点，该模型假设各个变量的均值为0，因此不对x和y进行中心化，但是对x进行归一化，而且归一化因子也是假设变量均值为0计算出来的该变量的标准差。在使用该模型进行测试的时候，需要首先对x和y进行中心化和归一化，此时因子是使用训练模型时候进行中心化和归一化的因子，然后再与系数相乘得到预测结果。在步骤(3)岭回归计算之后，进行如下对异常点处理：1)将MB-seq检测深度为0的位点定义为甲基化水平为0；2)结合MB-seq甲基化水平的观测值(MB level),甲基化CpG个数(MB mCG),MB-seq测序深度(MB depth),当前CpG侧翼+/-100bp区域的MB-seq检测到的平均甲基化水平(MB back),这200bp范围的甲基化CpG位点总数(MB mCG)，以及每一个CpG位点上下游100bp的基因组CpG密度、GC含量,CpG-OE值等对甲基化水平检测的影响，利用岭回归导入到模型中，并且机器学习得到某一胞嘧啶位点甲基化水平；3)将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回归的甲基化水平值小于0的位点自动归为甲基化水平为0。所述的相对甲基化信息包括：MB-seq甲基化水平的观测值MB level,甲基化CpG个数MB mCG,MB-seq测序深度MB depth,当前CpG侧翼+/-100bp区域的MB-seq检测到的平均甲基化水平MB back,这200bp范围的甲基化CpG位点总数MB mCG。

岭回归模型是采用岭回归实现的一种正则化线性回归。当多个预测因子含有非0系数并且呈现正态分布的时候，岭回归是理想的回归方法，岭回归对每一个预测因子影响小的模型尤其适用，并且它防止 线性回归模型系数由于共线性而导致无法模拟和高变异性。岭回归对共线预测因子的系数收缩并趋于零，例如,给出k相同的预测因子,都将获得相同的系数等于1/k的单个因子的回归系数值。因此,岭回归不会使某个因子消失，不能将某些因子摒弃来获得最优预测数据集。岭回归(2)估计解决回归问题(1)使用惩罚最小二乘法:

y＝μl_n+Xβ+e₁              (1) 

y＝(y₁，…，y_n)^T其中是观察表型的向量,l_n是一个n维列向量的，μ是一种常见的截距,是n×p矩阵的表示,β表示回归系数的向量,e₁是残差的向量和是残差的误差。

$\widehat{\mathrm{\β}}\left(\mathrm{ridge}\right)=\begin{array}{c}\arg \min \\ \mathrm{\β}\end{array}{\left|\right|y-\mathrm{X\β}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|\mathrm{\β}\left|\right|}_2^2---\left(2\right)$

其中

${\left|\right|y-\mathrm{X\β}\left|\right|}_2^2={\mathrm{\Σ}}_{i=1}^n{({\mathrm{\γ}}_i-x_i^T\mathrm{\β})}^2$

是(二次方程式)损失函数(即残差平方和),是X向量的第i个行。

${\left|\right|\mathrm{\β}\left|\right|}_2^2={\mathrm{\Σ}}_{j=1}^p{\mathrm{\β}}_j^2$

是基于β罚分,λ≥0是调优(罚分,正规化,或复杂化)参数,这些参数通过相对重要性决定经验误差和惩罚调节罚分的强度(即线性收缩)。λ值越大,收缩量越大。λ的值依赖于数据,通过数据驱动的方法(交叉验证)进行确定使用。

CpG密度、GC含量以及CpG-OE值三者计算方法分别为：

CpG密度：某一个CpG上下游各100bp范围内CpG个数除于201bp长度得到此CpG位点的CpG密度；

GC含量：某一个CpG上下游各100bp范围内C和G总数除于201bp长度得到此CpG位点的GC含量；

CpG-OE值：CpG上下游各100bp范围内CpG个数乘于210bp，然后除于C和G个数的乘积。

岭回归矫正MB-seq甲基化水平的系统，包括以下模块：

提取模块：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-OE值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信息；所述的相对甲基化信息包括：MB-seq甲基化水平的观测值MB level,甲基化CpG个数MB mCG,MB-seq测序深度MB depth,当前CpG侧翼+/-100bp区域的MB-seq检测到的平均甲基化水平MB back,这200bp范围的甲基化CpG位点总数MB mCG；

建模模块：根据基因组信息和甲基化信息，采用岭回归模型对真实甲基化水平RRBS level和回归参数建立回归模型；

回归模块：利用岭回归理论，并依据提取出来的基因组信息和甲基化信息，对基因组上的胞嘧啶位点进行回归以得到甲基化水平的模块。

优选的，提取模块提取信息后，设置阈值过滤模块，还将提取到的信息进行阈值过滤，过滤低质量碱基和序列，并过滤adapter污染序列。

优选的，在回归模块计算之前，还设置数据训练和测试模块，采用交叉验证评估模型进行数据训练和测试：

a).将预测特征变量和真实的甲基化水平分成训练和测试数据集；随机抽取50％的CpG位点作为训练数据,剩下的50％作为测试数据；

b).先使用训练数据训练模型；再计算预测甲基化水平值和RRBS测量的甲基化水平值之间的相关性系数；这个过程重复N次,N次的平均相关性系数用来表示模型的预测精度；

对于每个基因组元件，单独进行训练和岭回归测试；而对同时位于多个基因组元件的CpG位点，；取多个预测值的平均值；

c).甲基化水平的预测是全基因组范围的，并且对于RRBS原本就覆盖的位点，采取RRBS的观测值作为最终的甲基化水平；所有未被RRBS覆盖的CpG位点，一律认为其未被甲基化，并且不用于岭回归,甲基 化水平预测值小于0或者大于分别基于岭回归的原则规整到0和1。

优选的，N≥1000。

数据训练和测试模块对模型数据训练时：

a).当变量间存在共线性的时候，通过引入lambda表达式以解决最小二乘回归得到的系数不稳定，方差很大的问题；

b).当模型包含常数项时，岭回归函数对y进行中心化，以y的均值作为因子；对x进行中心化和归一化，以x中各个变量的均值和标准差作为因子；这样对x和y处理后，x和y的均值为0，这使得回归平面经过原点，即常数项为0；

c).当模型不包含常数项时，因为要强制通过原点，该模型假设各个变量的均值为0，因此不对x和y进行中心化，但是对x进行归一化，而且归一化因子也是假设变量均值为0计算出来的该变量的标准差。

数据训练和测试模块对模型进行测试的时候，需要首先对x和y进行中心化和归一化，此时因子是使用训练模型时候进行中心化和归一化的因子，然后再与系数相乘得到预测结果。

在回归模块计算之后，设置异常点处理模块，用于对异常点处理：

1)将MB-seq检测深度为0的位点定义为甲基化水平为0；

2)结合MB-seq甲基化水平的观测值(MB level),甲基化CpG个数(MB mCG),MB-seq测序深度(MB depth),当前CpG侧翼+/-100bp区域的MB-seq检测到的平均甲基化水平(MB back),这200bp范围的甲基化CpG位点总数(MB mCG)，以及每一个CpG位点上下游100bp的基因组CpG密度、GC含量,CpG-OE值等对甲基化水平检测的影响，利用岭回归导入到模型中，并且机器学习得到某一胞嘧啶位点甲基化水平；

3)将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回归的甲基化水平值小于0的位点自动归为甲基化水平为0。

Claims (8)

1.基于岭回归矫正MB-seq甲基化水平的方法，其特征在于，包括以下步骤：

(1)提取信息

(2)建模

(3)岭回归计算；

其中，所述的步骤(1)需要提取的信息有：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-OE值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信息；

所述的步骤(2)建模如下：

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\underset{j=0}{\overset{p}{\mathrm{\Σ}}}{w}_j{x}_{\mathrm{ij}})}^2+\mathrm{\λ}\underset{j=0}{\overset{p}{\mathrm{\Σ}}}{w}_j^2$

其中：

y:目标函数；为以RRBS高通量测序数据唯一比对结果中提取到的覆盖到的每个胞嘧啶的绝对甲基化信息；

x:回归变量矩阵；包括行、列；每行代表每个CpG变量；每列分别为每个变量的CpG密度、GC含量、CpG-OE值以及相对甲基化信息；

所述的相对甲基化信息包括：MB-seq甲基化水平的观测值MBlevel,甲基化CpG个数MB mCG,MB-seq测序深度MB depth,当前CpG侧翼+/-100bp区域的MB-seq检测到的平均甲基化水平MB back,这200bp范围的甲基化CpG位点总数MB mCG；

所述的步骤(3)岭回归计算具体是

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\underset{j=0}{\overset{p}{\mathrm{\Σ}}}{w}_j{x}_{\mathrm{ij}})}^2+\mathrm{\λ}\underset{j=0}{\overset{p}{\mathrm{\Σ}}}{w}_j^2$

对求导，结果为

2X^T(Y-XW)-2λW

令其为0，求得的值：

$\hat{w}={(X^{T}X+\mathrm{\λI})}^{-1}{X}^{T}Y$

输入新的回归变量矩阵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位点，一律认为其未被甲基化，并且不用于岭回归,甲基化水平预测值小于0或者大于分别基于岭回归的原则规整到0和1。

4.如权利要求3所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：N≥1000。

5.如权利要求3所述的基于岭回归矫正MB-seq甲基化水平的方法，其特征在于：在模型数据训练时：

a).当变量间存在共线性的时候，通过引入lambda表达式以解决最小二乘回归得到的系数不稳定，方差很大的问题；

b).当模型包含常数项时，岭回归函数对y进行中心化，以y的均值作为因子；对x进行中心化和归一化，以x中各个变量的均值和标准差作为因子；这样对x和y处理后，x和y的均值为0，这使得回归平面经过原点，即常数项为0；

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侧翼+/-100bp区域的MB-seq检测到的平均甲基化水平(MB back),这200bp范围的甲基化CpG位点总数(MB mCG)，以及每一个CpG位点上下游100bp的基因组CpG密度、GC含量,CpG-OE值等对甲基化水平检测的影响，利用岭回归导入到模型中，并且机器学习得到某一胞嘧啶位点甲基化水平；

3)将回归得到的甲基化水平超过1的位点自动归为甲基化水平为1，而回归的甲基化水平值小于0的位点自动归为甲基化水平为0。

8.岭回归矫正MB-seq甲基化水平的系统，其特征在于：包括以下模块：

提取模块：从参考基因组序列中提取基因组CpG密度、GC含量和CpG-OE值；从MB-seq高通量测序数据唯一比对结果中，提取已知基因组上每个胞嘧啶的相对甲基化信息；从RRBS高通量测序数据唯一比对结果中，提取覆盖到的每个胞嘧啶的绝对甲基化信息；所述的相对甲基化信息包括：MB-seq甲基化水平的观测值MB level,甲基化CpG个数MB mCG,MB-seq测序深度MB depth,当前CpG侧翼+/-100bp区域的MB-seq检测到的平均甲基化水平MB back,这200bp范围的甲基化CpG位点总数MB mCG；

建模模块：根据基因组信息和甲基化信息，采用岭回归模型对真实甲基化水平RRBS level和回归参数建立回归模型；

回归模块：利用岭回归理论，并依据提取出来的基因组信息和甲基化信息，对基因组上的胞嘧啶位点进行回归以得到甲基化水平的模块。