动物育种中育种值方差组分估计方法、装置、设备和介质
本发明涉及动物育种，尤其涉及一种动物育种中育种值方差组分估计方法、装置、设备和介质。
背景技术：
1、在猪、牛、鸡等家畜家禽育种以及作物基因组选择中，blup、gblup、ssgblup 和线性混合模型被广泛用于估计育种值和方差组分。方差组分估计通常采用 reml、em-reml 或ai-reml 方法，估计结果直接影响遗传力、显性方差、残差方差和育种值可靠性。
2、随着基因芯片和测序成本下降，单个育种群体中的基因型个体数量和标记数量持续增长。传统 reml 实现需要处理表型协方差矩阵、投影矩阵、关系矩阵或其他关系核及其逆矩阵。在大规模场景中，显式构建或求逆这些稠密矩阵会造成较高存储压力和求解压力。
3、已有规避亲缘关系矩阵求逆的 blup 模型求解方案能够在单次 blup 求解中避免关系矩阵逆，并得到育种值。然而，reml 方差组分估计不是一次 blup 求解，而是在每轮迭代中反复计算得分向量、迹项和 average information （平均信息）矩阵。即使底层blup 求解已经免逆，ai-reml 中仍存在投影矩阵作用于工作向量的计算瓶颈。
4、此外，em-reml 与 ai-reml 在工程行为上具有互补性但属于两条不同迭代线路：ai-reml 更新通常收敛较快，是本方案中的优先主路径；em-reml 更新由其期望最大化原理决定，通常需要更多收敛轮次，因而不作为常规前置步骤，仅在 ai-reml 难以完成收敛的少数极端情形下作为补充路径。已有技术通常将二者作为不同算法实现，难以在无逆关系核、随机迹估计和多右端批量求解框架下复用同一底层求解资源。因此，需要一种能够在同一 if 求解体系中分别支持 ai 快速主路径和 em 补充稳定路径的 reml 方差组分估计方案。
5、在 ai-reml 中，将方差组分方向导数矩阵作用于 py 后写成工作向量属于本领域常规数学表达。真正需要解决的是如何在不显式形成投影矩阵、表型协方差逆矩阵或关系矩阵逆的情况下，稳定得到投影矩阵作用于工作向量后的结果。本发明使用工作向量伪表型求解和伪残差提取解决该瓶颈。
技术实现思路
1、本发明提供一种动物育种中育种值方差组分估计方法、装置、设备和介质，要解决的技术问题是：在大规模动物育种和基因组遗传评估中，如何在不显式构建投影矩阵、表型协方差逆、关系矩阵逆以及稠密矩阵乘法的前提下，将 ai-reml 快速主路径与 em-reml补充稳定路径作为两条独立线路统一到同一无逆求解体系中，使二者共享无逆求解器、关系核算子、随机迹估计和伪残差提取机制，并使 rmom 初值能够同时服务于两条路径。
2、第一方面，本发明提供一种动物育种中育种值方差组分估计方法，所述方法包括：
3、获取育种群体的观测表型向量、固定效应设计矩阵、随机效应设计矩阵以及基因型标记矩阵；
4、根据所述观测表型向量、固定效应设计矩阵和随机效应设计矩阵，建立包括固定效应、至少一个随机效应和残差的混合模型，并根据所述基因型标记矩阵构建各方差组分方向的关系核矩阵；
5、根据所述混合模型对应的混合模型方程构建无逆if-mme求解器，所述无逆if-mme求解器通过对所述混合模型方程进行等价变换以规避关系核逆矩阵的显式构建，并以任意表型空间向量作为右端项输入，输出对应的固定效应解、随机效应解和残差；
6、初始化方差组分，并以所述观测表型向量作为右端项输入所述无逆if-mme求解器，根据所述无逆if-mme求解器的输出得到真实残差和随机效应预测值；
7、根据所述真实残差和随机效应预测值构造各方差组分方向的工作向量；
8、将所述工作向量作为伪表型输入所述无逆if-mme求解器，根据所述无逆if-mme求解器的输出得到固定效应伪解和随机效应伪解；
9、根据所述伪表型、固定效应伪解和随机效应伪解提取伪残差，所述伪残差用于表示所述reml投影矩阵对所述工作向量的作用结果；
10、根据所述工作向量和所述伪残差装配平均信息矩阵元素，并通过随机迹估计方法形成得分向量；
11、根据所述平均信息矩阵和所述得分向量求解小规模线性系统，得到方差组分更新步长，并根据所述更新步长更新当前方差组分；
12、重复执行以观测表型向量作为右端项输入无逆if-mme求解器至根据平均信息矩阵和得分向量更新当前方差组分的迭代过程直至满足收敛判据，输出估计得到的加性遗传方差组分、残差方差组分以及由所述加性遗传方差组分和残差方差组分计算的遗传力。
13、在一种可能的设计中所述混合模型表示为：
14、 （1）
15、 （2）
16、其中，为观测表型向量，为固定效应设计矩阵，为固定效应向量，为第个随机效应方向的设计矩阵，为第个随机效应方向的随机效应向量，为残差向量；为第个随机效应方向的关系核矩阵，为第个方差组分，为残差协方差矩阵，为正态分布。
17、在一种可能的设计中，所述无逆if-mme求解器中，关系核矩阵与向量的乘法通过关系核算子实现而无需显式求逆，所述关系核算子用于实现关系核矩阵与任意向量的乘法运算；所述无逆if-mme求解器对任意表型空间向量按如下公式构造右端项：
18、 （3）
19、其中，表示以向量为右端项的右端项向量，表示矩阵或向量的转置，为固定效应设计矩阵，为第个随机效应方向的设计矩阵；
20、当所述关系核由基因型标记矩阵隐式给出时，所述关系核算子对任意向量的作用结果通告如下公式计算：
21、 （12）
22、其中，为关系核矩阵，为中心化并标准化后的标记矩阵，为标准化系数。
23、在一种可能的设计中，根据所述真实残差和随机效应预测值，通过如下公式构造各方差组分方向的工作向量：
24、 （4）
25、其中，为reml投影矩阵，用于消除固定效应影响；为第个方差组分方向对应的导数矩阵，等于
v/
θ
i；
v为表型协方差矩阵，由随机效应协方差和残差协方差共同决定；
θ
i为待估计的方差组分向量；为观测表型向量。
26、在一种可能的设计中，所述伪残差通过如下公式提取：
27、 （5）
28、其中，为第个工作向量，为固定效应设计矩阵，为以为伪表型求解得到的固定效应伪解，为第个随机效应方向的设计矩阵，为以为伪表型求解得到的第个随机效应伪解；
29、所述伪残差与投影矩阵对所述工作向量的作用结果之间满足如下等价关系：
30、 （6）
31、其中，为残差协方差矩阵；
32、当残差协方差矩阵时，即同质残差方差情形下，等价关系表示为：
33、 （7）
34、其中，为残差方差组分。
35、在一种可能的设计中，所述平均信息矩阵元素通过如下公式装配：
36、 （8）
37、当残差协方差矩阵时，即同质残差方差情形下，所述平均信息矩阵元素通过如下公式装配：
38、 （9）
39、所述得分向量的第
i个元素由第一部分和第二部分构成，第一部分为基于真实表型求解得到的工作向量，第二部分为迹项
tr(pm
i
)；所述迹项
tr(pm
i
)通过hutchinson随机迹估计方法获得：生成满足零均值、单位方差条件的随机探针，将所述随机探针作为伪表型输入所述无逆if-mme求解器，根据得到的伪残差通过向量内积估计所述迹项；将所述第二部分从所述第一部分中减去，得到所述得分向量的第
i个元素。
40、在一种可能的设计中，所述方差组分通过如下公式进行更新：
41、 （11）
42、其中，为当前方差组分下的平均信息矩阵，为待估计的方差组分向量，为方差组分更新步长，为reml得分向量，为第轮迭代的当前方差组分，表示迭代轮次，为更新后的方差组分，为阻尼系数，用于控制更新幅度并将方差组分约束在正边界内；
43、当所述混合模型同时包括加性效应和显性效应时，所述混合模型表示为：
44、 （13）
45、 （14）
46、其中，为观测表型向量，为固定效应设计矩阵，为固定效应向量，为加性随机效应设计矩阵，为加性随机效应向量，为显性随机效应设计矩阵，为显性随机效应向量，为残差向量，为加性关系核矩阵，为显性关系核矩阵，为加性遗传方差组分，为显性遗传方差组分，为残差方差组分。
47、第二方面，本发明提供一种动物育种中育种值方差组分估计装置，所述装置包括：
48、数据获取模块，被配置为获取育种群体的观测表型向量、固定效应设计矩阵、随机效应设计矩阵以及基因型标记矩阵；
49、混合模型构建模块，被配置为根据所述观测表型向量、固定效应设计矩阵和随机效应设计矩阵，建立包括固定效应、至少一个随机效应和残差的混合模型，并根据所述基因型标记矩阵构建各方差组分方向的关系核矩阵；
50、求解器构建模块，被配置为根据所述混合模型对应的混合模型方程构建无逆if-mme求解器，所述无逆if-mme求解器通过对所述混合模型方程进行等价变换以规避关系核逆矩阵的显式构建，并以任意表型空间向量作为右端项输入，输出对应的固定效应解、随机效应解和残差；
51、迭代求解模块，被配置为初始化方差组分，并以所述观测表型向量作为右端项输入所述无逆if-mme求解器，根据所述无逆if-mme求解器的输出得到真实残差和随机效应预测值；
52、根据所述真实残差和随机效应预测值构造各方差组分方向的工作向量；
53、将所述工作向量作为伪表型输入所述无逆if-mme求解器，根据所述无逆if-mme求解器的输出得到固定效应伪解和随机效应伪解；
54、根据所述伪表型、固定效应伪解和随机效应伪解提取伪残差，所述伪残差用于表示所述reml投影矩阵对所述工作向量的作用结果；
55、根据所述工作向量和所述伪残差装配平均信息矩阵元素，并通过随机迹估计方法形成得分向量；
56、根据所述平均信息矩阵和所述得分向量求解小规模线性系统，得到方差组分更新步长，并根据所述更新步长更新当前方差组分；
57、重复执行以观测表型向量作为右端项输入无逆if-mme求解器至根据平均信息矩阵和得分向量更新当前方差组分的迭代过程直至满足收敛判据，输出估计得到的加性遗传方差组分、残差方差组分以及由所述加性遗传方差组分和残差方差组分计算的遗传力。
58、第三方面，本发明实施例提供一种电子设备，包括：至少一个处理器和存储器；所述存储器存储计算机执行指令；所述至少一个处理器执行所述存储器存储的计算机执行指令，使得所述至少一个处理器执行如上第一方面以及第一方面各种可能的设计所述的动物育种中育种值方差组分估计方法。
59、第四方面，本发明实施例提供一种计算机可读存储介质，所述计算机可读存储介质中存储有计算机执行指令，当处理器执行所述计算机执行指令时，实现如上第一方面以及第一方面各种可能的设计所述的动物育种中育种值方差组分估计方法。
60、第五方面，本发明实施例提供一种计算机程序产品，包括计算机程序，所述计算机程序被处理器执行时，实现如上第一方面以及第一方面各种可能的设计所述的动物育种中育种值方差组分估计方法。
61、本发明提供的动物育种中育种值方差组分估计方法、装置、设备和介质，其至少具有如下有益效果：
62、第一，本发明通过构建无逆if-mme求解器，将关系核逆矩阵的显式构建替换为关系核算子的矩阵-向量乘运算，并将ai-reml中的投影矩阵作用转化为工作向量伪表型求解后的伪残差提取，避免显式构建关系核逆矩阵、表型协方差逆矩阵和reml投影矩阵，大幅降低内存占用和计算复杂度。
63、第二，本发明复用同一个无逆if-mme求解器完成真实表型求解、工作向量伪表型求解以及随机探针伪表型求解，所述工作向量由真实表型求解结果构造，所述随机探针用于hutchinson随机迹估计中估计得分向量所需的迹项，从而降低重复矩阵构建和缓存压力。
64、第三，本发明的平均信息矩阵通过工作向量与伪残差的内积装配，即按照式（8）或式（9）计算平均信息矩阵元素，适合多右端项批量求解、cpu-gpu分层并行执行和大规模基因型流式计算。
65、第四，本发明可扩展到加性、显性、多核和多性状方差组分估计，只需为新增方差组分提供相应关系核算子和工作向量，即可沿用同一伪残差提取和平均信息矩阵装配机制。
66、第五，本发明通过阻尼系数α控制方差组分更新幅度并将方差组分约束在正边界内，确保迭代过程的数值稳定性。
67、第六，本发明与传统reml方法进行基准测试提供可复现实验接口，能够记录方差组分估计值、遗传力估计值、运行时间、峰值内存和收敛状态等指标。在10000个体、49656个snp的基准数据中，if-aireml在低遗传力和高遗传力场景下的总耗时分别为54.17s和57.49s，峰值内存分别为1.44gb和1.43gb；真实数据隐式rmom+if-aireml主路径在85362个体数据上收敛得到h2=0.34840368，总耗时8894.16s，峰值内存46.69gb。