 Methodology article
 Open Access
 Published:
A statistical measure for the skewness of X chromosome inactivation based on casecontrol design
BMC Bioinformatics volume 20, Article number: 11 (2019)
Abstract
Background
Skewed X chromosome inactivation (XCI), which is a nonrandom process, is frequently observed in both healthy and affected females. Furthermore, skewed XCI has been reported to be related to many Xlinked diseases. However, no statistical method is available in the literature to measure the degree of the skewness of XCI for casecontrol design. Therefore, it is necessary to develop methods for such a task.
Results
In this article, we first proposed a statistical measure for the degree of XCI skewing by using a casecontrol design, which is a ratio of two logistic regression coefficients after a simple reparameterization. Based on the point estimate of the ratio, we further developed three types of confidence intervals (the likelihood ratio, Fieller’s and delta methods) to evaluate its variation. Simulation results demonstrated that the likelihood ratio method and the Fieller’s method have more accurate coverage probability and more balanced tail errors than the delta method. We also applied these proposed methods to analyze the Graves’ disease data for their practical use and found that rs3827440 probably undergoes a skewed XCI pattern with 68.7% of cells in heterozygous females having the risk allele T active, while the other 31.3% of cells keeping the normal allele C active.
Conclusions
For practical application, we suggest using the Fieller’s method in large samples due to the noniterative computation procedure and using the LR method otherwise for its robustness despite its slightly heavy computational burden.
Background
X chromosome inactivation (XCI) is an epigenetic phenomenon. Under XCI, one of two X chromosomes in females is silenced during early embryonic development to achieve dosage compensation between two sexes [1]. As such, the genetic effect of two risk alleles in females is expected to be equivalent to that of one risk allele in males. Most of Xlinked genes undergo XCI and only about 15% of genes on X chromosome escape from XCI (XCIE) [2]. Both alleles in the genes under XCIE will be active, which are similar to autosomal genes. Generally, XCI has been treated as random (XCIR) where both maternal and paternal X chromosomes have equal chance to be inactivated, i.e. for an Xlinked gene, nearly 50% of the cells have one allele active while the remaining cells have the other allele active. However, recent studies have revealed that the skewed XCI (XCIS) is a biological plausibility and even a common feature in both healthy and affected females [3–5]. XCIS is a nonrandom process, which has been defined as a significant deviation from XCIR, for instance, the inactivation of one of the alleles in more than 75% of cells [6–8].
The mechanism of XCIS remains mysterious and XCIS in human may be likely caused by secondary selection [6, 9, 10]. Specifically, the initial choice of active X chromosome is considered as random. However, during body growth, when an Xlinked mutation affects cells proliferation or survival, there will be a larger or smaller proportion of cells with an active mutant allele. Due to the selection pressure, this type of secondary skewing varies in different tissues and is also associated with age. For heterozygous females, positive selection cells with mutant allele will lead to more severe expression of the disease, whereas negative selection cells with mutant allele can provide protection from deleterious effects [11, 12]. For example, in heterozygous females with a mutant FoxP3 allele, the XCIS against the mutant allele in specific tissues can prevent autoimmune disease, whereas the XCIS towards the mutant allele in breast epithelial cells can result in breast cancer [13]. On the other hand, some diseases, such as ovarian cancer, Rett syndrome, Klinefelter syndrome, and recurrent miscarriages, are reported to be related to XCIS [14–17]. Therefore, it is necessary to develop methods for measuring such XCI skewing.
Recently, there has been an increasing interest to incorporate the information on XCI into Xchromosome genetic association studies [18–23]. Clayton’s method first takes XCIR into account and treats males as homozygous females [18]. In this regard, two genotypes of males are coded as 0 or 2, while three genotypes of females are coded as 0, 1 or 2, respectively. This coding strategy also implies that the genetic effect of heterozygous genotype in females lies midway between two homozygous genotypes, which seems reasonable as in heterozygous females about half of cells express the mutant allele while the rest of cells express the normal allele. However, this method does not consider the XCIE and XCIS patterns. So, a resamplingbased method was proposed by maximizing the likelihood ratio (LR) over all the three biological patterns (XCIE, XCIR and XCIS), where the three genotypes of females are coded as 0, γ or 2 under XCIS [21]. Note that γ is an unknown parameter which is used to measure the degree of XCI skewing. For instance, γ=1 represents XCIR; γ=1.5 indicates XCIS where 75% of the cells have the mutant allele active, whereas the other 25% of the cells have the normal allele active. On the other hand, the detection of XCIS is either by measuring the level of methylation or by integrative analysis of whole exome and RNA sequencing data [24, 25]. Although Xu et al. has recently developed a statistical measure for the skewness of XCI based on family trios [26], there is still no statistical method available in the literature to measure the skewness of XCI for casecontrol design.
Therefore, in this article, we first showed that γ can be represented as a ratio of two logistic regression coefficients after a simple reparameterization, based on casecontrol data. We then obtained the point estimate of γ by the maximum likelihood estimates (MLEs) of these two regression coefficients. Further, we derived the confidence interval (CI) of γ by the delta method, the Fieller’s method and the LR method. We also applied all the proposed approaches to analyze the Graves’ disease data for their practical use.
Results
Statistical properties of confidence interval
Tables 1 and 2 list the estimated coverage probability (CP), left tail error (ML), right tail error (MR) (missing the true value of γ), ML/(ML+MR) and proportion of the discontinuous CIs (DP) of the LR, Fieller’s and delta methods under various simulation settings with N=500, ρ=0, and λ_{2}=1.5 and 2, respectively. From the tables, the LR and Fieller’s methods control the CP well when p=0.3. However, when p=0.1, both the LR and Fieller’s methods appear to overestimate the CP except for γ=0. Note that the CPs of the LR method are closer to the preset level than the Fieller’s method, which indicates the robustness property of the LR method for relatively small samples. Besides, the delta method generally has the worst CP under all the situations. When p=0.1, the delta method overestimates the CP for γ=0, while underestimates the CP for γ=1, 1.5 and 2, irrespective of λ_{2} being 1.5 or 2. When p increases from 0.1 to 0.3, the delta method overestimates the CP for γ=0, 0.5 and 1, while underestimates the CP for γ=1.5 and 2, regardless of λ_{2}=1.5 or 2. From the estimated ML, MR and ML/(ML+MR) values, we find that ML and MR of the deltatype CIs are not balanced since nearly all the values of ML/(ML+MR) are far away from 0.5, while the LR and Fieller’s methods have more balanced ML and MR than the delta method, especially when p=0.3. On the other hand, we see that the values of the DP for both the LR and Fieller’s methods are not over 3% under our simulation settings. Further, when p increases from 0.1 to 0.3 or λ_{2} changes from 1.5 to 2, DP will generally become smaller.
Tables 3 and 4 give the estimated CP, ML, MR, ML/(ML+MR) and DP of the LR, Fieller’s and delta methods when N=2000, ρ=0, and λ_{2}=1.5 and 2, respectively. When N increases from 500 to 2000, the LR method has similar performance with the Fieller’s method. It can be seen from the tables that the CPs of all the methods are more accurate. Both the LR and Fieller’s methods control the CP well, while the delta method still generally has the poor CP, especially when p=0.1. Note that when p=0.3, all the values of ML/(ML+MR) for the LR method and the Fieller’s method are around 0.5 regardless of λ_{2} being 1.5 or 2. But when p=0.1, the values of ML/(ML+MR) for the LR method and the Fieller’s method are deviated from 0.5, especially when γ_{0} are 0 and 2. Further, the Fieller’s method has slightly more balanced tail errors than the LR method when p=0.1. In addition, the delta method has the most unbalanced tail errors. We also find that the values of the DP generally decrease to be less than 2% when N increases to be 2000. In general, the LR method and the Fieller’s method control the CP well with the relatively balanced tail errors on the left and on the right. All the other results of CP, ML, MR, ML/(ML+MR) and DP with ρ=0.05 are given in Tables S1S4 [see Additional file 1], which are similar to those in Tables 1, 2, 3 and 4 except for N=500 and p=0.1, indicating that HardyWeinberg disequilibrium has limited effect on the results. Notice that under the scenario of N=500 and p=0.1, we observe that the CPs of all the methods in Additional file 1: Tables S1 and S2 are better than those in Tables 1 and 2, respectively. One possible explanation of this phenomenon is that the genotype frequency of AA in the control sample increases from 0.01 to 0.0145 when ρ changes from 0 to 0.05.
Sizes and powers
We also simulated the corresponding size and power for testing γ=γ_{0} [see Appendix B of Additional file 1]. The size results are given in Additional file 1: Tables S5–S8 and the power results are displayed in Figures S1S12 [see Additional file 1]. It can be seen that the LR method and the Fieller’s method control the size well except for N=500 and p=0.1, while the size of the delta method is either conservative or inflated. On the other hand, the power of the LR method and the Fieller’s method are close to each other, but the LR method is generally slightly more powerful than the Fieller’s method. However, the power of the delta method can be quite different from those of the LR and Fieller’s methods.
Application to Graves’ disease data
The GPR174 gene is located on X chromosome, which is associated with autoimmune thyroid disease, including Graves’ disease [27]. An X chromosome genomewide association study (GWAS) was conducted by Chu et al. to study the association between the GPR174 gene and Graves’ disease among Han population [27]. In this study, 14,141 single nucleotide polymorphisms on X chromosome were genotyped. Among them, rs3827440 is a nonsynonymous single nucleotide polymorphism within the GPR174 gene, with the minor allele frequency being 0.45 in this population, and thus is a functional variant of interest. Further, statistical analysis of both the GWAS data and the replication data showed that rs3827440 is statistically significantly associated with Graves’ disease. At rs3827440, there are two alleles T and C, where T is the susceptible allele which is associated with a higher expression level of the GPR174 gene. Several studies [7, 15] showed that XCIS is associated with autoimmune thyroid disease. So, we applied the LR, Fieller’s and delta methods to explore if rs3827440 undergoes an XCIS pattern. We only selected the females to estimate the degree of XCI skewing as well as its 95% CI. In the GWAS stage, 2242 females were sampled (1115 cases and 1127 controls). In the case group, the numbers of females with genotypes CC, TC, and TT are 163, 508, and 444, respectively. Those in the control group are 219, 541, and 367, respectively. In the replication stage, 6260 females were sampled with genotype counts 471, 1606, and 1298 in the case group, and 584, 1344, and 957 in the control group for CC, TC, and TT, respectively. The estimated allele frequency of T in females is 0.57 and 0.56 for the GWAS stage and the replication stage, respectively. We applied each of the proposed methods to the data in the GWAS stage and those in the replication stage. After that, we used our proposed methods to deal with the pooled data, by incorporating the stage as a covariate.
Table 5 gives the point estimate \(\hat {\gamma }\) and its 95% CIs, based on the LR, Fieller’s and delta methods. From the table, we observe that the LRtype CIs and the Fieller’s CIs are almost the same. The deltatype CIs are nested within the LRtype and Fieller’s CIs for the replication stage and the pooled data, which may be caused by the fact that the delta method underestimates the CP. We also find that the LRtype CI and the Fieller’s CI are asymmetrical around its point estimate in the GWAS stage but are nearly symmetrical around the point estimate in the replication stage and the pooled analysis, which is probably due to the larger sample size in the replication stage and the pooled dataset. In the GWAS stage, the point estimate \(\hat {\gamma }\) is 0.957. All of the three types of CIs contain 1 (XCIR). In the replication stage, the point estimate \(\hat {\gamma }\) is 1.513 and all the CIs do not contain 1. The results in the replication stage suggest the XCIS pattern at rs3827440 with 75.7% (1.513/2) of cells having the risk allele T active and the other 24.3% of cells having the normal allele C active. Note that the statistical results for both two stage data suggest different XCI patterns. One possible reason is that the variance of \(\hat {\gamma }\) is larger in the GWAS data and there may exist study heterogeneity between those two stages. The results for the pooled data give the point estimate \(\hat {\gamma }=1.373\) by adjusting the stage and all of the three types of CIs do not contain 1. This demonstrates that rs3827440 probably undergoes the XCIS pattern with 68.7% (1.373/2) of cells keeping the risk allele T active, while the other 31.3% of cells keeping the normal allele C active. However, this observation needs to be further confirmed by functional analysis of this variant.
Discussion
In this article, we proposed a statistical measure to estimate the degree of the skewness of XCI (i.e. γ). We first showed that γ can be expressed as a ratio of two logistic regression coefficients. Then, we constructed a ratio estimate \(\hat {\gamma }\) for γ and also derived three types of CIs (the LR, Fieller’s and delta methods) to evaluate its variation. The delta method is a simple and noniterative procedure but generally has poor statistical properties, which is probably caused by the skewness of \(\hat {\gamma }\). On the other hand, the LR method and the Fieller’s method are based on a simple reparameterization procedure and thus does not require the normality assumption of \(\hat {\gamma }\). The simulation results demonstrate that the LR method and the Fieller’s method have better performance than the delta method. On the other hand, note that the LRtype CI will be close to the Fieller’s CI when N is large. In this regard, the Fieller’s CI is preferential since it is a noniterative procedure. However, when N is relatively small, the LRtype CI is recommended for its robustness. In addition, our software SkewXCI is freely available at http://www.echobelt.org/web/UploadFiles/SkewXCI.html, which is implemented in R (http://www.rproject.org/, version 3.5.1).
Our proposed methods have several limitations. First, our methods assume that the genetic effect of the mutant allele among all the cells is additive on the disease. On the other hand, notice that Model (1) under XCI is different from genetic models (dominant, additive, and recessive) on autosomes or on X chromosome under XCIE. Specifically, genetic model defines the relationship between two alleles at a locus and usually varies from locus to locus [28–30]. However, when XCI occurs, only one allele is active at each locus in each cell and most of the loci share the same XCI pattern. As such, the magnitude of γ/2 is a measure of the proportion of cells with the mutant allele active among all the cells in heterozygous females. For instance, adrenoleukodystrophy has been previously viewed as an Xlinked recessive disorder where the female carriers are commonly thought to be normal or only mildly affected [31]. However, a recent study showed that the heterozygous females with adrenoleukodystrophy have a wide spectrum of clinical manifestations, ranging from mild to severe phenotypes, which is probably due to the various degree of XCIS towards the mutant allele. Second, we simply cut the estimated CI within the interval [ 0, 2] and this may lead to potential loss of information. However, if we incorporate this interval constraint into statistical inference, then the LR, Fieller’s and delta methods no longer follow a simple chisquare distribution or a standard normal distribution due to the boundary problem [32]. An alternative method is the Bayesian inference, where such constraint can be regarded as prior information. For instance, when no other information is available, we can choose an uniform prior distribution within the interval [ 0,2] for γ. Once the posterior distribution is derived, its percentiles or variance can be used to construct the corresponding CI. Third, note that the validity of our proposed measure is based on the assumption that there exists association between disease and allele A. Therefore, in GWAS, we can first screen the associated single nucleotide polymorphisms as candidate loci before making any inference about γ. If such association is not statistically significant, our proposed methods may not be reliable. In this situation, according to Fieller’s theorem, the Fieller’s CI and the LRtype CI can be discontinuous as shown in Tables 1, 2, 3, and 4, which is difficult to interpret.
Generally, the LR method and the Fieller’s method have accurate CP and control the ML and MR well, and hence are recommended in practical application. In future work, we will incorporate the information on the interval constraint into analysis so as to further improve the efficiency of the proposed methods. Moreover, we will generalize our methods to quantitative traits.
Conclusions
When the sample size is greater than 2000, the Fieller’s method has similar performance to the LR method and thus is preferential due to the noniterative computation procedure. However, the LR method is recommended otherwise because it has better statistical properties, especially in small samples.
Methods
Point estimation for γ
Consider an Xlinked diallelic locus with normal allele a and mutant allele A. We only select the females because XCI is unrelated to males. For females, suppose that aa, Aa and AA are three genotypes and let X={0, γ, 2} be the corresponding genotypic value, respectively, with γ∈ [ 0, 2]. For a casecontrol design, let Y=1 (0) denote that the female is affected (unaffected). Then, the association between Y and X can be expressed using a logistic regression model
where β_{0} is the intercept, β is the regression coefficient for X, z is a vector of covariates that need to be adjusted (e.g. age), and b^{T} is a vector of regression coefficients for z.
To estimate γ, we decompose the genotypic value X as X=γX_{1}+(2−γ)X_{2}, where X_{1}=I_{{G=Aa or AA}}, X_{2}=I_{{G=AA}}, G denotes the genotype of the female and I_{{.}} is the indicator function. It can be seen that X_{1} indicates if the genotype contains the mutant allele A and X_{2} represents if the genotype is the homozygote AA. As such, Model (1) becomes
Let β_{1}=βγ and β_{2}=β(2−γ). Then, the above model can be rewritten as
Further, due to this reparameterization, γ can be represented as
when β=(β_{1}+β_{2})/2≠0. γ can only be well defined in presence of the association between the disease and the allele A. Note that γ∈ [ 0, 2] means that β_{1} and β_{2} have the same sign. That is, the genetic effect of heterozygous genotype in females lies between those of two homozygotes, which is generally satisfied in real applications. From Eq. (3), we have γ=0 (2) if and only if β_{1}=0 and β_{2}≠0 (β_{1}≠0 and β_{2}=0), representing XCIS fully towards the normal (mutant) allele a (A), while γ=1 if and only if β_{1}=β_{2}≠0, which means XCIR. So, if we get the MLEs \(\hat {\beta }_{1}\) and \(\hat {\beta }_{2}\) of β_{1} and β_{2}, then \(\hat {\gamma }=2\hat {\beta }_{1}/(\hat {\beta }_{1}+\hat {\beta }_{2})\) is the MLE of γ by the invariance property of MLE.
Note that \(\hat {\beta }_{1}\) and \(\hat {\beta }_{2}\) can be easily calculated through the standard logistic regression procedure. Specifically, suppose that we collect N unrelated females from a homogeneous population. Then, the loglikelihood function of the sample can be written as
where y_{i}, x_{i1}, x_{i2} and z_{i} respectively are the values of Y, X_{1}, X_{2} and z of female i. Then, \(\hat {\beta }_{1}\) and \(\hat {\beta }_{2}\) are obtained by maximizing the above loglikelihood function, i.e.
where \(\hat {\beta _{0}}\) and \(\hat {\boldsymbol {b}}^{T}\) are the MLEs of β_{0} and b^{T}, respectively.
Confidence interval of γ based on delta method
Once the point estimate of γ is derived, we need to calculate the standard error or CI to evaluate its precision. Since \(\hat {\gamma }\) is also a ratio estimate, a natural idea is to use the first order Taylor series expansion of \(\hat {\gamma }\) and then obtain its asymptotic variance. Specifically, by the consistency of MLE, \(\hat {\beta }_{1}\) and \(\hat {\beta }_{2}\) are close to β_{1} and β_{2}, respectively, when N is large. Note that β=(β_{1}+β_{2})/2 and thus γ can be rewritten as β_{1}/β. Making a first order Taylor expansion of \(\hat {\gamma }\) around the point (β_{1}, β) and evaluating this at \(\left (\hat {\beta }_{1},\ \hat {\beta }\right)\), we have
where \(\hat {\beta }=\left (\hat {\beta _{1}} +\hat {\beta _{2}}\right)/2\). Taking variance from both sides, the above equation becomes
Notice that
and
where \(\text {Var}\left (\hat {\beta _{1}}\right)\), \(\text {Var}\left (\hat {\beta _{2}}\right)\) and \({\text {Cov}}\left (\hat {\beta _{1}},\hat {\beta _{2}}\right)\) are the elements of the variancecovariance matrix V of \(\hat {\beta _{1}}\) and \(\hat {\beta _{2}}\). Generally, V has no simple form when covariates are included in the model, but can be derived from the empirical Fisher’s information matrix \(\boldsymbol {\hat {I}}\) for (β_{0}, β_{1}, β_{2}, b^{T})^{T} [33]. For Model (2),
where U=(1, X_{1}, X_{2}, z) is the design matrix, X_{1}=(x_{11},x_{21},...,x_{N1})^{T}, X_{2}=(x_{12},x_{22},...,x_{N2})^{T}, z=(z_{1},z_{2},...,z_{N})^{T}, \(\boldsymbol {\hat {W}}={\text {diag}}\left (\hat {w}_{1},\hat {w}_{2},...,\hat {w}_{N}\right)\) is a diagonal matrix with diagonal elements
and
represents the estimated penetrances for female i. Once \(\boldsymbol {\hat {I}}\) is estimated, the partial information matrix \(\boldsymbol {\hat {I}}_{1}\) for β_{1} and β_{2} given β_{0} and b^{T} can be computed and thus \(\boldsymbol {V}=\boldsymbol {\hat {I}}^{1}_{1}\). If there is no covariate in the model, then V has the following form [see Appendix A of Additional file 1]
where n_{aa}, n_{Aa} and n_{AA} are the numbers of the females with aa, Aa and AA, respectively, and N=n_{aa}+n_{Aa}+n_{AA}; \(\hat {w}_{aa}\), \(\hat {w}_{aa}\) and \(\hat {w}_{aa}\) are the weighted elements for aa, Aa and AA, respectively, with
and
and
representing the estimated penetrances for aa, Aa and AA, respectively.
Replacing β and β_{1} by \(\hat {\beta }\) and \(\hat {\beta _{1}}\) in Eq. (4), we estimate the deltatype standard error [34]
As such, the deltatype CI \(\left (\gamma _{L}^{d},\gamma _{U}^{d}\right)\) at level (1−α) can be expressed as
where Z_{1−α/2} denotes the (1−α/2)quantile of the standard normal distribution. Note that the estimated CI may be out of the range of [0, 2] when the variation is large, which should be cut off. To test the null hypothesis H_{0}:γ=γ_{0} against the alternative hypothesis H_{1}:γ≠γ_{0}, we have
under H_{0}, where γ_{0} is an arbitrary constant between [ 0, 2], such as 1 (XCIR).
The delta method is a noniterative procedure and thus is easy to be implemented. However, the CI of a ratio estimate is generally skewed, while the deltatype CI is symmetrical [35, 36]. Therefore, it is necessary to propose the Fieller’s and likelihood ratio methods to overcome this shortcoming in the following sections.
Confidence interval of γ based on Fieller’s method
The Fieller’s method is another widely used noniterative approach for constructing CI for ratio estimate [37]. This type of CI can be asymmetrical around \(\hat {\gamma }\). To propose the Fieller’s CI, we first need to build a Wald test for testing γ=γ_{0}. Specifically, under γ=γ_{0}, we have β_{1}−γ_{0}β=0. Therefore, the Wald test for testing γ=γ_{0} can be written as
which follows a standard normal distribution. Then, the confidence limits \(\gamma _{L}^{f}\) and \(\gamma _{U}^{f} \left (\gamma _{L}^{f}<\gamma _{U}^{f}\right)\) for Fieller’s CI at level (1−α) can be found by solving the following equation
Rearranging the above equation yields a quadratic equation with respect to γ_{0}
where
and
Suppose Δ=E^{2}−4DF>0, then this equation must have two unequal roots with \(\gamma _{L}^{f}\) and \(\gamma _{U}^{f}\) being \(\left (E\pm \sqrt {\Delta }\right)/2D\). According to Fieller’s theorem, we know that D>0 implies Δ>0. In this situation, the Fieller’s CI is continuous and can be denoted by \(\left (\gamma _{L}^{f}, \gamma _{U}^{f}\right)\). Note that D>0 is equivalent to \(\left \hat {\beta }/\sqrt {\text {Var}\left (\hat {\beta }\right)}\right  >Z_{1\frac {\alpha }{2}}\). That is, there exists statistically significant association between the disease and the allele A at the significance level α. However, if there is no such association (i.e. D<0), the Fieller’s CI will be unbounded. For instance, if Δ<0, the Fieller’s CI will be (−∞,∞). If Δ>0, the Fieller’s CI will be \(\left (\infty, \gamma _{L}^{f}\right)\bigcup \left (\gamma _{U}^{f},\infty \right)\), which is the discontinuous CI. In real applications, it generally makes little sense to infer about γ if there is no association between the disease and the allele A according to its definition. In addition, the Fieller’s CI should also be restricted to the interval [0,2] when needed.
The Fieller’s method usually demonstrates better coverage probability than the delta method. Notice that the Fieller’s CI is based on the inversion of the Wald test. Since the LR test is expected to have more robust properties in small samples, so it is desirable to propose the LR method in the next section.
Confidence interval of γ based on likelihood ratio method
To obtain the LRbased CI, we first construct a likelihood ratio test for testing γ=γ_{0}. As mentioned above, we have derived the MLEs \(\hat {\beta }_{0}\), \(\hat {\beta }_{1}\), \(\hat {\beta }_{2}\) and \(\hat {\boldsymbol {b}}^{T}\) of β_{0}, β_{1}, β_{2} and b^{T} under H_{1}. To calculate the likelihood ratio test statistic λ, we further evaluate the likelihood function under H_{0}:γ=γ_{0}. If H_{0} holds, the genotypic value X equals 0, γ_{0} and 2 for aa, Aa and AA, respectively. In this regard, Model (1) is reduced to be a standard logistic model and the loglikelihood function under H_{0} can be written as
where x_{i} is the genotypic value of X of female i. Let \(\tilde {\beta }_{0}\), \(\tilde {\beta }\) and \(\tilde {\boldsymbol {b}}^{T}\) be the MLEs of β_{0}, β and b^{T} under H_{0}, respectively. Then,
and λ can be computed as
λ asymptotically follows a chisquare distribution with the degree of freedom being one \(\left (\text {i.e.}\, \chi _{1}^{2}\right)\).
Now, we introduce how to obtain the LRtype CI. For each given γ_{0}, we can calculate the corresponding value of the loglikelihood function l_{0} and \(\tilde {\boldsymbol {\theta }}=\left (\tilde {\beta }_{0}, \tilde {\beta }, \tilde {\boldsymbol {b}}^{T}\right)^{T}\) under H_{0}. So, l_{0} and \(\tilde {\boldsymbol {\theta }}\) are single variable functions with respect to γ_{0} and can be denoted by l_{0}=l_{0}(γ_{0}) and \(\tilde {\boldsymbol {\theta }}=\tilde {\boldsymbol {\theta }}(\gamma _{0})\), respectively. At the significance level α, the confidence limits \(\gamma _{L}^{l}\) and \(\gamma _{U}^{l} \left (\gamma _{L}^{l}<\gamma _{U}^{l}\right)\) of the LRtype CI is determined by the following equation with respect to γ_{0}
where q_{1−α} denotes the (1−α)quantile of the \(\chi _{1}^{2}\) distribution. Obviously, Eq. (5) has no closed form solutions and numerical method can be adopted. Note that θ=(β_{0},β,b^{T})^{T} are nuisance parameters in Eq. (5) which depend on γ_{0}. To solve this equation, it generally requires several iterations with different values of γ_{0}, and for each γ_{0}, the iterative maximization over the remaining parameters is also needed to determine \(\tilde {\boldsymbol {\theta }}\). This procedure is relatively timeconsuming. Therefore, to reduce the computational burden, borrowing the idea of Venzon and Moolgavkar [38], we can find the roots of Eq. (5) by solving the following system of nonlinear equations
which is easily implemented in the commonly used software (e.g. nleqslv package in R). Note that the above system differs only in the first equation from the system (with the first equation being replaced by \(\frac {\partial l_{0}}{\partial \gamma _{0}}(\gamma _{0}, \boldsymbol {\theta })=0\)) that defines the MLEs \(\hat {\gamma }\) and \(\hat {\boldsymbol {\theta }}=\left (\hat {\beta }_{0}, \hat {\beta }, \hat {\boldsymbol {b}}^{T}\right)^{T}\). Therefore, finding a root of such system almost has the same difficulty as that of finding the MLEs of Model (2) [38]. As such, this algorithm is generally more efficient.
On the other hand, based on the fact that the Wald test and the LR test are asymptotically equivalent in large samples, we know that the confidence limits of the Fieller’s CI and the LRtype CI should be close to each other. Therefore, we used the confidence limits of the Fieller’s CI as the initial values for γ_{0}. For example, when searching the lower limit, we chose the initial values for γ_{0} and θ as \(\gamma _{L}^{f}\) and \(\tilde {\boldsymbol {\theta }}\!\left (\gamma _{L}^{f}\right)\), respectively, where \(\tilde {\boldsymbol {\theta }}\!\left (\gamma _{L}^{f}\right)\) can be computed from the standard logistic regression procedure. Similarly, we used the same strategy to search the upper limit. The algorithm based on this choice of the initial values works well in most situations. However, in some scenarios, the Fieller’s CI and the LRtype CI may be very different. Thus, using the confidence limits of the Fieller’s CI as the initial values may cause that the algorithm does not converge. In this regard, we should directly solve the single variable function of Eq. (5). For example, we can use the bisection method to find the roots of Eq. (5) within the interval [ 0,2] (e.g. rootSolve package in R).
Like the Fieller’s CI, the LRtype CI can be unbounded when there is no association between the disease and the allele A. Specifically, when Equation (5) has no root, then the LRtype CI will be (−∞,∞). Otherwise, there will be two roots \(\gamma _{L}^{l}\) and \(\gamma _{U}^{l}\). If \(\gamma _{L}^{l}<\hat {\gamma }<\gamma _{U}^{l}\), then the LRtype CI is continuous and can be represented as \(\left (\gamma _{L}^{l},\gamma _{U}^{l}\right)\). If \(\hat {\gamma }\not \in \left (\gamma _{L}^{l},\gamma _{U}^{l}\right)\), then the LRtype CI will be \(\left (\infty, \gamma _{L}^{l}\right)\bigcup \left (\gamma _{U}^{l},\infty \right)\), which is the discontinuous CI. Similar to the delta and Fieller’s methods, the LRtype CI is also truncated by [ 0, 2] if necessary.
The LRbased CI and the Fieller’s CI can be asymmetrical which is an appealing choice, compared to the delta method. This is because the distribution of a ratio estimate is generally nonnormal with a heavy tail, especially when N is small. Additionally, it will be quite straightforward to incorporate covariates using the LR method.
Simulation settings
For simplicity, we assumed that there is no covariate included in the model in our simulation study. We incorporated the covariate into the real data analysis later. For a casecontrol design, we presumed that the genotype distribution in the case group and that in the control group of females follow trinomial distributions with probabilities (h_{0}, h_{1}, h_{2}) and (g_{0}, g_{1}, g_{2}), respectively, where h_{0} (g_{0}), h_{1} (g_{1}) and h_{2} (g_{2}) are the frequencies of aa, Aa and AA in the case (control) group, respectively. Let h_{0}/g_{0}=m, h_{1}/g_{1}=λ_{1}(h_{0}/g_{0})=λ_{1}m and h_{2}/g_{2}=λ_{2}(h_{0}/g_{0})=λ_{2}m, where λ_{1} and λ_{2} are the odds ratios for Aa and AA compared to aa in females, respectively. By h_{0}+h_{1}+h_{2}=1, we have m=1/(g_{0}+λ_{1}g_{1}+λ_{2}g_{2}), h_{0}=g_{0}/(g_{0}+λ_{1}g_{1}+λ_{2}g_{2}), h_{1}=λ_{1}g_{1}/(g_{0}+λ_{1}g_{1}+λ_{2}g_{2}) and h_{2}=λ_{2}g_{2}/(g_{0}+λ_{1}g_{1}+λ_{2}g_{2}). Thus, (h_{0}, h_{1}, h_{2}) is calculated by (g_{0}, g_{1}, g_{2}), λ_{1} and λ_{2}. The value of (g_{0}, g_{1}, g_{2}) is determined by the allele frequency of A (denoted by p=1−q) and the inbreeding coefficient ρ, i.e. g_{0}=q^{2}+ρpq, g_{1}=2(1−ρ)pq and g_{2}=p^{2}+ρpq. Specifically, ρ=0 implies that HardyWeinberg equilibrium holds in the control group and ρ≠0 indicates HardyWeinberg disequilibrium. Furthermore, from Model (2), β_{0}=log(m), β_{1}=log(λ_{1}) and β_{1}+β_{2}=log(λ_{2}). Then, γ=2log(λ_{1})/log(λ_{2}), i.e. \(\lambda _{1}=\lambda _{2}^{\gamma /2}\). As such, we defined the simulation settings as follows: p is fixed at 0.1 and 0.3, and ρ is set to be 0 and 0.05. The true value of γ is fixed at 0, 0.5, 1, 1.5 and 2. λ_{2} is assigned to be 1.5 and 2. We selected N as 500 (2000) where both the case and control groups have 250 (1000) females. The confidence level is fixed at (1−α)=95% and the number of replications k is 10,000.
We compared the performance of three types of CIs based on CP, ML, MR, ML/(ML+MR) and DP. CP is defined to be the proportion that the CI contains the true value of γ among k replications, regardless of whether the CI is continuous or not. ML and MR are calculated by \(\text {ML}=\#\left [(\gamma _{0}<\gamma _{L})\cap \left (\gamma _{L}\le \hat {\gamma }\le \gamma _{U}\right)\right ]/k\), and \(\text {MR}=\#\left [(\gamma _{0}>\gamma _{U})\cap \left (\gamma _{L}\le \hat {\gamma }\le \gamma _{U}\right)\right ]/k\), respectively, where # denotes the counting measure, and γ_{L} and γ_{U} are the confidence limits of the estimated CI. Note that \(\gamma _{L}\le \hat {\gamma }\le \gamma _{U}\) means that the CI is continuous. As such, we only consider the continuous CIs when estimating ML and MR, since it is impossible to distinguish between the left side and the right side if the CI is discontinuous. Further, DP is computed as \(1\#(\gamma _{L}\le \hat {\gamma }\le \gamma _{U})/k\). We believed that a good CI should control the CP well as well as have the balanced ML and MR. ML and MR are used together to measure the location of CI. If a balance between ML and MR is achieved, then ML/(ML+MR) is close to 0.5. On the other hand, note that the deltatype CI is always bounded. Therefore, the DP for the delta method will stay at 0. However, for the Fieller’s CI and the LRtype CI, small DP is desirable.
Abbreviations
 CI:

Confidence interval
 CP:

Coverage probability
 DP:

Proportion of discontinuous CIs
 GWAS:

Genomewide association study
 LR:

Likelihood ratio
 ML:

Left tail error
 MLE:

Maximum likelihood estimate
 MR:

Right tail error
 XCI; X chromosome inactivation; XCIE:

Escape from X chromosome inactivation
 XCIR; Random X chromosome inactivation; XCIS:

Skewed X chromosome inactivation
References
 1
Lyon MF. Gene action in the Xchromosome of the mouse (Mus musculus L.)Nature. 1961; 190:372–3.
 2
Berletch JB, Ma W, Yang F, Shendure J, Noble WS, Disteche CM, et al.Escape from X inactivation varies in mouse tissues. PloS Genet. 2015; 11:e1005079.
 3
Plenge RM, Stevenson RA, Lubs HA, Schwartz CE, Willard HF. Skewed Xchromosome inactivation is a common feature of Xlinked mental retardation disorders. Am J Hum Genet. 2002; 71:168–73.
 4
AmosLandgraf JM, Cottle A, Plenge RM, Friez M, Schwartz CE, Longshore J, et al.X chromosome–inactivation patterns of 1005 phenotypically unaffected females. Am J Hum Genet. 2006; 79:493–9.
 5
Busque L, Paquette Y, Provost S, Roy DC, Levine RL, Mollica L, et al.Skewing of Xinactivation ratios in blood cells of aging women is confirmed by independent methodologies. Blood. 2009; 113:3472–4.
 6
Minks J, Robinson WP, Brown CJ. A skewed view of X chromosome inactivation. J Clin Invest. 2008; 118:20–3.
 7
Chabchoub G, Uz E, Maalej A, Mustafa CA, Rebai A, Mnif M, et al.Analysis of skewed Xchromosome inactivation in females with rheumatoid arthritis and autoimmune thyroid diseases. Arthritis Res Ther. 2009; 11:R106.
 8
Renault NKE, Pritchett SM, Howell RE, Greer WL, Sapienza C, Ørstavik KH, et al.Human Xchromosome inactivation pattern distributions fit a model of genetically influenced choice better than models of completely random choice. Eur J Hum Genet. 2013; 21:1396–402.
 9
Brown CJ. Skewed Xchromosome inactivation: cause or consequence?J Natl Cancer Inst. 2010; 91:303–4.
 10
Belmont JW. Genetic control of X inactivation and processes leading to Xinactivation skewing. Am J Hum Genet. 1996; 58:1101–8.
 11
Medema RH, Boudewijn MT. The X factor: skewing X inactivation towards cancer. Cell. 2007; 129:1275–86.
 12
Deng X, Berletch JB, Nguyen DK, Disteche CM. X chromosome regulation: diverse patterns in development, tissues and disease. Nat Rev Genet. 2014; 15:367–78.
 13
Zuo T, Wang L, Morrison C, Chang X, Zhang H, Li W, et al.FOXP3 is an Xlinked breast cancer suppressor gene and an important repressor of the HER2/ErbB2 oncogene. Cell. 2007; 129:1275–86.
 14
Li G, Jin T, Liang H, Tu Y, Zhang W, Gong L, et al.Skewed Xchromosome inactivation in patients with esophageal carcinoma. Diagn Pathol. 2013; 8:56–62.
 15
Simmonds MJ, Kavvoura PK, Brand OJ, Newby PR, Jackson LE, Hargreaves CE, et al.Skewed X chromosome inactivation and female preponderance in autoimmune thyroid disease: an association study and metaanalysis. J Clin Endocrinol Metab. 2013; 99:127–31.
 16
Iitsuka Y, Bock A, Nguyen D, SamangoSprouse CA, Simpson JL, Bischoff FZ. Evidence of skewed Xchromosome inactivation in 47, XXY and 48, XXYY Klinefelter patients. Am J Med Genet Part A. 2001; 98:25–31.
 17
Sangha KK, Stephenson MD, Brown CJ, Robinson WP. Extremely skewed Xchromosome inactivation is increased in women with recurrent spontaneous abortion. Am J Hum Genet. 1999; 65:913–17.
 18
Clayton D. Testing for association on the X chromosome. Biostatistics. 2008; 9:593–600.
 19
Clayton D. Sex chromosomes and genetic association studies. Genome Med. 2009; 1:110–6.
 20
Hickey PF, Bahlo M. X chromosome association testing in genome wide association studies. Genet Epidemiol. 2011; 35:664–70.
 21
Wang J, Yu R, Shete S. Xchromosome genetic association test accounting for Xinactivation, skewed Xinactivation, and escape from Xinactivation. Genet Epidemiol. 2014; 38:483–93.
 22
Chen Z, Ng HKT, Li J, Liu Q, Huang H. Detecting associated singlenucleotide polymorphisms on the X chromosome in case control genomewide association studies. Stat Methods Med Res. 2017; 26:567–82.
 23
Ma L, Hoffman G, Keinan A. Xinactivation informs variancebased testing for Xlinked association of a quantitative trait. BMC Genomics. 2015; 16:241–9.
 24
Busque L, Mio R, Mattioli J, Brais E, Blais N, Lalonde Y, et al.Nonrandom Xinactivation patterns in normal females: lyonization ratios vary with age. Blood. 1996; 88:59–65.
 25
Szelinger S, Malenica I, Corneveaux JJ, Siniard AL, Kurdoglu AA, Ramsey KM, et al.Characterization of X chromosome inactivation using integrated analysis of wholeexome and mRNA sequencing. Plos One. 2014; 9:e113036.
 26
Xu SQ, Zhang Y, Wang P, Liu W, Wu XB, Zhou JY. A statistical measure for the skewness of X chromosome inactivation based on family trios. BMC Genet. 2018; 19:109.
 27
Chu X, Shen M, Xie F, Miao XJ, Shou WH, Liu L, et al.An X chromosomewide association analysis identifies variants in GPR174 as a risk factor for Graves’ disease. J Med Genet. 2013; 50:479–85.
 28
Dobyns WB, Filauro A, Tomson BN, Chan AS, Ho AW, Ting NT, et al.Inheritance of most Xlinked traits is not dominant or recessive, just Xlinked. Am J Med Genet Part A. 2004; 129:136–43.
 29
Biggar RJ, Bergen AW, Poulsen GN. Impact of X chromosome genes in explaining the excess risk of cancer in males. Am J Epidemiol. 2009; 170:65–71.
 30
Ober C, Loisel DA, Gilad Y. Sexspecific genetic architecture of human disease. Nat Rev Genet. 2008; 9:911–22.
 31
Lourenço CM, Simão GN, Santos AC, Marques JW. Sexspecific genetic architecture of human disease. Arq Neuropsiquiatr. 2012; 70:487–91.
 32
Molenberghs G, Verbeke G. Likelihood ratio, score, and Wald tests in a constrained parameter space. Am Stat. 2007; 61:22–7.
 33
Wood SN. Generalized Additive Models: An Introduction With R. 1st ed. London: Chapman & Hall Ltd; 2006.
 34
Oehlert GW. A note on the delta method. Am Stat. 1992; 46:27–9.
 35
Tin M. Comparison of some ratio estimators. J Am Stat Assoc. 1965; 60:294–307.
 36
Choquet D, L’ecuyer P, Léger C. Bootstrap confidence intervals for ratios of expectations. Acm T Model Comput S. 1999; 9:326–48.
 37
Fieller EC. Some problems in interval estimation. J R Stat Soc B. 1954; 16:175–85.
 38
Venzon DJ, Moolgavkar SH. A method for computing profilelikelihoodbased confidence intervals. J R Stat Soc CAppl. 1988; 37:87–94.
Acknowledgements
The authors would like to thank the editor and three anonymous reviewers for their valuable comments which highly improved the presentation of the article and the authors also appreciate Dr. Yuqiang Xue’s kind assistance.
Funding
This work was supported by the National Natural Science Foundation of China [81773544, 81373098, 81573207], the National and Guangdong University Students’ Innovation and Enterprise Training Project of China [201612121017], and the General Program of Applied Basic Research Programs of Yunnan Province [2017FB002]. All the funding supporters had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Availability of data and materials
The dataset supporting the conclusions of this article can be found in the article of Chu et al. [27].
Author information
Affiliations
Contributions
PW helped design the study, drafted the article, and conducted the simulation study. YZ helped conduct the simulation study and the data analysis, and design the study’s analytic strategy. BQW revised the article critically and helped interpret the results of the data analysis. JLL helped design the study, plot the figures and conduct the literature review. YXW helped interpret the results of the data analysis and prepare the introduction section and the discussion section. DP helped draft the article, helped the analysis and the interpretation of the data, and revised the article. XBW and WKF helped design the study, reviewed the whole paper and critically revised the article. JYZ helped design the study, supervised the field activities, and directed its implementation, including quality assurance and control. All authors read and approved this version of the manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file
Additional file 1
Appendices and Supplementary tables and figures. Appendix A Derivation of the closed form of V without covariates; Appendix B Size and power results for testing γ=γ_{0}; Tables S1S2 Estimated CP (%), ML (%), MR (%), ML/(ML+MR) and DP (%) of the twosided 95% CI when N=500, ρ=0.05, and λ_{2}=1.5 or 2 for the LR, Fieller’s and delta methods, respectively; Tables S3S4 Estimated CP (%), ML (%), MR (%), ML/(ML+MR) and DP (%) of the twosided 95% CI when N=2000, ρ=0.05, and λ_{2}=1.5 or 2 for the LR, Fieller’s and delta methods, respectively; Tables S5S8 Estimated size (%) for testing H_{0}: γ=γ_{0} with N=500 or 2000, and ρ=0 or 0.05 for the LR, Fieller’s and delta methods, respectively; Figures S1S6 Power comparison of the LR, Fieller’s and delta methods against γ. The simulation is based on 10,000 replicates with N=500, ρ=0 or 0.05, and γ_{0}=0, 1 or 2, respectively; Figures S7S12 Power comparison of the LR, Fieller’s and delta methods against γ. The simulation is based on 10,000 replicates with N=2000, ρ=0 or 0.05, and γ_{0}=0, 1 or 2, respectively. (PDF 117 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Wang, P., Zhang, Y., Wang, B. et al. A statistical measure for the skewness of X chromosome inactivation based on casecontrol design. BMC Bioinformatics 20, 11 (2019). https://doi.org/10.1186/s1285901825872
Received:
Accepted:
Published:
Keywords
 X chromosome inactivation
 Skewness
 Casecontrol design
 Confidence interval
 Graves’ disease