A statistical measure for the skewness of X chromosome inactivation for quantitative traits and its application to the MCTFR data

Background X chromosome inactivation (XCI) is that one of two chromosomes in mammalian females is silenced during early development of embryos. There has been a statistical measure for the degree of the skewness of XCI for qualitative traits. However, no method is available for such task at quantitative trait loci. Results In this article, we extend the existing statistical measure for the skewness of XCI for qualitative traits, and the likelihood ratio, Fieller’s and delta methods for constructing the corresponding confidence intervals, and make them accommodate quantitative traits. The proposed measure is a ratio of two linear regression coefficients when association exists. Noting that XCI may cause variance heterogeneity of the traits across different genotypes in females, we obtain the point estimate and confidence intervals of the measure by incorporating such information. The hypothesis testing of the proposed methods is also investigated. We conduct extensive simulation studies to assess the performance of the proposed methods. Simulation results demonstrate that the median of the point estimates of the measure is very close to the pre-specified true value. The likelihood ratio and Fieller’s methods control the size well, and have the similar test power and accurate coverage probability, which perform better than the delta method. So far, we are not aware of any association study for the X-chromosomal loci in the Minnesota Center for Twin and Family Research data. So, we apply our proposed methods to these data for their practical use and find that only the rs792959 locus, which is simultaneously associated with the illicit drug composite score and behavioral disinhibition composite score, may undergo XCI skewing. However, this needs to be confirmed by molecular genetics. Conclusions We recommend the Fieller’s method in practical use because it is a non-iterative procedure and has the similar performance to the likelihood ratio method. Supplementary Information The online version contains supplementary material available at 10.1186/s12863-021-00978-z.


Background
In genome-wide association study (GWAS), many human diseases have been found to be associated with X-chromosomal genes, such as autoimmune diseases [1,2], asthma [3], Duchenne muscular dystrophy [4,5], adrenoleukodystrophy [6], Wiskott-Aldrich syndrome [7] and some cancers [8][9][10][11][12]. However, development of methods for identifying association with genetic variants on X chromosome still lags behind that on autosomes due to the unique inheritance pattern of X chromosome [13]. The number of X chromosomes is different between males and females in mammals. There are two copies of X chromosome in mammalian females, one of which is paternal and the other is maternal, while mammalian males have only one maternal X chromosome. To compensate for this X chromosome dosage difference between sexes, one of two chromosomes in females is silenced during the early development of embryos, which is called X chromosome inactivation (XCI) [14][15][16][17][18]. Random XCI (XCI-R) is a process that either the paternal or the maternal allele at an X-chromosomal locus is randomly chosen to be silenced in all cells, which is common in most females [19]. However, skewed XCI (XCI-S) is also observed in a proportion of females, which is a non-random process and is defined as the observation of inactivation of the same allele in more than 75% cells [9,[20][21][22][23]. In addition, not all of the X-linked genes undergo XCI and the pseudo-autosomal region on both sex chromosomes does not require dosage compensation. In humans, over 15% X-linked genes have been shown to escape from XCI (XCI-E) [24,25].
In population genetics, there has been an increasing interest in the incorporation of the information on XCI into association analysis for qualitative traits [26][27][28][29][30] and quantitative traits [31][32][33][34], which may greatly improve the test power. For qualitative traits, Clayton [26] first took account of XCI in detecting the association between X-chromosomal markers and diseases by regarding males as homozygous females. However, the Clayton's method only considers the XCI-R pattern and does not incorporate the XCI-S and XCI-E patterns. So, Wang et al. [27] developed a resampling-based approach for case-control data simultaneously combining the information on three XCI patterns (XCI-R, XCI-S and XCI-E) by coding three genotypes in females as 0, γ and 2, where γ is an unknown parameter, takes possible values between 0 and 2, and can be used to measure the degree of the skewness of XCI. For X-linked quantitative trait loci (QTL), Zhang et al. [31] proposed a familybased association test, where the quantitative trait under study is required to follow a normal distribution. Although the involved variances of the trait value for males and females are assumed to be different, those for three genotypes in females are fixed to be the same. However, according to Ma et al. [32], XCI may lead to variance heterogeneity of the traits across different genotypes in females and the variance of the trait in heterozygous females is generally higher than that in homozygous females. So, based on only unrelated females, Ma et al. [32] suggested a test for X-linked association via inflated variance in heterozygous females, a weighted test for Xlinked association which considers different variances, and the combined test of these two tests using the Stouffer's Z-score method. Gao et al. [33] further developed the XWAS software toolset to facilitate GWAS on X chromosome, which includes the three test statistics proposed by Ma et al. [32]. Deng et al. [34] put forward a sex-specific Levene's test, and a generalized Levene's test based on a two-stage regression model accounting for sex-specific mean and variance effects, to test for association. The original Levene's test is robust to certain types of non-normal distribution, particularly when data are non-normal but symmetric [34], while the generalized Levene's test may not. It should be noted that the above methods for QTL only incorporate the XCI-R and XCI-E patterns and do not consider the XCI-S pattern. On the other hand, Wang et al. [35] has recently proposed a statistical measure available for the degree of XCI skewing for case-control data and developed three methods (likelihood ratio (LR), Fieller's and delta) to construct the corresponding confidence intervals (CIs). However, they are only applicable to qualitative traits and are not suitable for quantitative traits.
Therefore, in this article, we first extend the existing statistical measure for the degree of XCI skewing (i.e., γ) for qualitative traits [35] and make it accommodate quantitative traits. It is shown that the proposed γ is a ratio of two linear regression coefficients in the presence of association between the traits under study and the genotypes. We estimate the linear regression coefficients by incorporating the information on the variance heterogeneity across different genotypes in females and then obtain the point estimate of γ. Then, we extend the existing LR, Fieller's and delta methods for constructing the CIs of γ and make them suitable for quantitative traits. The simulation studies under various simulation settings are conducted to assess the performance of the proposed methods. We also apply the proposed methods to the Minnesota Center for Twin and Family Research (MCTFR) data for their practical use. Note that so far, we are not aware of any association study for the X-chromosomal markers in the MCTFR data, although there have been some previous association studies which only focused on autosomal markers [36][37][38][39][40][41][42][43].

Sizes and powers
The empirical type I error rates of the corresponding tests for the proposed LR, Fieller's and delta methods based on the sample size n = 1,000 and 2,000 are respectively given in Tables 1 and 2, where the additive effect size a = 0.1 and 0.3, the allele frequency p = 0.1 and 0.3, and the inbreeding coefficient ρ = 0. Under all the situations considered, the sizes of the proposed LR and Fieller's methods stay close to the pre-specified nominal level of 5%, irrespective of the values of n, a and p, which verifies their validity. However, the delta method has the inflated or conservative type I error rates in most scenarios. Additional file 1: Tables S1 and S2 show the sizes for the proposed LR, Fieller's and delta methods with ρ = 0.05 based on the sample size n = 1,000 and 2, 000, respectively, which are similar to those in Tables 1  and 2. This demonstrates that the Hardy-Weinberg disequilibrium almost has no effect on the sizes.
Note that the delta method does not control the sizes well. So, we only simulate the powers of the LR and Fieller's methods. Figures 1, 2 and 3 display the estimated powers for the LR and Fieller's methods against γ (γ ≠ γ 0 ) with a = 0.1 and 0.3, p = 0.1 and 0.3, n = 1,000, and ρ = 0 when γ 0 = 0, 1 and 2, respectively. Figures 4, 5 and 6 plot the corresponding estimated powers with a = 0.1 and 0.3, p = 0.1 and 0.3, n = 2,000, and ρ = 0 when γ 0 = 0, 1 and 2, respectively. The other power results are shown in Additional file 1: Figures S1-S14. It can be seen from these figures that the power of the LR method is almost the same as that of the Fieller's method. The powers of the LR and Fieller's methods gradually but asymmetrically become larger with |γ − γ 0 | increasing. When other parameters are unchanged, the powers with p = 0.3 are bigger than those with p = 0.1 (e.g., Fig. 1b vs. Fig. 1a, Fig. 1d vs. Fig. 1c). However, note that in σ 2 1 ¼ θð1−θÞa 2 þ 1:1, θ(1 − θ)a 2 attains its maximum 0.25a 2 when θ = 0.5 (i.e., γ = 1). The corresponding values of σ 2 1 for a = 0.1 and 0.3 are 1.1025 and 1.1225, respectively, which are not so different from each other. Furthermore, when γ = 0 or 2, θ(1 − θ)a 2 = 0, which is not related to the value of a. So, the powers with a = 0.1 and those with a = 0.3 are close to each other (e.g., Fig. 1a vs. Fig. 1c, Fig. 1b vs. Fig. 1d). When the sample size n is changed from 1,000 to 2,000, the LR and Fieller's methods are more powerful (e.g., Fig. 4 vs. Fig. 1). Finally, we find that the Hardy-Weinberg disequilibrium has little influence on the power results, e.g., by comparing Fig. 1 (ρ = 0) with Additional file 1: Figure S3 (ρ = 0.05).  Median of point estimate and statistical properties of confidence intervals Tables 3 and 4 show the estimated median of the point estimates of γ, CP, ML, MR, ML/(ML + MR), DP and EP of the two-sided 95% CIs of γ for the LR, Fieller's and delta methods against γ, with a = 0.1 and 0.3, p = 0.1 and 0.3, and ρ = 0 based on 10,000 replicates for n = 1,000 and 2,000, respectively. From these two tables, we find that in all the cases considered, the median ofγ maintains very close to the true value of γ. As for the CI, the LR and Fieller's methods have similar performance in the CP and the CPs of both methods are controlled around 95%, regardless of the values of a, p, γ and n. However, the CP of the delta method is underestimated or overestimated in most of the considered situations. The values of the ML/(ML + MR) for the LR and Fieller's methods generally stay close to 0.5, except for the cases of p = 0.1 and n = 1,000, and the situations of γ = 0 and 2, while the ML/(ML + MR) of the delta method always gets far way from 0.5. This indicates that the LR and Fieller's methods achieve more balance between ML and MR than the delta method. The LR and Fieller's methods have comparable performance in the DP and EP. The values of the DP of both methods are zero or close to zero, except for p = 0.1 and γ = 0, which is indicative of few discontinuous CIs to occur. However, the EP results of the LR and Fieller's methods show that there still exist a few CIs which are empty sets or reduced to be a point. On the other hand, the DP and EP of the delta method are zero for all the simulation settings. This is because the CI based on the delta method is always bounded and is a continuous interval. The ML/ (ML + MR), DP and EP of the LR and Fieller's methods appear not to be greatly affected by the values of a (0.1 or 0.3), while the LR and Fieller's methods perform worse in the ML/(ML + MR) and the DP when p = 0.1, compared to those with p = 0.3. When the sample size increases from 1,000 (Table 3) to 2,000 (Table 4)  of the EP with n = 2,000 are a little bigger than those with n = 1,000. This may be because γ= 0 and 2 are the endpoints of the interval [0, 2], which are the extreme cases. The corresponding results of the median ofγ, CP, ML, MR, ML/(ML + MR), DP and EP of the 95% CIs of γ for the LR, Fieller's and delta methods with ρ = 0.05 for n = 1,000 and 2,000 are given in Additional file 1: Tables S3 and S4, respectively. By comparing Table 3 with Additional file 1: Table S3 (or comparing Table 4 with Additional file 1: Table S4), we can see that the results in both tables are similar to each other, which means that the Hardy-Weinberg disequilibrium has no great effect on the point estimation and the interval estimation of γ.

Application to MCTFR data
We applied the proposed LR, Fieller's and delta methods to the data from the MCTFR GWAS of Behavioral Disinhibition for their practical use, and considered the following five quantitative traits: the nicotine composite score, alcohol consumption composite score, alcohol dependence composite score (DEP), illicit drug composite score (DRG) and behavioral disinhibition composite score (BD). The MCTFR data are made available for download from the database of Genotypes and Phenotypes (accession number: phs000620.v1.p1). In the MCTFR data, there are 2183 families (7377 subjects consisting of 3546 males and 3831 females), including 182 families with 1 member (182 subjects), 290 families with 2 members (580 subjects), 294 families with 3 members (882 subjects), 1352 families with 4 members (5408 subjects), and 65 families with 5 members (325 subjects). Among them, nuclear families are composed of the parents and two offspring who are monozygotic twins, full biological non-twin siblings, adopted siblings and mixed siblings with 1 biological offspring and 1 adopted offspring. Figure 7 shows more details of the family structure in the MCTFR data. Twelve thousand three hundred fifty-four single nucleotide polymorphisms (SNPs) on the X chromosome were genotyped. Note that our proposed methods are applicable in the presence of association between the SNPs and the quantitative traits of interest. So, we first conducted the association analysis for each locus and each trait. When only analyzing a single trait for all the 12,354 SNPs, the significance level was set to be α ′ = 0.05/12,354 = 4.047 × 10 −6 based on Bonferroni correction. When simultaneously analyzing multiple traits, Deng et al. [34] and McGue et al. [41] fixed the significance level at 1 × 10 −3 for their association analysis. Therefore, in this application, we also used this significance level for the association study when simultaneously considering multiple traits. Then, we calculated the point estimate and the corresponding CIs of the skewness of XCI at the 95% confidence level for all the SNPs which are associated with a single trait at the 4.047 × 10 −6 level or are simultaneously associated with two or more traits at the 1 × 10 −3 level. However, we found that all these traits, and the transformed traits (e.g., log(y + 1)) do not satisfy the normality assumption. As such, we used the existing Levene's test [34] to detect the association between the SNPs and the traits, which is robust to certain types of non-normal distribution.
The following quality control rules are used to filter the data. First, note that the proposed three methods for the interval estimation of γ only utilize unrelated females. On the other hand, although the adopted offspring in the nuclear families are biologically independent of their adopted parents, they might come from a subpopulation which is different from that of their parents. So, we deleted all the males in the data and all the offspring in the nuclear families, including the biological offspring and the adopted offspring. Second, genotyped female individuals with missing genotype rate over 10% were excluded. Third, the SNPs with missing genotype rate over 10% were deleted. Finally, we applied the PLINK software to carry out the HWE tests for SNPs [44] and the significance level is set to be 1 × 10 −4 [45]. The SNPs with the minor allele frequency (MAF) less than 5% or those out of HWE were also excluded. As such, a total of 1955 unrelated females and 11,355 SNPs were included in this application.
The Levene's test identified one SNP (rs17261621) which is only associated with the DRG trait at the 4.047 × 10 −6 level, two SNPs (rs792959 and rs17261621) which are associated with both the DRG and BD traits and three SNPs (rs4825722, rs4825726 and rs2196260) which are associated with both the DEP and BD traits at the 1 × 10 −3 level. The corresponding P values of the Levene's test and the HWE test together with the position, the MAF, the point estimates and the CIs of γ based on the LR, Fieller's and delta methods for these five SNPs are given in Table 5. For the DRG trait and the rs792959 locus, the point estimate of γ, and the 95% CIs of the LR, Fieller's and delta methods are 2, (1.0294, 2], (1.0293, 2] and [0, 2], respectively. For the BD trait and the rs792959 locus, the point estimate of γ, and the corresponding 95% CIs are 2, (1.0306, 2], (1.0304, 2] and [0, 2], respectively. The CIs of the LR and Fieller's methods for the DRG and BD traits are very similar and do not contain 1. Thus,γ being 2 indicates that at rs792959, 100% (2/2) of cells in heterozygous females have allele G active, and 0% of cells express allele A, which demonstrates the XCI-S towards allele G. However, the CIs of the delta method at rs792959 contain 1 (i.e., XCI-R). The conclusions drawn from the LR and Fieller's methods here are similar to those drawn from our simulation study. However, the truncated point estimateγ is 2, which is the right endpoint of the interval [0, 2]. This may be because the proposed LR and Fieller's methods require that the traits under study follow a normal distribution, while the DRG and BD traits are not normally distributed. Further, all the CIs for the other four SNPs contain 1, indicating random XCI. Particularly, for the BD trait and the rs4825722 locus, the CIs of the LR, Fieller's and delta methods are [0, 2], which provides no information on the XCI pattern.

Discussion
In this article, we extended the existing statistical measure for the degree of XCI skewing (i.e., γ) and the existing LR, Fieller's and delta methods for constructing the CIs of γ for qualitative traits [35], and made them suitable for quantitative traits. The proposed γ is a ratio of two linear regression coefficients in the presence of association between the traits under study and the genotypes. According to Ma et al. [32], XCI may cause variance heterogeneity of the traits across different genotypes in females. As such, we estimated the linear regression coefficients by incorporating the information on the variance heterogeneity and then obtained the point estimate of γ. The Fieller's and delta methods for calculating the CIs are simple and non-iterative procedures, while the LR method is an iterative one which needs more computing time. On the other hand, the hypothesis testing of the LR, Fieller's and delta methods was also investigated. We conducted extensive simulation studies (two different values of additive effect, two groups of allele Nð0; 1Þ is not so strictly correct. An-other possible reason why the delta method performs so poorly is that the first order Taylor expansion ofγ does not suffice. To investigate the performance of the delta method with higher order Taylor expansion, we used the second order Taylor expansion ofγ to calculate the asymptotic variance ofγ , which can be implemented in the "propagate" package in R software [46]. However, most of the estimated type I error rates for the delta method are still inflated or conservative, even though they appear to be controlled better than those in Tables 1 and 2 (data not shown for brevity). Therefore, in practical use, we recommend the Fieller's method because it is a non-iterative procedure and has the similar performance to the LR method. So far, we are not aware of any association study for the X-chromosomal SNPs in the MCTFR data. In fact, we also found that all the five traits in the MCTFR data are not normally distributed. On the other hand, when simultaneously analyzing multiple traits for the Xchromosomal SNPs, Deng et al. [34] fixed the significance level at 1 × 10 −3 for their association analysis. So, in our real data application, we used the existing Levene's test [34] to test for the association between the X-chromosomal SNPs and the five traits at the significance level of 1 × 10 −3 , which does not require the normality assumption for the traits. However, when only analyzing a single trait for all the 12,354 SNPs, the significance level is set to be α ′ = 0.05/12,354 = 4.047 × 10 −6 based on Bonferroni correction. One SNP (rs17261621) is shown to be only associated with the DRG trait at the 4.047 × 10 −6 level, two SNPs (rs792959 and rs17261621) are identified to be associated with both the DRG and BD traits, and three SNPs (rs4825722, rs4825726 and rs2196260) are found to be associated with both the DEP and BD traits at the 1 × 10 −3 level. In addition, we applied the proposed LR, Fieller's and delta methods to these five SNPs and calculated the CIs of the skewness of XCI at the 95% confidence level. The CIs based on the LR and Fieller's methods show that only rs792959 undergoes XCI-S. However, these conclusions need to be further confirmed by molecular genetics. On the other hand, the proposed LR and Fieller's methods require that the traits under study follow a normal distribution, while the DEP, DRG and BD traits are not normally distributed. Since we have no suitable data of this kind available, it is of future interest to apply the three proposed methods to datasets with traits following normal distributions and to further confirm their practical use.
Besides, the proposed methods have the following issues to discuss. First, to make the point estimate and the CIs of γ more interpretable, we simply use the interval [0, 2] to truncate the original point estimate and the original CIs, which may cause potential loss of information, and may also lead to the truncated CIs being empty sets when the original CIs lie outside [0, 2]. Fortunately, from our simulation study, the proportion of the CIs being empty sets or being reduced to be a point among all the simulation replications is all less than 2.7%. On the other hand, to incorporate the interval constraint of [0, 2] into statistical inference, we will develop a future Bayesian method to estimate the skewness of XCI by considering such constraint as prior information. Second, the proposed methods require the association between the traits and the SNPs being present. As such, in genome-wide association study, we could regard the screening of the associated SNPs as a preliminary step     before estimating the skewness of XCI. If such association is not statistically significant, the LR and Fieller's methods may result in the discontinuous CIs, which is difficult to interpret. Third, the normality assumption of the traits under study is needed in the proposed methods. In future, we will extend them to accommodate the traits not normally distributed. Finally, the proposed methods are only applicable to unrelated female subjects. Thus, we will extend the proposed methods and make them suitable for data with family or pedigree structure in future studies.

Conclusions
We recommend the Fieller's method in practical use because it is a non-iterative procedure and almost has the  same performance as the LR method. On the other hand, only rs792959, which is identified to be associated with both the DRG and BD traits, may undergo XCI-S, which needs to be confirmed by molecular genetics.

Notations and point estimate of γ
Consider an X-linked diallelic QTL. Let D and d represent the mutant allele and the normal allele at the QTL with the frequencies being p and q (p + q = 1), respectively. Note that XCI is unrelated to males and we only consider females here. Then, there are three possible genotypes dd, Dd and DD at the QTL in females. The corresponding frequencies are g 0 = q 2 + ρpq, g 1 = 2(1 − ρ)pq and g 2 = p 2 + ρpq, respectively, where ρ is the inbreeding coefficient. ρ = 0 means that Hardy-Weinberg equilibrium (HWE) holds in the female population, while ρ ≠ 0 denotes Hardy-Weinberg disequilibrium. Suppose that Y and G are the value of the quantitative trait under study and the genotype of a female subject, respectively. Notice that XCI may lead to variance heterogeneity of Y across different genotypes [32]. So, we assume that Y j G¼dd Nðμ 0 ; σ 2 0 Þ , Y j G¼Dd Nðμ 1 ; σ 2 1 Þ and Y j G¼DD Nðμ 2 ; σ 2 2 Þ: Further, let X 1 = I {G = Dd or DD} and X 2 = I {G = DD} . As such, X 1 denotes that this female carries at least one mutant allele D and X 2 indicates that she is a homozygote DD. Then, to construct the statistical measure of the skewness of XCI, we consider the following linear regression model where Z is a vector of covariates which need to be adjusted, β 0 is the intercept, β 1 and β 2 respectively are the regression coefficients of X 1 and X 2 , and b is a vector of the regression coefficients for Z. From Eq. (1), Under XCI-R, μ 1 should lie midway between μ 0 and μ 2 , which is β 0 þ Hence, for heterozygous females, any statistically significant deviation from such value can be regarded as an evidence of XCI-S. This is equivalent to that is far away from 0.5. Therefore, we define the following parameter γ to measure the skewness of XCI with β 1 + β 2 ≠ 0. And θ = γ/2, on the average, is indicative of the proportion of cells in a Dd female keeping the mutant allele D active. Thus, γ = 1 denotes XCI-R. 1 < γ ≤ 2 means the XCI-S towards D and 0 ≤ γ < 1 represents the XCI-S against D. For example, when γ = 1.5, then θ = 75%, which means that 75% cells have mutant allele D active and the other 25% cells express the normal allele d. Let β = (β 1 + β 2 )/2, then γ = β 1 /β. In this regard, we obtain β 1 = γβ and β 2 = (2 − γ)β, where β ≠ 0. Let X = γX 1 + (2 − γ)X 2 , then Eq. (1) becomes Here, the genotypic value X equals 0, γ and 2 for genotypes dd, Dd and DD, respectively, which implies that the definition of γ coincides with the coding strategy of Wang et al. [27] for XCI. On the other hand, from Eq. (2), we observe that γ can be well defined when the association between Y and the allele of interest is present (i.e., β ≠ 0).
Assume that we collect a sample of n independent females. Let n 0 , n 1 and n 2 be the number of females with genotypes dd, Dd and DD, respectively. So, n 0 + n 1 + n 2 = n. Let y ij , x ij1 , x ij2 and z ij denote the values of Y, X 1 , X 2 and Z of female j (j = 1, 2, ⋯, n i ), where i = 0, 1, 2 respectively correspond to genotypes dd, Dd and DD. According to Eq. (1), the log-likelihood function of the sample is Then, by maximizing the above equation, the maximum likelihood estimatesβ 0 ;β 1 ;β 2 ;σ 0 ;σ 1 ;σ 2 andb of β 0 , β 1 , β 2 , σ 0 , σ 1 , σ 2 and b can be derived. As such, according to Eq. (2), the point estimate of γ can be given as 2]. Then, the final point estimate of γ is cut off by 0 and 2, and denoted byγ.
Sinceγ is a ratio estimate, borrowing the idea of Wang et al. [35], we find that the standard error ofγ can be approximated by using the delta method. Specifically, take a first order Taylor expansion of γ around the point (β 1 , β) and evaluate it at ðβ 1 ;βÞ , which yieldsγ ≈ β 2 Þ are the elements of the variance-covariance matrix of β 1 andβ 2 , which can be computed by using the "glm" function in R software [46]. Respctively replacing β 1 and β byβ 1 andβ in Eq. (5) [35]. Under H 0 : γ = γ 0 , we have β 1 − γ 0 β = 0. Then, we can construct the following Wald test for testing H 0 : γ = γ 0 The confidence limits of the 100(1 − α)% CI based on the j > z α=2 are equivalent to each other, which mean that the association is present at the significance level of α. In addition, when Δ = 0 or A = 0, the CI is reduced to a point. As such, the 100(1 − α)% CI based on the Fieller's method is if Δ < 0 and A < 0 It should be noted that if Δ > 0, the above CI may be an empty set. When Δ > 0 and A < 0, the corresponding CI may be the union of two disjoint intervals, which is the discontinuous CI.

Simulation settings
For simplicity, we do not include any covariate in the model. The frequency p of allele D at the locus on X chromosome is fixed at 0.1 and 0.3. The inbreeding coefficient ρ is taken as 0 and 0.05 to respectively simulate the situation of HWE and that of Hardy-Weinberg disequilibrium. Let β 0 = 0.1 and β = 0.3. Further, β 1 = γβ and β 2 = (2 − γ)β are calculated from β = 0.3 and γ, where γ takes values of 0, 0.5, 1, 1.5 and 2. γ 0 is also assigned to be 0, 0.5, 1, 1.5 and 2 with γ = γ 0 to simulate the type I error rates of the proposed methods and γ ≠ γ 0 for simulating the their test powers. As mentioned in Ma et al. [32], the variance σ 2 1 of the trait value for heterozygous females is generally larger than σ 2 0 and σ 2 2 for homozygous females due to XCI. So, we set σ 2 0 ¼ σ 2 2 ¼ 1, and σ 2 1 ¼ θð1−θÞa 2 þ 1:1, where a is the additive effect of the QTL, θ = γ/2 is the inactivation ratio as mentioned before, and the variance caused by other factors is fixed to be 1.1. Here, a is set to be 0.1 and 0.3. The sample size n is selected to be 1,000 and 2,000. The genotype of each female is simulated according to the allele frequency p and the inbreeding coefficient ρ. Then, the trait value Y of this female given her genotype is generated by Y j G¼dd Nðβ 0 ; σ 2 0 Þ , Y j G¼Dd Nðβ 0 þ β 1 ; σ 2 1 Þ or Y j G¼DD Nðβ 0 þ β 1 þβ 2 ; σ 2 2 Þ. For each simulation setting, the simulations are conducted based on K = 10,000 replications and the significance level α is fixed at 5%. The simulation study is implemented in R software (version 3.2.5) [46].
Notice that the distribution of the point estimateγ may be asymmetric. So, we list the median ofγ 's over K replications to describe the central tendency of this skewed distribution. We assess the statistical properties of the CIs of γ by the following indexes. The coverage probability (CP) is the proportion that the CIs contain the true value γ among K replications, irrespective of the CI being continuous or discontinuous. DP and EP are the proportion of the discontinuous CIs and that of the CIs being an empty set or being reduced to be a point among K replications, respectively. Simulation study is also carried out to investigate the probabilities of the CI missing the true value γ on the left (ML) and on the right (MR), and the value of the ratio ML/(ML + MR), which is close to 0.5 when the balance between ML and MR is achieved. Here, ML ¼ #½ðγ<γ L Þ∩ðγ L ≤γ ≤ γ U Þ K and MR where # is the counting measure and (γ L , γ U ) 's are the continuous CIs. We only consider the continuous CIs when computing the ML and MR, because we cannot distinguish between the left side and the right side of the discontinuous CIs. SNP: Single nucleotide polymorphism; XCI: X chromosome inactivation; XCI-E: Escape from X chromosome inactivation; XCI-R: Random X chromosome inactivation; XCI-S: Skewed X chromosome inactivation