一种DNA甲基化芯片数据的扩展方法与流程

本发明属于DNA甲基化检测芯片的扩展
技术领域：
，更为具体地讲，涉及一种DNA甲基化芯片数据的扩展方法。
背景技术：
：作为人类基因组最为典型的表观遗传现象，DNA甲基化在多种关键生理活动中扮演重要角色。其甲基化状态与各种疾病，特别是癌症的发生密切相关。DNA甲基化检测方法最常用的方法是亚硫酸盐测序、亚硫酸氢微阵列和基于浓缩的方法。亚硫酸盐测序拥有最全面的基因组覆盖，但高测序深度使它变的非常昂贵。虽然亚硫酸氢微阵列和基础的浓缩的方法的成本相对较低，但是基于亚硫酸氢微阵列的平台被先验目标区域所约束，基于浓缩的方法具有相对低的分辨率和高敏感性实验偏差。在450K芯片检测已经广泛用于人体组织，尤其是在几千患者样本的DNA甲基化检测中，因此将450K芯片预测的覆盖范围的扩展将是非常有效的最大化使用现有的数据的方案。技术实现要素：本发明的目的在于克服现有技术的不足，提供一种DNA甲基化芯片数据的扩展方法，通过DNA甲基化芯片显著的扩大检测的覆盖范围，有效的最大化使用现有的数据。为实现上述发明目的，本发明一种DNA甲基化芯片数据的扩展方法，其特征在于，包括以下步骤：(1)、数据获取从公共数据库中获取现有的任意一种组织的每一个待预测CpG位点的上游和下游d1长度范围内基于DNA甲基化芯片测得的甲基化值1及此组织的DNA序列，以及基于全基因组亚硫酸氢钠测序法测得的每一个待预测CpG位点甲基化值2；(2)、整合DNA甲基化芯片测得的甲基化值1及此组织的DNA序列作为预测模型的特征对甲基化值1进行加权，并将此加权值作为预测模型的特征1；在DNA序列中，统计出每一个待预测CpG位点相邻区域d2长度范围内的1聚体和2聚体的DNA碱基对出现次数，作为预测模型的特征2；统计NpN的产生频率，作为预测模型的特征3；最后将特征1、特征2和特征3共同作为预测模型的输入特征；(3)、特征变换(3.1)、对输入特征进行标准化处理；设输入特征是一个n行p列的矩阵X＝(xi,j)n×p进行标准化变换：其中，表示第j列的平均值，sj表示第j列的标准差，对输入特征进行标准化处理后得到标准化矩阵Z；(3.2)、计算标准化矩阵Z的相关系数矩阵R：(3.3)、按照表示预设阈值，计算相关系数矩阵R的特征方程|R-λIp|＝0的特征根λk，Ip单位矩阵，确定出m个主成分分量；对每个特征根λk解方程Rbk＝λkbk，求得单位特征向量bk；(3.4)、将标准化矩阵Z转换为主成分Uk＝ZTbk，k＝1,2,…,m，Uk表示第k个主成分；(4)、构建随机森林回归模型先利用m个主成分分量和甲基化值2构建多个决策树，再将多个决策树按照其产生方式组合成随机森林回归模型；(5)、独立样本预测(5.1)、按照步骤(1)所述方法获取任意组织的每一个待预测CpG位点的上游和下游d1长度范围内基于DNA甲基化芯片测得的甲基化值及此组织的DNA序列；(5.2)、按照步骤(2)-(3)的方法整合DNA甲基化芯片测得的甲基化值及此组织的DNA序列得到输入特征，然后对输入特征进行特征变换得到m个主成分；(5.3)、将m个主成分输入步骤(4)训练好的随机森林模型进行预测得到每个待预测CpG位点的甲基化值即完成对DNA甲基化芯片数据的扩展。本发明的发明目的是这样实现的：本发明一种DNA甲基化芯片数据的扩展方法，通过预测DNA甲基化芯片未覆盖的CpG位点来实现DNA甲基化芯片数据的扩展。具体的讲，对于待预测的CpG位点，先基于DNA甲基化芯片测得的数据和DNA序列提取特征，然后进行特征变换并结合全基因组亚硫酸氢钠测序法测得的待预测点的甲基化值训练随机森林回归模型，最后使用训练好的随机森林回归模型对新数据进行预测。同时，本发明一种DNA甲基化芯片数据的扩展方法还具有以下有益效果：(1)、本发明充分利用了现有芯片数据，对于挖掘出与疾病相关的重要信息具有重大意义。(2)、我们的模型取得了令人满意的性能，它优于现有的常用方法EpiGRAPH。此外，在新的细胞类型的甲基化预测中我们展示了预测模型的泛化能力。由于公开发表的DNA甲基化芯片检测数据较多，因此该模型可以实现对多个组织的全基因组甲基化水平的预测。附图说明图1是一种DNA甲基化芯片数据的扩展方法流程图；图2是预测模型特征选择简图；图3是预测模型的相关系数及一致性示意图；图4是预测模型在H9上的预测结果图。具体实施方式下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。实施例为了方便描述，先对具体实施方式中出现的相关专业术语进行说明：CpG位点：CG碱基同时出现的位点；450K甲基化芯片：一种DNA甲基化芯片；bp：用来表示DNA长度的单位,也就是我们常说的碱基数；pearson相关系数：皮尔森相关系数；Spearman相关系数：斯皮尔曼相关系数；1聚体：物质单一状态时具有特定的性质或功能,不随多少改变，这里指A、C、G、T四种碱基；2聚体：相同或同一种类的物质,以成双的型态出现,可能具有单一状态时没有的性质或功能，这里指A/C/G/T两两出现的组合；NpN：与CpG相同的意思，N＝A/C/G/T，这里指A/C/G/T两两出现的组合；图1是一种DNA甲基化芯片数据的扩展方法流程图。在本实施例中，采用了450K甲基化芯片数据，并结合CpG位点周围500bp的序列中选取1聚体和2聚体的碱基特征和NpN比率，对全基因组的CpG位点甲基化水平进行了预测。如图1所示，本发明一种450K甲基化芯片数据的扩展方法，具体包括以下步骤：(1)、数据获取从公共数据库中获取胚胎干细胞(embryonicstemcells，ESC)中的H1的每一个待预测CpG位点的上游和下游d1＝1000bp长度范围内基于450K甲基化芯片测得的甲基化值1及此组织的DNA序列，以及基于全基因组亚硫酸氢钠测序法测得的每一个待预测CpG位点甲基化值2；图2(a)所示，待预测位点表示为空心点，450K甲基化芯片测得的位点为实心点；其中，每个CpG位点的甲基化水平表示为甲基化率；每个CpG位点的的甲基化值被描述成一个beta值从0(非甲基化)到1(完全甲基化)。(2)、整合450K甲基化芯片测得的甲基化值1及此组织的DNA序列作为预测模型的特征对甲基化值1进行加权，并将此加权值作为预测模型的特征1；在DNA序列中，统计出每一个待预测CpG位点相邻区域d2＝500bp长度范围内的1聚体和2聚体的DNA碱基对出现次数，作为预测模型的特征2；统计NpN的产生频率，作为预测模型的特征3；最后将特征1、特征2和特征3共同作为预测模型的输入特征；在本实施例中，对于每一个的CpG位点，收集相邻区域中1000个碱基对的450K芯片检测数据，将450K芯片检测的此1000个碱基对中的CpG位点的加权甲基化值作为一个特征；将一个CpG位点相邻区域中500bp范围内的1聚体和2聚体的DNA碱基对共20个特征，以及NpN(N＝A/C/G/T)共16个特征，如图2(b)所示，最后共同构成37个特征作为模型构建的输入特征。(3)、特征变换(3.1)、对输入特征进行标准化处理；设输入特征是一个n行p列的矩阵X＝(xi,j)n×p进行标准化变换：其中，表示第j列的平均值，sj表示第j列的标准差，对输入特征进行标准化处理后得到标准化矩阵Z；(3.2)、计算标准化矩阵Z的相关系数矩阵R：(3.3)、按照计算相关系数矩阵R的特征方程|R-λIp|＝0的特征根λk，Ip单位矩阵，确定出15个主成分分量；对每个特征根λk解方程Rbk＝λkbk，求得单位特征向量bk；(3.4)、将标准化矩阵Z转换为主成分Uk＝ZTbk，k＝1,2,…,15，Uk表示第k个主成分；(4)、构建随机森林回归模型先利用m个主成分分量和甲基化值2构建多个决策树，再将多个决策树按照其产生方式组合成随机森林回归模型；下面进行详细说明：(4.1)、构建决策树设训练样本向量的维度是F维，即训练集有F个属性。在构建开始之前选定一个参数f，满足f＜＜F，在构建每个内部节点的过程中，都需要从训练集中采用随机抽样的方法从他的所有F个属性选取f个属性，然后从f个属性中根据信息增益比，选出一个最优的属性充当分裂属性，进而是决策在此节点产生分裂。信息增益比的计算采用如下公式：其中S为全部样本集合，value(T)是属性T所有取值的集合，v是T的其中一个属性值，Sv是S中属性T的值为V的样例集合，|Sv|为Sv中所含样例数。Entropy(Sv)即表示信息增益，他的计算采用如下公式：其中，就是类别的总数，类别C是变量，它的取值是C1,C2,...,Cn，而每一个类别出现的概率分别是P(C1),P(C2),...,P(Cn)。(4.2)构建随机森林回归模型随机森林其实是由很多决策树组合而成的，但是其回归模型的性能往往依赖于组合是按照一种什么样的准则进行的，每一棵决策树的样本都取自一个训练集，都依赖于一个其产生方式的规则，这个规则往往用一个随机向量来表示。它的产生方式的递归描述如下：首先为第q棵决策树生成一个决定了其生成过程的随机向量θq，θq需要满足与之前生成的q-1个随机向量独立同分布。然后在原始训练样本X中利用随机抽样方法进行抽样，这棵决策树的生成过程中的数据都取自这次抽样的结果Xq。采用上述策略，随机向量θq和抽样结果Xq就能够生成第q棵决策树h(Xq,θq)。在对样本集经过q轮抽样之后，一共可以得到q棵决策树。当输入一个新样本时，随机森林中的每一颗决策树分别进行判断，最后对每颗决策树的结果取平均值。(4.3)、模型性能的评估；交叉检验法是十分常用的模型验证方法。其原理是将训练样本分成容量相同的w个子集，并对模型训练w次。在第u次(u＝1,2,L,w)训练时,要用除了第u个子集的所有子集训练模型，再用得到的模型对第u个子集计算误差。以w次误差的平均数值作为模型推广能力的近似数值。对于预测模型性能指标我们采用相关系数(Spearman相关系数和Spearman相关系数)、一致性、特异性(SP)、灵敏度(SE)、准确性(ACC)和马修相关系数(MCC)来进行评估。两个变量之间的Pearson相关系数定义为这两个变量的协方差与二者标准差积的商，即：其中，和μY分别是和Y的平均值，和σY分别是和Y的标准差，和Y分别代表拟合的甲基化值和实际WGBS记录的甲基化值。Spearman相关系数的计算公式为：其中，n为样本集大小，n行甲基化值和经过等级排序后的值为一致性是预测值与实际值之间的差值小于0.25的百分比。SE，SP，ACC和MCC的计算公式如下：其中TP表示预测正确的正样本(true-positive)；TN表示预测正确的负样本(true-negative)；FP表示预测错误的负样本(false-positive)；FN表示预测错误的正样本(false-negative)。通过使用3倍交叉验证测试20次，获得预测模型在每条染色体上的平均性能。预测值和真实值的Pearson相关系数和Spearman相关系数示于图3(a)，一致性示于图3(b)。在图3(a)和3(b)中结果最好的是在18号染色体上，相关系数为0.91(0.82)，一致性为88％；结果最差的是在Y染色体上，相关系数为0.84(0.73)，一致性为82％。(4.4)、方法对比；目前预测甲基化水平最常用的方法是EpiGRAPH，用我们的模型与其进行对比。把预测值经过二值化处理，当CpG位点的预测的甲基化值大于0.5的时候，我们认为的其甲基化状态为+1，反之当预测值小于0.5时，我们认为其甲基化之为-1。在结果最好与最坏的染色体上，用我们的模型与现有常用方法对比。如表1所示，我们构建的模型在18号染色体上，结果优于现有常用方法EpiGRAPH。如表2所示，我们构建的模型在Y染色体上，结果也优于现有常用方法EpiGRAPH。表1是预测模型在18号染色体上与现有方法EpiGRAPH的对比结果；Chr18SESPACCMCCRFinEpiGRAPH0.900.820.860.73RFinourmodel0.940.920.930.86表1表2是预测模型在Y染色体上与现有方法EpiGRAPH的对比结果；ChrYSESPACCMCCRFinEpiGRAPH0.960.260.820.33RFinourmodel0.920.730.850.74表2(5)、独立样本预测(5.1)、按照步骤(1)所述方法获取H9组织的每一个待预测CpG位点的上游和下游d1＝1000bp长度范围内基于DNA甲基化芯片测得的甲基化值及此组织的DNA序列；(5.2)、按照步骤(2)-(3)的方法整合450K甲基化芯片测得的甲基化值及此组织的DNA序列得到输入特征，然后对输入特征进行特征变换得到15个主成分；(5.3)将15个主成分输入步骤(4)训练好的随机森林模型进行预测得到每个待预测CpG位点的甲基化值即完成对DNA甲基化芯片数据的扩展。为了验证结果，下载对应的基于全基因组亚硫酸氢钠测序法测得的每一个待预测CpG位点甲基化值，在H9组织上预测的性能指标如图4所示，除了X染色体以外，预测值和真实值的Pearson相关系数为0.80±0.01，如图4(a)所示一致性为88％±1％，如图4(b)所示。预测的甲基化值在经过二进制处理之后，我们计算该预测结果的SE，SP，ACC和MCC,结果如图4(c)所示，得到的结果是令人满意的。尽管上面对本发明说明性的具体实施方式进行了描述，以便于本
技术领域：
的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本
技术领域：
的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。当前第1页1 2 3