A computational model for sex-specific genetic architecture of complex traits in humans: Implications for mapping pain sensitivity

Understanding differences in the genetic architecture of complex traits between the two sexes has significant implications for evolutionary studies and clinical diagnosis. However, our knowledge about sex-specific genetic architecture is limited largely because of a lack of analytical models that can detect and quantify the effects of sex on the complexity of quantitative genetic variation. Here, we derived a statistical model for mapping DNA sequence variants that contribute to sex-specific differences in allele frequencies, linkage disequilibria, and additive and dominance genetic effects due to haplotype diversity. This model allows a genome-wide search for functional haplotypes and the estimation and test of haplotype by sex interactions and sex-specific heritability. The model, validated by simulation studies, was used to detect sex-specific functional haplotypes that encode a pain sensitivity trait in humans. The model could have important implications for mapping complex trait genes and studying the detailed genetic architecture of sex-specific differences.


Background
Differences in males and females (sexual dimorphism) is ubiquitous in many biological aspects [1][2][3]. In humans, sexually dimorphic traits include those from morphological shapes and body size to brain development to disease susceptibility [4,5]. Substantial differences are also observed in sensitivities to pain and pain-killing drugs, and susceptibility to developing chronic pain between men and women [6][7][8]. All these sex-specific differences are due to varying expression of genes on the X/Y chromosome and autosomes, thought to result from differences in cellular and hormonal environments between the two sexes [9]. A growing body of research has been conducted to elucidate the genetic control of sexual dimorphism in various complex phenotypes by gene mapping approaches [4,5,10,11]. Despite these efforts, however, little is known about the genetic architecture underlying sex-related variation in a quantitative trait.
Since sex is easily determined, the effects of sex on morphological, developmental and pathological traits can be directly observed. However, characterizing the impacts of sex on the genetic architecture of these traits has been challenged by a lack of powerful statistical approaches. The motivation of this article is to develop a statistical and computational model that can systematically search for sex-specific genes contributing to quantitative variation and formulate testable hypotheses regarding the interplay between sexes and gene expression. Our model is principally different from those used in many previous studies that are aimed to detect sex-specific quantitative trait loci (QTLs) based on linkage or linkage disequilibrium analysis [2,4,5,12,13]. Our model will be founded on the statistical framework constructed by Liu et al. [14] to detect the effects and diversity of haplotypes constructed by single nucleotide polymorphisms (SNPs) that are genotyped at candidate genes or genome-wide [15]. Our model has been generalized to allow the test of sex differences in haplotype frequencies, allele frequencies and linkage disequilibria between different SNPs as well as additive and dominant effects of haplotypes on complex traits. It has power to identify sex-specific DNA sequence variants that encode complex phenotypes in men and women.

Notation
Suppose there is a diversity of haplotypes constructed by two SNPs each with two alleles designated as 1 and 0. Let p and q be the 1-allele frequencies for the first and second SNP, respectively. Thus, the 0-allele frequencies at these two different SNPs will be 1 -p and 1 -q. The two SNPs that are segregating in a natural human population form four haplotypes, [11], [10], [01], and [00], whose frequencies are constructed by allele frequencies and linkage disequilibrium (D) between the two SNPs, i.e., The parameters contained in equation (1) can be used to describe some important aspects of the genetic structure and diversity of a natural population. Thus, differences in genetic architecture between the two different sexes can be characterized by these sex-specific parameters. Because it is easy to derive the closed forms for estimating haplotype frequencies [14], we will estimate the linkage disequilibrium from the estimated haplotype frequencies.
Let and be the vectors of haplotype frequencies among males and females, respectively. All the genotypes for the two SNPs are consistent with diplotypes, except for the double heterozygote, 10/10, that belongs to a diplotype of either [11] [00] or [10] [01] ( Table 1). Assuming that the population is at Hardy-Weinberg equilibrium, the frequency of a diplotype is expressed as the product of the frequencies of the two haplotypes that construct it.      is the relative proportion of diplotype [11] [00] within the double heterozygote for sex k.
It can be seen from equation (3) or (4) that maximizing is equivalent to maximizing log and log individually in equation (5).

The EM algorithm
A closed-form solution for the EM algorithm [14] has been derived to estimate the unknown parameters that maximize the two sex-specific likelihoods of (5). The estimates of sex-specific haplotype frequencies are based on the log-likelihood function , whereas the estimates of sex-specific genotypic values of composite diplotypes and the residual variance are based on the loglikelihood function . These two different types of parameters can be estimated using a twostage hierarchical EM algorithm (see ref. [14] for a detailed implementation).

Model selection
According to equation (5), the summed likelihood across the sexes, , is formulated by assuming that haplotype [11] is a risk haplotype. However, a real risk haplotype is unknown from raw data (y k , S k ). An additional step for choosing the most likely risk haplotype should be implemented. The simplest way to do so is to calculate the likelihood values by assuming that any one of the four haplotypes can be a risk haplotype (Table 1). Thus, we obtain four possible likelihood values as follows: The largest likelihood value calculated is thought to correspond to the most likely risk haplotype. Under an optimal risk haplotype, we estimate sex-specific quantitative genetic parameters and .

Hypothesis tests
The genetic architecture of a quantitative trait is characterized by population (including haplotype frequencies, allele frequencies, and linkage disequilibria) and quantitative genetic parameters (including haplotype effects and mode of inheritance for haplotypes). The model proposed provides a meaningful way for estimating the genetic architecture of a trait and further testing sex-specific differences in genetic control.
After haplotype frequencies are estimated, allele frequencies and linkage disequilibrium between the two SNPs with each sex can be calculated as The influence of haplotypes on a quantitative trait is quantified in terms of the additive (a) and dominant genetic effects (d), and the mode of inheritance (ρ), which are estimated for each sex. Each of these population and quantitative genetic parameters can be tested when appropriate hypotheses are formulated.

Overall genetic control
Haplotype effects on the trait, i.e., the existence of functional haplotypes, in both male and female populations can be tested using the following hypotheses expressed as The log-likelihood ratio test statistic (LR) under these two hypotheses can be similarly calculated, where the L 0 and L 1 are the plug-in likelihood values under the null and alternative hypotheses of (8), respectively. Although the critical threshold for determining the existence of a functional haplotype can be based on empirical permutation tests, the LR may asymptotically follow a χ 2 distribution with four degrees of freedom, so that the threshold can be obtained from the

Sex-specific population genetic architecture
The male and female populations may be different in terms of population genetic parameters. Such sex-specific f k p k p k p k p k p k p k = + 11 00 00 00 10 01  At least one equality in p ΘΘΘΘ differences can be tested by formulating the following hypotheses for allele frequency at SNP 1, for allele frequency at SNP 2, and for the linkage disequilibrium between the two SNPs.
For each of the hypotheses (10)-(12), the LR values are calculated, which are each thought to asymptotically follow a χ 2 -distribution with one degree of freedom. Sex-specific differences in overall population genetic architecture can be tested with the null hypothesis H 0 : and D M = D F , with the corresponding LR value to be χ 2distributed with three degrees of freedom.

Sex-specific quantitative genetic architecture
Sex-specific differences in overall quantitative genetic architecture can be tested by formulating the hypotheses The LR value calculated under the null and alternative hypotheses is suggested to follow a χ 2 -distribution with two degrees of freedom. The rejection of the null hypothesis implies that the effects of the same haplotype are different between the two sexes. If there exists a sex-specific difference, the next step is to test whether this difference is due to the additive or dominant genetic effects, or both.

Sex-specific risk haplotypes
In the preceding sections, the same risk haplotype was assumed between the male and female populations. It is possible that the two sexes have different risk haplotypes. Let (j m = 2, 1, 0) and (j f = 2, 1, 0) be the genotypic values of composite diplotypes for the males and females constructed by a sex-specific rick haplotype. By reformulating the likelihood log of equation (5) based on sex-specific composite diplotypes, these genotypic values can be estimated with the EM algorithm. A best combination of risk haplotypes between the two sexes can be determined from the AIC values.

L-SNP model
It is possible that the two-and three-SNP models are too simple to characterize genetic variants for quantitative variation. We can develop a model that includes an arbitrary number of SNPs whose sequences are associated with the phenotypic variation. A key issue for the multi-SNP sequencing model is how to distinguish among 2 -1 different diplotypes for the same genotype heterozygous at loci. The relative frequencies of these diplotypes can be expressed in terms of haplotype frequencies.
Consider a a functional haplotype that contains L SNPs among which there exist linkage disequilibria of different orders. The two alleles, 1 and 0, at each of these SNPs are symbolized by r 1 , ..., r L , respectively. Let be the allele frequencies for these different SNPs within sex k. A haplotype frequency, denoted as , is decomposed into the following components: where D k 's are the linkage disequilibria of different orders among particular SNPs for sex k.
Sex-specific difference in terms of allele frequencies and linkage disequilibria between different SNPs as well as haplotype additive and dominance effects can be tested by formulating the corresponding hypotheses.

Pain genetics study
The model proposed was used to detect differences in the genetic architecture of pain sensitivity between men and women. Genetic and phenotypic data were from a pain genetics project in which 237 subjects (including 143 men and 94 women) from five different races were sampled for six SNPs at three candidate genes. As a demonstration of the utilization of the model, we will focus on two SNPs, OPRDT80G (with two alleles T and G) and OPRDT921C (with two alleles T and C), at the delta opioid receptor. Pain testing procedures followed Fillingim et al. [16]. The phenotypic values of traits were subtracted by the means for each race to remove the effect due to races. The two SNPs are highly significantly associated at p = 3.41 × 10 -5 for males and p = 9.63 × 10 -5 for females, with a normalized linkage disequilibrium of D' = 1.00, because alleles T from OPRDT80G and T from OPRDT921C as well alleles G from OPRDT80G and C from OPRDT921C tend to form the same haplotypes more frequently than at random. There is no sex-specific difference in allele frequencies at the two SNPs and their linkage disequilibrium.
By assuming that one of the haplotypes is a risk haplotype, we estimated the effects of each haplotype on a pain sensitivity trait, assessed with a baseline pressure pain threshold measured at the ulna, in the pooled male and female population. The genotypic values of composite diplotypes constructed by this risk haplotype and its non-risk haplotype counterpart were estimated and compared between different sexes. In both males and female, the three composite diplotypes do not display significant genetic differences in the pain trait studied, but the directions of the additive and dominance effects are different between the two sexes ( Table 2). In males, the non-risk haplotype tends to increase pressure pain thresholds, and it is overdominant to the risk haplotype, leading to increased pressure pain thresholds at a marginal significance level (p = 0.058) (Fig. 1). By contrast, in females, the non-risk haplotype tends to reduce pressure pain thresholds, and it also tends to be overdominant to the risk haplotype by reducing pressure pain thresholds. These discrepancies in both effect size and direction (Fig. 2) make the overall quantitative genetic architecture of the pain sensitivity trait significantly different between the two sexes (p = 1.49 × 10 -7 ) ( Table 2). Although the additive genetic effect displays a gene by sex interaction at the p = 0.03 significance level, a gene by sex interaction for the dominance effect is highly significant at p = 6.16 × 10 -7 . No significant difference was observed in inheritance mode between males and females.

Monte carlo simulation
Simulation studies were performed to test the statistical properties of the model proposed. Given a certain sample size (n), we simulated two SNPs by assuming different allele frequencies and linkage disequilibria between two sexes. The new model was used to analyze the simulated SNP and phenotypic data with the results tabulated in Table 3. Population genetic parameters including allele frequencies (p and q) and linkage disequilibrium (D) can well be estimated, with increasing precision when the sample size increases from 100 to 400. The power to detect the given sex-specific differences in allele frequencies and linkage disequilibrium is as high as 0.95 even with a modest sample size (100). Although the traditional model that does not implement sex-specific differences can provide precise estimates of these parameters, the estimates are generally between the true values of males and females. The estimation of quantitative genetic parameters including the additive (a) and dominance effects (d) needs the determination of an optimal risk haplotype. When all possible risk haplotypes were assumed for the simulated data, we found that the true risk haplotype gave the largest likelihood among the four possible cases. In general, quantitative genetic parameters can well be estimated, but the estimation precision increases dramatically with sample size and heritability. Although the additive effect can be obtained with reasonable precision at a modest heritability (0.1) with a modest sample size (100), the precise estimation of the dominance effect relies upon a larger heritability and sample size. Also, the given difference in the additive effect between two sexes can be detected with great power, even when both the sample size and heritability are small. But the same size of sex-specific difference in the dominance effect can be detected with the same power only when the sample size is 400 and heritability is 0.4 ( Table 3). The traditional model gave biased estimates of the sex-specific additive and dominance effects regardless of increasing sample size and heritability.
We conducted an additional simulation study, in which the data simulated under the assumption of no sex-specific differences in all genetic parameters were analyzed by the new and traditional model. As expected, both the models provide reasonable estimates of population and quantitative genetic parameters, with estimation precision increasing with increasing sample size and heritability (data not shown). This, in conjunction with the results in Table 3, suggests that the new model provides a general tool for study the genetic architecture of a complex trait, regardless of whether the genetic control of the trait is sexspecific.

Discussion
The genetic architecture of a quantitative trait is complex in terms of interactions between its underlying genetic factors and various environments including sex [1][2][3]. However, in many current studies, gene by sex interactions are often ignored simply because existing analytical models are not incorporated by environmental factors. While different phenotypic expressions of a trait between the two sexes can be easily measured [4,5], sex-specific discrepancy in the genetic control of the trait can be discerned only when a sophisticated model is used. There is strong evidence for sex-specific genetic influence [17,18] even for traits that display no sexual dimorphism [19]. Although quantitative genetic models are available to estimate sexspecific heritabilities due to aggregative effects of many genes [5,[19][20][21][22] or map sex-specific QTLs for phenotypic variation [4,5,[11][12][13]23], the model proposed in this article can dissect sex-specific genetic control at the DNA sequence levels.
Our model is founded on the conceptual framework for haplotyping a trait with single nucleotide polymorphisms (SNPs) formulated by Liu et al. [14]. Since haplotypes constructed by physically associated SNPs are thought to affect the expressivity of a complex trait [23,24], it is more precise to characterize such haplotype effects by incorporating gene by gene and gene by sex interactions. Lin and Wu [25]   lotype-haplotype interactions. The new model reported here can not only estimate sex-specific genetic parameters, but also provide a series of statistical procedures for testing sex-specific differences in the genetic architecture of quantitative variation. In a natural population, the structure and pattern of genetic variation can be studied by population genetic parameters, such as haplotype frequencies, allele frequencies and linkage disequilibria. Thus, the understanding of differences in these parameters between the two sexes help to infer the sex-specific genetic structure of a natural population and its evolutionary processes. As shown through simulation studies, our model is alert to discern sex-specific differences in basic population genetic parameters.
Liu et al. [14] assumed a so-called reference or risk haplotype that triggers an effect on complex traits in a different way from the other haplotypes. Thus, the combinations between the risk and non-risk haplotypes (composite diplotypes) will perform differently, depending on the type of combination, i.e., risk by risk, risk by non-risk and non-risk by non-risk. Liu et al. [14] proposed the concepts of the additive effect due to the substitution of the nonrisk (or risk) haplotype by the risk (or non-risk) haplotype and the dominance effect due to the interaction between the risk and non-risk haplotypes. These concepts have been integrated into the current model that allows the test of additive by sex and dominance by sex interaction effects. Simulation studies suggest that the new model displays adequate power to detect differences in these quantitative genetic parameters between the sexes, although the detection of sex-specific dominance effects needs a much larger sample size and/or heritability level.
Our model was used to analyze a real data set for pain genetics. Our analyses of a pain sensitivity trait-baseline pressure pain threshold measured at the ulna-by estimating genetic parameters and testing their sex-specific differences in a combination of male and female samples revealed that males and females have different population structure at two SNPs genotyped from a candidate gene (delta opioid receptor) for human pain and that haplo- types exert different genetic effects on the trait between the two sexes. A further test indicates that the risk haplotype [TT] detected by the model exemplifies sex-specific modes of inheritance in affecting the pain trait. While there are no differences among composite diplotypes in males, a significant additive effect was detected in females. Both additive and dominance effects due to the risk haplotype identified are different between the two sexes. Anholt and Mackay [2] described three major mechanisms that explain sex-specific difference in trait control, i.e., sex-specific effects (a gene affects only one sex), sex-biased effects (a gene affects both sexes but to different degrees), and sex-antagonistic effects (a gene affects both sexes but in opposite directions). In our example, the functional haplotype detected affects the pain trait in a sex-antagonistic effect manner, a mechanism thought to help the maintenance of genetic variation in natural populations [26].
In practice, failure to model sex-specific architecture may significantly hamper the ability to detect signals of functional genetic variants in genomewide screens. Although combining male and female data to increase sample size are tempting approaches to increase power, the estimates in this way will be biased from true sex-specific differences. Possible mechanisms that cause sex differences include parent-of-origin effects [27], linkage to or interaction with sex chromosomes, or differences arising from sex-specific hormonal environments. Our gene by sex interaction model that is incorporated by these mechanisms can be modified to consider interactions between genes and any other environments such as life style. Our interaction models should provide a more powerful tool to draw a detailed and precise picture of the genetic architecture of any complex traits that are important to human health.