CN107563147B - 一种估计基因组育种值的方法及装置 - Google Patents
一种估计基因组育种值的方法及装置 Download PDFInfo
- Publication number
- CN107563147B CN107563147B CN201710651753.2A CN201710651753A CN107563147B CN 107563147 B CN107563147 B CN 107563147B CN 201710651753 A CN201710651753 A CN 201710651753A CN 107563147 B CN107563147 B CN 107563147B
- Authority
- CN
- China
- Prior art keywords
- mpi
- vector
- matrix
- determining
- breeding value
- 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.)
- Active
Links
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明实施例提供一种估计基因组育种值的方法及装置，该方法包括：根据矩阵M‑1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的根据表型数据文件的记录条数、矩阵A‑1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值。本发明实施例可以实现计算机资源与运算速度的最优化，节约运算时间，提高工作效率。
Description
技术领域
本发明实施例涉及生物信息学技术领域，尤其涉及一种估计基因组育种值的方法及装置。
背景技术
在畜禽育种领域，数量遗传学理论的发展推动着育种方法技术的进步，而新的育种技术的实施与应用则迫切需要得到相应的分子生物学技术与计算机技术的支持。
2008年美国首先利用覆盖全基因组的标记信息对个体进行遗传评估。在此之前，用于个体遗传评估的信息主要包括系谱和个体表型记录。在实际育种工作中，无论使用何种个体遗传评估信息，一般使用最佳线性无偏预测法(best linear unbiasedprediction,BLUP)预测个体的育种值，继而选出遗传上优良的个体，用于繁育后代。Meuwissen等人(2001)提出利用全基因组标记信息进行基因组选择(genomic selection,GS)的方法。研究也表明，在奶牛育种中，实施基因组选择可比后裔测定方法节约92％的育种成本(Schaeffer,2006)。近年来，随着SNP芯片价格的不断下降，越来越多的遗传评估采用基因组选择进行个体遗传评估，与此同时，基因组选择中所包含的测定了SNP基因型信息的个体数也不断增加。例如，目前美国测定了SNP基因型的荷斯坦牛个体已达160多万头。基因组信息的激增使求解由BLUP方法构建的混合模型方程组(Mixed Model Equations,MME)所耗费的时间大大增加。
目前求解混合模型方程组的方法主要为迭代法，随着基因组信息的不断增加，混合模型方程组也越来越大，越来越复杂；相应地，每次迭代的耗时会相应增长。现有技术在利用迭代法计算基因组育种值时，采用的是串行计算方法，即依次顺序读取相应的待运算数据，计算效率较低，同时也造成计算机资源的浪费。
因此，亟需一种估计基因组育种值的方法，能够节约运算时间，充分利用计算机资源，提高工作效率。
发明内容
为解决现有的基因组育种值估计方法中存在的计算耗时过长、计算效率低的问题，本发明实施例提供一种估计基因组育种值的方法及装置。
第一方面，本发明实施例提供一种估计基因组育种值的方法，该方法包括：
根据矩阵M-1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的
根据表型数据文件的记录条数、矩阵A-1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和
根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值；
其中，所述矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，Wi为MPI进程MPIi读取的表型数据文件对应的系数矩阵，A-1为个体间加性遗传相关矩阵A的逆矩阵，n为大于1的整数，i＝1，2，3，…，n。
第二方面，本发明实施例提供一种估计基因组育种值的装置，该装置包括：
第一分进程单元，根据矩阵M-1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的
第二分进程单元，根据表型数据文件的记录条数、矩阵A-1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和
估计单元，根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值；
其中，所述矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，Wi为MPI进程MPIi读取的表型数据文件对应的系数矩阵，A-1为个体间加性遗传相关矩阵A的逆矩阵，n为大于1的整数，i＝1，2，3，…，n。
第三方面，本发明实施例还提供一种电子设备，该电子设备包括：存储器和处理器，所述处理器和所述存储器通过总线完成相互间的通信；所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行如下方法：根据矩阵M-1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的根据表型数据文件的记录条数、矩阵A-1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值；其中，所述矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，Wi为MPI进程MPIi读取的表型数据文件对应的系数矩阵，A-1为个体间加性遗传相关矩阵A的逆矩阵，n为大于1的整数，i＝1，2，3，…，n。
第四方面，本发明实施例还提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现如下方法：根据矩阵M-1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的根据表型数据文件的记录条数、矩阵A-1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值；其中，所述矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，Wi为MPI进程MPIi读取的表型数据文件对应的系数矩阵，A-1为个体间加性遗传相关矩阵A的逆矩阵，n为大于1的整数，i＝1，2，3，…，n。
本发明实施例提供的估计基因组育种值的方法和装置，通过将表型数据文件、M-1矩阵以及A-1矩阵进行分割并分配给多个MPI进程，结合MKL函数，实现并行化计算基因组育种值。在实际应用中，可以根据计算机使用情况设置MPI的进程数与MKL的函数运算时所用的线程数，控制并行计算所用资源，从而实现计算机资源与运算速度的最优化，节约了运算时间，提高了工作效率。
附图说明
图1为本发明实施例提供的估计基因组育种值的方法的流程示意图；
图2为本发明实施例提供的估计基因组育种值的算法流程示意图；
图3为本发明实施例提供的估计基因组育种值的装置的结构示意图；
图4为本发明实施例提供的电子设备的结构示意框图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
现代计算机，尤其是常用于执行运算任务的服务器，多是共享内存结构，且多具有多处理器多核的特点。Math Kernel Library(MKL)是广泛应用于科学、工程和金融领域的、优化的数学函数库。MKL中包含线性代数函数(向量-向量运算、向量-矩阵运算和矩阵-矩阵运算)、统计函数和数据拟合函数等。同时，其中很多函数都可以进行并行运算，加快运算速度。在MPI标准中，每个MPI进程都有各自的地址空间。在MPI-3标准中，每个MPI进程分配并维持自己的本地内存，并在共享内存区域显示部分内存。所有进程均可访问共享内存区域，因此可以减少进程中的数据交换。
本发明实施例在进行基因组育种值估计时，基于间接迭代技术的预条件共轭梯度迭代法的算法，利用MKL和MPI编程并行化处理待运算数据，以充分利用现代计算机多处理器多核的特点来加快混合模型方程组的求解。需要说明的是，本发明实施例中基因组育种值的估计还适用于利用系谱进行育种值估计。
本发明实施例中，通过根据所有个体遗传信息及表型数据文件，利用BLUP方法构建混合模型方程组，假设该混合模型方程组为Cx＝b，其中C为系数矩阵，b为右手项向量。使用预条件共轭梯度法求解x时，涉及四个向量，每次迭代对这四个向量进行更新。这四个向量分别为残差向量r，搜索方向向量d，解向量x和临时向量v。估计基因组育种值的过程也就是求解解向量x的过程。
图1为本发明实施例提供的估计基因组育种值的方法的流程示意图。如图1所示，该方法包括以下步骤：
步骤S101、根据矩阵M-1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的
具体地，矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，存储于二进制文件中，在进行基因组育种值估计时可以从所述二进制文件中直接读取。可选地，变换矩阵M为分块对角阵，为方便后续矩阵与向量的乘法运算，二进制文件以子块为单位进行存储，依次储存以下元素：子块的首行在M-1中的行号、子块的行数(子块为方阵，其行数等于列数)。当读取子块时，子块中的所有矩阵元素按行输出。使用辅助文件记录M-1中子块的个数及储存每个子块所用的字节数。
根据M-1矩阵的子块数和预设的MPI进程数n，将M-1按行切割为若n恰好被M-1矩阵的子块数整除，则每个MPI进程分到相同子块数的M-1矩阵元素；若n不能被M-1矩阵的子块数整除，则其熵所指代的剩余子块数的矩阵元素则全部分配给第一进程，所谓第一进程是指分到M-1首行矩阵元素的进程，即MPI1。
例如，若可以看到M-1有两个子块，二进制文件中的存储方式为：
第一行：1 2 1 4 2 1；
第二行：3 2 1 2 3 1。
从上述存储内容可知，第一个分块阵的首行在M-1中的行号为1，该分块阵的行数为2，矩阵元素为1 4 2 1；第二个分块阵的首行在M-1中的行号为3，该分块阵的行数也为2，矩阵元素为1 2 3 1。
若有2个MPI进程，则这2个进程在读取二进制文件中的M-1时，均匀分配的结果是：
进程MPI1分到为第一行：1 4 2 1；
进程MPI2分到的为第二行：1 2 3 1。
若可以看到M-1有3个子块，二进制文件中的存储方式为：
第一行：1 3 1 2 3 3 2 7 4 2 6
第二行：4 3 1 3 2 9 4 0 1 2 3
第三行：7 3 1 4 7 5 3 2 6 5 2
若有2个MPI进程，则这2个进程在读取二进制文件中的M-1时，均匀分配的结果是：
进程1分到的为第一行和第三行：1 2 3 3 2 7 4 2 6；1 4 7 5 3 2 65 2；
进程2分到的为第二行：1 3 2 9 4 0 1 2 3。
可以理解地是，上述举例仅为说明确定MPI进程MPIi读取的过程，在实际应用中，M-1矩阵一般为维数等于系数矩阵维数的对称阵。
步骤S102、根据待估计基因组的表型数据文件的记录条数、矩阵A-1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和
具体地，每条表型数据记录都对应系数矩阵C中的一个向量wi；将表型数据文件按记录条数和进程数均匀分割，各MPI进程同时读入分配到的表型数据文件，并计算读入记录所对应的wi，进而确定相应的系数矩阵Wi；A为个体间加性遗传相关矩阵，A可以利用系谱构建，也可以利用基因型信息构建，也可以同时利用系谱和基因型信息构建。类似矩阵M-1，A-1存储于二进制文件中，对于每个个体，在文件中依次写入：个体编号(从1开始连续编号)、与其他个体间非零的加性遗传相关系数的总个数、这些个体的编号、相应的遗传相关系数等数据，确定MPIi读取的的方法同相似，此处不再赘述。
步骤S103、根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值；
其中，所述矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，Wi为MPI进程MPIi读取的表型数据文件对应的系数矩阵，A-1为个体间加性遗传相关矩阵A的逆矩阵，n为大于1的整数，i＝1，2，3，…，n。
具体地，进程MPIi根据系数矩阵Wi和实现矩阵与向量的相差或向量与向量相乘的并行化计算，按照现有的算法流程，确定所述待估计基因组的基因组育种值。
本发明实施例提供的估计基因组育种值的方法，通过将大型分块矩阵均分为多个部分，从而利用计算机多进程，实现矩阵与向量，向量与向量之间的并行计算，可以充分利用计算机资源，加快计算速度，提高了工作效率。
在上述实施例的基础上，该方法中的根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值，具体包括：
根据MPI进程MPIi的与向量r的计算结果，确定所述基因组育种值的搜索方向向量d，其中，所述向量r为预设残差向量；
根据所述搜索方向向量d、MPI进程MPIi的系数矩阵Wi和利用间接迭代法，确定临时向量v；
根据所述临时向量v，确定所述基因组育种值。
具体地，右手项向量b是已知的(可以从二进制文件中直接读取得到)，预设残差向量r的初始值即为右手向量b。在计算过程中，根据具体情况，预设残差向量会相应地变化，具体会在后面具体说明，此处不再赘述。进程MPIi计算与向量r的乘积，各进程的计算结果，各进程将计算结果存储到共享内存中搜索方向向量d的对应部分。
MPI进程MPIi根据系数矩阵Wi和利用间接迭代法，分别计算系数矩阵Wi和搜索方向向量d的乘积，计算和搜索方向向量d的乘积，从而确定估计基因组育种值过程中的临时向量v，根据临时向量v，确定所述待估计的基因组育种值。
本发明实施例提供的估计基因组育种值的方法，通过利用计算机多进程，实现矩阵与向量，向量与向量乘积运算的并行计算，可以充分利用计算机资源，加快计算速度，提高了工作效率。
在上述各实施例的基础上，该方法中的根据所述搜索方向向量d、MPI进程MPIi的系数矩阵Wi和利用间接迭代法，确定所述临时向量v，具体包括：
根据MPI进程MPIi的系数矩阵Wi和所述搜索方向向量d的迭代计算结果，确定所述临时方向向量v的第一加和项；
根据MPI进程MPIi的和所述搜索方向向量d的迭代计算结果，确定所述临时方向向量v的第二加和项；
根据所述第一加和项和所述第二加和项，确定所述临时方向向量v。
具体地，利用间接迭代法进行计算的步骤如下：
其中，Nq为进程MPIi所分配的表型数据文件的记录条数，wj为所述表型数据文件中第j条数据所对应的向量，Rj为第j条数据所对应的残差方差-协方差矩阵，G0为加性遗传方差-协方差矩阵。根据MPI进程MPIi的系数矩阵Wi和所述搜索方向向量d的迭代计算结果，确定所述临时方向向量v的第一加和项根据MPI进程MPIi的和所述搜索方向向量d的迭代计算结果，确定所述临时方向向量v的第二加和项vd；之后，根据所述第一加和项和所述第二加和项，确定所述临时方向向量v。
本发明实施例提供的估计基因组育种值的方法，通过利用计算机多进程，实现矩阵与向量，向量与向量之间的并行计算，可以充分利用计算机资源，加快计算速度，提高了工作效率。
在上述各实施例的基础上，该方法中，所述MPI进程MPIi采用MKL函数进行相应的向量运算或矩阵运算。
图2为本发明实施例提供的估计基因组育种值的算法流程示意图。如图2所示，num指MPI进程数、conv为收敛阈值(一般为10-13或更小)、maxr为最大迭代次数，在利用计算机估计基因组育种值时，可以根据计算机使用情况设置MPI的进程数与MKL函数运算时所用的线程数。
步骤S201、在读入右手项向量b、预设的MPI进程数、收敛阈值、最大迭代次数的同时，为所有MPI进程创建共享内存，在该共享内存中存储r、d、v、b、x、c、q、α、β；
步骤S202、解向量x的初值设为0，残差向量r的初始值为向量b，c为用于指示收敛阈值的参数，q为用于指示迭代次数的参数，这两个参数的初始值均设为1。利用MKL的cblas_ddot计算||b||2。该步骤计算简单，可以使用多进程中的其中一个进程执行即可。
步骤S203、该步骤即为将矩阵M-1分为多个部分，所有进程同时读取M-1中的分块阵中的子块，然后利用MKL的cblas_dsymv函数计算读入的子块与向量r的乘积，r向量可从共享内存中读取，按照读取所在的行与向量d的对应关系，将所有进程的计算结果存储至共享内存中的向量d中；
步骤S204、利用MKL的cblas_ddot计算rTd。该步骤计算简单，可以使用多进程中的其中一个进程执行，之后将计算结果存入共享内存即可；
步骤S205、该步骤是对收敛阈值和迭代次数进行判断是否满足阈值条件。当收敛阈值c小于预设收敛阈值，或当前迭代次数q大于预设的最大迭代次数时，即可执行步骤S220：输出解向量；若不满足阈值条件，则进行下一次的迭代计算。该判断过程所有进程均需进行判断；
步骤S206、该步骤即将表型数据文件及A-1分为多个部分，分配给多个进程，实现并行化计算，具体可参考上述实施例，此处不再赘述；
步骤S207、该步骤利用MKL的cblas_ddot计算dTv的内积，使用多进程中的其中一个进程执行，之后将计算结果存入共享内存即可；
步骤S208、计算解向量x在搜索方向向量d上的强度α，使用多进程中的其中一个进程执行，之后将计算结果存入共享内存即可；
步骤S209、该步骤即为解向量x的更新过程。将上一次迭代计算得到的解向量x按行均分给各进程，将本次更新的αd向量按照同样的方式均分给各进程，所有进程同时进行加和运算。即，各进程均分x和αd向量长度的加和运算；
步骤S210、该步骤是为了校正残差向量，此处的判断条件q不限于被100整除，也可以被50整除，120整除，或者150整除等等。该判断过程所有进程均需进行判断；若q没有被100整除，则执行步骤S211；若q被100整除，则执行步骤S212；
步骤S211、r＝r-αv，类似步骤S209，所有进程均分向量r和向量αv向量长度的减法运算；该计算结果为r的近似值；
步骤S212、r＝b-Cx；Cx的计算过程类似步骤S206，此处不再赘述；之后，所有进程均分向量b和向量Cx向量长度的减法运算；该计算结果为r的准确值；
步骤S213、该步骤是对临时向量v进行更新，该步骤的计算方式同步骤S203，此处不再赘述；
步骤S214、该步骤类似步骤S207，利用MKL的cblas_ddot计算rTv的内积，使用多进程中的其中一个进程执行，之后将计算结果存入共享内存即可；
步骤S215、该步骤类似步骤S208，使用多进程中的其中一个进程执行，之后将计算结果存入共享内存即可；
步骤S216、更新f，使用多进程中的其中一个进程执行即可；
步骤S217、该步骤类似步骤S209，所有进程均分向量v和向量βd向量长度的加和运算；
步骤S218、该步骤类似步骤S207，利用MKL的cblas_ddot计算rTr的内积，使用多进程中的其中一个进程执行，之后将计算结果存入共享内存即可；
步骤S219，该步骤是在进行一次迭代计算之后，对c和q进行更新，使用多进程中的其中一个进程执行，之后将计算结果存入共享内存即可；
步骤S210，输出解向量x。
在实际应用中，利用上述计算流程，下面通过例子来具体说明。
实例1
某地区收集了8,085头奶牛的85,188条前三胎的产奶性状的数据，包括产奶量、乳脂量和乳蛋白量。系谱中共涉及19,975个个体。
若育种值估计模型为用随机回归模型进行育种值估计，模型如下：
yrsijkl＝HTDrsi+x(t)′rsklbrsj+z(t)′rsklursk+w(t)′rsklprsk+ersijkl,
其中，yrsijkl为场-测定日效应i，产犊月龄-季节j下，母牛k在第s胎的测定日记录l中性状r的观察值(前三胎产奶量、乳脂量和乳蛋白量的测定日记录)，HTDrsi为性状r第s胎场-测定日i效应(固定效应)，brsj为性状r第s胎产犊月龄-季节j效应下嵌套的固定回归系数向量(3-5月为第一季，6-8月为第二季，9-11为第三季，12-次年2月为第四季)，ursk为个体k性状r第s胎的加性遗传效应下嵌套的随机回归系数向量，prsk为个体k性状r第s胎的永久环境效应下嵌套的随机回归系数向量，ersijkl为残差。x(t)和z(t)使用四阶Legendre多项式，w(t)使用三阶Legendre多项式。
相应的矩阵形式为：
y＝Xb+Za+Wp+e，
式中y是表型观察值向量，b是场测定日固定效应和产犊月龄-季节的固定回归系数向量，a和p是加性遗传和永久环境效应的随机回归系数向量，X，Z和W分别是对应于b，a和p的设计矩阵，e是残差向量。其中，假设：
其中，A为根据系谱构建的分子亲缘相关矩阵，I为维数等于永久环境效应水平数的单位阵，C和P分别为加性效应和永久环境效应随机回归系数的方差-协方差矩阵。残差方差异质，即不同泌乳天数区间的残差方差-协方差不同，泌乳天数区间划分为：5-45d，46-155d，156-205d，206-305d。在一个泌乳天数区间内，R的表达式如下：
其中，R11，R22和R33分别对应于1-3胎胎次内的残差方差-协方差矩阵。
此时，相应的混合模型方程组为
相应于方法中示例的Cx＝b，显然
且
使用的变换矩阵的逆矩阵M-1为C阵中以个体进行分块的分块对角阵。
表1为使用不同进程数求解实例1的混合模型方程组所耗的时长表。从表1中可以看出，当进程数为2时，比使用单进程节约了36％的时间；当进程数为5时，比使用单进程节约了57％的时间。可以看出，本发明实施例所提出的方法对于加快求解混合模型方程组，节约运算时间是非常有效的。
表1使用不同进程数求解实例1的混合模型方程组所耗的时长
| 进程数 | 1 | 2 | 5 |
| 时间(s) | 244 | 155 | 105 |
实例2
现有1,000个体，每个个体均有逆回归育种值，同时，均有SNP基因型信息，SNP位点数为20,010。
若育种值估计模型为模型用GBLUP模型进行个体育种值估计，模型如下：
yi＝μ+ai+ei，
其中，yi为个体i的逆回归育种值，μ为群体均值，ai为个体的加性遗传效应，ei为残差。
相应的矩阵形式为：
y＝1μ+Za+e，
其中，y为逆回归育种值向量，μ为群体均值，a为育种值向量，Z为相应的设计矩阵，e为残差向量。其中，假设
其中，G为根据SNP标记构建的个体关系矩阵，为加性遗传方差，为残差方差。
此时，相应的混合模型方程组为
相应于方法中示例的Cx＝b，显然
且
使用的变换矩阵的逆矩阵M-1为C阵中以个体进行分块的分块对角阵。
表2为使用不同进程数求解实例2的混合模型方程组所耗的时长表。从表2中可以看出，当进程数为2时，比使用单进程节约了33％的时间；当进程数为5时，比使用单进程节约了67％的时间。
表2使用不同进程数求解实例2的混合模型方程组所耗的时长
| 进程数 | 1 | 2 | 5 |
| 时间(s) | 12 | 8 | 4 |
实例3
某地区收集了126,157头奶牛1,747,843条前三胎的产奶性状的数据，包括产奶量、乳脂量和乳蛋白量。系谱中共涉及272,172个个体。同时，其中7,556个个体具有SNP芯片信息，经过质控后，可利用的SNP位点数为72,507。
若育种值估计模型为用基于随机回归模型的基因组选择“一步”法进行育种值估计，模型如下：
yrsijkl＝HTDrsi+x(t)′rsklbrsj+z(t)′rsklursk+w(t)′rsklprsk+ersijkl，
其中，yrsijkl为场-测定日效应i，产犊月龄-季节j下，母牛k在第s胎的测定日记录l中性状r的观察值(前三胎产奶量、乳脂量和乳蛋白量的测定日记录)，HTDrsi为性状r第s胎场-测定日i效应(固定效应)，brsj为性状r第s胎产犊月龄-季节j效应下嵌套的固定回归系数向量(3-5月为第一季，6-8月为第二季，9-11为第三季，12-次年2月为第四季)，ursk为个体k性状r第s胎的加性遗传效应下嵌套的随机回归系数向量，prsk为个体k性状r第s胎的永久环境效应下嵌套的随机回归系数向量，ersijkl为残差。x(t)和z(t)用四Legendre多项式，w(t)用三阶Legendre多项式。
相应的矩阵形式为：
y＝Xb+Za+Wp+e，
式中y是表型观察值向量，b是场测定日固定效应和产犊月龄-季节的固定回归系数向量，a和p是加性遗传和永久环境效应的随机回归系数向量，X，Z和W分别是对应于b，a和p的设计矩阵，e是残差向量。其中，假设：
其中，I为维数等于永久环境效应水平数的单位阵，C和P分别为加性效应和永久环境效应随机回归系数的方差-协方差矩阵。残差方差异质，即不同泌乳天数区间的残差方差-协方差不同，泌乳天数区间划分为：5-45d，46-155d，156-205d，206-305d。在一个泌乳天数区间内，R的表达式如下：
其中，R11，R22和R33分别对应于1-3胎胎次内的残差方差-协方差矩阵。
H为个体关系矩阵，通过结合系谱与基因型信息构建，其逆矩阵的表达式为：
因此，相应的混合模型方程组为
相应于方法中示例的Cx＝b，显然
且
使用的变换矩阵的逆矩阵M-1为C阵中以个体进行分块的分块对角阵。
表3为使用不同进程数求解实例3的混合模型方程组所耗的时长表。从表中可以看出，当进程数为5时，比使用单进程节约了28％的时间；当进程数为10时，比使用单进程节约了56％的时间。
表3使用不同进程数求解实例3中混合模型方程组所耗的时长
| 进程数 | 1 | 5 | 10 |
| 时间(h) | 7.2 | 5.2 | 3.2 |
本发明实施例提供的估计基因组育种值的方法，将间接迭代法的预条件共轭梯度迭代计算基因组育种值的流程与MPI的并行编程相结合，统筹规划实现并行计算，从而实现计算机资源与运算速度的最优化。
图3为本发明实施例提供的估计基因组育种值的装置的结构示意图，如图3所示，该装置包括：第一分进程单元301、第二分进程单元302及估计单元303。其中，第一分进程单元301用于根据矩阵M-1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的第二分进程单元用于根据表型数据文件的记录条数、矩阵A-1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和估计单元303用于根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值；其中，所述矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，Wi为MPI进程MPIi读取的表型数据文件对应的系数矩阵，A-1为个体间加性遗传相关矩阵A的逆矩阵，n为大于1的整数，i＝1，2，3，…，n。
具体地，第一分进程单元301将矩阵M-1分配给多个MPI进程，第二分进程单元302将表型数据文件和矩阵A-1分配给多个MPI进程，估计单元303根据进程MPIi的系数矩阵Wi和进行基因组育种值的估计。本发明实施例提供的装置，是为了实现上述方法的，其功能具体可参考上述方法实施例，此处不再赘述。
在上实施例的基础上，该装置中的估计单元，具体包括：第一子估计单元，第二子估计单元和第三子估计单元。其中，第一子估计单元用于根据MPI进程MPIi的与向量r的计算结果，确定所述基因组育种值的搜索方向向量d，其中，所述向量r为预设残差向量；第二子估计单元用于根据所述搜索方向向量d、MPI进程MPIi的系数矩阵Wi和利用间接迭代法，确定临时向量v；第三子估计单元用于根据所述临时向量v，确定所述基因组育种值。
具体地，第一子估计单元根据与向量r，确定搜索方向向量d；第二子估计单元根据搜索方向向量d、Wi和确定临时向量v；第三子估计单元根据该向量v确定基因组育种值。本发明实施例提供的装置，是为了实现上述方法的，其功能具体可参考上述方法实施例，此处不再赘述。
在上述实施例的基础上，该装置中的第二子估计单元具体用于根据MPI进程MPIi的系数矩阵Wi和所述搜索方向向量d的迭代计算结果，确定所述临时向量v的第一加和项；根据MPI进程MPIi的和所述搜索方向向量d的迭代计算结果，确定所述临时方向向量v的第二加和项；根据所述第一加和项和所述第二加和项，确定所述临时向量v。
可选地，该装置中的估计单元在进行计算时，MPI进程MPIi采用MKL函数进行相应地向量运算或矩阵运算。
图4为本发明实施例提供的电子设备的结构示意框图，如图4所示，参照图4，所述电子设备，包括：处理器(processor)401、存储器(memory)402和总线403；其中，所述处理器401和所述存储器402通过所述总线403完成相互间的通信；所述处理器401用于调用所述存储器402中的程序指令，以执行上述各方法实施例所提供的方法，例如包括：根据矩阵M-1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的根据表型数据文件的记录条数、矩阵A-1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值；其中，所述矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，Wi为MPI进程MPIi读取的表型数据文件对应的系数矩阵，A-1为个体间加性遗传相关矩阵A的逆矩阵，n为大于1的整数，i＝1，2，3，…，n。
本发明实施例还提供一种计算机程序产品，所述计算机程序产品包括存储在非暂态计算机可读存储介质上的计算机程序，所述计算机程序包括程序指令，当所述程序指令被计算机执行时，计算机能够执行上述各方法实施例所提供的方法，例如包括：根据矩阵M-1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的根据表型数据文件的记录条数、矩阵A-1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值；其中，所述矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，Wi为MPI进程MPIi读取的表型数据文件对应的系数矩阵，A-1为个体间加性遗传相关矩阵A的逆矩阵，n为大于1的整数，i＝1，2，3，…，n。
本发明实施例还提供一种计算机可读存储介质，所述计算机可读存储介质存储计算机指令，所述计算机指令使所述计算机执行上述各方法实施例所提供的方法，例如包括：根据矩阵M-1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的根据表型数据文件的记录条数、矩阵A-1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值；其中，所述矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，Wi为MPI进程MPIi读取的表型数据文件对应的系数矩阵，A-1为个体间加性遗传相关矩阵A的逆矩阵，n为大于1的整数，i＝1，2，3，…，n。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
Claims (8)
1.一种估计基因组育种值的方法，其特征在于，包括：
根据矩阵M-1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的
根据表型数据文件的记录条数、矩阵A-1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和
根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值；
其中，所述矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，Wi为MPI进程MPIi读取的表型数据文件对应的系数矩阵，A-1为个体间加性遗传相关矩阵A的逆矩阵，n为大于1的整数，i＝1，2，3，…，n；
所述根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值，具体包括：
根据MPI进程MPIi的与向量r的计算结果，确定所述基因组育种值的搜索方向向量d，其中，所述向量r为预设残差向量；
根据所述搜索方向向量d、MPI进程MPIi的系数矩阵Wi和利用间接迭代法，确定临时向量v；
根据所述临时向量v，确定所述基因组育种值。
2.根据权利要求1所述的方法，所述根据所述搜索方向向量d、MPI进程MPIi的系数矩阵Wi和利用间接迭代法，确定临时向量v，具体包括：
根据MPI进程MPIi的系数矩阵Wi和所述搜索方向向量d的迭代计算结果，确定所述临时向量v的第一加和项；
根据MPI进程MPIi的和所述搜索方向向量d的迭代计算结果，确定所述临时向量v的第二加和项；
根据所述第一加和项和所述第二加和项，确定所述临时向量v。
3.根据权利要求1或2所述的方法，其特征在于，所述MPI进程MPIi利用MKL函数进行相应的向量运算或矩阵运算。
4.一种估计基因组育种值的装置，其特征在于，包括：
第一分进程单元，根据矩阵M-1的子块数和预设的MPI进程数n，确定MPI进程MPIi读取的
第二分进程单元，根据表型数据文件的记录条数、矩阵A-1中的个体数和所述MPI进程数n，确定MPI进程MPIi读取的系数矩阵Wi和
估计单元，根据MPI进程MPIi的系数矩阵Wi和确定待估计的基因组育种值；
其中，所述矩阵M-1为用于加速所述基因组育种值估计的变换矩阵M的逆矩阵，Wi为MPI进程MPIi读取的表型数据文件对应的系数矩阵，A-1为个体间加性遗传相关矩阵A的逆矩阵，n为大于1的整数，i＝1，2，3，…，n；
所述估计单元，具体包括：
第一子估计单元，根据MPI进程MPIi的与向量r的计算结果，确定所述基因组育种值的搜索方向向量d，其中，所述向量r为预设残差向量；
第二子估计单元，根据所述搜索方向向量d、MPI进程MPIi的系数矩阵Wi和利用间接迭代法，确定临时向量v；
第三子估计单元，根据所述临时向量v，确定所述基因组育种值。
5.根据权利要求4所述的装置，所述第二子估计单元，具体包括：
根据MPI进程MPIi的系数矩阵Wi和所述搜索方向向量d的迭代计算结果，确定所述临时向量v的第一加和项；
根据MPI进程MPIi的和所述搜索方向向量d的迭代计算结果，确定所述临时向量v的第二加和项；
根据所述第一加和项和所述第二加和项，确定所述临时向量v。
6.根据权利要求4或5所述的装置，其特征在于，所述MPI进程MPIi利用MKL函数进行相应地向量运算或矩阵运算。
7.一种电子设备，其特征在于，包括：
存储器和处理器，所述处理器和所述存储器通过总线完成相互间的通信；所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行如权利要求1至3任一所述的方法。
8.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1至3任一所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710651753.2A CN107563147B (zh) | 2017-08-02 | 2017-08-02 | 一种估计基因组育种值的方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710651753.2A CN107563147B (zh) | 2017-08-02 | 2017-08-02 | 一种估计基因组育种值的方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107563147A CN107563147A (zh) | 2018-01-09 |
| CN107563147B true CN107563147B (zh) | 2019-12-20 |
Family
ID=60975090
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710651753.2A Active CN107563147B (zh) | 2017-08-02 | 2017-08-02 | 一种估计基因组育种值的方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107563147B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108157293B (zh) * | 2018-02-07 | 2020-03-20 | 山东省农业科学院奶牛研究中心 | 一种基于系谱信息简化选择高生产性能a2a2纯合基因型奶牛的育种方法 |
| CN110211635B (zh) * | 2019-06-12 | 2020-07-21 | 北京康普森农业科技有限公司 | 用于畜禽基因组选择分析的方法及畜禽育种方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101539737B1 (ko) * | 2014-04-18 | 2015-07-28 | 주식회사 씨더스 | 유전체 정보와 분자마커를 이용한 여교잡 선발의 효율성 증진 기술 |
| CN104866904A (zh) * | 2015-06-16 | 2015-08-26 | 中电科软件信息服务有限公司 | 一种基于spark的遗传算法优化的BP神经网络并行化方法 |
| CN104965997A (zh) * | 2015-06-05 | 2015-10-07 | 浙江工业大学 | 一种基于植物功能与结构模型的作物虚拟育种方法 |
| CN105087571A (zh) * | 2015-09-09 | 2015-11-25 | 中国农业大学 | 一种筛查荷斯坦奶牛脊椎畸形综合征携带者的分子检测方法 |
| CN106779076A (zh) * | 2016-11-18 | 2017-05-31 | 栾图 | 基于生物信息的选育良种系统及其算法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20170091384A1 (en) * | 2015-09-28 | 2017-03-30 | International Business Machines Corporation | Systems and methods for fitting ld distributions at genomic scales |
-
2017
- 2017-08-02 CN CN201710651753.2A patent/CN107563147B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101539737B1 (ko) * | 2014-04-18 | 2015-07-28 | 주식회사 씨더스 | 유전체 정보와 분자마커를 이용한 여교잡 선발의 효율성 증진 기술 |
| CN104965997A (zh) * | 2015-06-05 | 2015-10-07 | 浙江工业大学 | 一种基于植物功能与结构模型的作物虚拟育种方法 |
| CN104866904A (zh) * | 2015-06-16 | 2015-08-26 | 中电科软件信息服务有限公司 | 一种基于spark的遗传算法优化的BP神经网络并行化方法 |
| CN105087571A (zh) * | 2015-09-09 | 2015-11-25 | 中国农业大学 | 一种筛查荷斯坦奶牛脊椎畸形综合征携带者的分子检测方法 |
| CN106779076A (zh) * | 2016-11-18 | 2017-05-31 | 栾图 | 基于生物信息的选育良种系统及其算法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107563147A (zh) | 2018-01-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Fernando et al. | A class of Bayesian methods to combine large numbers of genotyped and non-genotyped animals for whole-genome analyses | |
| US20180260709A1 (en) | Calculating device and method for a sparsely connected artificial neural network | |
| Jonas et al. | Does genomic selection have a future in plant breeding? | |
| Cleveland et al. | Practical implementation of cost-effective genomic selection in commercial pig breeding using imputation | |
| Hayashi et al. | A Bayesian method and its variational approximation for prediction of genomic breeding values in multiple traits | |
| Cavani et al. | Estimates of genetic parameters for reproductive traits in Brahman cattle breed | |
| CN107563147B (zh) | 一种估计基因组育种值的方法及装置 | |
| Masuda et al. | Application of supernodal sparse factorization and inversion to the estimation of (co) variance components by residual maximum likelihood | |
| Ning et al. | Efficient multivariate analysis algorithms for longitudinal genome-wide association studies | |
| Vandenplas et al. | Deflated preconditioned conjugate gradient method for solving single-step BLUP models efficiently | |
| Prakapenka et al. | GVCHAP: a computing pipeline for genomic prediction and variance component estimation using haplotypes and SNP markers | |
| Koivula et al. | Practical implementation of genetic groups in single-step genomic evaluations with Woodbury matrix identity–based genomic relationship inverse | |
| Strandén et al. | Bpop: an efficient program for estimating base population allele frequencies in single and multiple group structured populations | |
| Wu et al. | Parallel Markov chain Monte Carlo-bridging the gap to high-performance Bayesian computation in animal breeding and genetics | |
| Amuzu-Aweh et al. | Prediction of heterosis using genome-wide SNP-marker data: application to egg production traits in white Leghorn crosses | |
| Gonzalez-Peña et al. | Impact of pig insemination technique and semen preparation on profitability | |
| Tyrisevä et al. | Principal component approach in variance component estimation for international sire evaluation | |
| US11726757B2 (en) | Processor for performing dynamic programming according to an instruction, and a method for configuring a processor for dynamic programming via an instruction | |
| Faux et al. | A recursive algorithm for decomposition and creation of the inverse of the genomic relationship matrix | |
| Passafaro et al. | Would large dataset sample size unveil the potential of deep neural networks for improved genome-enabled prediction of complex traits? The case for body weight in broilers | |
| Silveira et al. | Regression trees in genomic selection for carcass traits in pigs | |
| Pitkänen et al. | From data to genomic breeding values with the MiX99 software suite | |
| CN115938477A (zh) | 多性状育种值的测定方法、装置、设备及存储介质 | |
| Freudenberg et al. | Accelerated matrix-vector multiplications for matrices involving genotype covariates with applications in genomic prediction | |
| Jamrozik et al. | Comparison of two computing algorithms for solving mixed model equations for multiple trait random regression test day models |
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 | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20210324 Address after: 100080 room 1703, 17 / F, 27 Zhongguancun Street, Haidian District, Beijing Patentee after: BEIJING NONGXIN INTERNET DATA TECHNOLOGY Co.,Ltd. Address before: 100193 No. 2 Old Summer Palace West Road, Beijing, Haidian District Patentee before: CHINA AGRICULTURAL University |
|
| TR01 | Transfer of patent right |