Improving genomic prediction of rhizomania resistance in sugar beet (Beta vulgaris L.) by implementing epistatic effects and feature selection

Experimental design and data preparation

The sugar beet population for this trial was developed by crossing two sugar beet lines and self-pollinating the resulting hybrids twice. This process resulted in a population of 155 S2 plants. These plants were genotyped using a customised SNP chip and subsequently self-pollinated. For each plant, 15 seeds were used as genotypes for this trial. Analysis of the SNP chip data revealed that each of the 155 genotypes was homozygous for resistance at both Rz1 and Rz2. Additionally, the population was expanded by including 15 seeds from a sugar beet line homozygous for resistance at Rz1 but not at Rz2.

Plants were grown for ten weeks in the greenhouse in soil infested with BNYVV, pathotype P. This variant of BNYVV contains five RNA strands⁴⁵ and is more aggressive than the variants with four RNA strands.¹⁷ After ten weeks, plants were removed from the soil and plant sap from lateral roots was extracted. Afterwards, the optical density (OD) value of each sample was measured using the double antibody sandwich enzyme-linked immunosorbent assay (DAS-ELISA). The OD values were measured after 60, 90, and 120 minutes using the Infinite F50^® (Tecan Group AG, Männedorf, Switzerland) at a wavelength of 405 nm. Harvest, sample preparation and conduction of the DAS-ELISA test followed the protocol described in.⁴⁶

Although DAS-ELISA is a commonly used tool to measure the concentration of BNYVV in samples from sugar beet,^7,10,23 it does not directly measure the virus concentration in form of the OD values.⁴⁷ To estimate the virus concentrations from the raw OD values as well as to reduce measurement errors for each 96-well plate, we transformed the non-normally distributed raw data to normally distributed data using an inverse logistic regression model.⁴⁸ The logistic regression model was derived using a serial dilution with 12 samples on each 96-well plate. The transformation followed the protocol in Ref. 48, with the adjustment that OD values were measured after 60, 90, and 120 minutes. At each of the three time points, the relationship between OD values and virus concentration was modelled using the logistic regression model in equation 1:

with OD_i being the OD value of sample i , C_i being the virus concentration of sample i, $\widetilde{\mathit{bc}}$ being the median of the buffer controls on each 96-well plate, tl being the technical limit of the machine, A describing the asymmetry of the curve, I being the relative virus concentration (C_i) at the inflection point if A=1, and S being the slope at the inflection point. Both A and S can be set freely with the lower limit of zero. Since the OD values were measured at three time points, each sample also provided three transformed values. Since the transformed data can be assumed to be normally distributed, the mean of the transformed data was calculated and used as the response variable to reduce technical errors during measurement at each time point. After transformation of the data, the mean of the transformed data was calculated for all plants descended from the same parent. Therefore, if no plants died during the trial, the mean from 15 plants was calculated as the phenotypic data point for the corresponding genotype.

After transformation and calculating the mean of the transformed values, the SNPs were prepared. SNPs with missing values were removed from the data set. Moreover, redundant SNPs were removed through linkage disequilibrium (LD) pruning. In this step, one of two SNPs that were correlated with more than an r² of 0.95 was removed to ensure that epistasis results were not confounded by LD.^33,49 Furthermore, SNPs with a major allele frequency of 0.95 or higher were removed as recommended in Refs. 50, 51 for genomic prediction studies. In this filtering step, it was ensured that only SNPs with a certain genomic variance in the population remained in the data set. After the final step of SNP filtering, 9,127 SNPs were kept in the data set. Finally, the remaining SNPs were recoded as 0 (homozygous major allele), 2 (homozygous minor allele), and 1 (heterozygous). This coding approach is recommended for analysing genotypic and additive genetic models.⁵² All filtering steps and the SNP recoding were performed using PLINK v1.90b6.10.⁵³

Genomic prediction and feature selection using single SNPs

After data were prepared, genomic prediction was performed using single SNPs. A total of 125 genotypes (representing 80% of the population) were randomly chosen as the training population, with 31 genotypes designated as the test population. This process, based on experimental designs used in previous studies involving feature selection,⁵⁴ was repeated ten times. Subsequently, genomic prediction was performed using random forest as it was recommended in the literature for genomic predictin in sugar beet.⁵⁵ To do so, the R package ranger, version 0.14.1, was used with default settings.⁵⁶

Prediction accuracy was assessed by predicting the test data set using a model derived from the training data set. Subsequently, the coefficient of determination (R²) was used to compare the predicted and the observed values in the test data set. The coefficient of determination is defined as the proportion of the explained variability of the total variability⁵⁷ and was used in previous studies as measure for the prediction accuracy.^58,59

To perform feature selection, the variable importance of each SNP was estimated using the function Boruta from the R package Boruta, version 7.0.0.⁶⁰ Boruta runs multiple iterations of random forest to assess the importance of each variable (in this context, SNPs). During these iterations, the importance of a variable is determined by measuring the decrease in prediction accuracy or increase in prediction error when the values of that variable are permuted or randomly shuffled. The more significant the decrease in accuracy, the higher the variable's importance. Boruta then compares the importance of each actual variable to the importance of randomly permuted versions (shadow variables).⁶⁰ This method has been used in earlier studies to evaluate SNP variable importance.^58,61,62 We chose the mean variable importance as estimate of the variable importance per SNP. We used a confidence level of 10−10, 2,000 trees and 200 as the maximum number of runs.

After the variable importance per SNP was estimated, feature selection was performed by carrying out genomic prediction with a random forest model that contained only the two SNPs with the highest variable importance. Subsequently, genomic prediction was performed using the three best SNPs, and so on. For each number of SNPs in the prediction model, the prediction accuracy was calculated. Since random forest can lead to different results in prediction accuracy even if the same data were used as training and test data set, prediction accuracy was estimated as the median prediction accuracy from ten repetitions. To prevent over-optimistic values for the prediction accuracy, variable importance was only estimated using the training data set.⁶³ In this way, genomic prediction using feature selection can be compared to genomic prediction using all available SNPs.^54,64

Performing the feature selection method as described above led to a variable importance per SNP as well as a prediction accuracy for each prediction model containing the i best SNPs for each of the ten random splits of the data set. In a next step, the prediction accuracy for each number of SNPs was defined as the median from these ten repetitions. Subsequently, the optimal number of SNPs was defined as the number of SNPs that maximised the median prediction accuracy from the ten repetitions.

Genomic prediction and feature selection using SNP pairs

In addition to genomic prediction with single SNPs, genomic prediction was also performed using SNP pairs. Since the genomic data set contained 9,127 SNPs after the data preparation and filtering, theoretically more than 41 million SNP pairs could be created out of these single SNPs. To reduce the resulting data set to a managable size, PLINK’s epistasis test was performed using default settings, testing the interaction term of each SNP pair for significance at a significance threshold of α=0.0001.⁶⁵ The selection of SNP pairs was performed with each training data set individually to prevent bias during the selection process.

After the SNP pairs were selected using PLINK’s epistasis test, the genotype of each SNP pair was defined using an additive-additive interaction model. In this way, the genotype of each SNP pair was defined as the product of both single SNPs.⁵² For instance, if any of the single SNPs was homozygous for the major allele (coded as 0 in the single SNPs), the genotype of the SNP pair was defined as 0. This was the case for five of the nine possible genotypes of a SNP pair. The combination of two heterozygous single SNPs would lead to a 1 for the SNP pair, the combination of a heterozygous SNP with a SNP that provides a homozygous minor allele would be 2, and the combination of two SNPs that provide homozygous minor alleles would be 4. The recoding of single SNPs as well as the resulting genotype of the SNP pair is summarised in Table 1.

After the SNP pairs were recoded, genomic prediction was performed with all SNP pairs that were selected from each training data set. Therefore, a prediction model was derived using random forest with each training data set (as described above with the ranger function with default settings), the phenotypic values of the corresponding test data set were predicted with the prediction model, and prediction accuracy was estimated as R² between the predicted and the observed values.

However, since PLINK’s epistasis test is based on linear regression and random forest is a method from machine learning, we developed an alternative for selecting SNP pairs from single SNPs based on machine learning methods. Therefore, we used the information about the variable importance of each single SNP as it was provided by the Boruta function and combined the single SNPs with the highest variable importance to all possible SNP pairs. Subsequently, genomic prediction was performed using all SNP pairs created with the best single SNPs and the resulting prediction accuracy was stored. This process was repeated for the three to 200 single SNPs with the highest variable importance. Finally, the number of single SNPs was determined where the resulting SNP pairs maximised the prediction accuracy. As with the other methods, this method was repeated with each training data set individually to avoid bias and afterwards the median from the ten repetitions was calculated for each number of analysed SNPs. An R script as well as the data from this trial are provided at https://github.com/tmlange/IFS\_SNPpairs.git⁴⁴ to give researchers the possibility to perform feature selection with SNP pairs based on the variable importance of single SNPs.

The feature selection with single SNPs and SNP pairs selected a certain number of best SNPs in each of the ten repetitions which were carried out independently from each other. The selected SNPs from each repetition were subsequently compared to assess their consistency across iterations. To evaluate the stability of the selected markers, a count was performed to determine the frequency of SNPs being chosen in these repetitions. SNPs that were selected at least 50% of the time were considered as robustly identified features. In this way, the results from the genomic prediction were used to detect SNPs that are associated with rhizomania resistance.

Genomic prediction and feature selection using single SNPs

The median prediction accuracies for all methods described are presented in Table 2. First, genomic prediction was conducted using all 9,127 single SNPs that remained after filtering. Genomic prediction with these single SNPs across ten random splits of the data set resulted in a median prediction accuracy of R² = 0.146.

In addition to genomic prediction with all SNPs, incremental feature selection was performed to optimise prediction accuracy by selecting a subset of the most informative SNPs. Figure 1 illustrates the number of SNPs in the prediction model on the X-axis and the corresponding median prediction accuracy from the ten repetitions on the Y-axis. It is evident that prediction accuracy increases steeply with the inclusion of the initial SNPs. However, just above R²=0.25, the prediction accuracy peaks and then gradually decreases, with the decline being less steep than the initial rise. Using this approach, the optimal set of SNPs for genomic prediction was identified. The prediction accuracy was maximised when 29 SNPs were included in the model, resulting in a median prediction accuracy of R²=0.267.

Figure 1. Median of the R² values from the ten repetitions of genomic prediction using random forest with the 2, … , 9,127 SNPs.

Genomic prediction and feature selection using SNP pairs

Besides genomic prediction and feature selection with single SNPs, similar approaches have been performed using SNP pairs. To perform genomic prediction with SNP pairs, PLINK’s epistasis test was performed with default settings for each training data set individually. After filtering via PLINK’s epistasis test, the resulting sample sizes ranged from 46,556 to 87,529 SNP pairs. Taking the 41 million theoretically possible SNP pairs into consideration, this is a reduction to 0.1% to 0.2%. When genomic prediction was performed with all SNP pairs that were left after filtering, the median prediction accuracy was R²=0.191.

Furthermore, feature selection was performed to identify an optimal subset of SNP pairs for genomic prediction. The best single SNPs, judged by their variable importance, were used to produce SNP pairs, and incremental feature selection was conducted using random forest with the SNP pairs derived from the top 3 to 200 single SNPs. Figure 2 displays the median prediction accuracy from the ten repetitions on the Y axis and the corresponding number of single SNPs that make up the SNP pairs on the X axis.

Figure 2. Median of the R² values from the ten repetitions of genomic prediction using random forest when the 3, … , 200 best single SNPs are combined to SNP pairs.

One can see that the prediction accuracy increases steeply when a small number of SNP pairs are included in the prediction model. Similar to Figure 1 that displays the prediction accuracy with the single SNPs, the prediction accuracy with the SNP pairs also forms a peak and decreases from there on. One can see that the peak height is slightly above R²=0.3 when SNP pairs are created using the 16 best SNPs.

Considering Figure 2, it appears that the number of SNP pairs has an effect on the resulting prediction accuracy. This suggests that the prediction accuracy might be affected by the number of SNP pairs that is selected via PLINK’s epistasis test. Therefore, we have analysed how many SNP pairs were selected via PLINK if the significance threshold was modified. However, the resulting prediction accuracy remained unchanged with different thresholds for PLINK’s epistasis test (tested for thresholds 10⁻², … , 10⁻⁸, data not shown). Consequently, we conclude that although the number of selected SNP pairs can be easily adjusted in PLINK’s epistasis test, this selection does not affect the resulting prediction accuracy.

Finally, we counted how often SNPs were selected in the ten repetitions of the feature selection. The four SNPs “SNP0425”, “SNP2484”, “SNP6428”, and “SNP7343” were selected at least 50% of the time in the feature selection in the ten repetitions. The SNPs “SNP2484” and “SNP6428” are located on chromosome 3. However, “SNP6428” is located on chromosome 2 and “SNP7343” is located on chromosome 5.

After identifying these four SNPs, SNP pairs were created as described in Table 1 with “SNP2484” as one of the two SNPs that were located on chromosome 3 together with “SNP6428” which was located on chromosome 2 as well as “SNP2484” together with “SNP7343” which was located on chromosome 5. Figure 3 displays the virus concentration depending on the four different genotypes resulting from the two SNP pairs.

Figure 3. The virus concentration of the plants in the trial based on the different genotypes of the SNP pair “SNP2484-SNP6428” on the left and the SNP pair “SNP2484-SNP7343” on the right.

The genotype of the SNP pair was defined as the additive-additive SNP-interaction model in Ref. 52.

In both graphs in Figure 3, it is evident that the genotype 0 (both SNPs homozygous for the major allele) resulted in the lowest virus concentrations, while genotype 4 (both SNPs homozygous for the minor allele) led to the highest virus concentrations. Analysing the effect of the genotypes on the virus concentration with an ANOVA produced highly significant p values for both SNP pairs (SNP2484-SNP6428: p = 6.6 • 10⁻⁸; SNP2484-SNP7343: p=1.1• 10⁻⁷) and R² values of R² = 0.2115 for SNP2484-SNP6428 and R² = 0.2066 for SNP2484-SNP7343. Thus, more than 20% of the total variability in the data can be explained with each SNP pair individually.