CN108509769B - 确定预定物种的基因表达和甲基化修饰调控的关系的方法 - Google Patents

Abstract

本发明公开了确定预定物种的基因表达和甲基化修饰调控的关系的方法，其包括：(1)以预定物种的父本、母本及其子代样本作为待测样本，并进行重亚硫酸盐全基因组甲基化测序；(2)确定子代样本中所有胞嘧啶位点以及包含胞嘧啶位点的读段；(3)确定所述包含胞嘧啶位点的读段中属于等位基因序列的读段；(4)确定每对目的读段的表观基因型；(5)将目的片段进行基因归类；(6)统计各基因包含的目的片段数目、三种表观基因型的比率，并获得各基因的表达量信息；(7)划分出多种候选基因组合；(8)进行皮尔逊相关性分析；(9)筛选目标基因组合；(10)进行多元线性回归。由此能够有效确定预定物种的基因表达和甲基化修饰调控的关系。

Description

确定预定物种的基因表达和甲基化修饰调控的关系的方法

技术领域

本发明涉及分子遗传技术领域，具体涉及表观遗传调控技术领域，更具体地，涉及确定预定物种的基因表达和甲基化修饰调控的关系的方法。

背景技术

基因表达受到多维度因素调节，十分复杂。而将单调控因素从多维因素调控网络剥离，构建其影响基因表达的数学模型，进而确定该单调控因素对于基因表达调控的影响，对于表观遗传调控研究意义重大。其中，DNA甲基化是非常重要的表观遗传调控因素，因而，构建某物种的DNA甲基化影响基因表达的数学模型，并确定该物种基因表达和甲基化修饰调控的关系，至关重要。

然而，目前，尚无确定预定物种基因表达和甲基化修饰调控的关系的报道。

发明内容

本发明旨在至少解决现有技术中存在的技术问题之一。为此，本发明的一个目的在于提出一种能够有效构建预定物种的DNA甲基化影响基因表达的数学模型，并确定该物种基因表达和甲基化修饰调控的关系的手段。

需要说明的是，本发明是基于发明人的下列发现和工作而完成的：

已知基因表达受到多维度因素调节，十分复杂，将单调控因素从多维因素调控网络剥离难度极大，而DNA甲基化又具有高度可变性，因而，将DNA甲基化从多维因素调控网络剥离，并构建基因表达与甲基化修饰调控关系的数学模型，难度非常大。

因而，发明人进行了一系列的科学研究和实验探索。最终，发明人在多种实验设计中，惊喜地发现了一种最为科学有效的确定预定物种基因表达和甲基化修饰调控的关系的方法：利用发明人首创的胞嘧啶位点表观基因型分型方法获得子代样本的胞嘧啶片段分型数据；基于胞嘧啶片段分型数据的数量进行多种分组，并针对各组，对各表观基因型与基因表达量做皮尔逊相关性分析，并筛选过滤各组的数据，找到在同一数据点时三种表观基因型同时与基因表达量极显著的分组，选择P值最小的一组数据；数据经过过滤筛选之后，将最显著的一组数据进行线性回归数学建模，以便得到可计算的基因表达受甲基化修饰调控关系的数学模型，从而可有效确定预定物种基因表达和甲基化修饰调控的关系。

进而，在本发明的一个方面，本发明提供了一种确定预定物种的基因表达和甲基化修饰调控的关系的方法。根据本发明的实施例，该方法包括以下步骤：

(1)以所述预定物种的父本、母本及其子代样本作为待测样本，并对所述待测样本进行重亚硫酸盐全基因组甲基化测序，以便分别获得父本、母本及其子代样本的测序数据，其中各测序数据均由多个读段组成；

(2)将所述子代样本的测序数据与参考基因组进行比对，以便确定所述子代样本中的所有胞嘧啶位点以及包含胞嘧啶位点的读段；

(3)针对所述子代样本的包含胞嘧啶位点的读段，对各胞嘧啶位点均进行位点上下游Call SNPs检测，以便确定所述子代样本的包含胞嘧啶位点的读段中属于等位基因序列的读段，并将每属于等位基因序列的两个读段作为一对目的读段；

(4)将所述子代样本的各目的读段分别与父本和母本的测序数据进行比对，以便确定每一对目的读段的两个读段的亲本来源，并基于每对目的读段的胞嘧啶甲基化信息，确定每对目的读段的表观基因型，

其中，所述目的读段的表观基因型分为三种：纯合甲基化基因型、杂合甲基化基因型和纯合非甲基化基因型，其中目的读段中分别来源于父本和母本的两个读段在同一位置上均发生胞嘧啶甲基化的，该目的片段为纯合甲基化基因型；两个读段在同一位置上均不发生胞嘧啶甲基化的，该目的片段为纯合非甲基化基因型；两个读段在同一位置上只有任意一个发生胞嘧啶甲基化的，该目的片段为杂合甲基化基因型；

(5)基于参考基因组序列，按照序列来源，将多对目的片段分别进行基因归类；

(6)统计各基因包含的目的片段的数目，确定每个基因包含的目的片段中三种表观基因型的比率，并获得各基因的表达量信息；

(7)基于设定候选基因包含的目的片段的数目的不同，对各基因进行多次分组，以便划分出多种候选基因组合；

(8)针对所述多种候选基因组合的每一种，将候选基因组合的所有基因的三种表观基因型比率分别与对应基因的基因表达量进行皮尔逊相关性分析；

(9)基于多种候选基因组合的皮尔逊相关性分析结果，筛选三种表观基因型比率同时与基因表达量显著相关且三个P值最小的一种候选基因组合作为目标基因组合；

(10)基于目标基因组合的表观基因型和基因表达量数据，进行多元线性回归，以便确定所述预定物种的基因表达和甲基化修饰调控的关系。

发明人惊奇地发现，利用该方法，通过皮尔逊相关性分析与多元线性回归，能够有效建立基因表达受甲基化修饰调控的关系的数学模型，从而能够有效确定所述预定物种的基因表达和甲基化修饰调控的关系。由此，也为从事DNA甲基化相关研究的科学工作者提供了一种有效的工具，通过该方法，可以实现对不同物种、不同样本进行基因表达受甲基化修饰调控关系的建模，确定各物种(待测样本)的基因表达和甲基化修饰调控的关系，从而能够有效促进表观遗传调控领域的技术研究。进一步，利用构建好的数学模型，可以有效了解该物种关键基因的甲基化修饰状况，进而能够通过5氮杂胞苷等可调节生物体内基因甲基化水平的试剂或其他手段，调节对应物种生物体内该关键基因的甲基化水平，进而达到调节该关键基因的表达，并最终调控生物体特定生理生化过程的目的。

根据本发明的一些实施例，步骤(2)进一步包括：利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96将比对结果进行染色体坐标排序、reads去重复处理。由此，比对结果更可靠，确定的胞嘧啶位点信息准确性高，有利于后续步骤的进行。

根据本发明的一些具体示例，利用GATK2对胞嘧啶位点附近5-15bp进行Call SNPs检测。由此，Call SNPs检测结果准确可靠，有利于后续步骤的进行。

根据本发明的实施例，利用SPSS 20进行所述多元线性回归。由此，建模结果准确性高，可重复性好。

其中，为方便理解，下面对本发明的方法的一些步骤进行详细说明：

首先，需要说明的是，所述预定物种优选为有参考基因组的物种，所述参考基因组为所述预定物种的已测序基因组或所述预定物种的近缘物种已测序基因组。例如，所述预定物种可以为人类、地山雀、金丝猴、非洲水稻、藏猪或复活草等有参考基因组的物种。

在步骤(1)中，确定待测样本后，对待测样本进行重亚硫酸盐全基因组甲基化测序，获得父母本和子代样本的基因组序列。根据本发明的一些实施例，在进行测序之前还可以进一步包括对待测样本进行基因组DNA提取的步骤；所述基因组DNA提取的方法不受特别限制，例如可以采用CTAB法。

而在得到各待测样本的基因组DNA后，可以进一步对基因组DNA进行质量、纯度和浓度检测，筛选获得合格的DNA样品进行下一步操作，若基因组DNA样品质量不合格需要重新提取基因组DNA。其中，基因组DNA质量检测优选使用琼脂糖凝胶电泳进行，可以根据电泳条带分析基因组DNA降解程度以及是否有RNA污染，操作方便。在检测基因组DNA质量合格后，优选进一步检测其纯度是否合格。根据本发明的一些实施例，优选的采用Nanodrop检测基因组DNA的纯度，当OD260/280为1.6～1.8时，认为基因组DNA的纯度合格，若基因组DNA的OD260/280不在上述范围内，则所述基因组DNA纯度不合格，需要重新提取。根据本发明的实施例，在筛选获得质量和纯度合格的基因组DNA后，优选的进一步采用Qubit对基因组DNA浓度进行精确定量，为后续操作提供依据。

其中，重亚硫酸盐全基因组甲基化测序的方法也不受特别限制(任何当前已知的技术或是未来的新方法均可)，只要能够有效获得待测样本的全基因组甲基化信息即可。根据本发明的一些具体示例，在获得用合格的父母本和子代样本基因组DNA样品后，采用重亚硫酸盐法构建上述合格的待测样本的基因组DNA样品测序文库，进而进行重亚硫酸盐全基因组甲基化测序。根据本发明的一些实施例，构建待测样本基因组DNA样品测序文库时可以加入比例为建库DNA起始量的1/1000的阴性对照lambda DNA。根据本发明的另一些实施例，可以按照以下步骤构建基因组DNA样品测序文库：首先随机打断基因组DNA至200～300bp获得DNA片段，所述随机打断基因组DNA优选的使用Covaris S220进行；然后对DNA片段进行平末端修复，所述平末端修复优选的使用内切酶和外切酶将黏性末端修复成为平末端；随后在平末端DNA片段后加尾巴A碱基获得带尾巴A的DNA片段，所述加尾巴A碱基步骤采用本领域常规的手段即可，无其他特殊要求；在带尾巴A的DNA片段上连接测序接头后进行Bisulfite处理；所述的Bisulfite处理优选的采用Zymo Research的EZ DNA MethylationGold Kit进行，所述Bisulfite处理的作用是将未发生甲基化的C变成U，在后续的PCR扩增过程后变为T，而甲基化的C则保持不变。在Bisulfite处理后进行PCR扩增即可获得基因组DNA样品测序文库。

根据本发明的实施例，进行重亚硫酸盐全基因组甲基化测序时，在获得基因组DNA样品测序文库后，可以进一步包括对文库进行质检筛选，以便获得合格基因组DNA样品测序文库。根据本发明的一些实施例，在质检文库之前优选的使用Qubit2.0对文库DNA浓度进行初步定量，然后根据文库DNA的浓度将其进行稀释，文库稀释后的浓度优选的为1ng/ul；在文库稀释后，优选的使用Agilent2100对文库的插入片段长度进行检测，筛选后的文库插入片段在320～520bp为合格；使用Q-PCR方法对文库的有效浓度进行准确定量，所述文库有效浓度>2nM为合格。

根据本发明的实施例，在获得合格基因组DNA样品测序文库后，即可对合格的DNA样品文库进行双末端Hiseq测序，以便获得父本、母本和子代样本的基因组序列，也即测序数据。

根据本发明的一些实施例，在步骤(2)中，将所述子代样本的测序数据与参考基因组进行比对时，所述参考基因组即为所述预定物种的基因组，如该预定物种未测序，则优选的选取该预定物种的已测序的近缘物种基因组作为参考基因组。

如前所述，步骤(2)可以进一步包括：利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96将比对结果进行染色体坐标排序、reads去重复处理。由此，比对结果准确，获得待测胞嘧啶位点后，有利于后续步骤(3)对已知胞嘧啶位点进行位点上下游Call SNPs检测，从而能够有效地区分子代的等位基因序列(由于序列经由重亚硫酸盐甲基化测序，所以胞嘧啶上的甲基化修饰状况已被记录)。

如前所述，步骤(4)中的所述目的读段的表观基因型分为三种：纯合甲基化基因型、杂合甲基化基因型和纯合非甲基化基因型，其中目的读段中分别来源于父本和母本的两个读段在同一位置上均发生胞嘧啶甲基化的，该目的片段为纯合甲基化基因型；两个读段在同一位置上均不发生胞嘧啶甲基化的，该目的片段为纯合非甲基化基因型；两个读段在同一位置上只有任意一个发生胞嘧啶甲基化的，该目的片段为杂合甲基化基因型。由此，在获得被区分的子代等位基因序列后，将已获得的被区分过子代等位基因序列与父母本基因组序列比对，即可有效完成胞嘧啶表观基因型分型，以便进行后续的步骤。其中，需要说明的是，该表观基因型分型方法是发明人首创。

根据本发明的实施例，在步骤(5)中，将多对目的片段分别进行基因归类包括：将所述多对目的片段(也即胞嘧啶甲基化分型片段)通过比对参考基因组序列进行注释，将不同片段归类于所属基因条目下，完成目的片段(也即胞嘧啶表观基因型分型的片段)与基因的mapping，最终不同基因获得的mapping成功的目的片段的数目是不同的。

根据本发明的实施例，在步骤(6)中，统计各基因包含的目的片段的数目，确定每个基因包含的目的片段中三种表观基因型的比率，包括：统计同一基因下不同类型的表观基因型的数量以及它们的比例，例如：某基因下mapping成功的目的片段共57个目的片段(有时也称为“胞嘧啶表观基因型分型片段”)中纯合甲基化基因型(有时也称为“表观基因型”)为30个(则，其比例＝30/57＝52.6％)，杂合甲基化基因型为15个(比例＝15/57＝26.3％)，而纯合非甲基化基因型为12个(比例＝12/57＝21.1％)。

根据本发明的一些实施例，在步骤(7)中，基于设定候选基因包含的目的片段的数目的不同，对各基因进行多次分组，其分组依据的目的片段的数目间隔不受特别限制，可以依据实际数据情况调整，例如可以分为候选基因包含的目的片段的数目≥5，≥10，≥15，≥20，≥25，≥30，≥35，≥40，≥45，≥50，≥55共11组。进而，在进行步骤(8)，即将各候选基因组合(即按照上述分组依据确定的各个分组)的所有基因的三种表观基因型比率分别与对应基因的基因表达量进行皮尔逊相关性分析后，再进行步骤(9)的候选基因组合筛选，以便筛选出三种表观基因型比率同时与基因表达量显著相关且三个P值最小的一种候选基因组合作为目标基因组合。例如，接上述的11组分组法，基于11组的皮尔逊相关性分析的结果，选择三种表观基因型同时与基因表达量显著相关且P值最小的一组(即一种候选基因组合)作为目标基因组合。

如前所述，可以利用SPSS 20进行所述多元线性回归。由此，可得该物种候选基因表达受甲基化修饰调控关系的数学模型方程：E＝aX+bY+cZ，E代表基因的表达量，X代表纯合甲基化表观基因型在整个基因中三种表观基因型中的比例，Y代表杂合甲基化表观基因型在整个基因中三种表观基因型中的比例，Z代表纯合非甲基化表观基因型在整个基因中三种表观基因型中的比例，a、b、c为相关系数。

此外，还需要说明的是，本发明的方法具有下列优点的至少之一：

1、现阶段还没有确定预定物种的基因表达和甲基化修饰调控的关系的方法的相关报道，而本发明第一次完成了基因表达受甲基化修饰调控关系的建模，利用本发明的方法能够有效确定预定物种的基因表达和甲基化修饰调控的关系。

2、本发明利用皮尔逊相关性分析与多元线性回归首次建立了基因表达受甲基化修饰调控关系的数学模型，该模型可以为从事DNA甲基化相关研究的科学工作者提供一种有效的工具，换言之，利用本发明的方法可以对不同物种的不同样本进行基因表达受甲基化修饰调控关系的建模，从而能够有效促进表观遗传学的研究。

3、本发明的方法采用已有用户友好度高的生物统计软件，操作简单，可以供不太了解生物统计的人使用，从而节省大量精力完成基因表达受甲基化修饰调控关系的数学建模。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1和图2显示了根据本发明的实施例，待测样本基因组DNA质量检测的琼脂糖凝胶电泳结果；以及

图3显示了根据本发明的实施例，获得毛白杨的基因表达和甲基化修饰调控的关系的数学模型经SPSS 20绘制的多元线性回归图形。

具体实施方式

下面将结合实施例对本发明的方案进行解释。本领域技术人员将会理解，下面的实施例仅用于说明本发明，而不应视为限定本发明的范围。实施例中未注明具体技术或条件的，按照本领域内的文献所描述的技术或条件(例如参考J.萨姆布鲁克等著，黄培堂等译的《分子克隆实验指南》，第三版，科学出版社)或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者，均为可以通过市购获得的常规产品，例如可以采购自Illumina公司。

一般方法：

本发明的确定预定物种的基因表达和甲基化修饰调控的关系的方法，一般包括以下步骤：

(1)以所述预定物种的父本、母本及其子代样本作为待测样本，并对所述待测样本进行重亚硫酸盐全基因组甲基化测序，以便分别获得父本、母本及其子代样本的测序数据，其中各测序数据均由多个读段组成；

(2)将所述子代样本的测序数据与参考基因组进行比对，以便确定所述子代样本中的所有胞嘧啶位点以及包含胞嘧啶位点的读段；

(3)针对所述子代样本的包含胞嘧啶位点的读段，对各胞嘧啶位点均进行位点上下游Call SNPs检测，以便确定所述子代样本的包含胞嘧啶位点的读段中属于等位基因序列的读段，并将每属于等位基因序列的两个读段作为一对目的读段；

(4)将所述子代样本的各目的读段分别与父本和母本的测序数据进行比对，以便确定每一对目的读段的两个读段的亲本来源，并基于每对目的读段的胞嘧啶甲基化信息，确定每对目的读段的表观基因型，

其中，所述目的读段的表观基因型分为三种：纯合甲基化基因型、杂合甲基化基因型和纯合非甲基化基因型，其中目的读段中分别来源于父本和母本的两个读段在同一位置上均发生胞嘧啶甲基化的，该目的片段为纯合甲基化基因型；两个读段在同一位置上均不发生胞嘧啶甲基化的，该目的片段为纯合非甲基化基因型；两个读段在同一位置上只有任意一个发生胞嘧啶甲基化的，该目的片段为杂合甲基化基因型；

(5)基于参考基因组序列，按照序列来源，将多对目的片段分别进行基因归类；

(6)统计各基因包含的目的片段的数目，确定每个基因包含的目的片段中三种表观基因型的比率，并获得各基因的表达量信息；

(7)基于设定候选基因包含的目的片段的数目的不同，对各基因进行多次分组，以便划分出多种候选基因组合；

(8)针对所述多种候选基因组合的每一种，将候选基因组合的所有基因的三种表观基因型比率分别与对应基因的基因表达量进行皮尔逊相关性分析；

(9)基于多种候选基因组合的皮尔逊相关性分析结果，筛选三种表观基因型比率同时与基因表达量显著相关且三个P值最小的一种候选基因组合作为目标基因组合；

(10)基于目标基因组合的表观基因型和基因表达量数据，进行多元线性回归，以便确定所述预定物种的基因表达和甲基化修饰调控的关系。

实施例1

参照本发明的确定预定物种的基因表达和甲基化修饰调控的关系的方法，按照上述的一般方法的步骤，确定毛白杨的基因表达和甲基化修饰调控的关系，具体如下：

其中，在本实例中，毛白杨L9作为杂交父本，LM50作为母本。

在进行测序之前进一步包括对待测样本进行基因组DNA提取的步骤；所述基因组DNA提取采用CTAB法；在得到各待测样本的基因组DNA后，进一步包括对基因组DNA进行质量、纯度和浓度检测，筛选获得合格的DNA样品进行下一步操作，若基因组DNA样品质量不合格需要重新提取基因组DNA。其中，基因组DNA质量检测使用琼脂糖凝胶电泳进行，根据电泳条带分析基因组DNA降解程度以及是否有RNA污染。结果见图1和图2，图1中样品为检测合格的，无降解、无RNA污染样品的电泳条带；图2中样品4号为严重降解；5号有重度RNA污染。本实施例采用筛选质量合格的如图1所示的基因组DNA进行下一步操作。在检测基因组DNA质量合格后，进一步检测其纯度是否合格。其中，本实施例采用Nanodrop检测基因组DNA的纯度，当OD260/280为1.6～1.8时，认为基因组DNA的纯度合格，若基因组DNA的OD260/280不在上述范围内，则所述基因组DNA纯度不合格，需要重新提取。在筛选获得质量和纯度合格的基因组DNA后，进一步采用Qubit对基因组DNA浓度进行精确定量，为后续操作提供依据。

在获得用合格的父母本和子代样本基因组DNA样品后，采用重亚硫酸盐法构建上述合格的待测样本的基因组DNA样品测序文库，进而进行重亚硫酸盐全基因组甲基化测序。构建待测样本基因组DNA样品测序文库时加入比例为建库DNA起始量的1/1000的阴性对照lambda DNA。

其中，构建基因组DNA样品测序文库的步骤如下：首先随机打断基因组DNA至200～300bp获得DNA片段，所述随机打断基因组DNA使用Covaris S220进行；然后对DNA片段进行平末端修复，所述平末端修复使用内切酶和外切酶将黏性末端修复成为平末端；随后在平末端DNA片段后加尾巴A碱基获得带尾巴A的DNA片段；在带尾巴A的DNA片段上连接测序接头后进行Bisulfite处理；所述的Bisulfite处理采用Zymo Research的EZ DNA MethylationGold Kit进行，所述Bisulfite处理的作用是将未发生甲基化的C变成U，在后续的PCR扩增过程后变为T，而甲基化的C则保持不变。在Bisulfite处理后进行PCR扩增即可获得基因组DNA样品测序文库。

进行重亚硫酸盐全基因组甲基化测序时，在获得基因组DNA样品测序文库后，进一步包括对文库进行质检筛选，以便获得合格基因组DNA样品测序文库。在质检文库之前进一步包括使用Qubit2.0对文库DNA浓度进行初步定量，然后根据文库DNA的浓度将其进行稀释，文库稀释后的浓度优选的为1ng/ul；在文库稀释后，使用Agilent2100对文库的插入片段长度进行检测，筛选后的文库插入片段在320～520bp为合格；使用Q-PCR方法对文库的有效浓度进行准确定量，所述文库有效浓度>2nM为合格。

在获得合格基因组DNA样品测序文库后，即可对合格的DNA样品文库进行双末端Hiseq测序，以便获得父本、母本和子代样本的基因组序列，也即测序数据。

步骤(2)进一步包括：利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96将比对结果进行染色体坐标排序、reads去重复处理。

步骤(3)中，通过GATK2(V3.2)对已知胞嘧啶位点附近10bp进行Call SNPs，从而区分等位基因序列(此时，由于序列经由重硫酸盐甲基化测序，所以胞嘧啶上的甲基化修饰状况已被记录)。

在步骤(5)中，将多对目的片段分别进行基因归类包括：将所述多对目的片段(也即胞嘧啶甲基化分型片段)通过比对参考基因组序列进行注释，将不同片段归类于所属基因条目下，完成目的片段(也即胞嘧啶表观基因型分型的片段)与基因的mapping，最终不同基因获得的mapping成功的目的片段的数目是不同的。例如，由Potri.007G001500的4个到Potri.011G014300的57个胞嘧啶表观基因型分型片段。

在步骤(6)中，统计各基因包含的目的片段的数目，确定每个基因包含的目的片段中三种表观基因型的比率，包括：统计同一基因下不同类型的表观基因型的数量以及它们的比例，例如：某基因下mapping成功的目的片段共57个目的片段(有时也称为“胞嘧啶表观基因型分型片段”)中纯合甲基化基因型(有时也称为“表观基因型”)为30个(则，其比例＝30/57＝52.6％)，杂合甲基化基因型为15个(比例＝15/57＝26.3％)，而纯合非甲基化基因型为12个(比例＝12/57＝21.1％)。

在步骤(7)中，基于设定候选基因包含的目的片段的数目的不同，对各基因进行多次分组。在本实施例中，基于候选基因包含的目的片段的数目≥5，≥10，≥15，≥20，≥25，≥30，≥35，≥40，≥45，≥50，≥55分为11组。

进而，在进行步骤(8)，即将各候选基因组合(即按照上述分组依据确定的各个分组)的所有基因的三种表观基因型比率分别与对应基因的基因表达量进行皮尔逊相关性分析后，再进行步骤(9)的候选基因组合筛选，以便筛选出三种表观基因型比率同时与基因表达量显著相关且三个P值最小的一种候选基因组合作为目标基因组合。其中，各分组的皮尔逊相关性分析结果见表1。

表1

表1中“**”代表P<0.01，“*”代表P<0.05。

基于表1的11组的皮尔逊相关性分析的结果可知，≥10和≥15时的两组数据，三种表观基因型都与基因表达量显著相关，此时，得到三种胞嘧啶表观基因型比率与基因表达量同时显著相关的数据组，选取P值最小的一组数据：≥10组数据，三种表观基因型在mapping基因中的比例都与基因表达量极显著相关，即P<0.01，而≥15组数据，只有纯合非甲基化表观基因型在mapping基因中所占比例与基因表达量极显著相关，而纯合甲基化表观基因型与杂合甲基化表观基因型在mapping基因中所占比例与基因表达量显著相关，即P<0.05。也即，在≥10和≥15时的两组数据中，≥10组数据有两个P值小于≥15组数据，因此，选择≥10组数据，也即选择候选基因包含的目的片段的数目≥10的分组——候选基因组合作为目标基因组合(见表1)。

进而，将候选基因包含的目的片段的数目≥10的候选基因组合的表观基因型和基因表达量数据，利用SPSS 20进行多元线性回归。由此，可得该物种候选基因表达受甲基化修饰调控关系的数学模型方程：E＝aX+bY+cZ，E代表基因的表达量，X代表纯合甲基化表观基因型在整个基因中三种表观基因型中的比例，Y代表杂合甲基化表观基因型在整个基因中三种表观基因型中的比例，Z代表纯合非甲基化表观基因型在整个基因中三种表观基因型中的比例，a、b、c为相关系数。

结果，得到可计算的基因表达受甲基化修饰调控关系的数学模型E＝-0.049X-0.024Y+0.036Z方程中E代表基因的表达量，X代表纯合甲基化表观基因型在整个基因中三种表观基因型中的比例，Y代表杂合甲基化表观基因型在整个基因中三种表观基因型中的比例，Z代表纯合非甲基化表观基因型在整个基因中三种表观基因型中的比例(如图3所示)。图3为数学模型经SPSS 20绘制的多元线性回归图形。

由此，通过上述获得的数学模型能够有效了解毛白杨关键基因(即候选基因)的甲基化修饰状况，进而，通过可调节毛白杨关键基因的甲基化水平的试剂(例如5氮杂胞苷，等)或其他手段，即可有效调节该关键基因的表达，从而达到调控毛白杨特定生理生化过程的目的。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。

Claims (4)

1.一种确定预定物种的基因表达和甲基化修饰调控的关系的方法，其特征在于，包括以下步骤：

(1)以所述预定物种的父本、母本及其子代样本作为待测样本，并对所述待测样本进行重亚硫酸盐全基因组甲基化测序，以便分别获得父本、母本及其子代样本的测序数据，其中各测序数据均由多个读段组成；

(2)将所述子代样本的测序数据与参考基因组进行比对，以便确定所述子代样本中的所有胞嘧啶位点以及包含胞嘧啶位点的读段；

(3)针对所述子代样本的包含胞嘧啶位点的读段，对各胞嘧啶位点均进行位点上下游Call SNPs检测，以便确定所述子代样本的包含胞嘧啶位点的读段中属于等位基因序列的读段，并将每属于等位基因序列的两个读段作为一对目的读段；

(4)将所述子代样本的各目的读段分别与父本和母本的测序数据进行比对，以便确定每一对目的读段的两个读段的亲本来源，并基于每对目的读段的胞嘧啶甲基化信息，确定每对目的读段的表观基因型，

其中，所述目的读段的表观基因型分为三种：纯合甲基化基因型、杂合甲基化基因型和纯合非甲基化基因型，其中目的读段中分别来源于父本和母本的两个读段在同一位置上均发生胞嘧啶甲基化的，该目的片段为纯合甲基化基因型；两个读段在同一位置上均不发生胞嘧啶甲基化的，该目的片段为纯合非甲基化基因型；两个读段在同一位置上只有任意一个发生胞嘧啶甲基化的，该目的片段为杂合甲基化基因型；

(5)基于参考基因组序列，按照序列来源，将多对目的片段分别进行基因归类；

(6)统计各基因包含的目的片段的数目，确定每个基因包含的目的片段中三种表观基因型的比率，并获得各基因的表达量信息；

(7)基于设定候选基因包含的目的片段的数目的不同，对各基因进行多次分组，以便划分出多种候选基因组合；

(8)针对所述多种候选基因组合的每一种，将候选基因组合的所有基因的三种表观基因型比率分别与对应基因的基因表达量进行皮尔逊相关性分析；

(9)基于多种候选基因组合的皮尔逊相关性分析结果，筛选三种表观基因型比率同时与基因表达量显著相关且三个P值最小的一种候选基因组合作为目标基因组合；

(10)基于目标基因组合的表观基因型和基因表达量数据，进行多元线性回归，以便确定所述预定物种的基因表达和甲基化修饰调控的关系。

2.根据权利要求1所述的方法，其特征在于，步骤(2)进一步包括：利用SAMTOOLSV0.1.18和PICARD-TOOLSV1.96将比对结果进行染色体坐标排序、reads去重复处理。

3.根据权利要求1所述的方法，其特征在于，利用GATK2对胞嘧啶位点附近5-15bp进行Call SNPs。

4.根据权利要求1所述的方法，其特征在于，利用SPSS 20进行所述多元线性回归。

Images