CN110610744B - 一种高效可并行运算且高准确性的全基因组选择方法 - Google Patents

一种高效可并行运算且高准确性的全基因组选择方法 Download PDF

Info

Publication number
CN110610744B
CN110610744B CN201910858021.XA CN201910858021A CN110610744B CN 110610744 B CN110610744 B CN 110610744B CN 201910858021 A CN201910858021 A CN 201910858021A CN 110610744 B CN110610744 B CN 110610744B
Authority
CN
China
Prior art keywords
value
optimal
mlm
calculating
flm
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
Application number
CN201910858021.XA
Other languages
English (en)
Other versions
CN110610744A (zh
Inventor
赵书红
尹立林
刘小磊
李新云
余梅
朱猛进
唐振双
许婧雅
殷东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan Shadow Gene Technology Co.,Ltd.
Original Assignee
Huazhong Agricultural University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Huazhong Agricultural University filed Critical Huazhong Agricultural University
Priority to CN201910858021.XA priority Critical patent/CN110610744B/zh
Publication of CN110610744A publication Critical patent/CN110610744A/zh
Application granted granted Critical
Publication of CN110610744B publication Critical patent/CN110610744B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics
    • G16B50/30—Data warehousing; Computing architectures

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Theoretical Computer Science (AREA)
  • Medical Informatics (AREA)
  • Biophysics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Genetics & Genomics (AREA)
  • Molecular Biology (AREA)
  • Analytical Chemistry (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Chemical & Material Sciences (AREA)
  • Bioethics (AREA)
  • Databases & Information Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及动植物育种及人类疾病预测技术领域，提供一种高效可并行运算且高准确性的全基因组选择方法。首先读取原始基因型文件和表型文件，构建新的基因型文件和表型文件，计算所有个体间的亲缘关系矩阵；然后提取新的表型文件中所有个体作为参考群，提取原始基因型文件中无表型数据的所有个体作为预测群；接着利用参考群数据进行全基因组关联分析，提取全基因组关联分析的结果特征；构建性状特异的模型库，采用交叉验证策略，依次优化最佳固定效应、最佳随机效应，从模型库中选取最优预测模型；最后，利用最优预测模型，计算预测群的基因组估计育种值。本发明能够快速准确且稳定地预测出个体基因组育种值，提升全基因组选择的准确性及效率。

Description

一种高效可并行运算且高准确性的全基因组选择方法
技术领域
本发明涉及动植物育种及人类疾病预测技术领域，特别是涉及一种高效可并行运算且高准确性的全基因组选择方法。
背景技术
随着覆盖整个基因组高密度单核苷酸多态性(SNP)基因分型技术的发展，全基因组选择(预测)作为基因组统计分析的强大工具，被广泛应用于植物和动物育种中复杂性状的遗传价值(种用价值)预测和评估，以及人类遗传学研究中。
现有的全基因组选择方法分为两类：一类是以全基因组最佳线性无偏预测GBLUP(Genomicbest linearunbiasedprediction)为代表的直接法，仅需构建个体间基因组关系矩阵，获取方差组分后可通过求解混合模型求得个体育种值；另一类以BayesB为代表的间接法，结合Bayes理论和隐马尔可夫迭代过程求取标记效应值，然后依照个体基因型对标记效应进行累加获得个体育种值。其中，直接法计算效率高，但由于其对性状遗传构建的假设简单，估计的育种值准确性较差；间接法对性状遗传构建的假设相对复杂，更符合性状遗传机制，具有更好的预测准确性，但由于其假设引入众多的未知参数，导致其参数求解过程极其复杂，计算效率较差，限制了间接法在实际预测中的应用。
发明内容
针对现有技术存在的问题，本发明提供一种高效可并行运算且高准确性的全基因组选择方法，能够快速准确且稳定地预测出个体基因组育种值，提升全基因组选择的准确性及计算效率。
本发明的技术方案为：
一种高效可并行运算且高准确性的全基因组选择方法，其特征在于，包括下述步骤：
步骤1：读取原始基因型文件和原始表型文件，提取原始基因型文件和原始表型文件中相同个体的基因型数据和表型数据，形成新的基因型文件和新的表型文件，并利用新的基因型文件计算所有个体间的亲缘关系矩阵G；
步骤2：提取新的表型文件中的所有个体作为参考群，提取原始基因型文件中无表型数据的所有个体作为预测群，得到参考群数据和预测群数据，将参考群随机分为M个规模相同的子参考群；其中，参考群数据包括参考群中每个个体的基因型数据和表型数据，预测群数据包括预测群中每个个体的基因型数据；
步骤3：利用参考群数据进行全基因组关联分析，提取全基因组关联分析的结果特征；构建性状特异的模型库，采用交叉验证策略，依次优化最佳固定效应、最佳随机效应，从模型库中选取最优预测模型；
步骤4：利用所述最优预测模型，计算所述预测群的基因组估计育种值。
所述步骤3包括下述步骤：
步骤3.1：使用L个线程并行执行相关系数计算；其中，使用第l∈{1,2,...,L}个线程执行相关系数计算包括：
步骤3.1.1：随机选取M-1个子参考群组合成测试群，将未被选取的子参考群作为验证群，利用测试群进行全基因组关联分析，采用选定的模型提取基因型中所有位点的显著P值为{P1l,P2l,...,Pnl,...,PNl}；其中，N为基因型中的位点总数，Pnl为第l个线程计算的基因型中第n个位点的显著P值，n∈{1,2,...,N}；
步骤3.1.2：用预设大小的窗口对所有位点按照基因组上的分布顺序进行划分，得到x个窗口内的位点；将每个窗口内的位点按照显著P值从小到大进行排序，选取每个窗口中显著P值最大的位点，形成位点集合X；
步骤3.1.3：利用测试群和GBLUP模型计算验证群的基因组估计育种值GEBV1，GEBV1只包含随机效应部分，计算GEBV1与验证群真实表型间的相关系数C0l；利用所述位点集合X测试固定效应模型FLM，从位点集合X中依次逐个不放回地取出位点并加入FLM中作为协变量，计算验证群的基因组估计育种值GEBV2，GEBV2只包含固定效应部分，并计算GEBV2与验证群真实表型间的相关系数集合为{Cf1l,Cf2l,...,Cfil,...,Cfxl}；利用所述位点集合X测试混合效应模型MLM，从位点集合X中依次逐个不放回地取出位点并加入MLM中作为协变量，计算验证群的基因组估计育种值GEBV3，GEBV3包含固定效应及随机效应两部分，并计算GEBV3与验证群真实表型间的相关系数集合为{Cm1l,Cm2l,...,Cmil,...,Cmxl}；其中，i∈{1,2,...,x}；
步骤3.1.4：若Cfil＞Cf,i-1,l且Cfil＞C0l，则位点集合X中第i个位点为FLM有效位点，第i个位点对应的窗口为FLM有效窗口，得到FLM有效窗口集合Fl；若Cmil＞Cm,i-1,l且Cmil＞C0l，则位点集合X中第i个位点为MLM有效位点，第i个位点对应的窗口为MLM有效窗口，得到MLM有效窗口集合Ml；
步骤3.2：计算{C01,C02,...,C0l,...,C0L}的均值为
Figure BDA0002198792600000031
计算{Cfi1,Cfi2,...,Cfil,...,CfiL}的均值为
Figure BDA0002198792600000032
{Cmi1,Cmi2,...,Cmil,...,CmiL}的均值为
Figure BDA0002198792600000033
i∈{1,2,...,x}，得到均值第一集合为
Figure BDA0002198792600000034
均值第二集合为
Figure BDA0002198792600000035
计算均值第一集合中元素的最大值为
Figure BDA0002198792600000036
均值第二集合中元素的最大值为
Figure BDA0002198792600000037
若
Figure BDA0002198792600000038
则选取最优预测模型为FLM；若
Figure BDA0002198792600000039
则选取最优预测模型为MLM；若
Figure BDA00021987926000000310
且
Figure BDA00021987926000000311
则选取最优预测模型为GBLUP；
步骤3.3：若最优预测模型为FLM，则对L个FLM有效窗口集合{F1,F2,...,Fl,...,FL}中的窗口进行计数，挑取出现次数大于或等于L×95％的FLM有效窗口作为终选FLM有效窗口；若最优预测模型为MLM，则对L个MLM有效窗口集合{M1,M2,...,Ml,...,ML}中的窗口进行计数，挑取出现次数大于或等于L×95％的MLM有效窗口作为终选MLM有效窗口；
步骤3.4：计算{Pn1,Pn2,...,Pnl,...,PnL}的指定值作为第n个位点的最终关联P值
Figure BDA00021987926000000312
得到所有位点的最终关联P值为
Figure BDA00021987926000000313
在终选FLM有效窗口中选取最终关联P值最大的位点作为FLM最佳协变量位点或在终选MLM有效窗口中选取最终关联P值最大的位点作为MLM最佳协变量位点；
步骤3.5：使用L个线程并行执行梯度下相关系数计算；其中，使用第l∈{1,2,...,L}个线程执行梯度下相关系数计算包括：
步骤3.5.1：若最优预测模型为GBLUP或MLM，则基于Vanraden算法初始化N×N的对角权重矩阵W＝diag(w1,w2,...,wN)＝diag(1,1,...,1)；对所述步骤3.1.1得到的所有位点的显著P值{P1l,P2l,...,Pnl,...,PNl}按照从小到大的顺序进行排序，得到排序后的显著P值序列{P1l',P2l',...,Pnl',...,PNl'}，将排序后的显著P值序列中前α％的元素对应的权重设置为放大倍数、将后(1-α％)的元素对应的权重保持不变仍为1，得到新的对角权重矩阵W'；其中，对α设置n1个梯度为{α1,α2,...,αp,...,αn1}、设置放大函数为logβP，对β设置n2个梯度{β1,β2,...,βq,...,βn2}，梯度βq下前αp％的元素中第k个元素对应的放大倍数为
Figure BDA0002198792600000041
结合Vanraden算法计算新的亲缘关系矩阵为T；若最佳预测模型为MLM，则加入所述MLM最佳协变量位点到MLM模型中；
步骤3.5.2：利用矩阵T及所述步骤3.1.1中的测试群计算所述步骤3.1.1中的验证群的基因组估计育种值GEBV4，并计算梯度{αp,βq}下GEBV4与所述步骤3.1.1中的验证群真实表型间的相关系数Cpql，继而得到所有梯度下的相关系数为{C11l,C12l,...,C1,n2,l,...,Cpql,Cp,q+1,l,...,Cp,n2,l,...,Cn1,1,l,Cn1,2,l,...,Cn1,n2,l}；
步骤3.6：计算{Cpq1,Cpq2,...,Cpql,...,CpqL}的均值为
Figure BDA0002198792600000042
得到均值第三集合为
Figure BDA0002198792600000043
计算均值第三集合中的最大值为
Figure BDA0002198792600000044
若
Figure BDA0002198792600000045
则无需对亲缘关系矩阵G进行优化；若
Figure BDA0002198792600000046
且{C11l,C12l,...,C1,n2,l,...,Cpql,Cp,q+1,l,...,Cp,n2,l,...,Cn1,1,l,Cn1,2,l,...,Cn1,n2,l|l∈{1,2,...,L}}中满足Cpql＞C0l的元素个数小于L×95％，则无需对亲缘关系矩阵G进行优化；若
Figure BDA0002198792600000047
且{C11l,C12l,...,C1,n2,l,...,Cpql,Cp,q+1,l,...,Cp,n2,l,...,Cn1,1,l,Cn1,2,l,...,Cn1,n2,l|l∈{1,2,...,L}}中满足Cpql＞C0l的元素个数大于或等于L×95％，则需要对亲缘关系矩阵G进行优化且最优前百分比为αp％、最优放大函数为
Figure BDA00021987926000000413
步骤3.7：若需要对亲缘关系矩阵G进行优化且最优前百分比为αp％、最优放大函数为
Figure BDA0002198792600000048
则将
Figure BDA0002198792600000049
按照最终关联P值从小到大排序，得到新的最终关联P值序列
Figure BDA00021987926000000410
对新的最终关联P值序列
Figure BDA00021987926000000411
中前αp％的元素中第k个元素对应的权重设置为放大倍数
Figure BDA00021987926000000412
后(1-αp％)的元素对应的权重不变且保持为1，形成对角权重矩阵，并结合Vanraden算法计算最佳性状特异的个体间亲缘关系矩阵G'。
所述步骤3.1.1中，所述选定的模型为GBLUP模型或固定效应模型FLM或混合效应模型MLM。
所述步骤3.4中，所述指定值为最大值或最小值或均值或中值。
所述步骤4包括下述步骤：
步骤4.1：若最优预测模型为FLM，则利用所述FLM最佳协变量位点及所述参考群计算所述预测群的基因组估计育种值；
步骤4.2：若最优预测模型为GBLUP且无需对亲缘关系矩阵G进行优化，则利用所述亲缘关系矩阵G及所述参考群计算所述预测群的基因组估计育种值；
步骤4.3：若最优预测模型为GBLUP且需要对亲缘关系矩阵G进行优化，则利用所述最佳性状特异的个体间亲缘关系矩阵G'及所述参考群计算所述预测群的基因组估计育种值；
步骤4.4：若最优预测模型为MLM且无需对亲缘关系矩阵G进行优化，则利用所述MLM最佳协变量位点、所述亲缘关系矩阵G及所述参考群计算所述预测群的基因组估计育种值；
步骤4.5：若最优预测模型为MLM且需要对亲缘关系矩阵G进行优化，则利用所述MLM最佳协变量位点、所述最佳性状特异的个体间亲缘关系矩阵G'及所述参考群计算所述预测群的基因组估计育种值。
本发明的有益效果为：
本发明在GBLUP模型的基础上融合Bayes理论的先验假设，结合全基因组关联分析(GWAS)信息，基于可并行的交叉验证，采用多元回归、网格搜索、二分法求极值等策略，针对不同性状自动选择最佳预测模型，准确筛选大效应标记以整合模型协变量，同时给予不同标记合适权重以构建个体关系矩阵，联合解析复杂性状遗传构建(GeneticArchitecture)，能够快速准确且稳定地预测出个体基因组育种值，提升全基因组选择的准确性及计算效率。
附图说明
图1为本发明的高效可并行运算且高准确性的全基因组选择方法的原理图。
具体实施方式
下面将结合附图和具体实施方式，对本发明作进一步描述。
如图1所示，本发明的高效可并行运算且高准确性的全基因组选择方法，其特征在于，包括下述步骤：
步骤1：读取原始基因型文件和原始表型文件，提取原始基因型文件和原始表型文件中相同个体的基因型数据和表型数据，形成新的基因型文件和新的表型文件，并利用新的基因型文件计算所有个体间的亲缘关系矩阵G。
本实施例中，采用R语言中的S4数据格式来建立磁盘与内存之间的数据映射；基因型文件和表型文件的读取及储存采用R CRAN::bigmemory软件包中的big.matrix格式，并提供plink格式[.ped/.map,.bed/.bim/.fam]、hapmap格式、VCF格式和numeric格式[0/1/2]四种数据格式的转换。检查、调整原始基因型文件和原始表型文件中的个体顺序，使其一致，并选留两个文件中同时存在的个体，获得新的表型文件及基因型文件。
步骤2：提取新的表型文件中的所有个体作为参考群，提取原始基因型文件中无表型数据的所有个体作为预测群，得到参考群数据和预测群数据，将参考群随机分为M个规模相同的子参考群；其中，参考群数据包括参考群中每个个体的基因型数据和表型数据，预测群数据包括预测群中每个个体的基因型数据。
步骤3：利用参考群数据进行全基因组关联分析，提取全基因组关联分析的结果特征；构建性状特异的模型库，采用交叉验证策略，依次优化最佳固定效应、最佳随机效应，从模型库中选取最优预测模型。
所述步骤3包括下述步骤：
步骤3.1：使用L个线程并行执行相关系数计算；其中，使用第l∈{1,2,...,L}个线程执行相关系数计算包括：
步骤3.1.1：随机选取M-1个子参考群组合成测试群，将未被选取的子参考群作为验证群，利用测试群进行全基因组关联分析，采用选定的模型提取基因型中所有位点的显著P值为{P1l,P2l,...,Pnl,...,PNl}；其中，N为基因型中的位点总数，Pnl为第l个线程计算的基因型中第n个位点的显著P值，n∈{1,2,...,N}；本实施例中，所述选定的模型为GBLUP模型或固定效应模型FLM或混合效应模型MLM；
步骤3.1.2：用预设大小的窗口对所有位点按照基因组上的分布顺序进行划分，得到x个窗口内的位点；将每个窗口内的位点按照显著P值从小到大进行排序，选取每个窗口中显著P值最大的位点，形成位点集合X；本实施例中，窗口大小为1MB；
步骤3.1.3：利用测试群和GBLUP模型计算验证群的基因组估计育种值GEBV1，GEBV1只包含随机效应部分，计算GEBV1与验证群真实表型间的相关系数C0l；利用所述位点集合X测试固定效应模型FLM(Fixed effect Linear Model)，从位点集合X中依次逐个不放回地取出位点并加入FLM中作为协变量，计算验证群的基因组估计育种值GEBV2，GEBV2只包含固定效应部分，并计算GEBV2与验证群真实表型间的相关系数集合为{Cf1l,Cf2l,...,Cfil,...,Cfxl}；利用所述位点集合X测试混合效应模型MLM(MixedeffectLinearModel)，从位点集合X中依次逐个不放回地取出位点并加入MLM中作为协变量，计算验证群的基因组估计育种值GEBV3，GEBV3包含固定效应及随机效应两部分，并计算GEBV3与验证群真实表型间的相关系数集合为{Cm1l,Cm2l,...,Cmil,...,Cmxl}；其中，i∈{1,2,...,x}；
步骤3.1.4：若Cfil＞Cf,i-1,l且Cfil＞C0l，则位点集合X中第i个位点为FLM有效位点，第i个位点对应的窗口为FLM有效窗口，得到FLM有效窗口集合Fl；若Cmil＞Cm,i-1,l且Cmil＞C0l，则位点集合X中第i个位点为MLM有效位点，第i个位点对应的窗口为MLM有效窗口，得到MLM有效窗口集合Ml。
步骤3.2：计算{C01,C02,...,C0l,...,C0L}的均值为
Figure BDA0002198792600000071
计算{Cfi1,Cfi2,...,Cfil,...,CfiL}的均值为
Figure BDA0002198792600000072
{Cmi1,Cmi2,...,Cmil,...,CmiL}的均值为
Figure BDA0002198792600000073
i∈{1,2,...,x}，得到均值第一集合为
Figure BDA0002198792600000074
均值第二集合为
Figure BDA0002198792600000075
计算均值第一集合中元素的最大值为
Figure BDA0002198792600000076
均值第二集合中元素的最大值为
Figure BDA0002198792600000077
若
Figure BDA0002198792600000078
则选取最优预测模型为FLM；若
Figure BDA0002198792600000079
则选取最优预测模型为MLM；若
Figure BDA00021987926000000710
且
Figure BDA00021987926000000711
则选取最优预测模型为GBLUP。
步骤3.3：若最优预测模型为FLM，则对L个FLM有效窗口集合{F1,F2,...,Fl,...,FL}中的窗口进行计数，挑取出现次数大于或等于L×95％的FLM有效窗口作为终选FLM有效窗口；若最优预测模型为MLM，则对L个MLM有效窗口集合{M1,M2,...,Ml,...,ML}中的窗口进行计数，挑取出现次数大于或等于L×95％的MLM有效窗口作为终选MLM有效窗口。
步骤3.4：计算{Pn1,Pn2,…,Pnl,…,PnL}的指定值作为第n个位点的最终关联P值
Figure BDA00021987926000000712
得到所有位点的最终关联P值为
Figure BDA00021987926000000713
在终选FLM有效窗口中选取最终关联P值最大的位点作为FLM最佳协变量位点或在终选MLM有效窗口中选取最终关联P值最大的位点作为MLM最佳协变量位点。其中，所述指定值为最大值或最小值或均值或中值。本实施例中，指定值为最大值，也即计算{Pn1,Pn2,…,Pnl,…,PnL}的最大值作为第n个位点的最终关联P值
Figure BDA00021987926000000714
步骤3.5：使用L个线程并行执行梯度下相关系数计算；其中，使用第l∈{1,2,…,L}个线程执行梯度下相关系数计算包括：
步骤3.5.1：若最优预测模型为GBLUP或MLM，则基于Vanraden算法初始化N×N的对角权重矩阵W＝diag(w1,w2,…,wN)＝diag(1,1,…,1)；对所述步骤3.1.1得到的所有位点的显著P值{P1l,P2l,…,Pnl,…,PNl}按照从小到大的顺序进行排序，得到排序后的显著P值序列{P1l',P2l',…,Pnl',...,PNl'}，将排序后的显著P值序列中前α％的元素对应的权重设置为放大倍数、将后(1-α％)的元素对应的权重保持不变仍为1，得到新的对角权重矩阵W'；其中，对α设置n1个梯度为{α1,α2,...,αp,...,αn1}、设置放大函数为logβP，对β设置n2个梯度{β1,β2,...,βq,...,βn2}，梯度βq下前αp％的元素中第k个元素对应的放大倍数为
Figure BDA0002198792600000081
结合Vanraden算法计算新的亲缘关系矩阵为T；若最佳预测模型为MLM，则加入所述MLM最佳协变量位点到MLM模型中；
步骤3.5.2：利用矩阵T及所述步骤3.1.1中的测试群计算所述步骤3.1.1中的验证群的基因组估计育种值GEBV4，并计算梯度{αp,βq}下GEBV4与所述步骤3.1.1中的验证群真实表型间的相关系数Cpql，继而得到所有梯度下的相关系数为{C11l,C12l,...,C1,n2,l,...,Cpql,Cp,q+1,l,...,Cp,n2,l,...,Cn1,1,l,Cn1,2,l,...,Cn1,n2,l}。
本实施例中，对α设置n1个梯度为{0.01,0.1,...,1,5}；对β设置n2个梯度{1.015,1.035,...,3,5,10}。
步骤3.6：计算{Cpq1,Cpq2,...,Cpql,...,CpqL}的均值为
Figure BDA0002198792600000082
得到均值第三集合为
Figure BDA0002198792600000083
计算均值第三集合中的最大值为
Figure BDA0002198792600000084
若
Figure BDA0002198792600000085
则无需对亲缘关系矩阵G进行优化；若
Figure BDA0002198792600000086
且{C11l,C12l,...,C1,n2,l,...,Cpql,Cp,q+1,l,...,Cp,n2,l,...,Cn1,1,l,Cn1,2,l,...,Cn1,n2,l|l∈{1,2,...,L}}中满足Cpql＞C0l的元素个数小于L×95％，则无需对亲缘关系矩阵G进行优化；若
Figure BDA0002198792600000087
且{C11l,C12l,...,C1,n2,l,...,Cpql,Cp,q+1,l,...,Cp,n2,l,...,Cn1,1,l,Cn1,2,l,...,Cn1,n2,l|l∈{1,2,...,L}}中满足Cpql＞C0l的元素个数大于或等于L×95％，则需要对亲缘关系矩阵G进行优化且最优前百分比为αp％、最优放大函数为logβq(P)。
步骤3.7：若需要对亲缘关系矩阵G进行优化且最优前百分比为αp％、最优放大函数为
Figure BDA0002198792600000091
则将
Figure BDA0002198792600000092
按照最终关联P值从小到大排序，得到新的最终关联P值序列
Figure BDA0002198792600000093
对新的最终关联P值序列
Figure BDA0002198792600000094
中前αp％的元素中第k个元素对应的权重设置为放大倍数
Figure BDA0002198792600000095
后(1-αp％)的元素对应的权重不变且保持为1，形成对角权重矩阵，并结合Vanraden算法计算最佳性状特异的个体间亲缘关系矩阵G'。
步骤4：利用所述最优预测模型，计算所述预测群的基因组估计育种值，具体包括下述步骤：
步骤4.1：若最优预测模型为FLM，则利用所述FLM最佳协变量位点及所述参考群计算所述预测群的基因组估计育种值；
步骤4.2：若最优预测模型为GBLUP且无需对亲缘关系矩阵G进行优化，则利用所述亲缘关系矩阵G及所述参考群计算所述预测群的基因组估计育种值；
步骤4.3：若最优预测模型为GBLUP且需要对亲缘关系矩阵G进行优化，则利用所述最佳性状特异的个体间亲缘关系矩阵G'及所述参考群计算所述预测群的基因组估计育种值；
步骤4.4：若最优预测模型为MLM且无需对亲缘关系矩阵G进行优化，则利用所述MLM最佳协变量位点、所述亲缘关系矩阵G及所述参考群计算所述预测群的基因组估计育种值；
步骤4.5：若最优预测模型为MLM且需要对亲缘关系矩阵G进行优化，则利用所述MLM最佳协变量位点、所述最佳性状特异的个体间亲缘关系矩阵G'及所述参考群计算所述预测群的基因组估计育种值。
显然，上述实施例仅仅是本发明的一部分实施例，而不是全部的实施例。上述实施例仅用于解释本发明，并不构成对本发明保护范围的限定。基于上述实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，也即凡在本申请的精神和原理之内所作的所有修改、等同替换和改进等，均落在本发明要求的保护范围内。

Claims (4)

1.一种高效可并行运算且高准确性的全基因组选择方法，其特征在于，包括下述步骤：
步骤1：读取原始基因型文件和原始表型文件，提取原始基因型文件和原始表型文件中相同个体的基因型数据和表型数据，形成新的基因型文件和新的表型文件，并利用新的基因型文件计算所有个体间的亲缘关系矩阵G；
步骤2：提取新的表型文件中的所有个体作为参考群，提取原始基因型文件中无表型数据的所有个体作为预测群，得到参考群数据和预测群数据，将参考群随机分为M个规模相同的子参考群；其中，参考群数据包括参考群中每个个体的基因型数据和表型数据，预测群数据包括预测群中每个个体的基因型数据；
步骤3：利用参考群数据进行全基因组关联分析，提取全基因组关联分析的结果特征；构建性状特异的模型库，采用交叉验证策略，依次优化最佳固定效应、最佳随机效应，从模型库中选取最优预测模型；
步骤4：利用所述最优预测模型，计算所述预测群的基因组估计育种值；
所述步骤3包括下述步骤：
步骤3.1：使用L个线程并行执行相关系数计算；其中，使用第l∈{1,2,...,L}个线程执行相关系数计算包括：
步骤3.1.1：随机选取M-1个子参考群组合成测试群，将未被选取的子参考群作为验证群，利用测试群进行全基因组关联分析，采用选定的模型提取基因型中所有位点的显著P值为{P1l,P2l,...,Pnl,...,PNl}；其中，N为基因型中的位点总数，Pnl为第l个线程计算的基因型中第n个位点的显著P值，n∈{1,2,...,N}；
步骤3.1.2：用预设大小的窗口对所有位点按照基因组上的分布顺序进行划分，得到x个窗口内的位点；将每个窗口内的位点按照显著P值从小到大进行排序，选取每个窗口中显著P值最大的位点，形成位点集合X；
步骤3.1.3：利用测试群和GBLUP模型计算验证群的基因组估计育种值GEBV1，GEBV1只包含随机效应部分，计算GEBV1与验证群真实表型间的相关系数C0l；利用所述位点集合X测试固定效应模型FLM，从位点集合X中依次逐个不放回地取出位点并加入FLM中作为协变量，计算验证群的基因组估计育种值GEBV2，GEBV2只包含固定效应部分，并计算GEBV2与验证群真实表型间的相关系数集合为{Cf1l,Cf2l,...,Cfil,...,Cfxl}；利用所述位点集合X测试混合效应模型MLM，从位点集合X中依次逐个不放回地取出位点并加入MLM中作为协变量，计算验证群的基因组估计育种值GEBV3，GEBV3包含固定效应及随机效应两部分，并计算GEBV3与验证群真实表型间的相关系数集合为{Cm1l,Cm2l,...,Cmil,...,Cmxl}；其中，i∈{1,2,...,x}；
步骤3.1.4：若Cfil＞Cf,i-1,l且Cfil＞C0l，则位点集合X中第i个位点为FLM有效位点，第i个位点对应的窗口为FLM有效窗口，得到FLM有效窗口集合Fl；若Cmil＞Cm,i-1,l且Cmil＞C0l，则位点集合X中第i个位点为MLM有效位点，第i个位点对应的窗口为MLM有效窗口，得到MLM有效窗口集合Ml；
步骤3.2：计算{C01,C02,...,C0l,...,C0L}的均值为
Figure FDA0002652121470000021
计算{Cfi1,Cfi2,...,Cfil,...,CfiL}的均值为
Figure FDA0002652121470000022
{Cmi1,Cmi2,...,Cmil,...,CmiL}的均值为
Figure FDA0002652121470000023
i∈{1,2,...,x}，得到均值第一集合为
Figure FDA0002652121470000024
均值第二集合为
Figure FDA0002652121470000025
计算均值第一集合中元素的最大值为
Figure FDA0002652121470000026
均值第二集合中元素的最大值为
Figure FDA0002652121470000027
若
Figure FDA0002652121470000028
则选取最优预测模型为FLM；若
Figure FDA0002652121470000029
则选取最优预测模型为MLM；若
Figure FDA00026521214700000210
且
Figure FDA00026521214700000211
则选取最优预测模型为GBLUP；
步骤3.3：若最优预测模型为FLM，则对L个FLM有效窗口集合{F1,F2,...,Fl,...,FL}中的窗口进行计数，挑取出现次数大于或等于L×95％的FLM有效窗口作为终选FLM有效窗口；若最优预测模型为MLM，则对L个MLM有效窗口集合{M1,M2,...,Ml,...,ML}中的窗口进行计数，挑取出现次数大于或等于L×95％的MLM有效窗口作为终选MLM有效窗口；
步骤3.4：计算{Pn1,Pn2,...,Pnl,...,PnL}的指定值作为第n个位点的最终关联P值
Figure FDA00026521214700000212
得到所有位点的最终关联P值为
Figure FDA00026521214700000213
在终选FLM有效窗口中选取最终关联P值最大的位点作为FLM最佳协变量位点或在终选MLM有效窗口中选取最终关联P值最大的位点作为MLM最佳协变量位点；
步骤3.5：使用L个线程并行执行梯度下相关系数计算；其中，使用第l∈{1,2,...,L}个线程执行梯度下相关系数计算包括：
步骤3.5.1：若最优预测模型为GBLUP或MLM，则基于Vanraden算法初始化N×N的对角权重矩阵W＝diag(w1,w2,...,wN)＝diag(1,1,...,1)；对所述步骤3.1.1得到的所有位点的显著P值{P1l,P2l,...,Pnl,...,PNl}按照从小到大的顺序进行排序，得到排序后的显著P值序列{P1l',P2l',...,Pnl',...,PNl'}，将排序后的显著P值序列中前α％的元素对应的权重设置为放大倍数、将后(1-α％)的元素对应的权重保持不变仍为1，得到新的对角权重矩阵W'；其中，对α设置n1个梯度为{α1,α2,...,αp,...,αn1}、设置放大函数为logβP，对β设置n2个梯度{β1,β2,...,βq,...,βn2}，梯度βq下前αp％的元素中第k个元素对应的放大倍数为
Figure FDA0002652121470000031
结合Vanraden算法计算新的亲缘关系矩阵为T；若最佳预测模型为MLM，则加入所述MLM最佳协变量位点到MLM模型中；
步骤3.5.2：利用矩阵T及所述步骤3.1.1中的测试群计算所述步骤3.1.1中的验证群的基因组估计育种值GEBV4，并计算梯度{αp,βq}下GEBV4与所述步骤3.1.1中的验证群真实表型间的相关系数Cpql，继而得到所有梯度下的相关系数为{C11l,C12l,...,C1,n2,l,...,Cpql,Cp,q+1,l,...,Cp,n2,l,...,Cn1,1,l,Cn1,2,l,...,Cn1,n2,l}；
步骤3.6：计算{Cpq1,Cpq2,...,Cpql,...,CpqL}的均值为
Figure FDA0002652121470000032
得到均值第三集合为
Figure FDA0002652121470000033
计算均值第三集合中的最大值为
Figure FDA0002652121470000034
若
Figure FDA0002652121470000035
则无需对亲缘关系矩阵G进行优化；若
Figure FDA0002652121470000036
且{C11l,C12l,...,C1,n2,l,...,Cpql,Cp,q+1,l,...,Cp,n2,l,...,Cn1,1,l,Cn1,2,l,...,Cn1,n2,l|l∈{1,2,...,L}}中满足Cpql＞C0l的元素个数小于L×95％，则无需对亲缘关系矩阵G进行优化；若
Figure FDA0002652121470000037
且{C11l,C12l,...,C1,n2,l,...,Cpql,Cp,q+1,l,...,Cp,n2,l,...,Cn1,1,l,Cn1,2,l,...,Cn1,n2,l|l∈{1,2,...,L}}中满足Cpql＞C0l的元素个数大于或等于L×95％，则需要对亲缘关系矩阵G进行优化且最优前百分比为αp％、最优放大函数为
Figure FDA0002652121470000038
步骤3.7：若需要对亲缘关系矩阵G进行优化且最优前百分比为αp％、最优放大函数为
Figure FDA0002652121470000041
则将
Figure FDA0002652121470000042
按照最终关联P值从小到大排序，得到新的最终关联P值序列
Figure FDA0002652121470000043
对新的最终关联P值序列
Figure FDA0002652121470000044
中前αp％的元素中第k个元素对应的权重设置为放大倍数
Figure FDA0002652121470000045
后(1-αp％)的元素对应的权重不变且保持为1，形成对角权重矩阵，并结合Vanraden算法计算最佳性状特异的个体间亲缘关系矩阵G'。
2.根据权利要求1所述的高效可并行运算且高准确性的全基因组选择方法，其特征在于，所述步骤3.1.1中，所述选定的模型为GBLUP模型或固定效应模型FLM或混合效应模型MLM。
3.根据权利要求1所述的高效可并行运算且高准确性的全基因组选择方法，其特征在于，所述步骤3.4中，所述指定值为最大值或最小值或均值或中值。
4.根据权利要求1所述的高效可并行运算且高准确性的全基因组选择方法，其特征在于，所述步骤4包括下述步骤：
步骤4.1：若最优预测模型为FLM，则利用所述FLM最佳协变量位点及所述参考群计算所述预测群的基因组估计育种值；
步骤4.2：若最优预测模型为GBLUP且无需对亲缘关系矩阵G进行优化，则利用所述亲缘关系矩阵G及所述参考群计算所述预测群的基因组估计育种值；
步骤4.3：若最优预测模型为GBLUP且需要对亲缘关系矩阵G进行优化，则利用所述最佳性状特异的个体间亲缘关系矩阵G'及所述参考群计算所述预测群的基因组估计育种值；
步骤4.4：若最优预测模型为MLM且无需对亲缘关系矩阵G进行优化，则利用所述MLM最佳协变量位点、所述亲缘关系矩阵G及所述参考群计算所述预测群的基因组估计育种值；
步骤4.5：若最优预测模型为MLM且需要对亲缘关系矩阵G进行优化，则利用所述MLM最佳协变量位点、所述最佳性状特异的个体间亲缘关系矩阵G'及所述参考群计算所述预测群的基因组估计育种值。
CN201910858021.XA 2019-09-11 2019-09-11 一种高效可并行运算且高准确性的全基因组选择方法 Active CN110610744B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910858021.XA CN110610744B (zh) 2019-09-11 2019-09-11 一种高效可并行运算且高准确性的全基因组选择方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910858021.XA CN110610744B (zh) 2019-09-11 2019-09-11 一种高效可并行运算且高准确性的全基因组选择方法

Publications (2)

Publication Number Publication Date
CN110610744A CN110610744A (zh) 2019-12-24
CN110610744B true CN110610744B (zh) 2020-10-23

Family

ID=68892703

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910858021.XA Active CN110610744B (zh) 2019-09-11 2019-09-11 一种高效可并行运算且高准确性的全基因组选择方法

Country Status (1)

Country Link
CN (1) CN110610744B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111524545B (zh) * 2020-04-30 2023-11-10 天津诺禾致源生物信息科技有限公司 全基因组选择育种的方法和装置
CN111883205B (zh) * 2020-07-14 2023-10-20 云南省烟草农业科学研究院 一种基于全基因组选择烟草有害成分预测方法及应用
CN111898807B (zh) * 2020-07-14 2024-02-27 云南省烟草农业科学研究院 一种基于全基因组选择烟叶产量预测方法及应用
CN111798920B (zh) * 2020-07-14 2023-10-20 云南省烟草农业科学研究院 基于全基因组选择烟草经济性状表型值预测方法及应用
CN113517020A (zh) * 2021-08-04 2021-10-19 华中农业大学 一种快速准确的动物基因组选配分析方法
CN115579057A (zh) * 2022-09-15 2023-01-06 中国科学院植物研究所 基于全基因组选择研究的水稻籽粒镉积累性状预测装置和预警系统
CN116469466B (zh) * 2023-04-11 2024-02-09 南京农业大学 一种高效预测菊花耐涝性的方法及其应用
CN117238363B (zh) * 2023-10-25 2024-04-16 青岛极智医学检验实验室有限公司 一种表型预测方法、预测系统、设备及介质
CN117831636B (zh) * 2024-03-04 2024-06-11 北京市农林科学院信息技术研究中心 利用融合模型实施基因组选择的方法、装置、设备及介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107446997A (zh) * 2017-07-21 2017-12-08 河北农业大学 与陆地棉纤维细度关联的snp分子标记及其应用
CN107619873A (zh) * 2017-08-30 2018-01-23 上海交通大学 基于关联分析和KASP开发waxy1基因内分子标记

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9672271B2 (en) * 2013-05-17 2017-06-06 Lawrence Sirovich Method for identifying and employing high risk genomic markers for the prediction of specific diseases
CN103914631A (zh) * 2014-02-26 2014-07-09 中国农业大学 一种基于snp芯片的综合基因组育种值估计方法及应用
CN107200775A (zh) * 2017-01-20 2017-09-26 华中农业大学 一种提高水稻柱头外露率的方法
CN108728574B (zh) * 2018-06-21 2021-05-28 贵州省油菜研究所 甘蓝型油菜籽粒密度性状的主效qtl位点、snp分子标记及应用
CN110211640B (zh) * 2019-06-05 2023-04-07 南通大学 一种基于gpu并行计算的复杂疾病基因互作关联分析方法
CN110211635B (zh) * 2019-06-12 2020-07-21 北京康普森农业科技有限公司 用于畜禽基因组选择分析的方法及畜禽育种方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107446997A (zh) * 2017-07-21 2017-12-08 河北农业大学 与陆地棉纤维细度关联的snp分子标记及其应用
CN107619873A (zh) * 2017-08-30 2018-01-23 上海交通大学 基于关联分析和KASP开发waxy1基因内分子标记

Also Published As

Publication number Publication date
CN110610744A (zh) 2019-12-24

Similar Documents

Publication Publication Date Title
CN110610744B (zh) 一种高效可并行运算且高准确性的全基因组选择方法
CN110211635B (zh) 用于畜禽基因组选择分析的方法及畜禽育种方法
Iwata et al. Bayesian association mapping of multiple quantitative trait loci and its application to the analysis of genetic variation among Oryza sativa L. germplasms
CN111524545B (zh) 全基因组选择育种的方法和装置
US20020119451A1 (en) System and method for predicting chromosomal regions that control phenotypic traits
Gianola et al. Genome-wide association studies with a genomic relationship matrix: a case study with wheat and arabidopsis
CN114360651A (zh) 一种基因组预测方法、预测系统及应用
CN115359845A (zh) 一种融合单细胞转录组的空间转录组生物组织亚结构解析方法
CN116580773A (zh) 基于集成学习的育种跨代表型预测方法与系统、电子设备
US20070172833A1 (en) Gene expression profile retrieving apparatus, gene expression profile retrieving method, and program
CN113742070A (zh) 一种低深度测序群体基因型填充计算内存优化方法
Hassani et al. Accuracy of prediction of simulated polygenic phenotypes and their underlying quantitative trait loci genotypes using real or imputed whole-genome markers in cattle
Jung et al. Identifying Differentially Expressed Genes in Meta‐Analysis via Bayesian Model‐Based Clustering
US6994965B2 (en) Method for displaying results of hybridization experiment
CN112102880A (zh) 品种鉴定的方法、其预测模型的构建方法和装置
Van Dongen et al. One-and two-sample tests for single-locus inbreeding coefficients using the bootstrap
CN115995262B (zh) 基于随机森林及lasso回归解析玉米遗传机理的方法
Alanoshahr et al. The impact of different genetic architectures on accuracy of genomic selection using three Bayesian methods
He et al. MINED: an efficient mutual information based epistasis detection method to improve quantitative genetic trait prediction
CN112837750B (zh) 一种动植物育种值预测方法及装置
He et al. Muse: A multi-locus sampling-based epistasis algorithm for quantitative genetic trait prediction
CN109817340B (zh) 疾病风险分布信息确定方法、装置、存储介质及设备
Tchakounte-Wakem A Comparison of Methods Taking into Account Asymmetry when Evaluating Differential Expression in Gene Expression Experiments
KUKUČKOVÁ et al. Genetic markers and biostatistical methods as appropriate tools to preserve genetic resources
CN112365926A (zh) 一种玉米自交系一般配合力和特殊配合力的全基因组预测方法及应用

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
TR01 Transfer of patent right

Effective date of registration: 20221018

Address after: 4th Floor, Building D5, Phase III, Wuhan Software New Town, Donghu New Technology Development Zone, Wuhan, Hubei

Patentee after: Wuhan Shadow Gene Technology Co.,Ltd.

Address before: 430070 No. 1 Lion Rock street, Hongshan District, Hubei, Wuhan

Patentee before: HUAZHONG AGRICULTURAL University