A Graphical Weighted Power Improving Multiplicity Correction Approach for SNP Selections

Controlling for the multiplicity effect is an essential part of determining statistical significance in large-scale single-locus association genome scans on Single Nucleotide Polymorphisms (SNPs). Bonferroni adjustment is a commonly used approach due to its simplicity, but is conservative and has low power for large-scale tests. The permutation test, which is a powerful and popular tool, is computationally expensive and may mislead in the presence of family structure. We propose a computationally efficient and powerful multiple testing correction approach for Linkage Disequilibrium (LD) based Quantitative Trait Loci (QTL) mapping on the basis of graphical weighted-Bonferroni methods. The proposed multiplicity adjustment method synthesizes weighted Bonferroni-based closed testing procedures into a powerful and versatile graphical approach. By tailoring different priorities for the two hypothesis tests involved in LD based QTL mapping, we are able to increase power and maintain computational efficiency and conceptual simplicity. The proposed approach enables strong control of the familywise error rate (FWER). The performance of the proposed approach as compared to the standard Bonferroni correction is illustrated by simulation and real data. We observe a consistent and moderate increase in power under all simulated circumstances, among different sample sizes, heritabilities, and number of SNPs. We also applied the proposed method to a real outbred mouse HDL cholesterol QTL mapping project where we detected the significant QTLs that were highlighted in the literature, while still ensuring strong control of the FWER.


INTRODUCTION
Linkage Disequilibrium (LD) analysis plays a fundamental role in gene mapping, as a tool for uncovering biological trait regulating genes. Many biological traits are influenced by genetic variants and hence it is possible to determine the rough genomic position of the causative variations through associations between SNPs and phenotype [1][2][3][4][5][6][7][8][9][10][11][12][13]. Among the popular SNP selection approaches, the single-locus association with solid multiplicity correction ability remains a powerful tool as associations can generally only be found over small distances [14]. Moreover, as tensof-thousands of SNPs for genome-wide association studies (GWAS) are under demand [5], the single SNP based LD analysis can at least provide a necessary initial screening selection to detect a subset of promising candidates for further exploration [15,16].
Despite the great progress which has already been made within LD based Quantitative Trait Loci (QTL) mapping, powerful and computationally efficient methods for largescale simultaneous testing of individual SNPs with strong control of the familywise error rate (FWER) are still lacking [17,18]. FWER is the most accepted error rate used to determine significance levels for large-scale testing problems *Address correspondence to this author at the Department of Mathematics and Statistics, Utah State University, Logan, UT 84322, USA; Tel: +001 435 797 0749; E-mail: Guifang.fu@usu.edu.
where the goal is to provide conclusive results. In some studies, researchers often control the False Discovery Rate (FDR) to obtain a large pool of potentially significant SNPs and then select only the most significant subset for validation due to cost restrictions. However , this rule can lead to unwanted results as the FDR is controlled only for all selected SNPs , and provides no promise of control for an arbitrarily selected subset of the significant SNPs. Thus, we recommend controlling the FWER (in place of the FDR) in exploratory scenarios where only the most promising results will be considered.
The Bonferroni correction, as one of the most widely used statistical procedures, is often employed to control the FWER when multiple tests are conducted. However, the Bonferroni correction is not favorable in large-scale testing because it substantially reduces the statistical power, hence decreasing the chances of detecting SNPs with real effects [19]. While permutation procedure has been widely employed to adjust for correlated tests, it is computationally expensive [15,20,21] and may mislead in the presence of family structure [22]. Moreover, the permutation approach was designed for only one test setting. In LD based QTL studies, the high likelihood of dependencies among SNPs and the two tests structure strongly demand a new multiplicity adjustment approach that maintains simplicity but is more powerful. hypothesis tests, one testing for the existence of a QTL for a given SNP (i.e. whether or not the QTL is associated with the phenotype), and the other testing for the strength of the LD between the SNP and the existing QTL (i.e. whether or not the QTL is successfully detected by the model). Although the existence of a significant QTL is the ultimate goal, the degree of LD between the QTL and SNP is also critical in guaranteeing the basic assumptions of the model. By Fu et al.'s assumptions [13], the unobservable QTL can be mapped by its association with the observable SNP through the conditional probability of the genotype of the QTL given the genotype of the SNP. Therefore, only QTLs that are not only significantly existing but also strongly linked with a SNP will be identified, i.e. identifying a significant QTL requires rejecting these two hypotheses simultaneously. Although the LD based QTL model has been successful in locating significant QTLs [13,23,24], the Bonferroni multiplicity correction approach used previously ignored two important issues. First, if the QTL existence test fails to reject, then the LD between the SNP and QTL is not identifiable in their mixture model. Second, the Bonferroni correction is too conservative for large-scale of SNPs.
In this article, we propose a new power improving multiplicity correction approach specially designed for the LD based QTL mapping, on the basis of graphical weighted-Bonferroni methods [25]. By introducing a logical structuring for the two tests involved for each SNP, i.e. higher order for QTL existence testing (primary) than the LD testing (secondary), the LD test will never be investigated if the primary test concludes that the QTL does not exist. By exploiting the priority ordering of the two hypotheses to adjust the p-values, our proposed approach can avoid the previously mentioned identifiability issue, and address the multiplicity correction for large-scale number of SNPs. None of the current LD based QTL methodologies directly overcome these challenges when performing these two tests [13,23,24]. Our proposed multiple correction approach with priority structuring has been shown to synthesize weighted Bonferroni-based closed testing procedures such as the weighted Bonferroni-Holm procedure, fixed sequence tests, gatekeeping procedures, and the fallback procedure into a powerful and versatile graphical approach [25], which we tailor here for the LD based QTL mapping.
In the following section we present the LD based QTL model and the two tests involved. Next, we describe in detail how we design the logical structuring to perform the multiplicity correction for the LD based QTL model. The significance of the power advantage of the proposed method over the Bonferroni correction is established through both simulations and one real QTL mapping project. Since sample size, heritability, and number of SNPs all determine the power of the method, we illustrate the power through heritability of 0.1 and 0.4, sample size small (100), medium (300), and large (500), and number of SNPs changing from 1, 10, 50, 100, 500, to 1,000. We end with a discussion of the results.

LD Based QTL Mapping Model
To map the rough location of the QTL regulating a certain biological trait, we apply the mixture model [13]. Under this model, QTL is detected by statistically modeling the genotypic variation through not only the association between phenotype and the putative QTL, but also the association between the putative QTL and SNP. Since the SNP genotype is observable, we can infer the probabilities of a putative QTL genotype by the conditional probability of QTL ( A) genotype given the SNP ( M) genotype, as long as there exists LD between the SNP and putative QTL [26].
The mixture model of [13] assumes each individual's phenotype i Y , n i , 1, = … , is a random variable from density denotes three distinct QTL genotypes. Each QTL genotype is assumed to induce a separate distribution of phenotypes. Typically, normal distributions are assumed for each . From these assumptions, the corresponding likelihood is expressed as [13] ), , where i l| is the conditional probability of individual i having QTL genotype l given their SNP genotypes , l is the phenotypic mean for QTL genotype l , is the common standard deviation for all genotypes, and ) , is the probability density of observations for individual i at QTL genotype l [13,26,28].  [13,27]. Hence, i l| is a function of p , q , and D . The EM algorithm is then applied to the likelihood in (1) to obtain maximum likelihood estimates for all parameters [13,27].

Two Hypothesis Tests
Through the likelihood in (1), the hypotheses vs = = : can be used to test if the QTL is significantly associated with phenotype Y ( i.e. existence of QTL ). Since all the unknown parameters in (1) were estimated by maximum likelihood estimates (MLEs), a log likelihood ratio statistic can be used to test the hypotheses in (2) [13]. The resulting test statistic ( 2 L ) is asymptotically distributed as a 2 2 under L H 0 for large enough samples [13,26].
Once the existence of a QTL is established, the test statistic used to judge whether or not the QTL is significantly associated with SNP is [13,29]: Here, 2 r is the square of the correlation coefficient between the SNP and QTL that has been used in most of the related literature, which has many good sampling properties [30,31]. Under D H 0 , 2 D is asymptotically distributed as 2 1 , from which the tail probability (p-value) of the observed level of association can be determined [10,23,26,29,32].
Because whether or not a QTL exists and whether or not the existing QTL is able to be successfully identified by the model are both critical, the QTL will not be interpreted as significant unless the two hypothesis tests (2) and (3) [33,34].
That is, the parameter i l| falls out of the model when the means are equal, as the , so that D which is computed from i l| can not be computed. As a result, testing 0 D under L H 0 is not meaningful. Secondly, a multiplicity correction is needed for simultaneous testing of all SNPs. Thirdly, the existence of a QTL underlying certain biological traits is the ultimate goal for the real application. Inspired by the idea of the graphical Bonferroni approach [25], we set the existence of QTL (2) to be the primary test and the LD test (3) to be the secondary test. If the primary test is not rejected, the secondary test will not be investigated. As a result, our proposed multiplicity correction approach increases the power, while preserving strong control of FWER and avoiding the identifiability issue of D under L H 0 .

Graphical Bonferroni Approach
The Graphical Bonferroni Approach (GBA) is a versatile and easily communicated general adjustment method for multiple testing [25]. Provided as a generalized framework, it must be specially tailored to each testing situation. Generally speaking, it is most powerful for situations where hypotheses can be partitioned into levels of importance such that the most important hypotheses are tested first and the lower level hypotheses are tested only if the higher level hypotheses show significant results.
All hypotheses of interest are depicted as nodes in a directed acyclic graph. Local significance thresholds for each node (hypothesis) dictate the local level at which each hypothesis is tested. Weighted edges between all nodes map the logical structuring of the designated testing approach. When a hypothesis is rejected, the weighted edges dictate the proportion of the locally assigned significance threshold that is passed from the rejected node to all connected nodes. Thus, the graph induces an iterative testing approach that is shown to result in a closed-test that admits a short-cut [25]. Further, Algorithm 1 of [25] provides a simple updating technique that performs the short-cut. Strong control of the FWER at level is proven to occur so long as three regularity conditions are met: 1) the sum of the local significance thresholds is no more than , 2) the sum of outgoing edge weights from each node are no larger than unity, and 3) no node has an edge connecting to itself [25].

Rejection Scheme
Since the ultimate goal is to check existence of QTL, the first interest is in testing (2) to see if the phenotype shows evidence of association with a latent QTL. Depending on the results of the test of (2), the testing for the given SNP will either end, or interest will be turned to testing   Fig. 1) cannot have a smaller adjusted p -value than its parent node ( L H 0 in Fig. 1).
As a result, for the single SNP analysis, either both hypotheses will be tested at level , or the testing will stop after L H 0 without considering D H 0 . Alternatively, the Bonferroni correction would test both hypotheses at /2 . Hence, the Bonferroni adjustment has less power, due to its smaller thresholds. Compared to the traditional Bonferroni, the only potential disadvantage of the GBA method is that it skips testing   Fig. 1, Fig. 2  Li H 0 and Di H 0 ) are rejected for any given SNP i . Second, the -level is split with a Bonferroni type allocation between the m top-level hypotheses while none of is initially provided to the m lower-level hypotheses. Upon rejection of a higher-level hypothesis, the lower-level child hypothesis receives all of the m / -level of the parent (edge weight of 1). If the lowerlevel child hypothesis is then also found to be significant, its threshold is then shared between all remaining higherlevel hypotheses (edge weights of 1/2).
The power advantage of our proposed GBA over the Bonferroni method is evident from the larger thresholds.
Where the Bonferroni method would test each hypothesis at the /6 -level, the GBA tests each hypothesis by thresholds that are no smaller than /3 . To demonstrate, assume that  Fig. 3. Notice in Fig. 3 the reconnecting of edge weights which previously attached to H . This demonstrates how edges determine not only the weight that will be passed, but also define the inheritance of edge weights. Fig. (2). Demonstration of the hierarchy of the GBA testing scheme for three SNPs. Fig. (3). Demonstration of the GBA testing scheme for three SNPs assuming that hypotheses and from the initial graph in (Fig. 2) are rejected.
Assume now that 1 0 D H of Fig. 3 can be rejected at the /3 -level. The graph updating (Fig. 4A) becomes more complicated with this rejection because the rejected hypothesis is both sending out and taking in edge weight from the same hypotheses (nodes). Specifically,   . (4). A: The updated graph from (Fig. 3) assuming the hypothesis of ( Fig. 3) is rejected at the -level. B: Graph resulting from the rejection of the hypothesis at the level.
In conclusion, we solve the three problems in the multiplicity testing scenario that exist in previous LD based QTL models. Firstly, the unidentifiable issue of D under the null hypothesis 3  concern that existence of QTL is the ultimate goal. In the remainder of this article, we show through simulation studies that the proposed GBA is more powerful for LD based QTL mapping than standard Bonferroni adjustments and thus leads to more scientific discovery while maintaining strong control of the FWER.

Power Simulation
We investigated a simulation study to quantify the power advantage of the proposed graphical Bonferroni approach (GBA) over the standard Bonferroni adjustment within the LD based QTL mapping model [13]. The QTL, phenotype, and SNPs were generated under the assumptions of the alternative hypotheses (described in Methodology section in (2) and (3) ). The QTL was generated using an assigned probability of 0.7 = q for the major allele. For each was used to code the QTL genotypes of aa, Aa, and AA, respectively. The normally distributed phenotype dependent on the value of the QTL is generated as . SNPs were then generated using the conditional probability of the SNP genotype given the value of the QTL genotype for each individual. In general, for an LD based QTL mapping model, researchers genotype the SNP first and then use the SNP to generate a QTL based on the conditional probability of QTL genotype given SNP genotype. However, for our purposes, we are interested in extending from single SNP mapping to multiple SNPs mapping. Therefore, we derive the conditional probability of SNP genotype given QTL genotype (see Table 1) from the Bayes Rule in Equation (6).
Sample sizes of 100 = n , 300 , and 500 were used to represent small, medium, and large sample sizes, respectively. The number of SNPs per simulation was set at 1 = m , 10, 50, 100, 500, and 1,000 to show the initial power under the single SNP scenario and the corresponding decreasing power trend as the number of SNPs increases. Finally, the heritability was set at two values, 0.1 = 2 H and 0.4, corresponding to high and low error variance [27]. The model error variance 2 was computed using the heritability and genetic variance of the QTL. Power estimates were averaged over 1,000 simulations.    Table 2 and depicted in Fig. 5 ) and a larger sample size ( 500 = n ) , the power of the multiplicity adjustment remains high even as the number of SNPs becomes large ( 1,000 = m ). However, in practice it is often expensive to collect so many sample measurements. It is worth mentioning that the power obtained from the GBA can achieve 80% for a large number of SNPs ( 1,000 = m ) but medium sample size ( 300 = n ). Moreover, even with a low heritability ( 0.1 = 2 H ), the power increase of the GBA over the Bonferroni adjustment allows for the possibility of maintaining the power level of the Bonferroni adjustment while decreasing the sample size of the study or increasing the number of SNPs, a great advantage for researchers. For example, under many cases, the power of our improved approach for 1,000 SNPs is similar or larger than the power of the Bonferroni adjustment for 500 SNPs.
Although the power increase of our proposed method improves moderately over the standard Bonferroni adjustment for the case of low heritability ( 0.1 = 2 H ) when the sample size is small ( 100 = n ), the power gains are still comparable to seminal results found by previous multiplicity improvements over their competitors [35]. All in all, our proposed GBA method shows a substantial increase in power over the Bonferroni adjustment under all 12 circumstances with the different combinations of sample size, number of SNPs, and heritability. We did not compare with permutation because it is computationally infeasible to compute these six settings using permutation. In addition, permutation approach is designed for the one testing structure, which is not the case here. Finally, we firmly believe that the GBA works better in the LD based QTL model because of the priority order setting, the logical consistency, uniformly better power in theory, and increase in computational efficiency.

Mouse HDL Cholesterol QTL Mapping Project
Epidemiological studies have consistently shown that the level of plasma high density lipoprotein (HDL) cholesterol is negatively correlated with the risks of coronary artery disease and gallstones [36,41]. Because of the inverse relationship between HDL and cardiovascular disease, there has been considerable interest in understanding genetic mechanisms contributing to variations in HDL levels. HDL levels vary considerably in different people, which are affected by interactions of multiple genes and environmental factors, and up to 70% of this variation in humans is genetically determined [37,42]. Because of the concordance between human QTLs regulating HDL and corresponding mouse loci and many easily controlled experimental advantages, mouse has become an animal model in HDL research. Numerous findings in HDL QTL associations are obtained from crosses between different inbred mouse strains. By crossing inbred strains that significantly differ in HDL levels and subsequently testing for association between HDL levels and genetic SNPs in the progeny, numerous significant QTLs involved in HDL have been identified in mouse [36,41,43,47].
Compared to the inbred mice strains with coarse mapping resolution, the QTL research on wild-caught and commercial stocks of outbred mice, as resources for genetic fine mapping, is far under developed. Zhang [48]. Three hundred 4-to-6-week-old male NMRI mice were purchased and individually housed with the same diet and environmental conditions. The blood samples of each mouse were measured by submandibular puncture after a 4-hr fast. Then plasma samples were frozen for measurement of HDL cholesterol. There were 10 mice removed because the standard deviation of individual blood pressure is greater than two. Another two mice were also discarded for their 99% identity of SNP genotypes. This caused the final sample size to be 288. A total of 581,672 high density SNP were initially genotyped by the Novartis Genomics Factory using the Mouse Diversity Genotyping Array [49]. In order to guarantee promising data for association mapping studies [50], only polymorphic SNPs with minor allele frequency greater than 2%, Hardy-Weinberg equilibrium 20 < 2 , and missing values less than 40% were retained. Moreover, identical SNPs within a 2Mb interval were collapsed. This left 44,428 unique SNP genotypes for their resulting analysis using three analysis methods, linear trend test, two way ANOVA, and EMMA [51]. From Zhang's work , adjustments for multiplicity at the genome-wide association level were made using a simulation approach [52] as well as the permutation approach [53]. They identified three loci as significant, with two loci on Chromosome1 (Chr1) and a single locus on Chromosome5 (Chr5) (see Fig. 3 of [48]). However, after a closer investigation, Zhang et al. reported that the significant findings in Mb182 of Chr1 is spurious.
We applied the introduced LD based QTL model [13] and the proposed GBA multiplicity correction approach to this outbred mouse HDL cholesterol genome data to compare our findings with the highly validated discoveries in current literature. Recalling the detailed adjustment structure of the GBA, it can be seen that the adjusted p -value obtained from GBA for the test of D H 0 will never be smaller than that of L H 0 . Hence, reporting the significant adjusted p -values for D H 0 is sufficient for demonstrating those SNPs that show strong evidence of linkage to a significantly existed QTL.
( Fig. 6)  supports two dramatically significant findings, on Chr1 at Mb173 and Chr5 at Mb125. These two significant discoveries are the same as the findings in current outbred mouse literature, compare to Fig. 3 of [48], but with an even stronger signal.
In Table 3 we notice that all significant QTLs detected from outbred mouse by our model are confirmed from reported findings obtained from inbred mouse crosses using very different approaches. Two QTLs have been reported coincident with candidate genes. Chr1 locus at Mb173, the highest peak in Fig. 6, is the major determinant of HDL, which has been detected as QTL Hdlq15 in inbred mouse strains multiple times (as referenced in Table 3). Combining mouse crosses with haplotype analysis for the HDL QTL located on Chr 1 locus at Mb173 reduced the list of candidates to a small amount. Numerous mouse crosses have linked HDL to this region, and Apoa2 has been identified as the gene underlying the QTL [38,40,41,43,46]; this gene has been highlighted in Nature Reviews Genetics [54]. Chr5 locus at Mb125, the second highest peak in Fig. 6, is located in the same locus as QTL Hdlq1 found by [46] and [44] (as referenced in Table 3). In addition, they conclude that Scarb1 (a well known gene involved in HDL metabolism) is the causal gene underlying Hdlq1 by haplotype analysis, gene sequencing, expression studies, and a spontaneous mutation [45,47].
There are a total of two significant detections from a large pool of 44,428 candidate SNPs from our model, with the findings confirmed by inbred mouse findings. Zhang et al. worked on exactly the same data, adjusting for multiplicity at the genome-wide association level using a simulation approach [52] as well as the permutation approach [53]. They made the same two significant detections with less significance level than our p-values, and they also find one spurious QTL. Therefore, our proposed approach brings a useful alternative approach for SNP selection literature.

DISCUSSION
Detecting significant genes that cause disease (for example the inverse relation between human cholesterol and cardiovascular disease) or regulate biological traits through LD based QTL mapping has been popular in many disciplines [1][2][3][4][5][6][7][8][9][10][11][12][13]. The new techniques can simultaneously consider tens of thousands of SNPs and hence bring big challenges to multiple testing. In addition, high dimensional biological traits, often reduced to multiple PC components, have been widely used and add yet another demand for a powerful and computationally efficient approach to adjust for multiple tests [13,[55][56][57].
These multiple tests require an adjustment on the resulting p -values in order to preserve control of the familywise error rate (FWER) at a pre-specified level . In some cases, follow up work on the significant findings may justify using the false discovery rate (FDR) as the error rate of interest. Typically however, the significant results are directly reported and therefore the FWER is the more desirable form of error rate to control. The current standard approach in LD based QTL mapping is to apply a Bonferroni adjustment to correct for multiplicity and preserve the FWER. As is well known, the Bonferroni correction is overly conservative for large numbers of tests, but the Fig. (6). The negative log of the GBA-adjusted -values for for each SNP in the mouse HDL cholesterol QTL mapping project. The red reference line corresponds to a 0.05 familywise error rate.  In this article, we tailored a multiple correction approach, based on graphical weighted-Bonferroni methods [25], which allows for the logical order among the two hypotheses to be structured into the multiplicity correction. As in the LD based QTL mapping model of [13], we need to test two hypotheses for each SNP, one with L H 0 about whether or not an association exists between QTL and phenotype, and the other with D H 0 about whether or not LD exists between SNP and QTL. Among these two tests, the QTL existence test has higher priority because the LD test will not be applicable if a QTL does not exist, and the existence of QTL is the ultimate goal in real applications. Although the logical structure of the two tests is known , none of the current LD based QTL literature considers this priority structure when performing these two tests [13,23,24,27]. By the structure of graphical weighted Bonferroni, the quantitative trait loci test and linkage disequilibrium test are integrated into a combined multiple testing correction framework [58]. As a result, GBA approach has more potential applications in QTL studies. For example, in a haplotype study, we can put QTL test, dominant or additive effect test into one multiple testing correction framework. Hence, if QTL test is not significant, we don't have to test dominant or additive effects.
The significance of the power advantage of the proposed method, established through simulations, and finally on real data, is such that we advocate its use whenever multiple tests are needed for the LD based QTL mapping design, where both L H 0 and D H 0 tests are considered.

CONFLICT OF INTEREST
The author(s) confirm that this article content has no conflict of interest.