ORIGINAL RESEARCH published: 06 February 2020 doi: 10.3389/fpls.2020.00025 Exploring Deep Learning for Complex Trait Genomic Prediction in Polyploid Outcrossing Species Laura M. Zingaretti 1*, Salvador Alejandro Gezan 2, Luis Felipe V. Ferrão 3, Luis F. Osorio 4, Amparo Monfort 1,5, Patricio R. Muñoz 3, Vance M. Whitaker 4 and Miguel Pérez-Enciso 1,6* 1 Centre for Research in Agricultural Genomics (CRAG) CSIC-IRTA-UAB-UB, Campus UAB, Barcelona, Spain, 2 School of Forest Resources and Conservation, University of Florida, Gainesville, FL, United States, 3 Blueberry Breeding and Genomics Lab, Horticultural Sciences Department, University of Florida, Gainesville, FL, United States, 4 IFAS Gulf Coast Research and Education Center, University of Florida, Wimauma, FL, United States, 5 Institut de Recerca i Tecnologia Agroalimentàries (IRTA), Barcelona, Spain, 6 ICREA, Passeig de Lluís Companys 23, Barcelona, Spain Genomic prediction (GP) is the procedure whereby the genetic merits of untested candidates are predicted using genome wide marker information. Although numerous Edited by: examples of GP exist in plants and animals, applications to polyploid organisms are still Xiyin Wang, North China University of Science scarce, partly due to limited genome resources and the complexity of this system. Deep and Technology, China learning (DL) techniques comprise a heterogeneous collection of machine learning Reviewed by: algorithms that have excelled at many prediction tasks. A potential advantage of DL for Kun Lu, GP over standard linear model methods is that DL can potentially take into account all Southwest University, China Khanh N. Q. Le, genetic interactions, including dominance and epistasis, which are expected to be of Taipei Medical University, Taiwan special relevance in most polyploids. In this study, we evaluated the predictive accuracy of *Correspondence: linear and DL techniques in two important small fruits or berries: strawberry and blueberry. Laura M. Zingaretti laura.zingaretti@cragenomica.es The two datasets contained a total of 1,358 allopolyploid strawberry (2n=8x=112) and Miguel Pérez-Enciso 1,802 autopolyploid blueberry (2n=4x=48) individuals, genotyped for 9,908 and 73,045 miguel.perez@uab.es single nucleotide polymorphism (SNP) markers, respectively, and phenotyped for five Specialty section: agronomic traits each. DL depends on numerous parameters that influence performance This article was submitted to and optimizing hyperparameter values can be a critical step. Here we show that Bioinformatics and interactions between hyperparameter combinations should be expected and that the Computational Biology, a section of the journal number of convolutional filters and regularization in the first layers can have an important Frontiers in Plant Science effect on model performance. In terms of genomic prediction, we did not find an Received: 22 October 2019 advantage of DL over linear model methods, except when the epistasis component Accepted: 10 January 2020 Published: 06 February 2020 was important. Linear Bayesian models were better than convolutional neural networks for Citation: the full additive architecture, whereas the opposite was observed under strong epistasis. Zingaretti LM, Gezan SA, Ferrão LF, However, by using a parameterization capable of taking into account these non-linear Osorio LF, Monfort A, Muñoz PR, effects, Bayesian linear models can match or exceed the predictive accuracy of DL. A Whitaker VM and Pérez-Enciso M (2020) Exploring Deep Learning for semiautomatic implementation of the DL pipeline is available at https://github.com/ Complex Trait Genomic Prediction in lauzingaretti/deepGP/. Polyploid Outcrossing Species. Front. Plant Sci. 11:25. Keywords: genomic prediction, genomic selection, polyploid species, deep learning, epistasis, complex traits, doi: 10.3389/fpls.2020.00025 strawberry, blueberry Frontiers in Plant Science | www.frontiersin.org 1 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids INTRODUCTION and, in contrast to DL, have not been designed to model non- additive genetic effects (such as dominance or epistasis); Deep learning (DL) techniques comprise a heterogeneous however, these effects can be incorporated in the model with collection of machine learning algorithms which have excelled appropriate parameterizations. at many prediction tasks, and this is a very active area of research One potential advantage of DL for GP over standard methods is (Min et al., 2017; Pattanayak, 2017; Namin et al., 2018). All DL that the whole genetic merit, including all non-additive effects, can algorithms employ multiple neuron layers and numerous potentially be predicted without the need to partition all effects. This architectures have been proposed: multiple layer perceptrons is an interesting property for clonally propagated outcrossing (MLPs), recurrent neural networks (RNNs), convolutional species, because genomes can be asexually reproduced from single neural networks (CNNs) (LeCun et al., 2015) and others. DL is plants once the desirable individual is found. It should also be a relatively straightforward to implement (https://keras.io/why- promising strategy in polyploids, although their complex genetic use-keras/), but optimum performance depends on an structure has delayed the availability of whole genome markers and adequate hyperparameter choice, which is not trivial and of specific analytic tools for, e.g. SNP calling (Slater et al., 2016; requires considerable computational resources (Young et al., Gezan et al., 2017; Bourke et al., 2018). A few studies have 2015; Chan et al., 2018). Although previous, limited evidence demonstrated the potential advantages of GS in allo and does not show a consistent advantage of DL over penalized linear autopolyploids (Gezan et al., 2017; Enciso-Rodriguez et al., 2018; methods for genomic prediction (GP) purposes (González- Nyine et al., 2018; Amadeu et al., 2019; de Bem Oliveira et al., 2019; Recio et al., 2014; Ma et al., 2017; Bellot et al., 2018; Juliana et al., 2019; Zingaretti et al., 2019), although its Montesinos-López et al., 2018a; Montesinos-López et al., implementation is still in its infancy. 2018b; Montesinos-López et al., 2019a), more efforts are When non-additive effects are investigated, there are two needed to fully understand the behavior and potential important points that need to be considered for higher ploidy constraints and capabilities of DL in GP scenarios. levels: i) there is a portion of the intra-locus allele interaction (i.e., Genomic selection (GS) is the breeding strategy consisting in dominance) that is passed to the progeny (particularly full-sibs), predicting complex traits using genomic-wide genetic markers. and ii) the definition of non-additive effects is more complex The idea was developed to overcome the limitations of marker- than in diploids as higher order interaction exist (Osborn et al., assisted selection (MAS) and was formalized by Meuwissen et al. 2003). Thus, methodologies that could model the whole genetic (2001). While MAS establishes a model with only the markers merit without restrictive assumptions could facilitate and with significant associations, genomic selection includes all, or improve the prediction for polyploid species, making DL an most available markers, for GP, irrespective of their effect and its attractive choice for genomic prediction. In practice, DL aims at significance. Due to the decrease in genotyping costs, genomic predicting the whole genetic merit, including interactions selection is becoming the standard tool in many plant and animal irrespective of their origin. breeding programs (Bernardo, 2008; González-Camacho et al., Among the polyploid species, strawberries (Fragaria x 2012; Crossa et al., 2013b; Meuwissen et al., 2013; Wiggans et al., ananassa) and blueberries (Vaccinium corymbosum) are 2017). There is an increasing number of successful applications of considered two of the most important soft fruit commodities. genomic selection in diploid and polyploid organisms where its use Considered a rich source of vitamins and minerals, fruit markets has generated important genetic gains by improving the accuracy of for both species have experienced a global increase in production breeding value prediction and dramatically reducing generation and consumption over the past decade (https://www.nass.usda. intervals (Crossa et al., 2013a; Castillo-Juárez et al., 2015; Duangjit gov/Publications/Todays_Reports/reports/ncit0619.pdf). To et al., 2016; Juliana et al., 2019; de Bem Oliveira et al., 2019). ensure that production and fruit quality meet the global In any scenario, GP poses statistical challenges since the demand, genetic improvement, and particularly GP, has a role number of markers is usually much larger than the number of to play in maximizing the utility, diversity, and yield of resources. individuals, i.e., the so-called large p (number of features) small n In this sense, previous experimental assessments performed in (sample size) paradigm (de los Campos et al., 2013; Pérez and De blueberry (Amadeu et al., 2019; de Bem Oliveira et al., 2019) and Los Campos, 2014). In this context, statistical methods require strawberry (Gezan et al., 2017) have proven the feasibility of either shrinkage, variable selection, or a combination of both incorporating genomic selection to either accelerate the pace or (Tibshirani, 1996). Most GP methods are based on linear models, improve the efficiency of breeding programs. From a genetic such as Genomic Best Linear Unbiased Prediction (GBLUP) standpoint, one important difference between both species is its (VanRaden, 2008), the Bayesian GP family (Meuwissen et al., inheritance pattern. Cultivated strawberry (Fragaria x ananassa) 2001; Pérez and De Los Campos, 2014), or LASSO (Tibshirani, is an allo-octoploid hybrid plant originated by cross between two 1996). In GBLUP, all marker effects are assumed to be normally wild octoploid species F. chiloensis and F.virginiana (Hancock distributed with equal variance and a homogeneous shrinkage is et al., 2008) both descendants of Fragaria diploid species; induced, whereas Bayesian models are more flexible and referred as allopolyploids, meiosis is mainly dictated by differential shrinkages and/or variable selection can be applied preferential pairing, exhibiting a diploid-like (or disomic) to distinct marker subsets. Note that these methods are linear segregation. In contrast, blueberry is a tetraploid organism Frontiers in Plant Science | www.frontiersin.org 2 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids originated from genome duplication within the same species. In 2015). After quality control, in which those markers with minor autopolyploids, the meiotic pairing is mainly described by forming allele frequencies (MAF) < 5% and with missing marker either random bivalents or multivalent during the division. Since the data > 5% were eliminated, 9,908 polymorphic SNP markers molecular mechanisms in auto and allopolyploids are quite were available. A total of five yield and fruit quality traits were complex, comparing new algorithms is a relevant issue to the evaluated in each trial: soluble solid content (brix), average fruit prospect of GP in these and other polyploid species. weight (AveWtT), total marketable yield (MktWtT), early In this study, we evaluated the performance of deep learning marketable yield (MktWtE), and percentage of culled for genomic prediction in two important horticultural species: (unmarketable) fruit (CullsTPer). Additional details for T2 and allo-octoploid strawberry and auto-tetraploid blueberry. We T4 can be found in Gezan et al. (2017). complement the empirical study with simulations to The blueberry population used in this study encompasses one understand better the impact of genetic architecture on DL cycle of the University of Florida blueberry breeding program's performance. Given the complexity of implementing DL, we recurrent selection and comprised 1,802 lines from 117 full-sib also provide a guideline on best practices for hyperparameter families. The population was originated from 146 parents that tuning and evaluate its importance in terms of predictive ability. presented superior phenotypic performance (cultivars and To facilitate reproducibility of these methods, a python-based advanced stage of breeding). Individuals were evaluated for five package for semiautomatic DL implementation, including auto yield and fruit quality-related traits: firmness, fruit size, weight, and allopolyploid organisms have been made available at https:// yield, and picking scar, which were collected during two github.com/lauzingaretti/deepGS/. production seasons. Phenotypes were pre-corrected for fixed year effects, as detailed in Amadeu et al. (2019) and de Bem Oliveira et al. (2019). A total of 73,045 SNPs was obtained using sequence capture by Rapid Genomics (Gainesville, FL), after MATERIALS AND METHODS aligning the reads against the high-quality “Draper” genome Plant Material and Genotypes assembly (Colle et al., 2019) as described in Benevenuto et al. Predictive performances were compared in two polyploid species (2019). Marker filtering followed these criteria: biallelic, mean (blueberry and strawberry), for a series of traits with presumably coverage > 40, minimum allele frequency > 0.01; maximum contrasting genetic architecture. A summary of both experimental missing data = 0.5%; minimum quality = 20. Also, individuals data sets is presented in the Table 1. with more than 50% missing data were removed, missing Regarding strawberry, we used 1,233 unique genotypes which genotypes were simply imputed with the mean. Tetraploid correspond to five advanced selection trials (T2, T4, T6, T8, and genotypes were called and the allele dosages were inferred with T10) from the strawberry breeding program at the University of the updog R package (Gerard et al., 2018). Standard genotype Florida, Institute of Food and Agricultural Sciences (USA). These calling with updog allows inferring genotypes according to the advanced trials were planted in five consecutive seasons and were number of allele copies, and genotypes can be coded say 0,1,2,3,4. given an even code starting with season 2013–2014 as T2 and In addition, as in de Bem Oliveira et al. (2019), here we ending with season 2018–2019 as T10. The number of lines in considered a set of “diploidized” genotypes that were obtained each trial was 217, 240, 236, 272, and 393 for T2, T4, T6, T8, and pooling all heterozygous genotypes in a single class, i.e., T10, respectively. Some of the genotypes in the last trial T10 were genotypes above 0,1,2,3,4 can be recoded as 0,1,1,1,2. The already tested in earlier trials, making the total number of rationale is that there can be incertitude on the exact number observations sum up to 1,358 (instead of 1,233). Plants were of allele copies in heterozygous genotypes. genotyped with the Axiom IStraw90 SNP array (Bassil et al., The GP methods evaluated in this study were assessed by true validation, which was obtained by splitting data into a training and a validation dataset. In the strawberry dataset, we considered that predicting performance of the last stage lines (T10) is the TABLE 1 | Summary of blueberry and strawberry experimental data sets used in most interest application for the industry and therefore the this paper. population was divided between training (T2, T4, T6, and T8 Strawberry (allopolyploid) Blueberry trials) and validation (T10) subsets with 965 and 393 lines, (autopolyploid) respectively. In the case of blueberry data, all samples were Ploidy 2n = 8x = 112 2n = 4x = 48 equally important and 30% of randomly sampled genotypes No. 1,358 (1,233 unique genotypes) 1,802 were assigned to the validation set. Predictive ability (PA) was observations defined as the correlation between observed and predicted No. SNPs 9,908 73,045 phenotypes in the validation set; prediction was computed Traits • Soluble solid content (brix) • Firmness from parameters estimated in the training dataset only. analyzed • Average fruit weight (AveWtT) • Fruit Size • Total marketable weight (MktWtT) • Weight • Early marketable yield (MktWtE) • Yield Genetic Structure and Heritability • Percentage of culled fruit (CullsTPer). • Scar Main Gezan et al. (2017) Amadeu et al. (2019) Inference reference and de Bem Oliveira Potential genetic structure was assessed by principal component et al. (2019) analysis (PCA) using all genotypes. Since genetic architecture Frontiers in Plant Science | www.frontiersin.org 3 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids may have an impact on GP performance and on the optimum The above parameterization assumes additivity of effects, GP model (Daetwyler et al., 2013), additive and non-additive although linear models can address non-linear relationships if genetic features were assessed by computing variance properly parameterized. Non-linear interactions can be modeled components from the model: by expressing g (Equation 2) in a general way, i.e., g = W w where W (centered and scaled) is a matrix of dummy variables that indicates y=m1+a+d+e+e (1) the number of copies of the reference allele ranging from 0 to the where the vector y represents the adjusted phenotype, m1 is the ploidy level (Slater et al., 2016; Enciso-Rodriguez et al., 2018). This intercept, a∼N(0, As2a ),d∼N(0, Ds2d ) and e∼N(0, Es2e ) are the model is, in principle, a good parameterization to account for non- additive, dominant and epistatic effects, respectively, and e ∼ linear interactions and we will refer to it as BRR general model (BRR-GM), since Bayes ridge regression was used. For more details, N(0, Is2e ) is the residual component. Matrices A and D were see Enciso-Rodriguez et al. (2018) and Amadeu et al. (2019). obtained using AGHmatrix package (Amadeu et al., 2016) for Non-linearity can also be managed by means of RKHS both strawberry (as diploid) and blueberry (autotetraploid) regression (Gianola et al., 2006) as an alternative to a linear species. For diploids, A and D were computed using regression for capturing complex interactions. This model VanRaden (2008) and Vitezica et al. (2013) methods, ZZ0 respectively. In fact, A = where Z is the matrix considers g in Equation 2 as N(0, Ks2g ) with K(xi , xi0 ) = exp( − 2oj pj (1 − pj ) that contains the centered individual genotype values and D = hjjxi − xi0 jj2 ), a kernel function where h is de bandwidth p MM0 is the dominance matrix, where the M elements parameter controlling how fast the covariance function drops 4 o − pj )2 j ½pj (1 with the distance between pairs of markers and jjxi − xi jj2 is the are − 2p2j , 2pj (1 − pj ), −2(1 − pj )2 for genotypes 0, 1, and 2, Euclidean distance between any two pairs of genotypes. This respectively. In the case of ploidy = 4, D was obtained as in parameterization induces a general matrix of genetic covariance Slater et al. (2016). The epistatic matrix (E) considered is the between markers. The key point here is that the kernel can model Hadamard product of additive effects, i.e. A ⊙ A (Henderson, non-linear relationships because it is a non-linear transformation 1984) Posterior distributions of genetic parameters were of the distances between the input variables. Empirical evidence obtained using Reproducing Kernel Hilbert Spaces (RKHS) confirms that it is an accurate approach to predict phenotypes of regression with BGLR package (Pérez and De Los Campos, complex traits (Gianola et al., 2008; de los Campos et al., 2009; de 2014). The additive, dominance and epistatic ratios were Los Campos et al., 2010). BRR-GM and RKHS were only estimated as: h ^2a = s2a =(s2a + s2 + s2e + s2e ), h ^2 = s2 =(s2a + s2 + s2e + s2e d d d d implemented for strawberry and simulated scenarios, since it ) and h^2e = s2e =(s2a + s2 + s2e + s2e ); where s2i the ith mean posterior was in strawberry where we found the trait with the largest d estimates of s2 as in Equation 1. We used both training and epistasis component, as described below. validation datasets combined in this stage, since this is purely a descriptive analysis and the values obtained are not employed in the later prediction stages. Deep Learning (Convolutional Neural Networks) Penalized Linear Methods DL has been described as a universal learning approach able to We compared the prediction performance of DL models with solve supervised, semi-supervised and unsupervised problems. two well-established linear methods: Bayesian Lasso (BL, Several DL architectures have been proposed, such as MLPs, Meuwissen et al., 2001) and Bayesian Ridge Regression (BRR, RNNs, CNNs, Generative Adversarial Networks (GANs) and Gianola, 2013). In these models, the trait can be expressed as: Reinforcement Learning (RL). Figure 1 shows a generic pipeline to evaluate DL in a GP context. y=m1+g+e (2) In our previous experiment (Bellot et al., 2018), CNNs were the best performing methods and therefore are the only ones discussed where m1 is the overall mean, g=Xb, X is the genotypes' matrix here. The advantage of CNNs in a GP context is that they can model and b is a vector of marker effects. In BRR, prior distributions of the correlation between adjacent input variables, that is, linkage marker effects b are N(0, Is2b ), whereas the prior distributions for disequilibrium between nearby SNPs. This is done via a l mathematical operation called convolution (Goodfellow et al., b in BL have a Laplace distribution, i.e., p(bj l, s2e ) = 2 exp( − 2se 2016). A typical CNN is made up of “convolutional layers”, l “pooling”, “flatten” and “dense” fully connected layers (Figure 2). jbj 2 ). Note that the Laplace distribution does not remove se In the “convolutional layer”, an operation called convolution is markers so, contrary to its frequentist counterpart, BL is not a performed along the input of predefined width and strides, which variable selection approach. Each model was fitted by using only are known as “kernel” and “filter” in the DL jargon, respectively. phenotypes from the training subset. The models were run using From a mathematical view, a convolution s(t) is a function that can the BGLR package (Pérez and De Los Campos, 2014) with a be defined as an “integral transform” (Widder, 1954): Gibbs sampler algorithm for a total of 6,000 cycles, discarding the first 1,000 samples for burn-in. sðt Þ = ð f *kÞðt Þ = ʃ kðt − xÞf ðxÞdx (3) Frontiers in Plant Science | www.frontiersin.org 4 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids FIGURE 1 | A generic deep learning (DL) pipeline for genomic prediction (GP) purposes. The general process includes the training and validation steps. In the training step, data are split into training and testing, DL hyperparameters are optimized by internal cross-validation with the test set and the model with the best predictive ability (PA) is chosen. In the validation step, the model PA is evaluated using a new set of data. FIGURE 2 | General CNN architecture employed in our workflow. The input layer is a SNP matrix of size n x p, where n is the number of training set and p, the number of SNPs. The convolutional layer consists on a nfilters convolution followed by a max-pooling layer with poolsize = 3 and an optional dropout. The outputs of max-pooling layer are joined together into one vector by flattening. All the neurons in the flatten layer are fully connected to the first dense layer. We tune the network using i dense layers with a variable number of hidden neurons in the respective hidden layers. The output of these dense layers is the prediction layer that uses linear function as activation. The neurons in convolutional and dense layers use relu, tanh or linear function as activations. where one of the functions (k or f) in Equation (3) must be a sðt Þ = ð f *kÞðt Þ = ox kðt − xÞf ðxÞ (4) kernel. Assuming that the kernel is represented by k, the convolution is the transformation of f (input data in the DL One of the main advantages of convolution networks is their context) into s(t). The operation is just the weighted sum of an capability to reduce the number of operations, i.e., the infinite number of copies f shifting over the kernel. The discrete hyperparameters to be estimated. As usual, an activation version of Equation 3 follows naturally as: function (generally non-linear) is applied after each convolution Frontiers in Plant Science | www.frontiersin.org 5 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids to produce the output layer. Finally, “pooling” layers reduces dimension product of individual genotypes in each pair. Total genetic and achieves a smoother representation, summarizing adjacent neurons value was the sum of effects across pairs of loci. by computing their maximum or mean. 3. Mixed: 80 individual additive SNPs and 60 epistatic SNP pairs were randomly chosen. Total genetic value was the sum Hyperparameter Optimization of effects across pairs of loci and individual additive loci. Since DL depends on numerous parameters that influence performance, optimizing hyperparameter can be a critical Allele substitution effects were sampled from a gamma unresolved step, which relies heavily on heuristics. Hence, it is distribution G(a = 1, b = 0.2). The trait was obtained adding surprising that many DL applications in GP have not paid the genetic value to an environmental normal residual. enough attention to this problem (Ma et al., 2017; Montesinos- Environmental variance was chosen such that broad-sense López et al., 2018b; Montesinos-López et al., 2019b). Several heritability was set to 0.50. For each genetic architecture, five approaches have been proposed for hyperparameter tuning (e.g., replicates were run. We compared BRR, BRR-GM, RKHS, and Bellot et al., 2018; Cho and Hegde, 2019; Le et al., 2019; DL. DL architectures were specifically optimized to each Rajaraman et al., 2019; Yoo, 2019). Here, DL architectures phenotypic trait, since no universal architecture is able to make were optimized using Talos (Autonomia Talos, 2019), which accurate predictions for all cases. works combining all parameters in a grid. Talos can choose the best model either maximizing the predictive accuracy or minimizing the error; the former criterion was employed here. RESULTS Since the approach can be expensive as the number of hyperparameters increases, a random search is the best strategy Population Structure and Genetic in practice. This rule evolves a list of CNN models for each Parameters phenotypic trait. We optimized the following hyperparameters No clear population structure was observed, neither in the (values considered within parentheses): activation function (relu, strawberry nor in the blueberry datasets (Figure S1). Note that tanh, linear), number of filters (16, 32, 64, 128), regularization genetic relationships between trials in strawberry data are rather (i.e., weight decay in DL terminology, 0, 0.1, 0.01, 0.001), learning uniform, irrespective of whether they are successive seasons or rate (0.1, 0.01, 0.001, 0.0025), number of neurons in fully not. This, together with the fact that little genotype by connected layer (4, 8, 12, 16), number of hidden layers environment (or year) interaction was observed (Gezan et al., (1,5,10), and dropout (0, 0.01, 0.1, 0.2). 2017), suggests a favorable scenario for GP. Talos output is the accuracy for each hyperparameter Heritability estimates in strawberry are slightly different from combination; we then used hyperparameter values as those obtained in the same material by Gezan et al. (2017) since independent variables and accuracy as target variable to run a here we used additional data and we removed genotypes tested random forest algorithm, which allowed us to compute the since here we used additional material and we removed hyperparameter value importance, measured as the decrease in genotypes tested more than once on different seasons. Gini's coefficient when adding the given hyperparameter. This Nevertheless, in agreement with previous results (Amadeu hyperparameter importance can be then used as guide to et al., 2019; de Bem Oliveira et al., 2019; Gezan et al., 2017) improve interpretability. The R package randomForest (Liaw narrow-sense heritabilities were moderate, ranging from 0.25 to and Wiener, 2002) was employed for this analysis. 0.35 for most strawberry (Figure 3) and blueberry (Figure 4) The DL algorithms used in this study were implemented in traits, except for strawberry average fruit weight (h2a = 0:58) The Keras (Chollet, 2015) and Tensorflow (Abadi et al., 2015) and degree of dominance found was quite low in general, especially in run on a GPU equipped Linux workstation. A generic script is strawberry. An exception was blueberry yield, where dominant publicly available at https://github.com/lauzingaretti/deepGS/. and epistatic variances were similar to the additive variance (Figure 4E). A remarkable case is percentage of culled fruit Simulation (CulsTPer) in strawberry, where the epistatic ratio (18%) was We studied the impact of genetic architecture on prediction only slightly smaller than the additive one (25%, Figure 3E). performance by simulation using the actual observed strawberry genotypes, assessing predictive performance with the same T10 Hyperparameter Importance strawberry genotypes (and genotypic data) as in the real CNN hyperparameters were optimized for each strawberry trait experiment, except that phenotypic responses were simulated. separately. Figure 5A shows the importance of each Three contrasting genetic architectures were considered: hyperparameter obtained from random forest by regressing the model predictive accuracies (obtained by an inner cross- 1. Additive: 200 randomly chosen SNPs were considered as validation) on all hyperparameter values combinations. causal loci. No dominance was simulated. Total individual Interestingly, the number of filters was overall the most genetic value was the sum of effects across loci. relevant factor, whereas other factors such as learning rate, whose 2. Epistatic: 100 epistatic pairs of SNPs were randomly sampled. importance has been claimed in the literature as critically important Epistasis was multiplicative by pairs, i.e., the genotype was the (Maas et al., 2013; Bawa and Kumar, 2019; Feng and Lu, 2019), Frontiers in Plant Science | www.frontiersin.org 6 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids FIGURE 3 | Posterior distributions of additive (blue), dominant (red), epistasis (green) fractions of variance in octoploid strawberry: (A) soluble solid content (brix); (B) early marketable yield (MktWtE); (C) total marketable yield (MktWtT); (D) average fruit weight (AveWtT); and (E) percentage of culled fruit (CullsTPer). Note the scale may vary along traits. FIGURE 4 | Posterior distributions of additive (blue), dominant (red), epistasis (green) fractions of variance in blueberry obtained with the tetraploid genotypes: (A) Firmness; (B) Scar; (C) Size; (D) Fruit Weight, and (E) Yield. Note the scale may vary along traits. played only a minor role. We also observed that the effect of each The final sets of hyperparameters for strawberry and hyperparameter depends on the layer, e.g. regularization or dropout blueberry phenotypes are indicated in Tables S1 and S2, were more important in first than in deep layers. respectively. Overall, our study shows that shallow In Figure 5A, the “trait” effect was excluded since it cannot be architectures are more competitive than deep architectures in controlled by the experimenter, although it was by far the most terms of PA, since the majority of models only included one influential variable. This is illustrated in Figure 5B, which shows the CNN layer. The number of filters -in combination with dropout- distribution of accuracies for each trait studied. Not only maximum has a large effect in the PA, but is highly dependent of the trait. accuracies varied across trait, the profiles were also extremely For instance, all optimal architectures for strawberry contain 128 different, usually multimodal. This suggests interactions between convolutions, but this is much more variable in the case of hyperparameter combinations, and it also indicates that trait– blueberry, with a range between 16 and 128 convolutional specific optimization should be performed whenever possible. operations. As for the fully connected layers, the situation is Figure 5 illustrates the kind of complex interactions that we less clear, and no obvious pattern is observed. We can highlight observed in hyperparameter optimization. For instance, Figures some characteristics though, for example, the number of hidden 5C, E, show the distinct influence of activation functions in fully connected layers is quite variable, but only a few neurons (4, percentage of culled fruit (Figure 5C) and brix (Figure 5E). 8, 12) are preferable in most of the architectures. As also reported Although “relu” activation function has been suggested as the by Waldmann (2018), combining weight decay and dropout activation of choice in recent DL literature (Maas et al., 2013; regularization is an efficient option to increase PA. Finally, the Pouladi et al., 2016), here we observed that linear or even best overlapping stride was 1 and optimum window size was 3 in sigmoid-like hyperbolic tangent (tanh) seemed to be a safer the convolutional layer, confirming Bellot et al. (2018) results. choice overall. It is relevant to note that interactions were clearly observed for some hyperparameters, such as the number Comparing Deep Learning With Bayesian of filters. For CulsTPer, either 16 or 128 filters resulted in optimum Penalized Linear Models accuracies, although they were also associated with the worst Figure 6A shows observed predictive abilities for each of the five hyperparameter combinations. In contrast, either 32 or 64 filters GP methods compared: BL, BRR, BRR-GM, RKHS, and CNNs in are to be preferred for average weight in strawberry (Figure 5F). strawberry. When averaged over traits in the strawberry species, Frontiers in Plant Science | www.frontiersin.org 7 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids FIGURE 5 | Hyperparameter influence on predictive accuracy in strawberry. Accuracy is defined as correlation between observed and predicted phenotypes by internal cross-validation. (A) Hyperparameter importance obtained from a random forest algorithm. nFilter: number of filters in the convolutional layer, activation_2, activation function in layer 2; reg_i, regularization in i-th layer; dropout_i, dropout rate in i-th layer; lr, learning rate; hidden_neurons, number of neurons in inter- mediate layers; hidden_layers, number of intermediate layers. (B) Distribution of accuracies along hyperparameter combinations for each phenotype. (C) Accuracies as a function of activation function for percentage of culls. (D) Accuracies as a function of number of filters for percentage of culls. (E) Distribution of accuracies as a function of activation for brix. (F) Distribution of accuracies as a function of number of filters for average fruit weight (AveWtT). PAs were 0.43, 0.43, 0.44, 0.44, and 0.44 for each of the five important differences between methods except in percentage of methods, respectively. By trait, the BRR-GM was best in AveWtT culled fruit. For this trait, CNN was ~20% better than any linear prediction, BL, BRR, and RKHS for MktWtE, RKHS, and BRR- model method. Interestingly, this trait was also the one with the GM for MkWtT, whereas CNN performed best in brix and largest epistatic component and exhibited a modest additive percentage of culled fruit. In all, nevertheless, there were no component (Figure 3E). Frontiers in Plant Science | www.frontiersin.org 8 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids FIGURE 6 | Predictive ability (PA) measured as correlation between observed and predicted phenotypes in the validation dataset in strawberry (A) and blueberry (B). Bayesian linear models (lasso and BRR) PAs in blueberry were computed with tetraploid genotypes, but were almost identical to those obtained with the diploidized ones. TABLE 2 | Posterior distribution means of variance component estimates (h ^2 ) As for the blueberry phenotypic traits, we found no and predictive ability (in simulated data using Bayes Ridge Regression (BRR), differences between GP methods BL and BRR (average PA = general model BRR (BRR-GM), Reproducing Kernel Hilbert Space regression 0.42), whereas CNNs were somewhat underperforming (average (RKHS), and Convolutional Neural Networks (CNN). PA = 0.40). The most remarkable result in blueberry is that CNN Replicate Architecture Genetic parame- Predictive ability (PA) performance was barely affected by the ploidy level employed to ter estimates build the genetic relationship matrix. In fact, the “diploid” option seemed more robust than the tetraploid one, except in fruit yield, ^ 2a h ^ 2d h ^ 2e h BRR BRR-GM RKHS CNN the only trait that was measured using a rating scale. 1 Additive 0.29 0.16 0.06 0.57 0.60 0.57 0.59 2 Additive 0.16 0.21 0.06 0.35 0.43 0.35 0.32 Simulation Study 3 Additive 0.26 0.25 0.05 0.52 0.58 0.51 0.51 Table 2 presents the main simulation results and Table S3, the 4 Additive 0.24 0.23 0.06 0.42 0.52 0.43 0.40 5 Additive 0.35 0.11 0.05 0.42 0.47 0.43 0.38 CNN architectures used for computing the PA in each replicate. 6 Mixed 0.17 0.19 0.06 0.33 0.44 0.33 0.30 Some interesting remarks can be made from these simulations. 7 Mixed 0.10 0.11 0.08 0.24 0.26 0.22 0.24 First, although biased, the variance component estimates do 8 Mixed 0.16 0.10 0.08 0.29 0.33 0.31 0.28 detect whether epistasis is important: h2e estimates are larger 9 Mixed 0.13 0.16 0.06 0.26 0.30 0.26 0.25 than the narrow-sense heritability in the presence of complete 10 Mixed 0.22 0.07 0.07 0.40 0.42 0.40 0.43 11 Epistatic 0.11 0.11 0.21 0.23 0.29 0.24 0.25 epistasis. Results are far less clear when only a fraction of loci 12 Epistatic 0.11 0.11 0.33 0.31 0.37 0.32 0.34 show epistasis. But the most relevant result is that, as we 13 Epistatic 0.12 0.09 0.23 0.34 0.38 0.35 0.32 hypothesized, predictive accuracies of CNN and additive 14 Epistatic 0.05 0.13 0.21 0.21 0.34 0.23 0.28 penalized methods were affected by genetic architecture. BRR 15 Epistatic 0.10 0.11 0.15 0.21 0.31 0.23 0.21 and RKHS were better than CNNs for the pure additive architecture, whereas the opposite was observed with pure DISCUSSION epistasis. However, BRR-GM, which accounts for non-linear relationships, was better than either CNNs or pure additive Supervised DL methods are examples of predictive modelling, linear models in most of the studied cases. consisting of approximating a mapping function (f) from input Frontiers in Plant Science | www.frontiersin.org 9 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids (X) to output (y) variables (Hornik et al., 1990). These problems (Ober et al., 2015; Gezan et al., 2017). This purpose can be include classification or regression tasks, to use the machine attained by different approaches. The simplest one is to learning jargon. Numerous successful applications of DL in incorporate a general matrix into the linear models made up of classification contexts have been published, e.g. pattern dummy variables. This model contains as many degrees of recognition (Drayer and Brox, 2014; Liang and Hu, 2015; Işin freedom as ploidy level per locus and allowing for any et al., 2016; Badrinarayanan et al., 2017) and natural language interaction structure between alleles (Enciso-Rodriguez et al., processing (NLP) (Deng and Liu, 2018). The DL implementation in 2018; Amadeu et al., 2019). RKHS models (Gianola et al., 2006; regression tasks is less abundant and the benefit of using these Gianola et al., 2008; de los Campos et al., 2009) are also able to methods remains uncertain (Bellot et al., 2018; Montesinos-López capture complex interaction patterns in a relatively et al., 2018a; Azodi et al., 2019). Most GP problems fall into the straightforward manner. Alternatively, a CNN can be regression task due to the complex nature of quantitative traits implemented using simply the raw data. Our analyses suggest (MacKay, 2009). So far, GP problems have been mainly addressed that DL can perform better than additive and RKHS models for using penalized linear models (De Los Campos et al., 2009; Crossa traits where the epistatic component is important and where et al., 2013a). More recently, the prediction of complex traits from narrow-sense heritability is low (e.g. percentage of culled fruit, genetic data is receiving attention from DL users (Ma et al., 2017; Figure 6). The simulation study performed in this work Montesinos-López et al., 2018a; Montesinos-López et al., 2019b). (Table 2) suggested that BRR including additivity, dominance The present work aim was to study the strengths and weaknesses of and the general dummy matrix described above can improve applying CNN to GP problems in polyploid species. CNN networks upon CNNs when the non-additive component is important, are attractive for addressing these problems, as they can although CNNs were better than strict additive linear models. accommodate situations where input variables are distributed Additional analyses with a wider range of phenotypic traits, along a space pattern, as with the case of SNPs. genetic structures and in larger datasets are needed to validate Implementing GS in polyploids is challenging. In allopolyploids, our results. genetic analyses have been traditionally implemented assuming An underlying goal of our work was to investigate the effect of diploidy, taking advantage of the fact that systems present disomic accounting for allele dosage in a GP context. Owing to the inheritance. High predictive performances have been observed in a complex nature of polyploids, genotype calling can be a challenge variety of allopolyploid species (e.g. cotton, strawberry, wheat) and and “diploidization”, i.e., considering a polyploid genome as traits (Gezan et al., 2017; Gapare et al., 2018; Juliana et al., 2019). diploid is usual (Bourke et al., 2018). Some studies have recently Recently, the importance of accounting for the contribution of investigated the effect of accounting the ploidy level in prediction subgenomes— potentially expressing epistatic effects— was accuracy in polyploids (Endelman et al., 2018; Nyine et al., 2018; considered in wheat, which shed light on the importance of Amadeu et al., 2019; Lara et al., 2019; Zingaretti et al., 2019). As accounting for this source of variation within the GP models in these previous results (de Bem Oliveira et al., 2019; Zingaretti (Santantonio et al., 2019). However, the scenario is even more et al., 2019), here we found that “diploidization”, in which all complex in autopolyploid species. Even with the recent advances in heterozygous genotypes are pooled, is as efficient and accurate as genotyping and sequencing technologies, the amount of genomic polyploid genotyping for prediction purposes, albeit it is trait information, and understanding, in most autopolyploid species is still dependent. Therefore, we conjecture that genomic selection, limited when compared to allopolyploid crops. One of the challenges particularly for low levels of ploidy, can pay off in polyploids is resolving the allelic dosage of individual locus (Bourke et al., 2018; even with simplified genotyping strategies. We need to be careful Gerard et al., 2018). From a quantitative genetics standpoint, we though as this approach may not be equally appropriate for all emphasize that polyploid species might present higher degrees of levels of ploidy and heterozygosity. For instance, this might be an complete and partial intra-locus interactions than diploids (Gallais, issue with sugarcane (with ploidy starting from 2n=20) as most 2003; Ferrão et al., 2018). Here, the interest of investigating DL individuals will be heterozygous. methods in polyploids is to take advantage of its non-linearity and less It is traditionally thought that DL requires extremely large restrictive assumptions for GP in comparison to the traditional linear datasets to be trained effectively (Alipanahi et al., 2015; Liang and model-based methods. Hu, 2015; Xiong et al., 2015). However, this and related works Previous studies (Ma et al., 2017; Bellot et al., 2018; (Ma et al., 2017; Bellot et al., 2018; Montesinos-López et al., Montesinos-López et al., 2018a; Montesinos-López et al., 2018a; Montesinos-López et al., 2019b) have shown that DL 2019b) have not shown clear advantages of DL over linear performance in GP is comparable to those of linear methods. model GP, as conventional models were competitive in terms Furthermore, the largest dataset analyzed so far with DL for of PA, with their added benefit of being faster and with more prediction (~100k individuals) did not show a consistent biological interpretability. However, DL could be better suited to advantage of DL (Bellot et al., 2018). Therefore, it seems that is explore non-linear components than linear models, especially the trait what really influences the success of DL and it appears when genotypes can be transmitted integrally, as occurs with not so critical the size of the dataset. This does not preclude, of asexual propagation. Certainly, the weak performance of classical course, that a large N is needed to advance in our knowledge on additive models in the presence of non-additive variance (e.g. best GP strategies. In fact, a larger N can be especially Figure 6 for percentage of culled fruit) confirms the relevance of recommended in clonally propagated species. It is well known developing methodologies that can incorporate non-linearity that an efficient breeding program tests a low number of crosses Frontiers in Plant Science | www.frontiersin.org 10 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids with a high number of genotypes in each of them. A cross would measure plant traits under say disease progression (Bergsträsser need to be tested if not much information is available though. et al., 2015), infrared thermography is able to scan temperature Numerous clonally propagated species of agricultural interest are and transpiration; NMR (nuclear magnetic resonance spectroscopy) polyploids, leading to high heterozygosity, non-linear and mass spectrometry (MS) are applied in plants metabolite interactions, and scarce prior knowledge about the crosses. In evaluation (Hong et al., 2016). These examples should be an ideal this scenario, as many cross-combinations as feasible should be scenario to neural networks as they involve imaging at high scale, produced to ensure the discovery and evaluation of the best complex, and heterogeneous datasets with multiple variables and genotypes (Grüneberg et al., 2009). The actual balance will outcome. In summary, we believe that the enormous amount of depend on the level of epistasis and dominance. If dominance data that can be automatically recorded revolutionizing plant is large, then the best clone would be within families with good breeding and the flexible nature of neural networks makes them performance; if dominance is low, this is not necessarily so. promising for meeting this future challenge. A drawback of DL models is that they lack biological or process interpretability and neither feature selection nor feature importance are obvious. In our opinion, GP algorithms are not DATA AVAILABILITY STATEMENT too useful for providing biological insight into the genetic basis of phenotypes; genome wide association studies should be more The strawberry dataset analyzed for this study can be obtained appropriate. In all, our results suggest that DL performance from V.M. Whitaker (vwhitaker@ufl.edu); the blueberry dataset improve as non-additive variance increases, a situation is usually can be obtained from P.R. Muñoz (p.munoz@ufl.edu) or Dryad encountered in fitness related traits. Digital Repository (accession number doi: 10.5061/ DL hyperparameter tuning is critical and difficult, especially in dryad.kd4jq6h) and https://gsajournals.figshare.com/articles/ terms of computational resources. Our analysis allows us to provide Supplemental_Material_for_de_Bem_Oliveira_et_al_2019/ some generic recommendations though. First, we and others (Bellot 7728365. et al., 2018; Montesinos-López et al., 2018a; Montesinos-López et al., 2019b) concluded that the predictive accuracy is mainly dependent of the trait, i.e., the architecture needs to be tuned for each trait AUTHOR CONTRIBUTIONS individually. Second, here we show that the popular relu activation MP-E conceived and supervised research. VW, PM, LO, SG, and function is not necessarily a universally valid activation function, LF contributed experimental data. LZ developed software and that interactions between hyperparameter combinations should be performed research. LZ and MP-E wrote the initial manuscript expected and that the number of convolutional filters and draft. All authors contributed to discussion and to writing the regularization in the first layers can have an important effect into final draft. the model performance (Figure 5). In general, we and other authors (Bellot et al., 2018; Waldmann, 2018) have reported that a shallow network is the best scenario in most cases. Nevertheless, DL can still FUNDING be attractive because it does not require feature engineering, a critical step in most machine learning methods. A further strength LZ was supported by a PhD grant from the Ministry of Economy of DL is its flexibility, e.g., it allows to define latent variables by using and Science (MINECO, Spain), by the MINECO grant autoencoder or embedding as a generative latent variable model. In AGL2016-78709-R to MP-E and from the EU through the addition, networks, even if shallow, can model complex BFU2016-77236-P (MINECO/AEI/FEDER, EU) and the relationships employing any non-linear activation function. “Centro de Excelencia Severo Ochoa 2016-2019” award SEV- Overall, there is no evidence that applying DL in GP 2015-0533. VW and LO were supported by the US Department applications necessarily improves the prediction accuracy upon of Agriculture/National Institute of Food and Agriculture that of classical linear model methods. PA depends on the trait Specialty Crop Research Initiative (SCRI) project ‘RosBREED: and is affected by many factors; no one algorithm is uniformly Combining disease resistance with horticultural quality in new better for all species and traits (Hu and Greene, 2018; Pérez- rosaceous cultivars' under Award Number 2014-51181-22378. Enciso and Zingaretti, 2019). PA usually decays if heritability is low or in the presence of high epistatic effects. Even under these ACKNOWLEDGMENTS conditions though, Bayesian models were better than CNNs in almost all cases (Tables S1, S2, Table 2). Even if performance of The VW lab acknowledges Dr. Sujeet Verma for curation of DL for GP is not outstanding, we cannot ignore that plant strawberry SNP data. breeding is based on both genotyping and phenotyping, and that high throughput phenotyping is critical for genomic dissection of complex traits (Cobb et al., 2013). Imaging and SUPPLEMENTARY MATERIAL computer vision can be employed to measure the physiological, growth, development, and other phenotypic properties of plants The Supplementary Material for this article can be found online with the advantage of being fast, non-invasive and a low-cost at: https://www.frontiersin.org/articles/10.3389/fpls.2020.00025/ strategy (Fahlgren et al., 2015), hyperspectral imaging is useful to full#supplementary-material Frontiers in Plant Science | www.frontiersin.org 11 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids REFERENCES Crossa, J. J., Beyene, Y., Kassa, S., Pérez, P., Hickey, J. M., Chen, C., et al. (2013a). Genomic prediction in maize breeding populations with genotyping-by- Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., et al. (2015). sequencing. G3. (Bethesda) 3, 1903–1926. doi: 10.1534/g3.113.008227 TensorFlow: large-scale machine learning on heterogeneous systems. Software Crossa, J., Pérez, P., Hickey, J., Burgueño, J., Ornella, L., Ceró N-Rojas, J., et al. available tensorflow. Org. 1, 1–19. doi: 10.1016/0076-6879(83)01039-3 (2013b). Genomic prediction in CIMMYT maize and wheat breeding Alipanahi, B., Delong, A., Weirauch, M. T., and Frey, B. J. (2015). Predicting the programs. Heredity (Edinb). 112, 48–60. doi: 10.1038/hdy.2013.16 sequence specificities of DNA- and RNA-binding proteins by deep learning. Daetwyler, H. D., Calus, M. P. L., Pong-Wong, R., de los Campos, G., and Hickey, Nat. Biotechnol. 33, 831–838. doi: 10.1038/nbt.3300 J. M. (2013). Genomic prediction in animals and plants: simulation of data, Amadeu, R. R., Cellon, C., Olmstead, J. W., Garcia, A. A. F., Resende, M. F. R., and validation, reporting, and benchmarking. Genetics 193, 347–365. doi: 10.1534/ Muñoz, P. R. (2016). AGHmatrix: R package to construct relationship matrices genetics.112.147983 for autotetraploid and diploid species: a blueberry example. Plant Genome 9. de Bem Oliveira, I., Resende, M. F. R., Ferrão, L. F. V., Amadeu, R. R., Endelman, J. B., doi: 10.3835/plantgenome2016.01.0009 Kirst, M., et al. (2019). Genomic prediction of autotetraploids; influence of Amadeu, R. R., Ferrão, L. F. V., Oliveira, I., de, B., Benevenuto, J., Endelman, J. B., relationship matrices, allele dosage, and continuous genotyping calls in et al. (2019). Impact of dominance effects on autotetraploid genomic phenotype prediction. G3. Genes Genomes Genet. 9, g3.400059.2019. doi: prediction. Crop Sci. 0, 0. doi: 10.2135/cropsci2019.02.0138 10.1534/g3.119.400059 Autonomio Talos [Computer software]. (2019). Retrieved from http://github.com/ de los Campos, G., Gianola, D., and Rosa, G. J. M. (2009). Reproducing kernel autonomio/talos. Hilbert spaces regression: a general framework for genetic evaluation1. J. Anim. Azodi, C. B., McCarren, A., Roantree, M., Campos, G. de los, and Shiu, S.-H. Sci. 87, 1883–1887. doi: 10.2527/jas.2008-1259 (2019). Benchmarking algorithms for genomic prediction of complex traits. De Los Campos, G., Naya, H., Gianola, D., Crossa, J., Legarra, A., Manfredi, E., bioRxiv 614479. doi: 10.1101/614479 et al. (2009). Predicting quantitative traits with regression models for dense Badrinarayanan, V., Kendall, A., and Cipolla, R. (2017). SegNet: A Deep Convolutional molecular markers and pedigree. Genetics 182, 375–385. doi: 10.1534/ encoder-decoder architecture for image segmentation. IEEE Trans. Pattern Anal. genetics.109.101501 Mach. Intell. 39, 2481–2495. doi: 10.1109/TPAMI.2016.2644615 de Los Campos, G., Gianola, D., Rosa, G. J., Weigel, K. A., and Crossa, J. (2010). Bassil, N. V., Davis, T. M., Zhang, H., Ficklin, S., Mittmann, M., Webster, T., et al. Semi-parametric genomic-enabled prediction of genetic values using (2015). Development and preliminary evaluation of a 90 K Axiom® SNP array reproducing kernel Hilbert spaces methods. Genet. Res. (Camb). 92, 295– for the allo-octoploid cultivated strawberry Fragaria × ananassa. BMC 308. doi: 10.1017/S0016672310000285 Genomics 16, 155. doi: 10.1186/s12864-015-1310-1 de los Campos, G., Hickey, J. M., Pong-Wong, R., Daetwyler, H. D., and Calus, M. P. L. Bawa, V. S., and Kumar, V. (2019). Linearized sigmoidal activation: a novel activation (2013). Whole-genome regression and prediction methods applied to plant and function with tractable non-linear characteristics to boost representation capability. animal breeding. Genetics 193, 327–345. doi: 10.1534/genetics.112.143313 Expert Syst. Appl. 120, 346–356. doi: 10.1016/j.eswa.2018.11.042 Deng, L., and Liu, Y. (2018). Deep Learning in Natural Language. (Springer) doi: Bellot, P., de Los Campos, G., and Pérez-Enciso, M. (2018). Can deep learning 10.1007/978-981-10-5209-5_11 improve genomic prediction of complex human traits? Genetics 210, 809–819. Drayer, B., and Brox, T. (2014). “Training deformable object models for human doi: 10.1534/genetics.118.301298 detection based on alignment and clustering,” in Lecture Notes in Computer Science Benevenuto, J., Ferrão, L. F. V., Amadeu, R. R., and Munoz, P. (2019). How can a (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in high-quality genome assembly help plant breeders? Gigascience 8, 1–4. doi: Bioinformatics). (Cham, Springer), 406–420. doi: 10.1007/978-3-319-10602-1_27 10.1093/gigascience/giz068 Duangjit, J., Causse, M., and Sauvage, C. (2016). Efficiency of genomic selection for Bergsträsser, S., Fanourakis, D., Schmittgen, S., Cendrero-Mateo, M., Jansen, M., tomato fruit quality. Mol. Breed. 36, 29. doi: 10.1007/s11032-016-0453-3 Scharr, H., et al. (2015). HyperART: non-invasive quantification of leaf traits Enciso-Rodriguez, F., Douches, D., Lopez-Cruz, M., Coombs, J., and de los using hyperspectral absorption-reflectance-transmittance imaging. Plant Campos, G. (2018). Genomic selection for late blight and common scab Methods 11, 1. doi: 10.1186/s13007-015-0043-0 resistance in tetraploid potato (Solanum tuberosum). G3. Genes Genomes Bernardo, R. (2008). Molecular markers and selection for complex traits in plants: Learning Genet. 8, 2471–2481. doi: 10.1534/g3.118.200273 from the last 20 years. Crop Sci. 48, 1649–1664. doi: 10.2135/cropsci2008.03.0131 Endelman, J. B., Carley, C. A. S., Bethke, P. C., Coombs, J. J., Clough, M. E., da Bourke, P. M., Voorrips, R. E., Visser, R. G. F., and Maliepaard, C. (2018). Tools Silva, W. L., et al. (2018). Genetic variance partitioning and Genome-Wide for genetic studies in experimental populations of polyploids. Front. Plant Sci. prediction with allele dosage information in autotetraploid Potato. Genetics 9, 513. doi: 10.3389/fpls.2018.00513 209, 77–87. doi: 10.1534/genetics.118.300685 Castillo-Juárez, H., Campos-Montes, G. R., Caballero-Zamora, A., and Montaldo, Fahlgren, N., Gehan, M. A., and Baxter, I. (2015). Lights, camera, action: high- H. H. (2015). Genetic improvement of pacific white shrimp [Penaeus throughput plant phenotyping is ready for a close-up. Curr. Opin. Plant Biol. (Litopenaeus) vannamei]: perspectives for genomic selection. Front. Genet. 6, 24, 93–99. doi: 10.1016/j.pbi.2015.02.006 93. doi: 10.3389/fgene.2015.00093 Feng, J., and Lu, S. (2019). “Performance Analysis of Various Activation Functions Chan, M., Scarafoni, D., Duarte, R., Thornton, J., and Skelly, L. (2018). “Learning in Artificial Neural Networks,” in Journal of Physics: Conference Series. (IOP). network architectures of deep CNNs under resource constraints,” in IEEE doi: 10.1088/1742-6596/1237/2/022030 Computer Society Conference on Computer Vision and Pattern Recognition Ferrão, L. F. V., Benevenuto, J., Oliveira, I., de, B., Cellon, C., Olmstead, J., et al. Workshops. (IEEE), 1784–1791. doi: 10.1109/CVPRW.2018.00222 (2018). Insights into the genetic basis of blueberry fruit-related traits using Cho, M., and Hegde, C. (2019). “Reducing the Search Space for Hyperparameter diploid and polyploid models in a GWAS context. Front. Ecol. Evol. 6, 107. doi: Optimization Using Group Sparsity,” in ICASSP, IEEE International 10.3389/fevo.2018.00107 Conference on Acoustics, Speech and Signal Processing - Proceedings (New Gallais, A. (2003). Quantitative genetics and breeding methods in autopolyploids York City: Institute of Electrical and Electronics Engineers Inc.), 3627–3631. plants. (Quae). doi: 10.1109/ICASSP.2019.8682434 Gapare, W., Liu, S., Conaty, W., Zhu, Q. H., Gillespie, V., Llewellyn, D., et al. Chollet, F. (2015). Keras: deep learning library for theano and tensorflow 7, T1. (2018). Historical datasets support genomic selection models for the prediction Available at: url: https://keras.io/k. of cotton fiber quality phenotypes across multiple environments. G3. Genes Cobb, J. N., DeClerck, G., Greenberg, A., Clark, R., and McCouch, S. (2013). Next- Genomes Genet. 8, 1721–1732. doi: 10.1534/g3.118.200140 generation phenotyping: requirements and strategies for enhancing our Gerard, D., Ferrão, L. F. V., Garcia, A. A. F., and Stephens, M. (2018). Genotyping understanding of genotype–phenotype relationships and its relevance to crop polyploids from messy sequencing data. Genetics 210, 789–807. doi: 10.1534/ improvement. Theor. Appl. Genet. 126, 867–887. doi: 10.1007/s00122-013-2066-0 genetics.118.301468 Colle, M., Leisner, C. P., Wai, C. M., Ou, S., Bird, K. A., Wang, J., et al. (2019). Gezan, S. A., Osorio, L. F., Verma, S., and Whitaker, V. M. (2017). An Haplotype-phased genome and evolution of phytonutrient pathways of experimental validation of genomic selection in octoploid strawberry. Hortic. tetraploid blueberry. Gigascience 8, 3. doi: 10.1093/gigascience/giz012 Res. 4, 16070. doi: 10.1038/hortres.2016.70 Frontiers in Plant Science | www.frontiersin.org 12 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids Gianola, D., Fernando, R. L., and Stella, A. (2006). Genomic-assisted prediction of Meuwissen, T. H. E., Hayes, B. J., and Goddard, M. E. (2001). Prediction of total genetic value with semiparametric procedures. Genetics 173, 1761–1776. doi: genetic value using genome-wide dense marker maps. Genetics 157, 1819– 10.1534/genetics.105.049510 1829. Gianola, D., van Kaam, J. B. C. H. M., and Toro, M. A. (2008). Reproducing kernel Meuwissen, T., Hayes, B., and Goddard, M. (2013). Accelerating improvement of hilbert spaces regression methods for genomic assisted prediction of livestock with genomic selection. Annu. Rev. Anim. Biosci. 1, 221–237. doi: quantitative traits. Genetics 178, 2289–2303. doi: 10.1534/genetics.107.084285 10.1146/annurev-animal-031412-103705 Gianola, D. (2013). Priors in whole-genome regression: the bayesian alphabet Min, S., Lee, B., and Yoon, S. (2017). Deep learning in bioinformatics. Brief. returns. Genetics 194, 573–596. doi: 10.1534/genetics.113.151753 Bioinform. 18, 851–869. doi: 10.1093/bib/bbw068 González-Camacho, J. M., de los Campos, G., Pérez, P., Gianola, D., Cairns, J. E., Montesinos-López, A., Montesinos-López, O. A., Hernández-Suárez, C. M., Mahuku, G., et al. (2012). Genome-enabled prediction of genetic values using Gianola, D., and Crossa, J. (2018a). Multi-environment genomic prediction radial basis function neural networks. Theor. Appl. Genet. 125, 759–771. doi: of plant traits using deep learners with dense architecture. G3. Genes Genomes 10.1007/s00122-012-1868-9 Genet. 8 (12), 3813–3828. doi: 10.1534/g3.118.200740 González-Recio, O., Rosa, G. J. M., and Gianola, D. (2014). Machine learning Montesinos-López, O. A., Montesinos-López, A., Crossa, J., Gianola, D., methods and predictive ability metrics for genome-wide prediction of complex Hernández-Suárez, C. M., and Martín-Vallejo, J. (2018b). Multi-trait, multi- traits. Livest. Sci. 166, 217–231. doi: 10.1016/j.livsci.2014.05.036 environment deep learning modeling for genomic-enabled prediction of plant Goodfellow, I., Bengio, Y., and Courville, A. (2016). Deep Learning (Cambridge: traits. G3. Genes Genomes Genet. 8, 3829–3840. doi: 10.1534/g3.118.200728 MIT Press Cambridge). Montesinos-López, O. A., Martín-Vallejo, J., Crossa, J., Gianola, D., Hernández- Grüneberg, W., Mwanga, R., Andrade, M., and Espinoza, J. (2009). Selection Suárez, C. M., Montesinos-López, A., et al. (2019a). A benchmarking between methods. Part 5: breeding clonally propagated crops. Plant Breed. Farmer deep learning, support vector machine and Bayesian threshold best linear Particip. 275–322. Available at: http://www.cabdirect.org/abstracts/ unbiased prediction for predicting ordinal traits in plant breeding. G3. Genes 20103075062.html. Genomes Genet. 9, 601–618. doi: 10.1534/g3.118.200998 Hancock, J. F., Sjulin, T. M., and Lobos, G. A. (2008). “Strawberries,” in Temperate Montesinos-López, O. A., Martín-Vallejo, J., Crossa, J., Gianola, D., Hernández- fruit crop breeding, vol. 393–437. (Wallingford, UK: Springer, Dordrecht). doi: Suárez, C. M., Montesinos-López, A., et al. (2019b). New deep learning 10.1007/978-1-4020-6907-9 genomic-based prediction model for multiple traits with binary, ordinal, and Henderson, C. R. (1984). Best linear unbiased prediction of performance and continuous phenotypes. G3. Genes Genomes Genet. 9, 1545–1556. doi: 10.1534/ breeding value. Biometrics, 172–192. Available at: papers3://publication/uuid/ g3.119.300585 627506AA-ACB7-491A-B468-9A3B5C2A52EC. doi: 10.2527/jas1985.601111x Namin, S. T., Esmaeilzadeh, M., Najafi, M., Brown, T. B., and Borevitz, J. O. Hong, J., Yang, L., Zhang, D., and Shi, J. (2016). Plant metabolomics: An (2018). Deep phenotyping: deep learning for temporal phenotype/genotype indispensable system biology tool for plant science. Int. J. Mol. Sci. 17, 767. classification. Plant Methods 14, 66. doi: 10.1186/s13007-018-0333-4 doi: 10.3390/ijms17060767 Nyine, M., Uwimana, B., Blavet, N., Hřibová, E., Vanrespaille, H., Batte, M., et al. Hornik, K., Stinchcombe, M., and White, H. (1990). Universal approximation of (2018). Genomic prediction in a multiploid crop: genotype by environment an unknown mapping and its derivatives using multilayer feedforward interaction and allele dosage effects on predictive ability in Banana. Plant networks. Neural Networks 3, 551–560. doi: 10.1016/0893-6080(90)90005-6 Genome 11, 1–16. doi: 10.3835/plantgenome2017.10.0090 Hu, Q., and Greene, C. S. (2018). Parameter tuning is a key part of dimensionality Ober, U., Huang, W., Magwire, M., Schlather, M., Simianer, H., and Mackay, T. F. C. reduction via deep variational autoencoders for single cell RNA (2015). Accounting for genetic architecture improves sequence based genomic transcriptomics. bioRxiv 385534. doi: 10.1101/385534 prediction for a drosophila fitness trait. PloS One 10, e0126880. doi: 10.1371/ Işin, A., Direkoǧlu, C., and Şah, M. (2016). “Review of MRI-based Brain Tumor journal.pone.0126880 Image Segmentation Using Deep Learning Methods,” in Procedia Computer Osborn, T. C., Chris Pires, J., Birchler, J. A., Auger, D. L., Chen, Z. J., Lee, H. S., Science (Amsterdam: Elsevier), 317–324. doi: 10.1016/j.procs.2016.09.407 et al. (2003). Understanding mechanisms of novel gene expression in Juliana, P., Poland, J., Huerta-Espino, J., Shrestha, S., Crossa, J., Crespo-Herrera, L., et al. polyploids. Trends Genet. 19, 141–147. doi: 10.1016/S0168-9525(03)00015-5 (2019). Improving grain yield, stress resilience and quality of bread wheat using Pérez, P., and De Los Campos, G. (2014). Genome-wide regression and prediction large-scale genomics. Nat. Genet. 51, 1530–1539. doi: 10.1038/s41588-019-0496-6 with the BGLR statistical package. Genetics 198, 483–495. doi: 10.1534/ Lara, L. A. D. C., Santos, M. F., Jank, L., Chiari, L., De Vilela, M. M., Amadeu, R. R., genetics.114.164442 et al. (2019). Genomic selection with allele dosage in Panicum maximum Jacq. Pérez-Enciso,, and Zingaretti, (2019). A guide on deep learning for complex trait G3. Genes Genomes Genet. 9, 2463–2475. doi: 10.1534/g3.118.200986 genomic prediction. Genes (Basel). 10, 553. doi: 10.3390/genes10070553 Le, N. Q. K., Huynh, T. T., Yapp, E. K. Y., and Yeh, H. Y. (2019). Identification of Pattanayak, S. (2017). “Unsupervised Learning with Restricted Boltzmann clathrin proteins by incorporating hyperparameter optimization in deep Machines and Auto-encoders,” in Pro Deep Learning with TensorFlow learning and PSSM profiles. Comput. Methods Prog. Biomed. 177, 81–88. (Berkeley, CA: Apress), 279–343. doi: 10.1007/978-1-4842-3096-1_5 doi: 10.1016/j.cmpb.2019.05.016 Pouladi, F., Salehinejad, H., and Gilani, A. M. (2016). Deep recurrent neural LeCun, Y., Bengio, Y., and Hinton, G. (2015). Deep learning. Nature 521, 436–444. networks for sequential phenotype prediction in genomics. arXiv Prepr. arXiv. doi: 10.1038/nature14539 1511.02554. Available at: https://arxiv.org/pdf/1511.02554.pdf [Accessed Liang, M., and Hu, X. (2015). “Recurrent convolutional neural network for object March 4, 2019]. recognition,” in Proceedings of the IEEE Computer Society Conference on Rajaraman, S., Jaeger, S., and Antani, S. K. (2019). Performance evaluation of deep Computer Vision and Pattern Recognition. (IEEE) 3367–3375. doi: 10.1109/ neural ensembles toward malaria parasite detection in thin-blood smear CVPR.2015.7298958 images. PeerJ 7, e6977. doi: 10.7717/peerj.6977 Liaw, A., and Wiener, M. (2002). Classification and regression by random forest. R Santantonio, N., Jannink, J. L., and Sorrells, M. (2019). A low resolution epistasis News 2, 18–22. mapping approach to identify chromosome arm interactions in allohexaploid Ma, W., Qiu, Z., Song, J., Cheng, Q., and Ma, C. (2017). DeepGS: predicting wheat. G3. Genes Genomes Genet. 9, 675–684. doi: 10.1534/g3.118.200646 phenotypes from genotypes using deep learning. bioRxiv 241414. doi: 10.1101/ Slater, A. T., Cogan, N. O. I., Forster, J. W., Hayes, B. J., and Daetwyler, H. D. 241414 (2016). Improving genetic gain with genomic selection in autotetraploid Maas, A. L., Hannun, A. Y., and Ng, A. Y. (2013). “Rectifier nonlinearities improve potato. Plant Genome 9, 1–15. doi: 10.3835/plantgenome2016.02.0021 neural network acoustic models,” in in ICML Workshop on Deep Learning for Tibshirani, R. (1996). Regression shrinkage and selection via the Lasso. J. R. Stat. Audio, Speech and Language Processing. Available at: https://www. Soc. Ser. B 58, 267–288. doi: 10.1111/j.2517-6161.1996.tb02080.x semanticscholar.org/paper/Rectifier-Nonlinearities-Improve-Neural- VanRaden, P. M. (2008). Efficient methods to compute genomic predictions. J. Dairy Network-Maas/367f2c63a6f6a10b3b64b8729d601e69337ee3cc [Accessed Sci. 91, 4414–4423. doi: 10.3168/jds.2007-0980 October 11, 2019]. Vitezica, Z. G., Varona, L., and Legarra, A. (2013). On the additive and dominant MacKay, T. F. (2009). Q & A: genetic analysis of quantitative traits. J. Biol. 8, 23. variance and covariance of individuals within the genomic selection scope. doi: 10.1186/jbiol133 Genetics 195, 1223–1230. doi: 10.1534/genetics.113.155176 Frontiers in Plant Science | www.frontiersin.org 13 February 2020 | Volume 11 | Article 25 Zingaretti et al. Exploring Deep Learning in Polyploids Waldmann, P. (2018). Approximate bayesian neural networks in genomic Performance Computing Environments - MLHPC ‘15,” (New York, New prediction. Genet. Sel. Evol. 50, 70. doi: 10.1186/s12711-018-0439-1 York, USA: ACM Press), 1–5. doi: 10.1145/2834892.2834896 Widder, D. V. (1954). The convolution transform. Bull. Am. Math. Soc. 60, 444– Zingaretti, M. L., Monfort, A., and Pérez-Enciso, M. (2019). pSBVB: a versatile 456. doi: 10.1090/S0002-9904-1954-09828-2 simulation tool to evaluate genomic selection in polyploid species. G3. Wiggans, G. R., Cole, J. B., Hubbard, S. M., and Sonstegard, T. S. (2017). Genomic (Bethesda) 9, 327–334. doi: 10.1534/g3.118.200942 Selection in dairy cattle: the USDA experience. Annu. Rev. Anim. Biosci. 5, 309–327. doi: 10.1146/annurev-animal-021815-111422 Conflict of Interest: The authors declare that the research was conducted in the Xiong, H. Y., Alipanahi, B., Lee, L. J., Bretschneider, H., Merico, D., Yuen, R. K. C., absence of any commercial or financial relationships that could be construed as a et al. (2015). RNA splicing. The human splicing code reveals new insights into potential conflict of interest. the genetic determinants of disease. Science 347, 1254806. doi: 10.1126/ science.1254806 Copyright © 2020 Zingaretti, Gezan, Ferrão, Osorio, Monfort, Muñoz, Whitaker and Yoo, Y. J. (2019). Hyperparameter optimization of deep neural network using Pérez-Enciso. This is an open-access article distributed under the terms of the Creative univariate dynamic encoding algorithm for searches. Knowledge-Based Syst. Commons Attribution License (CC BY). The use, distribution or reproduction in other 178, 74–83. doi: 10.1016/j.knosys.2019.04.019 forums is permitted, provided the original author(s) and the copyright owner(s) are Young, S. R., Rose, D. C., Karnowski, T. P., Lim, S.-H., and Patton, R. M. (2015). credited and that the original publication in this journal is cited, in accordance with “Optimizing deep learning hyper-parameters through an evolutionary accepted academic practice. No use, distribution or reproduction is permitted which algorithm. Proceedings of the Workshop on Machine Learning in High- does not comply with these terms. Frontiers in Plant Science | www.frontiersin.org 14 February 2020 | Volume 11 | Article 25