A Multiple-Trait Bayesian Variable Selection Regression Method for Integrating Phenotypic Causal Networks in Genome-Wide Association Studies

Bayesian regression methods that incorporate different mixture priors for marker effects are used in multi-trait genomic prediction. These methods can also be extended to genome-wide association studies (GWAS). In multiple-trait GWAS, incorporating the underlying causal structures among traits is essential for comprehensively understanding the relationship between genotypes and traits of interest. Therefore, we develop a GWAS methodology, SEM-Bayesian alphabet, which, by applying the structural equation model (SEM), can be used to incorporate causal structures into multi-trait Bayesian regression methods. SEM-Bayesian alphabet provides a more comprehensive understanding of the genotype-phenotype mapping than multi-trait GWAS by performing GWAS based on indirect, direct and overall marker effects. The superior performance of SEM-Bayesian alphabet was demonstrated by comparing its GWAS results with other similar multi-trait GWAS methods on real and simulated data. The software tool JWAS offers open-source routines to perform these analyses.

SEM-based GWAS (SEM-GWAS) methodology by applying SEM to linear mixed models for GWAS. They showed that while conventional GWAS methodology only provides overall SNP effects, SEM-GWAS can capture the complex causal relationships among traits and further decompose the overall SNP effects into direct and indirect effects.
The SEM-GWAS method proposed by Momen et al. (2018Momen et al. ( , 2019) is based on linear mixed models with a fixed substitution effect for the tested SNP and a random effect with covariances defined by a "genomic relationship matrix" computed from genotypes (VanRaden 2008) to account for genetic relatedness. Markers are usually implicitly assumed to affect all traits when the "genomic relationship matrix" is constructed in multi-trait analysis. However, this assumption is not biologically meaningful, especially in multi-trait analyses involving many traits. Cheng et al. (2018b) proposed a general class of multitrait Bayesian variable selection regression methods that use a broad range of mixture priors, e.g., multi-trait BayesCP, where each locus can affect any combination of traits, which allows us to more closely model the true biological mechanisms, e.g., pleiotropy (Cheng et al. 2018b).
The primary goal of this current research is to develop a multitrait Bayesian regression GWAS method that more closely resembles the underlying biological mechanisms including pleiotropy and causal structure among traits. In this paper, we develop and implement a new GWAS method called SEM-Bayesian alphabet, which integrates SEM to the multi-trait Bayesian variable selection methods, to incorporate the underlying biological mechanism. The term "Bayesian alphabet" denotes a collection of Bayesian regression models that differ in the priors adopted for marker effects (Gianola 2013). In this paper, we use SEM-BayesCP, a Bayesian variable selection method, to show the utility of the SEM-Bayesian alphabet. The performance of our proposed method is studied using real and simulated data.

MATERIALS AND METHODS
Multi-trait Bayesian regression model using mixture priors Assuming that individuals have all traits measured with a general mean as the only fixed effect, we write the multi-trait model for individual i from n genotyped individuals as: where y i is the vector of phenotypes of t traits for individual i, m is a vector of overall means for t traits, p is the number of genotyped loci, m ij is the genotype covariate at locus j for individual i (coded as 0,1,2), a j is the vector of marker effects of t traits for locus j, and e i is the vector of residuals of t traits for individual i. The fixed effects are assigned flat priors. The residuals, e i , are a priori assumed to be independently and identically distributed multivariate normal vectors with null mean and covariance matrix R, which is assumed to have an inverse Wishart prior distribution, W 21 S e ; n e ð Þ. Allowing each locus to affect any combination of traits, in a multiple-trait Bayesian variable selection method, e.g., multi-trait BayesCP (Cheng et al. 2018b), the vector of marker effects at locus j can be written as a where d jk is the indicator variable indicating whether the marker effect of locus j for trait k is zero or not, and b j is a priori assumed to be independently and identically distributed multivariate normal vectors with null mean and covariance matrix G, which is assumed to have an inverse Wishart prior distribution, W 21 t S b ; n b À Á . Given that a locus can have an effect on any combination of traits, we use numeric labels $1$; $2$; ⋯; $l$ to represent all 2 t possible combinations for d j , in which case the prior distribution for d j is: where P i is the prior probability that the vector d j corresponds to the vector labeled "i" and P P i ¼ 1. We assume the prior for P ¼ P 1 ; P 2 ; . . . P l ð Þis a uniform distribution.

Structural Equation Model
The linear SEM is composed of two parts: the measurement equation analyzing the relationship between the observable variables and latent variables, and the structural equation capturing the connections among latent variables (Anderson and Gerbing 1988). These two equations can be written as: where y i is the vector of observable variables for individual i, h i is a q · 1 vector of endogenous latent variables, j i is a r · 1 vector with exogenous latent variables, G 1 and G 2 are the matrix of unknown coefficients in structural equation, L is a t · q þ r ð Þ matrix of unknown structural coefficients, k 1 and k 2 are t · 1 and q · 1 vectors of residuals. The details of parameter estimation are discussed in Song and Lee (2012).
In our study, no latent variables are assumed and the sole observable variables are phenotypes. Thus only the causal relationship among observable variables, i.e., phenotypes, are fitted in the SEM model (also known as path analysis (Wright 1921)) as: where y i and L are defined as above, e i represents everything that is not explained by Ly i , and L is an t · t matrix of structural coefficients representing the causal structure recovered from the Inductive Causation (IC) algorithm as described in the next section.
To illustrate, we assume that the phenotypes of three traits for each individual (i.e., y 1 ; y 2 ; and y 3 for traits 1, 2, and 3) have the following causal relationship: where causal coefficient l ij represents that a 1-unit increase in trait i results in a l ij unit increase in trait j. Given the causal structure above, the L can be written as: (2) Searching causal structure As described above, fitting the SEM requires the causal structure among all traits to be known before analysis. To explore the widerange of possible underlying causal structures, we use the method from Valente et al. (2010) to discern the causal structure based on the posterior distribution of the residual covariance matrix. The reason we do not directly apply this method to phenotype data are that the covariance among phenotypes is likely confounded by genetic effects. The process of inferring causal structure is composed of three steps: 1. Fit the multi-trait BayesianCP model and obtain the posterior distribution of the residual covariance matrix. 2. Follow Valente et al. (2010) to derive the conditional independence relationship among traits based on the posterior distribution of the residual covariance matrix. In detail, we derive the residual partial correlation p(y i , y j |h), where h is a set of traits, to test whether trait y i is conditionally independent from y j . The highest posterior density (HPD) interval of 0.9 was used to make statistical decisions. If HPD interval of p(y i , y j |h) contains zero, y i and y j are regarded as conditionally independent on h. 3. Apply the IC algorithm (Pearl 2009) as described in the Appendix to the conditional independence relationship from step 2 to obtain the causal structure.

SEM-BayesCP
Assume e i ¼ m þ P p j¼1 m ij a j þ e i in equation (1) and follow assumptions in multi-trait BayesCP, the SEM-BayesCP model can be written as: Move Ly i from the right side to the left side of equation (3), and define L Ã ¼ I 2 L, where I is a t · t identity matrix and L is a t · t matrix of structural coefficients based on the discerned causal structure, the model becomes: To guarantee that the structural coefficient is identifiable, we assume that the residuals for each trait of individual i are independent with each other, which means the residual covariance matrix is diagonal (Wu et al. 2010;Momen et al. 2018). The vector of all non-zero elements in L, e.g., l ¼ l 12 ; l 13 ; l 23 ½ , is assumed to have a prior distribution: where 1 is a vector of ones, I is the identity matrix, and l 0 is a known mean for all elements in l. t 2 is a tuning parameter to adjust the sharpness degree of the prior (Gianola and Sorensen 2004). In this paper, we set l 0 ¼ 0 and t 2 ¼ 1. The priors for the remaining parameters are the same as in the section Multi-trait Bayesian regression model using mixture priors. Gibbs samplers are used to draw samples for all parameters. The full conditional distribution to draw samples for l is shown below. The derivations of the full conditional distributions of the remaining parameters of interest for Gibbs samplers are in Cheng et al. (2018b).
Full conditional distribution of L: We follow Gianola and Sorensen (2004) to obtain the full conditional distribution of L, with the difference between our derivation and Gianola and Sorensen (2004) being that we specify the causal structure with positions of parameters in the L. Let V denote all parameters except l in the SEM-BayesCP and use the causal structure L ¼ 0 0 0 l 12 0 0 l 13 l 23 0 0 @ 1 A as an example, the left hand side of equation (4), L Ã y i , can be written as: The conditional posterior distribution of l can be written as: (5) can be written as: Following the derivation in Gianola and Sorensen (2004) and the fact that L Ã j j ¼ 1 in a recursive system, the full conditional distribution of l is Decomposition of SNP effects In SEM-BayesCP, the marker effect for locus j, a j , is considered as the vector of direct marker effect of t traits. The indirect effect of locus Volume 10 December 2020 | SEM-Bayes-GWAS | 4441 j of t traits can be calculated as P t21 r¼1 L r a j . The overall effect of locus j on t traits is computed as P t21 r¼0 L r a j or I2L ð Þ 21 a j , which is the summation of both direct and indirect effect of locus j. For example, given a causal structure L ¼ 0 0 0 l 12 0 0 l 13 l 23 0 0 @ 1 A , the direct effect for locus j on three traits is a j ¼ a 1j a 2j a 3j 0 @ 1 A , and the indirect effect for locus j on three traits is calculated as and the overall effect of locus j on trait k is a 1j l 12 a 1j þ a 2j l 13 þ l 12 l 23 ð Þ a 1j þ l 23 a 2j þ a 3j 0 @ 1 A .

Inference of association based on genomic windows
Markers in a genomic window are usually highly correlated, indicating that any single marker may not show a strong association with the trait even though a causal variant exists in the window. In this paper, we make an inference of association based on genomic windows, because multiple markers inside a genomic window may jointly capture the signal from the causal variant (Fernando and Garrick 2013;Fernando et al. 2017).
To make an inference of association based on genomic windows, posterior distribution for the proportion of the genetic variance explained by markers in genomic window w, q w , is estimated from MCMC samples of overall, direct, and indirect marker effects as follows. For one MCMC sample of all marker effects on one trait, let a direct , a indirect , and a overall denote direct, indirect, and overall effects of all markers respectively.
The genetic value that is attributed to genomic window w is calculated as: where M w is a matrix of marker covariates in window w and a w;direct , a w;indirect , and a w;overall are the MCMC samples of direct, indirect, and overall marker effects for SNPs in window w. Then the variance explained by the genomic window w is defined as: Similarly, the total genetic variance is calculated as: The proportion of the genetic variance explained by direct, indirect, and overall marker effects in the genomic window w is calculated as: Given the MCMC samples of q w , the window posterior probability of association (WPPA) is calculated as the proportion of MCMC samples of q w that exceed a specific value T (Fernando and Garrick 2013;Chen et al. 2017;Lloyd-Jones et al. 2017). In this paper, associations are tested for non-overlapping windows of 100 SNPs, and genomic windows that explain over 1 N of the total genetic variance were deemed to be of potential interest (i.e., T ¼ 1 N , where N is the total number of windows).

Data analysis
Real data: The Rice Diversity Panel with 413 Oryza sativa individual accessions was used in the analysis. Three traits were considered, including plant height (PH), flowering time in Arkansas (FTA), and panicle number per plant (PN) in our GWAS. After removing the records with missing data for these three traits and genotype with minor allele frequency , 0:05, 370 individuals with 33,519 SNPs genotyped were included in our analysis. The phenotypic and genotypic data were publicly available for download from http:// www.ricediversity.org/. It has been shown that using a threshold of WPPA ¼ a to declare a significant genomic window restricts the proportion of false positives (PFP) to , 1 2 a (Fernando et al. 2017). A previous GWAS (Zhao et al. 2011) identified significantly associated SNPs in chromosome 6 for flowering time in Arkansas (FTA) using the same dataset. A threshold of WPPA = 0.8 and p-value = 5·10 26 in our GWAS analysis resulted in similarly significant signals. This result suggests that a WPPA of 0.8 and p-value = 5·10 26 are reasonable for declaring a significant genomic window.
Simulated data: To compare SEM-BayesCP with SEM-GWAS of Momen et al. (2018), we simulated data based on real genotypes from the Rice Diversity Panel. The simulation scenarios in Chen et al. (2017) were applied to simulate different genetic architectures. The QTL effects were generated from unit-gamma distribution (scale = 1) with three different shape parameters (g): fewer QTL with large effects (g ¼ 0:18), fewer QTL with small or large effects and many QTL with intermediate effects (g ¼ 3:0), and the intermediate case (g ¼ 1:48). In addition to the distribution of QTL effects, the number of QTL (n QTL ) may play an important role in GWAS, thus three numbers of QTL (n QTL ¼ 30; 90; 300) combined with the three shape parameters were used to create 9 scenarios. For each scenario, 50 replicated populations were simulated with QTL positions randomly sampled across the genome. Trait 1 was assumed to have a causal effect on trait 2 with causal structural coefficient l ¼ 1:0. The QTL effects were simulated under two scenarios (Figure 1): the QTL have direct effect on both trait 1 and trait 2 in scenario 1, where QTL only have direct effect on trait 1 in scenario 2. Half of the QTL were simulated following scenario 1, while the remaining followed scenario 2. Phenotypes for two traits were generated based on heritability of 0.5. In the simulated data analysis, the causal structure was assumed known to exclude the bias caused by searching causal structures. The structural coefficients were assumed known in SEM-GWAS.
We have implemented SEM-Bayesian alphabet in JWAS (Cheng et al. 2018a), an open-source, publicly available package for single-trait and multi-trait whole-genome analyses written in the freely available Julia language. More details can be found at https://reworkhow.github.io/JWAS.jl/latest/.

Simulated data result
The performances of SEM-BayesCP and SEM-GWAS were compared based on the AUC (Area under Receiver Operating Characteristic). Inference of association on genomic windows for SEM-GWAS was based on the minimum p-value (Begum et al. 2016), i.e., genomic windows containing at least one significant variant are declared as significant windows. To exclude the irrelevant AUC with low levels of specificity, only the partial area under the curve up until the false positive rate of 5% (pAUC5) (Chen et al. 2017;Ma et al. 2013) was calculated. For the convenience of comparison, all pAUC5 measurements were rescaled such that the pAUC5 of the random classifier is equal to 1. The R package ROCR (Sing et al. 2005) was used to obtain the pAUC5; the paired t-tests (p-value , 0.1) were used for comparing both the pAUC5 mean across all scenarios (overall mean comparison) and the pAUC5 mean for each level of n QTL and shape parameters g (marginal mean comparison).
The GWAS results based on overall, direct and indirect effect were shown in Table 1. For the overall effect result on trait 1, there is no significant difference between SEM-BayesCP and SEM-GWAS in both overall mean comparison and marginal mean comparison. The direct effect on trait 1 is the same as the overall effect on trait 1 since the direct effect and overall effect on trait 1 are equal based on the causal structure. For the overall effect on trait 2, the pAUC5 mean of SEM-BayesCP is significantly higher than that of SEM-GWAS in the overall mean comparison, and some marginal mean comparisons (e.g., n QTL ¼ 30 and g ¼ 0:18). For the direct effect on trait 2, though higher overall mean of pAUC5 is usually observed in SEM-BayesCP, there is no significant difference (p-value , 0.1) between SEM-BayesCP and SEM-GWAS in both overall mean comparison and marginal mean comparison. For the indirect effect result on trait 2, similar to the overall effect result on trait 2, the pAUC5 mean of SEM-BayesCP is significantly higher than that of SEM-GWAS in the overall mean comparison, and some marginal mean comparisons (e.g., n QTL ¼ 30 and g ¼ 0:18).

Real data result
Causal structure and structural coefficients: The causal structure among three traits is inferred by the IC algorithm from the estimated posterior distribution of the residual covariance matrix in the multitrait BayesCP model. Figure 2 shows three potential phenotypic causal structures among traits PH (y 1 ), FTA (y 2 ), and PN (y 3 ) recovered for the 0.9 HPD interval. The causal structure matrices for IC1 (L 1 ), IC2 (L 2 ), and IC3 (L 3 ) are: Figure 1 a 1 : the direct effect on trait y 1 ; a 2 : the direct effect on trait y 2 ; la 1 : the indirect effect on trait y 2 . The graph shows the simulated two QTL effect simulation scenarios. Scenario 1 represents that the QTL has direct effect on both trait 1 and trait 2, whereas scenario 2 represents that the QTL only have direct effect on trait 1.
n■ Table 1 Overall and marginal mean of rescaled pAUC5 of SEM-BayesCP and SEM-GWAS based on overall,direct and indirect effect  These three causal structures are fitted in the SEM-BayesCP model, and samples from posterior distributions for coefficients in these causal structures are obtained. The 90% credible intervals for structural coefficients in IC1, IC2, and IC3 are shown in Table 2. It is worth noting that the causal structures in Figure 2 provides the same set of marginal and conditional independencies, extra biological knowledge is required to further infer the causal structure. In this paper, the SEM-BayesCP model is proposed to incorporate (known) underlying causal structure among traits for GWAS. Thus, for the simplicity of our presentation, the causal structure is assumed known in Real Data Result section, e.g., IC1 is used to demonstrate the performance of SEM-BayesCP.
Decomposition of SNP effects: Direct, indirect, and overall SNP effects for all markers are estimated from SEM-BayesCP and SEM-GWAS. In SEM-BayesCP, direct SNP effects are assigned mixture priors, where each locus can affect any combinations of traits directly; samples from posterior distributions of indirect effects are obtained using joint samples from posterior distributions of L and direct SNP effects a j . In IC1, for trait PH, the overall SNP effect is equal to the direct SNP effect, because there is no intermediate trait.
For trait FTA, the overall SNP effect is composed of direct SNP effect and indirect SNP effect transmitted from PH. So the overall SNP effect for FTA is given by summing the direct SNP effect and indirect SNP effect. Similarly, for trait PN, the overall SNP effect is obtained by summing the direct SNP effect and indirect SNP effect transmitted from PH.
The results of GWAS from SEM-BayesCP and SEM-GWAS incorporating causal structure IC1 are shown in Figure 3. Significant signals are found only for trait FTA. SEM-BayesCP adopts a threshold of WPPA = 0.8 to declare a significant genomic window and SEM-GWAS adopts a threshold of p-value = 5 · 10 26 . The overall SNP effects are partitioned into direct and indirect effects, and GWAS are performed for the direct, indirect, and overall SNP effects separately for trait FTA. In Figure 3, the blue points, pink points and green points represents the significant genomic windows located in chromosome 1, chromosome 5 and chromosome 6. Window A contains SNPs from "id1000759" to "id1001229"; window B contains SNPs from "id1023967" to "id1024499"; window C contains SNPs from "id5013234" to "id5013920"; window D contains SNPs from "id6005814" to "id6006470".
For the overall effects, in the SEM-GWAS, window D achieved -log(p-value) 14.21; in the SEM-BayesCP, window C achieved WPPA 0.90, window D achieved WPPA 0.88, and window A achieved WPPA 0.82. For the direct effect, in the SEM-GWAS, window D achieved -log(p-value) 12.24; in the SEM-BayesCP, window C achieved WPPA 0.90, window D achieved WPPA 0.86, and window A achieved WPPA 0.80. For the indirect effect, in the SEM-GWAS, window B achieved -log(p-value) 14.65; in the SEM-BayesCP, although no window is identified as significant in SEM-BayesCP, a peak was observed at window B with WPPA 0.52. Further, for all three effects, the results from SEM-BayesCP and SEM-GWAS are correlated (the correlation between the WPPA from SEM-BayesCP and -log(p-value) from the SEM-GWAS is higher than 0.5). The correlation of indirect effect results from these two methods results achieved 0.70. Also, for both SEM-BayesCP and SEM-GWAS, the overall effect is more correlated with direct effect rather than indirect effect. The magnitudes for overall, direct, and indirect SNP effect in SEM-BayesCP are also shown in Figure 5. Though most large overall SNP effects consist of a large direct SNP effect and a relatively small indirect SNP effect, the indirect effect of some SNPs play an important role, e.g., the overall effect of SNP "id1024159", as shown in Figure 5, consists of a large indirect SNP effect and relatively small direct effect.
Genetic pleiotropy in SEM-BayesCP As shown in Figure 4, the posterior distribution of the parameter P is obtained, and markers show different levels of pleiotropy for direct SNP effects. In SEM-BayesCP, each SNP can have direct effects on any combination of traits, and the parameter P is used to estimate the proportion of SNPs having different levels of pleiotropy. Indirect SNP effects on one trait are transmitted from direct SNP effects on other intermediate traits. For example, in causal structure IC1, the indirect SNP effects on FTA is transmitted from direct SNP effects on the trait PH. A SNP having no direct effect on FTA may have overall effects on FTA if its direct effect on trait PH is non-zero. The proportions of markers affecting different combinations of traits through overall SNP effects can also be estimated. This result is shown in Figure 4, and different probabilities are observed for some cases between overall and direct SNP effects. More SNPs have effects on all traits simultaneously when overall SNP effects are considered compared to direct SNP effects (case 1). If only overall SNP effects are considered, n■ Table 2 The 90% credible interval for causal structural coefficients in the three causal structures some cases having non-zero probabilities for indirect SNP effects are hidden by the causal relationships among traits (cases 2-4). The same patterns for direct and overall SNP effects are observed in cases 5-8 because there is no causal relationship between trait FTA and PN.

DISCUSSION
The complex causal relationships among multiple traits are usually not considered in conventional multi-trait GWAS.
Here we propose the SEM-Bayesian alphabet method to incorporate pre-inferred causal structures among multiple traits into multi-trait Bayesian regression methods. SEM-Bayesian alphabet accounts for causal structures among traits, and has the potential advantage of estimating causal effects, providing genomic window-based inference, as well as providing a comprehensive understanding of the underlying biological mechanism.

GWAS
To show the potential utility of SEM-Bayesian alphabet, simulated data were used to compare SEM-BayesCP with SEM-GWAS. A wide variety of potential genomic architectures were constructed by the combination of different levels of skewness of gamma distribution for QTL effects (g) and different numbers of QTL(n QTL ) (Chen et al. 2017). BayesCP was also performed to estimate overall SNP effects on the same datasets, and similar results as those in the SEM-BayesCP were obtained (not shown in this paper). This is reasonable since the SEM-Bayesian alphabet model can be reduced to a model similar to Bayesian regression by reparameterization, indicating that the joint likelihood functions of SEM-Bayesian alphabet and Bayesian regression are similar. Compared to Bayesian regression, the SEM-Bayesian alphabet provides a more comprehensive understanding of the underlying biological mechanism by decomposing overall SNP effects into direct and indirect SNP effects (Figures 3, 4 and 5).
The comparison between SEM-BayesCP and SEM-GWAS using simulated data were shown in Table 1. As shown in our results, SEM-BayesCP has relatively the same or higher pAUC5 than SEM-GWAS in all simulation scenarios. In some scenarios, SEM-BayesCP has significantly higher pAUC5 than SEM-GWAS. For example, when one trait is affected by few QTL of large effects (e.g., n QTL = 30,shape g ¼ 0:18 ð Þ ), SEM-BayesCP has significantly higher pAUC5 than SEM-GWAS to infer indirect and overall effects. Though significant difference is not observed for direct Figure 3 GWAS results based on overall, direct and indirect SNP effects from SEM-BayesCP and SEM-GWAS incorporating IC1 causal structure for the trait flowering time at Arkansas (FTA). The horizontal dash line represents the threshold 0.8 or -log(5 · 10 26 ). X-axis represents the location of genomic windows along the 12 chromosomes; Y-axis represents window posterior probability of association (WPPA) for SEM-BayesCP and negative logarithm of the p-value (-log(p)) for SEM-GWAS. Colored points represent genomic windows with WPPA $ 0.8 or p-value # 5 · 10 26 . The blue points, pink points and green points represents the significant genomic windows located in chromosome 1, chromosome 5 and chromosome 6. Figure 4 Estimated proportion of markers affecting combinations of traits, P, from SEM-BayesCP incorporating IC1 causal structure. For each scenario, it is estimated for both direct SNP effects (left) and overall SNP effects (right). effect, higher overall mean of pAUC5 is usually observed in SEM-BayesCP.

Causal structure
The causal structure is assumed to be known in SEM-Bayesian alphabet, and it is usually discerned by three types of algorithms: the constraint-based algorithm, the score-based algorithms, and the hybrid algorithms. The IC algorithm (Pearl 2009;Valente et al. 2010) used in this paper is a typical constraint-based algorithm, which is based on conditional independence tests. The score-based algorithms apply the heuristic optimization techniques, which set an initial graph structure and assign an initial goodness-of-fit score to it, and then maximize the goodness-of-fit score to obtain the most possible causal structure. The hybrid algorithm is a hybrid of both the constraintbased and the score-based algorithms. It utilizes conditional independence tests to reduce the space of candidate causal structures, and uses network scores to identify the optimal structure among them (Scutari 2014). The causal structures inferred from these algorithms may be different. Note that different evaluation criteria may also result in different outcome causal structures. For example, in this paper, if we choose 0.99 instead of 0.9 HPD interval to search for causal structures, there will be no edge between the traits PH and PN.

Decomposition of SNP effects
In some previous analysis (Mi et al. 2010;Momen et al. 2018Momen et al. , 2019, the indirect SNP effect of locus j of t traits is obtained by multiplying the estimated L,L, and estimated direct SNP effects,â j , as P t21 r¼1L râ j . This is similar to using posterior means of causal structural coefficients and direct SNP effects for calculation of the indirect SNP effects. In our method, indirect SNP effects are estimated using joint samples from posterior distributions of L and a j . We compared these two approaches for indirect SNP effect estimation on real rice data, and found that the indirect effects estimated from these two approaches are slightly different. The SEM-BayesCP approach should be used in indirect SNP effect estimation due to the fact that L and a j may be highly dependent. CONCLUSION SEM-Bayesian alphabet provides more interpretation into biological mechanisms than Bayesian regression methods by decomposing the overall SNP effects into direct and indirect SNP effects. In SEM-Bayesian alphabet, posterior distributions of the overall, direct, and indirect SNP effects, as well as causal structure coefficients, are obtained, which are used to make inferences about these parameters. Compared to the typical GWAS method incorporating causal structure among multiple traits, such as SEM-GWAS, SEM-Bayesian alphabet obtains the posterior distributions for the proportion of variance attributed to a genomic region to detect causal loci (i.e., the use of WPPA). The level of gene pleiotropy, e.g., proportion of markers affecting different combinations of traits as shown in Figure 4, can also be further dissected into direct and indirect SNP effects. Also, with estimating structural coefficients, SEM-Bayesian alphabet still has relatively same or greater pAUC5 than SEM-GWAS in all scenarios of simulated data. In summary, SEM-Bayesian alphabet offers a more comprehensive understanding of the underlying biological mechanisms including pleiotropy and causal relationships among traits than conventional GWAS, as well as has a potential advantage in the GWAS inference than other GWAS considering complex causal effect among multiple traits.