CN106022005B - 一种连续性状和阈性状基因组育种值联合估计的贝叶斯方法 - Google Patents

Abstract

本发明公开了一种连续性状和阈性状基因组育种值联合估计的贝叶斯方法。该方法为基于线－阈模型的新贝叶斯方法，称为LT‑BayesCπ，用于连续性状和阈性状联合分析。使用模拟数据和第十四届QTL‑MAS国际研讨会公共数据验证LT‑BayesCπ，其基因组预测的准确性与基于单性状模型的BayesCπ、BayesTCπ进行比较，同时研究了其性能表现的影响因素。本发明结果表明，在所有情形下，LT‑BayesCπ对阈性状的基因组预测准确性比BayesTCπ显著增加，然而对连续性状的准确性与BayesCπ相当。

Description

一种连续性状和阈性状基因组育种值联合估计的贝叶斯方法

技术领域

本发明涉及一种连续性状和阈性状基因组育种值联合估计的贝叶斯方法。

背景技术

随着单核苷酸多态(SNP)芯片和基因型测序技术的发展，许多全基因组多态被用于动植物育种实践。基因组选择可不需个体本身表型信息，而使用全基因组标记信息，预测种畜遗传价值。在经典的基因组预测过程中，首先使用同时有基因型和表型信息的个体组成参群估计SNP效应，然后使用这些估计的效应构建预测方程，基于候选个体的基因型去计算他们的基因组育种值(GEBVs)。因此，在基因组选择中，一个合适的模型是准确预测基因组育种值的关键。

许多贝叶斯模型被提出去估计基因组育种值。在第一篇基因组选择论文中，两个贝叶斯方法(BayesA和BayesB)被提出去估计基因组育种值，被广泛用于后续的基因组选择研究。然而，BayesA和BayesB有两个缺点。一个缺点是，无论基因型或表型的数量，位点特异方差的完全条件后验分布的自由度相对于其先验分布只增加了1，如Gianola等人所指出，SNP效应的压缩率严重依赖于尺度参数。另一个缺点是，对于BayesB，某个位点SNP效应值为零的概率值π为已知。为了解决这两个缺点，Habier等人提出了两个新方法BayesCπ和BayesDπ。对于BayesCπ，一个共同方差赋于非零效应SNP替代位点特异方差，而对于BayesDπ，位点特异方差的先验逆卡方分布的尺度参数作为未知参数，其先验分布为Gamma(1，1)。BayesCπ和BayesDπ都把π作为一个未知参数，从数据推断获得。

Wang等人拓展常规的三种估计连续性状基因组育种值的贝叶斯方法(BayesA、BayesB和BayesCπ)为阈模型方法，相应地称为BayesTA、BayesTB和BayesTCπ，用于阈性状基因组育种值估计。模拟研究结果表明，对于阈性状，这三种新的BayesT方法表现均优于相应的常规贝叶斯方法。BayesTCπ在三种新方法中表现最好，被推荐为阈性状基因组选择的方法。

目前，单性状模型通常被用于基因组预测实践。理论上，对于遗传相关的多性状，联合基因组评估比单性状基因组预测准确性更高。然而，连续性状和阈性状的联合分析尚缺少。本发明提出一个基于线-阈模型的新贝叶斯方法，称为LT-BayesCπ，用于连续性状和阈性状联合分析。使用模拟数据和第十四届QTL-MAS国际研讨会公共数据验证LT-BayesCπ，其基因组预测的准确性与基于单性状模型的BayesCπ、BayesTCπ进行比较，同时研究了其性能表现的影响因素。

发明内容

本发明要解决的技术问题是提供一种连续性状和阈性状基因组育种值联合估计的贝叶斯方法。

为了解决上述技术问题，本发明采用的技术方案是，连续性状和阈性状基因组育种值联合估计的贝叶斯方法，包括以下步骤：

(1)模型

令y′₁＝{y_1，i}，i＝1，2，…，n为连续性状观察值向量，y′₂＝{y_2，i}，i＝1，2，…，n为阈性状观察值向量，l′＝{l_i}，i＝1，2，…，n为与阈性状相关联的潜在变量向量；线-阈模型为：

其中β₁、β₂为固定效应向量；g₁、g₂为SNP效应向量；e₁、e₂为随机残差向量；x₁、x₂为β₁、β₂的关联矩阵；Z是基因型指示矩阵，其中赋值0、1和2分别对应于基因型的11、12和22；

令v′＝[y′₁，l′]，给定β和g时，v服从如下分布：

其中β′＝[β′₁，β′₂]，g′＝[g′₁，g′₂]，

然后给定β，g和R_e，则抽样模型为：

其中

阈性状的潜在变量是观测不到的，其残差方差不可估；因此，在本研究中阈性状的残差方差固定为

(2)LT-BayesCπ的MCMC执行过程

a)先验分布

在本研究中，下面的先验分布是为了构建一个分级模型；

对于“固定”效应β：

β|β_min，β_max～∪(β_min，β_max)；

对SNP效应g：

每个SNP对于两个性状都为零效应和至少对于某一个性状为非零效应的概率分别为π和1-π；对于后者，SNP效应先验分布为二元正态分布，即：

其中

G₀和R_e的先验分布都是逆威夏特分布：

对于阈值t：

其中T＝{(t₁，...，t_k-1)|t_min≤t₁≤…≤t_k-1≤t_max}；

对于SNP效应为零的概率值π：

p(π)～∪(0，1)；

b)联合后验分布

参数向量增加潜在变量1和指示变量δ，表示为(Ω，l，δ)，其中Ω＝(β，g，G₀，R_e，t，π)；

参数向量的联合后验分布为：

p(Ω，l，δ|y₁，y₂)∝p(y₁，y₂|Ω，l，δ)p(Ω，l，δ)

＝p(y₁，l|Ω)P(y₂|Ω，l)p(Ω，δ)；

c)完全条件后验分布

潜在变量

潜在变量l_i的完全条件后验分布是以t_j-1到t_j为截断点的截断正态分布，其期望和方差如下：

d)“固定”效应和SNP效应

建立混合模型方程组(mixed model equations，MME)其中则θ_i的完全条件后验分布为：

其中

e)方差-协方差矩阵

SNP效应和残差效应的方差-协方差矩阵完全条件后验分布为：

其他参数的完全条件后验分布和吉布斯抽样过程与单性状方法BayesTCπ的相似。

本发明的有益效果是：

本发明提出一个基于线-阈模型的新贝叶斯方法，称为LT-BayesCπ，用于连续性状和阈性状联合分析。结果表明，在所有情形下，LT-BayesCπ对阈性状的基因组预测准确性比BayesTCπ显著增加，然而对连续性状的准确性与BayesCπ相当。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1是本发明实施例的标准情形下，随即选取一个重复的QTL模拟效应及SNP估计效应绝对值。面板Q_traitA和Q_trait B显示QTL模拟效应绝对值。面板Cpi_trait A，LTCpi_trait A，TCpi_trait B，LTCpi_trait B分别显示BayesCπ，LT-BayesCπ，BayesTCπ方法对性状A或性状B估计的SNP效应绝对值。

图2是本发明实施例中当QTL数量20-500时，三世代各方法基因组育种值准确性。

图3是本发明实施例中当性状A遗传力0.3-0.8时，三世代各方法基因组育种值准确性。

图4是本发明实施例中当阈性状发生率0.05-0.5时，三世代各方法基因组育种值准确性。

具体实施方式

本实施例提出一个基于线-阈模型的新贝叶斯方法，称为LT-BayesCπ，用于连续性状和阈性状联合分析。

1.方法

1.1模型

令y′₁＝{y_1，i}，i＝1，2，…，n为连续性状观察值向量，y′₂＝{y_2，i}，i＝1，2，…，n为阈性状观察值向量，l′＝{l_i}，i＝1，2，…，n为与阈性状相关联的潜在变量向量。线-阈模型为：

其中β₁、β₂为固定效应向量；g₁、g₂为SNP效应向量；e₁、e₂为随机残差向量；x₁、x₂为β₁、β₂的关联矩阵；Z是基因型指示矩阵，其中赋值0、1和2分别对应于基因型的11、12和22。令v′＝[y′₁，l′]，给定β和g时，v服从如下分布：

其中β′＝[β′₁，β′₂]，g′＝[g′₁，g′₂]，

然后给定β，g和R_e，则抽样模型为：

其中

阈性状的潜在变量是观测不到的，其残差方差不可估。因此，在本研究中阈性状的残差方差固定为

1.2LT-BayesCπ的MCMC执行过程

a)先验分布

在本研究中，下面的先验分布是为了构建一个分级模型。

对于“固定”效应β：

β|β_min，β_max～∪(β_min，β_max)。

对SNP效应g：

每个SNP对于两个性状都为零效应和至少对于某一个性状为非零效应的概率分别为π和1-π。对于后者，SNP效应先验分布为二元正态分布，即：

其中

G₀和R_e的先验分布都是逆威夏特分布：

对于阈值t：

其中T＝{(t₁，...，t_k-1)|t_min≤t₁≤…≤t_k-1≤t_max}。

对于SNP效应为零的概率值π：

p(π)～∪(0，1)。

b)联合后验分布

参数向量增加潜在变量1和指示变量δ，表示为(Ω，l，δ)，其中Ω＝(β，g，G₀，R_e，t，π)。

参数向量的联合后验分布为：

p(Ω，l，δ|y₁，y₂)∝p(y₁，y₂|Ω，l，δ)p(Ω，l，δ)＝p(y₁，l|Ω)p(y₂|Ω，l)p(Q，δ)。

c)完全条件后验分布

潜在变量

潜在变量l_i的完全条件后验分布是以t_j-1到t_j为截断点的截断正态分布，其期望和方差如下：

d)“固定”效应和SNP效应

建立混合模型方程组(mixed model equations，MME)其中则θ_i的完全条件后验分布为：

其中

e)方差-协方差矩阵

SNP效应和残差效应的方差-协方差矩阵完全条件后验分布为：

其他参数的完全条件后验分布和吉布斯抽样过程与单性状方法BayesTCπ的相似。

1.3模拟试验

模拟数据

为评估新提出的方法LT-BayesCπ，我们使用多性状基因组模拟软件GPOPSIM产生的模拟数据开展一系列模拟试验。

数据模拟开始于100个个体构成的基础群，接着是保持群体大小不变，经历没有世代重叠的1000个历史世代，标记为-999至0世代。在每个历史世代，50头公畜与50头母畜随即交配，每个交配产生2个后代(1公、1母)。所有位点标记在基础群中起始都是单一同态的，多态性位点标记在后续的世代中以1.25×10^-3的突变率增生，但会因遗传漂变降低多态，经过一些历史世代达到平衡状态。1000个历史世代后，另外再模拟六个世代，编号为1-6代。在1世代，群体大小由100扩大到1000，由0世代的50头公畜与50头母畜随机交配产生。1-5世代，每世代从500头公畜中随机选出50头公畜作为下一世代的父亲，500头母畜全作为下一世代的母亲。2-6世代，为了保持每一代群体大小数为1000，50头公畜随机与500头母畜交配，每头母畜产生2个后代(1公、1母)。1-2世代为本次试验的参考群，3-6世代为验证群。

模拟的基因组由5条染色体组成，总长5M(1M/染色体)。每条染色体上，均匀分布2000个标记，每两个标记之间包含一个潜在的QTL。最终真实的QTL从潜在的QTL中随机抽取。Haldane图距函数被用于计算相邻位点的重组率。

模拟两个遗传相关的性状(表示为性状A和性状B)，性状A是一个服从正态分布的连续性状，性状B是一个二分类性状，其潜在变量服从正态分布。所有抽样的QTL被分为三组，第一组、第二组和第三组。第一组QTL为多效性同时影响两个性状，而第二组和第三组QTL仅仅分别影响性状A、性状B。第一组中QTL的等位基因替代效应抽自一个根据两性状间遗传相关(r_AB)而变化的二元正态分布，第二组和第三组中QTL的等位基因替代效应抽自一元正态分布。除了r_AB＝0时，第一组、第二组和第三组的QTL数比值设定为0.0∶0.5∶0.5，两性状间相关(r_AB)为其他时，第一组、第二组和第三组的QTL数比值都设定为0.8∶0.1∶0.1。等位基因的替代效应重新尺度化使性状A和性状B的加性遗传方差分别等于2.0和1.0。环境相关设定为0.0，两个性状的环境效应独立地从一元正态分布抽取。

真实育种值(TBV)由所有QTL效应的和构成，加上环境效应，产生性状A的表型值和性状B的潜在变量值。1-6世代所有个体都模拟了基因型和真实育种值，但只有1-2世代的2000个体(参考群)赋予了表型值或潜在变量值。对于性状B，根据假设的发生率设定一个阈值，潜在变量值低于阈值的表型值为0，否则表型值为1。

我们首先模拟了一个标准情形，参数设定为：两性状遗传力和QTL数目60；性状A和B遗传相关r_AB＝0.5；性状B发生率0.30(即，30％的个体表型值为1)。为研究新方法基因组预测的影响因素，相对于标准情形改变有关参数，设置可选情形：遗传相关r_AB(0.0，0.2和0.80)；QTL数目(20，200和500)；连续性状A遗传力：(和0.8)；二分类性状B遗传力：(和0.5)；性状B的发生率(0.05，0.1和0.5)。对于每种情形，均模拟20个重复。

第十四届QTL-MAS国际研讨会公共数据

第十四届QTL-MAS国际研讨会公共数据(Szydlowski M，P：QTLMAS2010：simulated dataset.BMC Proc 2011，5Suppl 3(Suppl 3)：S3.)由连续的五个世代(F0-F4)的3226个个体组成。所有个体都模拟了基因型记录，然而只有F0-F3世代的2326个个体有两个性状的表型记录：一个数量性状Q和一个二分类性状B。有表型记录的个体(F0-F3世代)和没有表型记录的个体(F4世代)分别作为参考群和验证群。基因组由5条染色体构成，每条染色体长100Mbp，共模拟10031个两等位基因的SNP标记，没有任何缺失数据和判型错误。数量性状Q受37个QTL控制(其中，30个加性QTL，4个上位QTL和3个印记QTL)。30个加性QTL中，22个加性QTL影响二分类性状B，且不受其他QTL的控制。侠义遗传力(h²)，对于数量性状Q在公畜和母畜分别为0.52和0.39，而对于二分类性状B遗传力为0.48。两个性状间的遗传相关在公畜和母畜分别0.59和0.68。

SNP效应估计

三种贝叶斯方法都被用于参考群体进行SNP效应估计。新提出的LT-BayesCπ方法用于连续性状和阈性状的联合分析，而BayesCπ和BayesTCπ分别用于连续性状和阈性状的分析。对于每种模型，马尔科夫链都由50000个循环的吉布斯抽样产生，最初的30000个循环作为burn-in期被去掉。SNP效应的估计值由burn-in期之后所有SNP效应抽样值取平均获得。

基因组预测准确性

验证群中个体的基因组育种值是根据它们的标记基因型把所有标记效应求和获得。对于每个性状，真实育种值和基因组育种值的相关系数(r_TBV，GEBV)用于衡量准确性，真实育种值对基因组育种值的回归系数(b_TBV，GEBV)用于评估基因组预测的无偏性。然而，对于二分类性状，因为限定V_e＝1，所以基因组育种值的尺度与真实育种值的不一样。如此，b_TBV，GEBV必须使用校正回原始尺度，其中V_eo是模拟的真实残差方差。此外，实施T检验比较LT-BayesCπ与单性状方法BayesCπ或BayesTCπ之间的准确性差异。

2结果

2.1模拟数据

2.1.1标准情形下SNP效应的估计值

图1展示了在标准情形下随机选择一个重复的QTL模拟效应(图1Q_trait A和图1Q_trait B)，BayesCπ对性状A估计的SNP效应(图1Cpi_trait A)，LT-BayesCπ对性状A估计的SNP效应(图1LTCpi_trait A)，BayesTCπ对性状B估计的SNP效应(图1TCpi_trait B)和LT-BayesCπ对B性状估计的SNP效应(图1LTCpi_trait B)。连续性状A和阈性状B的QTL模拟效应范围分别为0-0.75和0-0.77，BayesCπ和LT-BayesCπ对连续性状A估计的SNP效应范围分别为0-0.55和0-0.70，BayesTCπ和LT-BayesCπ对阈性状B估计的SNP效应范围分别为0-0.14和0-0.21。所有方法都定位到了包含大效应QTL的片段。

2.1.2标准情形下基因组育种值的准确性

表1列出了在标准情形下3-6世代(验证群)基因组育种值的准确性，即GEBVs和模拟真实育种值的相关系数。对于所有的方法，其准确性如预期地随着世代下降，除了LT-BayesCπ对性状B的准确性第5世代略低于第6世代，可能是由于抽样误差。对于阈性状B，在所有世代中LT-BayesCπ均比BayesTCπ准确性高，都约提高了11％(表1)。然而对于连续性状A，LT-BayesCπ和BayesCπ比较，在所有世代中并没有提高。总之，对于性状B，LT-BayesCπ的准确性显著高于BayesTCx(P＜0.001)，对于性状A，其准确性与BayesCπ相当(P＞0.1)。

表1标准情况下3-6世代，三种方法获得的基因组育种值的准确性

注：20个重复的平均数±S.E.

2.1.3遗传相关对基因组预测的影响

表2展示了在四种不同遗传相关(0.00，0.20，0.50和0.80)情形下，三世代基因组预测的准确性。

对于连续性状A，无论两个性状的遗传相关性强弱，LT-BayesCπ和BayesCπ估计基因组育种值的准确性基本相同。表明对于连续性状，LT-BayesCπ和BayesCπ表现相当。同时，也显示随着两性状间遗传相关性的增加，LT-BayesCπ和BayesCπ准确性差异有减小趋势。

对于阈性状B，在两个性状没有遗传相关的情况下，LT-BayesCπ比BayesCπ的基因组预测准确性低5％。然而，随着遗传相关性的增加，LT-BayesCπ基因组预测的准确性急剧增加，这与BayesTCπ的表现不同，也与LT-BayesCπ、BayesCπ对连续性状A的表现不同。随着两性状间遗传相关性的增加，LT-BayesCπ相对于BayesTCπ的准确性改进也逐渐增加，特别当遗传相关系数达0.50和0.80时，前者比后者准确性分别显著提高0.116和0.201(P值＜0.001)。

表2四种不同遗传相关情形下三世代基因组育种值的准确性

注：***P值＜0.001，**P值＜0.01，*P值＜0.05

除估计基因组育种值之外，我们提出的LT-BayesCπ可更好的预测两性状间的遗传相关、残差相关和真实QTL比例(π值)。如表3所示，在所有不同遗传相关情形下，LT-BayesCπ获得π的估计值都是非常接近设定值0.006的。同时，在大部分情况下，遗传相关和残差相关的估计值也近似等于设定值，除了在遗传相关系数设定为0.80时，遗传相关和残差相关的估计值分别比设定值低0.13和高0.02。

表3四种不同遗传相关情形下LT-BayesCπ估计的遗传相关残差相关真实QTL比例

注：设定的残差相关和π分别是0和0.006

模拟的真实育种值对估计的基因组育种值的回归系数用于评价基因组预测的无偏性。如表4所示，对于连续性状A，BayesCπ的回归系数接近于1，LT-BayesCπ基因组预测无偏性稍差点。对于阈性状B，LT-BayesCπ和BayesTCπ产生了很大的偏差，尽管LT-BayesCπ比BayesTCπ表现得好很多，但在所有情况下他们得到的回归系数都比1.0大。

表4四种不同遗传相关情形下三世代模拟的真实育种值对估计的基因组育种值的回归系数

2.1.4QTL数目对基因组预测的影响

如图2所示，BayesCπ、BayesTCπ和LT-BayesCπ都对影响性状的QTL数目敏感，它们的基因组预测准确性随着QTL数目增加而迅速降低。当QTL数目从20增加到500时，BayesCπ、BayesTCπ和LT-BayesCπ(性状A)、LT-BayesCπ(性状B)的GEBVs准确性分别减少0.141，0.131，0.156和0.169。在相同情形下，我们新提出LT-BayesCπ方法，对于连续性状A，与单性状BayesCπ方法表现相当，而对于阈性状B，比单性状BayesTCπ方法表现好。

在QTL数目为20、60、200和500时，LT-BayesCπ和BayesCπ对连续性状A的准确性差异分别是0.002，0.010，0.012和0.017。然而，不同的是对于阈性状B，在所有的情形下，LT-BayesCπ的基因组预测准确性比BayesTCπ，分别高出0.113，0.116，0.095，0.075，表明LT-BayesCπ显著优于BayesTCπ(P＜0.001)

2.1.5连续性状A的遗传力对基因组预测的影响

图3显示了不同的性状A遗传力情况下三种方法在三世代基因组育种值准确性(遗传相关为0.5)。正如预期，随着性状A遗传力从0.3增加到0.8，BayesCπ和LT-BayesCπ(性状A)的基因组育种值准确性也持续增加。在所有情形下，BayesTCπ(性状B)的基因组育种值准确性都没有变化。然而，当性状A遗传力从0.3增加到0.8时，LT-BayesCπ(性状B)的基因组育种值准确性从0.581增加到0.632。这表明连续性状A遗传力的增加有利于提高LT-BayesCπ对阈性状B的基因组预测准确性。

2.1.6阈性状B的发生率对基因组预测的影响

图4显示了阈性状B不同发生率情况下各方法估计基因组育种值的准确性。当发生率从5％上升到50％时，LT-BayesCπ和BayesTCπ方法对阈性状B估计基因组育种值的准确性持续增加。然而，连续性状A的基因组预测准确性不受阈性状B发生率的影响。在阈性状B不同发生率情况下，LT-BayesCπ和BayesCπ对连续性状A估计基因组育种值准确性没有变化。

2.2第十四届QTL-MAS国际研讨会公共数据

为评估我们新提出的LT-BayesCπ方法的效果，我们也分析了第十四届QTL-MAS国际研讨会公共数据，该数据包含相关的连续性状Q和二分类性状B。对于每种贝叶斯方法，使用不同的随机数重复分析了10次。表5展示了验证群中基因组育种值的准确性和无偏性。对于连续性状Q，LT-BayesCπ方法比BayesCπ的基因组预测准确性有微小提高，虽然在我们的模拟试验中LT-BayesCπ和BayesCπ的准确性相当。同样地，和BayesTCπ相比，LT-BayesCπ方法对二分类性状B的估计基因组育种值准确性(P＜0.001)和无偏性都有明显改进，这和我们的模拟试验结果相似。

表5第十四届QTL-MAS国际研讨会公共数据验证群中三种方法的准确性和无偏性

注：不同随机数进行10次计算结果的平均数及标准误

3讨论

在传统的遗传评估中，多性状模型是为了利用遗传相关性状的信息来增加估计育种值准确性，已被广泛用于育种值估计实践。自从基因组选择概念于2001年被提出，许多模型，如GBLUP，BayesA，BayesB，BayesCπ，被开发用于基因组预测，并且大部分研究关注于单个连续性状基因组预测。最近，有若干研究将BayesA、BayesB、BayesCπ拓展用于单个阈性状。最近，一些研究将性状间的相关结构纳入考虑，多个连续性状进行联合基因组预测，结果表明可获得更准确的基因组预测。然而，连续性状与阈性状/二分类性状进行联合基因组预测的研究仍缺少。本实施例开发了新方法LT-BayesCπ去处理遗传相关的连续性状和阈性状的联合基因组预测。模拟数据和第十四届QTL-MAS国际研讨会公共数据的研究结果表明，在所有情形下，当使用LT-BayesCπ对连续性和阈性状进行联合分析，对于阈性状，LT-BayesCπ的准确性和无偏性都显著优于单性状方法BayesTCπ，而对于连续性状，其准确性和无偏性与单性状方法BayesCπ相似。

性状间遗传相关是多性状联合分析获得好效果的基本条件。两性状间遗传相关是由于多效QTL同时影响两个性状和/或影响不同性状QTL的连锁不平衡。在我们的模拟研究中，性状间的遗传相关主要由于共同QTL的多效性。模拟研究结果表明，性状间遗传相关越大，联合分析获益更多(表2)。然而，当遗传相关为零时，对于两个性状都是联合基因组预测的准确性低于分开的单个性状分析。Jia和Jannink也报道了该现象，对于没有遗传相关的两个连续性状，两性状模型比单性模型的表现差。原因可能是，在此情形下，多性状模型会导致相关性的非零估值，从而错误信息会在性状间共享。

在该模拟研究中，绝大多数情形下，连续性状没有从联合分析受益。一个原因是二分类性状的遗传力低(0.1)，不能提供足够的信息帮助连续性状提高准确性。其他的连续性状研究也获得类似结果。Jia和Jannink模拟了两个遗传力分别为0.1和0.5的连续性状，对于遗传力为0.5的性状，其准确性没有提高，而对于遗传力为0.1的性状，其准确性显著提高。本实施例中，对于第十四届QTL-MAS国际研讨会公共数据，两个性状的遗传力分别为和都在0.5附近，使用LT-BayesCπ进行联合分析，其连续性状和二分类性状的准确性都比单性状分析提高了，表明新提出的方法对于具有高遗传力的连续性状和二分类性状都有提高预测准确性的潜力。

控制性状的QTL数目影响贝叶斯方法基因组预测的准确性成为共识，我们的研究结果也再次证实此点。所有方法的准确性都随着QTL数目的上升而下降。当QTL数目上升，平均的每个QTL效应值变小，假如遗传方差固定，那么在该参考群体中估计SNP效应的准确性会降低。

当连续性状的遗传力上升，LT-BayesCπ的准确性不但对于连续性状增加，而且对于阈性状也增加。当二分类性状的遗传力上升，LT-BayesCπ的准确性不但对于二分类性状增加，而且对于连续性状也增加(图3)。研究结果表明，低遗传力性状可从相关的高遗传力性状获得信息，从而得到更高的预测准确性，Jia和Jannink和Guo等也观察到类似情况。该研究结果也与传统的遗传评估研究结论相一致，使用多性状模型的益处对于低遗传力性状更重要。

随着二分类性状的发生率向0.5增加，其基因组育种值的准确性也增加。原因在于小发生率的性状需要更大的参考群体去估计方差组分和获得足够的基因组育种值准确性。另一方面，连续性状基因组预测的准确性不受二分类性状发生率变化的影响(图4)。这可能由于LT-BayesCπ中的方差-协方差矩阵没有改变，从而对连续性状的影响是可忽略的。

4结论

本实施例的研究表明基于线-阈模型的LT-BayesCπ适用于连续性状和阈性状基因组育种值的联合预测。特别地，相对于单性状分析，使用LT-BayesCπ进行联合分析可显著提高二分类性状的准确性。两个性状间的遗传相关越大，获益会越多。连续性状或/和二分类性状的遗传力增大对两个性状准确性的提高都有利，尤其对于二分类性状。二分类性状发生率仅对二分类性状的预测准确性有影响。LT-BayesCπ是连续性状和阈性状联合分析的一个可选方法。

以上所述的本发明实施方式，并不构成对本发明保护范围的限定。任何在本发明的精神和原则之内所作的修改、等同替换和改进等，均应包含在本发明的权利要求保护范围之内。

Claims (1)

1.连续性状和阈性状基因组育种值联合估计的贝叶斯方法，包括以下步骤：

(1)模型

令y′₁＝{y_1，i}，i＝1，2，…，n为连续性状观察值向量，y′₂＝{y_2，i}，i＝1，2，…，n为阈性状观察值向量，l′＝{l_i}，i＝1，2，…，n为与阈性状相关联的潜在变量向量；线-阈模型为：

其中β₁、β₂为固定效应向量；g₁、g₂为SNP效应向量；e₁、e₂为随机残差向量；x₁、x₂为β₁、β₂的关联矩阵；z是基因型指示矩阵，其中赋值0、1和2分别对应于基因型的11、12和22；

令v′＝[y′₁，l′]，给定β和g时，v服从如下分布：

其中β′＝[β′₁，β′₂]，g′＝[g′₁，g′₂]，

然后给定β，g和R_e，则抽样模型为：

其中

阈性状的潜在变量是观测不到的，其残差方差不可估；因此，阈性状的残差方差固定为

(2)LT-BayesCπ的MCMC执行过程

a)先验分布

下面的先验分布是为了构建一个分级模型；

对于“固定”效应β：

β|β_min，β_max～∪(β_min，β_max)；

对SNP效应g：

每个SNP对于两个性状都为零效应和至少对于某一个性状为非零效应的概率分别为π和1-π；对于后者，SNP效应先验分布为二元正态分布，即：

其中

G₀和R_e的先验分布都是逆威夏特分布：

对于阈值t：

其中T＝{(t₁，...，t_k-1)|t_min≤t₁≤…≤t_k-1≤t_max}；

对于SNP效应为零的概率值π：

p(π)～∪(0，1)；

b)联合后验分布

参数向量增加潜在变量1和指示变量δ，表示为(Ω，l，δ)，其中Ω＝(β，g，G₀，R_e，t，π)；

参数向量的联合后验分布为：

p(Ω，l，δ|y₁，y₂)∝p(y₁，y₂|Ω，l，δ)p(Ω，l，δ)

＝p(y₁，l|Ω)p(y₂|Ω，1)p(Ω，δ)；

c)完全条件后验分布

潜在变量

潜在变量l_i的完全条件后验分布是以t_j-1到t_j为截断点的截断正态分布，其期望和方差如下：

d)“固定”效应和SNP效应

建立混合模型方程组其中则θ_i的完全条件后验分布为：

其中

e)方差-协方差矩阵

SNP效应和残差效应的方差-协方差矩阵完全条件后验分布为：

其他参数的完全条件后验分布和吉布斯抽样过程与单性状方法BayesTCπ的相似。