Large-scale multiple testing in genome-wide association studies via region-specific hidden Markov models

Candidate change point set

To effectively reduce the huge space of competing change points and save computation time, our ACP method needs a candidate change point set in advance. Here, we use a haplotype block partition method [23] to obtain the haplotype-block boundary points for each chromosome, which can be collected as the candidate change point set. Because the minimum length value of block L _min should be pre-specified in their haplotype block partition method, here, we let L _min be 300 for all our analysis.

Adaptive criterion-based partitioning procedure

Simultaneously inferring m _c and w _c can be regarded as a model selection problem. To select a desired model, the commonly used methods are established base on the criterion of minimizing the penalized negative maximum likelihood (e.g. BIC). However, many other existing criterions including BIC, assume that the observations are independent, which is not true in HMM. As a result, the effective sample sizes may be small owing to strong dependence among the observations, and the existing criterions may suffer from a failure of consistency. A data-driven penalized criterion was proposed in the Gaussian and least-squares regression model selection for independent observations [17, 24, 25]. Especially, [25] used this adaptive criterion for variable selection and clustering in Gaussian mixtures model and showed that this adaptive criterion outperforms other criterions (e.g. BIC) for small sample sizes. Following their work, we propose a data-driven penalized criterion for dependent observations in HMM.

where ${\widehat{\Psi }}_c$ is the maximum likelihood estimator of the parameters Ψ _c in HMMs for chromosome c using an EM algorithm, and the penalty function ${\text{pen}}_{L_c}(w_c,m_c)={\lambda }_c{D}_{(m_c,w_c)}$ is designed to avoid overfit problems. In this penalty function, λ _c >0 is a tuning parameter to be chosen depending on sample size $L_c=\sum _{r=1}^{R_c}{L}_{\mathit{\text{cr}}}$ and $D_{(m_c,w_c)}$ is the number of parameters in the model. Furthermore, in this paper, we have $D_{(m_c,w_c)}=\left(\right|w_c|+1)D_{m_c}$, where $D_{m_c}$ is the number of parameters that only depend on m _c . If we let ${\lambda }_c=\frac{ln\left(L_c\right)}{2}$, this penalty function becomes the penalty function of BIC for HMM.

Given a value of λ _c in the penalized criterion ${\text{pen}}_{L_c}(w_c,m_c)={\lambda }_c{D}_{(m_c,w_c)}$, we can find ${ŵ}_c$ and ${\widehat{m}}_c$ to minimize ${\text{crit}}_{{\lambda }_c}(m_c,w_c)$ of equation 1 by running Algorithm 1.

Algorithm 1

In step 1, we need to give pre-specified values of K _max and m _max in advance, where K _max denotes the maximum value of the number of change points for each chromosome, and m _max is the maximum value of the number of component. As we know, the number of true change points is usually far less than the number of the candidate change points in practical applications, so we can give a smaller value for K _max to save computation time. For m _max , Wei et al. [15] suggested that values of between four and six are usually chosen.

In step 2, following the methods of [26] and [27], we provide an optimal partitioning search method for change points to estimate ${ŵ}_{c,i}$ given λ _c and m _c , which is, in essence, a dynamic program (DP) algorithm. The detailed procedure about the optimal partitioning search method is shown in Additional file 1.

However, in practice, λ _c is unknown and needs to be calibrated and estimated from the data themselves. Slope heuristics [24] as well as its generalization, dimension jump method [17], are practical and effective calibration algorithms to estimate the optimal penalty ${\text{pen}}_{L_c,\mathit{\text{opt}}}(w_c,m_c)={\lambda }_{c,\mathit{\text{opt}}}{D}_{(m_c,w_c)}$. Here, we propose a sliding window-based dimension jump method to estimate λ_{c,o p t}, where the sliding window is used to avoid losing cases involving several successive jumps. When the width of the sliding window is 1, our proposed method becomes the dimension jump method of Wei et al. [17]. The following algorithm describes the detailed procedure for estimating λ_{c,o p t}.

Algorithm 2

At the end of Algorithm 2, we can obtain the estimation ${\widehat{\lambda }}_{c,\mathit{\text{opt}}}$ of the λ_{c,o p t}. Having ${\widehat{\lambda }}_{c,\mathit{\text{opt}}}$, we can then run Algorithm 1 to obtain ${\widehat{m}}_{c,\mathit{\text{opt}}}$ and ${ŵ}_{c,\mathit{\text{opt}}}$ as well as the desired optimal model by minimizing ${\text{crit}}_{{\lambda }_c}(m_c,w_c)$ of Equation (1). At the same time, we can get ${\widehat{\Psi }}_c$, the estimates of model parameters Ψ _c based on the optimal model, where c=1,…,C.

Theorem 1

Consider the multi-region HMMs defined in section 'Region-specific multi-HMM with change points and where the number of components known’. Let LIS₍₁₎,…,LIS_(L) be the ranked LIS values from all the regions of all chromosomes. Then, the RSPLIS procedure controls the global FDR at level α. In addition, the global FNR level of RSPLIS is β^∗+o(1), where β^∗ is the smallest FNR level among all valid FDR procedures at level α.

Simulations of the ACP method performance for model selection

which ensures that the true change point set {iL₀;i=1,…,4}⊂w⁰. To compare the performance of the two methods, we used sensitivity and specificity as measures. Sensitivity is defined as the average proportions of the true change points which are correctly identified as change points over the 100 times and the specificity is defined as the average proportions of the false change points which are not identified as the change points over the 100 times. We set K _max =8 and m _max =5 for our ACP method. At the same time, h takes the values of 2, and 20 for the window.

HMM-based simulations for comparing RSPLIS with PLIS in GWAS

For simplicity, in this simulation, suppose there are two chromosomes (c=2) in total, each of which consists of two stationary regions, and each region has 2000 SNPs (L _cr =2000,r=1,2,c=1,2). For each chromosome, we set K _max =4, m _max =3, and h=2 for our ACP method and gave the candidate change point set $w_c^0=\{300i;i=1,2,\dots ,\frac{L_c-300}{300}\}$, c=1,2. The purpose of this simulation is to compare RSPLIS with PLIS by finding disease-associated SNPs while controlling the FDR at a pre-specified level α=0.1 for the two chromosomes (combining chromosomes 1 and 2). We conducted simulation studies in the following two cases.

Case 2

In contrast to Case 1, to assess the performance of RSPLIS at the different disease risk levels, we varied the parameters of the z-value distribution while fixing the other parameters. We used the parameter settings in Table 5, where we varied the parameters of the z-value distribution by changing the value of υ₂ (υ₂=0.5,1,1.5,2). Furthermore, we let ${\mu }_{21}^{\left(c\right)}={\mu }_{11}^{\left(c\right)}+{\upsilon }_2$, ${\mu }_{11}^{\left(2\right)}={\mu }_{11}^{\left(1\right)}+0.5$ and varied the disease risk parameter ${\mu }_{11}^{\left(1\right)}$ from 0.5 to 2.0 with an increment of 0.5.

Chromosome	Region	Transition matrix	z-value distribution	z-value distribution
(c)	(r)		(Null)	(Non-null)
1	1	$\left(\begin{array}{ll}0.98& 0.02\\ 0.025& 0.975\end{array}\right)$	N(0,1)	$N({\mu }_{11}^{\left(1\right)},1)$
	2	$\left(\begin{array}{cc}0.98& 0.02\\ 0.025& 0.975\end{array}\right)$	N(0,1)	$N({\mu }_{21}^{\left(1\right)}={\mu }_{11}^{\left(1\right)}+{\upsilon }_2,1)$
2	1	$\left(\begin{array}{cc}0.98& 0.02\\ 0.025& 0.975\end{array}\right)$	N(0,1)	$N({\mu }_{11}^{\left(2\right)}={\mu }_{11}^{\left(1\right)}+0.5,1)$
	2	$\left(\begin{array}{cc}0.98& 0.02\\ 0.025& 0.975\end{array}\right)$	N(0,1)	$N({\mu }_{21}^{\left(2\right)}={\mu }_{11}^{\left(2\right)}+{\upsilon }_2,1)$

The simulation results are shown in Figures 1, 2, 3 and 4. From Figure 1 and Figure 3, we can obviously see that both RSPLIS and PLIS are well controlled at FDR level 0.1 asymptotically. Figure 2 and Figure 4 inform us that PLIS is dominated by RSPLIS for the power at Case 1 and Case 2, which indicates that our RSPLIS procedure is effective at dividing the chromosomes into smaller and more homogeneous regions by searching the change points. In addition, the difference in FNR levels (RSPLIS vs. PLIS) becomes smaller as ${\mu }_{11}^{\left(1\right)}$ increases for each model, which implies that RSPLIS is especially useful when the disease signals are weak.

To show that the higher power of RSPLIS is not gained to the detriment of a higher FDR level, we conducted a further simulation study. This study evaluated the sensitivities at different FDR levels for υ₁=0,0.15,0.30,0.45 in Case 1, and υ₂=0.5,1.0,1.5,2.0 in Case 2, where the sensitivities were calculated as the average proportions of correctly identified SNPs over the 100 replications. For the purpose of illustration, we have only listed the results for υ₁=0.15 of Case 1 and υ₂=1.0 of Case 2 in Figure 5 and Figure 6 respectively, because the other results were broadly similar. It is clear from Figures 5 and 6 that RSPLIS discovers more true disease-associated SNPs than PLIS at the same FDR level.

Genotype-based simulations for comparing RSPLIS with PLIS

This simulation evaluated the performance of selecting the relevant SNPs for RSPLIS and PLIS based on the genotype data set. In contrast to the simulation study in Subsection 'Application to the Daly data set’, we generated case-control genotype data with more realistic LD patterns. To this end, we constructed a genotype pool composed of genotypes from 60 samples for 23 chromosomes by randomly matching the pair of the known phased 120 haplotypes from the Illumina 550K. For simplicity, we only used SNPs selected from 2001-th SNP to 6000-th SNP of chromosome 7 and chromosome 15, respectively. Four SNPs were selected from each chromosome as the disease causal SNPs, each with a relative risk of 1.5. Specifically, the four SNPs, 400-th, 900-th, 1750-th, 3200-th, were chosen to be far away on chromosome 7, while the four other SNPs, 5600-th, 5604-th, 5608-th, 5612-th, were chosen based on their proximity to each other (i.e., separated by 3 SNPs) on chromosome 15.

where β₀=-9.5425 for a disease rate of 0.03, β _i = log(1.5), for i=1,…,8. We repeated the sampling procedure until we obtained 1000 cases and 1000 controls are obtained. The eight disease casual SNPs were then removed from the simulated data set, and the 39 SNPs that contained the three adjacent SNPs on each side of the eight disease-causal SNPs were regarded as the relevant SNPs. We assessed the performance of the testing procedure by the selection rate of relevant SNPs, where the percentages of the true positives (sensitivity) selected by the top M SNPs could be calculated easily. We set m _max =6, h=2, and K _max =5 for the ACP method.

We have plotted the average sensitivity curves for comparisons of RSPLIS vs. PLIS in Figure 7. It is apparent that our RSPLIS dominates PLIS in ranking the revelant SNPs. In summary, these results show that exploiting the heterogeneous chromosomal regions and searching the change points to find chromosomal regions that are more homogeneous has improved the precision of RSPLIS in that the number of false positives has been reduced while the statistical power has increased.

Model selection and estimation of HMM parameters

First, we used the K-Nearest neighbor method proposed by the R package [29] to impute the missing genotypes from the Daly et al. data [18]. Then, for each SNP, we obtained a p-value by conducting a χ²-test to assess associations between the allele frequencies and the disease status, furthermore, we get z-value by transforming p-value. We assumed that the null distribution is standard normal N(0,1) and the non-null distribution is a normal mixture$\sum _{i=1}^{m_c}{\xi }_{i}N({\mu }_{\mathit{\text{ri}}}^{\left(c\right)},{{\sigma }_{\mathit{\text{ri}}}^{\left(c\right)}}^2)$, where c=1 because there is only one chromosome in the Daly et al. data [18]. We used our ACP method to select the number of components and change points, where the parameters in HMMs were estimated by the EM algorithm. Thereafter, RSPLIS was used for multiple testing.

Data analysis