Regression-based approach for testing the association between multi-region haplotype configuration and complex trait

Background It is quite common that the genetic architecture of complex traits involves many genes and their interactions. Therefore, dealing with multiple unlinked genomic regions simultaneously is desirable. Results In this paper we develop a regression-based approach to assess the interactions of haplotypes that belong to different unlinked regions, and we use score statistics to test the null hypothesis of non-genetic association. Additionally, multiple marker combinations at each unlinked region are considered. The multiple tests are settled via the minP approach. The P value of the "best" multi-region multi-marker configuration is corrected via Monte-Carlo simulations. Through simulation studies, we assess the performance of the proposed approach and demonstrate its validity and power in testing for haplotype interaction association. Conclusion Our simulations showed that, for binary trait without covariates, our proposed methods prove to be equal and even more powerful than htr and hapcc which are part of the FAMHAP program. Additionally, our model can be applied to a wider variety of traits and allow adjustment for other covariates. To test the validity, our methods are applied to analyze the association between four unlinked candidate genes and pig meat quality.


Background
Haplotypes, the linear arrangement of alleles on the same chromosome inherited as a unit, provide a natural framework for testing the association between genetic markers and complex traits more efficiently than separate marker analysis [1]. There is strong evidence that several mutations in cis position within a single gene can interact to create a "super allele" that has a large effect on the observed phenotype. The biological explanation for these haplotype effects is that several mutations in a gene cause several amino acid changes in the ultimate protein product, and the joint effect of these amino acid changes can have a much larger influence on the function of the protein product than any single amino acid change. This emphasizes the importance of examining candidate genes by SNP haplotyping. Some studies focus on haplotypes within a given genomic region [2][3][4][5][6][7]. Because complex traits are presumed to be the results of interaction by a set of genes which may be located in different regions, some methods aim to test gene-gene interaction, and interactions of single markers from different unlinked regions [8][9][10][11]. Specifically, Becker et al [12] reported a method to deal with haplotype interaction in unlinked regions for a binomial trait. They find the best haplotype combination from the unlinked regions by permutation, which is a modification of Ge et al [13]. However, this method could only be applied to case-control association testing, and could not include other covariates.
Generalized Linear Model (GLM) is an extension of the general linear modeling process that allows models to be fitted for several kinds of traits, such as Gaussian, Poisson, Binomial, etc., and allows various covariates. Schaid et al [5] introduced score statistics, which are receiving increased attention because they require only computation of the null estimates and are asymptotically equivalent to Wald and likelihood ratio statistics under both null and Pitman alternative hypotheses. Some methods that use score tests based on GLM to test haplotype-trait association have demonstrated the validity and power of this statistic [6]. However, these methods only considered one genomic region. If considering multi-region multi-marker haplotype configurations, a severe multiple-testing problem will occur. To obtain uncorrected P-values for a specific marker combination, we use an unnested simulation introduced by Becker and Knapp [2], which is based on the algorithm proposed by Ge et al [13].
We propose an alternative approach that uses score statistics based on GLM to build the statistic T over which some of the unlinked regions are considered and some markers are chosen at the selected regions. Since the distribution of T is generally unknown and is generally not comparable, we replace T with P min which is inherited from the algorithm of Becker and Knapp [2]. This simulation method has already been validated by Manly [14] and Hoh et al [15], and has systematically been applied to some genetic data [2,12,16,17].

Simulation schemes
We conduct a simulation study to evaluate the power and type I error of the association test and to compare our approach with others. The haplotype data are generated in a way similar to that of Roeder et al. [18] and Tzeng et al [6].
In every simulation scheme, we consider two unlinked regions. We consider that markers are in strong linkage disequilibrium within each region, but markers from different regions are in linkage equilibrium. Therefore, we separately produce two regions by using a modified Hudson's MS program [19]. This program generates data under a coalescent model in which the recombination occurs uniformly over the region. The 4 samples sizes are 50, 100, 200 and 400, respectively. The scaled recombination rate, ρ = N e δ/bp, is set to 4 × 10 -3 for the recombination cold spots, and 100 times greater in the hot spots, with the effective population size N e is 1 × 10 4 . The scaled mutation for the entire region, 4N e μ/bp, is set to be 6 × 10 -4 .
Once the haplotypes have been generated, the first step is to restrict the disease or minor allele frequency. In this simulation, we set the allele frequency as 3 levels: 0.1, 0.3, and 0.5. We assume that the middle locus in every gene is the liability locus. Once a liability locus is chosen, a haplotype is defined as a segment of three adjacent SNPs in which the second SNP is the liability locus.
After randomly pairing haplotypes to form individual genotypes, we generate both continuous and binary trait values.

Continuous traits
For the Type I error test, we consider two simple models of quantitative traits simulated independently of the liability locus. Let model 1 include only an environmental effect e: Y = e. Let model 2 additionally incorporate a covariate Z: Y = γ × Z + e. In the models, e follows a standard normal distribution with mean 0 and variance 1, and Z is generated from a standard normal distribution. For assessing power and the effective selection of the best combination of markers, we also consider two models of quantitative traits simulated in association with the liability locus. Let model 3 decompose the trait value into MRHC effect and environmental effect e: Y = g + e, and model 4 additionally incorporate a covariate Z: In these models, g is the sum of all considered genes' effects. For the ith gene, g i has a discrete distribution and equals , , with probabilities q 2 , 2q(1-q) and (1-q) 2  Under all scenarios, we compute the global P value with 1500 permutation replicates for each simulated data set. Empirical significance levels and power were computed as the portion of simulated data sets for which the global P value was less than or equal to 0.05.

Comparison of three models htr, hapcc and HAPGLM
In order to check the validity and the accuracy of our HAPGLM approach, we first carry out simulations under the null hypothesis and compare it with hapcc and htr, which were implemented in the beta version of FAMHAP. htr performs a haplotype trend regression test proposed by Zaykin et al [7], and hapcc performs a χ 2 test for haplotypes proposed by Becker et al [12,16]. Here, we use model 5  The reason is that the disease individual prevalence is 0.10, and the percent for the disease is somewhat small, making it difficult to find the significant difference. However, when the sample size is more than 100, the power is near one for the three models. There are no significant differences among hapcc, htr and our method.

Three factors analysis for global test
To evaluate the test performance, we describe the results from our power and type I error study that use the above methods with various parameters. Type I error test includes 48 scenarios which include 4 models, 3 minor allele frequencies and 4 sample sizes. As shown in table 4, the type I error stabilizes in all the scenarios. Power test includes 48 scenarios which include 4 models, 3 disease allele frequencies and 4 same sizes. For the power calculations in table 5, the power is adversely affected by the small sample size and high disease allele frequency. Otherwise, if the sample size is at least to 100, the power is preserved. Therefore, we set sample size to test recombination affection and the specific MRHC testing as 100.

Recombination analysis for global test
In order to check the effect of recombination on the model, we first consider two different recombination levels at which the diversity of the haplotype is high and low. "High diversity" indicates that a minor or disease locus is located in the region of recombination hot spots and that the number of distinct haplotypes is 6-8; "low diversity" indicates that a minor or disease locus is located in a haplotype block and the number of distinct haplotype is 3-5.
In this simulation, we consider three SNPs in each region and we assume equal recombination level. The first recombination level has 2-4 different haplotypes. The haplotype distribution of the second recombination level consists of 6-7 theoretically possible haplotypes. From table 6 and 7, there are some differences between the two diversities on type I error and power, but there are no significant differences between the high and low diversities. Binary Traits Produce from model 1 above a given threshold Model 6 Binary Traits Produce from model 2 above a given threshold Model 7 Binary Traits Produce from model 3 above a given threshold Model 8 Binary Traits Produce from model 4 above a given threshold That is to say, the proposed method is not significantly affected by the recombination level.

Analysis for specific MRHC
The score test can easily compute the specific MRHC. In order to study the performance of the proposed method in detecting individual MRHCs, we set 3 disease allele frequencies and 4 model fits. The power for specific MRHC is presented in table 8. It is of interest that the proposed method is robust for detecting the specific MRHC.

Select the best combination of markers
In order to check the accuracy of our methods to select the best MRHC, the markers of the 2 genes in simulations are 123|456, and the marker 2 and 5 are as the liability loci. The frequency of the disease allele is 0.3. The result shows that for these 2 liability loci combination, the statistic T is the largest and the P value is the smallest. These are presented in table 9. The combination without these two loci presents low T and high P value.

Application to a pig meat quality dataset
Meat quality is very important in the pig meat production industry. Many candidate genes have been identified that could be used to improve this trait through marker assisted selection (MAS) [20]. The Heart Fatty Acid-Binding (H-FABP) gene encodes a type of cytosol protein that transports fatty acids from the cell membrane to other sites where 3-acyl-glyceride and phospholipids are synthesized and fatty acids are oxidized. Gerbens [21][22][23] discovered Msp I, Hae II and Hinf I polymorphisms of the H-FABP gene that is related to intramuscular fat content. Melanocortin-4 Receptor (MC4R) is believed to be a link between feed intake and body weight [24].
Polymorphism of the MC4R gene has been reported to be associated with back fat thickness [25]. Adipocyte Determination and Differentiation factor-1 (ADD1) can activate or restrain some genes in fat and glucose metabolism.
Research has suggested that the ADD1 gene can be used as a candidate gene for pork quality [26]. Calpastatin (CAST), which is an endogenous inhibitor (Ca 2+ dependent   cysteine proteinase), plays a central role in the regulation of calpain activity in cellsand is considered to be one of the major modulators of the calpains [27,28]. The CAST gene represents an excellent candidate gene for studying variation in pork quality. We aim to find association between multi-region haplotype effects from these candidate regions and meat quality.
Our data set is a sample which includes 93 unrelated fatteners from the following breeds/populations: 18 Meishan, 21 Sutai, 14 Yorkshire × Sutai, 16 Landrace × Sutai and 24 Duroc × Landrace × Yorkshire pigs. 8 polymorphic markers of the preceding genes in the populations have been reported [29,30], which are: 2 in ADD1, 1 in H-FABP, 1 in MC4R and 4 in CAST. We code these polymorphic markers as (A1, A2, H1, M1, C1, C2, C3, and C4). The χ 2 test of these polymorphic markers show that there are significant differences in 5 polymorphic markers (except for A1, A2, C4) between the five populations. We set up single locus models including sex and breed as environmental covariates, for every polymorphic marker, and use statistical software SAS macro GLM for calculations. The results show significant effects at 0.05 level for: A1, H1, C2, C3 and C4 in back fat thickness (BK); A1, C1 and C3 in meat color (MC); A1, H1, C1 and C4 in intramuscular fat content (IMF), and A1 in protein content. To apply our methods, we also incorporate two environmental covariates (sex and breed), and use back fat thickness, tenderness, drip loss, meat color, intramuscular fat content, pH 1 hour after slaughter, pH 24 hours after slaughter, and the content of protein as the dependent variable in the regression model. Table 10 illustrates the marker combination in which raw P values are lower at 0.01 level, Compared to single locus model, MRHC analysis can detect more markers which are significantly associated with traits.
Additionally, our methods are used to reconstruct the distinct MRHC from the above 4 genes which are in unlinked

Discussion
Presently many publications have proven that the genetic dissection of complex traits depends not only on the identification of genes involved in disease susceptibility but also on the elucidation of the synergistic role that genes play with other genes and with environmental factors [8][9][10][11][31][32][33]. Therefore, considering unlinked genomic regions simultaneously is desirable. There are two models hapcc and htr in FAMHAP program which can compute the haplotype interaction in unlinked region. For hapcc, Becker et al. [12] chose the usual χ 2 test statistic for contingency tables which can be applied only to case-control traits. For htr, the haplotype trend regression test proposed by Zaykin et al [7], chose F statistic and could be used for qualitative and quantitative traits, but won't allow other covariates in the model. Our proposed methods are based on score equations for GLMs which allow adjustment of covariates and can model qualitative and quantitative traits. For binary trait without covariate, the type I error and power comparison show that our model has the same power as hapcc and htr, and type I error is as expected. For a small number of markers, the run times for hapcc, htr and HAPGLM are approximately equal. Additionally, our model has obvious advantages. First, our model can be applied to analyze haplotype association across independent regions with adjusting of covariates for a wider variety of traits. Second, the score statistic ( ) of the individual MRHCs can be easily computed.
Our model adopted the simulation method proposed by Becker et al. [2,12,17], which can be computationally feasible to deal with the multiple marker combination. For our program, the evaluation of a single simulated data set with 15 markers in 3 regions will take no more than 10 seconds on average on two nodes with 3.0 GHz Intel with 512 MB main memory. In general, it will be possible to simultaneously consider about 600 to 1,200 hypotheses on a standard PC. Our program is very flexible to allow selection of loci and genes for analysis. However, in regions with too many possible haplotype combinations, our program runs out of memory. We need to consider more aggressive trimming parameters or other haplotype estimation algorithms. For example, we will improve our program by using the haplotype ancestral cluster idea to cluster rare haplotypes with similar ancestral haplotypes, which was used by Tzeng et al. [34]. Haplotypes of the entire block can be represented by a smaller set of SNPs which are referred to as tag SNPs [35]. In order to analyze more markers, it will be helpful to select tag markers at each region and to carry out the analysis on the set of these markers. Tag SNPs selection will save run-time, and we plan our further research along this path.
In this research, we proposed markers from different regions which are proposed to be in linkage equilibrium. z k 2   Markers from different regions can be in linkage disequilibrium, but the methods can allow such markers as if they were in the same region.
Since the current model assumes that the subjects are independent of each other (i.e., unrelated), it is critical to extend the current approach to account for the correlations between subjects, given their family data. Therefore, further studies are needed to address the impact and modeling strategies with regard to the assumptions in the model. This model is restricted to outcomes that can be placed in the generalized linear model framework. An extension to failure-time data could also be placed in the framework.

Conclusion
It is quite common that the genetic architecture of complex traits involves many genes and their interactions. Therefore, dealing with multiple unlinked genomic regions simultaneously is desirable. We developed a regression-based approach which can be applied to a wider variety of traits and allow adjustment for other covariates to assess the interactions of haplotypes that belong to different unlinked regions. Multiple marker combinations at each unlinked region are also considered. In addition, HAPGLM can be downloaded for free at: ftp:// public.sjtu.edu.cn/, user: ylhu0323, password: public.

Contribution to multi-region haplotype configurations
Consider R unlinked genomic regions, and m r observed markers for each of n unrelated individuals in region r. Further, we suppose that markers within each region are in strong linkage disequilibrium, while markers from different regions are in linkage equilibrium. What we wish to investigate is whether some of the genomic regions are associated with the phenotype of interest via some of the markers from each region.
Let G r be the multi-locus genotype of an individual at region r, and h r be a haplotype of region r. If the haplotypes and are compatible with G r , is then called a haplotype explanation of G r . Obviously, for a given G r , there may be several haplotype explanations.
First, we use the expectation-maximization (EM) algorithm [5] to obtain the maximum-likelihood estimates of the haplotype frequencies at each of the unlinked regions. After we pooled the rare haplotypes (with estimated frequencies <0.001) into a single group, we adopt the following formula ( [12]     The possible number for the combination is 2 6 -1 with 6 markers in 2 regions. Let Z = X e |X g and γ = (α|β). Then, the likelihood of trait y i for subject i, given the vector Z i , can be expressed as a GLM for exponential family data [36] according to where a, b, and c are known functions, and ϕ is the dispersion parameter. To implement the score statistics for different types of traits, we need only assume a distribution for the trait and to make the appropriate substitutions for the expected value of the trait, , the dispersion parameter a(ϕ), the ratio b"(η)/a(ϕ) and the link function. (see  table 1 of [5]).

Score test for incorporating contribution of MRHC
We derive score statistics to test the null hypothesis of no association between MRHC and trait, H 0 : α = 0. Let ζ denote the vector of nuisance parameters (β, μ, ϕ). The likelihood function for (α, ζ) on the basis of the data according to Tzeng et al [6] is ,....,  Where P(x g, i ) is the contribution of individual i to MRHCs.
The score function for α is the partial derivative of the likelihood equation (5), with respect to α.
S α is asymptotically distributed as multivariate normal [36]. We consider the generalized score test, which would ensure the asymptotic null χ 2 distribution even under model misspecification [5]. Define θ = (α, ζ) and let V α denote the variance of S α , the equation can be expressed according to Louis [37] and Tzeng et al [6] as: where, In appendix B, we show the above result.
With the above results, we can compute a global score statistic according to The score statistic is distributed asymptotically as χ 2 with degrees of freedom equal to the rank of V α .
Schaid et al [5] proved that the score function for α and the score function for haplotype probabilities are independent under the null hypothesis, so that the covariance between the two score functions is zero. Since the contribution to each MRHC is estimated from the haplotype frequency that is used to calculate the score statistic S α , the variance of the score statistic is not penalized by the use of estimated haplotype frequencies.
In this framework, we can readily compute score statistics for each MRHC according to [5]: where z k follows χ 2 (1) under the null hypothesis H 0 : α = 0.

Testing more than one hypothesis
If we select m markers in several genes, there would be 2 m -1 marker combinations. To test 2 m -1 combinations with associated raw P values, and declare the global P value the significance level for our analysis would lead to another multiple-testing problem. In order to avoid nested simulation, we use the method which Becker and Knapp [2] adapted from Ge et al [13]. The basic idea is that, to test B = 2 m -1 marker combinations, global P is estimated by the proportion of permutation samples with min B P t smaller than that in the observed data, where t is the simulation time. For each marker combination B ∈ B and for each permutation replicate i = 1,..., t, the raw P value of the ith permutation replicate is calculated as