{'fraction_non_alphanumeric': 0.03456725299728713, 'fraction_numerical': 0.019908987485779295, 'mean_word_length': 4.515202702702703, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improving the Efficiency of Genomic Selection

We investigate two approaches to increase the efficiency of phenotypic prediction from genome-wide markers, which is a key step for genomic selection (GS) in plant and animal breeding. The first approach is feature selection based on Markov blankets, which provide a theoretically-sound framework for identifying non-informative markers. Fitting GS models using only the informative markers results in simpler models, which may allow cost savings from reduced genotyping. We show that this is accompanied by no loss, and possibly a small gain, in predictive power for four GS models: partial least squares (PLS), ridge regression, LASSO and elastic net. The second approach is the choice of kinship coefficients for genomic best linear unbiased prediction (GBLUP). We compare kinships based on different combinations of centring and scaling of marker genotypes, and a newly proposed kinship measure that adjusts for linkage disequilibrium (LD). We illustrate the use of both approaches and examine their performances using three real-world data sets from plant and animal genetics. We find that elastic net with feature selection and GBLUP using LD-adjusted kinships performed similarly well, and were the best-performing methods in our study.


Introduction
The ever-increasing amount of genetic information available in plant and animal breeding is reflected in the development of sophisticated models for the prediction of quantitative traits from genome-wide markers (Heffner et al., 2009, Hayes et al., 2009, also known as genomic selection (GS). The markers are typically dense single-nucleotide polymorphisms (SNPs). Approaches to this problem have moved from models with simple specifications, such as ridge regression (Hoerl and Kennard, 1970) and the LASSO (Tibshirani, 1996), to models based on highly-structured hierarchical distributions or semiparametric approaches. Some examples are the Bayesian alphabet models (Gianola et al., 2009, Bayesian models with complex priors as in Guan and Stephens (2011), models based on reproducing kernel Hilbert spaces (RKHS) such as Bravo et al. (2009), or the Bayesian LASSO (Park andCasella, 2008, Vazquez et al., 2012). This complexity is motivated by the need to correctly model the genetic architecture of the trait under investigation while producing models that are easy to estimate even for large SNP profiles. We focus on two key aspects of these models: the inclusion of a preliminary step that removes SNPs that appear to be redundant, and the choice of kinship matrices to model the relatedness of the genotyped individuals.
The former is equivalent to feature selection (Koller and Sahami, 1996), and can also be achieved by shrinking SNP effects towards zero, either through the use of constraints (Zou and Hastie, 2005) or through appropriate prior distributions in a Bayesian setting (Meuwissen et al., 2001). We examine the effectiveness in GS of Markov blankets (Pearl, 1988), which have been extensively studied in graphical modelling. They provide a principled solution to feature selection problems, and can be implemented as a data pre-processing step prior to fitting the GS model. We implement Markov blanket feature selection within four GS models applied to three real-world data sets covering barley, rice and mouse genetics.
Kinship matrices were traditionally derived from pedigrees using a single definition, but with kinships now being calculated from SNP data many different definitions are available (Astle and Balding, 2009). We investigate four kinship matrices within genetic best linear unbiased prediction (GBLUP). These include a novel matrix introduced by Speed et al. (2012a) which adjusts for the bias introduced by differences in local linkage disequilibrium (LD), and has been shown to increase the precision of heritability estimates.

Markov Blankets and Feature Selection
The Markov Blanket of a variable of interest T , denoted as B(T ), is the minimal set of variables conditioned on which all other variables in the model are probabilistically independent of the target T (Pearl, 1988). The Markov blanket of a phenotype y in a GS model is the minimal set B(y) ⊂ X such that that is, the subset of SNPs B(y) that makes all other SNPs redundant as far as the trait y is concerned. Given this property, knowledge of only the SNPs in B(y) is enough to determine the probability distribution of y. Other SNPs become superfluous, either because they are not associated with the trait or because their effect is mediated by the SNPs in B(y). If B(y) were known, any GS model could be fitted using B(y) instead of the full SNP profile X with no loss of information, but in practice the need to estimate B(y) means that some information loss is possible. This two-stage approach contrasts with models such as BayesB (Meuwissen et al., 2001) and the LASSO (Tibshirani, 1996), which select significant SNP effects concurrently with model fitting and in a model-specific way. Markov blankets can be efficiently estimated from data through the use of conditional independence tests, such as parametric and non-parametric tests for partial correlation (Legendre, 2000, Hotelling, 1953 or mutual information (Scutari and Brogini, 2012). Tests in common use do not require any tuning parameter except for the type I error threshold α. The estimated B(y) will satisfy (1) only approximately because of type I and type II errors. The former arise from the noisiness inherent to the data and limited sample sizes, while the latter are typical of weak dependencies which will often be omitted from the Markov blanket.
Several computationally-efficient heuristic algorithms for Markov blanket estimation are available in literature, including Grow-Shrink (GS; Margaritis, 2003), Incremental Association (IAMB; Tsamardinos et al., 2003) and Hiton-MB (Aliferis et al., 2010). For instance, IAMB can be used to estimate the Markov blanket of a trait y as follows: 2. Forward Phase: until no change is made, (a) test each SNP X i for independence from y conditional on the current Markov blanket B(y); (b) admit into B(y) the SNP whose test returned the lowest p-value if that p-value is smaller than α.
3. Backward Phase: for each X i ∈ B(y), remove X i from B(y) if y is independent of X i conditional on B(y) \ X i .
As a result, conditional independence tests are performed in order of increasing complexity, thus ensuring that in practice only a small number of SNPs is used for each test. Compared to single-SNP analyses, such as those described in Macciotta et al. (2009) andSchulz-Streeck et al. (2011), feature selection with Markov blankets is computationally more expensive because of the use of conditional (y ⊥ ⊥ X i | B(y) \ X i ) instead of marginal (y ⊥ ⊥ X i ) independence tests. However, as shown in Section 4, Figure 5, conditional tests are more effective at discarding SNPs that carry essentially the same information about the trait and select subsets with more predictive power for the same size.

Kinship Estimation
In the past, pedigree information was used to specify kinships, but such information is often missing or inadequate. SNP-based methods for measuring kinships have become increasingly common and have the advantage of measuring the realised amount of genome sharing, as opposed to the expected value provided by pedigree-based methods (Astle andBalding, 2009, Forni et al., 2011) The SNP-based kinship of two individuals is usually based on the average over SNPs of the product of their genotypes, coded as 0, 1 and 2 according to the count of one of the two alleles. By design, it can only capture the additive components of kinship, and it has very low power in identifying non-additive ones. In the following, we denote this genotype matrix with X, with rows corresponding to individuals and columns to SNPs, and with X i its ith column.
In human genetics, kinship is commonly measured as the proportion of shared alleles at each locus (Morris and Cardon, 2007). This approach is also known as identical-by-state (IBS) kinship, and will be denoted by K 0 . Unlike other kinship matrices below, K 0 is always non-negative. However, it cannot be expressed in the form XX T , which leads to parameters directly interpretable as SNP effect sizes (see Section 3 for details).
Another choice, common in plant and animal genetics, is to centre the genotypes (Habier et al., 2007, VanRaden, 2008 and estimate the kinship matrix as where m is the number of markers and p i is a vector with every entry equal to the population allele fraction, usually estimated as the mean of X i /2. Centring improves interpretability, since kinship values can be interpreted as an excess or deficiency of allele sharing compared with random allocation of alleles, and so zero can be interpreted as "unrelated". However, the requirement to estimate the p i , usually from the same data set, can cause problems in some settings (Astle and Balding, 2009). One criticism of both the above choices is that the sharing of a rare allele between two individuals counts the same as the sharing of a common allele. One natural approach to giving more weight to the sharing of a rare allele is to standardise over SNPs, thus obtaining The (i, j) entry of K 2 can be interpreted as an average over SNPs of the correlation coefficient estimated from a single pair of individuals, i and j (Astle and Balding, 2009). A modification of K 2 has been recently proposed by Speed et al. (2012a), based on evidence that the effects of SNPs are sensitive to uneven LD across the genome. In particular, SNP effects are over-estimated in high-LD regions and under-estimated in low-LD regions due the uneven tagging of causal variants. The contributions of causal variants are picked up by a larger number of SNPs in high-LD regions compared to low-LD regions, thus introducing bias in the GS models and in turn in subsequent inference such as prediction or heritability estimation. To correct for this bias, SNPs can be re-weighted: where the weight vector w = [w 1 · · · w m ] solves and C i is a vector of squared correlations of SNP i with neighbouring SNPs. SNP effects are set to decay exponentially with physical distance, according to a decay rate λ whose value reflects the average LD for the data set. As a result we have that where d i j is the distance between SNPs i and j, is approximately constant as the weights offset differences in LD as measured by the squared correlations C i j . For computational reasons, the minimisation in (5) is performed separately on different chromosomes and, within each chromosome, on different regions chosen based on λ .

Materials and Methods
We explored the effects of the approaches outlined in Section 2 on the predictive power of GS models using three publicly-available real-world data sets including continuous phenotypic traits. The yield data from the AGOUEB project  consist of 227 UK winter barley varieties and 810 SNPs. The heterogeneous mouse population  from the Wellcome Trust Case Control Consortium (WTCCC) consists of 1940 SNP profiles and 12545 SNPs; among the recorded traits, we consider growth rate and weight. The rice data set from Zhao et al. (2011) consists of 413 varieties of Oryza sativa with 73808 SNPs; among the 34 recorded traits, we consider the number of seeds per panicle because of its low variability among the various subpopulations included in the original analysis.
All data sets have been preprocessed by removing SNPs with minor allele frequencies < 1% and those with > 20% missing data. The missing data in the remaining SNPs have been imputed using the impute R package (Hastie et al., 2012). Other widely used imputation methods in genetics, such as that implemented in MaCH (Li et al., 2010), were not available because of the absence of accurate SNP maps; the position of many SNPs is unknown, and only genetic distances (in cM) were available between mapped SNPs. Furthermore, we removed one SNP from each pair whose allele counts have correlation > 0.90 to increase the numerical stability of the models.
• Partial least squares (PLS) regression as implemented in the pls R package (Mevik et al., 2011).
• Genetic BLUP (GBLUP) implemented in the synbreed R package (Wimmer et al., 2012). It uses the linear mixed model where g are the random effects and Z is a design matrix that can be used for example to indicate the same genotype exposed to different environments. Any positive definite matrix can be used for K. Fixed effects can also be included in (8) in order to capture purely environmental effects (Heffner et al., 2009).
When K can be expressed in the form XX T , GBLUP can be shown to be equivalent to the Bayesian linear regression in which K determines the transformation X * of the SNP genotypes. For instance, the original X i are used when K = K 1 ; the scaledX i from (3) when K = K 2 ; and the weighted w iXi / ∑ w i from (4) when K = K 3 . This formulation of GBLUP results in a more natural interpretation of SNP effects, and is sometimes known as random regression BLUP (RR-BLUP). An overview of its properties can be found in Piepho (2009) and Piepho et al. (2012).
Markov blanket feature selection has been performed with the IAMB algorithm as implemented in the bnlearn R package (Scutari, 2010), using the exact Student's t test for Pearson's correlation with a type I error threshold of α = 0.15. Each GS model was fitted both using all the available SNPs and using only the SNPs included in the Markov blanket.
The different kinship matrices were investigated within GBLUP, as the other GS models do not include an explicit kinship term. K 1 and K 2 were computed using synbreed. For K 3 , we used the freely available LDAK software (Speed et al., 2012b). The LD decay rate was set to λ = 50cM for the AGOUEB data, λ = 0.2cM for the mouse data and λ = 100cM for the rice data. Such values were found, through experimentation, to ensure the LD adjustment was effective without affecting the genetic information present in the SNP profiles. K 0 was computed with PLINK (Purcell et al., 2007). All configurations of GS models and kinships were fitted once using all SNPs available after preprocessing the data and once using only those in the Markov blanket.
The predictive power of the GS models was measured with Pearson's correlation coefficient ρ between the observed trait values and the predictions obtained from 10-fold cross-validation. For each model, cross-validation was run 5 times.
Markov blankets, kinship matrices and GS models were fitted separately for each fold in each cross-validation run, and the resulting correlations averaged. The correlation between observed and fitted trait values is also reported as a measure of goodness of fit. Table 1 reports the observed correlations (ρ, i.e. the correlation between the observed and the fitted trait values) and the predictive correlations (ρ CV , i.e. the correlations obtained from cross-validation) for PLS, ridge regression, LASSO and the elastic net. The corresponding correlations arising from the subset of SNPs included in the Markov blankets are labelled ρ MB and ρ CV,MB , respectively.

Results
First of all, we note that for α = 0.15 Markov blankets only select a small number of SNPs, regardless of the dimension of the SNP profile. The average size of the Markov blankets obtained from cross-validation is 185 for the AGOUEB data, 543 (for growth rate) and 525 (for weight) for the mouse data, and 293 for the rice data. Of those SNPs, 136 (74%) appear in at least half of the crossvalidation folds for AGOUEB, 241 (46%) for the mouse data and weight, 276 (51%) for the mouse data and growth rate, but only 15 (5%) for the rice data. This can be attributed to the very low ratio between sample size and number of SNPs in the rice data (< 0.01) compared to the mouse (0.15) and AGOUEB (0.28) data. As expected, the dimension reduction is smaller in the case of the AGOUEB data because of the limited number of available SNPs, despite the extensive LD present in cultivated UK barley , Rostoks et al., 2006. On the other hand, only a small proportion of the original SNPs are retained for the mouse and rice data sets (about 4% and 0.4%, respectively). In each case, the number of SNPs included in the Markov blankets is smaller than the sample size, thus ensuring the regularity and numerical stability of the GS models.
The position of mapped SNPs within the respective genomes is shown in Figure 1 (AGOUEB), Figure 2 (rice), Figure 3 (mice, weight) and Figure 4 (mice, growth). For all but the AGOUEB data, we can see how the Markov blankets arising from cross-validation identify some regions as associated with the trait being modelled (e.g. SNPs in the range [0.1cM, 75.8cM] of chromosome 1 are included with high probability for both traits in the mice data set) while completely discarding other regions (e.g. [96.3cM, 108.3cM] in chromosome 2 and [70.6cM, 88.5cM] in chromosome 3). The positions of these regions may provide useful prior information in subsequent association studies and in targeting future genotyping efforts. In the case of the AGOUEB data, marker density is not high enough to identify regions with markedly different association levels.  Table 1: Correlation coefficients for PLS, ridge regression, LASSO and the elastic net: ρ is the correlation between observed and fitted trait values; ρ CV is the predictive correlation obtained from cross-validation; ρ MB and ρ CV,MB are the corresponding quantities obtained using only the SNPs in the Markov blanket. ∆ 1 = ρ MB − ρ and ∆ 2 = ρ CV,MB − ρ CV . The highest value for each quantity and data set is shown in bold. The average dimension of the Markov blanket over cross-validation is reported in parentheses for each data set and trait.
CHR 1 CHR 2 CHR 3 CHR 4 CHR 5 CHR 6 CHR 7     Table 2: Correlation coefficients obtained in GBLUP using the four kinship matrices defined in Section 2.2. The highest value for each quantity and data set is shown in bold. ρ and ρ CV are defined as in Table 1.
We observe no loss in the predictive power of the GS models after the Markov blanket feature selection. In fact, the increased numerical stability resulting from the reduced number of SNPs slightly improved the predictive power of the GS models. The average of ρ CV over the four analyses was 0.481, 0.498, 0.471 and 0.521 for PLS, ridge, LASSO and elastic net respectively, while the corresponding averages for ρ CV,MB were 0.502, 0.509, 0.487 and 0.520, all with an approximate standard deviation of 0.0057 computed as in Hooper (1958).
If we choose α < 0.15, we obtain Markov blankets that are too small to capture polygenic effects (results not shown). A possible explanation for this behaviour may be that large values of α allow Markov blankets to initially include SNPs that are weakly associated with the trait, to the point that they would be individually discarded. In addition, among them there may be sets of SNPs that are jointly significant due to epistasis, and such sets are retained in the Markov blanket.
Furthermore, Markov blankets outperform other subsamples of the same size. To show this, we generated for each data set 100 random subsets of SNPs of the same size as the corresponding Markov blanket. In addition, we also generated subsamples including the most significant SNPs from a single-SNP analysis under cross-validation. The same t test as in Markov blanket estimation was used to assess significance. Subsequently, we used them to fit the GS models and to compute the predictive correlations corresponding ρ CV,MB . As we can see from Figure 5, the Markov blanket always results in higher values of ρ CV,MB .
The elastic net consistently outperforms the other GS models both with and without the use of Markov blankets, except for the AGOUEB data set (in which ρ CV is essentially the same for ridge regression and the elastic net).
Overall, from Table 2 we see that the predictive performance of GBLUP improves as the kinship matrices progress from K 1 through to K 3 . K 0 , while not   Table 3: Correlation coefficients for GBLUP using the four kinship matrices defined in Section 2.2 and Markov blanket feature selection. ρ CV,MB is defined as in Table 1; ρ MB,KIN is the predictive correlation obtained from cross-validation with the use of Markov blankets but with the kinship matrices estimated from the full SNP profile. The highest value for each quantity and data set is shown in bold. ∆ 1 = ρ CV − ρ CV,MB and ∆ 2 = ρ CV − ρ MB,KIN , using the ρ CV reported in Table 2. being competitive with K 3 , outperforms at least one of K 1 and K 2 for all data sets but AGOUEB. The means of the four ρ CV values are 0.504 for K 0 , 0.493 for K 1 , 0.501 for K 2 , and 0.518 for K 3 , all with an approximate standard deviation of 0.0057. Thus, GBLUP with K 3 performs as well as the elastic net and outperforms PLS, ridge regression and the LASSO. Although the elastic net performed equally well with or without Markov blanket feature selection, that is not the case for GBLUP (Table 3). For all kinship matrices, the reduced size of the Markov blanket relative to the full SNP set detracts from the computation of kinship coefficients, leading to a substantial loss of predictive power. If all SNPs are available and can be used to compute the kinship matrices, then much but not all of this loss is restored.

Conclusions
We have shown that Markov blanket feature selection applied as a preliminary step in GS with a continuous trait is able to greatly reduce the size of the SNP set with no loss (and possibly a small gain) in the predictive power of PLS, ridge regression, LASSO and the elastic net. Among those models, the elastic net was the best performer, followed by ridge regression. If GS is to be performed repeatedly for the same phenotype, for example in successive generations of crops, Markov blanket feature selection opens the possibility of reducing costs by genotyping many fewer markers.
In the absence of a feature selection step, the LD-adjusted kinship matrix K 3 (Speed et al., 2012a) provides slightly better predictive power than the matrix with no LD adjustment K 2 (Astle and Balding, 2009) and the IBS kinship matrix K 0 produced by PLINK (Purcell et al., 2007). In turn, K 2 and K 0 appear superior to the matrix with neither LD adjustment nor standardising of SNPs K 1 (Habier et al., 2007). Using K 3 , GBLUP was competitive with the elastic net (both had mean ρ CV = 0.52 over the four datasets).
Markov blanket feature selection is not compatible with GBLUP because of the requirement for large numbers of SNPs to compute the kinship matrix. However, Markov blanket feature selection has only a small adverse effect on GBLUP if all SNPs are available for computing the kinship matrix.