CN110610744A - 一种高效可并行运算且高准确性的全基因组选择方法 - Google Patents
一种高效可并行运算且高准确性的全基因组选择方法 Download PDFInfo
- Publication number
- CN110610744A CN110610744A CN201910858021.XA CN201910858021A CN110610744A CN 110610744 A CN110610744 A CN 110610744A CN 201910858021 A CN201910858021 A CN 201910858021A CN 110610744 A CN110610744 A CN 110610744A
- Authority
- CN
- China
- Prior art keywords
- value
- optimal
- mlm
- calculating
- site
- 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.)
- Granted
Links
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}的均值为计算{Cfi1,Cfi2,...,Cfil,...,CfiL}的均值为{Cmi1,Cmi2,...,Cmil,...,CmiL}的均值为i∈{1,2,...,x},得到均值第一集合为均值第二集合为计算均值第一集合中元素的最大值为均值第二集合中元素的最大值为若则选取最优预测模型为FLM;若则选取最优预测模型为MLM;若且则选取最优预测模型为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值得到所有位点的最终关联P值为在终选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个元素对应的放大倍数为结合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}的均值为得到均值第三集合为计算均值第三集合中的最大值为若则无需对亲缘关系矩阵G进行优化;若且{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进行优化;若且{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%、最优放大函数为
步骤3.7:若需要对亲缘关系矩阵G进行优化且最优前百分比为αp%、最优放大函数为则将按照最终关联P值从小到大排序,得到新的最终关联P值序列对新的最终关联P值序列中前αp%的元素中第k个元素对应的权重设置为放大倍数后(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}的均值为计算{Cfi1,Cfi2,...,Cfil,...,CfiL}的均值为{Cmi1,Cmi2,...,Cmil,...,CmiL}的均值为i∈{1,2,...,x},得到均值第一集合为均值第二集合为计算均值第一集合中元素的最大值为均值第二集合中元素的最大值为若则选取最优预测模型为FLM;若则选取最优预测模型为MLM;若且则选取最优预测模型为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值得到所有位点的最终关联P值为在终选FLM有效窗口中选取最终关联P值最大的位点作为FLM最佳协变量位点或在终选MLM有效窗口中选取最终关联P值最大的位点作为MLM最佳协变量位点。其中,所述指定值为最大值或最小值或均值或中值。本实施例中,指定值为最大值,也即计算{Pn1,Pn2,…,Pnl,…,PnL}的最大值作为第n个位点的最终关联P值
步骤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个元素对应的放大倍数为结合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}的均值为得到均值第三集合为计算均值第三集合中的最大值为若则无需对亲缘关系矩阵G进行优化;若且{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进行优化;若且{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%、最优放大函数为则将按照最终关联P值从小到大排序,得到新的最终关联P值序列对新的最终关联P值序列中前αp%的元素中第k个元素对应的权重设置为放大倍数后(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 (5)
1.一种高效可并行运算且高准确性的全基因组选择方法,其特征在于,包括下述步骤:
步骤1:读取原始基因型文件和原始表型文件,提取原始基因型文件和原始表型文件中相同个体的基因型数据和表型数据,形成新的基因型文件和新的表型文件,并利用新的基因型文件计算所有个体间的亲缘关系矩阵G;
步骤2:提取新的表型文件中的所有个体作为参考群,提取原始基因型文件中无表型数据的所有个体作为预测群,得到参考群数据和预测群数据,将参考群随机分为M个规模相同的子参考群;其中,参考群数据包括参考群中每个个体的基因型数据和表型数据,预测群数据包括预测群中每个个体的基因型数据;
步骤3:利用参考群数据进行全基因组关联分析,提取全基因组关联分析的结果特征;构建性状特异的模型库,采用交叉验证策略,依次优化最佳固定效应、最佳随机效应,从模型库中选取最优预测模型;
步骤4:利用所述最优预测模型,计算所述预测群的基因组估计育种值。
2.根据权利要求1所述的高效可并行运算且高准确性的全基因组选择方法,其特征在于,所述步骤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}的均值为计算{Cfi1,Cfi2,...,Cfil,...,CfiL}的均值为{Cmi1,Cmi2,...,Cmil,...,CmiL}的均值为i∈{1,2,...,x},得到均值第一集合为均值第二集合为计算均值第一集合中元素的最大值为均值第二集合中元素的最大值为若则选取最优预测模型为FLM;若则选取最优预测模型为MLM;若且则选取最优预测模型为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值得到所有位点的最终关联P值为在终选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个元素对应的放大倍数为结合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}的均值为得到均值第三集合为计算均值第三集合中的最大值为若则无需对亲缘关系矩阵G进行优化;若且{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进行优化;若且{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%、最优放大函数为
步骤3.7:若需要对亲缘关系矩阵G进行优化且最优前百分比为αp%、最优放大函数为则将按照最终关联P值从小到大排序,得到新的最终关联P值序列对新的最终关联P值序列中前αp%的元素中第k个元素对应的权重设置为放大倍数后(1-αp%)的元素对应的权重不变且保持为1,形成对角权重矩阵,并结合Vanraden算法计算最佳性状特异的个体间亲缘关系矩阵G'。
3.根据权利要求2所述的高效可并行运算且高准确性的全基因组选择方法,其特征在于,所述步骤3.1.1中,所述选定的模型为GBLUP模型或固定效应模型FLM或混合效应模型MLM。
4.根据权利要求2所述的高效可并行运算且高准确性的全基因组选择方法,其特征在于,所述步骤3.4中,所述指定值为最大值或最小值或均值或中值。
5.根据权利要求2所述的高效可并行运算且高准确性的全基因组选择方法,其特征在于,所述步骤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'及所述参考群计算所述预测群的基因组估计育种值。
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 true CN110610744A (zh) | 2019-12-24 |
CN110610744B 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) |
Cited By (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111524545A (zh) * | 2020-04-30 | 2020-08-11 | 天津诺禾致源生物信息科技有限公司 | 全基因组选择育种的方法和装置 |
CN111798920A (zh) * | 2020-07-14 | 2020-10-20 | 云南省烟草农业科学研究院 | 基于全基因组选择烟草经济性状表型值预测方法及应用 |
CN111883205A (zh) * | 2020-07-14 | 2020-11-03 | 云南省烟草农业科学研究院 | 一种基于全基因组选择烟草有害成分预测方法及应用 |
CN111898807A (zh) * | 2020-07-14 | 2020-11-06 | 云南省烟草农业科学研究院 | 一种基于全基因组选择烟叶产量预测方法及应用 |
CN113517020A (zh) * | 2021-08-04 | 2021-10-19 | 华中农业大学 | 一种快速准确的动物基因组选配分析方法 |
CN116469466A (zh) * | 2023-04-11 | 2023-07-21 | 南京农业大学 | 一种高效预测菊花耐涝性的方法及其应用 |
CN117238363A (zh) * | 2023-10-25 | 2023-12-15 | 青岛极智医学检验实验室有限公司 | 一种表型预测方法、预测系统、设备及介质 |
WO2024056056A1 (zh) * | 2022-09-15 | 2024-03-21 | 中国科学院植物研究所 | 基于全基因组选择研究的水稻籽粒镉积累性状预测装置和预警系统 |
CN117831636A (zh) * | 2024-03-04 | 2024-04-05 | 北京市农林科学院信息技术研究中心 | 利用融合模型实施基因组选择的方法、装置、设备及介质 |
Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103914631A (zh) * | 2014-02-26 | 2014-07-09 | 中国农业大学 | 一种基于snp芯片的综合基因组育种值估计方法及应用 |
CN107200775A (zh) * | 2017-01-20 | 2017-09-26 | 华中农业大学 | 一种提高水稻柱头外露率的方法 |
CN107446997A (zh) * | 2017-07-21 | 2017-12-08 | 河北农业大学 | 与陆地棉纤维细度关联的snp分子标记及其应用 |
CN107619873A (zh) * | 2017-08-30 | 2018-01-23 | 上海交通大学 | 基于关联分析和KASP开发waxy1基因内分子标记 |
US20180046698A1 (en) * | 2013-05-17 | 2018-02-15 | Lawrence Sirovich | Intrinsic chromosomal linkage and disease prediction |
CN108728574A (zh) * | 2018-06-21 | 2018-11-02 | 贵州省油菜研究所 | 甘蓝型油菜籽粒密度性状的主效qtl位点、snp分子标记及应用 |
CN110211640A (zh) * | 2019-06-05 | 2019-09-06 | 南通大学 | 一种基于gpu并行计算的复杂疾病基因互作关联分析方法 |
CN110211635A (zh) * | 2019-06-12 | 2019-09-06 | 北京康普森农业科技有限公司 | 用于畜禽基因组选择分析的方法及畜禽育种方法 |
-
2019
- 2019-09-11 CN CN201910858021.XA patent/CN110610744B/zh active Active
Patent Citations (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20180046698A1 (en) * | 2013-05-17 | 2018-02-15 | Lawrence Sirovich | Intrinsic chromosomal linkage and disease prediction |
CN103914631A (zh) * | 2014-02-26 | 2014-07-09 | 中国农业大学 | 一种基于snp芯片的综合基因组育种值估计方法及应用 |
CN107200775A (zh) * | 2017-01-20 | 2017-09-26 | 华中农业大学 | 一种提高水稻柱头外露率的方法 |
CN107446997A (zh) * | 2017-07-21 | 2017-12-08 | 河北农业大学 | 与陆地棉纤维细度关联的snp分子标记及其应用 |
CN107619873A (zh) * | 2017-08-30 | 2018-01-23 | 上海交通大学 | 基于关联分析和KASP开发waxy1基因内分子标记 |
CN108728574A (zh) * | 2018-06-21 | 2018-11-02 | 贵州省油菜研究所 | 甘蓝型油菜籽粒密度性状的主效qtl位点、snp分子标记及应用 |
CN110211640A (zh) * | 2019-06-05 | 2019-09-06 | 南通大学 | 一种基于gpu并行计算的复杂疾病基因互作关联分析方法 |
CN110211635A (zh) * | 2019-06-12 | 2019-09-06 | 北京康普森农业科技有限公司 | 用于畜禽基因组选择分析的方法及畜禽育种方法 |
Non-Patent Citations (3)
Title |
---|
J314159的个人博客: "GWAS和Genomic prediction概念、原理及应用", 《HTTP://BLOG.SCIENCENET.CN/BLOG-292064-846224.HTML》 * |
刘小磊: "一种交替运用固定效应和随机效应模型优化全基因组关联分析的算法开发", 《中国博士学位论文全文数据库基础科学辑》 * |
尹立林: "全基因组选择模型研究进展及展望", 《畜牧兽医学报》 * |
Cited By (16)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111524545B (zh) * | 2020-04-30 | 2023-11-10 | 天津诺禾致源生物信息科技有限公司 | 全基因组选择育种的方法和装置 |
CN111524545A (zh) * | 2020-04-30 | 2020-08-11 | 天津诺禾致源生物信息科技有限公司 | 全基因组选择育种的方法和装置 |
CN111798920A (zh) * | 2020-07-14 | 2020-10-20 | 云南省烟草农业科学研究院 | 基于全基因组选择烟草经济性状表型值预测方法及应用 |
CN111883205A (zh) * | 2020-07-14 | 2020-11-03 | 云南省烟草农业科学研究院 | 一种基于全基因组选择烟草有害成分预测方法及应用 |
CN111898807A (zh) * | 2020-07-14 | 2020-11-06 | 云南省烟草农业科学研究院 | 一种基于全基因组选择烟叶产量预测方法及应用 |
CN111898807B (zh) * | 2020-07-14 | 2024-02-27 | 云南省烟草农业科学研究院 | 一种基于全基因组选择烟叶产量预测方法及应用 |
CN111883205B (zh) * | 2020-07-14 | 2023-10-20 | 云南省烟草农业科学研究院 | 一种基于全基因组选择烟草有害成分预测方法及应用 |
CN111798920B (zh) * | 2020-07-14 | 2023-10-20 | 云南省烟草农业科学研究院 | 基于全基因组选择烟草经济性状表型值预测方法及应用 |
CN113517020A (zh) * | 2021-08-04 | 2021-10-19 | 华中农业大学 | 一种快速准确的动物基因组选配分析方法 |
WO2024056056A1 (zh) * | 2022-09-15 | 2024-03-21 | 中国科学院植物研究所 | 基于全基因组选择研究的水稻籽粒镉积累性状预测装置和预警系统 |
CN116469466B (zh) * | 2023-04-11 | 2024-02-09 | 南京农业大学 | 一种高效预测菊花耐涝性的方法及其应用 |
CN116469466A (zh) * | 2023-04-11 | 2023-07-21 | 南京农业大学 | 一种高效预测菊花耐涝性的方法及其应用 |
CN117238363A (zh) * | 2023-10-25 | 2023-12-15 | 青岛极智医学检验实验室有限公司 | 一种表型预测方法、预测系统、设备及介质 |
CN117238363B (zh) * | 2023-10-25 | 2024-04-16 | 青岛极智医学检验实验室有限公司 | 一种表型预测方法、预测系统、设备及介质 |
CN117831636A (zh) * | 2024-03-04 | 2024-04-05 | 北京市农林科学院信息技术研究中心 | 利用融合模型实施基因组选择的方法、装置、设备及介质 |
CN117831636B (zh) * | 2024-03-04 | 2024-06-11 | 北京市农林科学院信息技术研究中心 | 利用融合模型实施基因组选择的方法、装置、设备及介质 |
Also Published As
Publication number | Publication date |
---|---|
CN110610744B (zh) | 2020-10-23 |
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 | |
CN116895334A (zh) | 用于估算或预测基因型和表型的方法和组成 | |
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) | 品种鉴定的方法、其预测模型的构建方法和装置 | |
Limpiti et al. | PCA-based informative SNP selection for analyzing population structure | |
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 | |
He et al. | Muse: A multi-locus sampling-based epistasis algorithm for quantitative genetic trait prediction | |
KUKUČKOVÁ et al. | Genetic markers and biostatistical methods as appropriate tools to preserve genetic resources | |
Tchakounte-Wakem | A Comparison of Methods Taking into Account Asymmetry when Evaluating Differential Expression in Gene Expression Experiments | |
Kawaguchi | Variable ranking by random forests model for genome-wide association study | |
KR20220099229A (ko) | 냉증 분류 방법 및 냉증 분류 시스템 |
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: 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 |
|
TR01 | Transfer of patent right |