CN103914632A - 一种快速估计基因组育种值的方法和应用 - Google Patents

Abstract

本发明属于生物信息学领域，提供了一种快速估计基因组育种值的方法，包括SNP芯片数据编辑；整理系谱文件；利用系谱信息构建A阵，从A阵中挑出具有SNP基因型的个体，根据个体间A阵元素构建子矩阵A₂₂；利用高密度SNP芯片信息构建G阵，利用公共开源LAPACK数学库中DGEMM过程快速构建G阵；构建H矩阵，调用公共开源LAPACK数学库中DGETRF和DGETRI矩阵求逆子程序，快速求得A阵、A₂₂阵以及G阵的逆矩阵，进而根据以上分块逆矩阵求解H逆矩阵；通过混合模型方程组求解个体的基因组育种值。本方法加快了利用综合系谱信息和基因组信息估计基因组育种值的计算速度，可推进基因组选择在动物育种领域的应用。

Description

一种快速估计基因组育种值的方法和应用

技术领域

本发明属于生物信息学技术领域，具体涉及一种利用系谱信息和基因组信息快速估计基因组育种值的方法和应用。 

背景技术

在畜禽育种领域，人工选择在品种培育和性状改良方面起着举足轻重的作用。自2007年来，随着各畜禽的全基因组高密度SNP芯片相继问世，基因组选择开始广泛用于畜禽育种领域，尤其奶牛育种。 

基因组选择的广泛应用与相应计算方法的发展密不可分。基于SNP效应方差的不同，基因组选择计算方法主要分为两大类：BLUP类方法和BAYES类方法。无论是BAYES类方法还是BLUP类方法，都是基于具有表型信息的SNP分型参考群个体，预测SNP分型候选个体，仅利用了参考群体和候选群体的基因组信息，而参考群体大小对基因组育种值估计准确性影响很大，参考群体越大，则基因组育种值估计准确性越高。但目前高密度SNP芯片成本仍然较高，这制约了基因组选择在更多畜禽分子育种中的应用。另一方面，基因组选择需要依赖传统的常规育种体系，需要根据常规育种估计的传统育种值作为新的表型，用来预测候选群体的基因组育种值，这造成了信息的双重利用，在BLUP方法中尤为明显。在传统育种值估计中，已经利用了基于系谱信息的个体亲缘关系，而在基因组育种值估计中，又利用了基于基因组信息的个体间亲缘关系，因此造成了个体关系的双重计算(double counting)，从而引起基因组育种值估计偏差。因此如果能够同时利用系谱信息和基因组信息，直接根据个体的表型测定值估计候选群体的基因组育种值，则不仅能提高基因组育种值估计的准确性，而且能够将常规育种体系和基因组选择体系融为一体，使畜禽育 种体系更完善。 

目前，综合基因组信息和系谱信息的估计基因组育种值方法仍存在实施上的困难，主要体现在(1)如何挑选无SNP基因型个体加入参考群体，不恰当的加入系谱信息不仅不能提高基因组育种值估计的准确性，而且带来计算时间的呈指数增加，满足不了实际育种的需要；(2)个体亲缘关系矩阵的快速构建和求逆。综合基因组信息和系谱信息能够更准确反映个体间亲缘关系，如何快速构建包含系谱信息和基因组信息，反映参考群体和候选群体个体间关系的H矩阵并求逆，不恰当的方法可能导致无法求解，将直接影响到基因组育种值评估的效率。因此，急需解决快速准确估计基因组育种值的问题。 

发明内容

育种值：种畜的种用价值，在数量遗传学中把决定某一数量性状的基因加性效应总和称为某一性状的个体育种值。 

基因组育种值：个体全基因组的SNP效应累加得到的育种值。 

参考群体：群体内个体具有SNP芯片基因型信息和表型数据信息，根据此参考群体可以估计整个基因组SNP标记效应，进而预测候选群体个体的基因组育种值。 

候选群体：由仅具有SNP芯片基因型信息的个体组成。 

针对现有技术不足，本发明的目的是提供一种快速估计基因组育种值的方法。 

为实现上述目的，本发明提供了一种快速估计基因组育种值的方法，包括以下步骤， 

步骤S1，SNP芯片数据编辑； 

步骤S2，整理系谱文件； 

步骤S3，利用系谱信息构建亲缘关系矩阵-A阵，从A阵中挑出具有SNP基因型的个体，根据个体间A阵元素构建子矩阵A₂₂； 

步骤S4，利用高密度SNP芯片信息构建基因组关系矩阵-G阵,利 用公共开源LAPACK数学库中DGEMM过程快速构建G阵，所述G阵根据构建，每个SNP位点两个等位基因用1、2进行编码，p_i代表第i个SNP第2个等位基因的频率，z的行数代表个体数，z的列数代表所用的SNP位点数； 

步骤S5，构建H矩阵，包含基于系谱信息的矩阵A和基于基因组信息的矩阵G，调用公共开源LAPACK数学库中DGETRF和DGETRI矩阵求逆子程序，快速求得A阵、A₂₂阵以及G阵的逆矩阵，进而根据以上分块逆矩阵，获得综合系谱信息和基因组信息的H逆矩阵； 

步骤S6,根据H矩阵，获取H矩阵中个体表型值； 

步骤S7，通过混合模型方程组求解个体的基因组育种值。 

优选的，所述步骤S1具体包括以下步骤， 

步骤S11，SNP芯片数据读取，读取SNP芯片数据，并用压缩格式存储文件以节省硬盘空间； 

步骤S12，缺失基因型填充，利用Beagle程序将芯片基因型中存在缺失的SNP标记或个体进行填充，提高芯片基因型检测质量； 

步骤S13，质量控制，所述质量控制参数为每个SNP标记检出率(call rate)和最小等位基因频率(MAF)、哈代-温伯格平衡检验和/或个体检出率(call rate)。 

优选的，步骤S2中所述整理系谱文件是根据SNP芯片数据个体，从整体系谱文件中挑出具有SNP芯片信息的个体，向上追溯父母系谱信息5-10代，并根据SNP信息对系谱文件父母及后代亲缘关系进行亲子鉴定并根据亲子鉴定结果调整原始系谱文件，主要表现为原始系谱文件亲子关系与亲子鉴定结果不一致时，系谱文件按亲子鉴定结果重新编排。系谱文件包含个体、母号、父号三个字段信息。 

优选的，所述步骤S4还包括根据A矩阵校正G矩阵，生成新的矩阵G*，以使基于系谱信息的A矩阵和基于基因组信息的G阵尺度相同或接近，降低估计偏差。 

本发明的另一目的是提供上述方法在动物育种方面的应用。 

本发明的有益效果：本发明提供一种快速估计基因组育种值的方法，在畜禽育种领域首次利用公共开源函数库LAPACK解决矩阵运算问题，从而加快了利用综合系谱信息和基因组信息估计基因组育种值的计算速度。可推进基因组选择在国内动物育种领域的应用，可更好地发挥基因组选择在动物育种领域的优势。 

附图说明

图1实施例1中本发明方法的流程图。 

具体实施方式

以下实施例用于说明本发明，但不用来限制本发明的范围。 

实施例1 

参见图1，本实施例提供了一种快速估计基因组育种值的方法，包括以下步骤， 

1)获取数据，数据来源于5974头中国荷斯坦牛母牛，出生于2004-2012年间，所有母牛进行了Illunima50K SNP芯片基因型测定，并对5个产奶性状产奶量、乳蛋白量、乳蛋白率、乳脂量、乳脂率进行了传统育种值估计。并对SNP芯片数据编辑，共47160SNP用于分析；根据5974头具有SNP基因型母牛系谱信息，向上追溯10代资料，包含130852个体，用于亲缘关系矩阵A阵构建；选取传统育种值可靠性高于0.45的个体5个产奶性状育种值作为表型，生成最终表型数据文件。 

2)利用系谱信息构建亲缘关系矩阵A阵。根据步骤1)中系谱追溯的130852个体构建A阵并求逆，选取其中的22106未进行基因型测定且是候选群体半同胞或母亲的母牛加入参考群体，从A逆中读取这22106个体和5974头具有SNP基因型母牛相应的值。 

A矩阵是所有动物个体之间的加性遗传相关矩阵。个体i和j之间的加性遗传相关是指在它们的基因组中具有同源相同基因的比列，或 者从个体i的基因组中随机抽取的一个基因与从个体j的基因组中随机抽取的一个基因同源相同的概率。 

构建A阵中的每一个元素采用以下的递推公式来计算： 

其中，s_i(s_j)和d_i(d_j)为个体i(j)的父亲和母亲。 

系谱文件的特点：所有个体按个体号，父亲号和母亲号列成一个三列表，且在个体一列中应包括所有在父和母列出现过的个体，在个体一列中应保证后代绝不会出现在其父母之前，一般可按出生日期排序，先出生的在前，为便于编写程序，个体应用自然数1开始连续编号。 

3)根据基因组信息构建个体关系矩阵G阵。根据5974头有基因型信息的母牛构建G阵，根据构建G阵，每个SNP位点两个等位基因用1、2进行编码，p_i代表第i个SNP等位基因2的频率(通过样本计算得出)，z的行数代表个体数，z的列数代表所用的SNP位点数，对于每个元素，若纯合子11,为0-2p_i；纯合子22，为2-2p_i；杂合子12或21，为1-2p_i。 

本步骤利用公共开源LAPACK数学库中过程DGEMM快速求解zz’，大大节省G矩阵构建时间。 

4)G矩阵校正。根据G矩阵与A矩阵尺度比例判断是否需要校正G阵，包括1)采用上述的构建方法，不再进行校正；2)采用上述的构建方法，然后利用公式G^*＝α+βG进行校正，生成新矩阵G*，其 中系数α、β来自方程组, 

Avg(diag(G))β+α＝Avg(diag(A)) 

Avg(offdiag(G))β+α＝Avg(offdiag(A))。 

5)构建H矩阵并求解H逆矩阵。 

对于实际应用中某些个体没有基因型记录的问题，利用新的关系矩阵构建，将系谱个体和基因型个体整合在一起，扩展到没有基因型而有表型记录的个体，新的矩阵被称为H矩阵，即： 

$H=\left[\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right]\left[\begin{array}{cc}{A}_{11}+A_{12}{A}_{22}^{-1}(G-A_{22})A_{22}^{-1}{A}_{21}& A_{12}{A}_{22}^{-1}G\\ G{A}_{22}^{-1}{A}_{21}& G\end{array}\right],$ 包含了A矩阵和G矩阵。 

基因组育种值估计中需要利用H逆矩阵信息，直接根据H矩阵求解会造成计算比较繁琐甚至当数据量大时无法求解。因此，跳过求H阵的繁琐过程，直接根据A逆、G逆或G*逆及A₂₂逆求解H逆矩阵，即 

$H^{-1}=A^{-1}+\left[\begin{array}{cc}0& 0\\ 0& G^{-1}-A_{22}^{-1}\end{array}\right]$

其中，A₂₂代表具有SNP基因型个体间的亲缘关系矩阵，是A矩阵的子矩阵，G矩阵代表基因组关系矩阵。 

即使如此简化，H^-1的快速求解仍然十分困难，制约了该方法在基因组选择中的高效应用，因此需要解决H^-1快速求解的问题。本方法调用公共开源LAPACK数学库中DGETRF和DGETRI矩阵求逆子程序，快速求得A阵、A₂₂阵以及G阵、G^*的逆矩阵。进而根据以上分块逆矩阵，获得综合系谱信息和基因组信息的H逆矩阵。 

6)建立混合模型方程组，利用共轭梯度迭代方法估计产奶量基因组育种值。混合模型方程组可表示为： 

(1)若e～N(0,Iσ_e ²)时， 

$\left[\begin{array}{cc}{1}^{\text{'}}1& 1^{\text{'}}Z\\ Z^{\text{'}}1& Z^{\text{'}}Z+A^{-1}k\end{array}\right]\left[\begin{array}{c}\mathrm{\μ}\\ \mathrm{\α}\end{array}\right]=\left[\begin{array}{c}{1}^{\text{'}}y\\ Z^{\text{'}}y\end{array}\right]$

(2)若e～N(0,Rσ_e ²)时， 

$\left[\begin{array}{cc}{1}^{\text{'}}{R}^{-1}1& 1^{\text{'}}{R}^{-1}Z\\ Z^{\text{'}}{R}^{-1}1& Z^{\text{'}}{R}^{-1}Z+A^{-1}k\end{array}\right]\left[\begin{array}{c}\mathrm{\μ}\\ \mathrm{\α}\end{array}\right]=\left[\begin{array}{c}{1}^{\text{'}}{R}^{-1}y\\ Z^{\text{'}}{R}^{-1}y\end{array}\right]$

其中，可通过迭代求解个体基因组育种值，(1)和(2)中A^-1用H^-1代替。 

7)重复步骤6)，估计乳蛋白量等其他4个产奶性状基因组育种值。实施例2应用真实数据评价本方法计算速度和准确性 

使用实施例1中数据，分别利用本发明计算构建H阵和求逆方法，与传统的对称矩阵求逆、R语言ginv函数比较计算时间。同时，与传统仅利用基因组信息的GBLUP方法比较预测基因组育种值准确性。将实施例1中5974头具有SNP基因型的541头2008年和2008后出生的母牛挑出，其传统育种值与基因组育种值相关作为标准衡量基因组育种值预测的准确性，相关越高说明基因组育种值估计越准确，说明本发明的计算方法效果越好。 

表1A、G、H阵构建和求逆及基因组育种值(GEBV)估计计算时间 

方法	H逆矩阵(分)1	A、G构建(分)2	GEBV估计时间(分)3
				实施例1	9	1	20
对称矩阵方法	75	120	360
				R语言函数ginv	5天	24小时	--

1：因为对称矩阵求逆方法和R语言函数ginv处理大数据时耗时过长，无法完成运算，仅比较了对三种方法对实施例1数据中5974个体构成的H逆矩阵比较； 

2：同1，为5974个体； 

3：5974头具有SNP基因型个体和22106头无基因型个体。 

由表1可以看出，在同样硬件配置下(服务器的基本配置为64核、128个CPU，8T硬盘以及512G内存)，本发明能够很短时间内完成。而传统的对称矩阵方法则无法在短时间内完成，A、G矩阵构建耗时120分钟，求H逆矩阵75分钟，计算时间是本发明10分钟的20倍，而 整个基因组育种值估计传统对称矩阵方法需要6个小时，本方法仅20分钟。同样，目前流行的R语言，仅求H逆矩阵就耗时5天，估计所有个体的基因组育种值则计算2个星期仍未完成，以此无法用于实际评估需要，仅能限于教学演示。 

表2基因组育种值估计准确性 

由表2可以看出，本发明由于综合利用了系谱信息和基因组信息，与传统的仅利用基因组信息的GBLUP方法相比，准确性提高幅度为12-16％，平均13％，能更准确估计个体的基因组育种值。同时，校正方法提高幅度不大，说明根据系谱信息构建的A矩阵和基于基因组信息的G矩阵尺度差别不大，基本不需要校正。 

实施例3仿真数据 

仿真数据是由2011年QTL-MAS Workshop利用LDSO软件模拟的一个远交群体，通过模拟历史群1000代，每代1000个体，当代群30代，每代150个体而构建的。分析所使用的数据是最后一代数据，共20个公畜家系，每个公畜与100个母畜交配，每个母畜产生15个后代。9990个SNP标记均匀分布在5条染色体上，每条染色体长度为1Morgen，均匀分布着1998SNP标记。 

仿真数据提供的具体数据信息为：系谱文件，SNP标记文件，以及表型文件，其中表型文件是每个15全同胞后代中，随机模拟其中10个体的表型值，而其他的5个个体为验证个体。同时提供3000后代的 真实育种值。 

对此仿真数据，计算方法同实施例1，选取本发明、GBLUP方法和贝叶斯类方法(BayesA,BayesB,BayesCpi)预测1000个后代的基因组育种值，并用基因组育种值与其真实育种值的相关系数作为基因组育种值预测的准确性，相关系数越高，说明基因组育种值估计越准确。同时将基因组育种值对真实育种值回归，回归系数衡量基因组育种值估计的无偏性，回归系数越接近1.0，说明基因组育种值估计越无偏。 

由表3可以看出，贝叶斯类方法由于假设SNP效应方差不同，更符合“微效多基因假说”遗传理论，因此准确性比同样只利用基因组信息的GBLUP方法高。但相比本发明，由于本发明整合了系谱信息和基因组信息，因此准确性最高，同时无偏性也最小。在计算时间方面，由于贝叶斯类方法需要借助吉布斯抽样技术估计SNP效应，因此耗时很长，三种方法都在4小时以上，并且计算时间会随着SNP标记数、样本大小呈指数性增长。而本发明和GBLUP在20分钟之内即完成运算，展示了实际应用潜力。 

表31000个后代基因组育种值估计准确性和无偏性 

虽然，上文中已经用一般性说明、具体实施方式及试验，对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。 

Claims (5)

1.一种快速估计基因组育种值的方法，其特征在于，包括以下步骤，

步骤S1，SNP芯片数据编辑；

步骤S2，整理系谱文件；

步骤S3，利用系谱信息构建亲缘关系矩阵-A阵，从A阵中挑出具有SNP基因型的个体，根据个体间A阵元素构建子矩阵A₂₂；

步骤S4，利用高密度SNP芯片信息构建基因组关系矩阵-G阵,利用公共开源LAPACK数学库中DGEMM过程快速构建G阵，所述G阵根据构建，每个SNP位点用1、2进行等位基因编码，p_t代表第i个SNP第2个等位基因频率，z的行数代表个体数，z的列数代表所用的SNP位点数；

步骤S5，构建H矩阵，包含基于系谱信息的矩阵A和基于基因组信息的矩阵G，调用公共开源LAPACK数学库中DGETRF和DGETRI矩阵求逆子程序，快速求得A阵、A₂₂阵以及G阵的逆矩阵，进而根据以上分块逆矩阵，获得综合系谱信息和基因组信息的H逆矩阵；

步骤S6,根据H矩阵，获取H矩阵中个体表型值；

步骤S7，通过混合模型方程组求解个体的基因组育种值。

2.根据权利要求1所述的方法，其特征在于，所述步骤S1具体包括以下步骤，

步骤S11，SNP芯片数据读取，读取SNP芯片数据，并用压缩格式存储文件以节省硬盘空间；

步骤S12，缺失基因型填充，利用Beagle程序将芯片基因型中存在缺失的SNP标记或个体进行填充，提高芯片基因型检测质量；

步骤S13，质量控制，所述质量控制参数为每个SNP标记检出率和最小等位基因频率、哈代-温伯格平衡检验和/或个体检出率。

3.根据权利要求1所述的方法，其特征在于，步骤S2中所述整理系谱文件是根据SNP芯片数据个体，从整体系谱文件中挑出具有SNP芯片信息的个体，向上追溯父母系谱信息5-10代，并根据SNP信息对系谱文件父母及后代亲缘关系进行亲子鉴定并根据亲子鉴定结果调整原始系谱文件，主要表现为原始系谱文件亲子关系与亲子鉴定结果不一致时，系谱文件按亲子鉴定结果重新编排。

4.根据权利要求1所述的方法，其特征在于，所述步骤S4还包括根据A矩阵校正G矩阵，生成新的矩阵G*，以使基于系谱信息的A矩阵和基于基因组信息的G阵尺度相同或接近，降低估计偏差。

5.根据权利要求1-4任一所述的方法在动物育种方面的应用。