CN102982253A - 一种多样本间甲基化差异检测方法及装置 - Google Patents

Abstract

本发明提供了一种多样本间甲基化差异的检测方法，该方法包括步骤：获得多组样本的基因组测序片段对应于参考基因组序列上的位置信息以及各个样本的甲基化数据；对在参考基因组序列上滑动的同一判断区间内的各组的甲基化数据进行分析，获取上述组基因组上存在甲基化差异的甲基化差异区域；对同一甲基化差异区域内的各组的甲基化数据进行分析，获取在甲基化差异区域内存在甲基化差异的具体组。本发明还提供了一种多样本间甲基化差异的检测装置。本发明的多样本间甲基化差异的检测方法及装置能查找并检测出多个样本基因组之间存在的甲基化差异性区域以及在该甲基化差异性区域存在甲基化差异的样本组。

Description

一种多样本间甲基化差异检测方法及装置

技术领域

本发明涉及表观遗传学领域，特别是涉及一种多样本间甲基化差异检测方法及装置。

背景技术

DNA甲基化已经成为表观遗传学和表观基因组学的重要研究内容。DNA甲基化是重要的表观遗传调控因子之一，对不同细胞、组织等甲基化修饰模式差异的研究，对于更好地解释组织与组织之间、个体与个体之间在表观修饰上存在差异的原因以及对于疾病易感人群，个体医疗甚至药物设计等方面的研究有着重大的意义。然而，相对于测定DNA甲基化谱的高通量实验技术的快速发展，从这些实验数据中查找与检测甲基化显著性差异区域的方法和装置的步伐却远远滞后。

目前已实现的甲基化显著性差异区域检测的方法，主要集中在两个样本之间的比较，例如，利用卡方检验或者t检验获得显著性差异区间，这种方法在一定程度上可以获得具有差异的区间，但是由于灵敏度较低，统计功效不高，不能实现多样本间甲基化显著性差异区域的查找与检测。另外就是利用基因芯片技术，但该技术具有检测结果不准确、实验要求高、费用昂贵的缺点。

因此，研究一种新的可以进行多样本间甲基化差异区域查找和甲基化差异检测并且具有较高灵敏度、特异度、准确性高，成本低的技术就成了亟待解决的问题。

发明内容

本发明主要解决的技术问题是提供一种多样本间甲基化差异检测方法及装置，能够准确、灵敏、快速地进行多样本间甲基化差异检测。

为解决上述技术问题，本发明采用的一个技术方案是：本发明提供了一种多样本间甲基化差异的检测方法，该方法包括步骤：获得多组样本的基因组测序片段对应于参考基因组序列上的位置信息以及各个样本的甲基化数据；对在参考基因组序列上滑动的同一判断区间内的各组的甲基化数据进行分析，获取上述组基因组上存在甲基化差异的甲基化差异区域；对同一甲基化差异区域内的各组的甲基化数据进行分析，获取在甲基化差异区域内存在甲基化差异的具体组。

根据本发明的一优选实施例，对在参考基因组序列上滑动的同一判断区间内的各组的甲基化数据进行分析，获取上述组基因组上存在甲基化差异的甲基化差异区域的步骤包括以下步骤：预设窗口长度，根据位置信息，从参考基因组序列起始端开始，以窗口长度作为判断区间的起始长度在参考基因组序列上设定判断区间；判断同一判断区间内的各组的甲基化数据是否符合方差分析前提；若符合方差分析前提，则将位置在同一判断区间内的各组的甲基化数据进行组间方差分析，若不符合方差分析前提，则将位置在同一判断区间内的各组的甲基化数据进行组间Kruskal-Wallis非参数检验，得到分析结果，根据分析结果判断上述组在该判断区间内是否存在甲基化差异；若存在甲基化差异，则延长该判断区间，判断位置在同一延长后的判断区间内的各组的甲基化数据是否符合方差分析前提，若符合则进行前述组间方差分析，若不符合则进行前述组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在当前判断区间内是否存在甲基化差异的步骤，重复本步骤直至判断得出上述组在当前判断区间内不存在甲基化差异，并输出该当前判断区间的信息作为甲基化差异区域；若不存在甲基化差异，则从上一判断区间的末端开始在参考基因组序列上以窗口长度作为判断区间的起始长度设定下一判断区间，判断位置在同一判断区间内的各组的甲基化数据是否符合方差分析前提，若符合则进行前述组间方差分析，若不符合则进行前述组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在该判断区间内是否存在甲基化差异，若存在甲基化差异则执行前述延长该判断区间、判断是否符合方差分析前提以及分析并判断上述组在当前判断区间内是否存在甲基化差异的步骤，若不存在甲基化差异则执行本步骤。

根据本发明的一优选实施例，的对同一甲基化差异区域内的各组的甲基化数据进行分析，获取在甲基化差异区域内存在甲基化差异的具体组的步骤包括以下步骤：通过最小显著性差异法对位置在同一甲基化差异区域内的且符合方差分析前提的各组的甲基化数据进行组间两两比较，通过Kruskal-Wallis Dunn方法对各组在甲基化差异区域内的且不符合方差分析前提的甲基化数据进行组间两两比较，得到比较结果；根据比较结果判断并输出在该甲基化差异区域内存在甲基化差异的具体组。

根据本发明的一优选实施例，对在参考基因组序列上滑动的同一判断区间内的各组的甲基化数据进行分析，获取上述组基因组上存在甲基化差异的甲基化差异区域的步骤包括以下步骤：在判断完上述组在当前判断区间是否存在甲基化差异之后，判断当前判断区间是否已达参考基因组序列末端，若是，则终止继续设定或延长判断区间以及分析并判断上述组在判断区间内是否存在甲基化差异的步骤。

根据本发明的一优选实施例，在判断得出上述组在当前判断区间内存在甲基化差异时，延长该判断区间的步骤具体为：将该判断区间延长一个预设步长。

根据本发明的一优选实施例，方差分析前提为：各组样本的数据是否具有独立性；同一判断区间内的各组的甲基化数据是否符合正态分布；同一判断区间内的各组的甲基化数据是否符合方差齐性。

本发明还提供了一种多样本间甲基化差异的检测装置，该装置包括：甲基化数据获取单元，用于获得多组样本的基因组测序片段对应于参考基因组序列上的位置信息以及各个样本的甲基化数据；甲基化差异区域获取单元，用于对在参考基因组序列上滑动的同一判断区间内的各组的甲基化数据进行分析，获取上述组基因组上存在甲基化差异的甲基化差异区域；甲基化差异组获取单元，用于对同一甲基化差异区域内的各组的甲基化数据进行分析，获取在甲基化差异区域内存在甲基化差异的具体组。

根据本发明的一优选实施例，甲基化差异区域获取单元包括：设置单元，用于预设窗口长度，根据位置信息，从参考基因组序列起始端开始，以窗口长度作为判断区间的起始长度在参考基因组序列上设定判断区间；判断单元，用于判断同一判断区间内的各组的甲基化数据是否符合方差分析前提；分析单元，若符合方差分析前提，分析单元用于将位置在同一判断区间内的各组的甲基化数据进行组间方差分析，若不符合方差分析前提，分析单元用于将位置在同一判断区间内的各组的甲基化数据进行组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在该判断区间内是否存在甲基化差异；控制单元，若存在甲基化差异，控制单元控制设置单元延长该判断区间，控制判断单元判断位置在同一延长后的判断区间内的各组的甲基化数据是否符合方差分析前提，控制分析单元在符合方差分析前提时进行前述组间方差分析，在不符合方差分析前提时进行前述组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在当前判断区间内是否存在甲基化差异的步骤，控制设置单元、判断单元及分析单元重复本步骤直至判断得出上述组在当前判断区间内不存在甲基化差异，并输出该当前判断区间的信息作为甲基化差异区域；若不存在甲基化差异，控制单元控制设置单元从上一判断区间的末端开始在参考基因组序列上以窗口长度作为判断区间的起始长度设定下一判断区间，控制判断单元判断位置在同一判断区间内的各组的甲基化数据是否符合方差分析前提，控制分析单元在符合方差分析前提时进行前述组间方差分析，控制分析单元在不符合方差分析前提时进行前述组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在该判断区间内是否存在甲基化差异的步骤，若存在甲基化差异则执行前述延长该判断区间、判断是否符合方差分析前提以及分析并判断上述组在当前判断区间内是否存在甲基化差异的步骤，若不存在甲基化差异则执行本步骤。

根据本发明的一优选实施例，甲基化差异组获取单元包括：比较单元，用于通过最小显著性差异法对位置在同一甲基化差异区域内的且符合方差分析前提的各组的甲基化数据进行组间两两比较，通过Kruskal-Wallis Dunn方法对各组在甲基化差异区域内的且不符合方差分析前提的甲基化数据进行组间两两比较，得到比较结果；输出单元，用于根据比较结果判断并输出在该甲基化差异区域内存在甲基化差异的具体组。

根据本发明的一优选实施例，甲基化差异区域获取单元还包括：终止单元，用于在分析单元判断完上述组在当前判断区间是否存在甲基化差异之后，判断当前判断区间是否已达参考基因组序列末端，若是，则终止设置单元、判断单元及分析单元。

根据本发明的一优选实施例，在分析单元判断得出上述组在当前判断区间内存在甲基化差异时，设置单元延长该判断区间具体为：将该判断区间延长一个预设步长。

根据本发明的一优选实施例，判断单元判断同一判断区间内的各组的甲基化数据是否符合方差分析前提具体为判断：各组样本的数据是否具有独立性；同一判断区间内的各组的甲基化数据是否符合正态分布；同一判断区间内的各组的甲基化数据是否符合方差齐性。

区别于现有甲基化显著性差异区域检测技术灵敏度较低、统计功效不高、结果不准确、实验要求高、费用昂贵的情况，本发明的多样本间甲基化差异检测方法及装置具有高灵敏度、特异性和快速的特点，能在解决大批量样本或者海量测序数据背景下准确查找并检测出样本基因组之间存在的甲基化差异性区域以及在该甲基化差异性区域存在甲基化差异的样本组，为进一步在表观遗传学方面的生物信息挖掘与研究提供基础，同时也为表观生物靶标，甚至药物设计等方面的研究提供一套启发性的方法。

附图说明

图1是本发明多样本间甲基化差异的检测方法一实施例的流程图；

图2是图1的多样本间甲基化差异的检测方法中的步骤S14的具体实现步骤的流程图；

图3是图1的多样本间甲基化差异的检测方法中的步骤S15的具体实现步骤的流程图；

图4是本发明多样本间甲基化差异的检测方法另一实施例的流程图；

图5是本发明多样本间甲基化差异的检测方法另一实施例的流程图；

图6是本发明多样本间甲基化差异的检测方法一个具体实施方式的流程图；

图7、图8是图6的多样本间甲基化差异的检测方法检测到的长度较短(≤600碱基)的甲基化显著性差异区域；

图9是图6的多样本间甲基化差异的检测方法检测到的中等长度(600-1500碱基)的甲基化显著性差异区域；

图10、图11是图6的多样本间甲基化差异的检测方法检测到的长度较长(＞1500碱基)的甲基化显著性差异区域；

图12是本发明提供的多样本间甲基化差异的检测装置一实施例的结构示意图。

具体实施方式

下面结合附图和实施例对本发明进行详细说明。

图1示出了本发明提供的多样本间甲基化差异的检测方法一实施例的流程图。本实施例中采用的测序技术为高通量测序技术，高通量测序技术可以为Illumina GA测序技术，也可以是现有的其它高通量测序技术。参考基因组序列可取于公共数据库(例如，美国国立生物技术信息中心，NCBI，该公共数据库可以通过如下途径获取http://www.ncbi.nlm.nih.gov/gene？term＝hvp)。

如图1所示，该多样本间甲基化差异的检测方法10包括步骤：

在步骤S13中，获得多组样本的基因组测序片段对应于参考基因组序列上的位置信息以及各个样本的甲基化数据。

当测序技术采用高通量测序技术时，可以通过任何一种短序列映射程序(如SOAP等映射程序)将高通量测序技术得到的样本的基因组测序片段比对到参考基因组序列上。根据比对的结果，获得了各个样本的基因组测序片段对应于参考基因组序列上的位置信息。

根据比对的结果，获得支持样本基因组测序片段中各个胞嘧啶的序列数量(对于MeDIP技术，则是各个胞嘧啶被覆盖的序列数量)。若采用MeDIP技术，则确定基因组测序片段对应于参考基因组序列上的具体位置，计算各个样本的基因组测序片段比对到参考基因组序列上后对参考基因组序列上各个胞嘧啶的深度覆盖情况，即直接计算覆盖在参考基因组序列上的各个胞嘧啶的测序片段的条数，并对计算获得的各个样本的深度覆盖情况进行标准化，以此来获得不同样本的甲基化数据。标准化的公式：sample_read_num_STD＝sample_read_num_ori *(sample_max_read_num/max_read_num)，其中，sample_read_num_STD表示标准化后样本的序列数；sample_read_num_ori表示样本实际序列数；sample_max_read_num表示样本对应文库的测序下机序列数；max_read_num表示待研究的各个样本文库中最大的测序下机序列数。标准化之后按照该公式重新计算各个样本胞嘧啶被覆盖的序列数量。

在步骤S14中，对在参考基因组序列上滑动的同一判断区间内的各组的甲基化数据进行分析，获取上述组基因组上存在甲基化差异的甲基化差异区域。图2是图1的多样本间甲基化差异的检测方法中的步骤S14的具体实现步骤的流程图。如图2所示，在本实施例中，步骤S14具体采用以下步骤实现：

步骤S141，预设窗口长度和预设步长，根据位置信息，从参考基因组序列起始端开始，以窗口长度作为判断区间的起始长度在参考基因组序列上设定判断区间，将预设步长作为判断区间的延长长度。

由于不同物种的甲基化差异情形不一定相同，因此窗口长度和步长的具体数值需要根据具体情况由用户决定的，通常初始窗口长度不会超过1000个碱基距离，步长最少长度为1个碱基长。

步骤S142，判断同一判断区间内的各组的甲基化数据是否符合方差分析前提。其中，方差分析前提具体为：各组样本的数据是否具有独立性；同一判断区间内的各组的甲基化数据是否符合正态分布；同一判断区间内的各组的甲基化数据是否符合方差齐性。

步骤S143，若符合方差分析前提则将位置在同一判断区间内的各组的甲基化数据进行组间方差分析，若不符合方差分析前提则将位置在同一判断区间内的各组的甲基化数据进行组间Kruskal-Wallis非参数检验，得到分析结果。

步骤S144，根据分析结果判断上述组在该判断区间内是否存在甲基化差异。

步骤S145，若存在甲基化差异，则延长该判断区间，判断位置在同一延长后的判断区间内的各组的所述甲基化数据是否符合方差分析前提，然后执行步骤S143及步骤S144，重复本步骤S145直至判断得出上述组在当前判断区间内不存在甲基化差异，并输出该当前判断区间的信息作为甲基化差异区域。

具体的，在判断得出上述组在当前判断区间内存在甲基化差异时延长该判断区间为将该判断区间延长一个预设步长。

步骤S146，若不存在甲基化差异，则从上一判断区间的末端开始在参考基因组序列上以所述窗口长度作为所述判断区间的起始长度设定下一判断区间，判断位置在同一判断区间内的各组的所述甲基化数据是否符合方差分析前提，然后执行步骤S143及步骤S144，若存在甲基化差异则执行前述S145步骤，若不存在甲基化差异则执行本步骤S146。

步骤S147，在判断完上述组在当前判断区间是否存在甲基化差异之后，判断当前判断区间是否已达所述参考基因组序列末端，若是，则终止继续设定或延长判断区间以及分析并判断上述组在判断区间内是否存在甲基化差异的步骤。

通过这整一个过程便可以得到上述组之间存在着甲基化显著差异的一系列区间的初步结果。根据本发明的另一实施例，接着再对这些区间进行相关的FDR(false discovery rate，错误发现率)过滤，最后得到最终的甲基化差异区域。

在步骤S15中，对同一甲基化差异区域内的各组的甲基化数据进行分析，获取在甲基化差异区域内存在甲基化差异的具体组。图3是图1的多样本间甲基化差异的检测方法中的步骤S15的具体实现步骤的流程图。如图3所示，在本实施例中，步骤S15具体采用以下步骤实现：

步骤S151，在获得甲基化差异区域后，若同一甲基化差异区域内的各组的甲基化数据满足方差分析前提则通过最小显著性差异法(Leastsignificant difference，简称LSD)对位置在同一甲基化差异区域内的各组的甲基化数据进行组间两两比较，若不满足方差分析前提则通过Kruskal-Wallis Dunn方法对各组在甲基化差异区域内的甲基化数据进行组间两两比较，得到比较结果。

步骤S152，根据比较结果判断并输出在该甲基化差异区域内存在甲基化差异的具体组。

如上所述，通过多重比较的方法确定了在这些甲基化差异区域中具体是哪些组和哪些组之间存在着差异。

特别的，在步骤S13之前，可以采用如下步骤对样本基因组进行预处理：

步骤Sp1，获得样本的基因组测序片段。具体可以采用MeDIP-seq(Methylated DNA Immunoprecipitation Sequencing，甲基化DNA免疫共沉淀)技术，即通过5′-甲基胞嘧啶抗体特异性富集样本基因组上发生甲基化的DNA片段，然后将这些DNA片段进行测序，获得样本的基因组测序片段。

步骤Sp2，将接头序列测序质量过低的基因组测序片段去除。该步骤可通过如下方式实现：预先设置样本接头序列的测序质量阈值(如5)和碱基数阈值(如3)，将接头序列中碱基的测序质量值低于测序质量阈值，且碱基的数量超过碱基数阈值的序列去除，例如，综合考虑测序条件和环境，将本实施例中10bp(碱基对)的接头序列中测序质量值低于5的碱基且个数大于3个的基因组测序片段去除。

步骤Sp3，将基因组测序片段中的样本接头序列与样本接头序列库进行比对，实现区分样本操作，并同时将样本接头序列从基因组测序片段中去除。具体包括如下步骤：

步骤Sp31，样本接头序列与样本接头序列库中序列进行完全匹配操作。

步骤Sp32，考虑到一系列实验过程中，样本接头序列可能出现降解情况，假设样本接序列降解1-2bp与样本接头序列库中序列对应部分进行完全匹配操作。

步骤Sp33，考虑到一系列实验过程中，样本接头序列发生碱基插入，本发明中允许样本序列仅有一个碱基的插入，在样本接头序列起始端进行完全匹配操作，当出现某碱基无法匹配时认为该碱基为插入碱基，跳过此碱基后继续严格的完全匹配操作。

步骤Sp34，考虑到一系列实验过程中，样本接头序列发生碱基缺失，本发明中允许样本序列仅有一个碱基的缺失，在样本接头序列中允许任何一个位置缺失一个碱基后，进行完全匹配操作。

完成步骤Sp31-Sp34后，按照步骤Sp31＞步骤Sp32＞步骤Sp33＞步骤Sp34的优先级顺序确定最终的样本接头序列的比对结果。而对于四步操作中四步均无比对结果，或者一个步骤同时比对到两个结果或仅有且同时Sp33、Sp34步骤比对出结果，则认为该比对结果是由于无法区分而判定为无效信息，并将相应的该条基因组测序片段去除。比对到同一样本接头序列的被认为是同一样本序列，从而实现样本区分的目的。最后去除每条有效的基因组测序片段中样本接头序列部分(长度范围大概8-11bp)。

图4示出了本发明提供的多样本间甲基化差异的检测方法一实施例的流程图。如图4所示，该多样本间甲基化差异的检测方法20包括步骤：S21、S23、S24、S25。其中步骤S23、S24和S25可以分别执行与图1所示的步骤S13、S14和S15相同或相似的技术内容，为简洁起见，这里不再赘述其技术内容。如图4所示，在步骤S23之前，执行步骤S21：对基因组测序片段进行过滤，以去除不合格的基因组测序片段。具体来说，步骤S21包括如下情况：

预先设置碱基的测序质量阈值和不合格碱基的比例阈值，其中，测序质量阈值和不合格碱基的比例阈值由具体测序技术及测序环境而定，例如，测序质量阈值设置为5，测序质量阈值低于5的碱基为不合格碱基，不合格碱基的比例阈值设置为50％，当基因组测序片段中碱基的测序质量值低于测序质量阈值，且不合格碱基的个数占整条序列碱基个数的比例超过比例阈值时，则认为该基因组测序片段是不合格序列并将其过滤掉；

当基因组测序片段的测序结果中不确定的碱基(如Illumina GA测序结果中的N)的个数超过整条序列碱基个数的10％，则认为该基因组测序片段是不合格序列并将其过滤掉；

与测序接头序列库进行比对，如果基因组测序片段中存在测序接头序列，则认为该基因组测序片段是不合格序列并将其过滤掉；

除样本接头序列外，与其它实验引入的的外源序列比对(如各种接头序列)，若序列中存在外源序列则认为该基因组测序片段是不合格序列并将其过滤掉。

本发明提供的多样本间甲基化差异检测方法，通过对基因组测序片段进行过滤，去除不合格的基因组测序片段，进一步降低了不合格基因组测序片段的影响，从而提高了检测分析的准确性。

图5示出了本发明提供的多样本间甲基化差异的检测方法一实施例的流程图。如图5所示，该多样本间甲基化差异的检测方法30包括步骤：S32、S33、S34、S35。其中步骤S33、S34、S35可以分别执行与图1所示的步骤S13、S14和S15相同或相似的技术内容，为简洁起见，这里不再赘述其技术内容。如图5所示，在步骤S33之前，执行步骤S32：“对多组样本的基因组测序片段与相应的参考基因组序列比对获得的结果进行筛选”。具体来说，步骤S32包括如下情况：

如唯一性的比较、比对长度的比较、错配数的比较，比对次数的比较等，筛选出每条序列比对结果最好及非常接近最好结果的比对信息，选用的筛选条件需视选用的比对软件、序列背景而定。最终仅保留筛选的比对结果为唯一的序列作为有效序列。

在下文的其他实施例中还将举例对前述步骤中的具体实现方式作进一步的详细介绍。

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合具体实施方式对本发明进行进一步详细说明，结果的分析中为简明起见也只用物种的3号染色体作为例子。图6示出了本发明提供的多样本间甲基化差异的检测方法一个具体实施方式的流程图。

样本：来自同种哺乳动物的三个不同品种，每个品种6个正常个体，每个个体提取8类脂肪组织和2类肌肉组织，共180个真实组织样本。

上机策略：由于皆是正常个体，无需设计对照组，只需设计一类文库，包含180个样本，超声打断样本的DNA片段后采用MeDIP技术沉降目标序列片段，然后采用Illumina GA高通量测序技术对这些目标序列片段进行测序，获得样本的基因组测序片段。

如图6所示，该多样本间甲基化差异的检测方法40包括：

步骤S41，接收高通量测序技术得到的基因组测序片段。

接收到基因组测序片段后，需对基因组测序片段进行过滤，以去除不合格的基因组测序片段。不合格的基因组测序片段包括：测序质量值低于5的碱基的个数超过整条序列碱基个数的50％的则认为是不合格序列；序列中测序结果中N的个数超过整条序列碱基个数的10％的则认为是不合格序列。

步骤S42，将基因组测序片段中的样本接头序列与样本接头序列库进行比对，实现区分样本操作，同时将接头序列从基因组测序片段中去除。

该步骤将接头序列中有测序质量低于5的碱基个数大于3个的序列去除，并执行与前述步骤Sp3相同或相似的内容，为简洁起见，此处不再赘述。此过程在序列下机时就已经完成，信息分析过程已经不需要再进行该步骤的区分。

步骤S43，与参考基因组序列进行比对并进行分析，获得各个样本的基因组测序片段对应于参考基因组序列上的位置信息以及各个样本基因组测序片段中胞嘧啶的甲基化数据。

采用SOAP(Short Oligonucleotide Analysis Package)映射程序，将高通量测序技术得到的基因组测序片段比对到参考基因组序列上。比对后，筛选每条基因组测序片段比对得最好的结果，即对于比到相同位置的基因组测序片段只选比对质量最好的，并且比对结果为唯一的基因组测序片段作为最终的有效序列。确定了样本的基因组测序片段在参考基因组序列上的具体位置之后，依照MeDIP技术的特点，计算各个样本基因组测序片段比对到参考基因组序列上后，对参考基因组序列上各个胞嘧啶的深度覆盖情况，并进行标准化，以此来获得不同样本的甲基化数据。标准化的方法请参考前述步骤S13的内容，为简洁起见，此处不再赘述。

步骤S44，对样本进行分组。

按照品种、组织和性别进行七类分析，各个分析中的具体分组如下：

(1)脂肪：按照不同的品种分成3组进行比较，用于检测脂肪组织在品种之间的甲基化显著性差异区域；

(2)脂肪：按照不同的脂肪组织分成8组进行比较，用于检测脂肪组织之间的显著性差异区域；

(3)脂肪：按照不同的性别分2组进行比较，用于检测脂肪组织在性别之间的甲基化显著性差异区域；

(4)肌肉：按照不同的品种分成3组进行比较，用于检测肌肉组织在品种之间的甲基化显著性差异区域；

(5)肌肉：按照不同的肌肉分成2组进行比较，用于检测肌肉组织之间的显著性差异区域；

(6)肌肉：按照不同的性别分2组进行比较，用于检测肌肉组织在性别之间的甲基化显著性差异区域。

(7)将2类肌肉组织分别与8类脂肪组织进行比较，检测肌肉组织与脂肪组织之间甲基化的显著性差异区域；

步骤S45，该步骤执行与前述步骤S141至步骤S147以及步骤S151至步骤S152相同或相似的内容，为简洁起见，此处不再赘述。如此，通过单因素方差分析以及Kruskal-Wallis检验最后便可以获得三个品种之间具有可信度高的甲基化显著性差异区域。

本实施例中只给出品种之间甲基化差异性区域的部分比较结果。

图7-11为第一种分类，即不同品种的脂肪组织之间甲基化显著性差异区域结果，其中图7、图8示出的是多样本间甲基化差异的检测方法400检测到的长度较短(≤600碱基)的甲基化显著性差异区域；图9示出的是多样本间甲基化差异的检测方法400检测到的中等长度(600-1500碱基)的甲基化显著性差异区域；图10、图11示出的是多样本间甲基化差异的检测方法400检测到的长度较长(＞1500碱基)的甲基化显著性差异区域。图中两条竖直虚线之间的区域即是检测到的差异区域，水平虚线则是深度的阈值，本实施方式中要求三个组中至少有一组的平均深度超过该阈值，该阈值为10，横坐标表示染色体上的位置，区间的两端是该差异性区域左右各500碱基长的范围，图中的三条线分别表示三个品种(即三个组)，图中的黑点，黑方框和黑菱形分别代表三个组所处区间中的胞嘧啶(C)和鸟嘌呤(G)二核苷酸位点，这是哺乳动物中主要的甲基化位点。由图可以看出，三个品种在所检测出来的区间相比于区间的两端中存在着明显的差异。由图9、图10和图11同样可以获得这种情形，这说明本发明的多样本间甲基化差异的检测方法400具有着很高敏感度以及特异性，对不同长度范围的差异性区间都具有很高的统计功效和敏感度。

至此，本实施例实现了利用生物信息学的方法高敏感度和特异性，快速地检测并识别各组间甲基化差异区间的目的。

虽然图7-11只对应于第一种分类分析的结果，但其它的分类分析所采用的方法和程序与第一种分类分析采用的方法和程序一样，由于第一种分类分析得出了结果，因此其它的分类分析也都能得出结果，在此对其它的分类分析结果就不一一赘述。

图12示出了本发明提供的多样本间甲基化差异的检测装置一实施例的结构示意图。如图12所示，该多样本间甲基化差异的检测装置50包括甲基化数据获取单元51、甲基化差异区域获取单元52和甲基化差异组获取单元53，其中：

甲基化数据获取单元51，用于获得多组样本的基因组测序片段对应于参考基因组序列上的位置信息以及各个样本的甲基化数据。

甲基化差异区域获取单元52，用于对在所述参考基因组序列上滑动的同一判断区间内的各组的所述甲基化数据进行分析，获取上述组基因组上存在甲基化差异的甲基化差异区域。

甲基化差异组获取单元53，用于对同一甲基化差异区域内的各组的甲基化数据进行分析，获取在甲基化差异区域内存在甲基化差异的具体组。

如图12所示，甲基化差异区域获取单元52包括设置单元521、判断单元522、分析单元523、控制单元524以及终止单元525。其中，

设置单元521，用于预设窗口长度和预设步长，根据位置信息，从参考基因组序列起始端开始，以所述窗口长度作为判断区间的起始长度在参考基因组序列上设定判断区间。

判断单元522，用于判断同一判断区间内的各组的甲基化数据是否符合方差分析前提。其中，方差分析前提的内容如步骤S142中所述，此处不再赘述。

分析单元523，若符合方差分析前提则分析单元522将位置在同一判断区间内的各组的甲基化数据进行组间方差分析，若不若符合方差分析前提则分析单元522将位置在同一判断区间内的各组的甲基化数据进行组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在该判断区间内是否存在甲基化差异。

控制单元524，若存在甲基化差异，控制单元524控制设置单元延长该判断区间，控制判断单元522判断位置在同一延长后的判断区间内的各组的甲基化数据是否符合方差分析前提，控制分析单元523在符合方差分析前提时进行前述组间方差分析，在不符合方差分析前提时进行前述组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在当前判断区间内是否存在甲基化差异的步骤，控制设置单元521、判断单元522及分析单元523重复本步骤直至判断得出上述组在当前判断区间内不存在甲基化差异，并输出该当前判断区间的信息作为甲基化差异区域，

若不存在甲基化差异，控制单元524控制设置单元521从上一判断区间的末端开始在参考基因组序列上以窗口长度作为判断区间的起始长度设定下一判断区间，控制判断单元522判断位置在同一判断区间内的各组的甲基化数据是否符合方差分析前提，控制分析单元523在符合方差分析前提时进行前述组间方差分析，控制分析单元523在不符合方差分析前提时进行前述组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在该判断区间内是否存在甲基化差异的步骤，若存在甲基化差异则执行前述延长该判断区间、判断是否符合方差分析前提以及分析并判断上述组在当前判断区间内是否存在甲基化差异的步骤，若不存在甲基化差异则执行本步骤。

其中，在判断得出上述组在当前判断区间内存在甲基化差异时，设置单元521延长该判断区间具体为：将该判断区间延长一个预设步长。

终止单元523，用于在分析单元523判断完上述组在当前判断区间是否存在甲基化差异之后，判断当前判断区间是否已达参考基因组序列末端，若是，则终止设置单元521、判断单元522及分析单元523。

在本发明提供的多样本间甲基化差异的检测装置另一实施例中，甲基化差异区域获取单元52还包括：过滤单元，用于对获得的甲基化差异区域进行FDR假阳性率过滤。

如图12所示，甲基化差异组获取单元53包括比较单元531及输出单元532。其中：

比较单元531，用于通过最小显著性差异法对位置在同一甲基化差异区域内的且符合方差分析前提的各组的甲基化数据进行组间两两比较，通过Kruskal-Wallis Dunn方法对各组在所述甲基化差异区域内的且不符合方差分析前提的甲基化数据进行组间两两比较，得到比较结果。

输出单元532，根据比较结果判断并输出在该甲基化差异区域内存在甲基化差异的具体组。

区别于现有技术的情况，本发明的多样本间甲基化差异检测方法及装置具有高灵敏度、特异性和快速的特点，能在解决大批量样本或者海量测序数据背景下准确查找并检测出样本基因组之间存在的甲基化差异性区域以及在该甲基化差异性区域存在甲基化差异的样本组，为进一步在表观遗传学方面的生物信息挖掘与研究提供基础，同时也为表观生物靶标，甚至药物设计等方面的研究提供一套启发性的方法。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (12)

1.一种多样本间甲基化差异的检测方法，其特征在于，所述方法包括步骤：

获得多组样本的基因组测序片段对应于参考基因组序列上的位置信息以及各个样本的甲基化数据；

对在所述参考基因组序列上滑动的同一判断区间内的各组的所述甲基化数据进行分析，获取上述组基因组上存在甲基化差异的甲基化差异区域；

对同一甲基化差异区域内的各组的所述甲基化数据进行分析，获取在所述甲基化差异区域内存在甲基化差异的具体组。

2.根据权利要求1所述的多样本间甲基化差异的检测方法，其特征在于，

对在所述参考基因组序列上滑动的同一判断区间内的各组的所述甲基化数据进行分析，获取上述组基因组上存在甲基化差异的甲基化差异区域的步骤包括以下步骤：

预设窗口长度，根据位置信息，从参考基因组序列起始端开始，以所述窗口长度作为判断区间的起始长度在参考基因组序列上设定判断区间；

判断同一判断区间内的各组的所述甲基化数据是否符合方差分析前提；

若符合方差分析前提，则将位置在同一判断区间内的各组的所述甲基化数据进行组间方差分析，若不符合方差分析前提，则将位置在同一判断区间内的各组的所述甲基化数据进行组间Kruskal-Wallis非参数检验，得到分析结果，根据分析结果判断上述组在该判断区间内是否存在甲基化差异；

若存在甲基化差异，则延长该判断区间，判断位置在同一延长后的判断区间内的各组的所述甲基化数据是否符合方差分析前提，若符合则进行前述组间方差分析，若不符合则进行前述组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在当前判断区间内是否存在甲基化差异的步骤，重复本步骤直至判断得出上述组在当前判断区间内不存在甲基化差异，并输出该当前判断区间的信息作为甲基化差异区域；

若不存在甲基化差异，则从上一判断区间的末端开始在参考基因组序列上以所述窗口长度作为所述判断区间的起始长度设定下一判断区间，判断位置在同一判断区间内的各组的所述甲基化数据是否符合方差分析前提，若符合则进行前述组间方差分析，若不符合则进行前述组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在该判断区间内是否存在甲基化差异，若存在甲基化差异则执行前述延长该判断区间、判断是否符合方差分析前提以及分析并判断上述组在当前判断区间内是否存在甲基化差异的步骤，若不存在甲基化差异则执行本步骤。

3.根据权利要求2所述的多样本间甲基化差异的检测方法，其特征在于，

所述的对同一甲基化差异区域内的各组的所述甲基化数据进行分析，获取在所述甲基化差异区域内存在甲基化差异的具体组的步骤包括以下步骤：

通过最小显著性差异法对位置在同一甲基化差异区域内的且符合方差分析前提的各组的所述甲基化数据进行组间两两比较，通过Kruskal-Wallis Dunn方法对各组在所述甲基化差异区域内的且不符合方差分析前提的所述甲基化数据进行组间两两比较，得到比较结果；

根据所述比较结果判断并输出在该甲基化差异区域内存在甲基化差异的具体组。

4.根据权利要求2所述的多样本间甲基化差异的检测方法，其特征在于，

对在所述参考基因组序列上滑动的同一判断区间内的各组的所述甲基化数据进行分析，获取上述组基因组上存在甲基化差异的甲基化差异区域的步骤包括以下步骤：

在判断完上述组在当前判断区间是否存在甲基化差异之后，判断当前判断区间是否已达所述参考基因组序列末端，若是，则终止继续设定或延长判断区间以及分析并判断上述组在判断区间内是否存在甲基化差异的步骤。

5.根据权利要求2所述的多样本间甲基化差异的检测方法，其特征在于，在判断得出上述组在当前判断区间内存在甲基化差异时，所述延长该判断区间的步骤具体为：将该判断区间延长一个预设步长。

6.根据权利要求2所述的多样本间甲基化差异的检测方法，其特征在于，所述方差分析前提为：

各组样本的数据是否具有独立性；

同一判断区间内的各组的所述甲基化数据是否符合正态分布；

同一判断区间内的各组的所述甲基化数据是否符合方差齐性。

7.一种多样本间甲基化差异的检测装置，其特征在于，所述装置包括：

甲基化数据获取单元，用于获得多组样本的基因组测序片段对应于参考基因组序列上的位置信息以及各个样本的甲基化数据；

甲基化差异区域获取单元，用于对在所述参考基因组序列上滑动的同一判断区间内的各组的所述甲基化数据进行分析，获取上述组基因组上存在甲基化差异的甲基化差异区域；

甲基化差异组获取单元，用于对同一甲基化差异区域内的各组的所述甲基化数据进行分析，获取在所述甲基化差异区域内存在甲基化差异的具体组。

8.根据权利要求7所述的多样本间甲基化差异的检测装置，其特征在于，其中，所述甲基化差异区域获取单元包括：

设置单元，用于预设窗口长度，根据位置信息，从参考基因组序列起始端开始，以所述窗口长度作为判断区间的起始长度在参考基因组序列上设定判断区间；

判断单元，用于判断同一判断区间内的各组的所述甲基化数据是否符合方差分析前提；

分析单元，若符合方差分析前提，所述分析单元用于将位置在同一判断区间内的各组的所述甲基化数据进行组间方差分析，若不符合方差分析前提，所述分析单元用于将位置在同一判断区间内的各组的所述甲基化数据进行组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在该判断区间内是否存在甲基化差异；

控制单元，若存在甲基化差异，所述控制单元控制所述设置单元延长该判断区间，控制所述判断单元判断位置在同一延长后的判断区间内的各组的所述甲基化数据是否符合方差分析前提，控制所述分析单元在符合方差分析前提时进行前述组间方差分析，在不符合方差分析前提时进行前述组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在当前判断区间内是否存在甲基化差异的步骤，控制所述设置单元、所述判断单元及所述分析单元重复本步骤直至判断得出上述组在当前判断区间内不存在甲基化差异，并输出该当前判断区间的信息作为甲基化差异区域；

若不存在甲基化差异，所述控制单元控制所述设置单元从上一判断区间的末端开始在参考基因组序列上以所述窗口长度作为所述判断区间的起始长度设定下一判断区间，控制所述判断单元判断位置在同一判断区间内的各组的所述甲基化数据是否符合方差分析前提，控制所述分析单元在符合方差分析前提时进行前述组间方差分析，控制所述分析单元在不符合方差分析前提时进行前述组间Kruskal-Wallis非参数检验，得到分析结果，并根据分析结果判断上述组在该判断区间内是否存在甲基化差异的步骤，若存在甲基化差异则执行前述延长该判断区间、判断是否符合方差分析前提以及分析并判断上述组在当前判断区间内是否存在甲基化差异的步骤，若不存在甲基化差异则执行本步骤。

9.根据权利要求8所述的多样本间甲基化差异的检测装置，其特征在于，其中，所述甲基化差异组获取单元包括：

比较单元，用于通过最小显著性差异法对位置在同一甲基化差异区域内的且符合方差分析前提的各组的所述甲基化数据进行组间两两比较，通过Kruskal-Wallis Dunn方法对各组在所述甲基化差异区域内的且不符合方差分析前提的所述甲基化数据进行组间两两比较，得到比较结果；

输出单元，用于根据所述比较结果判断并输出在该甲基化差异区域内存在甲基化差异的具体组。

10.根据权利要求8所述的多样本间甲基化差异的检测装置，其特征在于，所述甲基化差异区域获取单元还包括：

终止单元，用于在所述分析单元判断完上述组在当前判断区间是否存在甲基化差异之后，判断当前判断区间是否已达所述参考基因组序列末端，若是，则终止所述设置单元、判断单元及分析单元。

11.根据权利要求8所述的多样本间甲基化差异的检测装置，其特征在于，在所述分析单元判断得出上述组在当前判断区间内存在甲基化差异时，所述设置单元延长该判断区间具体为：将该判断区间延长一个预设步长。

12.根据权利要求8所述的多样本间甲基化差异的检测装置，其特征在于，所述判断单元判断同一判断区间内的各组的所述甲基化数据是否符合方差分析前提具体为判断：

各组样本的数据是否具有独立性；

同一判断区间内的各组的所述甲基化数据是否符合正态分布；

同一判断区间内的各组的所述甲基化数据是否符合方差齐性。

Images