A Simulation Method for Predicting Genotypes of Advanced Generations of Hybrid Combinations in Self-Pollinated Crops

Abstract

The present invention discloses a method for predicting and simulating the genotypes of advanced generations of hybrid combinations of self-pollinated crops, comprising the following steps: constructing a parental genotype data set and screening effective marker loci; calculating the linkage disequilibrium coefficient between marker loci, constructing a genomic linkage map, and determining the recombination rate between loci; using a random simulation algorithm to construct a gamete formation model and generate the genotypes of F1 generation individuals; based on the genotypes of the F1 generation, randomly simulating the selfing process to generate genotype data from F2 to Fn generations; establishing a mapping relationship between genotypes and phenotypes using a machine learning model to predict the phenotypic traits of advanced-generation selfing lines; the method of the present invention takes into account chromosomal exchange and linkage disequilibrium, can accurately simulate the gene recombination and segregation phenomena of the hybrid offspring of self-pollinated crops, and combines the gradient boosting decision tree model to achieve high-precision phenotypic prediction, providing an efficient auxiliary tool for the breeding of self-pollinated crops, which can significantly shorten the breeding cycle and improve the breeding efficiency.

Description

Technical Field

[0001] The present invention belongs to the technical field of cross-breeding, and particularly relates to a method for predicting and simulating the genotypes of high-generation hybrid combinations of self-pollinated crops. Background Art

[0002] Self-pollinated crops are an important type in agricultural production, and the selection of their hybrid combinations is a key means to improve crop yield and quality. Traditional methods for selecting self-pollinated crop hybrid combinations mainly rely on multi-generation self-pollination and field phenotypic selection, which require a large amount of time, manpower, and material resources. With the development of molecular biology and genomics technologies, molecular-assisted breeding methods based on DNA markers have gradually been applied to crop improvement, but there are still challenges in predicting the genotypes of high-generation inbred lines.

[0003] Existing prediction methods usually analyze based on actual sequencing data, cannot accurately simulate the gene recombination and segregation phenomena during the segregation process of hybrid offspring, lack consideration of chromosome exchange and linkage disequilibrium, resulting in low prediction accuracy. In addition, traditional methods cannot effectively predict the change trends of phenotypic traits after multi-generation self-pollination, making it impossible to optimize the allocation of resources in the breeding process.

[0004] Therefore, developing a method that can accurately simulate the genotypes of high-generation hybrid combinations of self-pollinated crops and predict trait performances is of great significance for improving breeding efficiency and shortening the breeding cycle. Summary of the Invention

[0005] The present invention provides a method for predicting and simulating the genotypes of high-generation hybrid combinations of self-pollinated crops. By constructing a parental genotype dataset of the self-pollinated crop hybrid combination, calculating the linkage disequilibrium coefficient between marker loci, constructing a gamete formation model using a random simulation algorithm, and establishing a mapping relationship between genotypes and phenotypes using a machine learning model, the purpose of accurately predicting the genotypes and phenotypic traits of high-generation inbred lines, improving breeding efficiency, and shortening the breeding cycle is achieved.

[0006] To achieve the above invention purpose, the specific technical solutions are as follows:

[0007] A method for predicting and simulating the genotypes of high-generation hybrid combinations of self-pollinated crops, the method comprising the following steps:

[0008] Step S1, constructing a parental genotype dataset of the self-pollinated crop hybrid combination, screening genomic marker loci for the parental genotype dataset, and obtaining effective marker loci;

[0009] Step S2, based on the parental genotype data, calculating the linkage disequilibrium coefficient between marker loci, constructing a genomic linkage map, and determining the recombination rate between loci;

[0010] Step S3, construct a gamete formation model using a random simulation algorithm, calculate the chromosome crossover probability according to the recombination rate, and generate the genotypes of the individuals in the

[0011] Step S4, based on the genotypes of the individuals in the generation, generate the genotype data of the individuals in the to generation through random simulation of the self-crossing process;

[0012] Step S5, establish a mapping relationship between genotypes and phenotypes using a machine learning model to predict the target trait performance of high-generation self-cross lines.

[0013] Further, in the said Step S1, the specific steps for screening genomic marker loci from the parental genotype dataset are:

[0014] Step S11, calculate the polymorphism information content value of each marker locus:

[0015] , where is the frequency of the th allele, and is the number of alleles;

[0016] Step S12, calculate the missing rate and minor allele frequency of the marker locus, and remove the marker locus with a missing rate greater than or a minor allele frequency less than ;

[0017] Step S13, sort the marker loci according to the genomic physical position, and when the distance between adjacent marker loci is less than the preset threshold, retain the marker locus with a higher PIC value;

[0018] The said preset threshold depends on the self-crossing crop species.

[0019] Further, calculating the linkage disequilibrium coefficient between marker loci in the said Step S2 includes:

[0020] Step S21, calculate the linkage disequilibrium coefficient and between loci : , where is the frequency of a specific allele combination at the two loci, and and are the frequencies of the respective alleles at loci and respectively;

[0021] Step S22, calculate the standardized linkage disequilibrium coefficient : ;

[0022] Step S23, construct a linkage map based on the value. When , it is considered that the two loci are strongly linked. When , it is considered moderately linked. When , it is considered weakly linked or unlinked.

[0023] Furthermore, the calculation formula for determining the recombination rate between loci in step S2 is:

[0024] , where is the recombination rate between locus and . is the physical distance between locus and . is the average recombination length of the inbred crop chromosome, which is adaptively adjusted according to the crop variety, and the value range is in cM.

[0025] Furthermore, the gamete formation model constructed by the random simulation algorithm in step S3 includes:

[0026] Step S31, based on the parental genotypes, randomly generate homologous chromosome exchange sites, and the probability of the exchange sites appearing is proportional to the recombination rate between the loci;

[0027] Step S32, simulate the gamete genotypes after chromosome exchange. For two parents and , for locus , if locus is within an odd exchange interval, the genotype of the hybrid individual at this locus comes from parent , expressed as: , if locus is within an even exchange interval, the genotype of the hybrid individual at this locus comes from parent , expressed as: , where and are the genotypes of parent and parent at locus ;

[0028] Step S33, randomly select the simulated parental gamete combinations to form generation hybrid individuals. The genotype of the individual is calculated as: , where and The genotypes of the paternal and maternal gametes at locus k, respectively.

[0029] Furthermore, the specific steps for generating high-generation individuals through random simulation of the self-crossing process in step S4 are as follows:

[0030] Step S41, perform self-crossing on individuals, simulate chromosome exchange during meiosis to produce gametes;

[0031] Step S42, randomly combine the produced gametes to form a population of individuals, the genotype of an individual at locus is calculated as: , where and are the genotypes of two gametes at locus ;

[0032] Step S43, continue to perform self-crossing on individuals, and similarly generate the genotypes of the generation, and so on until the

[0033] Step S44, calculate the homozygosity of each generation: , where is the total number of marker loci, is an indicator function that takes the value 1 when locus is homozygous and 0 otherwise; and represent the two states when locus is homozygous.

[0034] Furthermore, the establishment of the mapping relationship between genotypes and phenotypes using a machine learning model in step S5 includes:

[0035] Step S51, construct a training data set containing sample data with known genotypes and phenotypes, and divide the data into a training set and a validation set in a ratio of 7:3;

[0036] Step S52, use the Gradient Boosting Decision Tree model GBDT as the core prediction model. This model consists of multiple decision trees, and the calculation process of the model parameters is as follows:

[0037] Initialize the model , where, is the loss function; represents finding the value that minimizes the following expression, is the predicted value, is the true phenotypic value of the th sample;

[0038] For , is the number of decision trees, and is the residual of the th sample in the th iteration:

[0039] , is the partial derivative of the loss function with respect to the current predicted value, is the model after the (m - 1) th iteration, and is the genotype data of the th sample; Fit the residual to generate a decision tree to obtain the th leaf node region in the

[0040] th iteration; Calculate the optimal output value of each leaf node: ; Update the model: where is the learning rate, and its value range is

[0041] Step S53, the model optimization adopts -fold cross-validation method, , and determine the optimal hyperparameters through grid search, including the learning rate , the maximum depth of the tree, the number of decision trees and the minimum number of samples in the leaf node

[0042] Step S54, based on the high-generation genotype data generated in Step S4, use the optimized model to predict the target trait value: where is the genotype data matrix.

[0043] Furthermore, calculate the prediction accuracy, root mean square error, and coefficient of determination of the model;

[0044] Prediction accuracy , is the actual phenotypic value of the th sample, is the predicted target trait value of the

[0045] th sample, and is the number of samples. The closer the value of the prediction accuracy is to 1, the more accurate the prediction;

[0046] Root mean square error , and the smaller the value of the root mean square error, the more accurate the prediction; is the average value of the actual phenotypic values, and the closer the coefficient of determination value is to 1, the stronger the model's explanatory ability.

[0047] Furthermore, the GBDT model in step S52 further adopts feature screening and feature importance analysis, and uses the recursive feature elimination (RFE) method to screen important genomic marker sites. The recursive process is as follows:

[0048] At initialization, use all features, not less than 500;

[0049] Train the model and calculate the importance score of each feature ;

[0050] Remove the feature with the lowest importance score: ;

[0051] If the number of features is greater than the preset threshold, and the preset threshold is between 10% and 30% of the total number of features, continue to remove the feature with the lowest importance score; otherwise, obtain the feature subset ;

[0052] Calculate the importance index of each retained feature : , where is the importance score of feature in the th tree, is the total number of trees in the model;

[0053] Construct a spatial distribution map of genotype features, and visualize the distribution of marker sites with an importance index greater than 0.05 on the chromosome to locate potential functional regions.

[0054] Furthermore, by reducing the population size and the number of generations in the pure line selection process, resource optimization allocation is achieved; specifically: according to the predicted phenotypic values, sort the phenotypic values of the self-crossed offspring of each generation, and only select the top 20% of the individuals with the highest phenotypic values for the next generation of reproduction.

[0055] Compared with the prior art, the beneficial effects of the present invention are:

[0056] The present invention improves data quality and analysis efficiency by constructing a parental genotype dataset for inbred crop hybrid combinations and screening for effective marker loci; constructs a genomic linkage map based on the linkage disequilibrium coefficient, accurately calculates the recombination rate between loci, and through random simulation of the inbreeding process, and identifies key gene loci through feature screening and importance analysis, significantly improving the prediction accuracy; the method of the present invention provides an efficient and accurate calculation tool for inbred crop hybrid breeding, which can significantly shorten the breeding cycle and improve the breeding precision. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 1 It is a flowchart of a method for predicting and simulating the genotypes of advanced generations of inbred crop hybrid combinations of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0058] To make the objectives, technical solutions and advantages of the present invention clearer, the technical solutions in the present invention will be clearly and completely described below. Apparently, the described embodiments are some but not all of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art without making creative efforts based on the embodiments in the present invention belong to the scope of protection of the present invention.

[0059] As Figure 1 shown, it is a method for predicting and simulating the genotypes of advanced generations of inbred crop hybrid combinations of the present invention, and the method includes the following steps:

[0060] Step S1, construct a parental genotype dataset for inbred crop hybrid combinations, screen genomic marker loci for the parental genotype dataset, and obtain effective marker loci.

[0061] Parental genotype data can be obtained through technologies such as gene chips and high-throughput sequencing. For different inbred crops, such as rice, wheat, corn, etc., their genome sizes and structures are different, so data collection and collation need to be carried out according to the characteristics of different crops. Genomic marker locus screening is to improve the accuracy and calculation efficiency of subsequent analysis and remove marker loci with poor quality and low information content. The following details the steps of marker locus screening:

[0062] Step S11, calculate the polymorphism information content value of each marker locus:

[0063] , where is the frequency of the th allele, and is the number of alleles;

[0064] Step S12, calculate the missing rate and minor allele frequency of the marker locus, and remove those with a missing rate greater than or the minor allele frequency is less than of the marker locus;

[0065] Step S13, sort the marker loci according to the genomic physical positions, and when the distance between adjacent marker loci is less than a preset threshold, retain the marker locus with a higher PIC value;

[0066] The preset threshold depends on the self-pollinated crop species; for different self-pollinated crops, the setting of the preset threshold is different. For example, for crops with a smaller genome such as rice, the threshold for the distance between adjacent marker loci is set to 5 - 10 kb; for crops with a larger genome such as wheat, the threshold is set to 20 - 50 kb. This screening based on physical distance can effectively avoid redundant information caused by overly dense distribution of marker loci, while retaining sufficient genomic coverage to ensure the accuracy of subsequent analysis. The number of marker loci after screening can usually be reduced from hundreds of thousands initially to tens of thousands, greatly improving the computational efficiency while maintaining the integrity of genomic information.

[0067] Step S2, based on the parental genotype data, calculate the linkage disequilibrium coefficient between the marker loci, construct a genomic linkage map, and determine the recombination rate between the loci.

[0068] It should be noted that linkage disequilibrium (LD) refers to the phenomenon of non-random combination of alleles at different loci on the genome, which is an important concept in genetics and breeding research. By calculating the LD coefficient between loci, the linkage status of different regions on the genome can be understood, providing a basis for subsequent recombination rate calculation. The construction of a linkage map helps to identify linkage blocks (LD blocks) on chromosomes, and these regions often inherit as a whole during meiosis, which is of great significance for predicting the genotype segregation of hybrid offspring.

[0069] The calculation of the linkage disequilibrium coefficient between the marker loci in step S2 includes:

[0070] Step S21, calculate the linkage disequilibrium coefficient between and : , where is the frequency of a specific allele combination at the two loci, and are the frequencies of the respective alleles at loci and respectively;

[0071] Step S22, calculate the standardized linkage disequilibrium coefficient : ;

[0072] Step S23, based on values to construct a linkage map. When it is considered that two loci are strongly linked. When it is considered moderately linked. When it is considered weakly linked or unlinked.

[0073] The calculation formula for determining the recombination rate between loci in step S2 is:

[0074] , where is the recombination rate between locus and , is the physical distance between locus and , is the average recombination length of the inbred crop chromosome, which is adaptively adjusted according to the crop variety, and the value range is within cM.

[0075] The calculation of the recombination rate takes into account the non-linear relationship between physical distance and genetic distance, which conforms to the chromosome exchange law in the actual biological process. The setting of parameter needs to be adjusted according to the characteristics of different crops. For example, the value of rice is usually in the range of 20 - 30 cM, that of wheat is in the range of 40 - 60 cM, and that of maize is in the range of 30 - 50 cM.

[0076] A larger value indicates a lower chromosome exchange frequency and is applicable to crops with a lower recombination rate; a smaller value indicates a higher chromosome exchange frequency and is applicable to crops with active recombination.

[0077] Step S3, use the random simulation algorithm to construct a gamete formation model, calculate the chromosome exchange probability according to the recombination rate, and generate the genotypes of the

[0078] The random simulation algorithm for constructing a gamete formation model in step S3 includes:

[0079] Step S31, based on the parental genotypes, randomly generate homologous chromosome exchange sites, and the probability of the exchange sites appearing is proportional to the recombination rate between the loci;

[0080] Step S32, simulate the gamete genotypes after chromosome exchange. For two parents and , for locus , if locus If it is located within the odd crossover interval, the genotype of the hybrid individual at this locus is from the parent , expressed as: , if the locus is located within the even crossover interval, the genotype of the hybrid individual at this locus is from the parent , expressed as: , where and are the genotypes of the parents and the parent at the locus ;

[0081] Step S33, randomly select simulated parental gamete combinations to form generation hybrid individuals, The genotype of the individual is calculated as: , where and are the genotypes of the male and female gametes at locus k respectively.

[0082] The genotype of the F1 generation individuals is calculated by weighted average, reflecting the characteristics of heterozygous loci. In practical applications, genotypes are usually represented by 0, 1, 2, where 0 represents a homozygote (AA), 2 represents another homozygote (BB), and 1 represents a heterozygote (AB).

[0083] Therefore, the genotype value of the F1 generation hybrids at most loci is 1, indicating that they inherit one allele from each parent. This coding method facilitates subsequent numerical calculations and the application of machine learning models; it should be noted that the gamete formation model takes into account the randomness of chromosome crossover, making the simulation results closer to the genetic laws in the actual breeding process.

[0084] Step S4, based on the genotypes of the generation individuals, generate genotype data of the to generation individuals through random simulation of the self-crossing process.

[0085] The specific steps for generating high-generation individuals through random simulation of the self-crossing process in step S4 are as follows:

[0086] Step S41, self-cross the individuals, simulate the chromosome crossover during meiosis, and produce gametes;

[0087] Step S42, randomly combine the produced gametes to form a population of individuals, and the genotype of the individual at the locus is calculated as: , where and are the genotypes of two gametes at locus ;

[0088] Step S43, continue to self - cross the individuals, and similarly generate generation genotypes, and so on until generation;

[0089] Step S44, calculate the homozygosity of each generation: , where is the total number of marker loci, is an indicator function, which takes the value of 1 when the locus is homozygous, and 0 otherwise; and represent two states when the locus is homozygous.

[0090] As the number of self - cross generations increases, the homozygosity of individuals gradually increases. In theory, during the self - crossing process of self - pollinated crops, the homozygosity H should increase according to the formula , where n is the number of self - cross generations.

[0091] For example, the theoretical homozygosity of the F2 generation is 0.5, that of the F5 generation is 0.9375, and that of the F8 generation is 0.9922. The homozygosity obtained by simulating and calculating this method usually has good consistency with the theoretical value. However, considering the selection pressure and genotype adaptability differences in the actual breeding process, the real situation may deviate. By calculating the homozygosity of each generation, breeders can evaluate the degree of self - crossing.

[0092] Step S5, use a machine - learning model to establish a mapping relationship between genotypes and phenotypes, and predict the target trait performance of high - generation self - crossed lines.

[0093] The establishment of the mapping relationship between genotypes and phenotypes using a machine - learning model in the said step S5 includes:

[0094] Step S51, construct a training data set, including sample data with known genotypes and phenotypes, and divide the data into a training set and a validation set according to a ratio of 7:3;

[0095] The training data can be from the phenotypic data accumulated in previous breeding experiments or from genotype - phenotype association data in public databases; the data is divided in a ratio of 7:3 instead of the traditional 8:2 to ensure sufficient model training while providing enough validation samples to evaluate the generalization ability of the model.

[0096] For the case of a small sample size, the ratio can be appropriately adjusted, such as 6:4 or 5:5. In addition, during the data partitioning process, attention should be paid to maintaining the distribution consistency between the training set and the validation set. The stratified sampling method can be used to ensure that the two datasets have similar statistical characteristics.

[0097] Step S52: Use the Gradient Boosting Decision Tree model GBDT as the core prediction model. This model consists of multiple decision trees. The calculation process of the model parameters is as follows:

[0098] Initialize the model , where is the loss function; represents finding the value that minimizes the following expression, is the predicted value, is the th true phenotypic value of the sample;

[0099] For , is the number of decision trees. The residual of the th sample in the th round of iteration is:

[0100] , is the partial derivative of the loss function with respect to the current predicted value, is the model after the m - 1 th round of iteration, is the genotype data of the th sample; Fit the residual to generate a decision tree and obtain the th leaf node region

[0101] in the th round of iteration; Calculate the optimal output value of each leaf node: ; Update the model:

[0102] The GBDT model is suitable for genotype-phenotype prediction because it can capture the complex non-linear interaction between gene loci, which is consistent with the epistasis phenomenon in biology. Compared with traditional linear models such as multiple linear regression, GBDT has obvious advantages in dealing with high-dimensional features and complex relationships.

[0103] In this method, the loss function L can be selected according to different prediction targets. For continuous traits such as yield and plant height, the mean squared error (MSE) is usually selected as the loss function; for discrete traits such as disease resistance level, the logarithmic likelihood loss function can be selected. The learning rate η controls the contribution of each tree to the final model. A smaller learning rate helps to improve the model accuracy but requires more iteration times; a larger learning rate can accelerate the convergence speed but may lead to overfitting.

[0104] Step S53, the model optimization adopts the 5-fold cross-validation method, and determines the optimal hyperparameters through grid search, including the learning rate , the maximum depth of the tree, the number of decision trees and the minimum number of samples in a leaf node.

[0105] That is, the dataset is divided into 5 parts. Each time, 4 parts are used as the training set and 1 part is used as the validation set, and this is repeated 5 times so that each part of the data has the opportunity to be used as the validation set; during the grid search process, the value ranges of the hyperparameters include: the learning rate , the maximum depth of the tree , the number of decision trees , and the minimum number of samples in a leaf node . By combining different hyperparameter values for model training and evaluation, the parameter combination with the best performance on the validation set is selected as the final model parameters.

[0106] Step S54, based on the high-generation genotype data generated in Step S4, use the optimized model to predict the target trait value: , where is the genotype data matrix.

[0107] Calculate the prediction accuracy, root mean square error, and coefficient of determination of the model; the selection of model evaluation metrics should be based on the characteristics of the prediction task and the breeding objectives. The prediction accuracy (Accuracy) intuitively reflects the relative error between the predicted value and the actual value and is suitable for evaluating the prediction accuracy under different trait scales; the root mean square error (RMSE) provides an absolute measure of the prediction error and is suitable for comparing models of the same trait; the coefficient of determination (R²) measures the degree to which the model explains the variation of the dependent variable and is suitable for evaluating the overall goodness of fit of the model.

[0108] These three metrics should be considered comprehensively. For example, a model with a high R² but a high RMSE may capture the overall trend in the data but have a systematic bias; a model with a low RMSE but a low R² may predict accurately on average but cannot explain the variability of the data. For different breeding objectives, different evaluation metrics may need to be emphasized. For example, yield breeding may pay more attention to prediction accuracy, while quality breeding may pay more attention to the coefficient of determination.

[0109] Prediction accuracy , is the actual phenotypic value of the i-th sample, is the predicted target trait value of the i-th sample, is the number of samples. The closer the value of the prediction accuracy is to 1, the more accurate the prediction;

[0110] Root mean square error ; The smaller the root mean square error value, the more accurate the prediction;

[0111] Coefficient of determination , is the average value of the actual phenotypic values. The closer the value of the coefficient of determination is to 1, the stronger the model's explanatory ability.

[0112] The GBDT model in the step S52 further adopts feature screening and feature importance analysis, and uses the recursive feature elimination (RFE) method to screen important genomic marker sites. The recursive process is as follows:

[0113] At initialization, use all features, not less than 500;

[0114] Train the model and calculate the importance score of each feature ;

[0115] Remove the feature with the lowest importance score: ;

[0116] If the number of features is greater than the preset threshold, and the preset threshold is taken between 10% and 30% of the total number of features, continue to remove the feature with the lowest importance score; otherwise, obtain the feature subset

[0117] Calculate the importance index of each remaining feature : , where is the importance score of the feature in the th tree, is the total number of trees in the model;

[0118] Construct a spatial distribution map of the genotype features, and visualize the distribution of the marker sites with the importance index greater than 0.05 on the chromosome to locate potential functional regions.

[0119] Marker loci with an importance index greater than 0.05 are usually significantly associated with the target trait and may be located at or near the functional gene region controlling the trait. By visualizing the distribution of these important loci on the chromosome, potential QTL (quantitative trait locus) regions can be identified, providing clues for marker-assisted selection and gene function research. In practical applications, researchers can compare these loci with known gene annotation information to further verify the biological significance of the prediction results. For example, the analysis of rice yield traits may find that some important loci are concentrated on chromosomes 1, 3, and 6, which is consistent with the yield-related QTL regions found in previous studies.

[0120] Optimize resource allocation by reducing the population size and the number of generations in the pure-line selection process; specifically: according to the predicted phenotypic values, rank the phenotypic values of the selfed offspring in each generation and only select the top 20% of the individuals with the highest phenotypic values for reproduction in the next generation.

[0121] An overly high selection intensity (such as selecting the top 5%) can obtain more extreme phenotypic values, but it will greatly reduce the population genetic diversity and increase the risk of inbreeding depression; an overly low selection intensity (such as selecting the top 50%) cannot effectively improve the target trait.

[0122] A selection intensity of 20% can achieve a better balance between genetic gain and diversity maintenance; breeders can adjust this ratio appropriately according to specific crops and breeding goals. For example, for high-heritability traits such as stress resistance, a higher selection intensity can be adopted; for low-heritability traits such as yield, the selection intensity may need to be reduced to maintain sufficient genetic variation. In addition, this method can also achieve comprehensive selection of multiple traits. By constructing a comprehensive selection index, the predicted values of multiple target traits are weighted and integrated to achieve multi-objective breeding optimization.

[0123] To further illustrate the method and technical effects of the present invention and more clearly illustrate the implementation process of the present invention, the following takes the rice hybrid combination "Japonica A × Indica B" as an example to introduce in detail the specific application of the high-generation genotype prediction simulation method for self-crossing crop hybrid combinations.

[0124] Obtain the genomic data of the parents "Japonica A" and "Indica B" through high-throughput sequencing technology. The initial data contains approximately 250,000 SNP marker loci across the genome. Screen these marker loci:

[0125] Calculate the polymorphism information content (PIC) value of each marker locus. Taking the marker Chr01_10254562 (located on chromosome 1, with a physical position of 10,254,562 bp) as an example, the allele frequencies of this locus in the breeding population are P1 = 0.65 (A allele) and P2 = 0.35 (G allele) respectively. According to the PIC calculation formula: PIC = 0.455.

[0126] At the same time, calculate the deletion rate of the marker locus. It is found that the deletion rate of the Chr02_24513687 locus is 12.5%, which is greater than the preset threshold of 10%, so it is excluded; the minor allele frequency of the Chr03_5689421 locus is 3.2%, which is less than the preset threshold of 5%, and it is also excluded.

[0127] Further screen the marker loci with similar physical positions. In the rice genome, set the adjacent marker locus distance threshold to 5 kb. For example, the distance between the two loci Chr05_15687452 and Chr05_15689876 is 2,424 bp, which is less than the threshold. Their PIC values are 0.423 and 0.468 respectively. Retain the Chr05_15689876 locus with a higher PIC value.

[0128] After the above screening steps, a total of 42,568 high-quality marker loci are finally retained for subsequent analysis.

[0129] Based on the screened marker loci, calculate the linkage disequilibrium coefficient between loci. Taking the two loci Chr06_8945632 and Chr06_9012345 as an example:

[0130] In the population, the A allele frequency of the locus Chr06_8945632 , the T allele frequency of the locus Chr06_9012345 , the frequency of the AT combination at the two loci .

[0131] Linkage disequilibrium coefficient , the standardized linkage disequilibrium coefficient r² = 0.628; according to the r² value judgment, these two loci belong to moderate linkage (0.5 < r² ≤ 0.8). For the rice genome, set the average recombination length parameter l = 25 cM. For the above two loci, their physical distance = 66,713 bp. According to the recombination rate calculation formula: = 0.0825. This indicates that the recombination rate between these two loci is 8.25%, that is, there is about an 8.25% probability of exchange during meiosis.

[0132] By similar calculations, the recombination rates between all adjacent marker loci are obtained to construct a complete genomic linkage map.

[0133] A gamete formation model is constructed using a random simulation algorithm. Taking a region on chromosome 8 as an example, it contains 10 consecutive marker loci (simplified as 1 - 10).

[0134] Assume that the genotype of the parent "Japonica rice A" at these 10 loci is [0,0,0,0,0,0,0,0,0,0] (indicating all are AA homozygotes), and the genotype of the parent "Indica rice B" is [2,2,2,2,2,2,2,2,2,2] (indicating all are BB homozygotes).

[0135] According to the calculated recombination rate, chromosome exchange events are simulated; for example, during the formation of a gamete, randomly determine that the exchange site appears after locus 3 and locus 7:

[0136] If the locus is in an odd exchange interval: the genotype comes from the parent "Japonica rice A" and is 0; if the locus is in an even exchange interval: the genotype comes from the parent "Indica rice B" and is 2; then the simulated gamete genotype is: [0,0,0,2,2,2,2,0,0,0].

[0137] By simulating different combinations of exchange sites multiple times, multiple gametes are generated; randomly select two gametes to combine to form an F1 individual:

[0138] Gamete 1: [0,0,0,2,2,2,2,0,0,0];

[0139] Gamete 2: [0,0,2,2,2,2,0,0,0,0];

[0140] F1 individual genotype = (Gamete 1 + Gamete 2) / 2 = [0,0,1,2,2,2,1,0,0,0];

[0141] Among them, the locus with a value of 1 represents a heterozygous genotype (AB), the value of 0 represents a homozygous genotype (AA), and the value of 2 represents a homozygous genotype (BB).

[0142] In this way, 100 F1 individuals containing 42,568 marker loci across the whole genome are simulated.

[0143] Taking a specific F1 individual as an example, continue to simulate its self - crossing process. Assume that the genotype of this F1 individual on a certain chromosomal segment is [1,1,1,1,1,1,1,1,1,1] (all are heterozygous). During the meiosis of F1 to form gametes, chromosome exchange events are also simulated. For example, two gametes are generated:

[0144] Gamete 1: [0,0,0,2,2,2,2,0,0,0];

[0145] Gamete 2: [0,0,2,2,2,0,0,0,0,0];

[0146] These two gametes randomly combine to form F2 individuals, and their genotypes are:

[0147] F2 individual = (Gamete 1 + Gamete 2) / 2 = [0,0,1,2,2,1,1,0,0,0];

[0148] Calculate the homozygosity H of this F2 individual = 0.7; indicating that 70% of the loci have reached the homozygous state, which is slightly higher than the theoretical expected 50% homozygosity of the F2 generation. This difference is normal random fluctuation in actual breeding.

[0149] Continue to self - cross - simulate the F2 individuals to generate F3 - generation individuals. By the same method, simulate F4, F5... until F8 - generation individuals in turn. As the number of self - cross generations increases, the homozygosity gradually increases:

[0150] Average homozygosity of F2 generation: 0.53 (theoretical value 0.50);

[0151] Average homozygosity of F3 generation: 0.76 (theoretical value 0.75);

[0152] Average homozygosity of F4 generation: 0.88 (theoretical value 0.875);

[0153] Average homozygosity of F5 generation: 0.94 (theoretical value 0.9375);

[0154] Average homozygosity of F6 generation: 0.97 (theoretical value 0.9688);

[0155] Average homozygosity of F7 generation: 0.99 (theoretical value 0.9844);

[0156] Average homozygosity of F8 generation: 0.99 (theoretical value 0.9922).

[0157] Collected the genotype data and yield data of 500 known rice varieties as the training set, and divided them into 350 training samples and 150 validation samples according to the ratio of 7:3. Used the GBDT model for training, and initialized the model parameters:

[0158] Learning rate η = 0.05; maximum depth of decision tree = 5; number of decision trees M = 300; minimum number of samples in leaf nodes = 10; mean squared error (MSE) is used as the loss function, and hyperparameters are optimized through 5-fold cross-validation. The finally determined optimal parameters are: learning rate η = 0.03; maximum depth of decision tree = 7; number of decision trees M = 500; minimum number of samples in leaf nodes = 8.

[0159] The performance of the model was evaluated using the validation set: prediction accuracy Accuracy = 0.92; root mean squared error RMSE = 325.6 kg / ha; coefficient of determination R² = 0.87. This result indicates that the model has high prediction accuracy and can explain 87% of the yield variation.

[0160] The recursive feature elimination (RFE) method was used to screen important genomic marker loci. Initially, all 42,568 marker loci were used. Through iterative training and feature importance evaluation, finally 1,278 important marker loci were retained (about 3% of the total).

[0161] The importance index of each retained feature was calculated, and it was found that the importance indices of 38 marker loci were greater than 0.05, mainly distributed on chromosomes 1, 3, 6, 8, and 11. Visualizing these loci on the chromosomes, it was found that multiple loci clustered near the known QTL regions related to rice yield, such as the Gn1a (controlling grain number) on chromosome 6 and the sd1 (controlling plant height, indirectly affecting yield) gene region on chromosome 1.

[0162] Based on the prediction model, the yields of 1,000 individuals in the F2 generation were predicted, and they were ranked according to the predicted values. Only the top 20% (i.e., 200) individuals were selected to continue cultivation to the F3 generation.

[0163] Among the 200 individuals in the F3 generation, yield prediction and selection were carried out again, and the top 20% (i.e., 40) individuals were retained to continue cultivation to the F4 generation, and so on.

[0164] Finally, 5 high-yield pure lines were obtained in the F8 generation. The actual measured yield of the best strain "Rice High Yield No. 8" was 9,850 kg / ha, which was 12.3% higher than that of the control variety selected by the conventional breeding method.

[0165] Through the application of this method, the breeding cycle in rice breeding is shortened. The conventional breeding method takes 8 - 10 years to obtain stable high-yield varieties, while after applying this method, the breeding cycle is shortened to 5 years, improving the breeding efficiency. The resource utilization efficiency is enhanced. The traditional method requires cultivating and evaluating thousands of individuals from the F2 to F8 generations, while this method, through early prediction and screening, cultivates 1,000 individuals in the F2 generation, 200 in the F3 generation, 40 in the F4 generation, 8 in the F5 generation, and 5 in each of the F6 - F8 generations, totaling 1,258 individuals, saving about 75% of the field resources compared to the traditional method.

[0166] Among the finally selected 5 pure lines, 4 exhibit excellent trait combinations, with a success rate of 80%, far higher than the usually less than 30% success rate in traditional breeding methods. Through feature importance analysis, 6 new gene loci potentially related to rice yield were discovered, providing clues for further gene function research.

[0167] This example fully demonstrates the effectiveness and superiority of this method in actual breeding applications. It can not only accurately predict the phenotypic traits of high-generation inbred lines, but also significantly improve breeding efficiency, save breeding resources, and accelerate the cultivation process of excellent varieties.

[0168] It should be noted that this method is not only applicable to rice, but can also be extended to other important self-pollinated crops such as wheat, soybean, cotton, etc.; for different self-pollinated crops, relevant parameters need to be adjusted according to their genomic characteristics: wheat, as an allohexaploid with a high genomic complexity, should increase the strictness of marker site screening, adjust the adjacent marker site spacing threshold to 30 - 50 kb, and at the same time, the average chromosome recombination length parameter l should be set to 40 - 60 cM to reflect its lower recombination frequency; there are more repetitive sequences in the soybean genome, and when screening marker sites, the minor allele frequency threshold should be increased to 8% to reduce the influence of false positive sites, and its recombination length parameter l is usually set to 30 - 45 cM; for cotton, the differences between the A and D sub-genomes need to be considered, and the linkage disequilibrium parameters are set separately. The threshold can be appropriately reduced.

[0169] In addition, for the training of machine learning models for different crops, the GBDT model parameters also need to be adjusted according to the heritability of specific traits. For example, yield traits are more affected by the environment in wheat, and the learning rate should be reduced to 0.01 - 0.02, and the number of trees should be increased to 800 - 1000 to improve the model stability; it is estimated that this method can shorten the breeding cycle by about 40% in wheat breeding and increase the selection accuracy by 65% in soybean breeding.

[0170] The specific embodiments described above further elaborate on the purpose, technical solutions, and beneficial effects of the present invention. It should be understood that the above description is only for the specific embodiments of the present invention and is not used to limit the protection scope of the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention shall be included within the protection scope of the present invention.

Claims

1. A method for predicting and simulating the high-generation genotypes of self-fertilized crop hybrid combinations, characterized in that: The method comprises the following steps: Step S1, constructing a parental genotype dataset of a self-pollinated crop hybrid combination, and performing genomic marker site screening on the parental genotype dataset to obtain effective marker sites; The specific steps for screening genomic marker sites for parental genotype datasets are: Step S11, calculating the polymorphism information content of each marker site value: ,in For the The frequency of the alleles, is the number of alleles; Step S12, calculate the deletion rate and minor allele frequency of the marker locus, and eliminate those with a deletion rate greater than or the minor allele frequency is less than The marker site of Step S13, sorting the marker sites according to the physical positions of the genome, and retaining the marker sites with higher PIC values ​​when the distance between adjacent marker sites is less than a preset threshold; The preset threshold value depends on the type of self-pollinating crop, which is 5-10 kb for rice and 20-50 kb for wheat; Step S2, based on the parental genotype data, calculate the linkage disequilibrium coefficient between marker sites, construct a genome linkage map, and determine the recombination rate between sites; Calculation of linkage disequilibrium coefficients between marker loci includes: Step S21, calculate the location and Linkage disequilibrium coefficient : ,in is the frequency of a specific allele combination at two loci, and Site and The frequency of each allele; Step S22, calculating the standardized linkage disequilibrium coefficient : ; Step S23, based on The linkage map is constructed when When the two points are considered to be strongly linked, moderate linkage, When it is considered weakly linked or unlinked; The formula for determining the recombination rate between sites is: ,in For site and The recombination rate between For site and The physical distance between is the average recombinant length of the chromosome of the self-pollinated crop, which is adjusted adaptively according to the crop type and has a value range of cM; Step S3, using a random simulation algorithm to construct a gamete formation model, calculate the probability of chromosome exchange based on the recombination rate, and generate The genotype of the individuals of the generation; The stochastic simulation algorithm constructs the gamete formation model including: Step S31, based on the parental genotypes, randomly generate homologous chromosome exchange sites, and the probability of the exchange site appearing is proportional to the recombination rate between the sites; Step S32, simulating the gamete genotype after chromosome exchange, for the two parents and , for the site , if the point If the hybrid individual is located in an odd-numbered exchange interval, the genotype at this site comes from the parent , expressed as: , if the point If the hybrid individual is located in an even-numbered exchange interval, the genotype at this site comes from the parent , expressed as: ,in and Parents and parent At the location genotype; Step S33, randomly select simulated parent gametes to form The hybrid individuals The genotype of an individual is calculated as: ,in and are the genotypes of the paternal and maternal gametes at site k, respectively; Step S4, based on The genotypes of the individuals of the next generation are generated by random simulation of the selfing process. to Genotype data of the individuals of the previous generation; The specific steps of generating high-generation individuals through random simulation of the self-crossing process are: Step S41: Individuals self-fertilize, simulating The exchange of chromosomes during meiosis produces gametes; Step S42, random combination The gametes produced form Individual groups, Individual at the site The genotype is calculated as: ,in and For two gametes at the locus genotype; Step S43: Individuals continue to self-fertilize, similarly generating The genotypes are generated in this way. generation; Step S44, calculate the homozygosity of each generation: ,in is the total number of marker sites, is the indicative function, when the site The value is 1 when it is homozygous, otherwise it is 0; and Indicates the site Two states when homozygous; Step S5, using a machine learning model to establish a mapping relationship between genotype and phenotype, and predicting the target trait performance of high-generation inbred lines; The mapping relationship between genotype and phenotype established by machine learning model includes: Step S51, constructing a training data set, including sample data with known genotypes and phenotypes, and dividing the data into a training set and a validation set in a ratio of 7:3; Step S52, using the gradient boosting decision tree model GBDT as the core prediction model, the model consists of multiple decision trees, and the model parameter calculation process is: Initialize the model ,in, is the loss function; It means to find the expression that minimizes the following expression. value, is the predicted value, For the The true phenotypic value of each sample; for , is the number of decision trees, In the iteration The residual of the sample : , is the partial derivative of the loss function with respect to the current predicted value, is the model after the m-1th round of iteration, For the Genotype data of samples; fitting residuals Generate a decision tree and get In the iteration Leaf node area ; Calculate the optimal output value for each leaf node: ; Update the model: ,in is the learning rate, with a value range of ; Step S53, model optimization using The fold cross validation method, , determine the optimal hyperparameters through grid search, including the learning rate , maximum tree depth, number of decision trees and the minimum number of leaf node samples; Step S54, based on the high-generation genotype data generated in step S4, use the optimized model to predict the target trait value: ,in is the genotype data matrix.

2. The method for predicting and simulating high-generation genotypes of self-pollinated crop hybrid combinations according to claim 1, characterized in that: Calculate the prediction accuracy, root mean square error, and coefficient of determination of the model; Prediction Accuracy , is the actual phenotypic value of the i-th sample, is the predicted target trait value of the i-th sample, is the number of samples, and the closer the prediction accuracy value is to 1, the more accurate the prediction is; Root mean square error ; The smaller the root mean square error value, the more accurate the prediction; Coefficient of determination , It is the average value of the actual phenotypic value. The closer the coefficient of determination value is to 1, the stronger the explanatory power of the model is.

3. The method for predicting and simulating high-generation genotypes of self-pollinated crop hybrid combinations according to claim 2, characterized in that: The GBDT model in step S52 further adopts feature screening and feature importance analysis, and uses the recursive feature elimination (RFE) method to screen important genomic marker sites. The recursive process is: When initializing, use all Features, Not less than 500; Train the model and calculate the importance score for each feature ; Remove the feature with the lowest importance score: ; If the number of features If it is greater than the preset threshold, the features with the lowest importance score will be removed; Otherwise, get the feature subset ; The preset threshold is between 10% and 30% of the total number of features Calculate the importance index of each retained feature : ,in Features In the The importance scores in the tree, is the total number of trees in the model; Construct a spatial distribution map of genotype characteristics and convert the importance index The distribution of marker loci with a p-value greater than 0.05 on chromosomes was visualized to locate potential functional regions.

4. The method for predicting and simulating high-generation genotypes of self-pollinated crop hybrid combinations according to claim 3, characterized in that: By reducing the population size and the number of generations in the pure line selection process, optimal resource allocation can be achieved; specifically, the phenotypic values ​​of the self-pollinated offspring of each generation are ranked according to the predicted phenotypic values, and only the individuals with the top 20% of the phenotypic values ​​are selected for the next generation of reproduction.

Images