本发明公开了一种用于复杂性状遗传改良的基因组选择方法，包括以下步骤：(1)统计遗传模型的建立；(2)主效基因定位；(3)遗传参数估计；(4)基于个体全基因组育种值的选择。与现有技术相比，本发明的有益效果为：a)与GBLUP的模型假设相比，本申请提出的模型假设更符合生物学遗传规律，具有更高的预测力，即具有更高的基因组选择精度；b)本发明基于混合线性模型，具有很大的灵活性，便于模型扩展。
1.一种用于复杂性状遗传改良的基因组选择方法，其特征在于，包括以下步骤：(1)统计遗传模型的建立：
将基因分为主效基因和微效基因，建立全基因组主微基因全模型用于表型预测，所述全基因组主微基因全模型为简单加性遗传模型，假设群体大小为n，分子标记数为m，主效基因位点数为c，则所述全基因组主微基因全模型的矩阵形式可以表示为：其中，y为n×1的性状表型值向量，μ为群体均值向量；ai为第i个主效基因的加性效应，作固定效应；uk为第k个微效基因的加性效应，作随机效应且满足均值为零、方差为 的正态分布 xi和zk分别为ai和uk的系数向量；b和u分别为固定效应和随机效应向量，X和Z分别为对应的系数矩阵；ε为残差效应向量，作随机效应且服从均值为零、方差为的多元正态分布 其中In表示大小为n×n的单位矩阵，ε与u相互独立；
首先，对全基因组分子标记进行关联分析，对于第j个分子标记有：y＝μ+xjaj+ε，
其中，y为n×1的性状表型值向量，μ为群体均值向量；aj为该分子标记的效应值，作固定效应；xj为该分子标记对应的系数向量；ε为残差效应向量，作随机效应且服从多元正态分布其次，利用Henderson III方法在原假设H0：aj＝0下依次对全基因组分子标记位点显著性进行检验，当某一分子标记位点的F统计量值超过阈值时，则认为该标记效应显著，并将该位点作为主效基因选入所述全基因组主微基因全模型中；
在确定了主效基因位点后，对所述全基因组主微基因全模型进行参数估计，并估算个体全基因组育种值；
在得到所述全基因组主微基因全模型的参数估计值之后，依据个体全基因组育种值大小排序进行基因组选择。
2.如权利要求1所述的基因组选择方法，其特征在于，步骤(1)中，所述全基因组主微基因全模型为简单加性遗传模型时，y的期望值和方差‑协方差矩阵分别为：E(y)＝Xb，
其中，亲缘关系矩阵Ka＝ZZ ；如果用 表示微效基因总效应，此时所述全基因组主微基因全模型可以改写为：
3.如权利要求2所述的基因组选择方法，其特征在于，将个体全基因组育种值的估算分为两个部分，一个是作为固定效应的主效基因效应，用 估算；另一个是作为随机效应的微效基因效应，直接从整体角度得到其总体估计值，即微效基因总效应估计值 计算公式为：其中， 此时，只需要估算b、 和 并根据已知系数矩阵Z和X，即可求得并可由此计算个体全基因组育种值的估计值
4.如权利要求3所述的基因组选择方法，其特征在于，步骤(3)遗传参数估计时，利用最小范数二阶无偏估计法对随机效应的方差进行估算，得到方差 和 的估计值 和 再利用广义最小二乘法对固定效应的效应值进行估算，得到主效基因效应向量b的估计值
5.如权利要求2所述的基因组选择方法，其特征在于，步骤(1)中，所述全基因组主微基因全模型在简单加性遗传模型的基础上拓展到显性和上位性模型，拓展后为：其中，假设共定位到c个加性或显性主效基因位点以及e对上位性主效基因位点，ai和di分别表示主效基因位点i的加性和显性效应值，(aa)jj′、(dd)jj′、(ad)jj′和(da)jj′分别表示主效基因位点j和主效基因位点j′之间的加加上位性、显显上位性、加显上位性和显加上位a d aa dd ad da性效应值，作固定效应；用ξ、ξ、ξ 、ξ 、ξ 和ξ 分别表示六种遗传效应的微效基因总效应a d向量，作随机效应；用符号⊙表示向量与向量之间的点乘；x 和x分别表示由加性和显性两种不同的基因型编码方式得到的系数向量，对于具有两个等位基因A和a的位点j，第i个个体在该位点的加性编码系数 和显性编码系数 分别为，a d aa dd ad da a d aa dd ad da用b、b、b 、b 、b 、b 分别表示主效基因的六种遗传效应向量，X 、X、X 、X 、X 、X为对应系数矩阵。
6.如权利要求1所述的基因组选择方法，其特征在于，步骤(2)中，阈值大小通过置换检验来确定。
7.如权利要求6所述的基因组选择方法，其特征在于，在一型错误率为0.05时，至少需要进行1000次置换检验，即随机打乱1000次个体基因组信息与表型值的对应关系，获得
1000个置换样本，在每一个置换样本中，对所有分子标记依次做F统计检验并取其中最大的F值，根据1000个F最大值构建概率密度分布图，取原假设下达到95％置信度的F值作为显著性阈值。
8.如权利要求1所述的基因组选择方法，其特征在于，使用拟合优度和预测力评价建立的统计遗传模型，拟合优度衡量统计遗传模型拟合的能力，用 表示，公式如下：预测力则衡量统计遗传模型预测的能力，用 表示，公式如下：其中，ERESS和PRESS分别表示估计的、预测的残差平方，SS为总平方和； 和 表示第i个个体的表型估计值和预测值，和 分别表示第i个个体的残差估计值和残差预测值；
或者，在计算PRESS时，通过HAT矩阵中的中心杠杆值对ERESS进行修正后间接计算PRESS，从而减少重复验证的运算量；HAT矩阵，又称投射矩阵，是表型真实值和表型估计值T ‑1 T之间的一个转化矩阵，在固定模型中， H＝X(X X) X就被称为该固定模型的HAT矩阵，而HAT矩阵对角线上第i个元素hii被称为中心化杠杆值，HAT法中PRESS的计算公式如下：M F R F R ‑1 F T ‑1 ‑1 T ‑1
在混合线性模型假设下，H＝H+H (1‑H)，其中H ＝K(K+λI) ，H＝X(XV X) XV ，由此可以快速评价基因组选择精度，R F M
其中，H、H和H 分别表示用随机效应、用固定效应以及同时用固定效应和随机效应来预测个体表型值时HAT矩阵的形式；
K表示亲缘关系矩阵，V表示表型的方差‑协方差矩阵，I表示单位矩阵，
9.如权利要求1所述的基因组选择方法，其特征在于，步骤(4)中个体全基因组育种值的估计从训练群体拓展到训练群体与测试群体，拓展后为：假设训练群体个体数为n1，测试群体个体数为n2，分子标记数为m，主效位点数为c，且假设这两个群体具有相似的遗传背景，则原模型y＝Xb+Zu+ε可改写为：其中y1(n1×1)和y2(n2×1)分别为训练群体和测试群体的表型值向量，b和u分别为固定效应向量和随机效应向量，随机效应向量u服从均值为零、方差为 的多元正态分布Im表示大小为m×m的单位矩阵；X1和X2分别为训练群体和测试群体固定效应的系数矩阵，Z1和Z2分别为训练群体和测试群体随机效应的系数矩阵；ε1和ε2分别为训练群体和测试群体的残差效应向量，作随机效应，且分别服从均值为零、方差为 和的多元正态分布 和 分别表示大小为n1×n1和n2×n2的单位矩阵，ε1和ε2与u均相互独立；对应的期望值和方差‑协方差矩阵分别为：依据训练群体和测试群体大小可以将加性亲缘关系矩阵Ka划分为，在已知y1、X1、X2、Z1和Z2的前提下，并通过参数估计获得b、 和 的估计值，即可求得训练群体个体全基因组育种值的预测值对预测所得测试群体个体全基因组育种值的排序作为测试群体育种选择的依据。
一种用于复杂性状遗传改良的基因组选择方法
技术领域
[0001] 本发明涉及计算生物学技术领域，特别是涉及一种用于复杂性状遗传改良的基因组选择方法。
背景技术
[0002] 大部分人类疾病与农艺性状都是受到多基因和非遗传因素共同影响的复杂性状，对这些复杂性状的精确解析和预测对提高疾病诊断率以及提高作物品质都具有十分重要的意义。为了对复杂性状进行更精准地预测，Meuwissen等人(Meuwissen THE,Hayes BJ,and Goddard ME.Prediction of total genetic value using genome‑wide dense marker maps.Genetics,157(4):1819‑1829.(2001))最早提出了基因组选择(Genomic Selection，简称GS)，一种基于全基因组信息预测性状表型值的统计方法。与分子标记辅助选择(Marker‑assisted Selection，简称MAS)不同，基因组选择的最终目的不只是找到某个或某部分与性状关联的关键基因(或称主效基因)，而是要利用全基因组的信息(包括微效基因)去解释复杂性状的遗传机制，进而达到性状预测的目的。基因组选择主要分为以下三个步骤：首先，利用全基因组分子标记信息及表型信息，通过对训练群体的连锁或关联分析，建立性状育种值的预测模型；其次，根据建立的预测模型及测试群体的全基因组分子标记信息，对测试群体的个体表型值进行估算；最后，根据估算所得个体表型值进行择优选择。
[0003] 具体应用到农业领域，训练群体和测试群体往往来自同一批种质资源，在测序技术相当发达的今天，对大量种质资源的全基因组测序成本已经大大下降。但与此同时，高通量的农艺性状测量技术尚未发展成熟。通过基因组选择，我们只需从大量种质资源中选取部分种子进行种植及性状测量，再根据预测模型和剩余种质资源的基因组信息进行育种值估算，选择最具育种潜力的个体。可见，基因组选择无需对整个选育群体进行费时费力地表型测量，无论是对纯系育种还是杂交育种，都能大大提高品种选育效率，加快性状选择进程(Crossa J,Perez‑Rodriguez P,Cuevas J,Montesinos‑Lopez O,Jarquin D,de los Campos G,et al.Genomic Selection in Plant Breeding:Methods,Models,and Perspectives.Trends in Plant Science,22(11):961‑975.(2017))。基因组选择对于加深育种家们对种质资源的认识，提高个体育种值的预测精度以及为作物育种提供精准指导具有十分重要的意义。
[0004] 然而，现有基于基因组选择的预测模型和策略缺乏对复杂性状主效基因和微效基因的区分和整合，一是指忽略了主效基因和微效基因选择方式与效率上的差异，二是指未将主效基因效应和微效基因效应共同纳入育种值计算。Wang等人(Wang D,El‑Basyoni IS,Baenziger PS,Crossa J,Eskridge KM,and Dweikat I.Prediction of genetic values of quantitative traits with epistatic effects in plant breeding populations.Heredity,109(5):313‑319.(2012))直接用估算的主效基因效应值来计算个体育种值，忽略了微效基因对育种值的贡献。Xu等人(Xu SH,Zhu D,and Zhang QF.Predicting hybrid performance in rice using genomic best linear unbiased prediction.Proceedings of the National Academy of Sciences of the United States of America,111(34):12456‑12461.(2014))在GBLUP模型中将所有分子标记作为随机效应放入模型，并假设它们服从同一正态分布，以随机效应整体估算值作为个体育种值，用于预测水稻杂交组合在产量性状上的表现，忽视了主效基因对育种值的贡献。
发明内容
[0005] 本申请针对现有技术中存在的上述不足，提供了一种用于复杂性状遗传改良的基因组选择方法。
[0006] 一种用于复杂性状遗传改良的基因组选择方法，包括以下步骤：
[0007] (1)统计遗传模型的建立：
[0008] 将基因分为主效基因和微效基因，建立全基因组主微基因全模型用于表型预测，所述全基因组主微基因全模型为简单加性遗传模型，假设群体大小为n，分子标记数为m，主效基因位点数为c，则所述全基因组主微基因全模型的矩阵形式可以表示为：
[0009]
[0010] 其中，y为n×1的性状表型值向量，μ为群体均值；ai为第i个主效基因的加性效应，作固定效应；uk为第k个微效基因的加性效应，作随机效应且满足均值为零、方差为 的正态分布 假设；xi和zk分别为ai和uk的系数向量；b和u分别为固定效应和随机效应向量，X和Z分别为对应的系数矩阵；ε为残差效应向量，作随机效应且服从均值为零、方差为的多元正态分布 其中In表示大小为n×n的单位矩阵，ε与u相互独立；
[0011] (2)主效基因定位：
[0012] 首先，对全基因组分子标记进行关联分析，对于第j个分子标记有：
[0013] y＝μ+xjaj+ε，
[0014] 其中，y为n×1的性状表型值向量，μ为群体均值向量；aj为该分子标记的效应值，作固定效应；xj为该分子标记对应的系数向量；ε为残差效应向量，作随机效应且服从多元正态分布
[0015] 其次，利用Henderson III方法在原假设H0：aj＝0下依次对全基因组分子标记位点显著性进行检验，当某一分子标记位点的F统计量值超过阈值时，则认为该标记效应显著，并将该位点作为主效基因选入所述全基因组主微基因全模型中；
[0016] (3)遗传参数估计：
[0017] 在确定了主效基因位点后，对所述全基因组主微基因全模型进行参数估计，并估算个体全基因组育种值；
[0018] (4)基于个体全基因组育种值的选择：
[0019] 在得到所述全基因组主微基因全模型的参数估计值之后，依据个体全基因组育种值大小排序进行基因组选择。
[0020] 优选的，步骤(1)中，所述全基因组主微基因全模型为简单加性遗传模型时，y的期望值和方差‑协方差矩阵分别为：
[0021] E(y)＝Xb，
[0022]T
[0023] 其中，亲缘关系矩阵Ka＝ZZ ；如果用 表示微效基因总效应，此时所述全基因组主微基因全模型可以改写为：
[0024] y＝Xb+ξ+ε。
[0025] 更优选的，将个体全基因组育种值的估算分为两个部分，一个是作为固定效应的主效基因效应，用 估算；另一个是作为随机效应的微效基因效应，直接从整体角度得到其总体估计值，即微效基因总效应估计值 计算公式为：
[0026]
[0027] 其中， 此时，只需要估算b、 和 并根据已知系数矩阵Z和X，即可求得 并可由此计算个体全基因组育种值的估计值
[0028] 进一步优选的，步骤(3)遗传参数估计时，利用最小范数二阶无偏估计法对随机效应的方差进行估算，得到方差 和 的估计值 和 再利用广义最小二乘法对固定效应的效应值进行估算，得到主效基因效应b的估计值
[0029] 优选的，步骤(1)中，所述全基因组主微基因全模型为在简单加性遗传模型的基础上拓展到显性和上位性模型，拓展后为：
[0030]
[0031] 其中，假设共定位到c个加性或显性主效基因位点以及e对上位性主效基因位点，ai和di分别表示主效基因位点i的加性和显性效应值，(aa)jj′、(dd)jj′、(ad)jj′和(da)jj′分别表示主效基因位点j和主效基因位点j′之间的加加上位性、显显上位性、加显上位性和显加a d aa dd ad da上位性效应值，作固定效应；用ξ、ξ、ξ 、ξ 、ξ 和ξ 分别表示六种遗传效应的微效基因总a d
效应向量，作随机效应；用符号⊙表示向量与向量之间的点乘；x 和x分别表示由加性和显性两种不同的基因型编码方式得到的系数向量，对于具有两个等位基因A和a的位点j，第i个个体在该位点的加性编码系数 和显性编码系数 分别为，
[0032]
[0033] 用ba、bd、baa、bdd、bad、bda分别表示主效基因的六种遗传效应向量，Xa、Xd、Xaa、Xdd、ad daX 、X 为对应系数矩阵。
[0034] 优选的，步骤(2)中，阈值大小通过置换检验来确定。
[0035] 更优选的，在一型错误率为0.05时，至少需要进行1000次置换检验，即随机打乱1000次个体基因组信息与表型值的对应关系，获得1000个置换样本，在每一个置换样本中，对所有分子标记依次做F统计检验并取其中最大的F值，根据1000个F最大值构建概率密度分布图，取原假设下达到95％置信度的F值作为显著性阈值。
[0036] 优选的，使用拟合优度和预测力评价建立的统计遗传模型，
[0037] 拟合优度衡量统计遗传模型拟合的能力，用 表示，公式如下：
[0038]
[0039] 预测力则衡量统计遗传模型预测的能力，用 表示，公式如下：
[0040]
[0041] 其中，ERESS和PRESS分别表示估计残差平方和以及预测残差平方和，SS为总平方和； 和 表示第i个个体的表型估计值和预测值，和 分别表示第i个个体的残差估计值和残差预测值；
[0042] 或者，在计算PRESS时，通过HAT矩阵中的中心杠杆值对ERESS进行修正后间接计算PRESS，从而减少重复验证的运算量；HAT矩阵，又称投射矩阵，是表型真实值和表型估计值T ‑1 T之间的一个转化矩阵，在固定模型中， H＝X(XX) X就被称为该
固定模型的HAT矩阵，而HAT矩阵对角线上第i个元素hii被称为中心化杠杆值，HAT法中PRESS的计算公式如下：
[0043]
[0044] 在混合线性模型假设下，HM＝HF+HR(1‑HF)，其中HR＝K(K+λI)‑1，HF＝X(XTV‑1X)‑1XTV‑1，由此可以快速评价基因组选择精度。
[0045] 优选的，步骤(4)中个体全基因组育种值的估计从训练群体拓展到训练群体与测试群体，拓展后为：
[0046] 假设训练群体个体数为n1，测试群体个体数为n2，分子标记数为m，主效位点数为c，且假设这两个群体具有相似的遗传背景，则原模型y＝Xb+Zu+ε可改写为：
[0047]
[0048] 其中y1(n1×1)和y2(n2×1)分别为训练群体和测试群体的表型值向量，b和u分别为固定效应向量和随机效应向量，随机效应向量u服从均值为零、方差为 的多元正态分布 Im表示大小为m×m的单位矩阵；X1和X2分别为训练群体和测试群体固定效应的系数矩阵，Z1和Z2分别为训练群体和测试群体随机效应的系数矩阵；ε1和ε2分别为训练群体和测试群体的残差效应向量，作随机效应，且分别服从均值为零、方差为 和的多元正态分布 和 分别表示大小为n1
×n1和n2×n2的单位矩阵，ε1和ε2与u均相互独立；对应的期望值和方差‑协方差矩阵分别为：
[0049]
[0050]
[0051] 依据训练群体和测试群体大小可以将加性亲缘关系矩阵Ka划分为，
[0052]
[0053] 在已知y1、X1、X2、Z1和Z2的前提下，并通过参数估计获得b、 和 的估计值，即可求得训练群体个体全基因组育种值的预测值
[0054]
[0055] 对预测所得测试群体个体全基因组育种值的排序可以作为测试群体育种选择的依据。
[0056] 与现有技术相比，本发明的有益效果为：
[0057] a)与GBLUP的模型假设相比，本申请提出的模型假设更符合生物学遗传规律，具有更高的预测力，即具有更高的基因组选择精度；
[0058] b)本发明基于混合线性模型，具有很大的灵活性，便于模型扩展。
附图说明
[0059] 图1为比较MMIBLUP与GBLUP在分析三种不同遗传结构假设下模拟数据的表现结果图。
具体实施方式
[0060] 在混合线性模型框架下，我们认为复杂性状中的主效基因往往是涉及性状发生最核心环节的基因，基因效应较大，更符合固定效应的假设；相应地，微效基因往往涉及性状发生的修饰环节，基因效应较小，更符合随机效应的假设。但无论是主效基因还是微效基因，都对表型值产生了或多或少的影响，因此有必要共同纳入全模型育种值的估算，以提高基因组选择的精度。
[0061] 由此，我们提出了一种用于解析复杂性状遗传结构的基因组选择新策略，即针对复杂性状遗传特点，在基因组选择全模型中既包含GWAS鉴别的主效基因作固定效应，又包含GBLUP整体估算的微效基因作随机效应，主微效应既有区分又有整合，共同纳入个体全基因组育种值的估算。该方法包括：
[0062] (1)统计遗传模型的建立
[0063] 根据我们对复杂性状遗传结构的假设，建立全基因组主微基因全模型(major‑minor Integrated Best Linear Unbiased Prediction，简称MMIBLUP)用于表型预测。以简单加性遗传模型为例，假设群体大小为n，分子标记数为m，主效基因位点数为c，则所述全基因组主微基因全模型的矩阵形式可以表示为：
[0064]
[0065] 其中，y为n×1的性状表型值向量，μ为群体均值；ai为第i个主效基因的加性效应，作固定效应；uk为第k个微效基因的加性效应，作随机效应且满足均值为零、方差为 的正态分布 假设；xi和zk分别为ai和uk的系数向量；b和u分别为固定效应和随机效应向量，X和Z分别为对应的系数矩阵；ε为残差效应向量，作随机效应且服从均值为零、方差为 的多元正态分布 其中In表示大小为n×n的单位矩阵，ε与u相互独
立；
[0066] 在该模型假设下，y的期望值和方差‑协方差矩阵分别为：
[0067] E(y)＝Xb(2)
[0068]T
[0069] 其中，亲缘关系矩阵Ka＝ZZ ；如果用 表示微效基因总效应，此时所述全基因组主微基因全模型可以改写为：
[0070] y＝Xb+ξ+ε(4)
[0071] (2)主效基因定位
[0072] 主效基因定位是新策略下全模型建立的前提条件。因此首先对全基因组分子标记进行关联分析，对于第j个分子标记有，
[0073] y＝μ+xjaj+ε(5)
[0074] aj为该分子标记的效应值，作固定效应，xj为该分子标记对应的系数向量，其他模型的定义与公式(1)中相同。利用Henderson III方法在原假设H0：aj＝0下依次对全基因组分子标记位点显著性进行检验，当某一分子标记位点的F统计量值超过阈值，则认为该标记效应显著，并将该位点作为主效基因选入全模型中，该方法已成功应用于复杂性状QTL定位分析(Yang J,Zhu J,and Williams RW.Mapping the genetic architecture of complex traits in experimental populations.Bioinformatics,23(12):1527‑1536.(2007))。
[0075] 阈值大小通过置换检验来确定(Churchill GA,and Doerge RW.Empirical Threshold Values for Quantitative Trait Mapping.Genetics,138(3):963‑971.(1994))。在一型错误率为0.05时，至少需要进行1000次置换检验，即随机打乱1000次个体基因组信息与表型值的对应关系，获得1000个置换样本。在每一个置换样本中，对所有分子标记依次做F统计检验并取其中最大的F值。根据1000个F最大值构建概率密度分布图，取原假设下达到95％置信度的F值作为显著性阈值。
[0076] 我们假设根据Henderson III法F统计量值以及置换检验确定的阈值在全基因组范围确定c个主效基因位点。
[0077] (3)遗传参数估计
[0078] 在确定了c个主效基因位点后，对(1)中全模型进行参数估计。利用最小范数二阶无偏估计(MINQUE)法对随机效应的方差进行估算，得到方差估计值 和 再利用广义最小二乘法(GLS)对固定效应的效应值进行估算，得到主效基因效应估计值
[0079] 在新模型假设下，我们认为个体全基因组育种值(Genomic Estimated Breeding Value,简称GEBV)的估算应当包括两个部分，一个是作为固定效应的主效基因效应，用估算；另一个是作为随机效应的微效基因效应，直接从整体角度得到其总体估计值，即微效基因总效应估计值 计算公式为：
[0080]
[0081] 其中， 此时，只需要估算b、 和 并根据已知系数矩阵Z和X，即可求得 并可由此计算个体全基因组育种值的估计值
[0082] (4)基因组选择精度的评价
[0083] 拟合优度(Goodness of fit)和预测力(Predictability)是基因组选择精度(Accuracy)的重要评价指标，也是模拟中用于比较不同模型好坏的常用指标。拟合优度衡量模型拟合的能力，用 表示，预测力则衡量模型预测的能力，用 表示，它们分别可以用公式(7)和(8)求得，
[0084]
[0085]
[0086] 其中，ERESS和PRESS分别表示估计残差平方和以及预测残差平方和，SS为总平方和。 和 表示第i个个体的表型估计值和预测值， 和 分别表示第i个个体的残差估计值和残差预测值。在计算PRESS时，除了可以通过交叉验证法(Cross Validation,简称CV)直接计算外，还可以通过HAT矩阵中的中心杠杆值对ERESS进行修正后间接计算PRESS，从而减少重复验证的运算量。HAT矩阵，又称投射矩阵，是表型真实值和表型估计值之间的一个转T ‑1 T化矩阵。在固定模型中， H＝X(XX) X就被称为该固定模型的HAT
矩阵。而HAT矩阵对角线上第i个元素hii被称为中心化杠杆值，HAT法中PRESS的计算公式如公式(9)所示。
[0087]
[0088] HAT法拓展到混合线性模型下时，需根据不同育种值预测策略选择HAT矩阵。这里，R F M用H、H和H分别表示用随机效应、用固定效应以及同时用固定效应和随机效应来预测个体表型值时HAT矩阵的形式(表1)。
[0089] 表1混合线性模型中HAT矩阵在不同预测方法下的形式
[0090]
[0091] (5)基于个体全基因组育种值的选择
[0092] 基因组选择的一个最终目的是能够根据训练群体(Training Population，简称TRN)的基因型信息和表型信息建立模型，再依据模型及测试群体(Testing Population，简称TST)的基因型信息对测试群体的表型信息进行可靠预测，从而为育种家选择杂交组合及选育优势后代提供参考，因此有必要将个体全基因组育种值的公式拓展到测试群体。
[0093] 假设训练群体个体数为n1，测试群体个体数为n2，分子标记数为m，主效位点数为c，且假设这两个群体具有相似的遗传背景，则原模型y＝Xb+Zu+ε可改写为，
[0094]
[0095] 其中y1(n1×1)和y2(n2×1)分别为训练群体和测试群体的表型值向量，b和u分别为固定效应向量和随机效应向量，随机效应向量u服从均值为零、方差为 的多元正态分布 Im表示大小为m×m的单位矩阵；X1和X2分别为训练群体和测试群体固定效应的系数矩阵，Z1和Z2分别为训练群体和测试群体随机效应的系数矩阵；ε1和ε2分别为训练群体和测试群体的残差效应向量，作随机效应，且分别服从均值为零、方差为和 的多元正态分布 和 分别表示大小为
n1×n1和n2×n2的单位矩阵，ε1和ε2与u均相互独立；对应的期望值和方差‑协方差矩阵分别为：
[0096]
[0097]
[0098] 依据训练群体和测试群体大小可以将加性亲缘关系矩阵Ka划分为，
[0099]
[0100] 在已知y1、X1、X2、Z1和Z2的前提下，并通过参数估计获得b、 和 的估计值，即可求得训练群体个体全基因组育种值的预测值
[0101]
[0102] 由此可在只知道测试群体基因型信息的前提下，依赖训练群体和测试群体之间的亲缘关系矩阵对测试群体个体全基因组育种值进行预测。将群体内个体按照全基因组育种值大小从高到低排列，为育种选择提供了有价值的参考意见。
[0103] (6)显性和上位性遗传模型的拓展
[0104] 对于复杂性状遗传结构的解析仅仅依靠加性遗传模型是不够的。因此，在加性遗传模型的基础上，可以将全基因组主微基因全模型拓展到显性和上位性模型，拓展后为：
[0105]
[0106] 其中，假设共定位到c个加性或显性主效基因位点以及e对上位性主效基因位点，ai和di分别表示主效基因位点i的加性和显性效应值，(aa)jj′、(dd)jj′、(ad)jj′和(da)jj′分别表示主效基因位点j和主效基因位点j′之间的加加上位性、显显上位性、加显上位性和显加a d aa dd ad da上位性效应值，作固定效应；用ξ、ξ、ξ 、ξ 、ξ 和ξ 分别表示六种遗传效应的微效基因总a d
效应向量，作随机效应；用符号⊙表示向量与向量之间的点乘；x 和x分别表示由加性和显性两种不同的基因型编码方式得到的系数向量，对于具有两个等位基因A和a的位点j，第i个个体在该位点的加性编码系数 和显性编码系数 分别为，
[0107]
[0108] 用ba、bd、baa、bdd、bad、bda分别表示主效基因的六种遗传效应向量，Xa、Xd、Xaa、Xdd、ad daX 、X 为对应系数矩阵。
[0109] (7)其他
[0110] 与现有技术相比，本发明的有益效果为：
[0111] a)与GBLUP的模型假设相比，我们提出的模型假设更符合生物学遗传规律，具有更高的预测力，即具有更高的基因组选择精度；
[0112] b)本发明基于混合线性模型，具有很大的灵活性，便于模型扩展。
[0113] 下面结合具体实例对本发明作进一步阐释。
[0114] (1)材料
[0115] 我们通过模拟来测试基因组选择新策略MMIBLUP的表现，并与已有基因组选择方法GBLUP进行比较。假设模拟群体大小为n，分子标记数为m，标记间相互独立，MAF的范围在0.05‑0.5，得到了n个个体m个位点的基因型。为了更充分比较两种方法在不同遗传结构下的表现，我们共模拟了三种不同遗传假设类型下的复杂性状，包括主效基因遗传(Scenario I)、主效基因与微效基因共同遗传(Scenario II和Scenario III)和微效基因遗传(Scenario IV)。
[0116] 在第一种主效基因遗传(Scenario I)假设下，认为全基因组中只有小部分位点对表型有影响。因此在m个位点中随机选取ml个位点作为主效基因位点，并假设其服从同一正2
态分布 遗传率h分别取0.1、0.2、0.5和0.8来模拟性状的不同遗传水平。根据遗传率公式 求得相应的残差方差 并根据正态分布 模拟得
到个体残差大小。假设主效基因位点方差 ml占m的1％。
[0117] 在第二种主效基因与微效基因共同遗传(Scenario II和ScenarioIII)假设下，认为全基因组中除了有小部分主效基因参与性状遗传之外，还有大部分的微效基因共同参与性状遗传。因此，除了假设随机选取的ml个主效基因位点服从正态分布 之外，还假设剩余ms个微效基因位点服从另一个正态分布 根据主效基因变异大小占全基因组变异(包括主效基因和微效基因变异)大小的比例，可以分为占比中等的Scenario II和占比较小的ScenarioIII两种情况。在Scenario II下，假设ml占m的1％，且 ms＝m‑ml，且 在Scenario III下，同样假设ml占m的1％，且 ms＝m‑ml，但
此时，两种情况主效基因变异大小占全基因组变异大小的比例
分别约为0.50和0.25。此时的遗传率计算公式为
2
分别求出h分别取0.1、0.2、0.5和0.8时对应的残差方差大小
[0118] 在第三种微效基因遗传(Scenario IV)假设下，认为全基因组所有位点均为微效基因位点，假设所有位点的加性效应大小均服从同一正态分布 这里我们假设加性效应方差 再根据设定的遗传率大小及遗传率公式 求得残2
差效应的方差大小，遗传率h同样分别取0.1、0.2、0.5和0.8。
[0119] 在四种情况下，模拟的主效基因变异占全基因组变异(包括主效基因和微效基因)的比例依次递减。每组模拟数据重复20次，群体大小n＝2000，分子标记数m＝3000，并对表型值进行标准化，对基因型系数矩阵按列进行标准化。
[0120] (2)软件
[0121] 在已有软件QTXNetwork(http://ibi.zju.edu.cn/software/)的基础上，利用C++语言编写QTXNetwork‑GS模块。在QTXNetwork‑GS模块中，可以实现对二倍体群体三种不同遗传假设下复杂性状主效基因和微效基因的模拟，主效基因定位及其效应估计，基因组选择全模型的拟合和预测，以及拟合和预测能力的统计学评价。
[0122] (3)结果
[0123] 我们比较了MMIBLUP与GBLUP在分析三种不同遗传结构假设下模拟数据的表现。如图1所见，基于不同复杂性状遗传假设，四张图分别对应四种不同类型的模拟性状，即主效基因遗传(Scenario I)、主效基因与微效基因共同遗传(Scenario II&III)以及微效基因遗传(Scenario IV)，其中Scenario II中主效基因遗传变异占全基因组遗传变异的比例约为0.5，Scenario III中主效基因遗传变异占全基因组遗传变异的比例约为0.25。每张图中横坐标为遗传率依次为0.1、0.2、0.5和0.8。每组模拟数据重复20次，群体大小n＝2000，分子标记数m＝3000。比较新模型MMIBLUP与GBLUP模型基因组选择精度(即预测力)水平的高低。模拟结果表明，当存在主效基因遗传变异时，不论是否有微效基因遗传变异存在，即Scenario I、II和III三种情形下，MMIBLUP的预测力均高于GBLUP；当遗传假设为微效基因遗传时，即Scenario IV情形下，GBLUP的预测力高于MMIBLUP。
[0124] (4)其他
[0125] 最后，还需要特别注意的是，以上所举例子仅是本发明的具体实施例子。显然，本发明不仅仅限于以上实施例子，还可以有许多变通的情况。本领域的技术人员从本发明公开的内容直接推导出或联想到的所有变通情况，均认为是本发明的保护范围。
