CN104965997A - 一种基于植物功能与结构模型的作物虚拟育种方法 - Google Patents

Abstract

一种基于植物功能与结构模型的作物虚拟育种方法。本方法将结合生理和结构作物虚拟生长模型和数量遗传信息，通过模拟作物的生长发育过程，及繁殖过程，达到虚拟育种的目的，其步骤是：采集作物的生理数据、形态数据、数量基因数据和生长环境数据构建原始数据集；基于原始数据，通过计算机图形学的方法和作物生理生态学原理构建作物功能与结构模型；构建遗传模块，包含数量遗传信息，建立从数量遗传信息到相关模型参数的关联；加入对杂交过程中交叉、重组等遗传操作的模拟，实现对作物繁殖过程的模拟；通过对模型中个体的选择和杂交，实现虚拟育种；本发明的提出能够利用计算机技术实现作物的虚拟育种，从而为传统育种进程提供辅助和支持。

Description

一种基于植物功能与结构模型的作物虚拟育种方法

技术领域

本发明内容涉及植物生理生态学、植物可视化建模以及作物育种领域，是一种能够通过三维可视化模型的方式模拟作物育种过程的方法，可运用于对田间作物选育的辅助和支持。

背景技术

传统育种具有较为完善的育种理论和技术体系，而新技术的发展，包括生物技术、基因技术等，能促进了作物育种的成就与发展。目前新兴的辅助育种方式有分子标记辅助育种、细胞工程育种和基因工程育种等，将分子育种、基因育种等技术手段与常规育种方式结合起来，能够极大地促进育种水平的提高。然而，目前常规育种方式之外的其它新技术育种手段还都存在很大的局限，包括这方面研究需要大量技术积累，而且实施起来也需要特定的仪器和设备支持，这就造成了研究成本的增加，从而在发展过程中与实际育种过程脱离。而且即使是结合新技术的育种实验也需要较长周期的田间试验，尤其是在目标株型不确定的情况下缺乏针对性。另外，由于安全性评价和鉴定机制不够完善，相关新育种方式得到的粮食作物品种的安全性也受到非常大的质疑，这就让这些新技术的发展和推广受到限制。

作物功能与结构模型是一种比较新的可视化植物模型构建技术，它通过在计算机上以三维图形的方式展现作物的形态发育过程，并结合作物的生理过程和特定的环境因素对形态发育的调控作用，能够实现对作物的空间的结构在时间梯度上的变化。最近，植物模型研究者们开始将作物的数量遗传信息也加入到植物模型中，从而构建从作物数量基因、环境因素到主要生理进程，再到形态结构的发育动态这样完全的调控网络。

本发明主要植物功能与结构建模方法为基础，构建遗传模块，实现对作物繁殖过程中的遗传操作的模拟，并构建遗传信息到目标性状的关联，选择主要的遗传相关参数对生长动态进行调控影响，通过基因的差异性展现表现型的可塑性，最终模拟作物育种过程，实现作物的虚拟育种，可以结合现在分子育种等新的育种手段，对育种实践进行辅助。这种虚拟育种方法的提出尚属首次，具有较好的创新性和应用价值。

发明内容

本发明的目的是提供一种基于植物功能与结构建模的技术结合数量遗传信息模拟作物育种过程的方法，该方法可以较为真实地反映不同基因型个体的目标性状动态发育表现，从而通过计算机可视化的方式对模拟的作物个体进行虚拟选育，并快速得到不同标准的选育结果，最终为实际的田间育种过程提供辅助支持。

植物功能与结构建模的技术，是借助计算机可视化的手段，综合考虑植物生长的外因——气候因素、土壤因素、人工操作因素等，以及相关的内因——植物的主要生理进程、形态发育等，对植物的生长动态实现模拟的技术。通过在植物功能与结构模型中加入数量遗传信息，就能够构建“基因和环境-生理-形态”交互反馈网络，并以三维动态的方式进行展现；再加上对作物个体繁殖过程及其遗传操作的模拟，能够实现对育种过程的模拟。

本发明实现虚拟育种所采用的技术方案是：

一种基于功能与结构模型的作物虚拟育种方法运用包含生理功能、形态发育、数量遗传信息和特定环境因素的作物虚拟生长模型和人工选择操作结合的方法模拟作物的繁殖的育种过程，并包含以下步骤：

步骤1：数据采集

通过田间试验、分子试验、基因定位分析和文献搜索结合的方式，采集模型构建所需的数据，建立原始数据集，其中包含同一目标作物群体的各方面数据：作物生长环境数据、形态动态生长数据、主要生理进程数据、分子标记和遗传图谱数据、目标性状相关的数量基因数据等；

步骤2：功能与结构模型构建

首先，运用基于规则的植物建模方法，从计算机图形学技术出发，以作物生长原理为基础，采用扩展的XL(eXtended L-Systems)建模语言、Java编程语言以及实现图形替换的RGG(Relational Growth Grammar)语法规则，对作物器官形态进行重建、对生长过程进行模拟，建立作物结构模型(以水稻为例，即包含茎、叶、穗、籽粒等器官形态建成的可视化模型)；

其次，在结构模型的基础上，加入作物主要生理进程的模拟：运用LEAFC3光合作用模型，通过设定品种相关的光合参数，特定温度下的最大羧化作用速率V_m25(μmol m^-2s^-1)、光饱和时的潜在光合电子传递速率J_m25(μmol m^-2s^-1)、光合电子传输的活化能量E_j(J/mol)、CO₂的动力学参数K_c25(mol/mol)、O₂的动力学参数K_o25(mol/mol)等，以及品种无关的环境参数，比如空气温度T_a(℃)、空气中CO₂的浓度C_a(mol/mol)、相对湿度RH、叶片吸收的双向长短波辐射强度R_i(W/m²)、水平方向的风速Wspeed(m/s)等，模拟叶片内的光合相关的生化反应、气孔导度、叶片边界的物质和能量传输机制，运用Nikilov等人在1995年发表在Ecologial Modelling上(卷(期)：80(2-3)；页码：205-235)的算法计算植物叶片的CO₂短期稳态通量、水分及热量交换，从而模拟作物通过叶片光合作用产生同化物的速率，如下式所示：

A_n＝f(R_i,T_a,...,V_m25,J_m25,...)

这里记A_n为净光合速率(μmol m^-2s^-1)，则每个生长步长所积累的光合产量为：

$P_t=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(A_n\×a_i\×\mathrm{\Δ}{t}_d)$

这里，记P_t为光合产量(μmol)，a_i为第i张叶片的面积(m²)，n为作物植株个体叶片的数目，Δt_d为模型中一个生长步长(s)，即一天的时间长度，值为24×3600，另外，假设所有的同化物在分配到器官之前都汇集到一个同化物池中，记为AP(μmol)，则t时刻AP的增量即为：

ΔAP＝P_t-G_t

这里，记G_t为植物呼吸作用(包括生长呼吸和维持呼吸)消耗的同化物的量，即同化物池的瞬时增量ΔAP是时间t的光合产量P_t与消耗量的差；G_t由所有器官的生长量累加得来；特定器官的生长量由该器官生长函数结合源库模型计算得到：基于Yin等在2003年提出的beta生长函数，可以由下式计算：

$c_m=w_{\max }{\left(\frac{{2t}_e-t_m}{t_e(t_e-t_m)}\right)}^{\frac{t_m}{t_e-t_m}}$

记c_m为最大生长速率，t_m为生长速率最大的时刻，t_e为停止生长的时刻(即生长速率为0时)，也就是该器官的维度达到最大维度w_max(长度为m，面积为m²)的时刻(比如茎达到最大长度或者叶片达到最大面积)，进而可以计算在任意时刻t的特定器官的潜在生长速率gr_pot：

${\mathrm{gr}}_{\mathrm{pot}}=c_m\left(\frac{t_e-t}{t_e-t_m}\right){\left(\frac{t}{t_m}\right)}^{\frac{t_m}{t_e-t_m}}$

以水稻的株高为例，这里gr_pot即为t时刻株高的潜在生长速率，将所有器官的潜在生长速率累加，并乘以时间步长Δt_d，得到整个植株的库强度sd_tot：

sd_tot＝∑ss_potΔt_d

因此，可以根据特定器官的库强度占植株个体的库强度大小的比率，计算得到该器官生长步长内实际生长速率gr_real：

${\mathrm{gr}}_{\mathrm{real}}=\frac{{\mathrm{gr}}_{\mathrm{pot}}}{{\mathrm{sd}}_{\mathrm{tot}}}\mathrm{AP}({\mathrm{gr}}_{\mathrm{real}}\≤{\mathrm{gr}}_{\mathrm{pot}})$

这里的AP为当前时刻的同化物池大小；通过这些生长函数和分配函数的计算，能够实现整个植株在整个生长周期的发育动态展现；

接着，重点建立环境模型中光照模型，模拟太阳光源的位置和辐射度变化；虚拟的太阳光在天空中被分为直射光和散射光，包括他们在三维空间中的分布，以及在到达作物冠层后被冠层空间的遮挡等，从而实现叶片光合作用中光通量的大小；

最后，运用功能-结构建模技术，结合作物结构模型、生理模型和环境模型，以时间为轴线，通过基于XL语言和RGG的语法规则对作物器官的形成和生长做规则化解释，并实现规则和图形之间的迭代、替换，结合虚拟作物的同化物形成与分配，从而在结合生理功能、拓扑结构及光环境条件的基础上实现虚拟作物生长动态可视化，获得能够模拟作物生长过程的功能与结构模型；

步骤3：遗传模块的构建

为模型中个体增加遗传模块，该模块包含以数组形式存储的内在数量遗传属性，记为：

M＝{m₁,m₂,…,m_x,q₁,m_x+1,…,q_i,…,m_n}

D＝{d₁₂,d₂₃,…,d_xq1,…,d_(n-1)n}

这里，M表示其中一个染色体组上的分子标记基因型，其序列包括n个分子标记，m₁到m_n，其间根据特定位置信息放置了i个通过基因定位研究得到的数量性状基因座，q₁到q_i；序列每个位点取值为1或者-1，其中1表示该分子标记m或数量性状基因座q来自父本，-1表示来自母本，即父本中所有的分子标记基因型都为1，母本中都为-1；D表示相应的分子标记序列中相邻分子标记之间(或分子标记与数量位点之间)的遗传距离，d_(n-1)n表示第n-1分子标记与第n个分子标记之间的遗传距离(Morgan)，另有一个数组I_q存储数量性状基因座所处的位置信息，即在M中的下标：

I_q＝{I₁,I₂,…,I_i}

这里，从I₁到I_i的值表示第1个数量性状基因座到第i个基因座的下标；

通过上述两个数组的运算，能够得到数量性状基因座的基因型数组：

Q＝{x₁,x₂,…,x_i}

这里数组Q表示一个染色体组上数量基因的基因型，x_i即为上述M数组中的q_i，同时数量基因位点的效应值(用数组A表示)：

A＝{a₁,a₂,…,a_i}

a_i表示第i个位点上的数量性状基因座的加性效应值；另外根据实际情况，建立aa上位性效应数组(AA表示)和上位性基因型数组(Q_aa表示)，即一个位点的加性效应和另外一个位点的加性效应的互作效应及其相关位点基因型：

AA＝{aa₁,aa₂,…,aa_j}

Q_aa＝{x_aa1,x_aa2,…,x_aaj}

aa_j表示第j个加加上位性效应值，x_aaj表示第j个上位性效应位点的基因型；根据这些基因型信息和效应值信息，再加上目标性状的群体平均值μ，能够计算出特定品系个体的目标性状表现型值y：

$G\left(L_k\right)=\underset{i}{\overset{n}{\mathrm{\Σ}}}{x}_{i\left(L_k\right)}{a}_i+\underset{i}{\overset{n}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{n}{\mathrm{\Σ}}}{x}_{i\left(L_k\right)}{x}_{j\left(L_k\right)}{\mathrm{aa}}_{\mathrm{ij}}$

y(L_k)＝μ+G(L_k)

这里y(L_k)表示第k个品系个体的目标性状表现型值，由群体均值μ与该品系个体遗传效应值G(L_k)加和而得；而遗传效应值是所有位点加性效应a_i与该位点基因型x_i(Lk)乘积的累加，再加上所有上位性位点基因型x_i(Lk)x_j(Lk)与上位性效应值aa_ij的乘积；

这里计算得到的目标性状表现型值y(L_k)作为遗传相关的参数，代入到步骤2中的生长函数中，代替w_max的值，即同一性状生长函数中的最大生长维度w_max为遗传相关，取值由其个体基因型值和数量位点效应值决定；

以水稻(二倍体)的株高性状为例，结合两个染色体组上的基因型数据和加性及上位性效应值数据，能够得到株高的表现型值；将其用在主茎生长函数中作为最大长度w_max的值，从而控制该个体主茎在各个生长时刻的同化物竞争力和生长速率，并实现基因型差异到表现型差异的关联；

步骤4：模拟繁殖过程

首先，建立分子标记交叉(crossing over)算法：

Step 1：若还有未遍历的分子标记，则以未遍历的分子标记为入口，计算当前相邻分子标记(包括相邻的分子标记与数量基因)之间的交换率，由下式计算：

$r=\frac{1}{2}(1-e^{-2x})$

这里r是交换率，x是步骤3中遗传距离数组D中对应的分子标记之间的距离值；若全部分子标记都已经遍历，则执行Step 4；

Step 2：以计算出的交换率为概率条件，执行交换算法，若进入交换算法语句，则执行Step 3；若未进入交换算法语句，则执行Step 1；

Step 3：当前分子标记之后的所有分子标记在两个染色体组中进行交换；执行Step 1；

Step 4：执行结束，得到对应两组染色体组的目标分子标记基因型；

然后，主要运用染色体分离和重组(Recombination)算法模拟繁殖过程：

Step 1：父母本个体的染色体(以二倍体为例，记父本的两个染色体组为M₁₁和M₁₂；母本的记为M₂₁和M₂₂)，分别执行分子标记交叉算法，进行染色体组分离，得到交换后的的单倍染色体组M^’ ₁₁、M^’ ₁₂、M^’ ₂₁、M^’ ₂₂；

Step 2：以50％的概率取父本分离出的两组染色体组M^’ ₁₁、M^’ ₁₂中的一组，再同样以50％的概率取母本分离出的两组染色体组M^’ ₂₁、M^’ ₂₂中的一组，重组为子代一个个体染色体组的构成；

Step 3：通过生成的子代染色体组基因型，运用步骤3里的方法计算得到目标性状参数的取值，并应用的生长函数；

Step 4：如未达到所需的子代个体数目要求，则重新执行Step 1；如已达到子代个体数目要求，则本次繁殖结束，可以进行生长动态的模拟或者下一次的繁殖；

步骤5：进行虚拟选择育种

以选育DH群体为例，可以将运用作物功能与结构模型进行虚拟育种分为以下的几步：

Step 1：初始化模型群体，包括设置目标性状相关的数量遗传信息；

Step 2：在模型中群体生长可视化模拟的任意阶段，选择两个个体作为杂交的亲本；

Step 3：对选择的亲本执行步骤4所述的遗传操作，模拟繁殖过程；

Step 4：得到子1代群体：F1代，执行生长模拟；

Step 5：若要继续选育，执行Step 2；若选育完毕，则执行Step 6；

Step 6：群体染色体组分离，获得单倍体，并进行模拟加倍；

Step 7：得到目标群体：DH群体，执行生长模拟；

Step 8：通过模型可以得到目标群体的形态数据、生理数据、遗传数据等输出，虚拟选育过程结束；

除了DH群体的虚拟选育之外，上述步骤稍加修改即能实现其他群体的选育，比如重组自交系RIL群体等，总体过程是生长模拟→选择→繁殖→子代的生长模拟→继续选择→继续繁殖→新的子代……如此循环，按照特定的标准和目标来进行选择，就能得到模拟的后代群体，从而实现一般意义的虚拟育种。

本发明的有益效果是：

1)结合可视化模型的方法，首次实现对特定作物虚拟育种过程的模拟；

2)模型包含环境、生理、结构、遗传相关的数据，构建相互调控网络，比运用作物生理模型的方法更为系统，考虑因素更为全面；

3)通过有限的田间试验构建虚拟育种模型，能够在任意的生长时间以任意的育种策略对作物进行虚拟育种，得到模拟的育种结果，对不同的育种方式、育种标准进行比较，从而对育种实践进行辅助支持；

4)可以对育种后代的表现进行预测，节约大量的时间、人力和物力成本。

附图说明

图1作物虚拟育种模块图

图2以DH群体为例的虚拟育种流程图

具体实施方式

本发明构建作物功能与结构模型，并结合数量遗传信息和繁殖过程的模拟，能够实现作物的虚拟育种。结合附图，本方法的具体实施方式如下：

一种基于功能与结构模型的作物虚拟育种方法，运用包含生理功能、形态发育、数量遗传信息和特定环境因素的作物虚拟生长模型和人工选择操作结合的方法模拟作物的繁殖的育种过程，并包含以下步骤：

步骤1：数据采集

通过田间试验、分子试验、基因定位分析和文献搜索结合的方式，采集模型构建所需的数据，建立原始数据集，其中包含同一目标作物群体的各方面数据：作物生长环境数据、形态动态生长数据、主要生理进程数据、分子标记和遗传图谱数据、目标性状相关的数量基因数据等；

步骤2：功能与结构模型构建

首先，运用基于规则的植物建模方法，从计算机图形学技术出发，以作物生长原理为基础，采用扩展的XL(eXtended L-Systems)建模语言、Java编程语言以及实现图形替换的RGG(Relational Growth Grammar)语法规则，对作物器官形态进行重建、对生长过程进行模拟，建立作物结构模型(以水稻为例，即包含茎、叶、穗、籽粒等器官形态建成的可视化模型)；

其次，在结构模型的基础上，加入作物主要生理进程的模拟：运用LEAFC3光合作用模型，通过设定品种相关的光合参数，如25℃下最大羧化作用速率V_m25(μmol m^-2s^-1)、25℃下光饱和时的潜在光合电子传递速率J_m25(μmol m^-2s^-1)、光合电子传输的活化能量E_j(J/mol)、25℃下CO₂的动力学参数K_c25(mol/mol)、25℃下O₂的动力学参数K_o25(mol/mol)等，以及品种无关的环境参数，比如空气温度T_a(℃)、空气中CO₂的浓度C_a(mol/mol)、相对湿度RH、叶片吸收的双向长短波辐射强度R_i(W/m²)、水平方向的风速Wspeed(m/s)等，模拟叶片内的光合相关的生化反应、气孔导度、叶片边界的物质和能量传输机制，运用Nikilov等人在1995年发表在文献中的算法计算植物叶片的CO₂短期稳态通量、水分及热量交换，从而模拟作物通过叶片光合作用产生同化物的速率，如下式所示：

A_n＝f(R_i,T_a,...,V_m25,J_m25,...)

这里记A_n为净光合速率(μmol m^-2s^-1)，则每个生长步长所积累的光合产量为：

$P_t=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(A_n\×a_i\×\mathrm{\Δ}{t}_d)$

这里，记P_t为光合产量(μmol)，a_i为第i张叶片的面积(m²)，n为作物植株个体叶片的数目，Δt_d为模型中一个生长步长(s)，即一天的时间长度，值为24×3600，另外，假设所有的同化物在分配到器官之前都汇集到一个同化物池中，记为AP(μmol)，则t时刻AP的增量即为：

ΔAP＝P_t-G_t

这里，记G_t为植物呼吸作用(包括生长呼吸和维持呼吸)消耗的同化物的量，即同化物池的瞬时增量ΔAP是时间t的光合产量P_t与消耗量的差；G_t由所有器官的生长量累加得来；特定器官的生长量由该器官生长函数结合源库模型计算得到：基于Yin等在2003年提出的beta生长函数，可以由下式计算：

$c_m=w_{\max }{\left(\frac{{2t}_e-t_m}{t_e(t_e-t_m)}\right)}^{\frac{t_m}{t_e-t_m}}$

记c_m为最大生长速率，t_m为生长速率最大的时刻，t_e为停止生长的时刻(即生长速率为0时)，也就是该器官的维度达到最大维度w_max(长度为m，面积为m²)的时刻(比如茎达到最大长度或者叶片达到最大面积)，进而可以计算在任意时刻t的特定器官的潜在生长速率gr_pot：

${\mathrm{gr}}_{\mathrm{pot}}=c_m\left(\frac{t_e-t}{t_e-t_m}\right){\left(\frac{t}{t_m}\right)}^{\frac{t_m}{t_e-t_m}}$

以水稻的株高为例，这里gr_pot即为t时刻株高的潜在生长速率，将所有器官的潜在生长速率累加，并乘以时间步长Δt_d，得到整个植株的库强度sd_tot：

sd_tot＝∑ss_potΔt_d

因此，可以根据特定器官的库强度占植株个体的库强度大小的比率，计算得到该器官生长步长内实际生长速率gr_real：

${\mathrm{gr}}_{\mathrm{real}}=\frac{{\mathrm{gr}}_{\mathrm{pot}}}{{\mathrm{sd}}_{\mathrm{tot}}}\mathrm{AP}({\mathrm{gr}}_{\mathrm{real}}\≤{\mathrm{gr}}_{\mathrm{pot}})$

这里的AP为当前时刻的同化物池大小；通过这些生长函数和分配函数的计算，能够实现整个植株在整个生长周期的发育动态展现；

接着，重点建立环境模型中光照模型，模拟太阳光源的位置和辐射度变化；虚拟的太阳光在天空中被分为直射光和散射光，包括他们在三维空间中的分布，以及在到达作物冠层后被冠层空间的遮挡等，从而实现叶片光合作用中光通量的大小；

最后，运用功能-结构建模技术，结合作物结构模型、生理模型和环境模型，以时间为轴线，通过基于XL语言和RGG的语法规则对作物器官的形成和生长做规则化解释，并实现规则和图形之间的迭代、替换，结合虚拟作物的同化物形成与分配，从而在结合生理功能、拓扑结构及光环境条件的基础上实现虚拟作物生长动态可视化，获得能够模拟作物生长过程的功能与结构模型；

步骤3：遗传模块的构建

为模型中个体增加遗传模块，该模块包含以数组形式存储的内在数量遗传属性，记为：

M＝{m₁,m₂,…,m_x,q₁,m_x+1,…,q_i,…,m_n}

D＝{d₁₂,d₂₃,…,d_xq1,…,d_(n-1)n}

这里，M表示其中一个染色体组上的分子标记基因型，其序列包括n个分子标记，m₁到m_n，其间根据特定位置信息放置了i个通过基因定位研究得到的数量性状基因座，q₁到q_i；序列每个位点取值为1或者-1，其中1表示该分子标记m或数量性状基因座q来自父本，-1表示来自母本，即父本中所有的分子标记基因型都为1，母本中都为-1；D表示相应的分子标记序列中相邻分子标记之间(或分子标记与数量位点之间)的遗传距离，d_(n-1)n表示第n-1分子标记与第n个分子标记之间的遗传距离(Morgan)，另有一个数组I_q存储数量性状基因座所处的位置信息，即在M中的下标：

I_q＝{I₁,I₂,…,I_i}

这里，从I₁到I_i的值表示第1个数量性状基因座到第i个基因座的下标；

通过上述两个数组的运算，能够得到数量性状基因座的基因型数组：

Q＝{x₁,x₂,…,x_i}

这里数组Q表示一个染色体组上数量基因的基因型，x_i即为上述M数组中的q_i，同时数量基因位点的效应值(用数组A表示)：

A＝{a₁,a₂,…,a_i}

a_i表示第i个位点上的数量性状基因座的加性效应值；另外根据实际情况，建立aa上位性效应数组(AA表示)和上位性基因型数组(Q_aa表示)，即一个位点的加性效应和另外一个位点的加性效应的互作效应及其相关位点基因型：

AA＝{aa₁,aa₂,…,aa_j}

Q_aa＝{x_aa1,x_aa2,…,x_aaj}

aa_j表示第j个加加上位性效应值，x_aaj表示第j个上位性效应位点的基因型；根据这些基因型信息和效应值信息，再加上目标性状的群体平均值μ，能够计算出特定品系个体的目标性状表现型值y：

$G\left(L_k\right)=\underset{i}{\overset{n}{\mathrm{\Σ}}}{x}_{i\left(L_k\right)}{a}_i+\underset{i}{\overset{n}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{n}{\mathrm{\Σ}}}{x}_{i\left(L_k\right)}{x}_{j\left(L_k\right)}{\mathrm{aa}}_{\mathrm{ij}}$

y(L_k)＝μ+G(L_k)

这里y(L_k)表示第k个品系个体的目标性状表现型值，由群体均值μ与该品系个体遗传效应值G(L_k)加和而得；而遗传效应值是所有位点加性效应a_i与该位点基因型x_i(Lk)乘积的累加，再加上所有上位性位点基因型x_i(Lk)x_j(Lk)与上位性效应值aa_ij的乘积；

这里计算得到的目标性状表现型值y(L_k)作为遗传相关的参数，代入到步骤2中的生长函数中，代替w_max的值，即同一性状生长函数中的最大生长维度w_max为遗传相关，取值由其个体基因型值和数量位点效应值决定；

以水稻(二倍体)的株高性状为例，结合两个染色体组上的基因型数据和加性及上位性效应值数据，能够得到株高的表现型值；将其用在主茎生长函数中作为最大长度w_max的值，从而控制该个体主茎在各个生长时刻的同化物竞争力和生长速率，并实现基因型差异到表现型差异的关联；

步骤4：模拟繁殖过程

首先，建立分子标记交叉(crossing over)算法：

Step 1：若还有未遍历的分子标记，则以未遍历的分子标记为入口，计算当前相邻分子标记(包括相邻的分子标记与数量基因)之间的交换率，由下式计算：

$r=\frac{1}{2}(1-e^{-2x})$

这里r是交换率，x是步骤3中遗传距离数组D中对应的分子标记之间的距离值；若全部分子标记都已经遍历，则执行Step 4；

Step 2：以计算出的交换率为概率条件，执行交换算法，若进入交换算法语句，则执行Step 3；若未进入交换算法语句，则执行Step 1；

Step 3：当前分子标记之后的所有分子标记在两个染色体组中进行交换；执行Step 1；

Step 4：执行结束，得到对应两组染色体组的目标分子标记基因型；

然后，主要运用染色体分离和重组(Recombination)算法模拟繁殖过程：

Step 1：父母本个体的染色体(以二倍体为例，记父本的两个染色体组为M₁₁和M₁₂；母本的记为M₂₁和M₂₂)，分别执行分子标记交叉算法，进行染色体组分离，得到交换后的的单倍染色体组M^’ ₁₁、M^’ ₁₂、M^’ ₂₁、M^’ ₂₂；

Step 2：以50％的概率取父本分离出的两组染色体组M^’ ₁₁、M^’ ₁₂中的一组，再同样以50％的概率取母本分离出的两组染色体组M^’ ₂₁、M^’ ₂₂中的一组，重组为子代一个个体染色体组的构成；

Step 3：通过生成的子代染色体组基因型，运用步骤3里的方法计算得到目标性状参数的取值，并应用的生长函数；

Step 4：如未达到所需的子代个体数目要求，则重新执行Step 1；如已达到子代个体数目要求，则本次繁殖结束，可以进行生长动态的模拟或者下一次的繁殖；

步骤5：进行虚拟选择育种

以选育DH群体为例，可以将运用作物功能与结构模型进行虚拟育种分为以下的几步：

Step 1：初始化模型群体，包括设置目标性状相关的数量遗传信息；

Step 2：在模型中群体生长可视化模拟的任意阶段，选择两个个体作为杂交的亲本；

Step 3：对选择的亲本执行步骤4所述的遗传操作，模拟繁殖过程；

Step 4：得到子1代群体：F1代，执行生长模拟；

Step 5：若要继续选育，执行Step 2；若选育完毕，则执行Step 6；

Step 6：群体染色体组分离，获得单倍体，并进行模拟加倍；

Step 7：得到目标群体：DH群体，执行生长模拟；

Step 8：通过模型可以得到目标群体的形态数据、生理数据、遗传数据等输出，虚拟选育过程结束；

除了DH群体的虚拟选育之外，上述步骤稍加修改即能实现其他群体的选育，比如重组自交系RIL群体等，总体过程是生长模拟→选择→繁殖→子代的生长模拟→继续选择→继续繁殖→新的子代……如此循环，按照特定的标准和目标来进行选择，就能得到模拟的后代群体，从而实现一般意义的虚拟育种。

Claims (6)

1.一种基于功能与结构模型的作物虚拟育种方法，运用包含生理功能、形态发育、数量遗传信息和特定环境因素的作物虚拟生长模型和人工选择操作结合的方法模拟作物的繁殖的育种过程，并包含以下步骤：

步骤1：数据采集

通过田间试验、分子试验、基因定位分析和文献搜索结合的方式，采集模型构建所需的数据，建立原始数据集，其中包含同一目标作物群体的各方面数据：作物生长环境数据、形态动态生长数据、主要生理进程数据、分子标记和遗传图谱数据、目标性状相关的数量基因数据；

步骤2：功能与结构模型构建

21.首先，运用基于规则的植物建模方法，从计算机图形学技术出发，以作物生长原理为基础，采用扩展的XL(eXtended L-Systems)建模语言、Java编程语言以及实现图形替换的RGG(Relational Growth Grammar)语法规则，对作物器官形态进行重建、对生长过程进行模拟，建立作物结构模型；

22.其次，在结构模型的基础上，加入作物主要生理进程的模拟：运用LEAFC3光合作用模型，通过设定品种相关的光合参数，特定温度下的最大羧化作用速率V_m25(μmol m^-2s^-1)、光饱和时的潜在光合电子传递速率J_m25(μmol m^-2s^-1)、光合电子传输的活化能量E_j(J/mol)、CO₂的动力学参数K_c25(mol/mol)、O₂的动力学参数K_o25(mol/mol)，以及品种无关的环境参数，比如空气温度T_a(℃)、空气中CO₂的浓度C_a(mol/mol)、相对湿度RH、叶片吸收的双向长短波辐射强度R_i(W/m²)、水平方向的风速Wspeed(m/s)等，模拟叶片内的光合相关的生化反应、气孔导度、叶片边界的物质和能量传输机制，计算植物叶片的CO₂短期稳态通量、水分及热量交换，从而模拟作物通过叶片光合作用产生 同化物的速率，如下式所示：

A_n＝f(R_i,T_a,...,V_m25,J_m25,...) 

这里记A_n为净光合速率(μmol m^-2s^-1)，则每个生长步长所积累的光合产量为：

这里，记P_t为光合产量(μmol)，a_i为第i张叶片的面积(m²)，n为作物植株个体叶片的数目，Δt_d为模型中一个生长步长(s)，即一天的时间长度，值为24×3600，另外，假设所有的同化物在分配到器官之前都汇集到一个同化物池中，记为AP(μmol)，则t时刻AP的增量即为：

ΔAP＝P_t-G_t

这里，记G_t为植物呼吸作用(包括生长呼吸和维持呼吸)消耗的同化物的量，即同化物池的瞬时增量ΔAP是时间t的光合产量P_t与消耗量的差；G_t由所有器官的生长量累加得来；特定器官的生长量由该器官生长函数结合源库模型计算得到：基于Yin等在2003年提出的beta生长函数，可以由下式计算：

记c_m为最大生长速率，t_m为生长速率最大的时刻，t_e为停止生长的时刻(即生长速率为0时)，也就是该器官的维度达到最大维度w_max(长度为m，面积为m²)的时刻(比如茎达到最大长度或者叶片达到最大面积)，进而可以计算在任意时刻t的特定器官的潜在生长速率gr_pot：

以水稻的株高为例，这里gr_pot即为t时刻株高的潜在生长速率，将所有器官的潜在生长速率累加，并乘以时间步长Δt_d，得到整个植株的库强度sd_tot：

sd_tot＝∑ss_potΔt_d

因此，可以根据特定器官的库强度占植株个体的库强度大小的比率，计算得到该器官生长步长内实际生长速率gr_real：

这里的AP为当前时刻的同化物池大小；通过这些生长函数和分配函数的计算，能够实现整个植株在整个生长周期的发育动态展现；

接着，重点建立环境模型中光照模型，模拟太阳光源的位置和辐射度变化；虚拟的太阳光在天空中被分为直射光和散射光，包括他们在三维空间中的分布，以及在到达作物冠层后被冠层空间的遮挡等，从而实现叶片光合作用中光通量的大小；

23.最后，运用功能-结构建模技术，结合作物结构模型、生理模型和环境模型，以时间为轴线，通过基于XL语言和RGG的语法规则对作物器官的形成和生长做规则化解释，并实现规则和图形之间的迭代、替换，结合虚拟作物的同化物形成与分配，从而在结合生理功能、拓扑结构及光环境条件的基础上实现虚拟作物生长动态可视化，获得能够模拟作物生长过程的功能与结构模型；

步骤3：遗传模块的构建

为模型中个体增加遗传模块，该模块包含以数组形式存储的内在数量遗传属性，记为：

M＝{m₁,m₂,…,m_x,q₁,m_x+1,…,q_i,…,m_n}

D＝{d₁₂,d₂₃,…,d_xq1,…,d_(n-1)n}

这里，M表示其中一个染色体组上的分子标记基因型，其序列包括n个分子标记，m₁到m_n，其间根据特定位置信息放置了i个通过基因定位研究得到的数量性状基因座，q₁到q_i；序列每个位点取值为1或者-1，其中1表示该分子标记m 或数量性状基因座q来自父本，-1表示来自母本，即父本中所有的分子标记基因型都为1，母本中都为-1；D表示相应的分子标记序列中相邻分子标记之间(或分子标记与数量位点之间)的遗传距离，d_(n-1)n表示第n-1分子标记与第n个分子标记之间的遗传距离(Morgan)，另有一个数组I_q存储数量性状基因座所处的位置信息，即在M中的下标：

I_q＝{I₁,I₂,…,I_i}

这里，从I₁到I_i的值表示第1个数量性状基因座到第i个基因座的下标；

通过上述两个数组的运算，能够得到数量性状基因座的基因型数组：

Q＝{x₁,x₂,…,x_i}

这里数组Q表示一个染色体组上数量基因的基因型，x_i即为上述M数组中的q_i，同时数量基因位点的效应值(用数组A表示)：

A＝{a₁,a₂,…,a_i}

a_i表示第i个位点上的数量性状基因座的加性效应值；另外根据实际情况，建立aa上位性效应数组(AA表示)和上位性基因型数组(Q_aa表示)，即一个位点的加性效应和另外一个位点的加性效应的互作效应及其相关位点基因型：

AA＝{aa₁,aa₂,…,aa_j}

Q_aa＝{x_aa1,x_aa2,…,x_aaj}

aa_j表示第j个加加上位性效应值，x_aaj表示第j个上位性效应位点的基因型；根据这些基因型信息和效应值信息，再加上目标性状的群体平均值μ，能够计算出特定品系个体的目标性状表现型值y：

y(L_k)＝μ+G(L_k)

这里y(L_k)表示第k个品系个体的目标性状表现型值，由群体均值μ与该品系个 体遗传效应值G(L_k)加和而得；而遗传效应值是所有位点加性效应a_i与该位点基因型x_i(Lk)乘积的累加，再加上所有上位性位点基因型x_i(Lk)x_j(Lk)与上位性效应值aa_ij的乘积；

这里计算得到的目标性状表现型值y(L_k)作为遗传相关的参数，代入到步骤2中的生长函数中，代替w_max的值，即同一性状生长函数中的最大生长维度w_max为遗传相关，取值由其个体基因型值和数量位点效应值决定；

以水稻(二倍体)的株高性状为例，结合两个染色体组上的基因型数据和加性及上位性效应值数据，能够得到株高的表现型值；将其用在主茎生长函数中作为最大长度w_max的值，从而控制该个体主茎在各个生长时刻的同化物竞争力和生长速率，并实现基因型差异到表现型差异的关联；

步骤4：模拟繁殖过程

41.首先，建立分子标记交叉(crossing over)算法：

Step 411：若还有未遍历的分子标记，则以未遍历的分子标记为入口，计算当前相邻分子标记(包括相邻的分子标记与数量基因)之间的交换率，由下式计算：

这里r是交换率，x是步骤3中遗传距离数组D中对应的分子标记之间的距离值；若全部分子标记都已经遍历，则执行Step 4；

Step 412：以计算出的交换率为概率条件，执行交换算法，若进入交换算法语句，则执行Step 3；若未进入交换算法语句，则执行Step 1；

Step 413：当前分子标记之后的所有分子标记在两个染色体组中进行交换；执行Step 1；

Step 414：执行结束，得到对应两组染色体组的目标分子标记基因型；

42.然后，主要运用染色体分离和重组(Recombination)算法模拟繁殖过程：

Step 421：父母本个体的染色体(以二倍体为例，记父本的两个染色体组为M₁₁和M₁₂；母本的记为M₂₁和M₂₂)，分别执行分子标记交叉算法，进行染色体组分离，得到交换后的的单倍染色体组M’₁₁、M’₁₂、M’₂₁、M’₂₂；

Step 422：以50％的概率取父本分离出的两组染色体组M’₁₁、M’₁₂中的一组，再同样以50％的概率取母本分离出的两组染色体组M’₂₁、M’₂₂中的一组，重组为子代一个个体染色体组的构成；

Step 423：通过生成的子代染色体组基因型，运用步骤3里的方法计算得到目标性状参数的取值，并应用的生长函数；

Step 424：如未达到所需的子代个体数目要求，则重新执行Step 1；如已达到子代个体数目要求，则本次繁殖结束，可以进行生长动态的模拟或者下一次的繁殖；

步骤5：进行虚拟选择育种

以选育DH群体为例，可以将运用作物功能与结构模型进行虚拟育种分为以下的几步：

Step 51：初始化模型群体，包括设置目标性状相关的数量遗传信息；

Step 52：在模型中群体生长可视化模拟的任意阶段，选择两个个体作为杂交的亲本；

Step 53：对选择的亲本执行步骤4所述的遗传操作，模拟繁殖过程；

Step 54：得到子1代群体：F1代，执行生长模拟；

Step 55：若要继续选育，执行Step 2；若选育完毕，则执行Step 6；

Step 56：群体染色体组分离，获得单倍体，并进行模拟加倍；

Step 57：得到目标群体：DH群体，执行生长模拟；

Step 58：通过模型可以得到目标群体的形态数据、生理数据、遗传数据等输出，虚拟选育过程结束；

除了DH群体的虚拟选育之外，上述步骤稍加修改即能实现其他群体的选育，比如重组自交系RIL群体等，总体过程是生长模拟→选择→繁殖→子代的生长模拟→继续选择→继续繁殖→新的子代……如此循环，按照特定的标准和目标来进行选择，就能得到模拟的后代群体，从而实现一般意义的虚拟育种。

2.如权利要求1所述的基于功能与结构模型的作物虚拟育种方法，其特征在于：所述步骤1中的原始数据集，包含同一目标作物群体的生长环境、形态发育数据、生理进程数据、数量遗传信息数据等，其中数量遗传信息数据包括遗传图谱数据、分子标记基因型数据和目标性状相关的数量性状基因座位置和效应值数据，不同类型数据之间相互关联，互有影响。

3.如权利要求1所述的基于功能与结构模型的作物虚拟育种方法，其特征在于：所述步骤2中将LEAFC3光合作用模型、beta生长函数、同化物池、潜在生长速率应用到作物功能与结构模型中，构建植物生理-形态-环境的互反馈网络，并通过同化物池关联源和库，利用潜在生长速率计算总体库强度和特定器官对同化物的竞争力，利用同化物池和库强度大小来分配光合产物，以器官的最终维度大小为关键参数。

4.如权利要求1所述的基于功能与结构模型的作物虚拟育种方法，其特征在于：所述步骤3中，运用内置属性的方法存储作物个体虚拟分子标记序列和其它相关的遗传信息，通过遗传模型和虚拟的数量基因型计算表现型值，将表现型值作为该器官的最终维度代入到生长函数中计算潜在生长速率，竞争同化物，实现基因-生理-形态-环境之间的调控网络。

5.如权利要求1所述的基于功能与结构模型的作物虚拟育种方法，其特征 在于：所述步骤4中，在模型中通过模拟分子标记序列作为等位基因分离时的交叉、重组操作，实现对作物模拟个体杂交繁殖过程的模拟，其中通过分子标记之间的遗传距离信息计算相邻标记之间的交换率，模拟分离后的染色体组，再通过模拟染色体重组构建子代个体的虚拟染色体，即分子标记序列。

6.如权利要求1所述的基于功能与结构模型的作物虚拟育种方法，其特征在于：所述步骤5中，通过模型操作的方式可以在作物生长的任意时期，对虚拟作物个体进行选择、杂交、生长、再选择、再杂交等，如此往复，根据特定的选择策略和选择标准，可以进行虚拟育种。