{'fraction_non_alphanumeric': 0.047316742847469104, 'fraction_numerical': 0.012252412392077196, 'mean_word_length': 4.345814479638009, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Introducing Gaussian covariance graph models in genome-wide prediction

Several statistical models used in genome-wide prediction assume independence of marker allele substitution effects, but it is known that these effects might be correlated. In statistics, graphical models have been identified as a useful tool for covariance estimation in high dimensional problems and it is an area that has recently experienced a great expansion. In Gaussian covariance graph models (GCovGM), the joint distribution of a set of random variables is assumed to be Gaussian and the pattern of zeros of the covariance matrix is encoded in terms of an undirected graph G. In this study, methods adapting the theory of GCovGM to genome-wide prediction were developed (Bayes GCov, Bayes GCov-KR and Bayes GCov-H). In simulated and real datasets, improvements in correlation between phenotypes and predicted breeding values and accuracies of predicted breeding values were found. Our models account for correlation of marker effects and permit to accommodate general structures as opposed to models proposed in previous studies which consider spatial correlation only. In addition, they allow incorporation of biological information in the prediction process through its use when constructing graph G, and their extension to the multiallelic loci case is straightforward.


Introduction
Most of the Bayesian and classical models used in genome-wide prediction (Meuwissen et al., 2001) assume that marker allele substitution effects follow independent Gaussian distributions which induces a diagonal covariance matrix; however, some biological phenomena point to nonindependent effects. On one hand, the existence of linkage disequilibrium (LD) may create a spatial correlation of marker effects (Gianola et al., 2003;Yang and Tempelman, 2012). On the other hand, the complex interactions between regions of the genome and interactions of gene products in the metabolism also suggest that the assumption of independent effects may not be tenable. Thus, accounting for correlated marker allele substitution effects may increase the predictive performance of statistical models used in genome-wide prediction. Although it has been known that marker effects might be correlated, the problem of accounting for such a correlation has not been widely studied. So far, there have been few studies investigating this interesting problem. Gianola et al., (2003) described a series of frequentist and Bayesian models accounting for within chromosome correlated marker effects. Yang and Tempelman (2012) proposed a Bayesian antedependence model considering a nonstationary correlation structure of SNP effects. The two studies only considered correlations among nearby markers.
Covariance estimation is recognized as a challenging problem in statistics (Stein, 1975), especially in high dimensional problems under the "big small " condition where the sample covariance matrix is not of full rank (Rajaratnam et al., 2008). As a consequence, high dimensional covariance estimation using graphical models is a contemporary topic in statistics and machine learning. Regularization methods imposing sparsity on estimators through structural zeros in the covariance or inverse covariance matrix have gain attention during recent years, (Carvalho et al., 2007;Letac and Massan, 2007;Rajaratnam et al., 2008). In these models, the pattern of zeros of the covariance (covariance graph models) or precision matrix (concentration graph models) is defined by means of an undirected graph . The nodes of this graph represent the underlying random variables, and when the joint distribution of these variables is multivariate Gaussian, pairs of nodes not sharing an edge in are either, marginally independent (Gaussian covariance graph models) or conditionally (given all other variables) independent (Gaussian concentration graph models). This paper focuses on Gaussian covariance graph models (GCovGM). In statistics, the usefulness of these models in the analysis of high dimensional data exhibiting dependencies is well known (Carvalho et al., 2007;Rajaratnam et al., 2008); consequently, given the need for flexible statistical methods to account for correlated marker effects in genome-wide prediction, the introduction of GCovGM in this area seems promising.
To our knowledge, this is the first study adapting the theory of GCovGM to account for correlated SNP allele substitution effects in genome-wide prediction. The theory of GCovGM has been developed to estimate the covariance matrix of an observable -dimensional random vector using iid observations. In contrast, in genome-wide prediction, the problem involves predicting marker effects, estimating residual variance(s), and estimating the covariance matrix of an unobservable random vector (SNP effects) using one -dimensional vector with phenotypic information along with genomic information. Thus, the objective of this study was to develop methods that adapt the theory of GCovGM to genome-wide prediction in order to account for correlated marker allele substitution effects.

Materials and methods
This section is split into the following subsections. Firstly, due to the fact that GCovGM theory is not widely known in the realm of quantitative genetics, a brief introduction and details on the challenge encountered when adapting it to genome-wide prediction are presented. Then, statistical methods adapting GCovGM to genome-wide prediction are described along with some approaches to build the graph . Finally, datasets used to implement our methods are described.

Gaussian Covariance Graph Models
Here, the case of a known graph is considered. By known it is meant that the pattern of zeros in the covariance matrix is actually known or that is defined on the basis of domain-specific knowledge. Some basic concepts in graph theory are provided in supporting information (Appendix A); the reader not familiar with this topic is encouraged to read it before reading the rest of the paper. Hereinafter, the operator |•| represents the determinant when the argument is a matrix and cardinality when the argument is a set. Let , , … , be a set of vectors in ℝ identically and independently distributed (0, Σ), the target is to estimate Σ. The graph determines the null entries of Σ as explained before. Formally, the parameter space is defined as follows. Let = ( , ) be an undirected graph with vertex set and edge set , then Σ lies in the cone ℙ = { : ∈ ℙ + = 0 ℎ ( , ) ∉ }, where ℙ + is the space of positive definite matrices.
Maximum likelihood estimation is possible only when > and because of the constraints that it imposes when adapting GCovGM in genome-wide prediction (see supplementary material, Appendix B) this paper focuses on Bayesian approaches only.

Bayesian estimation
For natural exponential families (as in concentration graph models) a class of conjugate priors corresponding to the Diaconis-Ylvisaker prior (Diaconis and Yilvisaker, 1979) is frequently used. However, covariance graph models correspond to curved exponential families instead of natural exponential families. It is easily checked because (Σ) ∝ exp (− 2 (Σ −1 ) − 2 |Σ|) , Σ ∈ ℙ , does not have the form of a natural exponential family. Silva and Ghahramani (2009) introduced the family of conjugate priors known as inverse -Wishart ( ( , )) whose pdf has the following form: , … , ). Under this prior: This corresponds to a (̃,̃) distribution, ̃≔ + ,̃≔ + . An important issue that has to be considered now is for which values of and , , (•) is a valid density. To find sufficient conditions the modified Cholesky decomposition of Σ, Σ = ′, where is a lower triangular matrix with diagonal entries equal to one and is a strictly positive diagonal matrix, is used. Then, we have the following bijection {Σ } ≥ ,( , )∈ → ({ } > ,( , )∈ , ), and this transformation induces the density ) where = |{ : > , ( , ) ∈ }| ∀ = 1,2, … , − 1. From the mathematical point of view, the problem is to find sufficient conditions for the following integral to be finite: after some manipulations, it can be shown that these conditions are the following (Khare and Rajaratnam, 2011). 1) ∈ ℙ + , 2) + 2 > + 2 ∀ = 2,3, … , , where = |{ < : ( , ) ∈ }|. For covariance graph models, there is a block Gibbs sampler algorithm to draw samples from the posterior. This sampler is based on partitioning the covariance matrix as: and (•,•) denotes the Inverse Gamma (•,•) distribution. Using this result and permutations, the partition can be done for the random variables in every step. Hence, this is not a standard Gibbs sampler because partitions change in every step; however, convergence can be established using results from Asmussen and Glynn (2011).

The Khare-Rajaratnam family of flexible priors for decomposable graphs
When is decomposable and its vertices are ordered according to a perfect elimination scheme (Khare and Rajaratnam, 2012), there exists a wider family of more flexible priors developed by Khare and Rajaratnam (2011). The parameter of the ( , ) family is common for all ; however, for decomposable graphs a more flexible prior with pdf of the form ̅ , (Σ) ∝ ) , Σ ∈ ℙ , ∈ ℙ + , = ( 1 , 2 , … , ) can be used. In this prior density, every has its own shape parameter . The price paid for this extra flexibility is that the graph has to be decomposable. When considering the modified Cholesky decomposition of the covariance matrix, the density in terms of and is: Sufficient conditions for this to be a proper density are: ∈ ℙ + , > 2 + + 2 (Khare and Rajaratnam, 2011). This prior is conjugate because the posterior density, given by: is a ̅̃,̃( , ) density, where ̃ is as defined before and ̃× 1 = {̃} = + − 2 . Hereinafter, this family of priors will be denoted as GWKR( , ). If in addition to be decomposable the graph is also homogeneous, direct sampling from the posterior can be performed (this case is discussed later), otherwise MCMC methods are used to draw samples from the posterior. Details of a block Gibbs sampler and the proof of its convergence can be found in Khare and Rajaratnam (2011). The full conditional distributions used in Khare and Rajaratnam's Gibbs sampler (Khare and Rajaratnam, 2011) are the following. Let = ( , ) be a decomposable graph with its vertices ordered according to a perfect elimination scheme, let

Covariance graph models for homogeneous graphs
For covariance graph models, certain properties of the graph = ( , ) have appealing mathematical consequences on the estimation problem. Covariance graph models take advantage of the fact that homogeneous graphs admit a Hasse ordering of their nodes (see Appendix A). The importance of having a graph with this sort of ordering is summarized in the following theorem (Khare and Rajaratnam, 2011): Let = ( , ), be a homogeneous graph with a Hasse ordering of its nodes. Then, Σ = ′ ∈ ℙ ⟺ ∈ ℒ ⟺ −1 ∈ ℒ , that is, matrices and −1 preserve the pattern of zeros in Σ. When is homogeneous, direct samples from the posterior can be obtained by reparametrization in terms of T = L −1 . Let : = { } < ,( , )∈ , then, it follows that the random

Adapting GCovGM to genome-wide prediction
The linear regression model considered here is the following: = + (1) where ∈ ℝ is an observable random vector containing response variables (e.g., corrected phenotypes or deregressed BV), ∈ ℝ is an unknown vector of marker allele substitution effects, ∈ ℝ is a vector of residuals, × is a matrix whose entries correspond to one to one mappings from the set of genotypes to a subset of the integers for every individual at every locus where is map corresponding to the genotype of the ℎ individual for the ℎ marker. The distributional assumptions are: | , , 2~( , 2 ), |Σ~(0, Σ) and | 2~( 0, 2 ). In standard GCovGM, it is assumed that > 1 independent and identically distributed multivariate Gaussian random variables are observed. On the other hand, in genome-wide prediction, most of the times phenotypic data corresponds to a single -dimensional vector and the target is estimating Σ instead of [ ] ≔ = Σ ′ + 2 , that is, the target is estimating the covariance matrix of an unobservable -dimensional random variable using a single -dimensional vector containing phenotypes and the genomic information contained in . Zhang et al. (2013) proposed methods to estimate covariance matrices corresponding to the sum of a low rank symmetric matrix and a sparse matrix, but these methods require a sample size larger than one and do not estimate Σ direclty. These are the difficulties encountered when considering the theory of GCovGM as a mean to model correlated marker effects in genome-wide prediction.

A hierarchical Bayes formulation
The flexibility of Hierarchical Bayesian modeling permits to cope with the problem of adapting GCovGM to genome wide prediction; it provides a simple and principled solution. The residual variance is given the following conjugate prior: 2~( 2 , 2 ). Regarding the covariance matrix of marker effects, under the conventional GCovGM, Khare and Rajaratnam (2011) provided recursive equations to find the posterior mean in closed form for homogeneous graphs. However, as explained in the previous section, the target here is estimating the covariance matrix of an unobservable random vector; therefore, even for this sort of graphs sampling from the joint posterior distribution is required. To this end, the following simple but useful property permits the use of a Gibbs sampler. Notice that under model 1 it follows that ( | ) = ( | ). This property, and the conjugacy of the priors considered here (IGW and GWKR), imply that the full conditional of Σ pertains to the same family of the prior. Therefore, because it is possible to obtain samples from these families and all other full conditionals are standard distributions, a Gibbs sampler can be implemented (Robert and Casella, 2010). Under the model termed Bayes GCov: Σ|~( , ) which can be used for general graphs. Then, the joint posterior has the following form: and |~( * , * ), * ≔ + ′ , * ≔ + 1. If is decomposable and the conditional prior for is a GWKR( , ) distribution, then this variation of the model is referred to as Bayes GCov-KR. In this case, the full conditional distribution of Σ is GWKR( * , * ), where * ×1 ≔ { * } = 1 + − 2 . Finally, under the conditional prior ( , ), if the graph is homogeneous, the model is denoted as Bayes Cov-H just to stress the fact that a homogenous graph is being considered and therefore the model is reparametrized in terms of the modified Cholesky decomposition of . In this case, the Gibbs sampler is more efficient due to the fact that direct samples from the full conditional distribution of can be drawn.

Some criteria to define
One of the first steps to carry out analyses with our models is defining the graph , that is, defining the marginal covariance structure of marker allele substitution effects. To this end, some approaches based on genetic criteria are presented in this section. The first one is based on the idea of spatial correlation (Gianola et al. 2003, Yang andTempelman 2012). Using a physical or linkage map, a window is defined based on a given map distance, or a given number of markers and it is slid across each chromosome. The order of markers is dictated by the physical or linkage map. This strategy induces a differentially banded or a banded covariance matrix. A second approach is based on the use of biological information. Using tools such as gene annotation information, markers can be clustered based on their function using approximations like the one presented in Peñagaricano et al. (2015). This will create groups or sets of loci and there are two options: permit correlations among effects of markers in different blocks or not. Finally, linkage disequilibrium between loci can be used. In this case, one of the metrics used to assess LD is chosen and those pairs of loci having a metric greater than a predefined threshold will be neighbors in .

Data analyses
One of the main issues related to GCovGM is their computational burden. Consequently, to ensure computational tractability, two small datasets were simulated in order to implement the proposed models. A single genome formed by 5 chromosomes of 10 cM length each, with 1605 biallleic markers and 1000 biallelic QTL was simulated. This genome was created via a forward-intime approach using software QMSim (Sargolzaei and Schenkel, 2013). Phenotypic records were created as the sum of the breeding value and an error term. For dataset 1, QTL effects were drawn from independent zero mean Gaussian distributions and were scaled such that the additive genetic variance was equal to 50. On the other hand, for dataset 2, QTL allele substitution effects were simultaneously drawn from a multivariate Gaussian distribution with null vector mean and a banded covariance matrix with bandwidth of size 10. These effects were then scaled in order to have an additive variance of 50. Residuals were drawn from independent Gaussian distributions with null mean and variance equal to 50, consequently, heritability was 0.5. Only markers with a minor allele frequency larger than 0.08 were considered in the analyses and ten replicates of each dataset were simulated. The graph was based on windows defined by a fixed number of marker loci (6), which induces a decomposable-non-homogeneous graph; therefore, models Bayes GCov-KR and Bayes GCov were fitted. Bayes A, a Bayesian model assuming uncorrelated effects, which is frequently used in genome-wide prediction, was also fitted. Training sets were formed by individuals from generations zero and one, and validation sets were comprised of individuals from generation 2.
Models were also fit to a real dataset. This dataset corresponded to a subset of the multibreed Angus-Brahman herd described in Elzo et al. (2012) and Elzo et al. (2013), containing 102 animals genotyped with the Illumina GoldenGate Bovine3K BeadChip (Illumina, Inc., 2011). Three traits were analyzed: daily feed intake (DFI, ĥ 2 =0.31), ultrasonic measure of percent of intra-muscular fat (UPIMF, ĥ 2 = 0.53) and body weight taken the same day than UPIMF (UW, ĥ 2 = 0.54) corrected for contemporary group (year-pen), age of dam, age of individual, sex, Brahman fraction and heterozygosity of the individual. The training set was defined as the oldest 72 individuals and the validation set as the remaining ones. The graph was defined using windows of 5 adjacent markers and only within chromosome correlations were allowed.
Predictive performance was assessed using the following criteria. Pearson correlation of phenotypes and predicted breeding values in the validation set (predictive ability) for real and simulated data, and the Pearson correlation between true and predicted breeding values (accuracy) in training and validation sets for simulated data. In each analysis, 15000 MCMC samples (first 5000 were considered burn in) were obtained using the Gibbs samplers described before. Analyses were performed using in-house R scripts (R Core Team, 2015).

Results
Tables 1 and 2 present the summary of the performance of the models fitted to simulated and real datasets respectively. According to all criteria (APA, AAT, AAV), our models clearly outperformed Bayes A in the two simulated datasets, the differences being more marked in the case of independent QTL effects (dataset 1). In these datasets, the flexibility of the GWKR priors yielded a better predictive performance. Also, the performance of our methods tended to be less variable; Bayes A showed a smaller variation only for APA in dataset 2. In the real data set, all predictive abilities were low. For DFI, Bayes A had the highest predictive ability, while for UPIMF and UW it had the lowest. For these two traits, Bayes GCov had the best predictive performance.

General comments about the models
In this study, the theory of GCovGM was adapted to genome-wide prediction through hierarchical Bayesian modeling. This development permits to account for correlated marker allele substitution effects in a flexible way. This flexibility is due to the ability of our models to accommodate covariance structures arising from biological considerations such as information from metabolic pathways and not only from the assumption of spatial correlation as has been done in previous studies (Gianola et al., 2003;Yang and Tempelman, 2012). Thus, covariances between effects of markers which are not physically linked are permitted. Furthermore, the possibility of defining the graph using tools such as gene annotation provides a way to incorporate biological information in the prediction process.
Several approaches to define the graph based on biological principles were presented. These approaches involve the assumption of spatial correlation and the aforementioned use of existing bioinformatics tools to create "functional" sets of SNP whose effects are correlated. In general, the second strategy would induce graphs with no special properties. However, due to the theoretical and numerical advantages of decomposable graphs discussed previously, it is convenient to work with this sort of graphs whenever possible. To this end, in a submitted paper, we have proven two propositions and a corollary that provide conditions on the edges set and the ordering of markers, such that the induced graph is decomposable. For the sake of completeness, these propositions and the corollary are presented in Appendix B. Proposition 1 in Appendix B is the most general, but when is defined using biological information and subsets of different "functional" SNP sets are allowed to be correlated, its conditions are more difficult to satisfy. On the other hand, proposition 2 in Appendix B and its corollary are more restrictive in terms of the covariance structure, but they provide easier ways to order markers and define the edge set, that guarantee decomposability. Once the "functional" sets have been defined, if these conditions are not satisfied, these theoretical results provide a basis to find a decomposable super-graph containing the original graph . Such a supergraph is known as the cover of (Khare and Rajaratnam, 2012).
In GCovGM, the family of homogeneous graphs is the one with more attractive properties. This is why the implementation of Bayes GCov-H is easier and faster because direct sampling of is feasible. However, finding this kind of graphs is, in general, not an easy problem. An example of a homogeneous graph is a rooted tree where all nodes are children of a single parent (the root). Thus, under the approach of using biological information to define the graph , a homogeneous graph can easily be found as follows: The tree structure mentioned above is imposed to each "functional" set and no correlations between effects of markers in different sets are allowed. It also holds when each "functional" set is assumed to be a complete. All the strategies mentioned before might appear restrictive, but notice that assuming independent marker effects amounts to imposing a covariance structure as well. In fact it is a special case of our approach when the edge set is the empty set.
Here, the focus was on Bayesian models because under the GCovGM framework, they can deal with the "big p small n" setting. However, in Appendix C, a frequentist approach to find the empirical BLUP of is presented. This formulation is based on the EM algorithm (Dempster et al., 1977) combined with GCovGM theory which permits obtaining estimators of dispersion parameters Σ and 2 which are used to build the mixed model equations corresponding to model 1 whose solution yields the empirical BLUP of (Henderson,1963). This formulation involves a partition of data induced by the assumption that different groups (e.g., half-sib families) have different sets of marker effects. Such an assumption has also been considered by other authors like Gianola et al. (2003).
Even with the aid of bioinformatics, biochemistry and physiology to construct the graph , it may not reflect the actual underlying covariance structure, but important correlations might be captured resulting in an improvement of the accuracy of genome-wide prediction. Covariance model selection involves finding the pattern of zeros and estimating the non-zero elements of either the precision or the covariance matrix (Bickel and Levina 2008;Khare et al., 2013). Model selection in GCovGM has not been as well studied as its counterpart in Gaussian concentration graph models. There exist some frequentist methods that induce sparsity based on penalized likelihood approaches (Bien and Tibshirani, 2011) and others based on the idea of inducing sparsity in the parameter of the modified Cholesky decomposition of (Rothman et al., 2010). From the Bayesian perspective, some methods based on the Bayesian lasso have been proposed, e.g., Wang (2012), but their main limitation is the computational burden. Some of these models could be implemented in genome-wide prediction following approaches similar to those presented in this study.

Extension to multiallelic loci
Here, biallelic loci were considered, but in some cases multiallelic loci have to be dealt with. In the future, models could be fit using genotypes for actual genes instead of molecular markers. In such a case, there could be more than two alleles per locus. A similar situation occurs when fitting effects of haplotypes built from two or more consecutive markers (Meuwiseen et al., 2001;Calus et al., 2008). The methods developed here can be easily extended to the multiallelic case. If there are alleles at locus , then the corresponding columns of the design matrix are formed by defining − 1 variables as follows: is the genotype of the ℎ individual for the ℎ allele at locus and " − " represents an allele different from . The graph can be built based on the ideas discussed before, with extra considerations at the intra-locus level. For example, it could be assumed that effects of alleles of the same locus are all correlated.

Data analyses
In general, Bayes GCov and Bayes GCov-KR outperformed Bayes A except for DFI. Differences between our models and Bayes A were more marked when QTL effects were independent and models considered SNP effects (dataset 1). In the ideal scenario where the model considers the causal variants (the QTL) instead of markers, the benefit of accounting for marginal correlation is smaller as suggested by a smaller difference in the three criteria used to assess predictive performance. This behavior may suggest that when considering the causal variants instead of proxies as the SNP, models assuming independent effects yield an acceptable predictive performance even when the true covariance matrix is non-diagonal. Hopefully, this ideal scenario where the causal variants determining a phenotype, or at least most of them, are known will be reached in the near future. The largest difference between the method Bayes A and our methods (15.7%) was observed for AAV in dataset 1, while the smallest one (0.5%) was observed for AAT in dataset 2, in both cases, when comparing it with Bayes GCov-KR. Although Bayes GCov-KR had higher APA, AVA and ATA values in the simulated datasets, notice that the differences compared to Bayes G-Cov were small, being slightly larger in dataset 2; therefore, in these simulations the gain in fitting a more complex model which considers as many shape parameters as markers did not yield a notorious gain in accuracy or predictive ability. The gains in accuracy in the validation set observed in dataset 1 are larger than those found by Yang and Tempelman (2012) when comparing their antedependence models with their independent marker effects counterparts Bayes A and Bayes B, whereas gains in accuracy observed in dataset 2 were comparable (they found increments in accuracy of breeding values in the testing population up to 3%). The simulated data in Yang and Tempelman (2012) were similar to dataset 1, where it is expected that correlation among SNP arises from physical proximity to the same causal variants. They also considered a heritability value of 0.5. As to the real data analyses, first of all it has to be considered that these analyses are just a proof of concept and analyses involving larger populations with denser marker panels are needed. We did not consider a larger population and a denser SNP panel due to limited computational resources. For DFI, Bayes A had a higher PA whereas for UPIMF and UW Bayes GCov had the best PA followed by Bayes GCov-KR. For a mice population, Yang and Tempelman (2012) also found that a model assuming independence (Bayes B) outperformed one of their models accounting for correlation (ante-Bayes A) in terms of predictive ability. Although these traits have medium to high estimated heritabilities (Elzo et al., 2012;Elzo et al., 2013) PA were very low. Perhaps this is due to the low number of SNP considered here (2407). The study of Gianola et al. (2003) did not consider data analysis.

Final remarks
This paper introduces the theory of GCovGM in the context of genome-wide prediction which permits to account for correlated marker effects in a very flexible way in terms of the marginal covariance structure. Models developed here also allow incorporating biological information in the prediction process through its use when building graph .
Perfect elimination ordering. An ordering of a graph = ( , ) is defined to be a perfect elimination ordering if a triplet { , , } with > > such that ( , ) ∉ and ( , ), ( , ) ∈ does not exist. Subgraph. The graph ′ = ( ′, ′) is said to be a subgraph of graph = ( , ) if ′ ⊆ and ′ ⊆ . Induced subgraph. Consider the graph = ( , ) and a subset ⊆ . Define = ( × ) ∩ . The subgraph = ( , ) is defined to be a subgraph of induced by . Decomposable graph. An undirected graph = ( , ) is a decomposable graph if it does not contain a cycle of length greater than or equal to four as an induced subgraph. It turns out that decomposable graphs are characterized by the existence of a perfect elimination ordering of their vertices; therefore, a graph = ( , ) is decomposable iff its vertices admit a perfect elimination ordering.

Connected graph. A graph
is said to be connected if any pair of distinct vertices in are connected, that is, there exists a path between them. Directed edges. An edge is said to be directed if ( , ) ∉ whenever ( , ) ∈ . If ( , ) is a directed edge then is said to be a parent of and is said to be a child of v. Directed graph. A graph = ( , ) such that its edges are directed is defined as a directed graph. Directed acyclic graph. A directed acyclic graph (DAG) is a directed graph with no cycles. Tree. A tree is a connected graph with no cycle of length greater or equal than 3.

Rooted tree.
A rooted tree is a tree in which a particular node is distinguished from the others and designated the root of the tree. This node is the ancestor of all other nodes in the tree. An ancestor of a node in a rooted tree with root node is any node in the path from to . The edge set is defined as follows. If ̅ ≠ ̅ , → , and ∄ such that → → then put a directed edge from to . Hasse perfect vertex elimination scheme or Hasse ordering. Once the Hasse diagram of has being built, the nodes of are ordered in the following way. The ordering is descending starting from the root of the tree; therefore, nodes pertaining to equivalence classes on the top of the Hasse diagram are assigned the largest levels. Within every equivalence class with more than one node, the ordering is arbitrary. Hence, the ordering is not unique. Any ordering that gives an ancestor a higher level than any of its descendants in the Hasse diagram of is defined to be a Hasse perfect vertex elimination scheme or simply a Hasse ordering of the nodes of .

Maximum likelihood estimation of for general graphs, standard problem
If the sample size is larger than , then maximum likelihood estimation of Σ is feasible. After removing constant terms from the negative log-likelihood the following is the objective function to be minimized: * (Σ) = (Σ −1 ) + |Σ|, Σ ∈ ℙ , where is the sample covariance matrix. Notice that the objective involves Σ −1 instead of Σ. This objective function is not convex, which makes this minimization more difficult than the minimization problem for concentration graph models. One important feature of covariance graph models is that they correspond to curved exponential families instead of the well-studied exponential families as is the case of concentration graph models (Khare and Rajaratnam 2011), it poses a more challenging problem.
An iterative conditional fitting (ICF) algorithm to minimize * (Σ) was developed by Chaudhuri et al. (2007); however, because we are dealing with a non-convex optimization problem, convergence to a global or even a local minimum is not guaranteed.
The algorithm is based on the following partition of Σ: where Σ 11 is the 1,1 entry of Σ, . is the first column of Σ without the first entry and Σ −1,−1 is the submatrix of Σ resulting from deleting its first row and column. Using the standard rules for inversion by partitioning: Notice that knowing Σ, we can get ( .1 , Σ −1,−1 , 1 ) and vice versa; consequently, we have a one to one transformation. By using permutations, the same partition can be performed for every one of the random variables represented in graph . The high level of the algorithm is the following: 1) Partition Σ as ( .1 , Σ −1,−1 , 1 ), 2) minimize * (Σ) with respect to .1 treating as fixed the current values of Σ −1,−1 and 1 and 3) minimize * (Σ) with respect to 1 fixing the current values of Σ .1 and Σ −1,−1 . The same is repeated for the variables and it corresponds to one iteration of the algorithm.
The minimization problem is solved by minimizing the following quadratic form with respect to 1 (Chaudhuri et al.,2007): where 1 : = (Σ 1 ) (1, )∈ , .1 and −1,−1 are elements obtained after partitioning as Σ was partitioned in ( . 1) and 1 is a matrix of zeros such that: . = 1 . This is a standard problem and its solution is ̂1 = ( 1 ′ Σ −1,−1  Kauermann (1996) proposed to modify the objective function in order to make it a function of Σ, which makes the problem convex. The new objective function has the form: ̃( Σ) = (Σ −1 ) + |Σ|, Σ ∈ ℙ . A trick based on using the maximal cliques of is applied to solve this problem and the solution is known as the Kauermann's dual estimator. Under certain conditions, convergence of the ICF algorithm at least to a local stationary point can be proved (Drton et al., 2006).
where . Is the ℎ row of .

Maximum likelihood estimation of in genome-wide prediction
Unlike the Bayesian approach, envisaging a frequentist solution to the problem of adapting GCGM to genome-wide prediction is not straightforward and there is not a direct and principled method to cope with this problem. Therefore, some ad hoc assumptions have to be done in order to provide a frequentist formulation. The method proposed here involves two steps. The first one combines the EM algorithm (Dempster et al., 1977) with GCovGM to estimate covariance components. The second one involves plugging this estimates into mixed model equations corresponding to model 1 in order to obtain the empirical BLUP of (Henderson, 1963). According to the rationale of the EM-algorithm, we define as the augmented or missing data, then we find the maximizers of the augmented likelihood as if g were observable and finally we compute their expected values with respect to the distribution of the missing or augmented data given the observed data. As mentioned in the manuscript, maximization with respect to Σ is only possible if > . In model 1, we have a single -dimensional vector and the target is to estimate the residual variance and the covariance matrix of the -dimensional vector : therefore, in terms of the standard problem: =1.Thus, an ad hoc solution is to assume that data can be split into > groups such that each group has a different vector of marker effects, that is, = + , ∀ = 1,2, … , . Currently, as more and more animals are genotyped, for SNP panels of moderate density (e.g., 50K) the case > can be found. For example, this is the case of the Holstein population in the US. However, this is not the most common situation and it is important to notice that it does not imply that > which is the necessary condition to carry put maximum likelihood estimation of Σ. One of the simplest ways to split a population into groups is by considering families (e.g., half-sibs or full-sibs) as in Gianola et al. (2003). Currently, the requirement > will be met by very few populations when considering a relatively small number of markers and this is the reason to not considering the frequentist approach in the manuscript. Notwithstanding, in this appendix we provide a way to carry out maximum likelihood estimation of the dispersion parameter ≔ (Σ, 2 ) in a genome-wide prediction model based on multiple linear regression which later permits to obtain the empirical BLUP of .
It is also assumed that: 1 , … , are iid (0, Σ −1 ), 1 , … , are iid (0, 2 ) and ( , ′ ) = 0, ∀ 1 ≤ , ′ ≤ , where is the number of observations in group ; therefore, ∑ =1 = . Under these assumptions, the complete log-likelihood can be written as: The expected values of sufficient statistics for the covariance parameters taken with respect to the conditional distribution of the missing data given the observed data have to be found. The sufficient statistic for is ( , * ′ * ), * = − * * . Also, given , 1 , … , are independent with the following distributions: |~( −1 ′ 2 , −1 ), where ≔ ′ 2 + Σ −1 . Similarly, it follows that * |~( 2 −1 , 2 ( − 2 −1 )), where = * ′ ⨂Σ * + . Hence, Applying the Woodbury's identity, [ | ] can be alternatively expressed as: where ≔ Σ ′ + 2 . It does not require inversion of Σ, it requires inverting where Σ ( +1) is computed using methods explained before. For homogeneous graphs, function ℎ(•) has closed form after reparametrizing the objective function in terms of ( , ) as shown previously in this section. Once the algorithm converges and the maximum likelihood estimates of Ω and 2 are obtained, these are plugged in the mixed model equations corresponding to model 1 to obtain the empirical BLUP of (Henderson, 1963): orderings is reduced. Thus, in order to provide a simpler way to order markers, the following proposition only requires the existence of an ordering of the blocks and a structure on the edges set satisfying certain conditions that permit to find a perfect elimination ordering of markers.
Then the following ordering strategy (denoted by ) of marker loci is a perfect elimination ordering: once blocks have been ordered according to , markers are ordered in such a way that the smaller the index of a block the smaller the indices of the markers pertaining to that block. The ordering inside each block is done as follows: markers in are given the largest indices in block . In addition, under this ordering strategy, condition 2.2 is also necessary for to be a perfect elimination ordering whereas condition 2.1 is not.

Corollary to Proposition 2
Consider the "super graph" formed by regarding the blocks as super nodes and ℒ as a "super vertices set". Then, under conditions 2.1 and 1.4, if the "super vertices set" admits a perfect elimination ordering, the ordering defined in proposition 2 corresponds to a perfect elimination scheme.

Proofs of propositions 1 and 2, and corollary to proposition 2
Proof of proposition 1 Approaches based on marker's position induce either a banded or a differentially banded matrix and these structures are known to correspond to decomposable graphs. Similarly, F.1 induces block diagonal matrix which also corresponds to a decomposable graph. The only non-trivial part is to prove the statements made about conditions for F.2 to induce a decomposable graph. To prove it, we show that it is possible to find a perfect elimination ordering, and we use the fact that the existence of a perfect elimination ordering characterizes decomposable graphs.
On the other hand, if one block is linked to the other two, and these are the only existent links, condition 1.2 guarantees that . 1 does not hold. Until here, only linked blocks were considered. Notice that if at least one index of the triplet { , , } corresponds to a marker in an isolated (i.e. not linked) block, then . 1 does not happen.
In conclusion, under the four conditions of result 1 it is possible to find a perfect elimination ordering and therefore the graph induced under approach 3.b is decomposable. Now we proceed to prove necessity of conditions 1,2 and 3. To this end, it is proven that under 3.b, if these conditions do not hold, then the induced graph is not decomposable. Proof of necessity of conditions 1.1, 1.2 and 1.3 Hereinafter when the words "there is always" precede a condition, it is meant that under any ordering of markers, the graph always satisfies the condition. Suppose that 1.1 does not hold. It means that:  There is always at least one triplet of linked blocks { , ′ , ′′ }, 1 ≤ , ′ , ′′ ≤ |Ψ|, linked in the following way: ′ ≠ ′′ , ′ ≠ ′ ′′ , ′′ ≠ ′′ ′ , the three pairwise intersections , ′ and ′′ are non-empty, and there exist at least one triplet { , , } such that > > and ∈ ′′ ∩ , ∈ ′ or ∈ ′′ , ∈ ′ ∩ ′ , and ∈ ′′ , > > , or  There is always at least one triplet of linked blocks { , ′ , ′′ } and a triplet of indices { , , } such that > > , exactly one or two of the three sets { , ′ , ′′ } are empty and pertains to a non-empty intersection set. It is easy to notice that in any of the two events described in the first case, condition A.1 is satisfied. On the other hand, in the second case, if pertains to the (non-empty) intersection set of a block, say then we can always by pick ∈ ′′ and ∈ ′ such that ( , ) ∉ , then, because ∈ it follows that ( , ), ( , ) ∈ , that is, . 1 holds. Therefore, if 1.1 does not hold, a perfect elimination ordering does not exist and the graph is not decomposable. Similarly, if there is always at least one triplet of blocks { , ′ , ′′ } satisfying the converse of condition 1.2 then . 1 always holds, that is, a perfect elimination ordering does not exist and consequently the graph is not decomposable. Finally, if 1.3 does not hold then always exists at least a linked pair { , ′ } such that ∃ ∈ ′ , ∈ ′ ∩ , ∈ ′ , > > which immediately implies . 1. Proof of non-necessity of condition 1. 4 We refute the statement that condition 1.4 is necessary by constructing a counter-example. Suppose that a graph is comprised of two blocks. Let be the subset of containing edges comprised of one node pertaining to block 1 and one node from block 2. Also assume that block 1 contains nodes { , , , }, 1 2 ≔ { , , }, 2 1 ≔ { , } and ≔ {( , ), ( , ), ( , ), ( , ), ( , )}. Consider any ordering satisfying the following: ( ) = 1, ( ) = 2, ( ) = 3, ( ) = 4. Consider > > , ∈ 2 1 , if = ( ) then ( , ) ∉ ⇒ ( , ) ∉ , if = ( ) then ( , ) ∉ ⇒ = 1 = ( ) ∴ ∄ < and consequently . 1 does not hold. We have found a family of perfect elimination orderings and its existence implies that the graph is decomposable even though condition 4 does not hold. It proves that condition 3 is not a necessary condition. ∎ Proof of proposition 2 Proof of sufficiency The general approach is very similar to that used to prove sufficiency of conditions stated in result 1. Let Λ be the set of blocks linked to at least one block. Assume that markers have been ordered as stated in result 2, this ordering is denoted by . Let be the set of indices corresponding to markers in block under . Let { , , } be three arbitrary indices such that according to , > > , and ∈ , ∈ Λ be the subset of indices corresponding to markers in block whose effects are partially correlated with the effects of a subset of markers in another block 1 ≤ ≤ |Λ|. In the following, it is proven that under conditions 1.4, 2.1 and 2.2 this within-block ordering strategy yields a perfect elimination ordering. Four cases depending on the position of indices , , are considered. Case 1) The triplet { , , } pertains to the same block. Because each block corresponds to a complete graph, it follows that in this case condition . 1 is not satisfied. Case 2) ∈ , { , } ∈ ′ , { , ′ } ∈ Λ , ( ) > ( ′ ). If ( , ) ∉ , by condition 1.4, it means that one of the following mutually exclusive events must occur:  ∉ and ∈ ′  ∉ and ∉ ′  ∈ and ∉ ′ It follows that: ∉ ⇒ ( , ) ∉ , this suffices to show that condition . 1 cannot be attained under the first and second events. On the other hand: ∈ and ∉ ′ ⇒ ∉ ′ (∵ < under ) ⇒ ( , ) ∉ , which is enough to show that condition . 1 does not hold under the second scenario. Case 3) { , } ∈ , ∈ ′ , > ′ . In this case ( , ) ∈ because all blocks are complete; therefore, condition . 1 cannot be met. Case 4) ∈ , ∈ ′ , ∈ ′′ , ( ) > ( ′ ) > ( ′′ ). If ∉ ′′ it is clear that ( , ), ( , ) ∉ ∀ ∈ ∀ ∈ ′ disregarding the position of indices and within their corresponding blocks and the linkage relationship among the three blocks, thus . 1 does not hold and, therefore; in the following, only the case ∈ ′′ is considered. If each one of the three blocks is linked with the other two, then: ( , ) ∉ ⇒ ∉ or ∉ ′ ⇒ ( , ) ∉ or ( , ) ∉ (∵ condition 2.1). Similarly, if one block is linked to the other two, and these are the only existent links, condition 2.2 guarantees that . 1 does not hold. Consequently, under the two conditions of result 2 and condition 1.4, ordering is a perfect elimination ordering. Thus, the graph induced under approach 3.b with and edges set as defined in result 2 is decomposable.

Proof of necessity of condition 2.2
To prove necessity of condition 2.2, we show that if it is not satisfied, then is not a perfect elimination ordering. If there is at least one triplet of blocks { , ′ , ′′ } such that ( , ′ ), ( , ′′ ) ∈ , ( ′ , ′′ ) ∉ , ( ) < ( ′ ) < ( ′′) then by picking ∈ ′′ , ∈ ′ and ∈ ′ it follows that > > , and ( , ) ∉ , ( , ), ( , ) ∈ , that is, condition A.1 holds and as a consequence, is not a perfect elimination ordering, which implies necessity of condition 2.2. Proof of non-necessity of condition 2.1 By constructing a counter example, we refute the necessity of condition 2.1. Suppose that contains three blocks such that ( 1 ) = 1, ( 2 ) = 2, ( 3 ) = 3, 1 is only linked to 2 and 2 is only linked to 3 , 2 1 ≠ 2 3 and markers have been ordered according to . It is easy to see that in this graph, for any triplet { , , } such that according to , > > if ( , ) ∉ then ( , ) ∉ which implies that condition A.1 is never reached and consequently, is a perfect elimination ordering despite of the fact that condition 2.1 does not hold. Consequently, 2.1 is not a necessary condition.
∎ Proof of corollary to proposition 2 Let be the "super graph" defined in the corollary. It is assumed that conditions 2.1 and 1.4 hold. Consequently, we only need to check that condition 2.2 is satisfied. It follows immediately by noticing that condition 2.2 is nothing but the definition of a perfect elimination ordering of the blocks, i.e., a perfect elimination ordering of the "super vertices" of . ∎