Performance prediction of crosses in plant breeding through genotype by environment interactions

Performance prediction of potential crosses plays a significant role in plant breeding, which aims to produce new crop varieties that have higher yields, require fewer resources, and are more adaptable to the changing environments. In the 2020 Syngenta crop challenge, Syngenta challenged participants to predict the yield performance of a list of potential breeding crosses of inbreds and testers based on their historical yield data in different environments. They released a dataset that contained the observed yields for 294,128 corn hybrids through the crossing of 593 unique inbreds and 496 unique testers across multiple environments between 2016 and 2018. To address this challenge, we designed a new predictive approach that integrates random forest and an optimization model for G ×\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times $$\end{document} E interaction detection. Our computational experiment found that our approach achieved a relative root-mean-square-error (RMSE) of 0.0869 for the validation data, outperforming other state-of-the-art models such as factorization machine and extreme gradient boosting tree. Our model was also able to detect genotype by environment interactions that are potentially biologically insightful. This model won the first place in the 2020 Syngenta crop challenge in analytics.

www.nature.com/scientificreports/ effects and environment for genomic predictions. State-of-the-art machine learning models have also been used for crop yield prediction, including stepwise multiple linear regression 25 , neural networks 25,26 , convolutional neural networks 27,28 , recurrent neural networks 28 , multiple regression 26 , random forest 29 , weighted histograms regression 30 , and association rule mining and decision tree 31 .
In this paper, we propose a new model for predicting the yield performance of new hybrids based on historical data of other hybrids. This model integrates a random forest with a combinatorial optimization-based interaction-detection model and attempts to combine their strengths. The random forest model 32 is known for its capability to approximate general form nonlinear relationships among the variables. On the other hand, the interaction-detection model originated from a recently published algorithm 33 that has been shown to be particularly effective in detecting epistatic type of interactions. Our model extends that algorithm to the detection of genotype by environment interactions (G × E).
Our computational results using the 2020 Syngenta crop challenge data suggested that the proposed model can accurately predict the performance of untested cross combinations of inbreds and testers. Moreover, results of our prediction model can also reveal biologically meaningful insights, such as the best hybrids for specific environments.

Problem definition
Most of the effort in a breeding program is related to evaluating inbreds by crossing to another inbred known as a tester. According to the problem statement of the 2020 Syngenta crop challenge, "it is a plant breeder's job to identify the best parent combinations by creating experimental hybrids and assessing the hybrids' performance by 'testing' it in multiple environments to identify the hybrids that perform best. " While the yield performance of a hybrid is largely related to the parents, it is also affected by many factors that are hard to predict, such as heterosis and interactions between genotype and the environment.
The objective of the 2020 Syngenta crop challenge was to design a model for predicting the yield performance of a list of inbred-tester combinations based on historical datasets that included yield, genetic group, and pedigree information of hybrids collected in different environments over a number of years. If successful, this challenge will stimulate novel design of predictive models and algorithms for yield prediction of inbred-tester combinations and progeny testing of inbreds, which will help breeders make the most promising crosses without having to rely on large-scale trial-and-error that is expensive, labor intensive, and time consuming. The 2020 Syngenta crop challenge released the following dataset for commercial corn. training dataset.
• Yield: Historical yield performances were measured for 10,919 unique biparental hybrids. To provide realistic data without revealing proprietary information, actual yield values were anonymized to make the average and standard deviation of yields approximately 1.0 and 0.1, respectively. The range of the yields was from 0.0472 to 1.8001. • Genetic clusters: No genetic marker information was available, but the genetic clusters of 593 unique inbreds and 496 unique testers were provided. Syngenta grouped the inbreds and testers into some clusters according to their genetic similarities using internal methods. There were 14 inbred clusters and 13 tester clusters. • Environment: Out of a total of 593 × 496 = 294,128 possible combinations of inbred-tester crosses, the training data included 10,919 unique hybrids that were planted across 280 locations between 2016 and 2018, each year with a unique set of weather conditions. The information that we had for the environment is 280 location IDs and 3 years such that there were 599 unique location-weather combinations in the training set. The total number of unique hybrids-location-weather combinations was 155,765, some of which had multiple replications, so the total number of yield records was 199,476. However, this training dataset represents only 0.08% of all possible 593 × 496 × 280 × 3 = 247,067,520 hybrids-location-weather combinations.
test dataset. The test dataset includes a set of inbred-tester combinations whose yield performances need to be predicted. The environments in which these hybrids would be grown were not specified in the crop challenge. evaluation criteria. The evaluation criteria for the 2020 Syngenta crop challenge in analytics were "accuracy of the predicted values in the test set based on root mean squared error, simplicity and intuitiveness of the solution, clarity in the explanation, and the quality and clarity of the finalist's presentation at the 2020 INFORMS Conference on Business Analytics and Operations Research. " Our model won the first place in this competition. For this paper, we evaluated the proposed model in terms of prediction accuracy. Because we did not have access to the ground truth yield of the test dataset, we divided the given dataset to training and validation subsets using tenfold cross-validation (CV). Then, we used the average performance of the proposed model as the evaluation criteria.

Method
Data preprocessing. We defined the input variable X as one-hot coding of hybrid-location-weather combinations and the output variable y as the corresponding yield. To accommodate this definition, four types of training data were converted to binary using the one-hot coding preprocessing: inbred and tester indices, genetic cluster, location ID, and weather. For those hybrid-location-weather combinations with multiple replications, the average yield was used as the output data. As such, the training data has a dimension of 155,765 observations proposed model and algorithm. We proposed a hybrid model for this challenge, which combines random forest with G × E interaction detection techniques. The overview of the model is diagramed in Fig. 1. This model consists of three main components: a random forest model that captures the complex nonlinear relationship between input and output variables, a G × E interaction detection model that captures interactions among hybrid, location, and weather variables, and another random forest model that utilizes the interactions to augment the prediction performance of the first random forest model. Details of these components are described in the rest of this section.
Random forest model 1. Random forest 32 is an ensemble learning model that can be used for classification or regression by constructing a multitude of decision trees. To grow each tree, a random subset of features is selected along with replacement sampling (bootstrap sampling) used to select different subsets of the observations. Therefore, observations in the dataset that were not included in the bootstrapped samples are considered as out-of-bag observations, and the performance of the tree is evaluated by the average out-of-bag error. Due to the builtin component of cross-validation, the random forest is less prone to overfitting. The random forest model 1 takes the one-hot matrix X as input and predicts the corresponding yield performance ŷ as output. This model is sensitive to three hyperparameters: the number of trees should be large enough to stabilize the error rate and small enough to be tractable; the number of features controls tree correlation, and the node size (minimum size of terminal nodes) determines the complexity of the individual trees. A tenfold CV was used to partition dataset to training and validation subsets. For each fold, we used the training subset for training and parameter tuning. A fivefold CV over train partition for each fold was applied to tune the parameters. Table 1 gives the values of these hyperparameters using a fivefold CV over the whole dataset to get the best values that lead to good performance on the validation dataset.
G × E interactions model. The random forest model has the capability to approximate nonlinear relationships among the variables. It grows many classification trees by randomly selecting subsets of features. As such, this model is ineffective in discovering specific combinations of features that have the most significant interactions. Therefore, we also introduced a combinatorial optimization-based model to augment the random forest by strategically searching for G × E interactions.
The G × E interactions model was designed to detect interactions among specific hybrid, location, and weather variables. This model is built off of a recently published algorithm 33 , which was designed to detect genetic interactions in the form of epistases. The algorithm was found to be effective in detecting multiple interactions involving multiple variables. The G × E interactions model considers yield as a linear function of input variables and their interactions, shown as follows.
• β 0 , β j , and b k are the effects of baseline, variable j, and interaction k, respectively, and • ǫ i is random noise for observation i.
In this model, the interactions are defined by a matrix α , which has a dimension of K × p , where K is the number of interactions that the proposed model tries to decipher and p is the number of variables. Each column of this matrix corresponds to a variable and each row corresponds to an interaction. Moreover, each element of matrix α can take three possible values 0, 0.5, 1. If α k,j = 0 , then interaction k requires that variable j be 0 ( X i,j = 0 ) for any individual i to receive this effect. If α k,j = 1 , then interaction k requires that variable j be 1 ( X i,j = 1 ) for any individual i to receive this effect. If α k,j = 0.5 , then variable j is not involved in interaction k. Given matrix α , the matrix Z can be subsequently calculated to determine whether or not the individuals receive the interactions. The dimension of the binary matrix Z is n × K , with each row corresponding to one individual and each column corresponding to one interaction. If Z i,k = 1 , then individual i receives the interaction k, and Z i,k = 0 otherwise. This complex relationship can be captured mathematically as: individual i receives interaction k (Z i,k = 1) if and only if X i,j + α k,j � = 1 , or equivalently X i,j = α k,j , for each variable j. The key to model (1) is to find Z from a given training dataset (X Train , y Train ) , which requires the estimation of the number of interactions and the combination of variables that are involved in each interaction. When Z has been determined, model (1) reduces to a multiple linear regression that is easy to solve and interpret. Figure 2 illustrates an over-simplified example of G × E interactions on corn yield. The given training data gives the yield of n = 8 corn plants with all possible combinations of p = 3 variables: high-yield (1) or low-yield (0) gene, fertile (1) or infertile (0) soil, wet (1) or dry (0) weather. No random noise was added to simplify the illustration. The figure shows the solution to the model (Eq. 1). Matrix Z has three columns, indicating three interactions.
The rest of the solution indicates that the baseline yield is β 0 = 2 , the high yield gene, and wet weather contribute additional β 1 = 1 and β 3 = 2 , respectively, and the fertile soil has no additive effect ( β 2 = 0). In our model, a similar approach is used to detect interactions among hybrid, soil, and weather at a much larger scale with n = 155, 765 and p = 1, 399 . To overcome the computational challenges, we used a similar heuristic algorithm as in 33 , which had three desirable features: (1) it used cross-validation to avoid-overfitting; (2) it was able to find local optimal solutions efficiently; and (3) it could be parameterized to balance computation time and solution quality.
Random forest model 2. Although the interaction model can decipher the interactions between binary predictors, it cannot find more complex nonlinear function of interactions. Hence, we feed the results of the G × E model into another random forest to identify more complex nonlinear interactions. Random forest model 2 was designed to predict the residual prediction from random forest model 1. Let ŷ 1 and ŷ 2 denote the predictions from random forest models 1 and 2. The overall model output ŷ 1 +ŷ 2 will provide a more accurate prediction of yield, y, than ŷ 1 if ŷ 2 can be trained to estimate y −ŷ 1 .
To achieve this objective, we feed matrix Z from the G × E interactions model to random forest 2 to predict not only linear G × E interactions described in matrix Z but also more complex and nonlinear interactions. This model is trained using the residual of y −ŷ 1 to improve its accuracy. The tuned hyperparameters for the random forest model 2 are reported in Table 2. The same process as the random forest model 1 was applied to tune hyperparameters.
The proposed model combines the strengths of combinatorial optimization in identifying G × E interactions and random forest in producing accurate predictions using complex and nonlinear functions. As such, it is a trade-off between insight and accuracy. It will be shown in the computational experiments that this hybrid model produced more insightful and accurate predictions than using either model alone.

Quantitative results
In this section, we report the results of our computational experiments, which were designed to test the performance of the proposed algorithm with respect to other benchmark approaches.
(1)  35 was implemented in Python by the authors. • An extreme gradient boosting tree (XGBoost) 36 model was trained using the xgboost 36 package in R, which was an efficient and scalable implementation of gradient boosting framework. Three hyperparameters were tuned using fivefold cross validation (without data leakage): "nrounds", "eta", and "gamma". • A G × E interactions model 33 was implemented in MATLAB (Version 2018a), which used heuristic algorithms to detect multi-way and multi-effect epistasis (interactions between binary variables). It is equivalent to the G × E interactions model without integrating with the random forest models. • A random forest 32 was trained using the ranger 37 packages in R, which was an ensemble of decision trees and trains with the bagging method, equivalent to the random forest model 1 without the interaction model and the random forest model 2 in our proposed model. Three hyperparameters were tuned using fivefold cross-validation: the number of trees, number of features, and node size.  Three metrics were used for evaluating and comparing the predictive models' performances: RMSE, which presents the difference between predicted and observed values, Mean Absolute Error (MAE), which measures the average magnitude of the prediction errors, without considering their direction, and R 2 , the coefficient of determination defined as the proportion of the variance in the response variable that is explained by independent variables. Because the ground truth of the test dataset was never released, we partitioned the training dataset into training and validation subsets in a tenfold CV manner. For each fold, we tuned the parameters and trained the models using the training set, and then their performances were evaluated using the validation set. We made sure that no validation data was leaked in the model training process. The average RMSE, MAE, and R 2 values over ten partitions for the six algorithms are reported in Table 3. These results indicate that the proposed model outperformed other algorithms in all measures. Since the random forest model was part of our proposed model and it outperformed the first four machine learning algorithms, these results indicated the effectiveness of both the random forest method and our G × E interactions detection model.
The performance of the proposed model is also illustrated in Fig. 3, which plots the average predicted yields against actual observations for all inbreds and testers. The results suggest that our proposed model's prediction is close to the observation, both on average and in terms of probability density distributions. Table 3. Average RMSE, MAE, and R 2 of six algorithms for yield prediction. A 10-fold cross-validation on the training dataset was used for algorithm performance evaluation, since the ground truth yield of the test dataset was never released.  , respectively. The predicted and observed average yield for the 14 inbred clusters and 13 tester clusters are summarized in Table 4.
Genotype and environment interactions. The proposed model was able to provide not only accurate yield prediction but also genotype and environment interactions that could be biologically insightful. Figures 4  and 5 show the two-way and three-way interactions between variables, respectively. The results indicate that weather variables involve in more interactions following soil and genotype.   Figure 6, which can help breeders predict the most promising crosses. The average yields for four combinations of crosses are given in Table 5. These results appear to suggest that testers have a slightly higher weight in determining the yield performance of their progeny.

conclusion
We proposed a new model to address the 2020 Syngenta crop challenge, which combines random forest with an G × E interactions model to predict yield performance of inbreds and testers based on historical yield data in multiple years and environments. Random forest model has been found to be an effective and powerful machine learning model for prediction, yet it has its limitations in the degrees and types of interactions among the predictors. Based on a recently published algorithm for detecting multi-way and multi-effect epistatic effects, the G × E interactions model captures both linear and nonlinear interactions of the genotype by environment effects. The combination of random forest and the G × E interactions model was found to be effective in predicting yield performances of inbred-tester combinations in our computational study using tenfold validation, achieving a 0.0869 validation RMSE, 0.0648 validation MAE, and 0.3448 R-squared value, outperforming four other popular machine learning algorithms as the benchmark. Moreover, our proposed model was also more explainable than other machine learning models by yielding genotype by environment interactions. Results from our proposed model will be able to help breeders test progeny and identify the best parent combinations to produce new hybrids with improved yield performances. www.nature.com/scientificreports/

Data availability
The data analyzed in this study was provided by Syngenta for 2020 Syngenta crop challenge. We accessed the data through annual Syngenta crop challenge. During the challenge, September 2019 to January 2020, the data was open to the public. Researchers who wish to access the data may do so by contacting Syngenta directly (https :// www.ideac onnec tion.com/synge nta-crop-chall enge/chall enge.php).