A weighted AMMI algorithm for nonreplicated data

The objective of this work was to propose a weighting scheme for the additive main effects and multiplicative interactions (AMMI) model, as well as to assess the usefulness of this W-AMMI model in the study of genotype x environment interaction (GxE) and quantitative trait locus x environment interaction (QxE) for nonreplicated data. Data from the 'Harrington' x TR306 barley (Hordeum vulgare) mapping population, with 141 genotypes evaluated in 25 environments, were used to compare the results from the AMMI model with those of two proposed versions of the W-AMMI model: equal weights per row and equal weights per column. The proposed W-AMMI columns algorithm is viable to analyze data with heterogeneous variance, when there are no replicates available. The use of the AMMI and W-AMMI models, in the indicated cases, improves QTL detection, besides providing a sound interpretation of GxE and a better understanding of QxE, which allows obtaining valuable information on increasing productivities in different environments.


Introduction
The genotype x environment interaction (GxE) and the quantitative trait locus (QTL) x environment interaction (QxE) are common phenomena in multienvironmental trials (METs), and they represent a major challenge for breeders who intend to develop more adapted genotypes to different environmental conditions.The modelling strategies that have been used to understand GxE and QxE are based on fixed effect models, such as regression techniques (Rodrigues et al., 2011;Pereira et al., 2012aPereira et al., , 2012b)), as well as on singular-value decomposition techniques (SVD) (Gauch Jr., 1992;Paderewski et al., 2011;Paderewski & Rodrigues, 2014), and on mixed effects models (Alimi et al., 2012).
The additive main effects and the multiplicative interaction (AMMI) (Gauch Jr., 1992) is the most widely used model to understand GxE and QxE.However, when the error variance is heterogeneous throughout the environments, or when data are contaminated (when the presence of even a single outlier, if extreme, may lead to misinterpretations), the use of the AMMI model might not be appropriate (Rodrigues et al., 2014(Rodrigues et al., , 2016)).In cases when error variance is not homogeneous across environments, Rodrigues et al. (2014) proposed a generalization of the AMMI model, by which a weighted linear model is used to model the main effects, and a weighted low-rank SVD based algorithm is used to model the multiplicative effects and the weighted AMMI or W-AMMI.Rodrigues et al. (2014) also proposed a weighting scheme based on the inverse of the error variance, for cases when data are partially replicated.Romagosa et al. (1996) and Gauch et al. ( 2011) evaluated the usefulness of the AMMI methodology to study the QxE, and to identify potentially involved QTLs in the control of the interaction, in order to identify specifically adapted genotypes for each environment, and to a better understanding of them.The idea is to perform QTL scans on the AMMI predicted values, instead of the scans on the observed phenotypic data (for instance, yield), aiming at increasing the scores of the logarithm of odds (LOD) of the QTL detections.
The objective of this work was to propose a weighting scheme for the additive main effects and multiplicative interactions (AMMI) model, as well as to assess the usefulness of this W-AMMI model in the study of genotype x environment interaction (GxE) and quantitative trait locus x environment interaction (QxE) for nonreplicated data.

Materials and Methods
The AMMI model is a unimultivariate method that uses the analysis of variance to estimate the additive main effects for genotypes (G) and for environments (E).Thus, the singular value decomposition is applied to the residuals of the Anova, in order to estimate the multiplicative interaction terms.The AMMI model can be written as follows: in which: Y ijk is the observed phenotypic value of the genotype i, in the environment j and in the block k; μ is the grand mean; α i is the genotype i main effects as deviations from μ; β j is the environment j main effects as deviations from μ; θ jk is the effect of the environment j in the block k; λ n is the unique value for interaction of the principal component (IPC) on the axes n, that is, the singular value; γ ni and δ nj are the IPC scores of the genotype i and the environment j for the axis n, that is, the left and right singular vectors, respectively; N ≤ min(I -1, J -1), with I representing the number of genotypes, and J, the number of environments; and ε ijk is the experimental error with normal distribution.Depending on the number n of terms (axes or principal components) retained to describe the pattern of the interaction, the model is denoted by AMMI0, AMMI1, ..., AMMIF.In the AMMI0, no axis of interaction is considered; in the AMMI1, only the first axis of interaction is considered, and so forth until AMMIF, where all N axes of interaction are considered.
The matrix formulation of the AMMI model can be given by in which Y is the (I X J) two-way table of genotypic means across trials, or environments.Each column of Y represents the vector of genotypic means, as obtained from the phenotypic analysis of a corresponding trial, by an appropriate mixed model analysis that accounts for experimental design features and spatial trends.
The additive part of the model contains the term 1 I 1 T J μ, an intercept term, being an (I X J) matrix with the grand mean μ in all positions; α I J T 1 , , an (I X J) matrix of genotypic main effects, as deviations from the grand mean (equal rows); and, 1 , I J T β an (I X J) matrix of environmental main effects, as deviations from the grand mean (equal columns).The interaction part of the model, is approximated by the matrix product UDV T , with U being an (I X N) matrix whose columns contain the left singular vectors of the interaction; D an (N X N) diagonal matrix containing the singular values of Y*; and, V an (J X N) matrix whose columns contain the right singular vectors of interaction (Rodrigues et al., 2014).
When the two-way data table has missing data, or when the error variance is not constant across the environments, the cells of the table should be weighted differently in the model, as they account for less information (in the case of missing values), or less reliable information (in the case of a larger error variance).To account for the heterogeneity of the error variance, Rodrigues et al. (2014) proposed the W-AMMI model that replaces the standard Anova by a weighted linear model, and the standard SVD, by a weighted SVD.This approach is based on the expectation-maximization (EM) algorithm, by which the sum of squares of the difference between two consecutive interactions, X (t+1) and X (t) , is greater than a small value.For instance, 10 -9 , we compute: X (t+1) = SVD(W ʘ Y + (1 -W) ʘ X t ) in which: W is an (I X J) matrix with weights W i,j , 0≤ W i,j ≤1; 1 is a (I X J) matrix of ones in all positions; is the Hadamard product of matrices; t is the iteration number; and X is a low-rank approximation, with rank (X)=N.The results of this procedure are the U N , D N , and V N matrices, so that  Y ≈ U N D N V' N , and is the rank of the approximation that needs to be decided prior to the estimation of the model parameters (Rodrigues et al., 2014).This means that the weighted AMMI models are not nested: for instance, the weighted AMMI2 model will have different PC1 scores from the AMMI1 model, although the differences are small.
By applying the weighted SVD as described in the previous equation to the matrix  Y , and replacing it in the matrix formulation of the AMMI model, it will result in the W-AMMI model.This generalization of the AMMI model takes into account the differences in error variances across environments and eventually missing cells, and it can be applied to all data sets in which the AMMI model is appropriate.A requirement for the application of the W-AMMI model is that the error variance in each environment must be computed.Consequently, replicated data per environment is required, at least partially (Rodrigues et al., 2014).
When replicated data is not available, either because the number of breeding lines is large and replicated experiments would be too expensive, or because the original replicated data is not made available as open data to the scientific community, its statistical analysis is more difficult and less reliable.In this case, it is not possible to compute the error variances across environments as described above, and reported in Rodrigues et al. (2014).In a preliminary attempt to adapt the weighted AMMI model for nonreplicated data, two cases were considered here: in the first one _ the W-AMMI rows model _ the weights are designated as the inverse of the variance of the environments across genotypes, that is, equal weights per row; and, in the second one _ the W-AMMI columns _ the weights are designated as the inverse of the variance of the genotypes across environments, that is, equal weights per column.For instance, if the variance for a given genotype across environments, or for a given environment across genotypes, is high, the weight for that genotype (or environment), in the final model, should be lower.Therefore, we are not exactly following the idea of Rodrigues et al. (2014), for whom the weights are proportional to the inverse of the error variance.Instead, we are down-weighing the genotypes (W-AMMI rows), or the environments (W-AMMI columns), with a higher-interaction variance, that is, the predicted values will be more similar for different genotypes or environments, which will affect the LOD scores of the QTLs based on the predicted values.However, if the contribution from error dominates the contribution from GxE, this approach might be a good proxy for the error variance.

Results and Discussion
For QTL detection and for the analysis of QxE, we used the 'Harrington' x TR306 barley mapping population (Tinker et al., 1996) (Figure 1), which includes 141 genotypes and 25 environments, as well as the information on 140 phenotypes and 127 markers.Although the original phenotypic data might be replicated, we used here the mean values because the online repository contains this information only.
Following the procedure proposed by Eastment & Krzanowski (1982) for choosing the number of components in the AMMI model, three components were selected to be retained, and the AMMI3 was considered.The first IPC axis captured 19% of the sum of squares of the GxE; the second one, 15%; and the third one, 8%; the other axes are responsible for 58% of the interaction.Consequently, the AMMI model with two components explains 34% of the interaction sum of squares, and the AMMI model with three components explains 42% of this interaction.For a better visualization, the AMMI2 biplot with the first two IPC is depicted, accounting for 34% of the interaction sum of squares (Figure 2).
In order to minimize any possible result distortion, if the variance due to the environment is high, a weighted analysis of the GxE can be performed to produce more potentially reliable results, for the visualization and QTL detection.For this purpose, smaller weights can be given to the environments with higher variance.Biplots for the models W-AMMI rows, AMMI, and W-AMMI columns are presented in Figure 2.After comparing these biplots, it is possible to see that the angles formed between the environments and the position of some genotypes undergo an insignificant In the comparison of the AMMI and W-AMMI columns, the angle changes between the genotypes and the environments are more evident.For instance, the genotype 4 is recommended to the environment SK92c by the AMMI model.However, the same genotype 4 shows a strong and positive interaction with the change, almost imperceptible.A small reduction of the angle between the genotype 59 and the environment AB92c, when comparing the two biplots, illustrates this result.This occurs due to the variation caused by the weights in the rows of the W-AMMI model.environment AB92c, when analyzing the W-AMMI columns model.Therefore, the most recommended genotype to the environment SK92c now is the number 80.We also noticed that genotypes 25 and 92 had a high correlation with the environment MB93, in the W-AMMI columns, and that these genotypes were not recommended to this environment with the AMMI model.
The results regarding the QxE did not consider the rows of the W-AMMI model because its predicted values are very similar to the ones by the standard AMMI model, and, therefore, their results were very similar.QTL scans to each of the 25 environments of the mapping population are presented (Figures 3  and 4 ).These figures include the QTL scans for the raw data, the QTL scans for the AMMI3 predicted values (AQ analysis), and the QTL scans for the "columns" of predicted values of W-AMMI3 (W-AQ columns).Table 1 shows the QTL detections for all environments, including the chromosome, the LOD score, and the position of the detected QTLs in the 'Harrington' x TR306 Barley mapping population, for the raw data, for the predicted values by the standard AMMI model, and for the predicted values by the W-AMMI columns model.
By analyzing the magnitude of the LOD scores in the QTL scans, we noticed a visible improvement when using the AMMI columns predicted values, not only in the number of detected QTLs, but also on the LOD scores, in comparison to the QTL scans on the raw phenotypic data and on the AMMI predicted values (AQ analysis).As expected, based on the results from Gauch et al. ( 2011), the QTL scans on the raw phenotypic data provide lower LOD scores than the QTL scans on the AMMI predicted values.
For instance, in the environment AK93, which was the most significantly associated with the interaction, the chromosomes 2 and 3 show QTLs only when the scan was made using the W-AMMI model.Yet, the chromosome 2 had a high-LOD score in the environment AB92b when the QTL scan was obtained with the AMMI or W-AMMI columns predicted values.The highest-LOD score values were found for the W-AMMI columns predicted values, and the biggest value was approximately 9.46 on chromosome 3, for the environment AK93.
Table 1 shows the number of detected QTL per chromosome and the mean LOD scores for those the clear QTLs in chromosomes 3 and 7 were only detected 4 and 9 times, respectively, when considering the QTL scans of the raw data, but 20 and 16 when considering the AMMI predicted values, and 25 and 16 times, when considering the W-AMMI columns predicted values.This reinforces the idea that the QTLs obtained from the AMMI predicted values gain strength from other environments, and that the proposed W-AMMI columns improves even further the QTL detections.
The weights for the W-AMMI algorithm can be chosen in accordance with set requirements in Smith et al. (2001), Möhring & Piepho (2009), and Welham et al. (2010).It is worth mentioning that the weights in the W-AMMI algorithm used here require a (re)scaling that leads them to values between zero and one.In specific cases with little (error) variance heterogeneity (Gauch et al., 2011), the standard AMMI model is totally appropriate.By directing the approach that will be used in a given experiment -the AMMI or W-AMMI models -, the error or residual variance for each environment should be calculated and, after that, checked for its homogeneity in all environments.If the error variance in the environments is homogeneous, the results from the AMMI model will be similar to those in the W-AMMI approach.Therefore, the AMMI standard model strategy is already sufficient.However, when the error variations have a high heterogeneity among the environments, the use of the AMMI model is not recommended, and the W-AMMI algorithm should be used (Rodrigues et al., 2014).
The techniques presented here to detect and understand GxE and QxE are based on statistical principles, with applicability in microbial and plant populations studied in various environments, and they can be adapted to genetic studies on animals and humans (Gauch et al., 2011;Rodrigues et al., 2014).

Conclusions
1.The proposed W-AMMI column algorithm is viable to analyze data that shows an heterogeneous variance, when there are not available replicates .
2. The use of the AMMI and W-AMMI models in the indicated cases improves the detection of quantitative trait loci, and provides a better understanding of the interaction between quantitative trait loci and the environment, which allows breeders to obtain a valuable information to increase crop productivity in different environments.
Table 1.Number of detected QTL per chromosome, and the mean LOD scores (data inside partentheses) for these detections, using the raw phenotypic data, the predicted values with the standard AMMI model with three multiplicative terms, and the predicted values with the W-AMMI column model with three multiplicative terms.

Figure 1 .
Figure 1.Genetic map data for the 'Harrington' x TR306 barley mapping population.

Figure 2 .Figure 3 .
Figure 2 .Biplots for the 'Harrington' x TR306 barley mapping population, with the models AMMI2, W-AMMI2 (weighting by rows), and W-AMMI2 (weighting by columns).The points represent the genotypes and the arrows correspond to the environments.

Figure 4 .
Figure 4. Quantitative trait loci (QTL) scans for the last nine environments of the yield data for the 'Harrington' x TR306 barley mapping population: dotted grey lines represent the QTL scans for the observed phenotypic data; gray lines represent the QTL scans considering the AMMI3 predicted values; and black lines represent the QTL scans obtained with the predicted values in the W-AMMI3 columns.All analyses were based on composite interval mapping.The codes above each scan represent the environments.