Gene-Based Methods for Estimating the Degree of the Skewness of X Chromosome Inactivation

Skewed X chromosome inactivation (XCI-S) has been reported to be associated with some X-linked diseases, and currently several methods have been proposed to estimate the degree of the XCI-S (denoted as γ) for a single locus. However, no method has been available to estimate γ for genes. Therefore, in this paper, we first propose the point estimate and the penalized point estimate of γ for genes, and then derive its confidence intervals based on the Fieller’s and penalized Fieller’s methods, respectively. Further, we consider the constraint condition of γ∈[0, 2] and propose the Bayesian methods to obtain the point estimates and the credible intervals of γ, where a truncated normal prior and a uniform prior are respectively used (denoted as GBN and GBU). The simulation results show that the Bayesian methods can avoid the extreme point estimates (0 or 2), the empty sets, the noninformative intervals ([0, 2]) and the discontinuous intervals to occur. GBN performs best in both the point estimation and the interval estimation. Finally, we apply the proposed methods to the Minnesota Center for Twin and Family Research data for their practical use. In summary, in practical applications, we recommend using GBN to estimate γ of genes.


Introduction
X chromosome inactivation (XCI) is an important epigenetic phenomenon. Under the XCI, one of two X chromosomes in females is silenced in the early stage of embryonic development to ensure that the transcriptional dosage of X chromosomes in females and that in males are balanced [1]. Generally, there are three patterns of the XCI [2], random X chromosome inactivation (XCI-R), skewed X chromosome inactivation (XCI-S) [3][4][5][6], and escape from X chromosome inactivation (XCI-E) [7,8]. The XCI-R means that the paternal and maternal X chromosomes in females have the same probabilities to be inactive, i.e., for a locus on the X chromosome, approximately 50% of the cells inactivate one of the alleles, while the remaining 50% of the cells keep the other allele inactive. Under the XCI-E, the alleles on both the X chromosomes in females are expressed, which are similar to those at an autosomal locus. For humans, about 15-30% of the X-linked genes have been reported to undergo the XCI-E [7]. Finally, the XCI-S is defined as more than 75% of the cells in females inactivating the same allele [9]. For some extreme skewed cases, it is possible that more than 90% of the cells keep the same allele silenced [9,10]. As such, the difference in the number of the X chromosomes in females and males and the complexity of the XCI make the association tests for the X chromosomes more complicated than those for the autosomes.
The skewness of the XCI can reflect, or cause, biological consequences for females [9]. The clonal expansion of a somatic cell in females may lead to a cell population with extremely skewed XCI [9]. For some X-linked disorders, there is strong selection of the cells which keep the mutant allele inactive in the heterozygous carriers and, hence, assessing the degree of the skewness of the XCI is helpful in terms of being indicative of the carrier's disease status [11]. Further, the degree of the skewness of the XCI can determine the severity of certain X-linked diseases, such as haemophilia B [12,13]. On the other hand, even for the same mutant allele, the XCI-S in different tissues or cells may result in different clinical consequences. For example, in heterozygous females with a mutant FoxP3 allele, the XCI-S against the mutant allele in specific tissues can prevent autoimmune disease, while the XCI-S skewed towards the mutant allele in breast epithelial cells can cause breast cancer [14]. Besides this, studies have shown that some diseases, such as ovarian cancer, Rett syndrome, Duchenne muscular dystrophy and recurrent miscarriage, are also related to the XCI-S [15][16][17][18]. Therefore, in recent years, researchers have proposed some methods to test the association between the alleles at an X-chromosomal single nucleotide polymorphism (SNP) locus and traits [19][20][21][22][23][24][25][26]. For example, Wang et al. [23] developed a permutation-based test statistic which considers all the XCI patterns. For the XCI-R and the XCI-S, this method respectively codes three female genotypes (dd, Dd and DD) as 0, γ and 2 at an X-chromosomal SNP, with the major allele d and the minor allele D, where γ ∈ [0, 2] is an unknown genotypic value for heterozygous females, and respectively codes two male genotypes (d and D) as 0 and 2. Here, γ can be used to measure the degree of the XCI skewing. For instance, γ ∈ [0, 1) is indicative of the XCI-S skewed towards the minor allele D, γ = 1 means that the XCI pattern is the XCI-R, and γ ∈ (1, 2] indicates the XCI-S skewed towards the major allele d. For the XCI-E, three female genotypes are coded as 0, 1 and 2, and two male genotypes are coded as 0 and 1, respectively. However, the X-chromosomal association tests mentioned above are only applicable to a single SNP and common variants, and are not suitable for genetic regions or genes with multiple SNPs and rare variants. Rare variants refer to the variants with a minor allele frequency (MAF) less than 1%, and those with MAF ≥ 1% are called common variants [27,28]. Over the past few years, genome-wide association studies have identified many common variants associated with complex traits, but these variants usually explain only a small part of the estimated heritability for a given trait. On the other hand, it has been shown that rare variants play a key role in influencing traits [29]. Single-variant tests often have low test power when applied to the rare variants. Therefore, many statistical methods had been presented, which focus on testing the cumulative effect of rare variants in genetic regions or SNP sets (such as genes), including the burden test and the variance-component tests [27,[30][31][32][33]. The burden test collapses all the rare variants in a genetic region into a single burden variable, and then regresses the trait on the burden variable to test the cumulative effect of the rare variants in that region [27]. The variance-component tests, such as the sequence kernel association test (SKAT), do not directly aggregate the variants in the modeling process, but aggregate the association between the variants and the trait through a kernel matrix [33]. Another method, SKAT-O, proposed by Lee et al. [34], has the advantages of both the burden and SKAT tests, but the time cost is higher than the previous two methods. All these methods have one thing in common, i.e., increasing the weights of rare variants' contributions and decreasing the weights of common variants' contributions. However, for a trait-related gene, the relative influence of rare and common variants is not known [35]. Therefore, Iuliana et al. [35] put forward several multi-locus association tests, such as the adaptive sum test, which consider the effects of both common and rare variants on the trait, and these methods are more powerful when the genes simultaneously contain rare and common variants. Note that these multi-locus association tests are all based on genetic regions or genes on autosomes, and may not be directly applied to the X chromosomes. Therefore, Clement et al. [36] improved the traditional burden test, SKAT and SKAT-O methods and suggested three gene-based X-chromosomal association tests. However, these methods only take account of the XCI-R and XCI-E patterns. What is more, the FxSKAT method, proposed by Asuman et al. [37], is not only applicable to pedigree data, but also takes the XCI-E into account during the analysis process.
Except for testing the association between the genes on the X chromosome and the traits under study, it is also important to develop methods to measure the corresponding degree of the skewness of the XCI (denoted as γ). At present, researchers have put forward several methods to estimate γ for a single SNP, which can simultaneously get the point estimates and the confidence intervals (CIs) of γ. Specifically, Xu et al. [38] proposed a statistical index for estimating γ based on family trios (both parents and their daughter), which can be represented as the ratio of two relative risks in females, and derived the corresponding CI with the likelihood ratio (LR) test. Wang et al. [39] used the ratio of two regression coefficients of a logistic regression to estimate γ, and obtained the CIs with the LR, Fieller's and delta methods, respectively. Li et al. [40] further extended the methods of Wang et al. so that they can accommodate quantitative traits. However, the above-mentioned methods are all constructed for a single SNP, and are not suitable for genetic regions or genes containing multiple SNPs. Furthermore, when applied to rare variants, they perform poorly. In addition, it should be noted that the delta method cannot control the coverage probability (CP) well, and the LR and Fieller's methods have similar performance in the interval estimation, while the Fieller's method is computationally efficient. Thus, the Fieller's method is recommended in practice. However, both the LR and Fieller's methods may yield unbounded CIs when the denominators in the ratios used to estimate γ are close to 0. Fortunately, the penalized Fieller's (PF) method, which was proposed by Wang et al. [41], can be used to conduct the ratio estimation and always get the bounded CIs by choosing an appropriate penalty parameter. However, it has not been applied to the estimation of the degree of the skewness of the XCI yet. On the other hand, the above-mentioned methods do not consider the constraint condition of γ ∈ [0, 2], and simply cut off the point estimates and the CIs within [0, 2], which may result in extreme point estimates (0 or 2) and empty sets or noninformative CIs (i.e., [0,2]). In contrast, the Bayesian methods can effectively utilize the prior information of each unknown parameter in the analysis, and have been widely used in statistical genetics [42]. Therefore, in this paper, we borrow the idea of the burden test, aggregate all the variants in a gene under study into a burden variable by selecting the appropriate weights, and then estimate the mean degree of the skewness of the XCI over all the SNPs in the gene based on the burden variable. We first propose the point estimate and the penalized point estimate of γ for the gene, and then derive its CIs based on the Fieller's and PF methods, respectively. Then, by considering the constraint condition of γ ∈ [0, 2], we propose the Bayesian methods to obtain the point estimates and the credible intervals of γ. Specifically, after getting enough samples drawn from the posterior distribution of γ, we calculate the mode of the samples as the point estimate of γ and the highest posterior density interval (HPDI) as the credible interval of γ [43]. We conduct extensive simulation studies to compare the performances of the proposed point estimation methods and the interval estimation methods for γ. Finally, we demonstrate the practical utility of the proposed methods by applying them to the Minnesota Center for Twin and Family Research (MCTFR) data.

Notations
Suppose that we only collect n female subjects, because male subjects provide no information on the XCI skewing. Consider an X-linked trait (quantitative or qualitative) and let y i represent the trait value of the ith female (i = 1, 2, . . . , n), then Y = (y 1 , y 2 , . . . , y n ) T is the vector of the trait values for all the females. Assume that a gene which contains J SNPs is associated with this trait, and let d j and D j denote the major allele and the minor allele at the jth SNP (j = 1, 2, . . . , J), respectively. Let G ij be the genotype at the jth SNP of the ith female (i.e., G ij = d j d j , D j d j or D j D j ). If we use γ ∈ [0, 2] to measure the mean degree of the skewness of the XCI for all the SNPs in the gene, then g ij = 0, γ and 2 can be used to denote the genotypic values for genotypes d j d j , D j d j and D j D j , respectively. As such, G i = g i1 , g i2 , . . . , g iJ T is the vector of the genotypic values at the J SNPs of the ith female. Therefore, we consider the association between the gene and the trait based on the following generalized linear model where h(·) is a link function; Z i = (Z i1 , Z i2 , . . . , Z im ) T is the vector of m covariates of the ith female, which are needed to be adjusted, and Z = (Z 1 , Z 2 , . . . , Z n ) T is an n × m covariate matrix; µ i = E(y i |G i , Z i ) is the conditional mean of the ith female's trait value given G i and Z i ; β 0 is the intercept, β = β 1 , β 2 , . . . , β J T is the vector of the regression coefficients of G i , Based on the idea of the burden test [27], we aggregate all the SNPs in the gene into a burden variable and let X i = J ∑ j=1 ω j g ij , where ω j is a weight for the jth SNP. Here we assume that ω j is a function with respect to the MAF at the jth SNP (denoted as MAF j ), i.e., ω j = Beta MAF j , 0.5, 0.5 [35]. So, model (1) can be rewritten as where β c is the regression coefficient of X i . Next, we consider two variables g ij = 1 means that the ith female contains at least one minor allele at the jth SNP, and g (2) ij = 1 denotes that the female is a homozygote D j D j at the jth SNP. Through simple transformations, we can get g ij = γg (1) ij + (2 − γ)g (2) ij , and X i can be expressed as ij and X To estimate the mean degree of the XCI skewing for the gene (i.e., γ), we substitute X i = γX (1) i + (2 − γ)X (2) i into model (2) and get For quantitative traits, h(·) is the identity function, and model (3) can be written as where ε i is the random error and follows N 0, σ 2 . In this case, the unknown parameters are , and the corresponding likelihood function of the sample is As for qualitative traits, h(·) is the logit function, and model (3) is written as The unknown parameters are θ 2 = β 0 , β c , γ, b T T and the likelihood function is where y i = 1 and 0 respectively indicate that the ith female is a case and a control, As such, we obtain β c = β (1) c + β (2) c /2 and γ can be expressed as By assuming that the degree of the skewness of the XCI at the jth SNP is γ j , γ satisfies, under a certain condition (the proof is given in Appendix A), .j − g (2) .j γ j ∑ J j=1 ω j g (1) .j − g (2) .j , where g (1) .j = n ∑ i=1 g (1) ij is the number of the females who contain at least one minor allele at the jth SNP, and g (2) .j = n ∑ i=1 g (2) ij is the number of the females whose genotypes at the jth SNP are D j D j . So, γ is the weighted mean of the γ j 's for all the SNPs in the gene with the weights being ω j g (1) .j − g (2) .j / J ∑ j=1 ω j g (1) .j − g (2) .j . When there are rare variants at some SNPs or when the variation of the γ j 's in the gene is large, γ is still well defined for the whole gene. On the other hand, from Equation (5), γ can be well defined if there is an association between the gene and the trait (i.e., β c = β (1) c + β  (2) c = 0, indicating that all the γ j 's are 2 and the XCI-S is completely skewed towards the major allele for each SNP. However, γ = 1 means that on the average, the gene undergoes the XCI-R or the XCI-E. After obtaining the estimates of β (1) c and β (2) c , respectively denoted byβ (1) c andβ (2) c which can be derived by the maximum likelihood method, the point estimate of γ can be expressed asγ = 2β

Point Estimate and CI of γ by Fieller's Method
Note that γ should take the possible values from the interval [0, 2]. So, the original estimateγ = 2β (1) c / β (1) c +β (2) c needs to be cut off in [0, 2] and the resulting estimate is denoted byγ GF . Further, we utilize the Fieller's method to get the CI of γ. Specifically, borrowing the idea of Wang et al. [39], we haveβ c = β (1) c +β (2) c /2, To construct the CI of γ, we first establish a Wald test under the null hypothesis H 0 : γ = γ 0 , where γ 0 is a pre-specified value (e.g., 1, which means that on the average, the gene undergoes the XCI-R or the XCI-E). As such, we have β (1) c − γ 0 β c = 0, and the Wald test statistic isβ (1) Therefore, the 100(1 − α)% CI of γ can be derived by solving the following equation where Z 1−α/2 is the (1 − α/2) upper quantile of the standard normal distribution. Rearrange the above equation with respect to γ 0 into a quadratic equation where , the CI of γ will degenerate to be a point. The CI of γ for other cases is as follows It should be noted that even in the case of ∆ > 0, the CI of γ obtained by the Fieller's method may still be an empty set. And in the case of ∆ > 0 and A < 0, the CI may be composed of two parts, which is the discontinuous interval.

Penalized Point Estimate and CI of γ by PF Method
As mentioned above, we constructγ =β (1) c /β c as the point estimate of γ, wherê β c = β (1) c +β (2) c /2. However, if the denominatorβ c is very close to 0,γ will tend to the infinity. The CI of γ based on the Fieller's method before the truncation is usually unbounded. To deal with this issue in the ratio estimate and borrow the idea of Wang et al. [41], we propose the following PF method to obtain the penalized point estimate of γ and the corresponding CI. Consider the penalized log-likelihood function of β c as follows: where λ > 0 is a penalty parameter and is taken to be Z 2 1−α/2 /4 as suggested by Wang et al. [41] because the CI obtained by the PF method is always bounded with λ = Z 2 1−α/2 /4. By maximizing the function pl, we have the penalized denominator β c =β c /2 + sign β c β2 c /4 + λV ar β c , where sign(·) is the signum function. Further, we can getV ar If we replaceβ c by β c to obtain the point estimate γ =β (1) c / β c , then γ is a biased estimate of γ. To reduce this bias, we need to correct the numeratorβ (1) c by β (1) c =β (1) c + γ β c −β c . Correspondingly, we can getV ar β (1) c ,β c + 4(1 − ξ) 2 γ 2V ar β c andĈ ov β (1) c , β c =Ĉ ov β (1) c ,β c − 2ξ(1 − ξ) γV ar β c . After ob-taining the corrected denominator β c and the corrected numerator β (1) c ,γ * = β (1) c / β c truncated by [0, 2] is the penalized point estimate of γ, which is denoted byγ GPF . The construction process of the corresponding CI ofγ GPF is the same as the Fieller's method, except for respectively replacingβ c ,β (1) c ,V ar β c ,V ar β (1) c andĈ ov β (1) c ,β c by β c , β (1) c , c , β c in Equation (6). However, it should be noted that although the CI of γ based on the PF method is always bounded when λ = Z 2 1−α/2 /4, it may be out of [0, 2] and we need to truncate it by [0, 2].

Point Estimate and Credible Interval of γ by Bayesian Method
Note that the point estimates (γ GF andγ GPF ), and the corresponding CIs mentioned above, are cut off in the interval [0, 2] and cannot directly incorporate the information on γ ∈ [0, 2]. Therefore, in this subsection, we introduce the Bayesian method to give the point estimate and the credible interval of γ by considering the prior information of γ ∈ [0, 2]. Specifically, we have the posterior distribution of the unknown parameter θ . as follows f θ . Y, X (1) , X (2) , where f (θ . ) is the joint prior distribution of θ . ; when the traits are quantitative, θ . = θ 1 and L . (θ . ) = L 1 (θ 1 ); when the traits are qualitative, θ . = θ 2 and L . (θ . ) = L 2 (θ 2 ). However, in general, we cannot get the analytical solutions of f θ . Y, X (1) , X (2) , Z . Therefore, it is not feasible to directly sample from the posterior distribution. Fortunately, there are several algorithms for sampling from an approximate distribution of the posterior distribution, such as the Hamiltonian Monte Carlo (HMC) algorithm which can be implemented by the "rstan" package in R [43]. On the other hand, according to Annis et al. [43], the correlation between the parameters has little influence on the HMC algorithm. To simplify the operations, and improve the sampling efficiency, we assume that the unknown parameters in θ . are independent of each other, and use the HMC algorithm to sample from the approximate posterior distribution of θ . . In other words, we choose the prior distribution for each unknown parameter separately. The prior distributions of the parameters in θ . are selected as follows. To reduce the influence of the selection of the prior distributions on the results, for nuisance parameters β 0 , β c and b (there is an additional nuisance parameter σ when the trait is quantitative), we choose the weak prior distributions [44]. Specifically, we assume that the prior distributions of β 0 and β c are both N 0, 10 2 , and that of b is MV N 0, diag 10 2 , 10 2 , . . . , 10 2 . For quantitative traits, we also specify the prior distribution of σ to be an exponential distribution, i.e., σ ∼ exp(1). As for the parameter γ of interest, which is used to measure the mean degree of the skewness of the XCI over all the SNPs in the gene and satisfies the constraint condition of γ ∈ [0, 2], we consider two possible prior distributions. The first one is the truncated normal distribution, with both parameters being 1 and the values ranging from 0 to 2, and the probability density function of the prior distribution is where φ(·) is the probability density function of the standard normal distribution. In this way, γ not only satisfies the constraint condition of γ ∈ [0, 2], but also the probability of γ being close to 1 is the highest, which is consistent with the literature [2], i.e., most of the SNPs on the X chromosome undergo the XCI-R. Meanwhile, the selected truncated normal distribution of γ also avoids that the probability of γ taking the extreme value (0 or 2) is too low, which may be more suitable for practical applications. The second prior distribution of γ is a uniform distribution, i.e., γ ∼ U(0, 2).
After specifying the prior distributions of all the unknown parameters, we can get enough samples of γ through the HMC algorithm, and then calculate the mode of the samples as the point estimate of γ, and the highest posterior density interval (HPDI) as the credible interval of γ. Here, we denote the Bayesian methods with the truncated normal prior and the uniform prior as GBN and GBU, and the point estimates obtained by these two methods are denoted asγ GBN andγ GBU , respectively.

Simulation Settings
We conducted extensive simulation studies to evaluate the performances of the proposed point estimation and interval estimation methods. The number of female subjects (i.e., the sample size n) is set to be 500 and 2000. Consider a gene associated with the trait under study and the number of the SNPs in the gene (i.e., J) is fixed at 100, i.e., we assume that all the 100 SNPs are associated with the trait. Meanwhile, we define η as the proportion of rare variants among the 100 SNPs. To explore the effect of η on the proposed methods, we set η = 0, 0.4 and 1, which correspond to the cases of all the 100 SNPs only including common variants, the 100 SNPs simultaneously containing common and rare variants, and all the 100 SNPs only consisting of rare variants, respectively. Among them, the MAFs for common variants are sampled from U(0.01, 0.5), while the MAFs for rare variants are randomly simulated from U(0.005, 0.01) [45][46][47]. We generate the genotypes of n female subjects by referring to the ideas of Wang et al. [45], Basu et al. [46], and Turkmen et al. [47]. We first generate a latent vector V = (V 1 , V 2 , . . . , V 100 ) T from the multivariate normal distribution with the mean vector being 0 and the elements of the variance-covariance matrix satisfying Var(V j ) = 1 and Corr V j , V k = ρ |j−k| (j, k = 1, 2, . . . , 100) [45,47], where the linkage disequilibrium among the SNPs is taken into consideration. For simplicity, we set ρ = 0.5 in our simulation studies. Once V is generated, it is then transformed to 0 (major allele) or 1 (minor allele) determined by the corresponding MAFs. This process is repeated twice, and two simulated vectors of length 100 are put together to form the genotypes at the 100 SNPs for a female subject. After simulating the genotypes of n female subjects, we have an n × 100 genotypic value matrix G = (G 1 , G 2 , . . . , G n ) T with the elements being 0, 1 or 2, and then we replace the elements of G equal to 1 with γ to simulate the information on the XCI-S. Note that to simplify the simulation and better evaluate the performances of our proposed methods (e.g., the calculation of the mean squared errors (MSEs) of the point estimates requires a single true value of γ for each replicate; the details are given later), we set the degrees of the XCI skewing γ j 's at all the 100 SNPs to be the same in the simulation study (i.e., γ j = γ, j = 1, 2, . . . , 100).
We only consider a covariate Q, which is generated from the standard normal distribution. For the quantitative trait, we simulate the trait value y i of the ith female according to the following model where ε i is the random error, which is generated from the standard normal distribution; β 0 is the intercept and δ is the regression coefficient of the covariate Q, and both the parameters are set to be 0.5 [36]; β j = e log 10 MAF j /2 is the regression coefficient of the genotypic value g ij at the jth SNP [33,34,36], where e is the tuning parameter and is used to avoid the effect of a SNP being too large or too small [36]. To highlight the effects of rare variants on the trait, we set e = 1.5 when the jth SNP has a rare variant, otherwise e = 1.1. Further, notice that the directions of the effects of different SNPs on the trait may be different. Therefore, we consider the proportion of the SNPs with positive effects among the 100 SNPs (denoted by τ) and set τ to be 0.6 and 1, indicating that the effect directions of some SNPs are positive and some are negative, and all the SNP effects are positive, respectively. We do not simulate the case of τ = 0 (i.e., all the SNP effects are negative) because the results with τ = 0 are similar to those with τ = 1. As for the qualitative trait, the model for generating the affection status y i of the ith female is as follows All of the parameters are the same as when simulating the quantitative trait, except that we need to set the case-control ratio to be 1:1.
After simulating the genotypes and the trait values, we use model (4) to obtain the estimates of β (1) c and β andM AF j is the estimate of the MAF at the jth SNP. Then, we get the point estimateγ GF , the penalized point estimateγ GPF , and the CIs of γ derived by the Fieller's and the PF methods. As for the Bayesian methods, the HMC algorithm is implemented through the "sampling" function in the R package "rstan". We set 8 chains for the parallel sampling in the simulation. For each chain, we extract 10,000 samples, and the first 5000 are used for warm-up. So, we finally get 40,000 samples. To ensure the convergence, the target acceptance rate is set to be 0.99. The above simulation steps are all implemented in the R software (version 4.1.1, http://r-project.org, accessed on 5 January 2022). For each simulation setting, the number of the replicates is fixed to be 500, and for each replicate, the true value of γ is sampled from the uniform distribution U(0, 2). To evaluate the accuracy and the robustness ofγ GBN ,γ GBU ,γ GPF andγ GF , we calculate the MSEs of these point estimates. Here, (γ s − γ s ) 2 /500, where γ s represents the true value of γ andγ s is the point estimate in the sth replicate (s = 1, 2, . . . , 500). Note thatγ GBN andγ GBU are always between 0 and 2, so we only compute the proportions ofγ GPF andγ GF taking the extreme values (0 or 2), respectively. Meanwhile, scatter plots are used to show the relationship between the four point estimates and the true values of γ. To compare the performances of the GBN, GBU, PF and Fieller's methods in the interval estimation, we calculate the CP as well as the mean, the median, the standard deviation and the interquartile range of the widths of the 95% HPDIs or CIs (denoted by W mean , W median , W sd and W iqr ), respectively. For the PF and Fieller's methods, we also compute the proportions of the empty sets (EP), the noninformative intervals (NP), and the discontinuous intervals (DP) to further compare the effectiveness of these two methods, where the noninformative interval means the CI being [0, 2]. However, it should be noted that the GBN and GBU methods avoid the cases of empty sets, noninformative intervals, and discontinuous intervals occurring. In addition, we draw the scatter plots between the interval widths of the four proposed methods and the true values of γ.

Simulation Results
The proportions of the extreme values (0 or 2) forγ GPF andγ GF are shown in Table 1. It can be seen from the table that the proportions of the point estimates equal to 0 are the same for bothγ GPF andγ GF , while the proportion of the point estimates equal to 2 for γ GPF is reduced. This is because before the truncation, bothγ * = β (1) c / β c andγ =β (1) c /β c always have the same sign, andγ * is bounded. Specifically, whenγ * andγ are negative, γ GPF andγ GF are both 0. On the other hand, whenγ * andγ are positive, compared withγ, the proportion ofγ * being greater than 2 decreases. Further, from Table 1, with the increase of the sample size or the trait changing from qualitative to quantitative, the proportions of the extreme values forγ GPF andγ GF both become less. Next, let us take a look at the effects of the proportion of the rare variants (η) and the proportion of the SNPs with the positive effects (τ) among all the SNPs on the proportions of the extreme values forγ GPF andγ GF . Under the situation that the trait is quantitative and τ = 0.6 (i.e., the effect directions of some SNPs are positive and some are negative), the proportions of the extreme values (0 and 2) forγ GPF andγ GF with η = 0 (all the SNPs only include common variants) are less than those with η = 1 (all the SNPs only consist of rare variants), irrespective of the sample size (n). As for the qualitative trait, when n = 2000 and τ = 0.6, the proportion of the extreme values equal to 0 forγ GPF and the proportions of the extreme values (0 and 2) for γ GF with η = 0 are smaller than those with η = 1, while the proportion of the extreme values equal to 2 forγ GPF with η = 0 (12.8%) is larger than that with η = 1 (10.4%). When the trait is qualitative, n = 500 and τ = 0.6, the results are similar to those with n = 2000, except that the proportion of the extreme values equal to 2 forγ GF with η = 0 (20.0%) and that with η = 1 (19.2%) are very close to each other. In addition, the proportions of the extreme values (0 or 2) forγ GPF andγ GF . have no obvious trends for other cases of different values of η and τ. The MSEs of the four point estimates (γ GBN ,γ GBU ,γ GPF andγ GF ) are listed in Table 2. From Table 2, we can see that the MSEs ofγ GBN andγ GBU are smaller than those ofγ GPF andγ GF , and the MSE ofγ GBN is the smallest. When the sample size increases or the trait turns from qualitative to quantitative, the MSEs of these four point estimates decrease significantly. In general, the MSEs of the four point estimates gradually become larger when η changes from 0, 0.4 to 1 (i.e., higher proportion of rare variants) and other parameters are kept unchanged, except for the case when the trait is quantitative, n = 500 and τ = 1. On the other hand, the MSEs of the four point estimates with τ = 0.6 (i.e., the effect directions of some SNPs are positive and some are negative) are smaller than those with τ = 1 (i.e., all the SNP effects are positive), when other parameters are fixed.   Figures 1 and 2 are the scatter plots of the four point estimates against the true values of γ for the quantitative trait with n = 500, and τ = 0.6 and 1, respectively. In each figure, subplots (a)-(d) (four subplots in the first row) are respectively the scatter plots ofγ GBN ,γ GBU ,γ GPF andγ GF with η = 0; subplots (e)-(h) (four subplots in the second row) and subplots (i)-(l) (four subplots in the third row) are the corresponding scatter plots with η = 0.4 and 1, respectively. By comparing the four subplots in the same row of each figure, we find that the two point estimates (γ GBN andγ GBU ) obtained by the Bayesian methods are closer to the true values of γ, and both perform better thanγ GPF andγ GF . On the other hand, note that the distribution of the true value of γ is U(0, 2), and it can be seen from the figures that the distributions ofγ GBN andγ GBU are more uniform, while the distributions ofγ GPF andγ GF are skewed towards the extreme values (0 and 2). Meanwhile, by respectively comparing subplots (a), (e) and (i) forγ GBN with subplots (b), (f) and (j) forγ GBU , there is a little greater dispersion forγ GBU thanγ GBN . In addition, from subplots (c), (g) and (k) forγ GPF and subplots (d), (h) and (l) forγ GF , we observe that there exist many extreme point estimates forγ GPF and γ GF (represented by the blue points). Moreover, the scatter plots forγ GPF andγ GF provide the additional information that most of the extreme point estimates generally occur when the true values of γ are less than 0.5 or greater than 1.5. Further, by comparing the subplots in different rows of each figure when τ = 0.6 (Figures 1, S1, S3 and S5), i.e., η changing from 0, 0.4 to 1, the dispersions of the four point estimates generally increase, indicating that, in general, the MSEs of the four point estimates become larger, which are consistent with the results in Table 2. The numbers of the blue points in subplots (c) and (d) with η = 0 are much less than those in subplots (k) and (l) with η = 1, respectively. However, for those figures with τ = 1 (Figures 2, S2, S4 and S6), there is no obvious trend for the number of the blue points. Compared to Figure 1 (τ = 0.6), the agreements between the four point estimates and the true values of γ in Figure 2 (τ = 1) are worse, which can also be seen in other figures (Figures S1, S3 and S5 vs. Figures S2, S4 and S6, respectively). Observing Figure 2, we find that the four point estimation methods may underestimate γ when τ = 1. Finally, these four point estimation methods perform better for the quantitative trait than for the qualitative trait (Figures 1, 2, S1 and S2 vs. Figures S3-S6, respectively), and when the sample size increases (Figures S1, S2, S5 and S6 vs. Figures 1, 2, S3 and S4, respectively).     Table  3, we observe that the EPs of the PF method are generally smaller than, or equal to, those of the Fieller's method, except for the quantitative trait with = 500, = 0.4 or 1, and = 1, and the qualitative trait with = 500 or 2000, = 0.4 or 1, and = 1. However, the NPs of the PF method are always smaller than, or equal to, those of the Fieller's method. Note that when we use the PF and Fieller's methods to calculate the CIs of , we need to truncate the CIs by the interval [0, 2]. As such, compared to the Fieller's method, the PF method can get shorter CIs, which means that the PF method reduces the possibility of the truncated CIs being the noninformative intervals. On the other hand, if the CIs before the truncation are disjoint from the interval [0, 2] , the PF method will increase the possibility that the truncated CIs are empty sets, which is the reason why the PF method may have bigger EPs than the Fieller's method in some scenarios. In addition, all the DPs of the PF method are equal to 0. This is because we consider the penalty parameter =    Table 3, we observe that the EPs of the PF method are generally smaller than, or equal to, those of the Fieller's method, except for the quantitative trait with n = 500, η = 0.4 or 1, and τ = 1, and the qualitative trait with n = 500 or 2000, η = 0.4 or 1, and τ = 1. However, the NPs of the PF method are always smaller than, or equal to, those of the Fieller's method. Note that when we use the PF and Fieller's methods to calculate the CIs of γ, we need to truncate the CIs by the interval [0, 2]. As such, compared to the Fieller's method, the PF method can get shorter CIs, which means that the PF method reduces the possibility of the truncated CIs being the noninformative intervals. On the other hand, if the CIs before the truncation are disjoint from the interval [0, 2], the PF method will increase the possibility that the truncated CIs are empty sets, which is the reason why the PF method may have bigger EPs than the Fieller's method in some scenarios. In addition, all the DPs of the PF method are equal to 0. This is because we consider the penalty parameter λ = Z 2 1−α/2 /4, and the CIs derived by the PF method are always continuous. With increase of the sample size, the NPs of the PF and Fieller's methods and the DPs of the Fieller's method become smaller. Moreover, under the same simulation settings, the NPs of both methods, and the DPs of the Fieller's method, for the quantitative trait are less than those for the qualitative trait. Under the situation that τ = 0.6, when η changes from 0, 0.4 to 1 and other parameters are kept unchanged, the EPs of both methods have no obvious trends, while the NPs of both methods and the DPs of the Fieller's method generally become larger. As for τ = 1, when η changing from 0, 0.4 to 1 and other parameters being fixed, the EPs of the PF method appear larger except for the quantitative trait and n = 2000, while the DPs of the Fieller's method are relatively stable, and the NPs of the PF and Fieller's methods show a trend of first increasing and then decreasing on most occasions. On the other hand, when other parameters are fixed, the EPs and NPs of the PF and Fieller's methods with τ = 0.6 are smaller than those with τ = 1 in most cases, and the DPs of the Fieller's method with τ = 0.6 are larger than or equal to those with τ = 1. The CPs, W mean and W median of the GBN, GBU, PF and Fieller's methods are displayed in Table 4, and the corresponding W sd and W iqr are given in Table 5. Table 4 demonstrates that, for the quantitative trait, the CPs of the GBN, GBU and Fieller's methods are controlled around 95%. However, when n = 500, η = 1 and τ = 1, the CP of the PF method is underestimated (87.8%). As the sample size increases to 2000 and other parameters remain unchanged, the CP of the PF method is 96.6%. For the qualitative trait, when n = 500, the CPs of the GBN, GBU and PF methods are underestimated in most situations. With the increase of the sample size to 2000, the CPs of these three methods generally increase to be around 95%, but the CPs when η = 1 and τ = 1 are still underestimated. Thus, for this simulation setting, we conduct an additional simulation study with larger sample sizes (3000 and 4000), and the corresponding results are presented in Table S1. It is shown in Table S1 that the CPs of these three methods are closer to 95% when the sample size continues to increase. This is explainable by the fact that qualitative traits generally require larger samples to achieve the same CPs than quantitative traits. In addition, we can see from Table 4 that the Fieller's method has higher CPs under various simulation settings for the qualitative trait. However, according to Table 3, when the trait is qualitative, the NPs of the Fieller's method are relatively high, which means that many CIs obtained by the Fieller's method are the noninformative intervals (i.e., [0,2]). This may explain why the CPs of the Fieller's method are on the high side. Further, from Tables 4 and 5, the W mean , W median , W sd and W iqr of the GBN and GBU methods are smaller than those of the PF and Fieller's methods in most situations. The GBN method has the smallest W mean , W median and W iqr in most cases, and it also has the smallest W sd under all the simulated settings. As can be seen from Table 4, when the trait is qualitative and n = 500, the W median 's of the Fieller's method are all 2, which indicates that in this case, more than half of the CIs based on the Fieller's method are the noninformative intervals. This is consistent with the results of the NPs in Table 3. When the sample size increases, or the trait turns from qualitative into quantitative, the W mean 's and W median 's of the four interval estimation methods greatly decrease. However, for the W sd and W iqr , there are different trends in some situations. For example, when the trait is qualitative, the W sd 's and W iqr 's of the four methods become larger in most cases as the sample size increases. Note that the widths of the intervals obtained by the four methods are closer to 2 and the corresponding variation will be smaller when n = 500. With the sample size increasing, the widths of the intervals gradually decrease and the corresponding variation appears larger, which may cause the bigger W sd and W iqr .  In the case of τ = 0.6, the four methods have larger W mean 's and W median 's in most cases when η changes from 0, 0.4 to 1, while for the scenario of τ = 1, the four methods show a trend of first increasing and then decreasing on most occasions, except that the W mean 's and W median 's of the Fieller's method are gradually larger for the qualitative trait. When the trait is quantitative and τ = 0.6, the W sd 's and W iqr 's of the four methods become larger with η increasing from 0, 0.4 to 1, irrespective of the sample size. When the trait is qualitative, n = 500 and τ = 0.6, as η is bigger, the W sd 's and W iqr 's of the four methods generally are smaller, while when n = 2000, the W sd 's of the four methods and the W iqr 's of the GBN and GBU methods are relatively stable, and the W iqr 's of the PF and Fieller's methods generally become larger. For the quantitative trait with n = 500 and τ = 1, with the increase of η, the W sd 's and W iqr 's of the GBN, GBU and Fieller's methods appear smaller and those of the PF method are larger in most situations, while in the case of n = 2000, the four methods usually have larger W sd 's and W iqr 's. When the trait is qualitative and τ = 1, with η increasing, the W sd 's and W iqr 's of the GBN and GBU methods present a tendency of first decreasing and then increasing on most occasions, while those of the PF method are larger in most cases, and those of the Fieller's method become smaller, irrespective of the sample size. On the other hand, when other parameters are fixed, the W mean 's and W median 's of the four methods with τ = 0.6 are smaller than those with τ = 1, except for the W mean 's of the GBN, GBU and PF methods and the W median 's of the GBN and GBU methods for the quantitative trait with n = 500 and η = 1, and those for the qualitative trait with n = 500 or 2000, and η = 1. Under the scenarios where η is kept unchanged, the W sd 's and W iqr 's of the GBN, GBU and Fieller's methods with τ = 0.6 are generally larger than those with τ = 1 for the quantitative trait with n = 500, and the qualitative trait with n = 500 or 2000, while there are different trends for the quantitative trait with n = 2000. In addition, the W sd 's and W iqr 's of the PF method with τ = 0.6 generally are smaller than those with τ = 1, when other parameters are fixed. Figures 3, 4 and S7-S12 are the scatter plots of the widths of the 95% HPDIs or CIs obtained by the four interval estimation methods (GBN, GBU, PF and Fieller) against the true values of γ under different simulation settings. We can clearly observe the distributions of the widths of the HPDIs or CIs through these figures. For example, Figures 3 and 4 are the scatter plots of the widths of the HPDIs or CIs against the true values of γ for the quantitative trait with n = 500, and τ = 0.6 and 1, respectively. In each figure, subplots (a)-(d) (four subplots in the first row) are respectively the scatter plots for the GBN, GBU, PF and Fieller's methods with η = 0; subplots (e)-(h) (four subplots in the second row) and subplots (i)-(l) (four subplots in the third row) are the corresponding scatter plots with η = 0.4 and 1, respectively. It can be seen from the four subplots in the same row of each figure that the distributions of the widths of the HPDIs for the GBN and GBU methods are similar, and both have smaller dispersions than those of the CIs for the PF and Fieller's methods. Furthermore, these figures display that the distributions of the interval widths for the PF and Fieller's methods are greatly more skewed towards 2 than the GBN and GBU methods. We respectively compare subplots (a), (e) and (i) for the GBN method with subplots (b), (f) and (j) for the GBU method and find that the dispersions of the widths of the HPDIs for the GBN method are slightly smaller than the GBU method. Additionally, subplots (c), (g) and (k) for the PF method, and subplots (d), (h) and (l) for the Fieller's method, show that the PF and Fieller's methods may yield empty sets or noninformative intervals (displayed by the blue points), and the Fieller's method may also get discontinuous intervals (shown by the orange points). By comparing the subplots in different rows of each figure (Figures 3 and S7) when the trait is quantitative and τ = 0.6, the dispersions of the widths of the HPDIs or CIs become slightly larger as η changing from 0, 0.4 to 1, and it can also be seen from Figure 3 that the distributions of the interval widths are a little more skewed towards 2. On the other hand, when the trait is qualitative with τ = 0.6 ( Figures S9 and S11), there are no obvious trends in the dispersions of the interval widths, except that their distributions are more skewed towards 2. However, under the situation that τ = 1 (Figures 4, S8, S10 and S12), the points in these figures become less discrete in most cases when η increases, and the overall widths of the four interval estimation methods also somewhat decrease, except for the scenarios where the trait is quantitative and n = 2000, and the trait is qualitative and n = 500. Further, by comparing the figures for different values of τ (Figures 3, S7, S9 and S11 vs. Figures 4, S8, S10 and S12, respectively), it can be found that the overall widths of the HPDIs or the CIs obtained by the four interval estimation methods with τ = 0.6 are generally smaller than those with τ = 1, except for those with η = 1. Lastly, as the trait turns from qualitative into quantitative ( Figures S9-S12 vs. Figures 3, 4, S7 and S8, respectively) or the sample size increases (Figures 3, 4, S9 and S10 vs. Figures S7, S8, S11 and S12, respectively), the performances of the four interval estimation methods are greatly improved.

Application to MCTFR Data
The MCTFR Genome-Wide Association Study of Behavioral Disinhibition is a familybased epidemiological study of substance abuse and related psychopathology. The dataset can be made available from the database of Genotypes and Phenotypes with accession numbers 86747-6 and 95621-5 (https://www.ncbi.nlm.nih.gov/projects/gap/cgibin/study.cgi?study_id=phs000620.v1.p1, accessed on 5 January 2022). The dataset includes 2183 families and 7377 participants (3831 female subjects and 3546 male subjects). Among them, only 5960 subjects have both the phenotypic and genotypic data, while others do not have phenotypic data or do not have genotypic data. There are five quantitative traits included in this dataset: the nicotine composite score, the alcohol consumption composite score, the alcohol dependence composite score, the illicit drug composite score and the non-substance use related behavioral disinhibition composite score. To avoid the influence of family structure on the results, we exclude offspring from the real data application. At the same time, we only need the information of female subjects, so we also exclude male subjects from the analysis. Meanwhile, 12,354 SNPs on the X chromosome are included in the dataset. We use the following quality control

Application to MCTFR Data
The MCTFR Genome-Wide Association Study of Behavioral Disinhibition is a familybased epidemiological study of substance abuse and related psychopathology. The dataset can be made available from the database of Genotypes and Phenotypes with accession numbers 86747-6 and 95621-5 (https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi? study_id=phs000620.v1.p1, accessed on 5 January 2022). The dataset includes 2183 families and 7377 participants (3831 female subjects and 3546 male subjects). Among them, only 5960 subjects have both the phenotypic and genotypic data, while others do not have phenotypic data or do not have genotypic data. There are five quantitative traits included in this dataset: the nicotine composite score, the alcohol consumption composite score, the alcohol dependence composite score, the illicit drug composite score and the non-substance use related behavioral disinhibition composite score. To avoid the influence of family structure on the results, we exclude offspring from the real data application. At the same time, we only need the information of female subjects, so we also exclude male subjects from the analysis. Meanwhile, 12,354 SNPs on the X chromosome are included in the dataset. We use the following quality control criteria to filter the SNPs [48,49]: (1) genotype call rate being less than 99%, (2) MAF being smaller than 1 × 10 −5 , (3) individual call rate being below 99%, and (4) the p value of the Hardy-Weinberg equilibrium test being less than 1 × 10 −6 . Finally, 1994 female subjects and 12,342 SNPs on the X chromosome are utilized to conduct real data analysis. Since we estimate the degree of the skewness of the XCI based on genes, we first need to find the genes which each SNP belongs to. Based on the GRCH38 (Genome Reference Consortium Human Genome Build 38, https://uswest.ensembl.org/, accessed on 25 February 2022) reference, we use the "getBM" function in the R package "biomaRt" to match the SNPs to the genes on the X chromosome [45]. As such, we find 733 matched genes, while there are some genes containing only a single SNP in the dataset. As there have been several methods available to estimate the degree of the skewness of the XCI for a single SNP, we exclude genes consisting of only one SNP. Therefore, only 493 genes are included in the subsequent analysis.
Note that estimating γ requires the genes on the X chromosome to be associated with the traits. So, we need to test if the associations between the genes and the traits exist before using our proposed methods to estimate the degree of skewness of the XCI. Notice that the five traits in the MCTFR dataset do not follow normal distributions; therefore, we use the rank-based inverse normal transformation to transform the trait data [50]. Further, to adjust the effects of other variables, we incorporate two covariates, age and year of birth, into the application [48]. Due to the fact that we only use female subjects, we still apply the adaptive sum test proposed by Iuliana et al. [35] to test for the association between each gene and each trait. Unlike other multi-locus association analysis methods, when there are both rare and common variants in a gene, the adaptive sum test still maintains high test power. We set the significance level to be α = 0.05/(5 × 493) = 2.03 × 10 −5 based on the Bonferroni correction. After identifying the genes associated with the traits, we calculate the four point estimates of γ (γ GBN ,γ GBU ,γ GPF andγ GF ), and then use the GBN, GBU, PF and Fieller's methods to obtain the corresponding HPDIs or CIs.
We finally identify only one gene, TMEM47, statistically significantly associated with the alcohol dependence composite score (p value = 2.32 × 10 −6 ). There are two SNPs (rs10522027 and rs5928615) included in the gene. The estimated MAFs of these two SNPs are 0.1407 and 0.0998, respectively, which means that both SNPs only contain common variants. It has been confirmed that TMEM47 is located in the NC_000023.11 region and includes three exons. Studies have shown that the gene is expressed in the bladder, adipose and 23 other tissues and found that the overexpression of TMEM47 may induce resistance in patients to certain chemotherapy drugs [51,52]. The four point estimates (γ GBN ,γ GBU ,γ GPF andγ GF ) of γ for the gene are 0.4703, 0.4547, 0.4816 and 0.4847, and the 95% HPDIs or CIs derived by the GBN, GBU, PF and Fieller's methods are (0.0023, 1.2380), (0.0337, 1.3083), (0.0562, 1.2410) and (0.0557, 1.3896), respectively. That is to say, the point estimates are all less than 0.5, while the 95% HPDIs or CIs all contain 1, which means that the XCI pattern for TMEM47 on the alcohol dependence composite score may be the XCI-R or the XCI-E. By comparing the interval widths of these four interval estimation methods, we find that the width of the CI obtained by the PF method is the shortest, followed by the HPDI obtained by the GBN method, and the longest is the CI yielded by the Fieller's method.

Discussion
In this paper, we propose four point estimates (γ GBN ,γ GBU ,γ GPF andγ GF ) and four interval estimation methods (GBN, GBU, PF and Fieller) of the degree of the skewness of the XCI for a gene (i.e., γ). Among the point estimates,γ GF is constructed by truncating the ratio of the two regression coefficients by the interval [0, 2]. And,γ GPF is obtained by choosing the penalty parameter λ = Z 2 1−α/2 /4, and respectively correcting the denominator and the numerator, which is also truncated by [0,2]. Both theγ GBN andγ GBU are developed, based on the Bayesian theory, by considering the prior information of γ ∈ [0, 2], and the corresponding prior distributions of γ are respectively a truncated normal distribution and a uniform distribution. Use ofγ GBN andγ GBU can avoid the extreme point estimates of γ (0 or 2) occurring. Among the interval estimation methods, the Fieller's method has been widely used to construct the CIs of a ratio estimate. The PF method can always get the bounded CIs by choosing an appropriate penalty parameter. The GBN and GBU methods calculate the HPDIs of the samples randomly chosen from the approximate posterior distributions of γ as the credible intervals, which can avoid empty sets, noninformative intervals (i.e., [0,2]) and discontinuous intervals to occur. We conducted extensive simulation studies to compare their performances, by simulating different types of traits (quantitative and qualitative), different sample sizes (n = 500 and 2000), different proportions of rare variants among all the SNPs considered (η = 0, 0.4 and 1), and different proportions of the SNPs with positive effects among all the SNPs considered (τ = 0.6 and 1). The simulation results showed that there may exist some extreme point estimates forγ GPF andγ GF , especially when the sample size is small or the proportion of rare variants is high. The least MSE, in most situations, is derived fromγ GBN , and the MSEs ofγ GBN andγ GBU are smaller than those ofγ GPF andγ GF . As for the interval estimation, the CIs derived by the Fieller's method may be empty sets, noninformative intervals and discontinuous intervals. Although the PF method can avoid discontinuous intervals, the resulting CIs can be empty sets and noninformative intervals. In addition, most of the CPs of the GBN and GBU methods can be controlled around 95%, and a larger sample size is required only when the trait is qualitative and all the SNPs are rare variants. For qualitative traits, the CPs of the PF method appear a little low when the sample size is relatively small. However, the CPs of the Fieller's method seem to be well controlled, which is due to the large proportion of noninformative intervals. The GBN method has the smallest W mean , W median and W iqr in most situations, and the least W sd under all the simulation settings. Therefore, we recommend usingγ GBN and the GBN method to estimate the degree of the XCI skewing in practical applications.
On the other hand, concerning the simulation settings and the simulation results, we further discuss the following issues. Firstly, we consider the influence of the proportion of rare variants (η) and the proportion of the SNPs with positive effects (τ) among all the SNPs in the gene under study on the estimation results. When τ = 0.6 and other parameters are fixed, the proportions of the extreme values (0 and 2) forγ GPF andγ GF with η = 0 are generally less than those with η = 1, while they have no obvious trends for other cases of different values of η and τ. In general, the MSEs of the four point estimates generally become larger as η changes from 0, 0.4 to 1 and other parameters are kept unchanged. The four point estimates with τ = 0.6 always have smaller MSEs than τ = 1. The changing trends of the EPs, NPs and DPs of the PF and Fieller's methods with the increase of η are related to τ. Furthermore, the EPs and NPs of the PF and Fieller's methods with τ = 0.6 generally are smaller than τ = 1, while the DPs of the Fieller's method with τ = 0.6 are larger than or equal to those with τ = 1. On the other hand, in the case of τ = 0.6, the four interval estimation methods have larger W mean 's and W median 's in most cases with η changing from 0, 0.4 to 1, while for the scenario of τ = 1, those of the four methods show a trend of first increasing and then decreasing on most occasions. The changing tendencies of the W sd 's and W iqr 's of the four methods, with η increasing, are affected by the trait type, n and τ. When other parameters are kept unchanged, the W mean 's and W median 's of the four methods with τ = 0.6 are smaller than those with τ = 1 in most cases. Besides this, the findings, by comparing the W sd 's and W iqr 's of the GBN, GBU and Fieller's methods for τ = 0.6 with those for τ = 1, are related to the trait type and n, while the W sd 's and W iqr 's of the PF method with τ = 0.6 are generally smaller than those with τ = 1. Secondly, to better evaluate the performances of the proposed methods, we set the degrees of the XCI skewing at all the SNPs in the gene to be the same in our simulation studies. For example, when we calculate the MSEs of the point estimates and the CPs of the HPDIs or the CIs, a single true value of γ for each replicate is required. However, note that there may be different degrees of the XCI skewing at different SNPs, and, actually, we can also consider this issue in our simulation studies, although we have no appropriate evaluation indexes to assess the performances of the proposed methods for this situation. Finally, when we simulate quantitative traits, the random error ε i is generated from the standard normal distribution, where the standard deviation (σ) is equal to 1. To further illustrate the effect of different values of σ on the estimation results, we conducted additional simulation studies with n = 2000 and assume that ε i follows N(0, 4), where σ = 2. The corresponding results are presented in Tables S2-S4 and Figures S13-S16. As can be seen from these tables and figures, the Bayesian methods still have obvious advantages in both the point estimation and the interval estimation. Further, the four point estimation methods, and the four interval estimation methods with σ = 2, perform worse than σ = 1.
We applied the proposed methods to the MCTFR data and identified a gene, TMEM47, which is statistically significantly associated with the alcohol dependence composite score. However, although the four point estimates of γ for the gene TMEM47 on the alcohol dependence composite score are all smaller than 0.5, the corresponding 95% HPDIs or CIs all contain 1, which means that the XCI pattern for this gene may not be the XCI-S. Further, we observed that the width of the CI obtained by the PF method is the shortest, followed by the HPDI obtained by the GBN method, and the longest was the CI yielded by the Fieller's method. However, it should be noted that the PF method may not control the CP well (e.g., Table S3).
Last, but not least, there are still some issues in our proposed methods which need to be discussed. Firstly, we would like to further discuss the effect of the truncation by the interval [0, 2] on the point estimation and the interval estimation of γ. When we use thê γ GPF andγ GF to estimate γ, both of them are truncated by [0,2]. If the point estimates before the truncation (γ * andγ) lie outside [0, 2],γ GPF andγ GF become the extreme values (0 or 2). Correspondingly, when using the PF and Fieller's methods to construct the CIs of γ, it is easy to obtain empty sets or noninformative intervals. On the contrary, the Bayesian methods can avoid extreme point estimates, empty sets and noninformative intervals by specifying the appropriate prior distributions of γ and making full use of the constraint condition of γ ∈ [0, 2]. In addition, the extreme point estimate of 0 (2) means that the XCI is completely skewed towards the minor alleles (major alleles) at all the SNPs in a gene. However, these phenomena are not common in practice [2]. Meanwhile, it should be noted that empty sets and noninformative intervals are not informative, and the discontinuous CIs are also not useful, because the discontinuous CIs cannot be clearly explained in practice. Secondly, since the XCI patterns at different SNPs may be different, our estimatedγ is just the mean degree of the skewness of the XCI over all the SNPs in the gene under study, and we cannot obtain the degree of the skewness of the XCI for each SNP in this gene. Meanwhile, in the process of estimating γ, the target allele is the minor one at each SNP, and it is not possible to distinguish the disease allele from the normal allele at each SNP. Therefore, we can only identify whether or not the XCI of the gene is skewed towards the minor alleles, but it is not possible to know whether the XCI is skewed towards the disease alleles or the normal alleles. Thirdly, the proposed Bayesian methods need to specify the prior distributions of all the unknown parameters in advance, and the selection of the prior distributions may have a certain impact on the results. For simplicity, we only considered two possible prior distributions for γ, and one prior distribution for each of the other unknown parameters. However, the prior distributions of these parameters are usually unknown, and we cannot guarantee that the weak prior distributions we used are the most appropriate. We provide an R package named GEXCIS, which is publicly available at https://github.com/Meng-KaiLi/GEXCIS (accessed on 30 April 2022), and can be used to estimate the degree of the skewness of the XCI for genes through the proposed methods in this paper. This R package also allows researchers to specify the prior distribution of each unknown parameter from their own research backgrounds. Fourthly, the Bayesian methods use the HMC algorithm for the sampling, which is not affected by the correlation between unknown parameters. Therefore, to improve computational efficiency, we assumed that all the unknown parameters are independent. However, the Bayesian methods, taking the correlation between the parameters into account, should have better performance, which is our future work. Fifthly, if the HPDIs or CIs we get contain 1, which means that the XCI pattern for the gene is the XCI-R or the XCI-E, our proposed methods cannot distinguish them. Therefore, in our future study, we will consider including males' information to distinguish the XCI-R from the XCI-E. Finally, the proposed methods are only applicable to independent female subjects, and we will extend them in future so that they could accommodate the family data.

Conclusions
We propose four point estimates and four interval estimation methods to estimate γ of genes. Among the four point estimates,γ GF may have the extreme point estimates, andγ GPF can only reduce the occurrence of the extreme point estimates equal to 2, whilê γ GBN andγ GBU can avoid the extreme point estimates occurring. As for the four interval estimation methods, the Fieller's method may derive empty sets, discontinuous intervals and noninformative intervals, and the PF method can avoid the occurrence of discontinuous intervals and get less noninformative intervals, while the GBN and GBU methods do not yield these three types of the intervals. However, it should be noted that through these proposed methods, we cannot obtain the degree of the skewness of the XCI for each SNP in the gene, and cannot know whether the XCI is skewed towards the disease alleles or the normal alleles. In summary, the point estimates obtained by the GBN method always have the least MSE, and the HPDIs of the GBN method generally have the shortest width and the lowest variation, so we recommend using the GBN method in practical applications.

Conflicts of Interest:
The authors declare no conflict of interest.