CN113223606A - 一种用于复杂性状遗传改良的基因组选择方法 - Google Patents

一种用于复杂性状遗传改良的基因组选择方法 Download PDF

Info

Publication number
CN113223606A
CN113223606A CN202110522399.XA CN202110522399A CN113223606A CN 113223606 A CN113223606 A CN 113223606A CN 202110522399 A CN202110522399 A CN 202110522399A CN 113223606 A CN113223606 A CN 113223606A
Authority
CN
China
Prior art keywords
effect
value
gene
genome
model
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110522399.XA
Other languages
English (en)
Other versions
CN113223606B (zh
Inventor
徐海明
张齐心
刘臣涛
朱天能
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN202110522399.XA priority Critical patent/CN113223606B/zh
Publication of CN113223606A publication Critical patent/CN113223606A/zh
Application granted granted Critical
Publication of CN113223606B publication Critical patent/CN113223606B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations

Landscapes

  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Molecular Biology (AREA)
  • Medical Informatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • Analytical Chemistry (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Genetics & Genomics (AREA)
  • Chemical & Material Sciences (AREA)
  • Physiology (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明公开了一种用于复杂性状遗传改良的基因组选择方法，包括以下步骤：(1)统计遗传模型的建立；(2)主效基因定位；(3)遗传参数估计；(4)基于个体全基因组育种值的选择。与现有技术相比，本发明的有益效果为：a)与GBLUP的模型假设相比，本申请提出的模型假设更符合生物学遗传规律，具有更高的预测力，即具有更高的基因组选择精度；b)本发明基于混合线性模型，具有很大的灵活性，便于模型扩展。

Description

一种用于复杂性状遗传改良的基因组选择方法
技术领域
本发明涉及计算生物学技术领域，特别是涉及一种用于复杂性状遗传改良的基因组选择方法。
背景技术
大部分人类疾病与农艺性状都是受到多基因和非遗传因素共同影响的复杂性状，对这些复杂性状的精确解析和预测对提高疾病诊断率以及提高作物品质都具有十分重要的意义。为了对复杂性状进行更精准地预测，Meuwissen等人(Meuwissen THE,Hayes BJ,and Goddard ME.Prediction of total genetic value using genome-wide densemarker maps.Genetics,157(4):1819-1829.(2001))最早提出了基因组选择(GenomicSelection，简称GS)，一种基于全基因组信息预测性状表型值的统计方法。与分子标记辅助选择(Marker-assisted Selection，简称MAS)不同，基因组选择的最终目的不只是找到某个或某部分与性状关联的关键基因(或称主效基因)，而是要利用全基因组的信息(包括微效基因)去解释复杂性状的遗传机制，进而达到性状预测的目的。基因组选择主要分为以下三个步骤：首先，利用全基因组分子标记信息及表型信息，通过对训练群体的连锁或关联分析，建立性状育种值的预测模型；其次，根据建立的预测模型及测试群体的全基因组分子标记信息，对测试群体的个体表型值进行估算；最后，根据估算所得个体表型值进行择优选择。
具体应用到农业领域，训练群体和测试群体往往来自同一批种质资源，在测序技术相当发达的今天，对大量种质资源的全基因组测序成本已经大大下降。但与此同时，高通量的农艺性状测量技术尚未发展成熟。通过基因组选择，我们只需从大量种质资源中选取部分种子进行种植及性状测量，再根据预测模型和剩余种质资源的基因组信息进行育种值估算，选择最具育种潜力的个体。可见，基因组选择无需对整个选育群体进行费时费力地表型测量，无论是对纯系育种还是杂交育种，都能大大提高品种选育效率，加快性状选择进程(Crossa J,Perez-Rodriguez P,Cuevas J,Montesinos-Lopez O,Jarquin D,de losCampos G,et al.Genomic Selection in Plant Breeding:Methods,Models,andPerspectives.Trends in Plant Science,22(11):961-975.(2017))。基因组选择对于加深育种家们对种质资源的认识，提高个体育种值的预测精度以及为作物育种提供精准指导具有十分重要的意义。
然而，现有基于基因组选择的预测模型和策略缺乏对复杂性状主效基因和微效基因的区分和整合，一是指忽略了主效基因和微效基因选择方式与效率上的差异，二是指未将主效基因效应和微效基因效应共同纳入育种值计算。Wang等人(Wang D,El-Basyoni IS,Baenziger PS,Crossa J,Eskridge KM,and Dweikat I.Prediction of genetic valuesof quantitative traits with epistatic effects in plant breedingpopulations.Heredity,109(5):313-319.(2012))直接用估算的主效基因效应值来计算个体育种值，忽略了微效基因对育种值的贡献。Xu等人(Xu SH,Zhu D,and ZhangQF.Predicting hybrid performance in rice using genomic best linear unbiasedprediction.Proceedings of the National Academy of Sciences of the UnitedStates of America,111(34):12456-12461.(2014))在GBLUP模型中将所有分子标记作为随机效应放入模型，并假设它们服从同一正态分布，以随机效应整体估算值作为个体育种值，用于预测水稻杂交组合在产量性状上的表现，忽视了主效基因对育种值的贡献。
发明内容
本申请针对现有技术中存在的上述不足，提供了一种用于复杂性状遗传改良的基因组选择方法。
一种用于复杂性状遗传改良的基因组选择方法，包括以下步骤：
(1)统计遗传模型的建立：
将基因分为主效基因和微效基因，建立全基因组主微基因全模型用于表型预测，所述全基因组主微基因全模型为简单加性遗传模型，假设群体大小为n，分子标记数为m，主效基因位点数为c，则所述全基因组主微基因全模型的矩阵形式可以表示为：
Figure BDA0003064545300000021
其中，y为n×1的性状表型值向量，μ为群体均值；ai为第i个主效基因的加性效应，作固定效应；uk为第k个微效基因的加性效应，作随机效应且满足均值为零、方差为
Figure BDA0003064545300000022
的正态分布
Figure BDA0003064545300000023
假设；xi和zk分别为ai和uk的系数向量；b和u分别为固定效应和随机效应向量，X和Z分别为对应的系数矩阵；ε为残差效应向量，作随机效应且服从均值为零、方差为
Figure BDA0003064545300000024
的多元正态分布
Figure BDA0003064545300000025
其中In表示大小为n×n的单位矩阵，ε与u相互独立；
(2)主效基因定位：
首先，对全基因组分子标记进行关联分析，对于第j个分子标记有：
y＝μ+xjaj+ε，
其中，y为n×1的性状表型值向量，μ为群体均值向量；aj为该分子标记的效应值，作固定效应；xj为该分子标记对应的系数向量；ε为残差效应向量，作随机效应且服从多元正态分布
Figure BDA0003064545300000031
其次，利用Henderson III方法在原假设H0：aj＝0下依次对全基因组分子标记位点显著性进行检验，当某一分子标记位点的F统计量值超过阈值时，则认为该标记效应显著，并将该位点作为主效基因选入所述全基因组主微基因全模型中；
(3)遗传参数估计：
在确定了主效基因位点后，对所述全基因组主微基因全模型进行参数估计，并估算个体全基因组育种值；
(4)基于个体全基因组育种值的选择：
在得到所述全基因组主微基因全模型的参数估计值之后，依据个体全基因组育种值大小排序进行基因组选择。
优选的，步骤(1)中，所述全基因组主微基因全模型为简单加性遗传模型时，y的期望值和方差-协方差矩阵分别为：
E(y)＝Xb，
Figure BDA0003064545300000032
其中，亲缘关系矩阵Ka＝ZZT；如果用
Figure BDA0003064545300000033
表示微效基因总效应，此时所述全基因组主微基因全模型可以改写为：
y＝Xb+ξ+ε。
更优选的，将个体全基因组育种值的估算分为两个部分，一个是作为固定效应的主效基因效应，用
Figure BDA0003064545300000034
估算；另一个是作为随机效应的微效基因效应，直接从整体角度得到其总体估计值，即微效基因总效应估计值
Figure BDA0003064545300000035
计算公式为：
Figure BDA0003064545300000036
其中，
Figure BDA0003064545300000037
此时，只需要估算b、
Figure BDA0003064545300000038
和
Figure BDA0003064545300000039
并根据已知系数矩阵Z和X，即可求得
Figure BDA00030645453000000310
并可由此计算个体全基因组育种值的估计值
Figure BDA00030645453000000311
进一步优选的，步骤(3)遗传参数估计时，利用最小范数二阶无偏估计法对随机效应的方差进行估算，得到方差
Figure BDA00030645453000000312
和
Figure BDA00030645453000000313
的估计值
Figure BDA00030645453000000314
和
Figure BDA00030645453000000315
再利用广义最小二乘法对固定效应的效应值进行估算，得到主效基因效应b的估计值
Figure BDA0003064545300000041
优选的，步骤(1)中，所述全基因组主微基因全模型为在简单加性遗传模型的基础上拓展到显性和上位性模型，拓展后为：
Figure BDA0003064545300000042
其中，假设共定位到c个加性或显性主效基因位点以及e对上位性主效基因位点，ai和di分别表示主效基因位点i的加性和显性效应值，(aa)jj′、(dd)jj′、(ad)jj′和(da)jj′分别表示主效基因位点j和主效基因位点j′之间的加加上位性、显显上位性、加显上位性和显加上位性效应值，作固定效应；用ξa、ξd、ξaa、ξdd、ξad和ξda分别表示六种遗传效应的微效基因总效应向量，作随机效应；用符号⊙表示向量与向量之间的点乘；xa和xd分别表示由加性和显性两种不同的基因型编码方式得到的系数向量，对于具有两个等位基因A和a的位点j，第i个个体在该位点的加性编码系数
Figure BDA0003064545300000043
和显性编码系数
Figure BDA0003064545300000044
分别为，
Figure BDA0003064545300000045
用ba、bd、baa、bdd、bad、bda分别表示主效基因的六种遗传效应向量，Xa、Xd、Xaa、Xdd、Xad、Xda为对应系数矩阵。
优选的，步骤(2)中，阈值大小通过置换检验来确定。
更优选的，在一型错误率为0.05时，至少需要进行1000次置换检验，即随机打乱1000次个体基因组信息与表型值的对应关系，获得1000个置换样本，在每一个置换样本中，对所有分子标记依次做F统计检验并取其中最大的F值，根据1000个F最大值构建概率密度分布图，取原假设下达到95％置信度的F值作为显著性阈值。
优选的，使用拟合优度和预测力评价建立的统计遗传模型，
拟合优度衡量统计遗传模型拟合的能力，用
Figure BDA0003064545300000046
表示，公式如下：
Figure BDA0003064545300000047
预测力则衡量统计遗传模型预测的能力，用
Figure BDA0003064545300000051
表示，公式如下：
Figure BDA0003064545300000052
其中，ERESS和PRESS分别表示估计残差平方和以及预测残差平方和，SS为总平方和；
Figure BDA0003064545300000053
和
Figure BDA0003064545300000054
表示第i个个体的表型估计值和预测值，
Figure BDA0003064545300000055
和
Figure BDA0003064545300000056
分别表示第i个个体的残差估计值和残差预测值；
或者，在计算PRESS时，通过HAT矩阵中的中心杠杆值对ERESS进行修正后间接计算PRESS，从而减少重复验证的运算量；HAT矩阵，又称投射矩阵，是表型真实值和表型估计值之间的一个转化矩阵，在固定模型中，
Figure BDA0003064545300000057
H＝X(XTX)-1XT就被称为该固定模型的HAT矩阵，而HAT矩阵对角线上第i个元素hii被称为中心化杠杆值，HAT法中PRESS的计算公式如下：
Figure BDA0003064545300000059
在混合线性模型假设下，HM＝HF+HR(1-HF)，其中HR＝K(K+λI)-1，HF＝X(XTV-1X)-1XTV-1，由此可以快速评价基因组选择精度。
优选的，步骤(4)中个体全基因组育种值的估计从训练群体拓展到训练群体与测试群体，拓展后为：
假设训练群体个体数为n1，测试群体个体数为n2，分子标记数为m，主效位点数为c，且假设这两个群体具有相似的遗传背景，则原模型y＝Xb+Zu+ε可改写为：
Figure BDA00030645453000000510
其中y1(n1×1)和y2(n2×1)分别为训练群体和测试群体的表型值向量，b和u分别为固定效应向量和随机效应向量，随机效应向量u服从均值为零、方差为
Figure BDA00030645453000000511
的多元正态分布
Figure BDA00030645453000000512
Im表示大小为m×m的单位矩阵；X1和X2分别为训练群体和测试群体固定效应的系数矩阵，Z1和Z2分别为训练群体和测试群体随机效应的系数矩阵；ε1和ε2分别为训练群体和测试群体的残差效应向量，作随机效应，且分别服从均值为零、方差为
Figure BDA00030645453000000513
和
Figure BDA00030645453000000514
的多元正态分布
Figure BDA00030645453000000515
Figure BDA00030645453000000516
Figure BDA00030645453000000520
和
Figure BDA00030645453000000517
分别表示大小为n1×n1和n2×n2的单位矩阵，ε1和ε2与u均相互独立；对应的期望值和方差-协方差矩阵分别为：
Figure BDA00030645453000000518
Figure BDA00030645453000000519
依据训练群体和测试群体大小可以将加性亲缘关系矩阵Ka划分为，
Figure BDA0003064545300000061
在已知y1、X1、X2、Z1和Z2的前提下，并通过参数估计获得b、
Figure BDA0003064545300000062
和
Figure BDA0003064545300000063
的估计值，即可求得训练群体个体全基因组育种值的预测值
Figure BDA0003064545300000064
Figure BDA0003064545300000065
对预测所得测试群体个体全基因组育种值的排序可以作为测试群体育种选择的依据。
与现有技术相比，本发明的有益效果为：
a)与GBLUP的模型假设相比，本申请提出的模型假设更符合生物学遗传规律，具有更高的预测力，即具有更高的基因组选择精度；
b)本发明基于混合线性模型，具有很大的灵活性，便于模型扩展。
附图说明
图1为比较MMIBLUP与GBLUP在分析三种不同遗传结构假设下模拟数据的表现结果图。
具体实施方式
在混合线性模型框架下，我们认为复杂性状中的主效基因往往是涉及性状发生最核心环节的基因，基因效应较大，更符合固定效应的假设；相应地，微效基因往往涉及性状发生的修饰环节，基因效应较小，更符合随机效应的假设。但无论是主效基因还是微效基因，都对表型值产生了或多或少的影响，因此有必要共同纳入全模型育种值的估算，以提高基因组选择的精度。
由此，我们提出了一种用于解析复杂性状遗传结构的基因组选择新策略，即针对复杂性状遗传特点，在基因组选择全模型中既包含GWAS鉴别的主效基因作固定效应，又包含GBLUP整体估算的微效基因作随机效应，主微效应既有区分又有整合，共同纳入个体全基因组育种值的估算。该方法包括：
(1)统计遗传模型的建立
根据我们对复杂性状遗传结构的假设，建立全基因组主微基因全模型(major-minor Integrated Best Linear Unbiased Prediction，简称MMIBLUP)用于表型预测。以简单加性遗传模型为例，假设群体大小为n，分子标记数为m，主效基因位点数为c，则所述全基因组主微基因全模型的矩阵形式可以表示为：
Figure BDA0003064545300000071
其中，y为n×1的性状表型值向量，μ为群体均值；ai为第i个主效基因的加性效应，作固定效应；uk为第k个微效基因的加性效应，作随机效应且满足均值为零、方差为
Figure BDA0003064545300000072
的正态分布
Figure BDA0003064545300000073
假设；xi和zk分别为ai和uk的系数向量；b和u分别为固定效应和随机效应向量，X和Z分别为对应的系数矩阵；ε为残差效应向量，作随机效应且服从均值为零、方差为
Figure BDA0003064545300000074
的多元正态分布
Figure BDA0003064545300000075
其中In表示大小为n×n的单位矩阵，ε与u相互独立；
在该模型假设下，y的期望值和方差-协方差矩阵分别为：
E(y)＝Xb(2)
Figure BDA0003064545300000076
其中，亲缘关系矩阵Ka＝ZZT；如果用
Figure BDA0003064545300000077
表示微效基因总效应，此时所述全基因组主微基因全模型可以改写为：
y＝Xb+ξ+ε(4)
(2)主效基因定位
主效基因定位是新策略下全模型建立的前提条件。因此首先对全基因组分子标记进行关联分析，对于第j个分子标记有，
y＝μ+xjaj+ε(5)
aj为该分子标记的效应值，作固定效应，xj为该分子标记对应的系数向量，其他模型的定义与公式(1)中相同。利用Henderson III方法在原假设H0：aj＝0下依次对全基因组分子标记位点显著性进行检验，当某一分子标记位点的F统计量值超过阈值，则认为该标记效应显著，并将该位点作为主效基因选入全模型中，该方法已成功应用于复杂性状QTL定位分析(Yang J,Zhu J,and Williams RW.Mapping the genetic architecture of complextraits in experimental populations.Bioinformatics,23(12):1527-1536.(2007))。
阈值大小通过置换检验来确定(Churchill GA,and Doerge RW.EmpiricalThreshold Values for Quantitative Trait Mapping.Genetics,138(3):963-971.(1994))。在一型错误率为0.05时，至少需要进行1000次置换检验，即随机打乱1000次个体基因组信息与表型值的对应关系，获得1000个置换样本。在每一个置换样本中，对所有分子标记依次做F统计检验并取其中最大的F值。根据1000个F最大值构建概率密度分布图，取原假设下达到95％置信度的F值作为显著性阈值。
我们假设根据Henderson III法F统计量值以及置换检验确定的阈值在全基因组范围确定c个主效基因位点。
(3)遗传参数估计
在确定了c个主效基因位点后，对(1)中全模型进行参数估计。利用最小范数二阶无偏估计(MINQUE)法对随机效应的方差进行估算，得到方差估计值
Figure BDA0003064545300000081
和
Figure BDA0003064545300000082
再利用广义最小二乘法(GLS)对固定效应的效应值进行估算，得到主效基因效应估计值
Figure BDA0003064545300000083
在新模型假设下，我们认为个体全基因组育种值(Genomic Estimated BreedingValue,简称GEBV)的估算应当包括两个部分，一个是作为固定效应的主效基因效应，用
Figure BDA0003064545300000084
估算；另一个是作为随机效应的微效基因效应，直接从整体角度得到其总体估计值，即微效基因总效应估计值
Figure BDA0003064545300000085
计算公式为：
Figure BDA0003064545300000086
其中，
Figure BDA0003064545300000087
此时，只需要估算b、
Figure BDA0003064545300000088
和
Figure BDA0003064545300000089
并根据已知系数矩阵Z和X，即可求得
Figure BDA00030645453000000810
并可由此计算个体全基因组育种值的估计值
Figure BDA00030645453000000811
(4)基因组选择精度的评价
拟合优度(Goodness of fit)和预测力(Predictability)是基因组选择精度(Accuracy)的重要评价指标，也是模拟中用于比较不同模型好坏的常用指标。拟合优度衡量模型拟合的能力，用
Figure BDA00030645453000000812
表示，预测力则衡量模型预测的能力，用
Figure BDA00030645453000000813
表示，它们分别可以用公式(7)和(8)求得，
Figure BDA00030645453000000814
Figure BDA00030645453000000815
其中，ERESS和PRESS分别表示估计残差平方和以及预测残差平方和，SS为总平方和。
Figure BDA00030645453000000816
和
Figure BDA00030645453000000817
表示第i个个体的表型估计值和预测值，
Figure BDA00030645453000000818
和
Figure BDA00030645453000000819
分别表示第i个个体的残差估计值和残差预测值。在计算PRESS时，除了可以通过交叉验证法(Cross Validation,简称CV)直接计算外，还可以通过HAT矩阵中的中心杠杆值对ERESS进行修正后间接计算PRESS，从而减少重复验证的运算量。HAT矩阵，又称投射矩阵，是表型真实值和表型估计值之间的一个转化矩阵。在固定模型中，
Figure BDA00030645453000000820
H＝X(XTX)-1XT就被称为该固定模型的HAT矩阵。而HAT矩阵对角线上第i个元素hii被称为中心化杠杆值，HAT法中PRESS的计算公式如公式(9)所示。
Figure BDA0003064545300000091
HAT法拓展到混合线性模型下时，需根据不同育种值预测策略选择HAT矩阵。这里，用HR、HF和HM分别表示用随机效应、用固定效应以及同时用固定效应和随机效应来预测个体表型值时HAT矩阵的形式(表1)。
表1混合线性模型中HAT矩阵在不同预测方法下的形式
Figure BDA0003064545300000092
(5)基于个体全基因组育种值的选择
基因组选择的一个最终目的是能够根据训练群体(Training Population，简称TRN)的基因型信息和表型信息建立模型，再依据模型及测试群体(Testing Population，简称TST)的基因型信息对测试群体的表型信息进行可靠预测，从而为育种家选择杂交组合及选育优势后代提供参考，因此有必要将个体全基因组育种值的公式拓展到测试群体。
假设训练群体个体数为n1，测试群体个体数为n2，分子标记数为m，主效位点数为c，且假设这两个群体具有相似的遗传背景，则原模型y＝Xb+Zu+ε可改写为，
Figure BDA0003064545300000093
其中y1(n1×1)和y2(n2×1)分别为训练群体和测试群体的表型值向量，b和u分别为固定效应向量和随机效应向量，随机效应向量u服从均值为零、方差为
Figure BDA0003064545300000094
的多元正态分布
Figure BDA0003064545300000095
Im表示大小为m×m的单位矩阵；X1和X2分别为训练群体和测试群体固定效应的系数矩阵，Z1和Z2分别为训练群体和测试群体随机效应的系数矩阵；ε1和ε2分别为训练群体和测试群体的残差效应向量，作随机效应，且分别服从均值为零、方差为
Figure BDA0003064545300000096
和
Figure BDA0003064545300000097
的多元正态分布
Figure BDA0003064545300000098
Figure BDA0003064545300000099
Figure BDA00030645453000000912
和
Figure BDA00030645453000000910
分别表示大小为n1×n1和n2×n2的单位矩阵，ε1和ε2与u均相互独立；对应的期望值和方差-协方差矩阵分别为：
Figure BDA00030645453000000911
Figure BDA0003064545300000101
依据训练群体和测试群体大小可以将加性亲缘关系矩阵Ka划分为，
Figure BDA0003064545300000102
在已知y1、X1、X2、Z1和Z2的前提下，并通过参数估计获得b、
Figure BDA0003064545300000103
和
Figure BDA0003064545300000104
的估计值，即可求得训练群体个体全基因组育种值的预测值
Figure BDA0003064545300000105
Figure BDA0003064545300000106
由此可在只知道测试群体基因型信息的前提下，依赖训练群体和测试群体之间的亲缘关系矩阵对测试群体个体全基因组育种值进行预测。将群体内个体按照全基因组育种值大小从高到低排列，为育种选择提供了有价值的参考意见。
(6)显性和上位性遗传模型的拓展
对于复杂性状遗传结构的解析仅仅依靠加性遗传模型是不够的。因此，在加性遗传模型的基础上，可以将全基因组主微基因全模型拓展到显性和上位性模型，拓展后为：
Figure BDA0003064545300000107
其中，假设共定位到c个加性或显性主效基因位点以及e对上位性主效基因位点，ai和di分别表示主效基因位点i的加性和显性效应值，(aa)jj′、(dd)jj′、(ad)jj′和(da)jj′分别表示主效基因位点j和主效基因位点j′之间的加加上位性、显显上位性、加显上位性和显加上位性效应值，作固定效应；用ξa、ξd、ξaa、ξdd、ξad和ξda分别表示六种遗传效应的微效基因总效应向量，作随机效应；用符号⊙表示向量与向量之间的点乘；xa和xd分别表示由加性和显性两种不同的基因型编码方式得到的系数向量，对于具有两个等位基因A和a的位点j，第i个个体在该位点的加性编码系数
Figure BDA0003064545300000111
和显性编码系数
Figure BDA0003064545300000112
分别为，
Figure BDA0003064545300000113
用ba、bd、baa、bdd、bad、bda分别表示主效基因的六种遗传效应向量，Xa、Xd、Xaa、Xdd、Xad、Xda为对应系数矩阵。
(7)其他
与现有技术相比，本发明的有益效果为：
a)与GBLUP的模型假设相比，我们提出的模型假设更符合生物学遗传规律，具有更高的预测力，即具有更高的基因组选择精度；
b)本发明基于混合线性模型，具有很大的灵活性，便于模型扩展。
下面结合具体实例对本发明作进一步阐释。
(1)材料
我们通过模拟来测试基因组选择新策略MMIBLUP的表现，并与已有基因组选择方法GBLUP进行比较。假设模拟群体大小为n，分子标记数为m，标记间相互独立，MAF的范围在0.05-0.5，得到了n个个体m个位点的基因型。为了更充分比较两种方法在不同遗传结构下的表现，我们共模拟了三种不同遗传假设类型下的复杂性状，包括主效基因遗传(ScenarioI)、主效基因与微效基因共同遗传(Scenario II和Scenario III)和微效基因遗传(Scenario IV)。
在第一种主效基因遗传(Scenario I)假设下，认为全基因组中只有小部分位点对表型有影响。因此在m个位点中随机选取ml个位点作为主效基因位点，并假设其服从同一正态分布
Figure BDA0003064545300000114
遗传率h2分别取0.1、0.2、0.5和0.8来模拟性状的不同遗传水平。根据遗传率公式
Figure BDA0003064545300000115
求得相应的残差方差
Figure BDA0003064545300000116
并根据正态分布
Figure BDA0003064545300000117
模拟得到个体残差大小。假设主效基因位点方差
Figure BDA0003064545300000118
ml占m的1％。
在第二种主效基因与微效基因共同遗传(Scenario II和ScenarioIII)假设下，认为全基因组中除了有小部分主效基因参与性状遗传之外，还有大部分的微效基因共同参与性状遗传。因此，除了假设随机选取的ml个主效基因位点服从正态分布
Figure BDA0003064545300000119
之外，还假设剩余ms个微效基因位点服从另一个正态分布
Figure BDA00030645453000001110
根据主效基因变异大小占全基因组变异(包括主效基因和微效基因变异)大小的比例，可以分为占比中等的Scenario II和占比较小的ScenarioIII两种情况。在Scenario II下，假设ml占m的1％，且
Figure BDA0003064545300000121
ms＝m-ml，且
Figure BDA0003064545300000122
在Scenario III下，同样假设ml占m的1％，且
Figure BDA0003064545300000123
ms＝m-ml，但
Figure BDA0003064545300000124
此时，两种情况主效基因变异大小占全基因组变异大小的比例
Figure BDA0003064545300000125
分别约为0.50和0.25。此时的遗传率计算公式为
Figure BDA0003064545300000126
分别求出h2分别取0.1、0.2、0.5和0.8时对应的残差方差大小
Figure BDA0003064545300000127
在第三种微效基因遗传(Scenario IV)假设下，认为全基因组所有位点均为微效基因位点，假设所有位点的加性效应大小均服从同一正态分布
Figure BDA0003064545300000128
这里我们假设加性效应方差
Figure BDA0003064545300000129
再根据设定的遗传率大小及遗传率公式
Figure BDA00030645453000001210
求得残差效应的方差大小，遗传率h2同样分别取0.1、0.2、0.5和0.8。
在四种情况下，模拟的主效基因变异占全基因组变异(包括主效基因和微效基因)的比例依次递减。每组模拟数据重复20次，群体大小n＝2000，分子标记数m＝3000，并对表型值进行标准化，对基因型系数矩阵按列进行标准化。
(2)软件
在已有软件QTXNetwork(http://ibi.zju.edu.cn/software/)的基础上，利用C++语言编写QTXNetwork-GS模块。在QTXNetwork-GS模块中，可以实现对二倍体群体三种不同遗传假设下复杂性状主效基因和微效基因的模拟，主效基因定位及其效应估计，基因组选择全模型的拟合和预测，以及拟合和预测能力的统计学评价。
(3)结果
我们比较了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。
(4)其他
最后，还需要特别注意的是，以上所举例子仅是本发明的具体实施例子。显然，本发明不仅仅限于以上实施例子，还可以有许多变通的情况。本领域的技术人员从本发明公开的内容直接推导出或联想到的所有变通情况，均认为是本发明的保护范围。

Claims (9)

1.一种用于复杂性状遗传改良的基因组选择方法，其特征在于，包括以下步骤：
(1)统计遗传模型的建立：
将基因分为主效基因和微效基因，建立全基因组主微基因全模型用于表型预测，所述全基因组主微基因全模型为简单加性遗传模型，假设群体大小为n，分子标记数为m，主效基因位点数为c，则所述全基因组主微基因全模型的矩阵形式可以表示为：
Figure FDA0003064545290000011
其中，y为n×1的性状表型值向量，μ为群体均值向量；ai为第i个主效基因的加性效应，作固定效应；uk为第k个微效基因的加性效应，作随机效应且满足均值为零、方差为
Figure FDA0003064545290000012
的正态分布
Figure FDA0003064545290000013
xi和zk分别为ai和uk的系数向量；b和u分别为固定效应和随机效应向量，X和Z分别为对应的系数矩阵；ε为残差效应向量，作随机效应且服从均值为零、方差为
Figure FDA0003064545290000014
的多元正态分布
Figure FDA0003064545290000015
其中In表示大小为n×n的单位矩阵，ε与u相互独立；
(2)主效基因定位：
首先，对全基因组分子标记进行关联分析，对于第j个分子标记有：
y＝μ+xjaj+ε，
其中，y为n×1的性状表型值向量，μ为群体均值向量；aj为该分子标记的效应值，作固定效应；xj为该分子标记对应的系数向量；ε为残差效应向量，作随机效应且服从多元正态分布
Figure FDA0003064545290000016
其次，利用Henderson III方法在原假设H0：aj＝0下依次对全基因组分子标记位点显著性进行检验，当某一分子标记位点的F统计量值超过阈值时，则认为该标记效应显著，并将该位点作为主效基因选入所述全基因组主微基因全模型中；
(3)遗传参数估计：
在确定了主效基因位点后，对所述全基因组主微基因全模型进行参数估计，并估算个体全基因组育种值；
(4)基于个体全基因组育种值的选择：
在得到所述全基因组主微基因全模型的参数估计值之后，依据个体全基因组育种值大小排序进行基因组选择。
2.如权利要求1所述的基因组选择方法，其特征在于，步骤(1)中，所述全基因组主微基因全模型为简单加性遗传模型时，y的期望值和方差-协方差矩阵分别为：
E(y)＝Xb，
Figure FDA0003064545290000021
其中，亲缘关系矩阵Ka＝ZZT；如果用
Figure FDA0003064545290000022
表示微效基因总效应，此时所述全基因组主微基因全模型可以改写为：
y＝Xb+ξ+ε。
3.如权利要求2所述的基因组选择方法，其特征在于，将个体全基因组育种值的估算分为两个部分，一个是作为固定效应的主效基因效应，用
Figure FDA0003064545290000023
估算；另一个是作为随机效应的微效基因效应，直接从整体角度得到其总体估计值，即微效基因总效应估计值
Figure FDA0003064545290000024
计算公式为：
Figure FDA0003064545290000025
其中，
Figure FDA0003064545290000026
此时，只需要估算b、
Figure FDA0003064545290000027
和
Figure FDA0003064545290000028
并根据已知系数矩阵Z和X，即可求得
Figure FDA0003064545290000029
并可由此计算个体全基因组育种值的估计值
Figure FDA00030645452900000210
4.如权利要求3所述的基因组选择方法，其特征在于，步骤(3)遗传参数估计时，利用最小范数二阶无偏估计法对随机效应的方差进行估算，得到方差
Figure FDA00030645452900000211
和
Figure FDA00030645452900000212
的估计值
Figure FDA00030645452900000213
和
Figure FDA00030645452900000214
再利用广义最小二乘法对固定效应的效应值进行估算，得到主效基因效应向量b的估计值
Figure FDA00030645452900000215
5.如权利要求2所述的基因组选择方法，其特征在于，步骤(1)中，所述全基因组主微基因全模型在简单加性遗传模型的基础上拓展到显性和上位性模型，拓展后为：
Figure FDA00030645452900000216
其中，假设共定位到c个加性或显性主效基因位点以及e对上位性主效基因位点，ai和di分别表示主效基因位点i的加性和显性效应值，(aa)jj′、(dd)jj′、(ad)jj′和(da)jj′分别表示主效基因位点j和主效基因位点j′之间的加加上位性、显显上位性、加显上位性和显加上位性效应值，作固定效应；用ξa、ξd、ξaa、ξdd、ξad和ξda分别表示六种遗传效应的微效基因总效应向量，作随机效应；用符号⊙表示向量与向量之间的点乘；xa和xd分别表示由加性和显性两种不同的基因型编码方式得到的系数向量，对于具有两个等位基因A和a的位点j，第i个个体在该位点的加性编码系数
Figure FDA0003064545290000031
和显性编码系数
Figure FDA0003064545290000032
分别为，
Figure FDA0003064545290000033
用ba、bd、baa、bdd、bad、bda分别表示主效基因的六种遗传效应向量，Xa、Xd、Xaa、Xdd、Xad、Xda为对应系数矩阵。
6.如权利要求1所述的基因组选择方法，其特征在于，步骤(2)中，阈值大小通过置换检验来确定。
7.如权利要求6所述的基因组选择方法，其特征在于，在一型错误率为0.05时，至少需要进行1000次置换检验，即随机打乱1000次个体基因组信息与表型值的对应关系，获得1000个置换样本，在每一个置换样本中，对所有分子标记依次做F统计检验并取其中最大的F值，根据1000个F最大值构建概率密度分布图，取原假设下达到95％置信度的F值作为显著性阈值。
8.如权利要求1所述的基因组选择方法，其特征在于，使用拟合优度和预测力评价建立的统计遗传模型，
拟合优度衡量统计遗传模型拟合的能力，用
Figure FDA0003064545290000034
表示，公式如下：
Figure FDA0003064545290000035
预测力则衡量统计遗传模型预测的能力，用
Figure FDA0003064545290000036
表示，公式如下：
Figure FDA0003064545290000037
其中，ERESS和PRESS分别表示估计的、预测的残差平方，SS为总平方和；
Figure FDA0003064545290000038
和
Figure FDA0003064545290000039
表示第i个个体的表型估计值和预测值，
Figure FDA00030645452900000310
和
Figure FDA00030645452900000311
分别表示第i个个体的残差估计值和残差预测值；
或者，在计算PRESS时，通过HAT矩阵中的中心杠杆值对ERESS进行修正后间接计算PRESS，从而减少重复验证的运算量；HAT矩阵，又称投射矩阵，是表型真实值和表型估计值之间的一个转化矩阵，在固定模型中，
Figure FDA00030645452900000312
H＝X(XTX)-1XT就被称为该固定模型的HAT矩阵，而HAT矩阵对角线上第i个元素hii被称为中心化杠杆值，HAT法中PRESS的计算公式如下：
Figure FDA0003064545290000041
在混合线性模型假设下，HM＝HF+HR(1-HF)，其中HR＝K(K+λI)-1，HF＝X(XTV-1X)-1XTV-1，由此可以快速评价基因组选择精度。
9.如权利要求1所述的基因组选择方法，其特征在于，步骤(4)中个体全基因组育种值的估计从训练群体拓展到训练群体与测试群体，拓展后为：
假设训练群体个体数为n1，测试群体个体数为n2，分子标记数为m，主效位点数为c，且假设这两个群体具有相似的遗传背景，则原模型y＝Xb+Zu+ε可改写为：
Figure FDA0003064545290000042
其中y1(n1×1)和y2(n2×1)分别为训练群体和测试群体的表型值向量，b和u分别为固定效应向量和随机效应向量，随机效应向量u服从均值为零、方差为
Figure FDA0003064545290000043
的多元正态分布
Figure FDA0003064545290000044
Im表示大小为m×m的单位矩阵；X1和X2分别为训练群体和测试群体固定效应的系数矩阵，Z1和Z2分别为训练群体和测试群体随机效应的系数矩阵；ε1和ε2分别为训练群体和测试群体的残差效应向量，作随机效应，且分别服从均值为零、方差为
Figure FDA0003064545290000045
和
Figure FDA0003064545290000046
的多元正态分布
Figure FDA0003064545290000047
Figure FDA0003064545290000048
Figure FDA0003064545290000049
和
Figure FDA00030645452900000410
分别表示大小为n1×n1和n2×n2的单位矩阵，ε1和ε2与u均相互独立；对应的期望值和方差-协方差矩阵分别为：
Figure FDA00030645452900000411
Figure FDA00030645452900000412
依据训练群体和测试群体大小可以将加性亲缘关系矩阵Ka划分为，
Figure FDA00030645452900000413
在已知y1、X1、X2、Z1和Z2的前提下，并通过参数估计获得b、
Figure FDA00030645452900000414
和
Figure FDA00030645452900000415
的估计值，即可求得训练群体个体全基因组育种值的预测值
Figure FDA00030645452900000416
Figure FDA00030645452900000417
对预测所得测试群体个体全基因组育种值的排序作为测试群体育种选择的依据。
CN202110522399.XA 2021-05-13 2021-05-13 一种用于复杂性状遗传改良的基因组选择方法 Active CN113223606B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110522399.XA CN113223606B (zh) 2021-05-13 2021-05-13 一种用于复杂性状遗传改良的基因组选择方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110522399.XA CN113223606B (zh) 2021-05-13 2021-05-13 一种用于复杂性状遗传改良的基因组选择方法

Publications (2)

Publication Number Publication Date
CN113223606A true CN113223606A (zh) 2021-08-06
CN113223606B CN113223606B (zh) 2022-05-27

Family

ID=77095405

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110522399.XA Active CN113223606B (zh) 2021-05-13 2021-05-13 一种用于复杂性状遗传改良的基因组选择方法

Country Status (1)

Country Link
CN (1) CN113223606B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114898809A (zh) * 2022-04-11 2022-08-12 中国科学院数学与系统科学研究院 适用复杂性状的基因-环境交互的分析方法及存储介质
CN116863998A (zh) * 2023-06-21 2023-10-10 扬州大学 一种基于遗传算法的全基因组预测方法及其应用
CN117238363A (zh) * 2023-10-25 2023-12-15 青岛极智医学检验实验室有限公司 一种表型预测方法、预测系统、设备及介质

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101338341A (zh) * 2008-06-16 2009-01-07 南京农业大学 不结球白菜高Vc含量基因的分子标记方法
EP2449488A1 (en) * 2009-07-02 2012-05-09 BIOCRATES Life Sciences AG Method for normalization in metabolomics analysis methods with endogenous reference metabolites
CN103229693A (zh) * 2013-02-01 2013-08-07 南京林业大学 一种杉木优良个体的选择方法
CN103451183A (zh) * 2013-09-18 2013-12-18 中国科学院遗传与发育生物学研究所 与小麦小穗数主效基因位点紧密连锁的分子标记及其获取方法和应用
WO2016089815A1 (en) * 2014-12-01 2016-06-09 Danisco Us Inc Fungal host strains, dna constructs, and methods of use
CN109101786A (zh) * 2018-08-29 2018-12-28 广东省农业科学院动物科学研究所 一种整合显性效应的基因组育种值估计方法
CN109524059A (zh) * 2018-12-28 2019-03-26 华中农业大学 一种快速稳定的动物个体基因组育种值评估方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101338341A (zh) * 2008-06-16 2009-01-07 南京农业大学 不结球白菜高Vc含量基因的分子标记方法
EP2449488A1 (en) * 2009-07-02 2012-05-09 BIOCRATES Life Sciences AG Method for normalization in metabolomics analysis methods with endogenous reference metabolites
CN103229693A (zh) * 2013-02-01 2013-08-07 南京林业大学 一种杉木优良个体的选择方法
CN103451183A (zh) * 2013-09-18 2013-12-18 中国科学院遗传与发育生物学研究所 与小麦小穗数主效基因位点紧密连锁的分子标记及其获取方法和应用
WO2016089815A1 (en) * 2014-12-01 2016-06-09 Danisco Us Inc Fungal host strains, dna constructs, and methods of use
CN109101786A (zh) * 2018-08-29 2018-12-28 广东省农业科学院动物科学研究所 一种整合显性效应的基因组育种值估计方法
CN109524059A (zh) * 2018-12-28 2019-03-26 华中农业大学 一种快速稳定的动物个体基因组育种值评估方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
ZHIHONG ZHU 等: "Statistical method for mapping QTLs for complex traits based on two backcross populations", 《CHINESE SCIENCE BULLETIN》 *
ZHIWU ZHANG 等: "Mixed linear model approach adapted for genome-wide association studies", 《NATURE GENETICS》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114898809A (zh) * 2022-04-11 2022-08-12 中国科学院数学与系统科学研究院 适用复杂性状的基因-环境交互的分析方法及存储介质
CN116863998A (zh) * 2023-06-21 2023-10-10 扬州大学 一种基于遗传算法的全基因组预测方法及其应用
CN116863998B (zh) * 2023-06-21 2024-04-05 扬州大学 一种基于遗传算法的全基因组预测方法及其应用
CN117238363A (zh) * 2023-10-25 2023-12-15 青岛极智医学检验实验室有限公司 一种表型预测方法、预测系统、设备及介质
CN117238363B (zh) * 2023-10-25 2024-04-16 青岛极智医学检验实验室有限公司 一种表型预测方法、预测系统、设备及介质

Also Published As

Publication number Publication date
CN113223606B (zh) 2022-05-27

Similar Documents

Publication Publication Date Title
CN113223606B (zh) 一种用于复杂性状遗传改良的基因组选择方法
Cros et al. Genomic selection prediction accuracy in a perennial crop: case study of oil palm (Elaeis guineensis Jacq.)
Ly et al. Whole-genome prediction of reaction norms to environmental stress in bread wheat (Triticum aestivum L.) by genomic random regression
Fernandes Júnior et al. Genomic prediction of breeding values for carcass traits in Nellore cattle
Piepho et al. BLUP for phenotypic selection in plant breeding and variety testing
Ukrainetz et al. Assessing the sensitivities of genomic selection for growth and wood quality traits in lodgepole pine using Bayesian models
EP3326093B1 (en) Improved computer implemented method for predicting true agronomical value of a plant
Rio et al. Genomic selection efficiency and a priori estimation of accuracy in a structured dent maize panel
Misztal et al. Emerging issues in genomic selection
Lara et al. Temporal and genomic analysis of additive genetic variance in breeding programmes
CN111883206A (zh) 一种拟合非加性效应的基因组估计育种值的方法
Garrick et al. Genomic prediction and genome-wide association studies in beef and dairy cattle.
Guillaume et al. Estimation by simulation of the efficiency of the French marker-assisted selection program in dairy cattle (Open Access publication)
Van Eeuwijk et al. Mixed model approaches for the identification of QTLs within a maize hybrid breeding program
Maenhout et al. Prediction of maize single-cross hybrid performance: support vector machine regression versus best linear prediction
Quezada et al. Genomic breeding values’ prediction including populational selfing rate in an open-pollinated Eucalyptus globulus breeding population
You et al. Genomic cross prediction for linseed improvement
Das et al. A Bayesian framework for functional mapping through joint modeling of longitudinal and time‐to‐event data
Amaya Martínez et al. Genetic evaluations in cattle using the single-step genomic best linear unbiased predictor
Salehi et al. Assessment of parametric and non-parametric methods for prediction of quantitative traits with non-additive genetic architecture
CN108470112A (zh) 新配杂交组合表型的预测方法
Atefi et al. Accuracy of genomic prediction under different genetic architectures and estimation methods
Santos et al. Genome-wide prediction of maize single-cross performance, considering non-additive genetic effects
RU2729364C1 (ru) Способ оценки селекционного материала по генотипу для подбора желательных родительских пар, идентификации желательных генотипов из популяций и оценки селекционных линий
CN116863998B (zh) 一种基于遗传算法的全基因组预测方法及其应用

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant