基于DNA甲基化水平的木本植物叶片表型特征和光合特性预测模型的构建方法及预测方法

摘要

本发明提供了基于DNA甲基化水平的木本植物叶片性状和光合特性模型的构建方法及预测方法，属于生物分析技术领域。本发明基于随机森林选取体现地理位置差异的重要特征变量，筛选得到7个叶片特征变量，确定最优聚类数目，利用改进的FCM聚类算法得到每组聚类叶片样本；根据变量间相关性和梯度提升树得的酶切组合重要性，获得对每组聚类叶片样本中重要酶切组合；以所述酶切组合的DNA甲基化水平作为回归变量，基于高斯径向基函数构建LS‑SVM回归预测模型；输入重要酶切组合的DNA甲基化水平来准确地预测叶形状因子，叶面积和净光合速率。本方法用于预测木本植物表型特征和光合特性，同时筛选优良性状木本植物个体。

权利要求

1.一种基于DNA甲基化水平的木本植物的叶片表型特征和光合特性预测模型的构建方法，包括以下步骤：

1)收集全国自然分布的同一物种木本植物的叶片，得到叶片代表样本；

2)测定所述叶片代表样本的表型特征和光合特性，得到叶片的光合特性数据和叶片表型特征数据；

所述叶片表型特征包括叶面积、叶长度、叶宽度、叶周长、叶长宽比和叶形状因子；

所述光合特性包括净光合速率、气孔导度、CO2浓度和水分利用率；

3)测定所述叶片代表样本中酶切片段的DNA甲基化状态，得到酶切片段的DNA甲基化水平，计算得到每个叶片代表样本的全基因组平均DNA甲基化水平；

4)以所述叶片的光合特性数据、叶片表型特征数据和叶片的全基因组的平均DNA甲基化水平为候选变量，采用随机森林模型筛选对地理位置产生差异的重要特征变量，得到叶片表型特征数据和全基因组平均DNA甲基化水平；

5)将所述叶片表型特征数据和全基因组平均DNA甲基化水平作为特征变量，利用NbClust软件包中组内平方误差和26个聚类评估指标确定叶片代表样本的最优聚类组数；

6)将所述叶片代表样本的最优聚类组数输入到改进的模糊C-均值聚类算法中，计算得到每组聚类样本的隶属度矩阵；

所述改进的模糊C-均值聚类算法具体步骤如下：

a.给定所述最优聚类组数c，n是样本数据个数，设定迭代停止阈值为ε＝10-5，设置最大迭代次数t＝300，随机初始化隶属矩阵U，令迭代计数器t＝0；

设有限集合X＝{x1,x2,...,xn}，且X中的元素有m个特征变量，X表示为n×m的矩阵如下：其中，m表示特征变量的数目，n代表叶片代表样本个数；

把矩阵X的n个样本分为c组，其中2≤c≤n，分成的c组的模糊聚类矩阵U为：所述矩阵U中，μij表示样本xj与类别i的隶属关系，且0≤μij≤1,c个聚类中心为：

选用最小误差平方和作为聚类准则，聚类分析的目标函数如式(1)所示：加上约束条件如式(2)所示公式：

求解得式(3)所示公式：

b.根据式(3)计算更新模糊聚类矩阵和聚类中心矩阵；

c.若P(t)-P(t-1)＜ε，则停止计算，并输出模糊聚类矩阵U和聚类中心矩阵P，否则令t＝t+

1，转向步骤(2)直至输出矩阵U和矩阵P；

7)基于所述步骤6)得到的每组聚类样本的隶属度矩阵，计算每组聚类样本中全基因组平均DNA甲基化水平与叶片各特征变量的相关性和梯度提升树得到的酶切组合重要性，得到每一组聚类样本的重要酶切片段组合；

所述叶片各特征变量包括叶面积、净光合速率和叶形状因子；

8)以所述步骤7)得到的重要酶切片段组合的DNA甲基化水平作为回归变量，利用高斯径向基函数建立LS-SVM回归预测模型，得到叶片表型特征和光合特性模型如式(9)所示；

步骤2和3之间没有时间顺序的限制。

2.根据权利要求1所述的构建方法，其特征在于，所述步骤4)中采用随机森林模型筛选为选择特征变量的Mean Decrease Accuracy和Mean Decrease Gini的均值在15以上的变量作为重要特征变量。

3.根据权利要求1所述的构建方法，其特征在于，利用高斯径向基函数建立LS-SVM回归预测模型的方法，包括以下步骤：在SVM中，假设样本训练集为T＝{ti|ti＝(xi,yi),xi∈Rn, 样本训练集T中，xi为训练样本集中第i个样本的的输入变量，yi为训练样本集中第i个样本的的输出变量，R代表实数域，n代表输入样本个数，回归函数为式(4)中，w和b为回归参数， 为特征映射，x为训练样本集的输入变量；

而在LS-SVM中给问题转化为求解：

式(5)中，γ为正则化参数，ξ为松弛因子，In×1＝(1,1,...,1)′，In×1＝(1,1,...,1)′为这样的一个矩阵；

构造Lagrange函数如下：

其中，α为拉格朗日乘子；下求L(w,b,ξ,α)的鞍点，即最优点；

消去式(7)中的w，ξ，可得：

式(8)中，Ω＝(xix′j),i,j＝1,2,...,n，E为n×n的单位矩阵；

通过解式(8)得到α和b，则最小二乘支持向量机的估计函数为：其中，k(xi,x′j)为核函数，选取Gauss径向基核函数在进行回归前已对输入变量进行了标准化，并通过参数寻优，令γ＝10，σ＝1。

4.根据权利要求1所述构建方法，其特征在于，所述步骤3)中叶片代表样本中酶切片段的DNA甲基化状态的测定方法，包括以下步骤：

31)采用EcoRI/HpaII和EcoRI/MspI限制性内切酶对所述叶片代表样本的基因组DNA进行酶切，得到的酶切片段；

32)将所述酶切片段依次进行预扩增和选择性扩增，根据得到的选择性扩增产物进行分型，根据分型结果得到酶切片段的DNA甲基化状态。

5.根据权利要求4所述构建方法，其特征在于，所述分型是将所述选择性扩增产物进行电泳，将得到的电泳条带在二进制字符矩阵中进行评分，用“0”表示条带缺失，用“1”表示条带的存在；CNG(1，0)代表半甲基化状态，CG(0，1)代表全甲基化状态，(1，1)代表无甲基化状态，(0，0)代表未知甲基化状态；全基因组的甲基化水平计算公式如式(10)所示：全基因组的DNA甲基化水平＝(半甲基化状态的位点+全甲基化状态位点+未知甲基化状态的位点)/(半甲基化状态的位点+全甲基化状态位点+未知甲基化状态的位点+无甲基化状态的位点)    式(10)。

6.根据权利要求1所述构建方法，其特征在于，所述叶片代表样本的数量为200个以上。

7.根据权利要求1所述构建方法，其特征在于，所述叶形状因子的计算公式如式(11)所示：

8.一种基于DNA甲基化水平的木本植物的叶片表型特征和光合特性预测模型的预测方法，其特征在于，用权利要求1～7任意一项所述方法中的重要酶切片段组合的DNA甲基化水平输入权利要求1～7任意一项所述方法构建的叶片表型特征和光合特性模型中预测叶形状因子、叶面积和叶净光合速率；

所述重要酶切片段组合的DNA甲基化水平为所述权利要求1～7任意一项所述构建方法中得到的重要酶切片段组合按照全基因组的甲基化水平计算公式计算得到的DNA甲基化水平。

9.根据权利要求8所述预测方法，其特征在于，酶切组合的DNA甲基化水平对叶型因子，叶面积和净光合速率的影响程度绘制边际效应图；

所述边际效应图的绘制方法为将重要酶切片段组合的DNA甲基化水平作为输入变量，然后调用gbm中自带的绘制边际效应图的函数绘制较重要的酶切片段对每个叶片的叶形状因子、叶面积和叶净光合速率的边际效应图。

说明书

基于DNA甲基化水平的木本植物叶片表型特征和光合特性预

测模型的构建方法及预测方法

技术领域

[0001] 本发明属于生物信息领域，具体涉及基于DNA甲基化水平的木本植物叶片表型特征和光合特性预测模型的构建方法及预测方法。

背景技术

[0002] DNA甲基化最常发生于二核苷酸内的胞嘧啶的第五位碳原子上，是一种重要的表观遗传学修饰。DNA甲基化抑制转座因子，假基因，重复序列和个体基因的表达，在许多生物学过程如基因表达，胚胎发育，细胞增殖，分化和染色体稳定性中起关键作用。在植物中，DNA甲基化发生于CG、CHG、CHH位点(H代表除鸟嘌呤之外的其他碱基)，并且在高等植物细胞中，被甲基化的胞嘧啶最多可以占到总胞嘧啶数的50％。DNA甲基化控制植物的生长和发育，特别涉及基因表达和DNA复制的调控。因此，研究植物的DNA甲基化有助于我们了解DNA甲基化调节植物生长发育的方式，具有现实意义。

[0003] 以往对DNA甲基化的研究主要集中在定性分析，缺少定量研究，且对于植物的DNA甲基化研究也主要集中在草本植物上，对木本植物的研究还较少。例如，2015年，Wanneng Yang等人基于线性回归研究了全基因组DNA甲基化如何影响水稻的叶片性状(叶宽，叶长，叶面积等)。此外，2015年，Dong Ci等人报道了利用甲基化敏感多态性技术(MSAP)对小叶杨自然群体的DNA甲基化修饰位点进行全基因组扫描，获得甲基化多态性位点；同时利用主成分分析(PCA)和STRUCTURE软件对小叶杨自然群体的表观群体遗传结构进行解析。发现DNA甲基化位点上的基因可能在叶片发育和光合作用调控中起重要作用，对植物性状(包括叶)相关联形态和光合作用有一定影响。但是该研究仍旧停留在定性分析上，缺少对DNA甲基化分析的定量研究。同时，目前对林木植物关于DNA甲基化的研究依赖于全基因组的DNA甲基化的扫描，研究成本较高，且数据量庞大导致结果准确性较差。

发明内容

[0004] 有鉴于此，本发明的目的在于提供基于DNA甲基化水平的木本植物叶片表型特征和光合特性预测模型的构建方法及预测方法，预测准确率高。

[0005] 为了实现上述发明目的，本发明提供以下技术方案：

[0006] 本发明提供了基于DNA甲基化水平的木本植物叶片表型特征和光合特性预测模型的构建方法，包括以下步骤：

[0007] 1)收集全国自然分布的同一物种木本植物的叶片，得到叶片代表样本；

[0008] 2)测定所述叶片代表样本的光合特性和表型特征，得到叶片的光合特性数据和叶片表型特征数据；

[0009] 所述表型特征包括叶面积、叶长度、叶宽度、叶周长、叶长宽比和叶形状因子；

[0010] 所述光合特性包括净光合速率、气孔导度、CO2浓度和水分利用率；

[0011] 3)测定所述叶片代表样本中酶切片段的DNA甲基化状态，得到酶切片段的DNA甲基化水平，计算得到每个叶片代表样本的全基因组平均DNA甲基化水平；

[0012] 4)以所述叶片的光合特性数据、叶片表型特征数据和叶片的全基因组的平均DNA甲基化水平为候选变量，采用随机森林模型筛选对地理位置产生差异的重要特征变量，得到叶片表型特征数据和全基因组平均DNA甲基化水平；

[0013] 5)将所述叶片表型特征数据和全基因组平均DNA甲基化水平作为特征变量，利用NbClust软件包中组内平方误差和26个聚类评估指标确定叶片代表样本的最优聚类组数；

[0014] 6)将所述叶片代表样本的最优聚类组数输入到改进的模糊C-均值聚类算法中，计算得到每组聚类样本的隶属度矩阵；

[0015] 所述改进的模糊C-均值聚类算法具体步骤如下：

[0016] a.给定所述最优聚类组数c，n是样本数据个数，设定迭代停止阈值为105，设置最大迭代次数t 300，随机初始化隶属矩阵U，令迭代计数器t 0；

[0017] 设有限集合X{x1,x2,...,xn}，且X中的元素有m个特征变量，X表示为n m的矩阵如下：

[0018]

[0019] 其中，m表示特征变量的数目，n代表叶片代表样本个数；

[0020] 把矩阵X的n个样本分为c组(2 c n)，分成的c组的模糊聚类矩阵U为：

[0021]

[0022] 所述矩阵U中，ij表示样本xj与类别i的隶属关系，且 c个聚类中心为：

[0023]

[0024] 选用最小误差平方和作为聚类准则，聚类分析的目标函数如式(1)所示：

[0025]

[0026] 加上约束条件如式(2)所示公式：

[0027]

[0028] 求解得式(3)所示公式：

[0029]

[0030] b.根据式(3)计算更新模糊聚类矩阵和聚类中心矩阵；

[0031] c.若P(t)P(t1)，则停止计算，并输出模糊聚类矩阵U和聚类中心矩阵P，否则令t t 1，转向步骤(2)直至输出矩阵U和矩阵P；

[0032] 7)基于所述步骤6)得到的每组聚类样本的隶属度矩阵，计算每组聚类样本中全基因组平均DNA甲基化水平与叶片各特征变量的相关性和梯度提升树得到的酶切组合重要性，得到每一组聚类样本的重要酶切片段组合；

[0033] 所述叶片各特征变量包括叶面积、净光合速率和叶形状因子；

[0034] 8)以所述步骤7)得到的重要酶切片段组合的DNA甲基化水平作为回归变量，利用高斯径向基函数建立LS-SVM回归预测模型，得到叶片表型特征和光合特性模型如式(9)所示；

[0035]

[0036] 步骤2和3之间没有时间顺序的限制。

[0037] 优选的，所述步骤4)中采用随机森林模型筛选为选择选择特征变量的Mean Decrease Accuracy和Mean Decrease Gini的均值在15以上的变量作为重要特征变量。

[0038] 优选的，利用高斯径向基函数建立LS-SVM回归预测模型的方法，包括以下步骤：

[0039] 在SVM中，假设样本训练集为 样本训练集T中，xi为训练样本集中第i个样本的的输入变量yi为训练样本集中第i个样本的输出变量，R代表实数域，n代表输入样本个数，回归函数为

[0040] f(x)w(x)b                 ⑷

[0041] 式(4)中，w和b为回归参数，(x)为特征映射，x为训练样本集的输入变量；

[0042] 而在LS-SVM中给问题转化为求解：

[0043]

[0044] s.t.y (x)w bIξ                ⑸

[0045] 式(5)中，为正则化参数，ξ为松弛因子，In1(1,1,...,1)，In1(1,1,...,1)为 这样的一个矩阵；

[0046] 构造Lagrange函数如下：

[0047]

[0048] 其中α为拉格朗日乘子；下求L(w,b,ξ,α)的鞍点，即最优点；

[0049]

[0050] 消去式(7)中的w，ξ，可得：

[0051]

[0052] 式(8)中，Ω(xixj),i,j 1,2,...,n，E为n n的单位矩阵。

[0053] 通过解式(8)得到α和b，则最小二乘支持向量机的估计函数为：

[0054]

[0055] 其中，k(xi,xj)为核函数，选取Gauss径向基核函数k(xi,xj)

[0056] 在进行回归前已对输入变量进行了标准化，并通过参数寻优，令10，1。

[0057] 优选的，所述步骤3)中叶片代表样本中酶切片段的DNA甲基化状态的测定方法，包括以下步骤：

[0058] 31)采用EcoRI/HpaII和EcoRI/MspI限制性内切酶对所述叶片代表样本的基因组DNA进行酶切，得到的酶切片段；

[0059] 32)将所述酶切片段依次进行预扩增和选择性扩增，根据得到的选择性扩增产物进行分型，根据分型结果得到酶切片段的DNA甲基化状态。

[0060] 优选的，所述分型是将所述选择性扩增产物进行电泳，将得到的电泳条带在二进制字符矩阵中进行评分，用“0”表示条带缺失，用“1”表示条带的存在；CNG(1，0)代表半甲基化状态，CG(0，1)代表全甲基化状态，(1，1)代表无甲基化状态，(0，0)代表未知甲基化状态；全基因组的甲基化水平计算公式如式(10)所示：

[0061] 全基因组的DNA甲基化水平＝(半甲基化状态的位点+全甲基化状态位点+未知甲基化状态的位点)/(半甲基化状态的位点+全甲基化状态位点+未知甲基化状态的位点+无甲基化状态的位点)式(10)。

[0062] 优选的，所述叶片代表样本的数量为200个以上。

[0063] 优选的，所述叶形状因子的计算公式如式(11)所示：

[0064]

[0065] 本发明提供了基于DNA甲基化水平的木本植物叶片表型特征和光合特性预测模型的预测方法，用重要酶切片段组合的DNA甲基化水平输入所述方法构建的叶片表型特征和光合特性模型中预测叶形状因子、叶面积和叶净光合速率；

[0066] 所述重要酶切片段组合的DNA甲基化水平为所述所述构建方法中得到的重要酶切片段组合按照全基因组的甲基化水平计算公式计算得到的DNA甲基化水平。

[0067] 优选的，酶切组合的DNA甲基化水平对叶型因子，叶面积和净光合速率的影响程度绘制边际效应图。

[0068] 所述边际效应图的绘制方法为将重要酶切片段组合的DNA甲基化水平作为输入变量，然后调用gbm中自带的绘制边际效应图的函数绘制较重要的酶切片段对每个叶片的叶形状因子、叶面积和叶净光合速率的边际效应图。

[0069] 本发明提供的基于DNA甲基化水平的木本植物叶片表型特性和光合特性预测模型的构建方法，将叶片代表样本在聚类前，使用随机森林进行光合特性和表型特征选择，得到对所述木本植物地理分布差异有较大影响的特征变量。叶片代表样本聚类时，未使用传统的聚类方法(K-Means聚类，PAM聚类等)，而是使用的改进的模糊C-均值聚类算法，所得输出结果具有概率意义，保留了木本植物个体间差异的不确定性。同时在建立预测模型前，首先筛选较重要的酶切组合，降低了预测模型和实际操作的复杂度，使预测更准确。同时本发明首次建立了林木表观群体遗传学研究策略，解析了林木群体表观遗传结构，提出DNA甲基化可能会影响小叶杨自然群体的表型特征和光合特性。

[0070] 进一步的，本发明提供的方法，利用筛选的重要酶切片段组合的DNA甲基化水平对表型特征和光合特性做出边际效应图，实现定量研究DNA甲基化对表型特征和光合特性的影响。

附图说明

[0071] 图1为实施例1中采用森林随机模型筛选小叶杨特征变量图；

[0072] 图2为实施例1中基于SSE and 26聚类评估指数确定小叶杨样本最优聚类个数图；

[0073] 图3为实施例1中使用改进的FCM算法确定每类样本的个体图；

[0074] 图4为实施例1中各特征变量之间相关性聚类情况图；

[0075] 图5为实施例1中第一组聚类小叶杨样本中酶切片段与三个变量的重要性；

[0076] 图6为实施例1中第一组聚类小叶杨样本中酶切片段与三个变量的重要性；

[0077] 图7为实施例2中重要酶切组合的DNA甲基化水平对叶面积的边际效用的影响，图7-1为第一组聚类小叶杨样本；图7-2为第二聚类小叶杨样本；

[0078] 图8为实施例2中重要酶切组合的DNA甲基化水平对净光合速率的边际效用的影响，图8-1为第一组聚类小叶杨样本；图8-2为第二聚类小叶杨样本；

[0079] 图9为实施例2中重要酶切组合的DNA甲基化水平对叶片形状因子的边际效用的影响，图9-1为第一组聚类小叶杨样本；图9-2为第二聚类小叶杨样本；

[0080] 图10为实施例2中两个杨树亚群的叶片形态。

具体实施方式

[0081] 本发明提供了基于DNA甲基化水平的木本植物叶片表型特征和光合特性预测模型的构建方法，包括以下步骤：

[0082] 1)收集全国自然分布的一种木本植物的叶片，形成叶片代表样本；

[0083] 2)测定所述叶片代表样本的光合特性和表型特征，得到叶片的光合特性数据和叶片表型特征数据；

[0084] 所述表型特征包括叶面积、叶长度、叶宽度、叶周长、叶长宽比和叶形状因子；

[0085] 所述光合特性包括净光合速率、气孔导度、CO2浓度和水分利用率；

[0086] 3)测定所述叶片代表样本中酶切片段的DNA甲基化状态，得到酶切片段的DNA甲基化水平，计算得到每个叶片代表样本的全基因组平均DNA甲基化水平；

[0087] 4)以所述叶片的光合特性数据、叶片表型特征数据和叶片的全基因组的平均DNA甲基化水平为候选变量，采用随机森林模型筛选对地理位置产生差异的重要特征变量，得到叶片表型特征和全基因组平均DNA甲基化水平；

[0088] 5)将叶片表型特征数据和全基因组平均DNA甲基化水平作为特征变量，利用NbClust软件包中组内平方误差和26个聚类评估指标确定叶片代表样本的最优聚类组数；

[0089] 6)将所述叶片代表样本的最优聚类数输入到改进的模糊C-均值聚类算法中，计算得到每组聚类样本的隶属度矩阵；

[0090] 所述改进的模糊C-均值聚类算法具体步骤如下：

[0091] a.给定所述最优聚类数c，n是样本数据个数，设定迭代停止阈值为105，设置最大迭代次数t 300，随机初始化隶属矩阵U，令迭代计数器t 0；

[0092] 设有限集合X{x1,x2,...,xn}，且X中的元素有m个特征变量，X表示为n m的矩阵如下：

[0093]

[0094] 其中，m表示特征变量的数目，n代表叶片代表样本数目；

[0095] 把矩阵X的n个样本分为c类(2 c n)，其模糊聚类矩阵U为：

[0096]

[0097] 其中，ij表示样本xj与类别i的隶属关系，且

[0098] c个聚类中心为：

[0099]

[0100] 选用最小误差平方和作为聚类准则，聚类分析的目标函数为：

[0101]

[0102] 加上约束条件得：

[0103]

[0104] 求解得：

[0105]

[0106] b.根据式(3)计算更新模糊聚类矩阵和聚类中心矩阵；

[0107] c.若P(t)P(t1)，则停止计算，并输出模糊聚类矩阵U和聚类中心矩阵P，否则令t t 1，转向步骤(2)直至得到模糊聚类矩阵U和聚类中心矩阵P；

[0108] 7)基于每组聚类样本的隶属度矩阵，计算每组聚类样本中全基因组平均DNA甲基化水平与叶片各特征变量的相关性和梯度提升树得到的酶切组合重要性，得到每一组聚类样本的重要酶切片段组合；

[0109] 所述叶片各特征变量包括叶面积、净光合速率和叶形状因子；

[0110] 8)以重要酶切片段组合的DNA甲基化水平作为回归变量，利用高斯径向基函数建立LS-SVM回归预测模型，用重要酶切片段组合的DNA甲基化水平输入LS-SVM回归预测模型中预测叶形状因子、叶面积和叶净光合速率。

[0111] 本发明收集全国自然分布的同一物种木本植物的叶片，得到叶片代表样本。

[0112] 在本发明中，叶片的收集地优选陕西、青海、河北、河南、宁夏、山西、北京、内蒙古。所述木本植物的种类没有特殊限制，本发明提供的方法可适用所有的木本植物。所述木本植物优选为杨树属，最优选为小叶杨。木本植物的叶片样本的数量优选为1000个以上。从采集的木本植物的叶片中挑选叶片代表样本。所述挑选的标准是能够覆盖自然小叶杨群体的整个地理分布。所述叶片代表样本的数量为200～500个。

[0113] 得到叶片代表样本后，本发明测定所述叶片代表样本的光合特性和表型特征，得到叶片的光合特性数据和叶片表型特征数据；所述叶片表型特征包括叶面积、叶长度、叶宽度、叶周长、叶长宽比和叶形状因子；所述光合特性包括净光合速率、气孔导度、CO2浓度和水分利用率。

[0114] 在本发明中，所述测量叶面积、叶长度、叶宽度、叶周长和叶宽度比五个表型性状时，采用了激光叶片面积测量装置测定。本发明对激光叶片面积测量装置没有特殊限制，采用本领域所熟知的叶面积测量仪即可。本发明实施例中，叶面积测量仪为便携式激光叶片面积测量装置(CI-202)。在本发明中，所述叶形状因子的的公式如式(11)所示：

[0115]

[0116] 在本发明中，所述光合特性采用便携式气体交换系统(Li-6400xt,LiCor)仪器测定。为了获得最大瞬时净光合速率，将光合光子通量密度(PPFD)设置为1600，CO2浓度设置为400。净光合速率、气孔导度、细胞间CO2浓度和水分利用效率(WUE)在净光合速率下均有记录。

[0117] 本发明测定所述叶片代表样本中酶切片段的DNA甲基化状态，得到酶切片段的DNA甲基化水平，计算得到每个叶片代表样本的全基因组平均DNA甲基化水平。

[0118] 在本发明中，所述叶片代表样本中酶切片段的DNA甲基化状态的测定方法，优选包括以下步骤：

[0119] 31)采用EcoRI/HpaII和EcoRI/MspI限制性内切酶分别对所述叶片代表样本的基因组DNA进行酶切，得到的酶切片段；

[0120] 32)将所述酶切片段依次进行预扩增和选择性扩增，根据得到的选择性扩增产物进行分型，根据分型结果得到酶切片段的DNA甲基化状态。

[0121] 本发明优选提取所述叶片代表样本的基因组DNA。

[0122] 在本发明中，所述叶片代表样本的基因组DNA的提取方法没有特殊限制，采用本领域所熟知的基因组DNA提取方法即可。在本发明实施例中，所述叶片代表样本的基因组DNA的提取采用试剂盒法提取得到。所述试剂盒采用DNA植物小型试剂盒(基根中国，上海)。

[0123] 得到所述叶片代表样本的基因组DNA后，本发明采用EcoRI/HpaII和EcoRI/MspI限制性内切酶分别对所述叶片代表样本的基因组DNA进行酶切，得到的酶切片段。

[0124] 在本发明中，所述EcoRI/HpaII和EcoRI/MspI限制性内切酶的来源没有特殊限制，采用本领域所熟知的酶来源即可。本发明对所述酶切的方法没有特殊限制，采用本领域所熟知的EcoRI/HpaII和EcoRI/MspI限制性内切酶的酶切方法即可。

[0125] 在本发明中，所述酶切片段进行预扩增和选择性扩增，得到选择性扩增PCR产物。本发明对所述预扩增和选择性扩增没有特殊限制，采用Dong Ci等人发表的Variation in genomic methylation in natural populations of Populus simonii is associated with leaf shape and photosynthetic traits(Journal of Experimental Botany,Vol.67,No.3pp.723–737,2016)。

[0126] 得到的选择性扩增产物后，本发明对所述选择性扩增产物进行分型，根据分型结果得到酶切片段的DNA甲基化状态。

[0127] 在本发明中，所述选择性扩增产物为EcoRI/HpaII限制性内切酶酶切产物条带时用H表示；所述选择性扩增产物为EcoRI/MspI限制性内切酶酶切产物条带时用M表示。若选择性扩增的PCR产物进行电泳后，HM均有带，则说明两者切割处都未发生甲基化，即CCGG序列；若H有带，M无带时，是半甲基化也有说是外甲基化，即5‘mCCGG序列；而H无带，M有带时，是全甲基化，也说是内甲基化，即5’CmCGG序列。HpaII不能切割任何双链胞嘧啶全甲基化，只能切割单链外甲基化。至于MspI，它可以切割内侧甲基化位点，既可以是全甲基化，也可以是半甲基化。

[0128] 得到的选择性扩增产物后，本发明对所述选择性扩增产物进行分型，根据分型结果得到酶切片段的DNA甲基化状态。

[0129] 在本发明中，所述分型优选将所述选择性扩增产物进行电泳，将得到的电泳条带在二进制字符矩阵中进行评分，用“0”表示条带缺失，用“1”表示条带的存在；(CNG(1，0))代表半甲基化状态，(CG(0，1)代表全甲基化状态，(1，1)代表全甲基化状态，(0，0)代表未知甲基化状态)全基因组的甲基化水平计算公式为式(1，0)所示：

[0130] 全基因组的DNA甲基化水平＝(半甲基化状态的位点+全甲基化状态位点+未知甲基化状态的位点)/(半甲基化状态的位点+全甲基化状态位点+未知甲基化状态的位点+无甲基化状态的位点)     式(10)。

[0131] 得到叶片的光合特性数据、叶片表型特征数据和叶片的全基因组的平均DNA甲基化水平后，本发明以所述叶片的光合特性数据、叶片表型特征数据和叶片的全基因组的平均DNA甲基化水平为候选变量，采用随机森林模型筛选对地理位置产生差异的重要特征变量，得到叶片表型特征和全基因组平均DNA甲基化水平。

[0132] 在本发明中，所述采用随机森林模型筛选对地理位置产生差异的重要特征变量的方法，是将Mean Decrease Accuracy和Mean Decrease Gini的均值在15以上的变量作为重要特征变量。本申请以小叶杨为样本，得到的重要特征变量是所述叶片的叶片表型特征数据和叶片的全基因组的平均DNA甲基化水平。

[0133] 得到叶片表型特征和全基因组平均DNA甲基化水平后，本发明将叶片表型特征数据和全基因组平均DNA甲基化水平作为特征变量，利用NbClust软件包中组内平方误差结合26个聚类评估指标确定叶片代表样本的最优聚类组数。

[0134] 在本发明中，使用组内平方误差(SSE)评估时，选择斜率最小的聚类数，且使用26个聚类评估指标评估时，选择组内标准数(Number Criteria)最大的聚类数为最优聚类组数。

[0135] 本发明将所述叶片代表样本的最优聚类数输入到改进的模糊C-均值聚类算法中，计算得到每组聚类样本的隶属度矩阵；

[0136] 所述改进的模糊C-均值聚类算法具体步骤如下：

[0137] a.给定所述最优聚类组数c，n是样本数据个数，设定迭代停止阈值为105，设置最大迭代次数t 300，随机初始化隶属矩阵U，令迭代计数器t 0；

[0138] 设有限集合X{x1,x2,...,xn}，且X中的元素有m个特征变量，X表示为n m的矩阵如下：

[0139]

[0140] 其中，m表示特征变量的数目，n代表叶片代表样本数目；

[0141] 把矩阵X的n个样本分为c类(2 c n)，其模糊聚类矩阵U为：

[0142]

[0143] 其中，ij表示样本xj与类别i的隶属关系，且

[0144] c个聚类中心为：

[0145]

[0146] 选用最小误差平方和作为聚类准则，聚类分析的目标函数为：

[0147]

[0148] 加上约束条件得：

[0149]

[0150] 求解得：

[0151]

[0152] b.根据式(3)计算更新模糊聚类矩阵和聚类中心矩阵；

[0153] c.若P(t)P(t1)，则停止计算，并输出模糊聚类矩阵U和聚类中心矩阵P，否则令t t 1，转向步骤(2)直至得到模糊聚类矩阵U和聚类中心矩阵P。

[0154] 得到每组聚类样本的隶属度矩阵后，本发明基于每组聚类样本的隶属度矩阵，计算每组聚类样本中全基因组平均DNA甲基化水平与叶片各特征变量的相关性，以及根据梯度提升树得到的酶切组合重要性。

[0155] 在本发明中，确定酶切组合重要性的方法如下：

[0156] 首先，分析各变量相关性的方法如下：

[0157] 为了分析各变量的相关关系及组合影响情况，将相关系数转化为距离度量，采用的方法为d 1 |r|，其中d为度量距离，r为相关系数；r的公式如下：

[0158]

[0159] 公式中Cov(X,Y)为X，Y的协方差，D(X)、D(Y)分别为X、Y的方差。)。通过多尺度自主重抽样，可以得到每份数据进行层次聚类的p值，以此来评估层次聚类的不确定性。我们采用这种方式对51个变量(叶面积、叶长度、叶宽、叶周长、叶长宽比、叶形状因子、净光合速率、气孔导度、CO2浓度和水分利用率以及41个酶切组合的DNA甲基化水平)进行系统聚类，得到了各变量间的相关关系，相关性较强的变量用红色方框标出。

[0160] 其次，根据梯度提升树得到的酶切组合重要性，具体方法如下：

[0161] 利用R语言中的gbm包训练梯度提升树，

[0162] 参数的设置如下：

[0163] distribution＝'gaussian'，

[0164] n.trees＝10000,

[0165] shrinkage＝0.01,

[0166] interaction.depth＝5,

[0167] bag.fraction＝0.5,

[0168] cv.folds＝10。

[0169] Gbm包根据训练模型计算出各个输入变量对相应变量的重要性，得分越高，对响应变量影响越大。这样就得到每一组聚类样本的重要酶切片段组合，我们选择重要性得分在2以上的酶切片段组合。

[0170] 在本发明中，利用R语言中的gbm包中绘制边际效应图的函数并设置参数绘制，训练梯度提升树，得到重要影响的酶切组合后，我们绘制边际效应图。

[0171] 所述绘制过程中各个参数的设置如下：

[0172] distribution＝'gaussian'，

[0173] n.trees＝10000,

[0174] shrinkage＝0.01,

[0175] interaction.depth＝5,

[0176] bag.fraction＝0.5,

[0177] cv.folds＝10。

[0178] 根据所得到的边际效应图，就可以分析每一类小叶杨中，DNA甲基化水平叶片面积，净光合速率叶形状因子的边际效应，从而分析DNA甲基化对小叶杨叶片性状及光合特性的影响最后，在得到重要酶切片段组合后，本发明以重要酶切片段组合的DNA甲基化水平作为回归变量，利用高斯径向基函数建立LS-SVM回归预测模型，得到如式(9)所示的叶表型特征模型；

[0179]

[0180] 在本发明中，在SVM中，假设样本训练集为 回归函数为

[0181] f(x)w(x)b                  ⑷

[0182] 其中，w，b为回归参数，(x)为特征映射。

[0183] 而在LS-SVM中给问题转化为求解：

[0184]

[0185] s.t.y(x)w bIξ                ⑸

[0186] 其中，为正则化参数，ξ为松弛因子，In1(1,1,...,1)；

[0187] 构造Lagrange函数如下：

[0188]

[0189] 其中α为拉格朗日乘子；下求L(w,b,ξ,α)的鞍点，即最优点；

[0190]

[0191] 消去式(7)中的w，ξ，可得：

[0192]

[0193] 式(8)中，Ω(xixj),i,j 1,2,...,n，E为n n的单位矩阵。

[0194] 通过解式(8)得到α和b，则最小二乘支持向量机的估计函数为：

[0195]

[0196] 其中，k(xi,xj)为核函数，选取Gauss径向基(RBF)核函数 在进行回归前已对输入变量进行了标准化，并通过参数寻优，令10，1。

[0197] 本发明提供了基于DNA甲基化水平的木本植物叶片表型特征和光合特性预测模型的预测方法，用重要酶切片段组合的DNA甲基化水平输入所述方法构建的叶片表型特征和光合特性模型中预测叶形状因子、叶面积和叶净光合速率；

[0198] 所述重要酶切片段组合的DNA甲基化水平为所述构建方法中得到的重要酶切片段组合按照全基因组的甲基化水平计算公式计算得到的DNA甲基化水平。

[0199] 在本发明中，酶切组合的DNA甲基化水平对叶型因子，叶面积和净光合速率的影响程度绘制边际效应图。

[0200] 所述边际效应图的绘制方法为将重要酶切片段组合的DNA甲基化水平作为输入变量，然后调用gbm中自带的绘制边际效应图的函数绘制较重要的酶切片段对每个叶片的叶形状因子、叶面积和叶净光合速率的边际效应图。通过做出边际效用图，我们可以得到酶切组合的DNA甲基化水平对叶型因子，叶面积和净光合速率的影响。可以直接找到影响较大的酶切片段，进一步研究其上的基因位点，减少时间与经济成本。

[0201] 下面结合实施例对本发明提供的一种基于DNA甲基化水平预测木本植物的叶片形状和光合特性模型的构建方法和预测方法进行详细的说明，但是不能把它们理解为对本发明保护范围的限定。

[0202] 实施例1

[0203] 下面以在全国采集到的235例小叶杨个体为例。

[0204] 实验数据的采集：

[0205] 采用Dong Ci等人发表的Variation in genomic methylation in natural populations of Populus simonii  is associated with leaf shape and photosynthetic traits(Journal of Experimental Botany,Vol.67,No.3pp.723–737,
2016)中公开的酶切片段甲基化数据计算得到235个小叶杨样本个体DNA基因组甲基化水平，结果如表1所示(CC表示Chicheng County：赤城县,ZJK:张家口FX:陕西富县；LY:Linyou County麟游县；LX:Langao County岚皋县,LC:Luochuan County洛川县,GQ:Gaoling Count高陵县；HZ:Huzhu County互助县,XH:Xinghai County兴海县；W:Dulan County都兰县,MY:
Menyuan County门源县,SX:Song County嵩县,YC:Yichuan County伊川县,JL:Zhongning County中宁县,NM:Baotou City包头市,NW:Ningwu County宁武县,TRT:Taoranting Park北京陶然亭公园)。同时采用便携式激光叶片面积测量装置(CI-202)测定叶片表型，叶宽(简称width)，周长(简称perim)，叶面积(简称area)，叶形状因子(简称fact)，叶长度(简称length)，叶长宽比(简称ratio)和平均DNA甲基化水平(简称Dmavg))；叶片的光合特性采用便携式气体交换系统(li-6400 xt licor)仪器测定。为了获得最大瞬时净光合速率，将光合光子通量密度(PPFD)设置为1600，CO2浓度设置为400。净光合速率、气孔导度、细胞间CO2浓度和水分利用效率(WUE)在净光合速率下均有记录。光合特性数据结果见表2(地名标注同上)。

[0206] 1.使用改进FCM算法对小叶杨样本进行分类。

[0207] 1.1基于随机森林选取能够体现地理位置差异的重要特征变量。

[0208] 图1为基于随机森林选取能够体现地理位置差异的重要特征变量。由图1可知，小叶杨样本的差异主要取决于七个变量，例如叶宽(简称width)，周长(简称perim)，叶面积(简称area)，叶形状因子(简称fact)，叶长度(简称length)，叶长宽比(简称ratio)和平均DNA甲基化水平(简称Dmavg)。

[0209] 1.2确定小叶杨样本最优聚类个数。图2为确定小叶杨样本最优聚类个数。从图2中可以看出，使用组内平方误差(SSE)评估时，选择斜率最小的聚类数为2，当使用26个聚类评估指标评估时，选择组内标准数(Number Criteria)最大的聚类数为2。因此，综合SSE与26个指标，推荐的最优聚类个数为2。

[0210] 1.3使用改进的FCM算法确定每类样本的个体。从图3中可以看出，第一类样本包含139例小叶杨个体，第二类样本包含97例小叶杨个体。

[0211] 2.预测模型回归变量的选择

[0212] 2.1变量相关性的研究

[0213] 由图4可以看出，发现有四组酶切组合的DNA甲基化水平相关性极强，比如H44E4和H47E4,H80E1和H80E14,H80E9和H86E9，并且第四组包含15个酶切组合(H65E7,H63E9,H86E7,H60E3,H60E1,H46E4,H31E1,H60E15,H65E8,H47E6,H63E11,H65E6,H34E3,H46E10,H46E11)。这些结果显示DNA甲基化在杨树基因组内是区域特异性的，并且强相关的酶切组合的DNA甲基化水平可能代表在这些区域上DNA甲基化的修饰过程是相似的。

[0214] 2.2基于机器学习算法(梯度提升树)所得的酶切组合重要性选择回归变量，由于该机器学习方法存在一定误差，所以每次学习后选出的重要酶切片段会略有不同。因此为了减小误差，本文通过多次学习后取出现较多次数的酶切片段，并按照重要性程度大小进行排序。

[0215] 由图5可知，第一类小叶杨样本中，酶切组合H31E1,H65E12,H80E14,H65E5,H80E2,H46E12,H80E12,H65E4和H80E1的DNA甲基化水平对叶型因子有较大影响，酶切组合H31E3,H46E11,H60E3,H63E10,H44E4,H80E2,H65E5,H63E11和H60E2的DNA甲基化水平对叶片面积影响较大，酶切组合H60E3,H31E3,H63E12,H46E4,H65E7,H86E16,H82E5,H80E1和H63E11的DNA甲基化水平对净光合速率影响较大。

[0216] 由图6可知，第二类小叶杨样本中，酶切组合H60E15,H60E1,H63E10,H65E6,H34E5,H65E7,H80E13,H82E5和H60E2的DNA甲基化水平对叶型因子有较大影响，酶切组合H80E13,H60E2,H82E5,H60E15,H63E11,H86E7,H80E9,H86E16,H65E6和H44E4的DNA甲基化水平对叶片面积影响较大，酶切组合H65E8,H60E1,H80E2,H86E7,H46E11,H44E15,H63E12,H82E5和H46E4的DNA甲基化水平对净光合速率影响较大。

[0217] 3.以筛选出的酶切组合作为回归变量，基于LS-SVM预测小叶杨表型特征和光合特性(叶型因子，叶片面积，净光合速率)的值，同时测定小叶杨表型特征和光合特性数据，结果见表3～表8。

[0218] 第一类小叶杨样本中：对叶型因子，叶片面积，净光合速率的预测准确率分别为96.26％,94.42％和96.88％。第二类小叶杨样本中，对叶型因子，叶片面积，净光合速率的预测准确率分别达到81.27％,92.1％和95.8％。

[0219] 实施例2

[0220] DNA甲基化水平与小叶杨叶片表型特征、光合特性的关系。

[0221] 为了探讨DNA甲基化水平对表型性状(叶片形状因子、叶片面积和净光合速率)的影响，我们分析了边际效用的数值。

[0222] 在第一类小叶杨样本中，酶切组合H31E3,H60E3,H63E10,H44E4,H80E2,H65E5,H63E11和H60E2的DNA甲基化水平对叶片面积的边际效应明显(如图7-1)。叶片面积的会随着酶切组合H60E3,H44E4和H60E2的DNA甲基化水平的升高而降低，随着H31E3,H63E10,H80E2,H65E5和H63E11的DNA甲基化水平的提高而提高。第二类小叶杨样本中，酶切组合H80E13,H60E2,H82E5,H60E15,H63E11,H86E7,H80E9,H86E16,H65E6和H44E4的DNA甲基化水平对叶片面积的边际效应十分明显(参考图7-2),并且，除了酶切组合H80E13,H60E15,H63E11和H86E7，其他酶切组合DNA甲基化水平的挺高会使叶片面积变小。

[0223] 对于净光合速率来说(如图8所示)，第一类小叶杨样本中，酶切组合H60E3,H31E3,H63E12,H46E4,H65E7,H86E16,H82E5,H80E1和H63E11的DNA甲基化水平的边际效应明显。其中，净光合速率会随着H31E3,H63E12和H65E7的DNA甲基化水平的提高而降低，随着H60E3,H46E4,H86E16,H82E5,H80E1和H63E11的DNA甲基化水平的提高而增大。在第二类样本中，酶切组合H65E8,H80E2和H86E7的DNA甲基化水平的边际效应明显，净光合速率会随着酶切组合H65E8和H86E7的DNA甲基化水平的提高而降低，随着酶切组合H80E2的DNA甲基化水平的提高而增大。

[0224] 从图9中可以看出，在第一类小叶杨样本中，酶切组合H31E1,H65E12,H80E14,H65E5,H80E2,H46E12,H80E12,H65E4和H80E1的DNA甲基化水平对叶形因子的边际效应十分明显，酶切组合H31E1,H65E12,H80E14,H65E5,H46E12和H65E4的DNA甲基化水平越高，叶形因子越小，但是叶形因子会随酶切组合H80E2,H80E12和H80E1的DNA甲基化水平的提高二提高。第二类小叶杨样本中，酶切组合H82E5和H60E2的DNA甲基化水平越高，叶形因子越大。但叶形因子会随H60E15,H60E1,H63E10,H65E6,H34E5,H65E7和H80E13的DNA甲基化水平的提高而降低。

[0225] 通过分析两类小叶杨样本，发现DNA甲基化水平的提高可能会使叶形因子减小。杨属的种具有表型的多样性。而DNA甲基化对杨属可塑性的贡献比较大。此时叶型因子可作为表型多样性的重要参考因子。叶型因子的值越接近1，则叶片的形状越接近圆形。小叶杨的叶片主要有以下两种(如图10)。通过计算他们的叶型因子，我们发现对于第一种叶片形态，叶型因子较大，在0.698-0.7853。第二种叶片形态的小叶杨的叶型因子一般小于0.5，并且其叶型因子值越大，叶片曲率越小。通过长期的实验和研究，我们发现第二种类型的小叶杨抗逆性较强。因此我们推断得到DNA甲基化可能是这一类小叶杨的抗逆性增强的原因。

[0226] 本发明为小叶杨在DNA甲基化作用下的生长发育提供了理论依据。

[0227] 通过对小叶杨的两个亚种群的分析，我们发现DNA甲基化水平的提高可能会降低叶片形状因子。杨属的种具有表型的多样性。而DNA甲基化对杨属可塑性的贡献比较大。此时叶型因子可作为表型多样性的重要参考因子。叶型因子的值越接近1，则叶片的形状越接近圆形。小叶杨的叶片主要有以下两种。通过计算他们的叶型因子，我们发现对于第一种叶片形态，叶型因子较大，在0.698-0.7853左右。第二种叶片形态的小叶杨的叶型因子一般小于0.5.并且其叶型因子值越大，叶片曲率越小。通过长期的实验和研究，我们发现第二种类型的小叶杨抗逆性较强。因此我们推断得到DNA甲基化可能是这一类小叶杨的抗逆性增强的原因。

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

[0229]

[0230]

[0231] 表4第二组聚类叶片样本的净光合速率的测定值和真实值

[0232]

[0233] 表5第二组聚类叶片样本的叶形状因子的测定值和真实值

[0234]特征 Pn的测定值 预测值 特征 Pn的测定值 预测值
CC27 0.09 0.10210262 FX63 0.65 0.61308217
CC40 0.15 0.15696717 FX7 0.1 0.11401913
FX104 0.05 0.0656351 FX76 0.6 0.56567484

[0235] 表6第一组聚类叶片样本的叶形状因子的测定值和真实值

[0236]特征 测定值 预测值 特征 测定值 预测值
CC10 0.63 0.6242989 FX61 0.59 0.5846074
CC11 0.61 0.6045739 FX64 0.7 0.6869076
CC12 0.6 0.5941265 FX65 0.64 0.636514

[0237] 表7第一组聚类叶片样本的叶片面积的测定值和真实值

[0238]特征 测定值 预测值 特征 测定值 预测值
CC10 26.98 27.45045 GQ54 10.78 12.71664
CC11 17.85 19.23766 GQ6 11.85 13.69531
CC12 15.07 16.25292 HZ19 18.23 19.44564

[0239] 表8第一组聚类叶片样本的净光合速率的测定值和真实值

[0240]特征 测定值 预测值 特征 测定值 预测值
CC10 18.1 17.869152 FX64 15.33 15.324247
CC11 17.07 16.932437 FX65 10.58 10.968846
CC12 17.03 17.026104 FX66 12.21 12.541189