一种用于无家系低遗传力品种的基因组选择方法及应用与流程

1.本发明涉及遗传育种技术领域，具体而言，涉及一种用于无家系低遗传力品种的基因组选择方法及应用。

背景技术：

2.粮食安全的重要基础是优质种质资源的利用。然而种业整体呈橄榄型分布，腰部的包括玉米、大豆、生猪、肉牛、肉羊以及诸多蔬菜品种的种子质量不高，特别是畜牧业，优质种源大部分来自国外。优质种质的遗传基础缺少深度研究，严重影响了种质资源利用。以水稻为例，除了少数不育基因、抗病和抗虫基因外，众多基因亟待发掘。所以难以从基因编辑水平提出优质的选育方案，常规育种大多靠经验育种，导致选择概率小，难育出突破性品种。
3.目前的基因组选择计算模型主要分以blup代表的直接法以及bayes代表的间接法。这两种方法均假设snps之间相互独立，且snps效应是同分布。但是实际上有效应的snps往往富集在qtls上，存在很明显的相关性。一方面，snps之间的相关性会造成估计的误差，另一方面，与性状不相关的snps也会稀释有效应的snps的作用。目前模型对snps数量以及参考体数量要求很高。若希望使用低密度snps进行基因组选择育种，snps之间非等效性与snps之间的相关性是无法规避的问题。因此，急需一种用于无家系低遗传力品种的基因组选择方法，在假设snps之间存在某种相关性且snps非等效的条件下建立了一种对snps效应值估计的模型。

技术实现要素：

4.为了解决上述问题，本发明的目的是提供一种用于无家系低遗传力品种的基因组选择方，通过建立一种全基因组选择的计算模型，为基因组选择育种提供一种弱假设模型cr-elastic net，解决全基因组选择育种现有模型对snps数量以及无家系参考体数量要求高的问题，进而实现通过更符合snps变异规律的假设实现对动物育种值估计准确度的提升；以及为无家系低遗传力体提供一种优于先进水平的模型，达到使用少量个体即可有效估计无家系体育种值的目的，推动种业快速发展。
5.为了实现上述技术目的，本技术提供了一种用于无家系低遗传力品种的选育方法，包括以下步骤：
6.基于snps的非等效条件，通过校正snps之间的相关性关系，评估snps效应值，对无家系低遗传力品种进行选育。
7.优选地，在评估snps效应值的过程中，根据非等效条件和相关性关系，构建用于评估snps效应值的cr-elastic net模型，通过设置常量微调惩罚量和模型代价函数，对snps效应值进行评估。
8.优选地，在构建cr-elastic net模型的过程中，elastic net模型的表达式为：
[0009][0010]
式中，w为snps效应值向量；w
t
为向量w的转置；||w||为w的转置共轭矩阵与矩阵w的积的最大特征根的平方根值，||w||＝w
t
w；为w的欧几里得范数的平方,为w里每个元素的平方和；yi为第i个观测的表型值；xi为第i个观测的基因型为一个n
×
1的向量，n为snps数量，λ为常量微调惩罚量，n为样本量，ρ为一个0-1之间的常数，其中，当ρ为0时模型代价函数等价于岭回归代价函数，ρ为1时模型代价函数等价于lasso回归代价函数。
[0011]
优选地，根据经验设定一个snps之间相关性的参数p，当snps之间相关系数低于1-p则认为其相互独立，高于1-p则认为其相关。通过模型拟合设置在该相关性参数下模型解释残差比例最高的λ.在以上p与λ下，获取snps的效应值。
[0012]
优选地，在评估snps效应值的过程中，通过重抽样，从总体中有放回抽取样本容量为n的m个子集，对m个子集通过弹性网络拟合，得到m组snps效应值wk，snps的分布服从原分布，snps真实效应值依概率收敛到{wk}的平均值。
[0013]
本发明公开了一种用于无家系低遗传力品种的选育应用，用于实现无家系低遗传力品种的选育系统，包括：
[0014]
第一分析模块，用于假设所述snps之间的相关性关系阈值；
[0015]
第二分析模块，用于获取在第一分析模块下最优的常量微调惩罚量；
[0016]
评估模块，用于评估所述snps效应值；
[0017]
选育模块，用于根据所述snps效应值，对无家系低遗传力品种进行选育。
[0018]
本发明公开了以下技术效果：
[0019]
本发明提供的cr-elastic net模型可用于无家系低遗传力品种的选育，有利推动种业发展。
附图说明
[0020]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0021]
图1是本发明所述的方法流程图。
[0022]
图2是在不同参考体数目下，cr-elastic net与gblup相对于传统pblup在50个模拟性状预测能力存在稳定提升的性状数。
[0023]
图3是在50个模拟性状中使用不同参考体数目cr-elastic net相对gblup预测能力提升超过5％，10％，20％的性状数目。
[0024]
图4是在真实数据集(半滑舌鳎)中cr-elastic net与gblup的平均交叉验证预测相关性对比。
[0025]
图5是在真实数据集(猪)中cr-elastic net与gblup在不同参考体数目下的平均交叉验证预测相关性对比。
具体实施方式
[0026]
为使本技术实施例的目的、技术方案和优点更加清楚，下面将结合本技术实施例中附图，对本技术实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本技术一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本技术实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本技术的实施例的详细描述并非旨在限制要求保护的本技术的范围，而是仅仅表示本技术的选定实施例。基于本技术的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本技术保护的范围。
[0027]
如图1所示，本发明提供了一种用于无家系低遗传力品种的选育方法，包括以下步骤：
[0028]
基于snps的非等效条件，通过校正snps之间的相关性关系，评估snps效应值，对无家系低遗传力品种进行选育。
[0029]
进一步优选地，本发明在评估snps效应值的过程中，本发明根据非等效条件和相关性关系，构建用于评估snps效应值的cr-elastic net模型，通过设置常量微调惩罚量和模型代价函数，对snps效应值进行评估。
[0030]
进一步优选地，本发明在构建cr-elastic net模型的过程中，本发明提到的elastic net模型的表达式为：
[0031][0032]
式中，w为snps效应值向量；w
t
为向量w的转置；||w|为w的转置共轭矩阵与矩阵w的积的最大特征根的平方根值，||w||＝w
t
w；为w的欧几里得范数的平方,即为w里每个元素的平方和；yi为第i个观测的表型值；xi为第i个观测的基因型为一个n
×
1的向量，n为snps数量，λ为常量微调惩罚量，n为样本量，ρ为一个0-1之间的常数，其中，当ρ为0时模型代价函数等价于岭回归代价函数，ρ为1时模型代价函数等价于lasso回归代价函数。进一步优选地，本发明在设置常量微调惩罚量的过程中，本发明提到的常量微调惩罚量越大l2-正则化惩罚项的影响越大，当常量微调惩罚量为0时无惩罚。
[0033]
进一步优选地，本发明在模型代价函数等价于岭回归代价函数时，保留所有snps。
[0034]
进一步优选地，本发明在当模型代价函数等价于lasso回归代价函数时，如果snps之间存在相关性，则随机选择其中一个snps。
[0035]
进一步优选地，本发明在评估snps效应值的过程中，本发明通过重抽样，从总体中有放回抽取样本容量为n的m个子集，对m个子集通过弹性网络拟合，得到m组snps效应值wk，snps的分布服从原分布，snps真实效应值依概率收敛到{wk}的平均值。
[0036]
本发明公开了一种用于无家系低遗传力品种的选育应用，用于实现无家系低遗传力品种的选育系统，包括：
[0037]
第一分析模块，用于假设所述snps之间的相关性关系阈值；
[0038]
第二分析模块，用于获取在第一分析模块下最优的常量微调惩罚量；
[0039]
评估模块，用于评估所述snps效应值；
[0040]
选育模块，用于根据所述snps效应值，对无家系低遗传力品种进行选育。
[0041]
本发明还公开了一种计算机程序，通过计算机程序实现了用于无家系低遗传力品
种的选育方法，进而将该程序嵌入到智能设备中，用于在假设snps的非等效条件以及snps之间的相关性关系下，生成snps效应值，并根据snps效应值，实现对无家系低遗传力品种进行智能选育。
[0042]
为了更好的理解本发明的含义，下面对本发明所涉及的术语解释如下：
[0043]
snp：单核苷酸多态性，single nucleotide polymorphism的缩写，表示在基因组水平上由单个核苷酸变异引起的dna序列多态性。
[0044]
参考体：在基因组选择中，具有表型和基因型的体，表型通过人工测量获得，基因型通过基因组测序获得，该体用于训练模型。
[0045]
候选体：在基因组选择中，只有基因型而没有表型的体，通常为选育过程中的备选个体，测定基因型后便可结合参考体估算这些个体的育种值，通过育种值的高低进行优选优育。
[0046]
cr-elastic net：可重复抽样弹性网络，是一种基于重抽样理论与弹性网络以及最小均方差估计理论，直接估计snps效应，结合分子标记以达到估算个体育种值的方法。
[0047]
gblup：基因组最佳线性无偏预测，genomic best linear unbiased prediction的缩写，是一种利用覆盖基因组的高密度分子标记构建基因组亲缘关系矩阵g矩阵估算个体育种值的方法。
[0048]
pblup:系谱最佳线性无偏预测，pedigree-based best linear unbiased prediction的缩写，是一种基于系谱估算个体育种值的方法。
[0049]
gebv：基因组估计育种值，genomic estimated breeding value的缩写，通过基因组选择方法在基因组水平上估算出的育种值。
[0050]
tbv：模拟数据真实育种值，定义为个体基因组上所有qtls效应的加和。
[0051]
mse准则：minimization criterion of mean-square error是一种优化的误差准则，就是使未知量与已知量的均方误差达到最小化，在这种条件下来确定所需的未知量。
[0052]
本发明为全基因组选择的全新计算模型(cr-elastic net)，示例主要包括如下步骤：
[0053]
1)模拟数据采集
[0054]
模拟数据由qmsim2生成10000个个体，个体遗传力为0.1的50个性状。个体染体数目5条，染体长度为100cm，每条染体上包含1000个snp以及20个qtl，snps与qtl随机分布于染体上，qtl效应值满足指数分布。
[0055]
2)真实数据采集
[0056]
2.1半滑舌鳎数据集
[0057]
采用舌鳎芯片38ksnps标记位点，舌鳎抗哈维氏弧菌表型(遗传力0.16)，数据集包含a.1527尾半滑舌鳎抗哈维氏弧菌表型与基因型；b.44尾半滑舌鳎及其子代感染存活率。
[0058]
2.2猪数据集
[0059]
采用猪60ksnps标记位点，遗传力为0.07的性状t1，数据集个体数目为3534。
[0060]
3)模型预测效果验证
[0061]
利用模拟数据与真实数据的全部snp使用gblup与cr-elastic net进行交叉验证，并比较gblup与cr-elastic net模型预测准确性。
[0062]
本发明所提供的cr-elastic net可用于低遗传力下优质个体的筛选，预测准确性
相对于gblup有显著提升。最低有效参考体数量远低于gblup。因此可加快优质种质的选育进程。
[0063]
下面结合实例对本发明进行详细说明。
[0064]
实例1：模拟数据生成以及cr-elastic net与gblup模型预测准确性计算。(linux环境)
[0065]
1、模拟数据生成：
[0066]
模拟数据集：模拟遗传力为0.1。下图为示例，模拟体具备2个世代，每个世代10000个个体，雌雄比例为1：1,仅第一代有基因型。qtl数目100个，共5条染体，每条染体长度100cm,每条染体上随机分布1000个snp位点，qtl效应满足指数分布。重复50次。生成文件“lm_mrk{1..50}.txt”；“p1_mrk{1..50}.txt”；“p1_data{1..50}.txt”.
[0067]
通过r语言整理为plink(map；ped)格式：
[0068]
通过plink1.9将输出文件转为.raw(linux)；
[0069]
获得文件geno{1..50}.raw；pheno{1..50}.csv；
[0070]
julia语言cr-elasticnet理论与实现；
[0071]
整理好的基因型数据“geno{1..50}.raw”,表型数据“sim{1..50}pheno.csv”，验证集文件“cross.val.txt”,训练集文件julia:“cross.train{500；1000；2000}.csv”；
[0072]
r:验证集文件:“valname.bin”；训练集文集：“list_trainnames.bin”；
[0073]
1.1、导入依赖包；
[0074]
1.2、计算；
[0075]
1.3、r语言gblup模型与实现：
[0076]
y＝xb+za+e
[0077]
其中，y为一个n
×
1的表型向量，n是样本大小；x为固定效应，b为固定效应的向量；z作为随机效应是一个n
×
n的亲缘关系矩阵，e为残差，a为随机效应拟合值为一个n
×
1的向量。若数据为模拟数据则x为全1向量，若存在固定效应则为固定效应矩阵。
[0078]
r语言gblup使用r程序包“sommer”实现。
[0079]
2、有效参考体数量判断：
[0080]
从模拟数据10000个有基因型个体中随机抽取1000个个体作为验证集，从剩余个体中随机抽取{500，1000，2000，3000，4000，5000}个体作为训练集，重复20次，分别通过cr-elasticnet与gblup计算训练集与验证集gebv，计算训练集gebv与tbv(定义为qtls的和)的pearson相关性。因模拟数据为无家系体，pblup方法预测准确性为0。以模型(gblup；cr-elasticnet)劣于blup预测效果作为原假设，对50个模拟性状进行符号检验。置信水平为α＝0.05，使用bonferroni校正.如表1所示50个模拟数据集中cr-elastic net与gblup相对pblup有效提升预测准确性表型数。其统计结果图如图2所示。
[0081]
表1
[0082]
[0083]
由表1与图2可以得出结论，在cr-elastic net在低参考体数目(500)下稳定相对于pblup方法预测准确度提升(阳性结果》80％)，而gblup在3000个参考体数目以上稳定相对于pblup方法预测准确度提升(阳性结果》80％)。
[0084]
3.数据预测准确性对比：
[0085]
从模拟数据10000个有基因型个体中随机抽取1000个个体作为验证集，从剩余个体中随机抽取{500，1000，2000}个个体作为训练集，重复20次，分别通过cr-elastic net与gblup计算训练集与验证集gebv，计算验证集gebv与tbv(定义为qtls的和)的pearson相关性。其中将不存在预测效果的性状预测准确性定义为0，预测准确度定义为acc＝r(gebv,tbv)/h.将上述对50个性状的平均准确度提升(upgrade rate)如表2所示50个模拟数据集cr-elastic net相对gblup预测准确性提升。在表2所述数据中cr-elastic net相对于gblup在不同参考体数下提升超过5％，10％，20％的性状数目，如图3所示
[0086]
表2
[0087]
[0088]
[0089][0090]
通过表2与图3可以得出，在500，1000，2000个个体作为训练集下，cr-elastic net模型预测效果显著优于gblup，在500个体训练集下平均提升47.5％，其中42个表型预测准确性提升超过10％，36个表型预测准确性提升高于20％。在1000个体训练集下平均提升36.41％，其中43个表型预测准确性提升超过10％，38个表型预测准确性提升高于20％。在2000个体训练集下平均提升8.8％，其中16个表型预测准确性提升超过10％，6个表型预测准确性提升高于20％。
[0091]
示例2.半滑舌鳎数据集预测准确性对比。
[0092]
1、交叉验证：
[0093]
使用r语言包sommer所提供的gblup算法，结合整理好的基因型数据genotype.csv,表型数据phonetype.csv,系谱数据pedigree.csv以及交叉验证分组文件k-cross.csv对全基因组重测序的参考体和候选体共1527尾半滑舌鳎个体进行基因组选择计算，其中gblup存在固定效应，cr-elastic net将固定效应转化为哑变量。计算验证集gebv与表型值之间的pearson相关，如表3所示cr-elastic net与gblup对1572尾舌鳎交叉验证预测准确性表。对数据进行成对t检验，差异极显著(p＝0.001846056；双尾)，结果如图4所示。
[0094]
表3
[0095][0096]
2对子代存活率验证
[0097]
分别通过cr-elastic net与gblup估算44尾半滑舌鳎候选体的gebv，用这些个体繁育后，其gebv与子代存活率见表4所示44尾半滑舌鳎候选个体子代家系gebv和感染存活率表(gblup and cr-elastic net)。
[0098]
表4
[0099]
[0100][0101]
从表4可以计算出对子代gebv与感染存活率之间的pearson相关系数cr-elastic net(0.7936)》gblup(0.71),提升12％。由此可知cr-elastic net模型是一种相对准确的估计模型，
[0102]
示例3.对猪遗传力为0.07的性状t1预测准确性对比：
[0103]
从3534头猪60ksnps数据集中选取有表型t1的个体共2804头猪，从中随机抽取{500，1000，2000}头猪作为训练集，其余作为验证集，重复50次。使用训练集分别使用cr-elasticnet与gblup模型计算验证集gebv与表型的pearson相关性。其中将无预测效果的性状预测准确性定义为0。预测准确度定义为acc＝r(gebv,pheno)/h.将上述对50个重复的平均准确度以及相对提升(upgrade rate)，如表5所示猪遗传力为0.07性状t1二分法交叉验证预测准确性表。对cr-elastic net与gblup在使用不同参考体数目预测遗传力为0.07
性状的预测准确度对比，如图5所示。
[0104]
表5
[0105][0106]
从表5和图5中可以得出cr-elastic net模型相对于gblup在猪的低遗传力表型真实数据集中也有很显著的提升。进一步证明了模拟数据得到的结论。
[0107]
上述模拟数据示例与两个真实数据示例结果表明，本发明使用的cr-elastic net模型在低遗传力性状的无家系种质选育中具备优良的预测效果。在模拟数据的多个性状在500(平均提升47.4％)，1000(平均提升36.4％)，2000(平均提升8.9％)个体作为训练集下，cr-elastic net模型预测效果显著优于gblup，有效提升预测准确度的表型占比高且相对需要的有效参考体数目低。在真实舌鳎数据下交叉验证预测准确性(提升34％)，对后代的预测准确性也由(0.71提升至0.793)，效果极佳。在猪数据集下对低遗传力的验证也存在显著提升。两个真实数据集佐证了模拟数据的结论综上所述，cr-elastic net模型在常用参考体数目(500-2000个)对无家系，低遗传力的性状预测效果达到国际顶尖水平且对最低参考体数目要求很低(》500)。可以在后续种质选育中推广应用。
[0108]
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0109]
在本发明的描述中，需要理解的是，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。
[0110]
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

技术特征：

1.一种用于无家系低遗传力品种的基因组选择方法，其特征在于，包括以下步骤：基于snps的非等效条件，通过校正所述snps之间的相关性关系，评估snps效应值；根据所述snps效应值，对无家系低遗传力品种进行选育。2.根据权利要求1所述一种用于无家系低遗传力品种的选育方法，其特征在于：在评估所述snps效应值的过程中，根据所述非等效条件和所述相关性关系，构建用于评估所述snps效应值的cr-elastic net模型，通过设置超参数(假定snps之间的相关性系数)，获取常量微调惩罚量和模型代价函数，对所述snps效应值进行评估。3.根据权利要求2所述一种用于无家系低遗传力品种的选育方法，其特征在于：在构建cr-elastic net模型的过程中，所elastic net模型的表达式为：式中，w为snps效应值向量；w
t
为向量w的转置；||w||为w的转置共轭矩阵与矩阵w的积的最大特征根的平方根值，||w||＝w
t
w；为w的欧几里得范数的平方,为w里每个元素的平方和；y
i
为第i个观测的表型值；x
i
为第i个观测的基因型为一个n
×
1的向量，n为snps数量，λ为常量微调惩罚量，n为样本量，ρ为一个0-1之间的常数，其中，当ρ为0时模型代价函数等价于岭回归代价函数，ρ为1时模型代价函数等价于lasso回归代价函数。经验设定p值为0.3，即当snps之间的相关系数高于0.7(1-0.3)则认为其相关。4.根据权利要求3所述一种用于无家系低遗传力品种的选育方法，其特征在于：根据经验设定一个snps之间相关性的参数p，当snps之间相关性系数低于1-p则认为其相互独立，高于1-p则认为其相关。通过模型拟合设置在该相关性参数下模型解释残差比例最高的λ.在以上p与λ下，获取snps的效应值。5.根据权利要求4所述一种用于无家系低遗传力品种的选育方法，其特征在于：在评估所述snps效应值的过程中，通过重抽样，从总体样本中有放回抽取样本容量为n的m个子集，对m个子集通过弹性网络拟合，得到m组snps效应值w
k
，snps的分布服从原分布，snps真实效应值依概率收敛到{w
k
}的平均值。6.一种用于无家系低遗传力品种的选育应用，其特征在于，用于实现无家系低遗传力品种的选育系统，包括：第一分析模块，用于假设所述snps之间的相关性关系阈值；第二分析模块，用于获取在第一分析模块下最优的常量微调惩罚量；评估模块，用于评估所述snps效应值；选育模块，用于根据所述snps效应值，对无家系低遗传力品种进行选育。

技术总结

本发明公开了一种用于无家系低遗传力品种的基因组选择方法及应用，包括：基于SNPs的非等效条件，通过校正SNPs之间的相关性关系，评估SNPs效应值；根据SNPs效应值，对无家系低遗传力品种进行选育；本发明提供的用于评估SNPs效应值的方法可用于无家系低遗传力品种的选育，有利推动种业发展。有利推动种业发展。有利推动种业发展。