一种植物基因组差异甲基化区域的检测方法

摘要

本发明提供了植物基因组差异甲基化区域的检测方法，属于基因组学和生物信息学技术领域。通过植物全基因组甲基化测序分别获得处理组和对照组全基因组上单位点甲基化数据；分别读取处理组和对照组全基因组上单位点甲基化信息，分区段储存；采用以200bp为一个窗口，以50bp为步长滑动窗口扫描存储的全基因组每一区段的甲基化信息，将扫描得到的甲基化信息分区段储存到python字典中；对存储在python字典中的全基因组每一区段的甲基化信息进行筛选，根据得到的有效全基因组上单位点甲基化信息统计每个窗口中各区段的甲基化信息。本发明利用python字典把整个染色体以50bp为区段储存，方便后续的运算，节省了大量时间。

权利要求

1.一种植物基因组差异甲基化区域的检测方法，其特征在于，包括以下步骤：

1)通过植物全基因组甲基化测序获得处理组全基因组上单位点甲基化数据和对照组全基因组上单位点甲基化数据；

2)分别读取处理组全基因组上单位点甲基化信息和对照组全基因组上单位点甲基化信息，将读取的两组全基因组上单位点甲基化信息分区段储存，得到全基因组每一区段的甲基化信息；

3)以200bp为一个窗口，以50bp为步长滑动窗口扫描所述全基因组每一区段的甲基化信息，将扫描得到的全基因组每一区段的甲基化信息分区段储存到python字典中；

所述两组全基因组上差异甲基化区域的甲基化信息包括染色体，区域，甲基化类型，每一种甲基化类型对应的胞嘧啶位点数量，每一种甲基化类型对应的胞嘧啶位点数量所支持的甲基化reads数量、每一种甲基化类型对应的胞嘧啶位点数量所不支持的甲基化reads数量和每一种甲基化类型对应的胞嘧啶位点数量所支持与不支持甲基化reads数量的总和；

其中区域号的设置方法为以步长50bp为一个区域，从0开始编号，每移动一个区域顺延一个编号；所述甲基化类型包括CG类型、CHG类型或CHH类型；所述全基因组上差异甲基化区域的甲基化信息的储存形式为3层字典储存；最外侧键的名称为染色体；从外向内数第二层键的名称为区域；从外向内数第三层键的名称为甲基化类型；最内层键为所述区域内每一种甲基化类型对应的胞嘧啶位点数量所支持的甲基化reads数量、所述区域内每一种甲基化类型对应的胞嘧啶位点数量所不支持的甲基化reads数量和所述区域内每一种甲基化类型胞嘧啶位点数量的总和；

4)对步骤3)中存储在python字典中的全基因组每一区段的甲基化信息进行预筛选、Fisher检验、FDR校正和再筛选，得到有效全基因组上单位点甲基化信息；

5)根据所述有效全基因组上单位点甲基化信息，统计每个窗口中各区段的甲基化信息。

2.根据权利要求1所述检测方法，其特征在于，步骤4)中所述预筛选的条件为每个窗口中不低于6个胞嘧啶位点。

3.根据权利要求1所述检测方法，其特征在于，步骤4)中Fisher检验的条件为PValue小于0.05。

4.根据权利要求1所述检测方法，其特征在于，步骤4)中FDR 校正的条件为Q-Value值小于0.05。

5.根据权利要求1所述检测方法，其特征在于，步骤1)中根据所述处理组全基因组上单点甲基化数据，统计全基因组的CG类型、CHG类型或CHH类型的甲基化水平和总甲基化水平。

6.根据权利要求5所述检测方法，其特征在于，步骤4)再筛选的条件包括两种情况：对对照组和处理组均发生甲基化的区域而言，与所述总甲基化水平相比，存在40％的甲基化水平的显著差异和处理组或对照组至少一个的甲基化水平在10％以上；

对对照组无而处理组有或对照组有而处理组无的差异甲基化区域而言，与所述全基因组的CG类型、CHG类型和CHH类型的甲基化水平相比，CG类型和CHG类型的甲基化水平差异要达到30％，CHH类型的甲基化水平差异要达到20％。

说明书

一种植物基因组差异甲基化区域的检测方法

技术领域

[0001] 本发明属于基因组学和生物信息学技术领域，具体涉及一种植物基因组差异甲基化区域的检测方法。

背景技术

[0002] 进入21世纪以来，受极端天气的影响，我国由干旱导致的经济损失越来越大。据统计，2015年我国因干旱造成的直接经济损失高达579亿人民币。在如此严重的干旱情况下，研究植物的抗旱机制就显得尤为的重要和迫切。1942年，Waddington提出“epigenetics”一词，它在随后被Wolff等人定义为：DNA序列没有变化但是基因表达发生了可遗传变化的一门学科。在随后的研究中，发现了DNA甲基化，组蛋白甲基化，乙酰化等表观修饰现象。DNA甲基化是指在DNA甲基化酶的催化下，将S-腺苷-甲硫氨酸(SAM)上的甲基连接到胞嘧啶(C)的第五位碳原子上。DNA甲基化作为表观遗传学中的重要一环，在现在的研究中越来越被科学家们重视。目前的研究发现，不管是在动物还是植物中，DNA甲基化大部分都发生在胞嘧啶C中。根据胞嘧啶C后面的碱基，可以把DNA的甲基化分为三种类型，即：CG，CHG，CHH(H＝A，C，T)。例如：CGT属于CG类型，CTG碱基属于CHG类型，CTT属于CHH类型。在过去的十年间，DNA甲基化检测技术在不停的进步，发展出包括甲基化敏感多态性扩增技术(MSAP)，高分辨率熔解曲线法(HighResolutionMelting，HRM)和DNA甲基化测序的“金标准”--Bisulfitesequencing(Bs-Seq)在内的DNA甲基化检测技术。

[0003] DNA甲基化参与了植物生长发育中的各个过程，包括种子的萌发，植物的开花，果实的成熟，DNA的复制，转录，翻译以及翻译后修饰等。而DNA甲基化最让人关注的地方在于它和基因表达的关系。DNA甲基化在植物的抗逆过程中，也发挥着重要的作用。早在1990年科学家就发现，DNA甲基化参与了植物对外界环境的响应。2010年，Tan等人利用甲基化敏感多态性(MSAP)技术发现，当玉米在经受盐胁迫后，一些逆转录转座子，热激蛋白和蛋白激酶的基因甲基化水平发生变化，进一步发现，玉米蛋白磷酸酶2C(zmPP2C)的第一个内含子的甲基化水平的升高会使得zmPP2C表达下降。2013年，Zhong等人用重亚硫酸盐测序技术(BS-Seq)精准地得到了番茄的全基因组甲基化水平—甲基化组。2015年，Garg等人利用BS-Seq发现，在干旱条件下水稻的甲基化水平上升。在水稻染色体水平上，甲基化在着丝粒的附近呈现一个高水平状态，这很可能是因为在着丝粒附近有着大量的重复序列，而DNA甲基化在重复序列能维持在高水平状态。随后发现在干旱条件下，水稻基因上游的甲基化水平会发生明显的变化。值得注意的是，在离转录起始位点300～400bp处，DNA甲基化率出现一个峰值，从峰值处开始到转录起始位点，DNA甲基化水平急剧下降。差异甲基化区域(DMR)指的是在染色体的相同位置，DNA甲基化水平有显著差异的区域。分析这些区域可以很好了解在逆境胁迫下体内发生甲基化变化的通路以及相关基因。Garg发现，在干旱处理下，有15525个DMR甲基化水平上升，而有7100个DMR的甲基化水平下降；在盐处理下，有7584个DMR的甲基化率上升，而有7430个DMR的甲基化率下降。把DMR相关的基因进行GO和KEGG分析后发现，这些基因主要涉及的是植物的代谢过程，对刺激物的反应以及对非生物逆境的响应过程。进一步分析这些基因的表达和DNA甲基化的关系可以发现，染色体上有的区域DNA甲基化和基因表达呈负相关关系，有的区域DNA甲基化和基因表达呈正相关的关系。然而DNA甲基化的研究仅仅限于少数几个物种，与植物响应干旱的关系和机制也不清晰，急需更多的研究和发现。

[0004] 目前，申请公开号为CN102061337A，一种组织特异性差异甲基化区域tDMR检测方法的专利公开了一种差异甲基化区域的检测方法，虽然该方法对多种基因组来源具有较好的适用性，但是检测分析过程耗时较长，不利于实现快速分析检测。

发明内容

[0005] 有鉴于此，本发明的目的在于提供一种植物基因组差异甲基化区域的检测方法的方法，所述方法不仅能实现统计差异甲基化区域信息，而且提高分析速度，提高工作效率。

[0006] 本发明提供的一种植物基因组差异甲基化区域的检测方法，包括以下步骤：

[0007] 1)通过植物全基因组甲基化测序获得处理组全基因组上单位点甲基化数据和对照组全基因组上单位点甲基化数据；

[0008] 2)分别读取处理组全基因组上单位点甲基化信息和对照组全基因组上单位点甲基化信息，将读取的两组全基因组上单位点甲基化信息分区段储存，得到全基因组每一区段的甲基化信息；

[0009] 3)采用以200bp为一个窗口，以50bp为步长滑动窗口扫描所述全基因组每一区段的甲基化信息，将扫描得到的全基因组每一区段的甲基化信息分区段储存到python字典中；

[0010] 4)对步骤3)中存储在python字典中的全基因组每一区段的甲基化信息进行预筛选、Fisher检验、FDR校正和再筛选，得到有效全基因组上单位点甲基化信息；

[0011] 5)根据所述有效全基因组上单位点甲基化信息，统计每个窗口中各区段的甲基化信息。

[0012] 优选的，步骤4)中所述预筛选条的条件为每个窗口中不低于6个胞嘧啶位点。

[0013] 优选的，步骤4)中Fisher检验的条件为PValue小于0.05。

[0014] 优选的，步骤4)中FDR校正的条件为Q-Value值小于0.05。

[0015] 优选的，步骤1)中根据所述处理组全基因组上单点甲基化数据，统计全基因组的CG类型、CHG类型或CHH类型的甲基化水平和总甲基化水平。

[0016] 优选的，步骤4)再筛选的条件包括两种情况：

[0017] 对对照组和处理组均发生甲基化的区域而言，与所述总甲基化水平相比，存在40％的甲基化水平的显著差异和处理组或对照组至少一个的甲基化水平在10％以上；

[0018] 对对照组无而处理组有或对照组有而处理组无的差异甲基化区域而言，与全基因组的CG类型、CHG类型和CHH类型的甲基化水平相比，CG类型和CHG类型的甲基化水平差异要达到30％，CHH类型的甲基化水平差异要达到20％。

[0019] 优选的，步骤3)中全基因组上差异甲基化区域的甲基化信息包括染色体，区域，甲基化类型，每一种甲基化类型对应的胞嘧啶位点数量，每一种甲基化类型对应的胞嘧啶位点数量所支持的甲基化reads数量、每一种甲基化类型对应的胞嘧啶位点数量所不支持的甲基化reads数量和所述区域内每一种甲基化类型胞嘧啶位点数量的总和；

[0020] 其中区域号的设置方法为以步长50bp为一个区域，从0开始编号，每移动一个区域顺延一个编号；

[0021] 所述甲基化类型包括CG类型、CHG类型或CHH类型。

[0022] 优选的，步骤3)中所述全基因组上差异甲基化区域的甲基化信息的储存形式为3层字典储存；

[0023] 最外侧键的名称为染色体；

[0024] 从外向内数第二层键的名称为区域；

[0025] 从外向内数第三层键的名称为甲基化类型；

[0026] 最内层键为所述区域内每一种甲基化类型对应的胞嘧啶位点数量所支持的甲基化reads数量、所述区域内每一种甲基化类型对应的胞嘧啶位点数量所不支持的甲基化reads数量和所述区域内每一种甲基化类型胞嘧啶位点数量的总和。

[0027] 本发明提供的植物基因组差异甲基化区域的检测方法，通过植物全基因组甲基化测序分别获得处理组和对照组全基因组上单位点甲基化数据；分别读取处理组和对照组的全基因组上单位点甲基化信息，将读取的两组全基因组上单位点甲基化信息分区段储存；采用以200bp为一个窗口，以50bp为步长滑动窗口扫描储存到python字典中全基因组每一区段的甲基化信息，将扫描得到的甲基化信息分区段储存到python字典中；对存储在python字典中的全基因组每一区段的甲基化信息进行预筛选、Fisher检验、FDR校正和再筛选，得到有效全基因组上单位点甲基化信息；统计有效全基因组上单位点甲基化信息每个窗口中各区段的甲基化信息。利用以200bp为一个窗口，以50bp为步长滑动窗口的方法扫描，并将扫描的甲基化差异区域信息储存到python字典，可显著提高筛选不同抗性处理间基因组水平上DMR的效率，比对照方法(华大专利)快8倍以上。实验证明：计算15G全基因组甲基化数据量时，是利用python的字典存储方式，每50bp储存一次，方便了后续的运算。计算0～200bp这个区域的时候，只需要输出第0/1/2/3这四次的结果就行了。本发明提供的方法耗时极大缩短，所获得DMR数量上相差不大，显然在效率上很有优势，且数据量越大，优势越明显，甚至可节省几个月的时间。

附图说明

[0028] 图1为滑窗法扫描全基因组上单位点甲基化数据示意图；

[0029] 图2为小立碗藓三种类型的DNA甲基化水平在基因组上的比例和基因结构上的分布情况；其中图2-A为在M01样品中小立碗藓三种类型的DNA甲基化比例；图2-B为在M02样品中小立碗藓三种类型的DNA甲基化比例；图2-C为在M01样品中小立碗藓三种类型的DNA甲基化分布；图2-D为在M02样品中小立碗藓三种类型的DNA甲基化分布；

[0030] 图3为干旱处理下三种类型的DMR统计结果，Hyper表示DNA甲基化上升，Hypo表示DNA甲基化下降；Total表示下降和上升的两类总和；

[0031] 图4为干旱处理下受到甲基化调控的基因参与的主要生物学过程；

[0032] 图5为本发明提供的方法与已发表方法的比较耗时结果柱形图。

具体实施方式

[0033] 本发明提供的一种植物基因组差异甲基化区域的检测方法，包括以下步骤：

[0034] 1)通过植物全基因组甲基化测序获得处理组全基因组上单位点甲基化数据和对照组全基因组上单位点甲基化数据；

[0035] 2)分别读取处理组全基因组上单位点甲基化信息和对照组全基因组上单位点甲基化信息，将读取的两组全基因组上单位点甲基化信息分区段储存到python字典中，得到全基因组每一区段的甲基化信息；

[0036] 3)采用以200bp为一个窗口，以50bp为步长滑动窗口扫描所述全基因组每一区段的甲基化信息，将扫描得到的全基因组每一区段的甲基化信息分区段储存到python字典中；

[0037] 4)对步骤3)中存储在python字典中的全基因组每一区段的甲基化信息进行预筛选、Fisher检验、FDR校正和再筛选，得到有效全基因组上单位点甲基化信息；

[0038] 5)根据所述有效全基因组上单位点甲基化信息，统计每个窗口中各区段的甲基化信息。

[0039] 本发明通过植物全基因组甲基化测序获得处理组全基因组上单位点甲基化数据和对照组全基因组上单位点甲基化数据。

[0040] 在本发明中，处理组植物的处理方式包括非正常生长条件的胁迫处理，例如高温、低温、干旱、洪涝、高盐、低盐、高光照或低光照等处理方法。在本发明实施例中，为了举例说明本发明所要保护的分析方法的实施步骤，采用干旱处理方法为例加以说明。所述对照组植物为正常条件下生长得到的植物。

[0041] 本发明对植物的种类不做任何具体限定，本发明提供的分析方法适合所有植物的基因组甲基化DMR分析。在本发明实施例中，为了举例说明本发明所要保护的分析方法的实施步骤，是以小立碗藓(Physcomitrella patens)为例加以说明。

[0042] 本发明对所述植物全基因组甲基化测序的方法没有特殊限制，采用本领域所熟知的植物全基因组甲基化测序的方法即可。在本发明实施例中，所述植物全基因组甲基化测序的方法优选为BS-Seq方法。所述植物全基因组甲基化测序前优选提取得到植物全基因组DNA。本发明对植物全基因组DNA的提取方法没有特殊限制，采用本领域所熟知的基因组DNA的抽取方法即可，例如手提法或试剂盒法。

[0043] 本发明对所述植物全基因组甲基化测序的结果有要求，包括以下内容：每个样本得到22G高质量reads数，总质量数达到88.7Gbp，转化率达到99.5％以上；覆盖基因组81％～91％的胞嘧啶位点，基因组测序深度达到38×～68×。

[0044] 单位点DNA甲基化水平的获取返回的高质量数据(去掉接头和低质量的数据cleanreads)与小立碗藓的基因组(https://phytozome.jgi.doe.gov/pz/portal.html#！info？alias＝Org_Ppatens)进行比对，得到全基因组单个胞嘧啶位点的甲基化信息(保证每个C位点至少6条reads能比对上)。用于统计单个位点甲基化水平的软件是bsmap-2.74(允许3个gap)，得到的bam(或者sam)文件用于后续统计分析得到单个胞嘧啶位点的甲基化状态。

[0045] 得到处理组全基因组上单位点甲基化信息和对照组全基因组上单位点甲基化信息后，本发明分别读取处理组全基因组上单位点甲基化信息和对照组全基因组上单位点甲基化信息，将读取的两组全基因组上单位点甲基化信息分区段储存，得到全基因组每一区段的甲基化信息。

[0046] 在本发明中，根据所述处理组全基因组上单点甲基化数据，统计全基因组的CG类型、CHG类型或CHH类型的甲基化水平和总甲基化水平，以便用于后续筛选中作为比较的基础。CG类型的甲基化率为CG类型的甲基化数量/CG总的数量；CHG为CHG类型的甲基化数量/CHG总的数量；CHH为CHH类型的甲基化数量/CHH总的数量。在统计全基因组三种类型甲基化水平时，甲基化率小于0.1的胞嘧啶位点认定为未甲基化的位点。

[0047] 在本发明中，所述分区段储存到python字典所用的模块为python模块。所述python模块包括argparse模块、os模块和scipy模块。所述argparse模块是用来输入命令行参数。os是python内置模块之一，用来python和linux服务器的交互。scipy是python模块中专门用来统计运算的一个模块，在本发明中用来进行fisher检验。所述区段的长度为50bp。

[0048] 得到全基因组每一区段的甲基化信息后，本发明采用以200bp为一个窗口，以50bp为步长滑动窗口扫描所述全基因组每一区段的甲基化信息，将扫描得到的全基因组每一区段的甲基化信息分区段储存到python字典中。

[0049] 在本发明中，所述两组全基因组上差异甲基化区域的甲基化信息优选包括染色体，区域，甲基化类型，每一种甲基化类型对应的胞嘧啶位点数量，每一种甲基化类型对应的胞嘧啶位点数量所支持的甲基化reads数量、每一种甲基化类型对应的胞嘧啶位点数量所不支持的甲基化reads数量和每一种甲基化类型对应的胞嘧啶位点数量所支持与不支持甲基化reads数量的总和。其中区域号的设置方法为以步长50bp为一个区域，从0开始编号，每移动一个区域顺延一个编号。所述甲基化类型包括CG类型、CHG类型或CHH类型。所述全基因组上差异甲基化区域的甲基化信息的储存形式为3层字典储存；最外侧键的名称为染色体；从外向内数第二层键的名称为区域；从外向内数第三层键的名称为甲基化类型；最内层键为所述区域内每一种甲基化类型对应的胞嘧啶位点数量所支持的甲基化reads数量、所述区域内每一种甲基化类型对应的胞嘧啶位点数量所不支持的甲基化reads数量和所述区域内每一种甲基化类型胞嘧啶位点数量的总和。所述字典形式为{chr:{binOrder:{type:[support,unsupport,total]}}}，例{chr1:{2:{CG:1580,265,25}}}代表的含义为：在一号染色体中，第二个区域(50bp为一个区域，从1开始)中，CG类型的胞嘧啶位点总共有25个，其中这25个位点支持甲基化的reads总和为1580，不支持甲基化的reads总数为265。

[0050] 存储后，本发明对存储在python字典中的全基因组每一区段的甲基化信息进行预筛选、Fisher检验、FDR校正和再筛选，得到有效全基因组上单位点甲基化信息。

[0051] 在本发明中，所述预筛选、Fisher检验、FDR校正和再筛选的过程也是以窗口为单位进行的。所述窗口的长度包括200bp。

[0052] 在本发明中，所述预筛选条的条件优选为每个窗口中不低于6个胞嘧啶位点。Fisher检验的条件优选为PValue小于0.05。FDR校正的条件优选为Q-Value值小于0.05。所述再筛选的条件优选包括两种情况：a.对对照组和处理组均发生甲基化的区域而言，与所述总甲基化水平相比，存在40％的甲基化水平的显著差异和处理组或对照组至少一个的甲基化水平在10％以上；b.对对照组无而处理组有或对照组有而处理组无的差异甲基化区域而言，与全基因组的CG类型、CHG类型和CHH类型的甲基化水平相比，CG类型和CHG类型的甲基化水平差异要达到30％，CHH类型的甲基化水平差异要达到20％。

[0053] 本发明根据所述有效全基因组上单位点甲基化信息，统计每个窗口中各区段的甲基化信息。

[0054] 在本发明中，所述从染色体的起始位置开始(1)，1～200为第一个窗口，50～250为第二个窗口，100～300为第三个窗口，以此类推。其中第一个窗口对应上面字典中binorder中的0～3，第二窗口对应上面字典中binorder中的1-4，以此类推。例如：在统计第一个窗口中的甲基化信息时，只需要统计上述字典中binorder中键为0～3的所有信息就行，这样就能快速统计出每个窗口的甲基化信息。

[0055] 得到每个窗口中各区段的甲基化信息，本发明统计DMR的分布和DMR相关的基因，对相关基因进行功能注释，得到干旱诱导了小立碗藓大量次级代谢、蛋白转运、细胞壁代谢和信号传导相关的基因发生DNA甲基化变化。本发明对所述功能注释的方法没有特殊限制，采用本领域所熟知的功能注释方案即可。

[0056] 下面结合实施例对本发明提供的一种植物基因组差异甲基化区域的检测方法的方法进行详细的说明，但是不能把它们理解为对本发明保护范围的限定。

[0057] 实施例1

[0058] 1、小立碗藓的培养以及干旱处理

[0059] 使用BCDAT培养基培养小立碗藓，取原丝体材料打磨继代，在14小时光照(光照强度500μmol·m-2·s-1)/10小时黑暗，25℃培养室中培养40天。将植物材料取出并进行干燥复水处理。约干燥处理2小时，失水量达到80％左右时，取样。BCDAT培养基：0.1％的TES，1％的储存液B，1％储存液D，1％储存液C，5mmol/L的酒石酸铵，1mmol/L的CaCl2·2H2O，0.8％的琼脂粉(固体)。高压灭菌后，室温保存。

[0060] 2、基因组DNA抽提

[0061] 选取了对照组即处理前(M01)和干旱处理(M02)两个材料作为样品，使用CTAB法提取基因组DNA。提取得到基因组DNA质量通过Nanodrop-2000及琼脂糖凝胶电泳测试。CTAB提取液：2％十六烷基三乙基溴化铵(cetytriethylammonium bromide，CTAB)，20mmol/L EDTA(pH值8.0)，100mmol/L Tris-HCl(pH值为8.0)，0.2％巯基乙醇，1.4mol/L NaCl，1％PVP。

[0062] 3、小立碗藓全基因组甲基化测序

[0063] 基因组DNA样品送北京百迈客生物科技有限公司用BS-Seq方法进行DNA甲基化水平检测。每个样本得到22G高质量reads数，总质量数达到88.7Gbp，转化率达到99.5％以上；覆盖基因组81％-91％的胞嘧啶位点，基因组测序深度达到38×-68×。

[0064] 为了探究在基因组层面2个样本的DNA甲基化水平是否发生了明显变化，统计了DNA甲基化在全基因组的分布。结果显示，在M01中，CHH类型的甲基化最多，占81.54％；其次是CG，占9.33％；最后是CG，占9.13％。单个位点DNA甲基化水平的获取返回的高质量数据(去掉接头和低质量的数据clean  reads)与小立碗藓的基因组(https://phytozome.jgi.doe.gov/pz/portal.html#！info？alias＝Org_Ppatens)进行比对，得到全基因组单个胞嘧啶位点的甲基化信息(保证每个C位点至少6条reads能比对上)。用于统计单个位点甲基化水平的软件是bsmap-2.74(允许3个gap)，得到的bam(或者sam)文件采用python模块(methratio.py)统计分析得到单个胞嘧啶位点的甲基化状态。

[0065] 4、全基因组水平上的甲基化水平变化分析

[0066] 根据得到的胞嘧啶单个位点的甲基化信息，统计出全基因组的甲基化水平。

[0067] CG类型的甲基化率为CG类型的甲基化数量/CG总的数量，CHG和CHH的计算方法与CG类型的甲基化率的计算方法一样。在统计全基因组甲基化水平时，甲基化率小于0.1的胞嘧啶位点认定为未甲基化的位点。

[0068] 结果：小立碗藓三种类型的DNA甲基化水平统计结果(图2)显示：在小立碗藓中，CHH类型的DNA甲基化占绝大部分(大于80％)而且甲基化水平比较均匀，而CG类型、CHG类型的DNA甲基化比较少，但是都处于高水平的DNA甲基化水平。在基因body区域的甲基化水平普遍大大低于基因外区域。这个结果和水稻，拟南芥等其他高等植物中报道的不一样，在其他植物中，占比最多的是CG类型的DNA甲基化，CHH是最少的，而且在其他植物中的CHH类型的DNA甲基化大部分处于低甲基化的水平，另外其他植物基因body区域有比较高的甲基化水平。

[0069] 实施例2

[0070] 把实施例1中得到的对照组(M01)和处理组(M02)全基因组上所有的染色体上的单个胞嘧啶位点的甲基化信息储存在python字典中。为了后面能快速找到染色体上某个区域的甲基化信息，存储时分区段(每个区域50bp)去储存胞嘧啶位点的信息。字典的具体形式为{chr:{binOrder:{type:[support,unsupport,total]}}}，比如{chr1:{2:{CG:1580,265,25}}},含义为：在一号染色体中，第二个区域(50bp为一个区域，从1开始)中，CG类型的胞嘧啶位点总共有25个，其中这25个位点支持甲基化的reads总和为1580，不支持甲基化的reads总数为265。

[0071] 比对对照组和处理组全基因组上所有的染色体上的单个胞嘧啶位点的甲基化信息，找出在基因组的相同位置发生DNA甲基化变化的区域，得到DMR区域信息文件。

[0072] 结果见图3，图3为干旱处理下三种类型的DMR统计结果，Hyper表示DNA甲基化上升，Hypo表示DNA甲基化下降；Total表示下降和上升的两类总和。以上结果显示，当受到干旱处理时，CHH，CHG，CG类型的DMR分别有15535，4578和2266个。CHH类型的DMR数目要远远大于CG，CHG类型的DMR数目。

[0073] 采用滑窗法来统计DMR区域的甲基化信息。所述滑窗法原理(见图1)，具体方法如下：以200bp为一个窗口，以50bp为步长的滑动窗口来扫描处理组和对照组的全基因组上的胞嘧啶位点甲基化信息(包括染色体，区域，甲基化类型，每一种甲基化类型对应的胞嘧啶位点数量，每一种甲基化类型对应的胞嘧啶位点数量所支持的甲基化reads数量、每一种甲基化类型对应的胞嘧啶位点数量所不支持的甲基化reads数量和每一种甲基化类型对应的胞嘧啶位点数量所支持与不支持甲基化reads数量的总和。其中区域号的设置方法为以步长50bp为一个区域，从0开始编号，每移动一个区域顺延一个编号)，统计每个区段(50bp)的reads数，对统计的胞嘧啶位点的信息以窗口为单位依次进行预筛选、Fisher检验、FDR校正和再筛选，得到有效全基因组DMR区域的信息。预筛选的条件为每个窗口保证有6个胞嘧啶位点。Fisher检验的条件为PValue小于0.05。用FDR(False discovery rate)校正P-Value后，Q-Value值要小于0.05。再筛选的条件包括40％的甲基化水平的显著差异以及至少一个样本的甲基化水平在10％以上；对对照组无而处理组有或对照组有而处理组无的差异甲基化区域而言，与全基因组的CG类型、CHG类型和CHH类型的甲基化水平相比，CG类型和CHG类型的甲基化水平差异要达到30％，CHH类型的甲基化水平差异要达到20％。

[0074] 将有效全基因组DMR区域的信息进行扫描，从染色体的起始位置开始，1～200为第一个窗口，50～250为第二个窗口，100～300为第三个窗口，以此类推。其中第一个窗口对应上面字典中binorder中的0-3，第二窗口对应上面字典中binorder中的1～4，以此类推。例如：需要统计第一个窗口中的甲基化信息，只需要统计上述字典中binorder中键为0～3的所有信息就行，快速统计出每个窗口的甲基化信息，将各窗口中特定甲基化类型的数量相加，能够快速得到全基因组DMR的数量。

[0075] 计算0～200bp这个区域的时候，只需要输出第0/1/2/3这四次的结果就行了，主要利用python字典计算把整个染色体以50bp为一小段储存，方便了后续的运算，节省了大量时间。计算15G全基因组甲基化数据量时，与已有方法(公开的专利CN201010557131)DMR分析速度比较，本发明耗时极大缩短，所获得DMR数量上相差不大，显然在效率上很有优势(图5)。当数据量大时越发突出优势，将几个月的时间节省了。

[0076] 实施例4

[0077] 本发明利用在线数据库进行生物学过程分析，具体见网站http://wego.genomics.org.cn/和http://kobas.cbi.pku.edu.cn/。经过分析得到干旱处理下受到甲基化调控的基因参与的主要生物学过程(见图4)。

[0078] 以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。