CN106779076A - 基于生物信息的选育良种系统及其算法 - Google Patents

Abstract

本发明涉及生物技术领域，具体公开了一种基于生物信息的选育良种系统及其算法，整合了全基因组选择中个体间的连锁不平衡和连锁分析关系信息，更加全面有效地利用了个体的基因组数据信息和系谱信息，可更好地发挥基因组选择在农业育种领域的优势，预测候选群体个体的全基因组育种值并计算其与仿真数据中该个体的真实育种值之间的相关系数作为全基因组育种值预测的准确性指标，能准确预测以实现高效率特异性育种此外本发明提供的方法也有助于检测和确定与某些特异性状相关联的功能基因，有助于实现特异性育种。

Description

基于生物信息的选育良种系统及其算法

技术领域

本发明涉及生物技术领域，具体涉及了一种基于生物信息的选育良种系统及其算法。

背景技术

在农业生物(作物、畜禽和水产生物)的产业链中,良种选育占有至关重要的地位。传统的育种方法大多基于遗传统计原理并结合生产实际,利用候选个体本身和(或)亲属的性状记录估计得到其育种值，并利用育种值为指标通过进行多世代的人工选育以实现品种培育和性状改良。然而这种方法往往存在着育种周期长,人力物力消耗大等缺点。如法国的“大约克”种猪就花了三十多年时间才选育成功。在乳牛育种领域通过传统后裔测定方式验证一头候选种公牛需要5-6年时间，大约投入5万美元。

基于候选个体的遗传性状由其所携带的基因来决定的原理，通过利用候选个体的基因信息在基因层面来直接评估而非仅仅利用系谱信息和记录性状的表型数据来考量候选个体的遗传性状的优劣，将能够大大缩短育种时间及成本投入。利用候选个体的基因信息的育种方法主要包括两种。第一种方法称为标记辅助选择(Marker AssistedSelection,MAS)方法，其试图通过部分基因组信息，即确定控制表型性状的基因信息来实现对目标性状的选择。该方法被有效利用的前提条件是与控制表型性状的基因相关的数量性状位点(QTL)能够确定并被精确定位。然而在实际育种应用中，决定某一经济性状的往往有很多基因，而每个基因的效应都不明显，从而造成QTL检测和定位的困难。即便存在有对经济性状贡献较大的主效基因，在构成表型性状的所有遗传变异中应用于MAS的遗传标记只能捕获主效基因所带来的那部分变异，而小效应累加起来所带来的变异却被忽视了。所以MAS方法只能有限地应用在经济性状的控制基因明确且控制基因的效应显著的性状选择中。

为了捕获构成表型的所有遗传变异，需要在候选个体整个基因组水平上检测影响目标性状的所有QTL来进行选择利用，即在候选个体全基因组范围内进行标记辅助选择，这就是第二种方法：本发明所涉及的全基因组选择方法。该方法利用个体整个基因组高密度遗传标记图谱以保证影响目标性状的每一个QTL都与高密度全基因组标记图谱中的至少一个遗传标记处于连锁不平衡(Linkage Disequilibrium，LD)状态。通过这种连锁不平衡信息来实现对所有影响目标性状的QTL的捕捉，再通过LDMIP方法来实现良种选育。该方法比目前全基因组选择的另一常用方法GBLUP有更高的准确率。

全基因组遗传标记还可以提供系谱中个体间更加精确的亲缘遗传关系，即连锁分析(Linkage Analysis，LA)信息。利用个体及其父母的全基因组信息可以计算个体等位基因从父母处遗传的遗传几率，可以用来示踪等位基因在系谱中一代一代的传递。这种连锁分析信息对于例如检测和确定与疾病相关的功能基因位置方面非常重要:利用遗传标记可以示踪附近对疾病有效应的基因。这样如果标记附近有功能基因，并且样本足够大，携带不同标记等位基因的个体的表现就会显著不同。在品种选育中，人们除了追求高产优质等特性之外，能够有效抵御疾病的抗病能力也往往是重要的选育特性。通过连锁分析信息便可以检测在一个家系中等位基因与疾病的传递是否相关，有助于选育抗病能力强的个体。

发明内容

为解决上述技术问题，我们提出了一种基于生物信息的选育良种系统及其算法，其目的：准确预测以实现高效率特异性育种。

为达到上述目的，本发明的技术方案如下：

一种基于生物信息的选育良种系统，具体实施步骤如下：

一、采用全基因组选择方法,首先获取参考群体和候选群体个体的基因组数据，其中，采用SNP遗传标记来获取基因型数据；

二、获取系谱文件并进行预处理，向上追溯父母系谱信息至少5代及以上，系谱文件需要包含个体编号、父亲编号和母亲编号至少三段信息，在系谱文件中筛选具有基因型数据的基因分型个体，向上追溯父母系谱信息至少5代及以上以获得GA系谱；通过基因分型个体的基因型数据信息对系谱文件中父母及后代亲缘关系进行亲子鉴定以核查GA系谱文件，当根据基因型数据获得的亲子鉴定结果与原始GA系谱文件所显示的亲子关系不一致时，应以基因型数据获得的亲子鉴定结果为基准重新编排原始GA系谱文件；

三、在完成基因型数据的预处理并得到了GA系谱文件后，通过LDMIP方法来完成，具体的LDMIP方法，如下:

A、GA系谱中未基因分型的祖先个体，推知其基因型数据；

B、根据基因型数据和GA系谱信息，计算个体等位基因从父母处遗传的遗传几率；

C、通过LDMIP程序的计算结果可以建立起GA系谱中参考群体和候选群体之间的LDLA关系矩阵，这一关系矩阵既包括参考群体和候选群体之间基于基因型信息的个体亲缘关系，又包括参考群体和候选群体之间基于系谱信息的血缘关系，得到LDLA关系矩阵之后，对其求逆，然后通过混合模型方程组求解候选群体个体的全基因组育种值；

D、同时利用个体等位基因从父母处遗传的遗传几率计算可以示踪等位基因在系谱中世代传递信息，如可以用于检测和确定与疾病相关的功能基因，有助于选育抗病能力强的个体。

优选的，采用SNP遗传标记来获取基因型数据具体方式如下：通过SNP芯片获取参考群体和候选群体的基因型数据并进行预处理，首先基因分型个体的基因型数据如有缺失，对所缺失的基因型予以填充；然后通过SNP标记检出率、最小等位基因频率、Hardy-Weinberg平衡检验、个体检出率等参数优化数据质量。

一种基于生物信息的选育良种系统的算法，对于基因型数据的每个SNP位点的两个等位基因用0、1进行编码从而得到每个个体的每个SNP位点基因型编码，其中，如第i个个体的第j个SNP位点的编码gij,则对于纯合子“00”的基因型编码为“0”；对于纯合子“11”编码为“2”；对于杂合子“10”或“01”编码为“1”；然后通过LDMIP程序根据GA系谱首先预测祖先个体基因型编码信息，然后计算每个基因座位的等位基因从父母处遗传的遗传几率，利用等位基因遗传几率可以得到连锁分析信息，通过LDMIP程序还可以同时得到连锁不平衡信息，与连锁分析信息共同组建GLDLA关系矩阵，其中，GLDLA关系矩阵的行数和列数相同，均为GA系谱中个体的数目。

优选的，构建GLDLA关系矩阵要涉及到构建基于系谱信息的血缘关系矩阵－A矩阵，和基于SNP信息的亲缘关系矩阵－G矩阵，构建A矩阵的元素依据以下规则来递推确定，即对于系谱中一个体i和它的父亲s和母亲d，其与系谱中另一个体j的关系为:

一、如果它的父亲和母亲都已知，则

aij＝aji＝0.5(ajs+ajd)

aii＝1+0.5(asd)；

二、如果它的父亲s已知，母亲d未知，则

aij＝aji＝0.5(ajs)

aii＝1；

三、如果它的母亲d已知，父亲s未知，则

aij＝aji＝0.5(ajd)

aii＝1；

四、如果它的父亲和母亲都未知，则

aij＝aji＝0

aii＝1；

G矩阵的构建依据如下方法：

其中Nm为SNP位点数目，对于个体i和个体j，X中的元素Xij由以下方法确定：

其中gij为前述第i个个体的第j个SNP位点的基因型编码，pj为第j个SNP位点的第一个等位基因的基因频率，通过基因型数据样本计算得出；

在A矩阵和G矩阵都构建好之后，利用如下方法构建GLDLA关系矩阵：

GLDLA＝D*A*D+Δ*G*Δ；

而D和Δ分别是两个对角矩阵，其中D矩阵的对角元；

Δ矩阵的对角元

通过GLDLA关系矩阵，对其求逆，利用逆矩阵以及参考群体个体的表型数据可以建立如下混合模型方程组，

其中λ＝σ_e ²/σ_a ²，通过迭代求解混合模型方程组即可预测候选群体个体的全基因组育种值。

通过上述技术方案，预测候选群体个体的全基因组育种值并计算其与仿真数据中该个体的真实育种值之间的相关系数作为全基因组育种值预测的准确性指标，能准确预测以实现高效率特异性育种。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明所公开的一种基于生物信息的选育良种系统及其算法的流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

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

如图1所示，一种基于生物信息的选育良种系统，采用全基因组选择方法,首先获取参考群体和候选群体个体的基因组数据。对于基因型数据目前普遍采用SNP遗传标记来获取，对于牛、猪、鸡等禽畜物种已有商用的SNP芯片可以应用。通过SNP芯片获取参考群体和候选群体的基因型数据并进行预处理：首先基因分型个体的基因型数据如有缺失，对所缺失的基因型予以填充。然后通过SNP标记检出率、最小等位基因频率、Hardy-Weinberg平衡检验、个体检出率等参数优化数据质量。

获取系谱文件并进行预处理，向上追溯父母系谱信息至少5代及以上。系谱文件需要包含个体编号、父亲编号和母亲编号至少三段信息。在系谱文件中筛选具有基因型数据的基因分型个体，向上追溯父母系谱信息至少5代及以上以获得基因分型-祖先系谱(Genotyped-Ancestor pedigree,GA系谱)。通过基因分型个体的基因型数据信息对系谱文件中父母及后代亲缘关系进行亲子鉴定以核查GA系谱文件。当根据基因型数据获得的亲子鉴定结果与原始GA系谱文件所显示的亲子关系不一致时，应以基因型数据获得的亲子鉴定结果为基准重新编排原始GA系谱文件。

在完成基因型数据的预处理并得到了GA系谱文件后，通过LDMIP方法来完成:

一、GA系谱中未基因分型的祖先个体，推知其基因型数据

二、根据基因型数据和GA系谱信息，计算个体等位基因从父母处遗传的遗传几率。

通过LDMIP程序的计算结果可以建立起GA系谱中参考群体和候选群体之间的连锁不平衡-连锁分析关系矩阵(Linkage Disequilibrium Linkage Analysis,LDLA关系矩阵)。这一关系矩阵既包括参考群体和候选群体之间基于基因型信息的个体亲缘关系(连锁不平衡关系)，又包括参考群体和候选群体之间基于系谱信息的血缘关系(连锁分析关系)。得到LDLA关系矩阵之后，对其求逆，然后通过混合模型方程组求解候选群体个体的全基因组育种值。

同时利用个体等位基因从父母处遗传的遗传几率计算可以示踪等位基因在系谱中世代传递信息，如可以用于检测和确定与疾病相关的功能基因，有助于选育抗病能力强的个体。

利用仿真数据计算综合利用连锁不平衡和连锁分析关系信息。仿真数据利用QMSim模拟软件生成，通过模拟历史群2000代，每代200个个体，当代群8代，每代随机选取25头公畜与随机选取的250头母畜随机交配，每头母畜产生2个后代。模拟的基因组含有5条染色体，每条染色体的长度均为1摩尔根(Morgan),均匀分布10000个SNP标记，总共有50000个SNP标记。在当代群中，从第4代开始记录基因型信息。

通过模拟产生仿真系谱数据、仿真基因型数据文件，并且可以得到全部个体的真实育种值。然后产生两个性状的仿真表现型数据，两个性状的遗传力分别为0.5和0.1。对于基因分型的个体，选择最后一代500个个体作为候选群体，不产生表现型数据。而前四代共2000个个体为参考群体，通过模拟仿真表现型数据。

对此仿真数据，首先对数据文件进行预处理。对于系谱数据，根据2500个基因分型个体向上追溯10代系谱得到基因分型个体-祖先(GA)系谱包含5045个个体。对于基因型数据，每个SNP位点的两个等位基因用0、1进行编码从而得到每个个体的每个SNP位点基因型编码。如第i个个体的第j个SNP位点的编码gij,则对于纯合子“00”的基因型编码为“0”；对于纯合子“11”编码为“2”；对于杂合子“10”或“01”编码为“1”。

然后利用LDMIP程序根据GA系谱首先预测祖先个体基因型编码信息，然后计算每个基因座位的等位基因从父母处遗传的遗传几率。利用等位基因遗传几率可以得到连锁分析信息。通过LDMIP程序还可以同时得到连锁不平衡信息，与连锁分析信息共同组建连锁不平衡-连锁分析(GLDLA)关系矩阵。矩阵的行数和列数相同，均为GA系谱中个体的数目。矩阵的元素说明了个体间连锁不平衡-连锁分析关系。

构建GLDLA关系矩阵要涉及到构建基于系谱信息的血缘关系矩阵，A矩阵，和基于SNP信息的亲缘关系矩阵，G矩阵。构建A矩阵的元素依据以下规则来递推确定，即对于系谱中一个体i和它的父亲s和母亲d，其与系谱中另一个体j的关系为:

A、如果它的父亲和母亲都已知，则

aij＝aji＝0.5(ajs+ajd)

aii＝1+0.5(asd)；

B、如果它的父亲s已知，母亲d未知，则

aij＝aji＝0.5(ajs)

aii＝1；

C、如果它的母亲d已知，父亲s未知，则

aij＝aji＝0.5(ajd)

aii＝1；

D、如果它的父亲和母亲都未知，则

aij＝aji＝0

aii＝1；

G矩阵的构建依据如下方法：

其中Nm为SNP位点数目，对于个体i和个体j，X中的元素Xij由以下方法确定：

其中gij为前述第i个个体的第j个SNP位点的基因型编码，pj为第j个SNP位点的第一个等位基因的基因频率，通过基因型数据样本计算得出。

在A矩阵和G矩阵都构建好之后，利用如下方法构建GLDLA关系矩阵：

GLDLA＝D*A*D+Δ*G*Δ

而D和Δ分别是两个对角矩阵。其中D矩阵的对角元

Δ矩阵的对角元

利用连锁不平衡-连锁分析(GLDLA)关系矩阵，对其求逆，利用逆矩阵以及参考群体个体的表型数据可以建立如下混合模型方程组，

其中λ＝σ_e ²/σ_a ²，通过迭代求解混合模型方程组即可预测候选群体个体的全基因组育种值。

为检验本发明方法的实施效果，预测候选群体个体的全基因组育种值并计算其与仿真数据中该个体的真实育种值之间的相关系数作为全基因组育种值预测的准确性指标。准确性越高说明本发明的实施效果越好。同时作为比较，分别只利用血缘关系矩阵A矩阵并求逆，然后以A-1代替混合模型方程组中的GLDLA-1求解混合模型方程组，所得的育种值则为利用系谱数据ABLUP的预测育种值。如果只利用亲缘关系矩阵G矩阵并求逆，以G-1代替混合模型方程组中的GLDLA-1求解混合模型方程组所得的育种值则为普通全基因组GBLUP的预测育种值。将ABLUP和GBLUP的预测育种值分别与个体的真实育种值相关，得出ABLUP和GBLUP的育种值预测准确性。下表为三种方法的育种值预测准确性对比结果：

性状遗传力h2	GLDLA方法	GBLUP方法	ABLUP方法
				0.5	0.805	0.740	0.620
0.1	0.563	0.531	0.408

可以看到采用同时应用连锁不平衡-连锁分析关系的GLDLA方法的准确性要高于普通基因组选择应用的GBLUP方法，更大大高于传统的ABLUP方法。

上述中的技术术语如下：

表型性状：决定一个生物个体可能不同于其他生物个体的形态、生理、生化和行为等数量性状。

遗传标记：是已知在染色体上位置的一种基因或DNA序列，可被用于鉴定生物个体或物种。其可被描述为可观测变异。

数量形状位点(QTL)：指占据一特定染色体区域的微效多基因群。

单核苷酸的多态性(SNP):是指在基因组上单个核苷酸的变异，包括转换、颠换、缺失和插入，形成的遗传标记。

等位基因:又称对偶基因，是一些占据染色体的基因座的可以复制的脱氧核糖核酸。

育种值：是指种畜的种用价值。在数量遗传学中把决定某一数量性状的基因加性效应值定义为该量性状的育种值。

全基因组育种值：是指在全基因组范围内通过基因组中标记效应累加得到的育种值。

以上所述的仅是本发明的优选实施方式，应当指出，对于本领域的普通技术人员来说，在不脱离本发明创造构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。

Claims (4)

1.一种基于生物信息的选育良种系统，其特征在于，具体实施步骤如下：

一、采用全基因组选择方法,首先获取参考群体和候选群体个体的基因组数据，其中，采用SNP遗传标记来获取基因型数据；

二、获取系谱文件并进行预处理，向上追溯父母系谱信息至少5代及以上，系谱文件需要包含个体编号、父亲编号和母亲编号至少三段信息，在系谱文件中筛选具有基因型数据的基因分型个体，向上追溯父母系谱信息至少5代及以上以获得GA系谱；通过基因分型个体的基因型数据信息对系谱文件中父母及后代亲缘关系进行亲子鉴定以核查GA系谱文件，当根据基因型数据获得的亲子鉴定结果与原始GA系谱文件所显示的亲子关系不一致时，应以基因型数据获得的亲子鉴定结果为基准重新编排原始GA系谱文件；

三、在完成基因型数据的预处理并得到了GA系谱文件后，通过LDMIP方法来完成，具体的LDMIP方法，如下:

A、GA系谱中未基因分型的祖先个体，推知其基因型数据；

B、根据基因型数据和GA系谱信息，计算个体等位基因从父母处遗传的遗传几率；

C、通过LDMIP程序的计算结果可以建立起GA系谱中参考群体和候选群体之间的LDLA关系矩阵，这一关系矩阵既包括参考群体和候选群体之间基于基因型信息的个体亲缘关系，又包括参考群体和候选群体之间基于系谱信息的血缘关系，得到LDLA关系矩阵之后，对其求逆，然后通过混合模型方程组求解候选群体个体的全基因组育种值；

D、同时利用个体等位基因从父母处遗传的遗传几率计算可以示踪等位基因在系谱中世代传递信息，如可以用于检测和确定与疾病相关的功能基因，有助于选育抗病能力强的个体。

2.根据权利要求1所述的一种基于生物信息的选育良种系统，其特征在于，采用SNP遗传标记来获取基因型数据具体方式如下：通过SNP芯片获取参考群体和候选群体的基因型数据并进行预处理，首先基因分型个体的基因型数据如有缺失，对所缺失的基因型予以填充；然后通过SNP标记检出率、最小等位基因频率、Hardy-Weinberg平衡检验、个体检出率等参数优化数据质量。

3.权利要求1所述的一种基于生物信息的选育良种系统的算法，其特征在于，对于基因型数据的每个SNP位点的两个等位基因用0、1进行编码从而得到每个个体的每个SNP位点基因型编码，其中，如第i个个体的第j个SNP位点的编码gij,则对于纯合子“00”的基因型编码为“0”；对于纯合子“11”编码为“2”；对于杂合子“10”或“01”编码为“1”；然后通过LDMIP程序根据GA系谱首先预测祖先个体基因型编码信息，然后计算每个基因座位的等位基因从父母处遗传的遗传几率，利用等位基因遗传几率可以得到连锁分析信息，通过LDMIP程序还可以同时得到连锁不平衡信息，与连锁分析信息共同组建GLDLA关系矩阵，其中，GLDLA关系矩阵的行数和列数相同，均为GA系谱中个体的数目。

4.根据权利要求3所述的一种基于生物信息的选育良种系统的算法，其特征在于，构建GLDLA关系矩阵要涉及到构建基于系谱信息的血缘关系矩阵－A矩阵，和基于SNP信息的亲缘关系矩阵－G矩阵，构建A矩阵的元素依据以下规则来递推确定，即对于系谱中一个体i和它的父亲s和母亲d，其与系谱中另一个体j的关系为:

一、如果它的父亲和母亲都已知，则

aij＝aji＝0.5(ajs+ajd)

aii＝1+0.5(asd)；

二、如果它的父亲s已知，母亲d未知，则

aij＝aji＝0.5(ajs)

aii＝1；

三、如果它的母亲d已知，父亲s未知，则

aij＝aji＝0.5(ajd)

aii＝1；

四、如果它的父亲和母亲都未知，则

aij＝aji＝0

aii＝1；

G矩阵的构建依据如下方法：

其中Nm为SNP位点数目，对于个体i和个体j，X中的元素Xij由以下方法确定：

其中gij为前述第i个个体的第j个SNP位点的基因型编码，pj为第j个SNP位点的第一个等位基因的基因频率，通过基因型数据样本计算得出；

在A矩阵和G矩阵都构建好之后，利用如下方法构建GLDLA关系矩阵：

GLDLA＝D*A*D+Δ*G*Δ；

而D和Δ分别是两个对角矩阵，其中D矩阵的对角元；

Δ矩阵的对角元

通过GLDLA关系矩阵，对其求逆，利用逆矩阵以及参考群体个体的表型数据可以建立如下混合模型方程组，

其中λ＝σ_e ²/σ_a ²，通过迭代求解混合模型方程组即可预测候选群体个体的全基因组育种值。