Likelihood Ratio Test for Excess Homozygosity at Marker Loci on X Chromosome.

The assumption of Hardy-Weinberg equilibrium (HWE) is generally required for association analysis using case-control design on autosomes; otherwise, the size may be inflated. There has been an increasing interest of exploring the association between diseases and markers on X chromosome and the effect of the departure from HWE on association analysis on X chromosome. Note that there are two hypotheses of interest regarding the X chromosome: (i) the frequencies of the same allele at a locus in males and females are equal and (ii) the inbreeding coefficient in females is zero (without excess homozygosity). Thus, excess homozygosity and significantly different minor allele frequencies between males and females are used to filter X-linked variants. There are two existing methods to test for (i) and (ii), respectively. However, their size and powers have not been studied yet. Further, there is no existing method to simultaneously detect both hypotheses till now. Therefore, in this article, we propose a novel likelihood ratio test for both (i) and (ii) on X chromosome. To further investigate the underlying reason why the null hypothesis is statistically rejected, we also develop two likelihood ratio tests for detecting (i) and (ii), respectively. Moreover, we explore the effect of population stratification on the proposed tests. From our simulation study, the size of the test for (i) is close to the nominal significance level. However, the size of the excess homozygosity test and the test for both (i) and (ii) is conservative. So, we propose parametric bootstrap techniques to evaluate their validity and performance. Simulation results show that the proposed methods with bootstrap techniques control the size well under the respective null hypothesis. Power comparison demonstrates that the methods with bootstrap techniques are more powerful than those without bootstrap procedure and the existing methods. The application of the proposed methods to a rheumatoid arthritis dataset indicates their utility.

Introduction

Association analysis is a useful tool to map disease loci by using markers on autosomes based on family data and case-control data [1-9]. There has been an increasing interest of exploring the association between diseases and markers on X chromosome and the effect of the departure from Hardy-Weinberg equilibrium (HWE) on association analysis on X chromosome [10-17]. Note that there are two hypotheses of interest regarding the X chromosome: (i) the frequencies of the same allele at a locus in males and females are equal and (ii) the inbreeding coefficient in females is zero (without excess homozygosity) in X-specific quality control [18, 19]. As such, excess homozygosity in females and significantly different minor allele frequencies between males and females are used to filter X-linked variants [20, 21]. The inbreeding coefficient is generally estimated by functions of excess homozygosity [22, 23], which may be caused by population substructure, consanguineous mating or factors like null alleles [24, 25]. Overall and Nichols developed an approach to distinguish population substructure and consanguinity by using multilocus genotype data [24]. On the other hand, Zheng et al. proposed two test statistics to test for the equality of the frequencies of the same allele in males and females and the zero inbreeding coefficient in females on X chromosome, respectively [14]. However, they only focused on association analysis on X chromosome and the type I error rates and powers of these two test statistics have not been studied yet. Further, there is no existing method to simultaneously detect both of the issues till now.

Therefore, in this article, we first combine two test statistics proposed in zheng et al. [14] and suggest Z 0 to simultaneously test for (i) the equality of the frequencies of the same allele in males and females and (ii) the zero inbreeding coefficient on X chromosome based on the collected sample. For the purpose of improving the test power for both (i) and (ii), a novel likelihood ratio test on X chromosome is proposed. We write out the likelihood functions of the collected sample under the null hypothesis and alternative hypothesis at a single locus on X chromosome, respectively. Next, we obtain the maximum likelihood estimates (MLEs) of the unknown parameters by expectation-maximization (EM) algorithms [26] and construct the corresponding likelihood ratio test (LRT0) statistic to test for both (i) and (ii). If the null hypothesis is statistically rejected, we further conduct two hypothesis testing issues to find the underlying reasons why the null hypothesis is violated by proposing another two likelihood ratio tests LRT1 (for the equality of the frequencies of the same allele in males and females) and LRT2 (for excess homozygosity). Note that the size of LRT0 and LRT2 is conservative from our simulation study. As such, we use parametric bootstrap techniques to evaluate the validity and performance of LRT0 and LRT2, which are respectively denoted by LRT0 and LRT2. Moreover, we explore the effect of population stratification on the proposed tests. In addition, the root mean squared error (RMSE) and bias are used to assess the accuracy of the MLEs of the unknown parameters. Finally, the application of the proposed methods to a rheumatoid arthritis (RA) dataset indicates their utility.

Materials and Methods

Background and notations

Consider a biallelic marker locus on X chromosome with alleles M 1 and M 2. Let p and p be the frequencies of M 1 in males and females, respectively. As such, the frequencies of M 2 in males and females are q = 1 − p and q = 1 − p , respectively. In females, let ρ be the inbreeding coefficient, which is generally nonnegative [27-29]. Thus, the frequencies of three genotypes M 1 M 1, M 1 M 2 and M 2 M 2 in females can be expressed as follows:

To this end, there is no excess homozygosity in females when ρ = 0; excess homozygosity exists when ρ > 0. Note that p ≠ p may be true on X chromosome. So, we construct the null hypothesis denoted by H 0: p = p and ρ = 0 to test for both of the hypotheses (i) and (ii). If the null hypothesis is violated, we need to investigate which one of p ≠ p and ρ > 0 is true. As such, we have other two hypothesis testing issues with the null hypothesis being H 01: p = p and H 02: ρ = 0, respectively. It should be noted that X chromosome has the problem of X chromosome inactivation and dosage compensation [30], but we do not consider them in this section. The corresponding discussion can be found later (see the Discussion section).

Assume that n 1 and n 0 represent the numbers of males with alleles M 1 and M 2 in a collected sample, respectively; n 2, n 1 and n 0 denote the numbers of females with genotypes M 1 M 1, M 1 M 2 and M 2 M 2, respectively. Then, N = n 1 + n 0 and N = n 2 + n 1 + n 0 are respectively the numbers of males and females in the sample, and N = N + N is the sample size.

Existing methods Z 1 and Z 2 for H 01 (equality of the frequencies of the same allele in males and females) and H 02 (zero inbreeding coefficient), respectively

Zheng et al. [14] proposed the test statistic to test for H 01: p = p , where and are the estimates of p and p , , and are the estimates of the variances of and under H 01, respectively, with . Under H 01, Z 1 asymptotically follows the chi-square distribution with one degree of freedom when the sample size is large enough.

Weir and cockerham [31] introduced the disequilibrium coefficient in females . In other words, testing for Δ = 0 is equivalent to testing for ρ = 0. Hence, zheng et al. [14] further developed the following test statistic to test for H 02: ρ = 0, where , , and . Under H 02, Z 2 approximately follows the chi-square distribution with one degree of freedom when N is large enough. It should be noted that the test Z 2 has nothing to do with male individuals and thus only needs female individuals.

Z 0 test for both hypotheses (i) and (ii) of interest regarding the X chromosome

Zheng et al. [14] showed that, under H 0: p = p and ρ = 0, Z 1 and Z 2 are independent. However, they did not propose the corresponding test statistic for H 0. As such, we suggest the test statistic to test for H 0: p = p and ρ = 0. Under H 0, Z 0 asymptotically follows the chi-square distribution with the degrees of freedom being 2. Moreover, it should be noted that we can use to estimate the inbreeding coefficient ρ.

Likelihood ratio test for both hypotheses (i) and (ii) of interest regarding the X chromosome

To construct a likelihood ratio test (LRT) for H 0: p = p and ρ = 0, we give the likelihood function of the sample as follows: where θ = (p , p , ρ). Firstly, we use the following EM algorithm to estimate the unknown parameters p , p and ρ under the alternative hypothesis (H 1: p ≠ p or ρ > 0). Suppose that Y = (Y 1, Y 2, Y 3, Y 4, Y 5) = (n 1, n 0, n 2, n 1, n 0) denotes the observed data. (Y 1, Y 2, Y 3, Y 4, Y 5) can be augmented by splitting the third cell into two cells W 1 and W 2, which are unobservable random variables such that Y 3 = W 1 + W 2 for female homozygote M 1 M 1 and W 1 and W 2 follow the binomial distributions with success probabilities and , respectively, and by splitting the fifth cell into two cells W 3 and W 4, where Y 5 = W 3 + W 4 for female homozygote M 2 M 2 and W 3 and W 4 follow the binomial distributions with success probabilities and , respectively. Thus, the likelihood function of complete data (n 1, n 0, w 1, w 2, n 1, w 3, w 4) is: where the normalizing constant is omitted for brevity.

At the E-step, the Q function at iteration (k + 1) is constructed as where θ ( is the estimate of θ at iteration k.

At the M-step, the estimated value θ ( of θ at iteration (k + 1) can be obtained by maximizing the Q function with respect to θ. Therefore, the MLEs of p , p and ρ at iteration (k + 1) are respectively

Note that the MLE of p is the same for all the iterations, which is also the same as zheng et al. [14]. In the above expressions, where . Given the initial value θ (0) of θ, the above-mentioned two steps continue until the convergence criterion is satisfied. For example, the absolute differences between the estimates of the parameters at two consecutive iterations are all less than 10−7. The value of θ obtained at the last iteration is taken as the MLE of θ under H 1.

Note that p = p and ρ = 0 under H 0. Let p = p = p , the pooled allele frequency of M 1. Then, L(θ) in Eq (1) can be rewritten as

Thus, the MLE of p under H 0 is , the estimated pooled allele frequency of M 1. Let . Then, we can construct the following LRT to test for H 0 which asymptotically follows a chi-square distribution with the degrees of freedom being 2 when the null hypothesis holds.

Likelihood ratio test for equality of frequencies of the same allele in males and females

Once the null hypothesis (H 0: p = p and ρ = 0) is rejected based on the result of Eq (6), we further need to consider the following two tests H 01: p = p and H 02: ρ = 0. Note that under the null hypothesis H 01: p = p = p, ρ may not be zero and we need to estimate it. Let ϕ = (p, ρ) and q = 1 − p. Thus, the corresponding likelihood function of complete data is

We use the following EM algorithm to estimate ϕ under H 01. The corresponding formulas at iteration (k + 1) are as follows where and are respectively the MLEs of p and ρ at iteration (k + 1), and . E (w 1|n 2), E (w 2|n 2), E (w 3|n 0) and E (w 4|n 0) in the above expressions are similar to E (w 1|n 2), E (w 2|n 2), E (w 3|n 0) and E (w 4|n 0) in Eqs (2)–(5), just replacing , and in Eqs (2)–(5) by , and , respectively. Let . Then, we propose the following test statistic LRT1 to test for the null hypothesis H 01: p = p , which approximately follows a chi-square distribution with the degree of freedom being 1 under H 01.

Likelihood ratio test for inbreeding coefficient being zero

Note that under the null hypothesis H 02 : ρ = 0, p and p may be different from each other and we need to estimate them separately. Let ψ = (p , p ) and L(θ) in Eq (1) can be rewritten as

Then, the MLEs of p and p are and , respectively, which are the same as zheng et al. [14]. Let . As such, we develop the following test statistic to test for H 02 : ρ = 0 which asymptotically follows a chi-square distribution with the degree of freedom being 1 under H 02. Just like the Z 2 test statistic, LRT2 only uses female individuals in the sample because the terms based on male individuals in the numerator and the denominator of the fraction are the same, which can be reduced.

Likelihood ratio tests via parametric bootstrap for H 0 and H 02

It should be noted from our simulation results (see the Results section) that the simulated type I error rates of LRT0 and LRT02 respectively for H 0 and H 02 are too conservative. On the other hand, several studies showed that the likelihood ratio tests may typically not follow a chi-square distribution asymptotically [31, 32], and hence their exact distributions can be obtained by Monte Carlo simulation [33]. Accordingly, we make use of parametric bootstrap techniques to evaluate the size and power of these two methods. For convenience, we denote these methods via parametric bootstrap by LRT0 and LRT2, respectively. We begin by describing the implementation steps for LRT0 as follows:

For a collected sample of size N with N males and N females, calculate the value of LRT0;

Compute the estimated pooled allele frequency based on the sample as follows: ;

Based on , calculate the frequencies of three genotypes M 1 M 1, M 1 M 2 and M 2 M 2 in females under H 0 in the following: , and , respectively, where ;

According to and , regenerate the alleles of N males; based on , and , regenerate the genotypes of N females;

Calculate the value of LRT0 based on the new N males and N females, denoted by ;

Repeat Steps 4 and 5 B times, which results in B test statistics , , …, ;

The P-value of the original LRT0 can be estimated as

For LRT2, we can conduct the steps similar to those mentioned above. Firstly, after obtaining the value of LRT2, calculate the frequencies of three genotypes M 1 M 1, M 1 M 2 and M 2 M 2 in females under H 02 in the following: , and , respectively, with . The alleles of N males stay the same as the original sample and only regenerate the genotypes of N females according to , and . Then, carry out the similar procedures of Steps 4–7 and we can obtain the the estimated P-value of LRT2.

Software implementation

We have written the XHWE software with R (http://www.r-project.org), which includes the eight test statistics: LRT0, LRT0, LRT1, LRT2, LRT2, Z 0, Z 1 and Z 2. The R package named XHWE is available on CRAN (http://cran.r-project.org/web/packages/XHWE/). The initial values of p , p , p and ρ in the EM algorithms are taken to be n 1/N , (2n 2 + n 1)/(2N ), (n 1 + 2n 2 + n 1)/(N + 2N ) and 0.02, respectively. The convergence criterion is that the absolute differences between the estimates of the parameters at two consecutive iterations are all less than 10−7 for the LRT-type statistics. The default maximum number of iterations is 1000. The input data file is the standard pedigree data. The XHWE software only uses the founders with genotypes available in it and will analyze marker loci one by one. The software outputs the values of all the test statistics and the corresponding P-values. Also, the XHWE software outputs the estimates of all the parameters under both the null and alternative hypotheses for each test statistic. The parameter estimates under the alternative hypothesis for the LRT-type test statistics are the same. However, under the respective null hypotheses of the LRT-type test statistics, the estimates may be different. It should be noted that the estimates of p and p under the null hypothesis of H 02 in this article and those in zheng et al. [14] are the same, respectively. The output results will be automatically saved in the text file named “results.txt”.

Simulation settings

Simulation study is conducted to assess the performance of the proposed LRT0, LRT0, LRT1, LRT2, LRT2 and Z 0 test statistics and to compare them with the existing Z 1 and Z 2 under various simulation settings which are similar to those in zheng et al. [14]. The allele frequency p in males takes two values: 0.3 and 0.5. When p is fixed, the value of p in females is taken as p = p + ϵ, where ϵ = 0, ±0.04 and ±0.05. The inbreeding coefficient ρ in females is set at 0 to 0.1 in increment of 0.05. The sample size is taken as 800 and 1200 with the ratio r = N : N being 2:1, 1.5:1, 1:1, 1:1.5 and 1:2. As mentioned earlier, when p = p and ρ = 0, the size of all the eight test statistics is simulated; when p = p and ρ > 0, the size of LRT1 and Z 1 is gotten; when p ≠ p and ρ = 0, the size of LRT2, LRT2 and Z 2 is obtained. Otherwise, we simulate the corresponding powers. In addition, it should be noted that for the fixed sample size (800 or 1200) simulated above, the powers of all the three test statistics LRT2, LRT2 and Z 2 for H 02 : ρ = 0 are not so large, from our simulation results below. On the other hand, these three test statistics only use female individuals. As such, we further obtain the sample size N required for LRT2 to gain 80% simulated power and then simulate the size and powers of LRT2, LRT2 and Z 2 under this sample size. To investigate how population structure affects the proposed methods, we also consider the following population stratification model with two subpopulations in our simulation study. p = 0.3 (0.5), p = p + ϵ, ϵ = 0, ±0.04 and ±0.05 in the first (second) subpopulation and the ϵ values are respectively denoted by ϵ 1 and ϵ 2. Assume that ρ = 0 in each subpopulation, and the ratio of each subpopulation constructing the population is set to 0.5. The sample size is taken to be 1800, where each individual is a female or a male with equal probability. Note that under population stratification, the null hypothesis H 0: p = p and ρ = 0 is generally not true. Thus, we use the population stratification model to study the powers of the proposed methods. The significance level is fixed at 5% and 10000 replications are simulated under each simulation setting. For LRT0 and LRT2 via parametric bootstrap, B is set to be 1000. Finally, to compare the efficiency of the parameter estimates of the proposed EM algorithms with those in zheng et al. [14] for each simulation setting, we use the RMSEs and biases to assess the accuracy of the parameter estimates, where and , and β is the parameter which needs to estimate.

Results

Simulation results

Table 1 lists the simulated size of LRT0, LRT0, LRT1, LRT2, LRT2, Z 0, Z 1 and Z 2 under H 0 : p = p = p and ρ = 0 with N = 800 and 1200 and p = 0.3 and 0.5 for different values of r = N : N . According to the table, the size of LRT1, Z 0, Z 1 and Z 2 is close to the nominal 5% level, while the size of LRT0 and LRT2 is too conservative. However, after the parametric bootstrap technique, LRT0 and LRT2 stay close to the nominal 5% level.

Table 1

Simulated size (in %) of LRT0, LRT0, LRT1, LRT2, LRT2, Z 0, Z 1 and Z 2 under H 0 : p = p = p and ρ = 0 with N = 800 and 1200 for different values of r and p.

N	r	p	LRT0	LRT0b	LRT1	LRT2	LRT2b	Z 0	Z 1	Z 2
800	2:1	0.3	3.01	4.81	5.02	1.91	4.87	4.83	5.22	4.83
	2:1	0.5	2.98	4.97	5.02	2.28	4.99	5.13	5.19	5.13
	1.5:1	0.3	2.92	4.81	4.74	2.19	4.75	4.99	4.93	4.99
	1.5:1	0.5	3.07	4.93	4.83	2.22	4.98	5.02	4.99	5.02
	1:1	0.3	3.17	5.05	4.74	2.59	5.36	4.82	4.81	4.82
	1:1	0.5	3.30	4.99	5.33	2.40	5.20	5.16	5.34	5.16
	1:1.5	0.3	3.05	5.18	5.09	2.40	5.11	5.09	5.28	5.09
	1:1.5	0.5	3.39	5.18	5.03	2.49	5.34	5.19	5.07	5.19
	1:2	0.3	3.13	4.89	4.77	2.18	4.81	5.13	4.89	5.13
	1:2	0.5	3.12	4.85	4.65	2.32	5.13	5.23	4.78	5.23
1200	2:1	0.3	3.42	5.31	4.84	2.35	5.07	5.47	4.97	5.47
	2:1	0.5	3.45	5.38	4.76	2.48	5.15	5.38	5.01	5.38
	1.5:1	0.3	2.91	4.84	4.84	2.30	5.12	4.83	5.02	5.16
	1.5:1	0.5	3.36	5.31	5.29	2.45	5.35	5.38	5.42	5.38
	1:1	0.3	2.88	4.77	5.05	2.15	4.73	4.73	5.18	4.73
	1:1	0.5	3.04	5.07	5.22	2.01	4.79	4.94	5.30	4.94
	1:1.5	0.3	2.93	4.97	4.75	2.25	4.98	4.98	4.87	4.98
	1:1.5	0.5	3.08	4.83	5.06	2.24	4.86	4.87	5.12	4.87
	1:2	0.3	3.13	4.98	4.83	2.39	5.31	5.03	4.92	5.03
	1:2	0.5	3.24	5.06	4.86	2.54	5.38	5.05	4.91	5.05

Fig 1 gives the simulated powers of the eight test statistics against r under H 1 : p ≠ p and ρ > 0 for different values of ρ (0.05 and 0.1) and N (800 and 1200), having p = 0.3 and p = 0.35. It is shown in the figure that LRT0 is more powerful than LRT0 and Z 0, and LRT0 and Z 0 have the similar performance in power (Fig 1a-1d in the first row), regarded of the inbreeding coefficient ρ, the sample size N and the ratio r. LRT1 and Z 1 have almost the same performance in power (Fig 1e-1h in the second row). LRT2 has much more power than LRT2 and Z 2, and LRT2 is a little less powerful than Z 2 (Fig 1i-1l in the third row). The powers of LRT1 and Z 1 are not so affected by the different values of r, while LRT0, LRT0, Z 0, LRT2, LRT2 and Z 2 are more and more powerful with the number of female individuals increasing (r changing from 2:1 to 1:2) when other parameters are fixed. We also find that the powers of LRT0, LRT0, Z 0, LRT2, LRT2 and Z 2 appear great reaction to the different values of ρ when N is fixed. Specially, their powers under ρ = 0.1 (subplots in the second and fourth columns, respectively) are much larger than those under ρ = 0.05 (subplots in the first and third columns, respectively). However, the powers of LRT1 and Z 1 are almost not influenced by ρ. Further, it can be seen in Fig 1 that LRT0, LRT0 and Z 0 with two degrees of freedom (subplots in the first row) are much more powerful than LRT1, Z 1, LRT2, LRT2 and Z 2 with one degree of freedom (subplots in the second and third rows). This is because the true model is p ≠ p and ρ > 0. In addition, when the sample size changes from 800 (subplots in the first and second columns) to 1200 (subplots in the third and fourth columns), all the test statistics are much more powerful.

Fig 1

Simulated powers of LRT0, LRT0, LRT1, LRT2, LRT2, Z 0, Z 1 and Z 2 against r = N : N under H 1 : p ≠ p and ρ > 0 based on 10000 replicates with p = 0.3 and p = 0.35.

In the first column: ρ = 0.05 and N = 800; in the second column: ρ = 0.1 and N = 800; in the third column: ρ = 0.05 and N = 1200; in the fourth column: ρ = 0.1 and N = 1200. In the first row, the powers of LRT0, LRT0 and Z 0 for H 0 : p = p and ρ = 0; in the second row, the powers of LRT1 and Z 1 for H 01 : p = p ; in the third row, the powers of LRT2, LRT2 and Z 2 for H 02 : ρ = 0.

Fig 2 displays the simulated size/powers of the eight test statistics against r under H 02 : ρ = 0 for different values of p , having p = 0.3 and N = 1200. The results in the third row of the figure are the size of LRT2, LRT2 and Z 2, while those in the first and the second rows of the figure are the powers of LRT0, LRT0 and Z 0, and those of LRT1 and Z 1, respectively. It is shown in the figure that the size of LRT2 and Z 2 maintains close to the nominal 5% level, while LRT2 is too conservative. As for the tests for H 01 : p = p , LRT1 and Z 1 almost have the same simulated power just like Fig 1. On the other hand, the powers of LRT0, LRT0, Z 0, LRT1 and Z 1 are not so affected by the ratio r. However, their powers are greatly influenced by the absolute difference |ϵ| = |p − p |. Specifically, their powers under p = 0.25 and p = 0.35 are much larger than those under p = 0.26 and p = 0.34. In addition, when the simulation setting is fixed, LRT1 and Z 1 with one degree of freedom are a little more powerful than LRT0, LRT0 and Z 0 with two degrees of freedom, because the true model is p ≠ p and ρ = 0. By comparing Fig 2d (ρ = 0), Fig 1c (ρ = 0.05) and Fig 1d (ρ = 0.1) under N = 1200, p = 0.3 and p = 0.35, LRT0, LRT0 and Z 0 are more and more powerful with ρ increasing.

Fig 2

Simulated size/powers of LRT0, LRT0, LRT1, LRT2, LRT2, Z 0, Z 1 and Z 2 against r = N : N under H 02 : ρ = 0 based on 10000 replicates with p = 0.3 and N = 1200.

In the first column: p = 0.25; in the second column: p = 0.26; in the third column: p = 0.34; in the fourth column: p = 0.35. In the first row, the powers of LRT0, LRT0 and Z 0 for H 0 : p = p and ρ = 0; in the second row, the powers of LRT1 and Z 1 for H 01 : p = p ; in the third row, the size of LRT2, LRT2 and Z 2 for H 02 : ρ = 0.

Figs A–G in S1 File show the corresponding results under other simulation settings with p ≠ p and ρ > 0, which are similar to those in Fig 1. Figs H and I in S1 File plot the corresponding results under p = p and ρ > 0, and Figs J–L in S1 File give the corresponding results under p ≠ p and ρ = 0. The more details refer to S1 File.

Table 2 shows the simulated size of LRT2, LRT2 and Z 2 for H 02 : ρ = 0 for different values of p , having N = 0 under the sample sizes N required for LRT2 to obtain 80% simulated power. Table 3 lists the simulated powers under these sample sizes for different values of p , having ρ = 0.05 and 0.1. From Table 2, we can see that the type I error rates of LRT2, LRT2 and Z 2 are close to the nominal significance level of 5%. It is shown in Table 3 that the power of LRT2 attains to about 80%, and the difference in power between LRT2 and Z 2 is about 10%.

Table 2

Simulated size (in %) of LRT2, LRT2 and Z 2, having N = 0 and ρ = 0.

N f	p f	LRT2	LRT2b	Z 2
2500	0.20	2.27	5.18	4.84
	0.25	2.27	4.96	4.89
	0.30	2.43	5.15	5.13
	0.35	2.31	4.99	5.14
	0.40	2.25	4.96	4.70
	0.45	2.43	5.05	4.93
	0.50	2.54	5.03	5.06
	0.55	2.38	5.26	4.86
	0.60	2.49	5.11	5.10
650	0.20	2.14	4.98	4.67
	0.25	2.03	4.62	4.82
	0.30	2.53	4.91	5.20
	0.35	2.45	5.05	5.11
	0.40	2.18	4.72	4.54
	0.45	2.04	4.84	4.65
	0.50	2.48	4.98	5.14
	0.55	2.30	4.74	5.12
	0.60	2.27	5.01	4.77

Table 3

Simulated powers (in %) of LRT2, LRT2 and Z 2, having N = 0.

N f	ρ	p f	LRT2	LRT2b	Z 2
2500	0.05	0.20	68.2	79.0	69.7
	0.05	0.25	69.1	79.4	69.9
	0.05	0.30	69.4	79.9	70.2
	0.05	0.35	69.8	80.2	70.3
	0.05	0.40	70.0	80.6	70.4
	0.05	0.45	69.3	79.9	69.6
	0.05	0.50	70.6	80.8	71.5
	0.05	0.55	71.2	81.0	71.5
	0.05	0.60	70.2	80.6	70.5
650	0.10	0.20	68.1	78.9	70.3
	0.10	0.25	69.2	79.9	70.8
	0.10	0.30	69.9	80.9	71.0
	0.10	0.35	70.3	80.5	71.1
	0.10	0.40	71.1	81.1	71.8
	0.10	0.45	71.9	81.9	72.5
	0.10	0.50	71.1	81.6	72.1
	0.10	0.55	70.7	81.4	71.3
	0.10	0.60	70.8	81.9	71.6

Tables A–J in S1 File list the RMSEs and biases of the estimates of p , p , the pooled allele frequency p and ρ for different values of p , p , ρ, r and N. It should be noted that the estimate of p based on the EM algorithm is the same as zheng et al. [14]. Further, the estimates and of p and p based on the EM algorithms have the similar RMSEs and biases as those from zheng et al. [14], respectively. However, when we focus on the estimate of ρ, we find that although the biases of and based on the EM algorithms are larger than in zheng et al. [14] for some cases, the RMSEs of and are smaller than for all the simulation settings.

Table 4 displays the simulated size/powers of LRT0, LRT0, LRT1, LRT2, LRT2, Z 0, Z 1 and Z 2 under the population stratification model. When ϵ 1 = ϵ 2 = 0, the size of LRT1 and Z 1 is obtained. Further, note that the ratios of two subpopulations in the whole population are equal. As such, ϵ 1 = −ϵ 2 will also cause the size of LRT1 and Z 1. Under other simulation settings, we get the powers of the eight test statistics. To investigate whether or not the population stratification model causes excess homozygosity, we save the values of 10000 ρ estimates for each estimation method (, or ). Then, calculate the corresponding mean and standard deviation (SD), which are also listed in Table 4. The results show that the population stratification model indeed leads to the positive inbreeding coefficient (i.e., excess homozygosity), which is consistent with Overall and Nichols [24]. The mean values ( and ) using the EM algorithm are a little larger than proposed in zheng et al. [14], while and have less standard deviation. On the other hand, the size of LRT1 and Z 1 is close to the nominal significance level of 5%. The power of LRT0 is larger than LRT0 and Z 0, and LRT0 and Z 0 have the similar powers, irrespective of the ϵ 1 and ϵ 2 values. LRT1 and Z 1 have almost the same powers. LRT2 is much more powerful than LRT2 and Z 2, and the power of LRT2 is a little smaller than Z 2. If ϵ 1 is fixed and ϵ 2 is changed, the ρ estimate increases with ϵ 2 increasing, and hence LRT2, LRT2 and Z 2 are more and more powerful; if ϵ 2 is fixed and ϵ 1 is changed, the ρ estimate decreases with ϵ 1 increasing, and hence LRT2, LRT2 and Z 2 are less and less powerful. This may be caused by p being taken to be 0.3 and 0.5 in the first and second subpopulations, respectively.

Table 4

Mean and standard deviation (SD) of ρ estimates over 10000 replications, and simulated size/powers (in %) of LRT0, LRT0, LRT1, LRT2, LRT2, Z 0, Z 1 and Z 2 under population stratification model.

ϵ		ρ^_1		ρ^_01		ρ^_z		Power
ϵ 1 a	ϵ 2 b	Mean	SD	Mean	SD	Mean	SD	LRT0	LRT0b	LRT1	LRT2	LRT2b	Z 0	Z 1	Z 2
-0.05	-0.05	0.043	0.031	0.046	0.032	0.042	0.034	72.3	78.6	72.1	24.1	37.0	72.6	71.5	25.0
	-0.04	0.050	0.031	0.052	0.032	0.049	0.034	67.8	74.4	63.4	31.0	43.9	67.7	63.0	31.6
	0.00	0.067	0.032	0.068	0.032	0.067	0.033	54.8	63.4	22.7	51.7	65.2	55.3	22.5	52.5
	0.04	0.088	0.034	0.088	0.034	0.087	0.034	65.0	72.3	4.7	73.6	82.8	65.6	4.7	74.5
	0.05	0.093	0.035	0.093	0.035	0.093	0.035	70.9	76.3	4.1	78.7	86.2	71.3	4.1	79.0
-0.04	-0.05	0.039	0.029	0.041	0.030	0.037	0.033	60.5	69.0	63.2	19.2	28.5	60.9	63.2	19.8
	-0.04	0.045	0.030	0.046	0.031	0.043	0.033	57.6	63.8	50.5	25.6	37.6	57.6	49.9	26.2
	0.00	0.061	0.033	0.062	0.033	0.060	0.034	42.2	51.1	17.0	44.5	57.7	42.6	17.0	45.0
	0.04	0.081	0.033	0.081	0.033	0.081	0.034	57.4	64.0	4.2	65.9	77.9	58.1	4.3	66.4
	0.05	0.088	0.033	0.088	0.033	0.088	0.033	65.2	72.5	5.0	74.7	85.2	65.8	5.0	75.7
0.00	-0.05	0.029	0.027	0.030	0.027	0.025	0.033	21.8	29.0	23.0	11.0	19.0	22.5	22.9	12.1
	-0.04	0.030	0.027	0.031	0.027	0.026	0.033	18.9	24.6	17.3	12.5	20.0	19.6	17.2	12.9
	0.00	0.043	0.029	0.043	0.029	0.041	0.032	17.7	24.0	4.3	23.1	34.9	18.2	4.3	23.9
	0.04	0.058	0.032	0.058	0.032	0.058	0.033	42.3	50.6	19.0	40.0	54.1	43.1	19.0	40.6
	0.05	0.063	0.032	0.063	0.032	0.063	0.033	49.4	57.4	23.3	46.2	59.5	50.1	23.5	46.9
0.04	-0.05	0.019	0.023	0.020	0.023	0.012	0.032	6.7	9.4	6.8	4.8	10.0	7.6	6.9	6.1
	-0.04	0.022	0.024	0.022	0.025	0.014	0.034	6.3	9.0	4.7	7.1	11.7	7.0	4.7	7.7
	0.00	0.031	0.028	0.031	0.028	0.026	0.034	18.0	23.7	17.5	11.8	20.7	18.9	17.5	12.8
	0.04	0.041	0.031	0.041	0.031	0.039	0.034	52.1	61.9	51.8	21.2	32.0	52.7	52.1	21.5
	0.05	0.046	0.031	0.047	0.031	0.045	0.034	61.7	69.5	58.6	27.1	38.9	62.1	58.8	27.6
0.05	-0.05	0.018	0.023	0.018	0.023	0.008	0.034	4.1	7.2	4.9	4.5	8.4	5.2	5.0	5.5
	-0.04	0.020	0.023	0.020	0.023	0.011	0.034	6.0	8.4	5.2	4.2	10.1	7.5	5.2	5.9
	0.00	0.028	0.028	0.028	0.028	0.022	0.035	22.1	28.2	24.1	10.4	16.9	22.9	24.2	11.2
	0.04	0.036	0.028	0.037	0.028	0.034	0.032	59.5	65.7	61.6	17.0	27.1	60.0	61.9	17.7
	0.05	0.042	0.031	0.043	0.030	0.040	0.033	68.1	74.1	67.2	22.6	32.1	68.3	67.6	23.4

a ϵ in the first subpopulation.

b ϵ in the second subpopulation.

Application to RA data

We apply the proposed methods to the RA dataset from North American Rheumatoid Arthritis Consortium for studying their practicability, which is available from Genetic Analysis Workshop 15. In this dataset, there are 1217 families. Note that many individuals’ genotypes are missing. On the other hand, to obtain a sample of which all the individuals are independent, we only select the available founders in this dataset, which results in a sample composed of 369 founders (N = 112 and N = 257) in the analysis. 293 SNP markers on X chromosome for each founder are included in this application. The significance level is fixed at α = 5%. Table 5 gives the corresponding results based on the P-values of LRT0, LRT1, LRT2, Z 0, Z 1 and Z 2. From Table 5, LRT0 identified 6 loci which Z 0 did not identify, and Z 0 identified 4 additional loci. One locus is detected by LRT1 that is not found by Z 1, and 4 additional loci are detected by Z 1. There are 12 loci identified by LRT2, which can not be identified by Z 2, and only 2 additional loci are identified by Z 2. However, there exist multiple testing problems because we simultaneously analyze 293 loci. So, Bonferroni correction is used (α′ = 0.05/293 = 1.71 × 10−4) and there is no statistically significant result to occur. The more details can be found in Tables K–M in S1 File.

Table 5

LRT0, LRT1, LRT2, Z 0, Z 1 and Z 2 results of application to rheumatoid arthritis data at 5% level.

A. Contingency table showing LRT0b and Z0 results at 5% level.
	PZ0 < 0.05	PZ0 ≥ 0.05	Total
PLRT0b < 0.05	11	6	17
PLRT0b ≥ 0.05	4	272	276
Total	15	278	293
B. Contingency table showing LRT1 and Z1 results at 5% level.
	PZ1 < 0.05	PZ1 ≥ 0.05	Total
PLRT1 < 0.05	9	1	10
PLRT1 ≥ 0.05	4	279	283
Total	13	280	293
C. Contingency table showing LRT2b and Z2 results at 5% level.
	PZ2 < 0.05	PZ2 ≥ 0.05	Total
PLRT2b < 0.05	14	12	26
PLRT2b ≥ 0.05	2	265	267
Total	16	277	293

To investigate the computational efficiency of the XHWE software, we implement the code with the default arguments for this dataset (1217 families and 293 SNPs), on a HP 2311f personal computer (Microsoft Windows 7 Enterprise (Service Pack 1), 4GB of RAM and 3.40 GHz Intel(R) Core(TM) i7 Duo processor) and record its computational time. This process needs 977 seconds. Therefore, on the average, the running time for a single SNP is about 3.3 seconds. For the genome-wide case, for example, one would analyze 200000 SNP markers on X chromosome for the family sample of the type mentioned above, which would lead to 1600000 tests for the hypotheses with running time being about 7.6 days on the personal computer of this type.

Discussion

The existing Z 1 and Z 2 tests were respectively proposed to test for H 01 : p = p and H 02 : ρ = 0. However, we find that there is no simulation study conducted to assess the validity of Z 1 and Z 2 and their performance [14]. Further, there is no existing method to simultaneously test for H 0 : p = p and ρ = 0. Therefore, in this article, we first combine these two test statistics and suggest Z 0 = Z 1 + Z 2 to test for the equality of the frequencies of the same allele in males and females and the zero inbreeding coefficient on X chromosome based on the collected sample, because Z 1 and Z 2 are independent of each other. What’s more, for the purpose of improving the test power, we propose several LRT-type test statistics. Firstly, we write out the likelihood functions under H 0 : p = p and ρ = 0 and H 1 : p ≠ p or ρ > 0 at a single SNP locus on X chromosome, respectively. Then, we obtain the MLEs of the male allele frequency, the female allele frequency and the inbreeding coefficient by the EM algorithms, where we use the RMSE and bias to assess the accuracy of the MLEs of these unknown parameters and construct the corresponding likelihood ratio test (LRT0) statistic under the null hypothesis H 0. If H 0 is statistically rejected, we further develop two LRT-type test statistics LRT1 and LRT2 respectively for H 01 : p = p and H 02 : ρ = 0. Note that LRT0 and LRT2 are too conservative from the simulated results. So, we use parametric bootstrap techniques and propose the LRT0 and LRT2 test statistics. We simulate the data under different parameter settings. Simulation results show that the proposed bootstrap-based methods LRT0 and LRT2, LRT1, Z 0 and the existing Z 1 and Z 2 control the type I error rates well under the respective null hypothesis. Power comparison demonstrates that LRT0 is more powerful than both LRT0 and Z 0. Under ρ > 0, LRT2 has much more power than LRT2 and Z 2, and LRT2 is a little less powerful than Z 2. In addition, LRT1 and Z 1 almost have the same power under p ≠ p .

As for the parameter estimates, the estimate of p based on the EM algorithm is the same as that in zheng et al. [14]. Further, the estimates and of p and the pooled allele frequency p based on the EM algorithms have the RMSE and bias similar to those from zheng et al. [14], respectively. However, although the biases of and based on the EM algorithms are larger than from zheng et al. [14] for some cases, the RMSEs of and are smaller than for all the simulation settings. In addition, the population stratification model indeed causes excess homozygosity, which is consistent with Overall and Nichols [24]. The mean values ( and ) using the EM algorithm are a little larger than proposed in zheng et al. [14], while and have less standard deviation.

Note that ρ = 0 and ρ > 0 in the null and alternative hypotheses of the likelihood ratio test LRT0 or LRT2, respectively, which causes the “boundary” problem and that the corresponding likelihood ratio test is not expected to follow a χ 2 distribution [31, 33]. This may be the reason why the size of LRT0 and LRT2 is too conservative. Therefore, we use parametric bootstrap techniques to obtain the exact distributions of LRT0 and LRT2 in this article.

Due to the presence of the X chromosome inactivation (XCI) and dosage compensation (DC), association analysis and excess homozygosity tests on X chromosome are more complicated than those on autosomes [34]. In the presence of XCI, only one allele from a pair of alleles in females is expressed [35]. Consequently, if considering a locus with two alleles M 1 and M 2, the effect of the M 1 allele in males should be equivalent to the difference between M 2 M 2 and M 1 M 1 homozygous females. As such, when we conduct analyses based on allele-counting, we must either count each allele twice in males or equivalently count each allele in females as 0.5, reflecting a “dosage compensation” for X inactivation [34]. It should be noted that LRT2, LRT2 and Z 2 for H 02 : ρ = 0 are not affected by XCI and DC because they only use female individuals in the collected sample. Similarly, Z 1 for H 01 : p = p is also not influenced by XCI and DC because it estimates the allele frequencies and the corresponding variances in males and females, respectively. Thus, Z 0 = Z 1 + Z 2 is still valid when XCI and DC exist. To investigate the effect of XCI and DC on LRT0, LRT0, LRT1 and LRT1, where LRT1 is the bootstrap version of LRT1, we carry out simulation study under several simulation settings in the presence of XCI and DC. The simulation settings and simulation results are listed in Table 6. It is shown in the table that the size of LRT2, Z 0, Z 1 and Z 2 stays close to the nominal 5% level and the size of LRT2 is still conservative. However, LRT0 and LRT1 without bootstrap cannot control the size well. Fortunately, the type I error rates of LRT0 and LRT1 with bootstrap are very close to 5%. Furthermore, LRT0 is more powerful than Z 0 almost for all the cases and LRT1 and Z 1 almost have the same performance in power. Therefore, in the presence of XCI and DC, LRT0, Z 1 and LRT2 are recommended. Finally, LRT0 and LRT2 can deal with samples of small size. However, LRT0 and LRT2 are based on the parametric bootstrap techniques, which are more computationally intensive.

Table 6

Simulated size/powers (in %) of LRT0, LRT0, LRT1, LRT1, LRT2, LRT2, Z 0, Z 1 and Z 2 based on 10000 Monte Carlo replications and 1000 bootstrap replications under X chromosome inactivation and dose compensation, having p = 0.3 and the ratio N : N = 1: 1.

N	ρ	p f	LRT0	LRT0b	LRT1	LRT1b	LRT2	LRT2b	Z 0	Z 1	Z 2
800	0.00	0.30	6.5	5.0	10.6	4.9	2.3	5.2	4.8	4.8	4.9
	0.05	0.30	17.2	13.3	10.8	4.9	16.2	25.9	13.0	5.0	17.0
	0.10	0.30	44.1	38.5	11.4	5.2	49.3	63.0	41.2	5.3	50.5
	0.00	0.34	31.7	27.3	42.2	29.2	2.5	5.4	23.5	29.7	5.2
	0.05	0.34	41.8	36.5	41.2	27.9	15.8	26.3	31.6	28.4	16.4
	0.10	0.34	65.1	59.5	41.0	28.2	50.9	64.0	58.5	28.5	52.0
	0.00	0.35	43.6	38.7	55.5	41.3	2.2	4.9	33.6	41.9	5.0
	0.05	0.35	54.3	48.9	55.3	41.5	15.8	26.0	42.0	41.9	16.7
	0.10	0.35	72.1	67.3	53.7	40.2	49.9	63.1	64.8	40.6	50.9
1200	0.00	0.30	6.7	4.8	10.9	5.1	2.2	4.8	5.1	5.0	5.5
	0.05	0.30	22.4	18.0	10.8	5.1	23.0	34.3	18.4	5.0	24.3
	0.10	0.30	60.9	54.8	11.2	5.5	67.3	78.5	58.6	5.5	68.3
	0.00	0.34	42.7	37.8	55.0	40.7	2.2	4.9	32.5	41.3	4.9
	0.05	0.34	57.3	52.1	54.4	40.8	22.1	33.6	46.6	41.0	22.9
	0.10	0.34	81.0	76.8	53.0	39.3	67.5	79.3	76.0	39.5	68.7
	0.00	0.35	59.5	53.8	70.6	57.4	2.3	4.9	47.5	57.9	5.0
	0.05	0.35	71.2	66.2	70.3	57.5	22.7	34.3	59.7	57.8	23.8
	0.10	0.35	88.1	84.7	68.5	55.4	67.7	78.8	83.6	56.0	68.7

Supporting Information.

Tables A–J, root mean squared errors (RMSE) and biases of estimates of p , p and ρ based on EM algorithm and zheng et al. [14] under different simulation settings. Tables K–M, LRT0, LRT0, Z 0, LRT1, Z 1, LRT2, LRT2, and Z 2 results of application to rheumatoid arthritis data, respectively. Figs A–L, simulated size/powers of LRT0, LRT0, LRT1, LRT2, LRT2, Z 0, Z 1 and Z 2 against r = N : N based on 10000 replicates under different simulation settings.

(PDF)

Click here for additional data file.