BACKGROUND
mapping disease susceptibility loci for complex diseases has proven to be challenging.
the etiology of a complex disease is multifactorial, with both genetic and environmental risk factors and interactions among these risk factors.
genome-wide association  studies are an alternative to traditional family-based linkage studies for mapping disease susceptibility loci for complex diseases.
the most common design to date for gwa studies is to genotype a dense panel  of single-nucleotide polymorphisms  for a large collection of unrelated cases and controls.
snps are then tested one marker at a time for association most often using contingency table analysis, i.e., a x2test or fisher's exact test with one or two degrees of freedom, depending on the hypothesis being tested.
the resulting set of p-values is then adjusted for multiple comparisons, usually a bonferroni correction for the number of tests performed.
single marker analysis involves the marginalization of effects over all genetic and environmental backgrounds.
consequently, effect sizes for single markers tend to be small , necessitating sample sizes of thousands to obtain reasonable levels of power at genome-wide significance levels  <cit> .
multi-locus methods may have increased power if they account for the correlation of loci due to linkage disequilibrium.
however, multi-locus methods generally have been observed to have decreased power because they tend to produce test statistics with larger degrees of freedom  <cit> .
therefore, in designing a powerful and efficient multi-locus method for generating testable hypotheses based on genome-wide snp data, we are motivated to reduce the degrees of freedom.
alternatively, consider a biochemical pathway consisting of many genes, which may be located on different chromosomes and hence may segregate independently.
this motivation raises the general question of how to test a set of markers that may be biologically related but not correlated via linkage disequilibrium.
possibilities include multivariate techniques such as the f test or multiple regression.
however, these techniques can also lose power due to large degrees of freedom.
a third issue for association analysis is low power for rare variants that medical sequencing is expected to uncover.
to overcome low power for variants with frequencies < 5%, specialized methods for grouping variants have been described  <cit> , but such methods were designed only for analysis of rare variants for which linkage disequilibrium is expected to be negligible.
furthermore, the two most commonly used estimators of linkage disequilibrium, r <dig> and d', are biased upwards when allele frequencies are low or sample sizes are small  <cit> .
in response to all of these motivations, we propose a flexible multi-locus method based on unions of multiple snps that yields a test that always has only one degree of freedom.
the method implicitly accounts for linkage disequilibrium and is appropriate for simultaneous analysis of common and rare variants.
we demonstrate by simulation that our method has valid control of the false positive error rate while yielding more power than traditional single marker analysis.
we also demonstrate by simulation that our method is sensitive to any source of differential correlation among markers.
we then analyze two publicly available data sets for parkinson disease  for which no significant associations were detected by single marker analysis  <cit> .
using our new method, we found evidence supporting susceptibility to pd at three loci.
methods
study samples
parkinson disease  is the second most common neurodegenerative disorder after alzheimer disease.
for the discovery sample, we obtained publicly available data from the progeni and genepd studies of familial pd  <cit> .
briefly,  <dig>  cases and controls were genotyped using illumina humancnv <dig> version1_c beadchips  and the illumina infinium ii assay protocol.
for the replication sample, we obtained publicly available data from a study of sporadic pd from the snp resource at the national institute of neurological disorders and stroke human genetics resource center dna and cell line repository http://ccr.coriell.org/ninds/ as well as clinical data  <cit> .
the original genotyping was performed in the laboratory of drs.
andrew singleton and john a. hardy , bethesda, md usa.
briefly,  <dig>  snps were genotyped using the illumina infinium i and humanhap <dig> assays on  <dig> cases and  <dig> neurologically normal controls.
the data made available had already been processed for quality control  <cit> .
the framework of logical unions
to illustrate the construction of a union based on classical  logic, consider two events.
let p be the probability of event x and p be the probability of event y. the probability of the union of events x and y, p , is p  = p  + p  - p .
to estimate the probability of the intersection p , we need a measure of correlation.
one commonly used measure of correlation is the pairwise correlation coefficient, given by r=covvarvar=p−ppp)p).
after rearrangement and substitution, p=p+p−pp−rp)p).
in the genetic context, let a represent the major allele and b represent the minor allele at a snp.
with respect to the minor allele, dominant coding of the three genotypes  is given by  and recessive coding is given by .
additive coding of the number of copies of the minor allele, i.e.,, can be achieved by summing over dominant and recessive variables since  +  =  and thus need not be separately considered  <cit> .
r <dig> is a commonly used measure of linkage disequilibrium.
thus, unions implicitly account for correlation between markers, including coupling  or repulsion  phase.
it is critical to note that this use of the correlation coefficient encompasses gametic phase disequilibrium as well as other sources of correlation such as epistasis and natural selection.
also note that, by definition, p  = p, demonstrating the equivalence of analysis of a union of one marker with traditional single marker analysis.
there are many possible coding schemes one could implement within this framework.
in this study, we construct unions over consecutive, non-overlapping sets of k markers.
the choice of non-overlapping sets was made to facilitate calculating the genome-wide testing burden.
let xij be a binary indicator variable for the ith individual and the jth union.
with two markers  under dominant coding for the minor allele,
xij={0for genotypea1a1a2a21for genotypesa1b1a2a <dig> b1b1a2a <dig> a1a1a2b <dig> a1a1b2b <dig> a1b1a2b <dig> a1b1b2b <dig> b1b1a2b <dig> b1b1b2b21for genotypesa1b1?2?
<dig> b1b1?2?
<dig> ?1?1a2b <dig> ?1?1b2b2?for genotypesa1a1?2?
<dig> ?1?1a2a <dig> ?1?1?2?
<dig>
in which ?
represents missing data.
under dominant coding, a value of  <dig> indicates the absence of at least one b allele  and a value of  <dig> indicates the presence of at least one b allele .
similarly, with two markers under recessive coding for the minor allele,
xij={0for genotypesa1a1a2a <dig> a1a1a2b <dig> a1b1a2a <dig> a1b1a2b21for genotypesa1a1b2b <dig> a1b1b2b <dig> b1b1a2a <dig> b1b1a2b <dig> b1b1b2b21for genotypesb1b1?2?
<dig> ?1?1b2b2?for genotypesa1a1?2?
<dig> a1b1?2?
<dig> ?1?1a2a <dig> ?1?1a2b <dig> ?1?1?2?
<dig>
under recessive coding, a value of  <dig> indicates the absence of at least one bb genotype  and a value of  <dig> indicates the presence of at least one bb genotype .
any individual for which xij = ?
is removed from analysis of the jth union.
note that this choice of implementation relies solely on counting individuals with pre-specified multi-locus genotypes, does not require estimation of pi, pj, or r, and does not assume hardy-weinberg equilibrium.
for case-control data, we test for association using fisher's exact test.
for the jth union, the  <dig> ×  <dig> contingency table is constructed by counting the numbers of individuals for which xij is  <dig> or  <dig> among cases and controls.
the test has only one degree of freedom, regardless of k. we control for multiple comparisons using a bonferroni correction for the number of unions tested.
the number of unions tested is given by ⌈#snpsk⌉, in which ┌x┐ is the ceiling function and returns the smallest integer not less than x. this implementation of multiplicity control preserves more power as k increases and as the testing burden is consequently reduced.
given  <dig>  snps in the discovery sample and accounting for both dominant and recessive coding over unions of size one to five, the genome-wide significance level was  <dig> ×∑k=15⌈ <dig> k⌉= <dig> ×10− <dig>  all union testing was performed in r  <cit> .
other analyses
imputation was performed using mach, version  <dig> .
<dig> http://www.sph.umich.edu/csg/abecasis/mach/download/.
for the reference panel, we retrieved the combined hapmap phase ii+iii raw genotype files for the ceu sample http://hapmap.ncbi.nlm.nih.gov/downloads/genotypes/latest_phaseii+iii_ncbi_b36/forward/non-redundant/.
we filtered the  <dig> , <dig> autosomal ceu snps based on the inclusion of unrelated individuals only, a minor allele frequency ≥  <dig> , a snp missingness rate ≤ 5%, and an individual missingness rate ≤ 5%.
we inferred haplotype phases for the reference data using the settings --rounds  <dig> --states  <dig>  we conditioned imputation on the maximum-likelihood estimates of the crossover map and the error rate map.
we retained all imputed genotype calls for which the posterior probability ≥  <dig> .
we performed haplotype analysis using plink version  <dig>  http://pngu.mgh.harvard.edu/purcell/plink/ <cit> .
briefly, we compared each haplotype against all of the other haplotypes.
RESULTS
simulation analysis
we first investigated the validity and power of our proposed method, exploring over sizes of unions of snps ranging from one to five.
note that analysis of a union of one snp is identical to single marker analysis.
under the null hypothesis of no association, analysis of unions has the expected per comparison error rate and is therefore valid .
for small effect sizes , analysis of unions is increasingly more powerful than single marker analysis as the frequency of the minor multi-locus genotype decreases .
if a union consists of both risk-increasing and risk-decreasing predictors, effects will cancel and power will be lost .
also, if a union consists of too many predictors with no effect on the outcome, then power to detect a predictor with an effect within the union will decrease .
analysis of unions of correlated predictors can be more powerful than analysis of unions of independent predictors .
epistasis  can either increase or decrease power, depending on the directions of the effect sizes for epistatic effects vs. marginal effects .
for epistasis to be detectable, minor genotype frequencies must be large so that the joint genotype counts are reasonably large.
the union test is a valid test of association of genotype frequencies if the correlation between markers is matched in cases vs. controls .
the union test is also a powerful test of differential correlation among markers, even in the absence of differences in genotype frequencies .
thus, the union test can simultaneously detect differences in genotype frequencies, epistasis, haplotype structure, natural selection, population structure, and any other process affecting correlation among markers.
comparison to other methods
we compared the union test to a logistic kernel machine-based test, a set-based test for common variant analysis that adaptively estimates the degrees of freedom given the correlational structure of the markers in the set  <cit> .
both tests were valid at the  <dig>  significance level .
the union test was more powerful than the kernel machine-based test for multiple independent markers with nonzero effects .
the union test was also more powerful than the kernel machine-based test when the minor genotype frequency was below ~ <dig>  .
the kernel machine-based test was more powerful than the union test for markers with opposing effects , dilute signal , correlated predictors , and epistasis .
the kernel machine-based test was conservative in the presence of differential correlation between markers under the null hypothesis and had no power to detect differential correlation under the alternative hypothesis .
we also compared the union test to the collapsing method  <cit> , a set-based test for rare variant analysis that collapses the data into a x <dig> test with one degree of freedom.
the collapsing method is based on aggregating genotypes across all loci in the set such that "an individual is coded as  <dig> if a rare allele is present at any of the variant sites and as  <dig> otherwise"  <cit> .
thus, the collapsing method is in fact identical to the union test  by construction and power is equivalent for both tests.
the two main differences between the collapsing method and the union test are that the former assumes that minor allele frequencies are low, i.e., ≤  <dig> , and that all markers are independent, i.e., there is no linkage disequilibrium   <cit> .
the collapsing method is thus a special case of the union test, with the union test applicable across the entire range of frequencies and across the entire range of linkage disequilibrium.
real data analysis
data processing for quality control for both samples is depicted in figure  <dig>  for the discovery sample, we retained  <dig> of the  <dig>  cases,  <dig> of the  <dig>  controls, and  <dig>  of the  <dig>  snps.
for the replication sample, we retained  <dig> of the  <dig> cases,  <dig> of the  <dig> controls, and  <dig>  of the  <dig>  snps.
to investigate the possibility of population stratification in the discovery and replication samples, we estimated the variance inflation factor of the genomic control method  <cit> .
we estimated an inflation factor of  <dig>  for the discovery sample and  <dig>  for the replication sample, indicating a negligible inflation of the false positive error rate due to population stratification .
under dominant coding, we identified no regions correlated with pd case/control status at a genome-wide significance level of  <dig>  × 10- <dig> .
we identified seven independent loci at suggestive levels of association .
of these, we replicated the region on chromosome 4p <dig> including the gene gak.
this region was discovered using unions of one, two, or three markers but was replicated only in the single marker analysis .
using haplotype analysis, we found that the association in both samples was driven by the same low-frequency haplotype  with an odds ratio of  <dig>  .
a in order, the snps are rs <dig> and rs <dig>
to directly compare single marker analysis, haplotype analysis, and our multi-locus method, we examined the first multi-locus union in table  <dig>  which consisted of snps rs <dig> and rs <dig>  using the replication sample .
by  <dig> ×  <dig> contingency table analysis of genotypes, rs <dig> was significantly associated with pd  but rs <dig> was not .
similarly, by  <dig> ×  <dig> contingency table analysis under dominant coding, rs <dig> was significantly associated with pd  but rs <dig> was not .
haplotype analysis revealed significant association with pd for the haplotype consisting of the minor allele at rs <dig> and the major allele at rs <dig>  i.e., haplotype tt vs. haplotypes cc and ct with haplotype tc being unobserved .
under dominant coding,  <dig> ×  <dig> contingency table analysis of the union did not yield significant association, because the union tested multi-locus genotype cctt vs. the other eight genotypes, effectively attenuating the signal of haplotype tt .
in general, for a union over k snps, the number of joint genotypes is 3k and  <dig> ×  <dig> contingency table analysis is more powerful than  <dig> × 3k contingency table analysis.
more importantly, as k increases, some joint genotypes likely will be unobserved, leading to cells in the contingency table having counts of zero.
collapsing a  <dig> × 3k contingency table into a  <dig> ×  <dig> contingency table is one way to address this problem of data sparsity.
in our implementation of dominant coding over multiple loci, the reference multi-locus genotype consists solely of diplotypes of the most common haplotype.
this example  illustrates the trade-off of gaining power by reducing the degrees of freedom vs. losing power by attenuation of signal.
under recessive coding, we identified no regions correlated with pd case/control status at a genome-wide significance level of  <dig>  × 10- <dig> .
at suggestive significance levels, we identified two potentially interesting loci at chromosomes 3p <dig> and 22q <dig> .
at chromosome 3p <dig>  snp rs <dig> was associated with decreased risk in the discovery sample but with increased risk in the replication sample.
we investigated this directional inconsistency locus using haplotype analysis .
we identified a haplotype block of eight snps including rs <dig>  in the discovery sample, the haplotype with the strongest effect carried the minor allele at rs <dig> and was associated with decreased risk.
in the replication sample, the haplotype with the strongest effect carried the major allele and was associated with increased risk.
thus, the effects were directionally consistent at the haplotypic level with respect to rs <dig> but were directionally inconsistent with respect to marginal effects at this snp.
at chromosome 22q <dig>  the marginal effects for snp rs <dig> were directionally inconsistent  and different haplotypes drove association in the two samples .
a in order, the snps are rs <dig>  rs <dig>  rs <dig>  rs <dig>  rs <dig>  rs <dig>  rs <dig>  and rs <dig>
a in order, the snps are rs <dig> and rs <dig>
finally, we examined execution time for the union test.
algorithmically, the union test consists of one hypothesis test per union, regardless of the number of markers in the union.
for a union of k markers, the union test is approximately k- <dig> times faster than the corresponding single marker test, because the computational burden of taking the union is less than the computational burden of performing one fisher's exact test.
discussion
recent attempts at elucidating the underlying genetic structure of common, complex diseases have been hampered by the use of methods that do not account for correlation in the data, are overly stringent, and are underpowered in the presence of allelic or locus heterogeneity and rare variants.
compared to other multi-locus methods designed to address these issues, unions offer many advantages, including accounting for linkage disequilibrium, epistatic interactions within  unions, and haplotype structure, with no requirement for the estimation of haplotypes or phase.
union testing can be viewed as a generalization of grouping schemes for analysis of independent rare variants  <cit>  that allows for arbitrary patterns of correlation among common and/or rare variants.
an additional advantage is that the counting algorithm implicitly accounts for correlation without explicit and potentially biased estimation of linkage disequilibrium measures and without the need of permutation testing to establish correct type i error rates.
union tests can be more powerful than single marker tests because they make fuller use of the information  in a multivariable framework.
power may also be gained by collapsing multi-locus genotype information into a single degree of freedom test  <cit> .
the union test is a composite test of association of genotype frequencies and differential correlation among markers.
union tests are also haplotype tests and are sensitive to processes such as natural selection and population structure that may lead to differential correlation among markers.
previous analyses did not detect substantial population structure in either pd data set  <cit> .
however, given that our multi-locus method is more powerful than single marker analysis under many genotypic configurations, it is possible that our method is more sensitive to cryptic relatedness or residual population stratification than is single marker analysis.
to better address this issue, regression with a covariate for parental ancestry provides a framework compatible with unions  <cit> .
our proposed method is very flexible.
the choice of which snps to combine into unions belongs to the investigator.
we chose to combine consecutive snps in non-overlapping sets in order to take advantage of linkage disequilibrium, so that union testing may be viewed as a generalization of haplotype testing, while minimizing the correlation between tests.
with an appropriate grouping scheme, unions may be used to address allelic heterogeneity, which occurs when the phenotype of interest is caused by different mutations within the same gene  <cit> .
the role of allelic heterogeneity in complex disease is exemplified by the several, individually rare mutations in nod <dig> associated with crohn's disease  <cit> .
alternatively, one could focus on rare variants and group snps with minor allele frequencies < 5%.
rare variants have been predicted to play an important role in complex diseases  <cit> , as has been illustrated for colorectal cancer  <cit>  and cholesterol metabolism  <cit> .
unions may also be used to address genetic heterogeneity, which occurs when multiple genes cause the same disease.
although we tested unions using contingency table analysis, it is straightforward to adapt our coding for use with generalized linear models and thus to incorporate covariates into the analysis.
our method is also flexible with respect to genetic models, such that coding for snps may be any combination of dominant, additive, or recessive  <cit> .
depending on the implementation , union testing can also reduce the testing burden compared to single marker analysis.
however, instead of non-overlapping windows as we used herein, a sliding window approach could be employed, at the cost of a higher testing burden.
there are two main disadvantages of unions.
one disadvantage is that the snp and/or haplotype within a union driving association cannot be inferred.
if genome-wide analysis is intended as exploratory and hypothesis-generating, then this disadvantage is not a concern.
the second disadvantage is that there are scenarios under which statistical power may be lost, depending on how collapsing is achieved.
as with other collapsing methods, the power of union testing decreases if the size of the union is too large, whether or not the included variants are causal  <cit> .
this limitation is more pronounced for markers with large minor allele frequencies.
similarly, the power of union testing decreases if both risk-increasing and risk-decreasing effects are included within a union.
one approach to address both attenuation and cancellation of signal is to use a mixture model with one component for risk-increasing effects, one component for risk-decreasing effects, and one component for null effects.
we illustrated union testing by reanalyzing genome-wide genotype data for parkinson disease.
as expected, our findings represent a superset of those originally described  <cit> .
snp rs <dig> is  <dig> kb upstream of the gene ngly <dig> , which encodes an enzyme thought to participate in proteasomal degradation of misfolded glycoproteins  <cit> .
snp rs <dig> is intronic in the gene gak , which regulates clathrin-mediated membrane trafficking  <cit> .
in the original analysis of these data, gak was found to be associated under an additive model  <cit> , although our analysis revealed that the association at this locus resulted from a dominant effect.
linkage disequilibrium levels are moderate to strong over a  <dig> kb region including gak as well as five other genes, leading to poor resolution at this locus.
snps rs <dig> and rs <dig> are intronic in the gene syn <dig> .
synapsins are essential for assembly of synaptic vesicles and modulate neurotransmitter release, with expression of syn <dig> being neuron-specific  <cit> .
as with other genome-wide association studies, the power of our analysis was limited by small effect sizes and small sample sizes, particularly a small replication sample.
we have developed a unified framework using unions of multi-locus genotypes to powerfully, flexibly, and efficiently analyze genome-wide genotype data, candidate gene data, or medical sequencing data.
implementation with standard statistical packages  is straightforward.
we replicated as significantly affecting pd risk a single, low-frequency haplotype in the region of the gak gene on chromosome  <dig> and identified two other regions with suggestive evidence of multiple low-frequency risk haplotypes.
CONCLUSIONS
recent attempts at elucidating the genetic architecture of complex traits have been hampered by the use of methods that do not account for correlation in the data, are overly stringent, and are underpowered in the presence of allelic or locus heterogeneity and rare variants.
to address these issues, we developed a powerful and efficient multivariable method based on unions of variants.
our method produces multivariable test statistics with reduced degrees of freedom compared to haplotype-based methods, implicitly accounts for linkage disequilibrium, and reduces the testing burden.
our method also generalizes collapsing methods as previously described for analysis of rare variants.
thus, we provide a unified framework for the analysis of one or more variants with any pattern of linkage disequilibrium and with any minor genotype frequency.
our method is suitable for genome-wide genotype data, candidate gene data, exome sequencing data, and whole genome sequencing data.
using our new method, we found evidence supporting susceptibility to pd at three loci.
authors' contributions
ds designed the study, analyzed the data, and drafted the manuscript.
lkv conceived of the study and contributed to the interpretation of the results and to the writing of the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
table comparing the validity of the union test and the kernel machine-based test.
see the legend for figure 1a for details.
click here for file
additional file 2
table comparing the power of the union test and the kernel machine-based test for multiple independent markers with moderate effects.
see the legend for figure 1b for details.
click here for file
additional file 3
table comparing the power of the union test and the kernel machine-based test for multiple independent markers with larger effects.
see the legend for figure 1c for details.
click here for file
additional file 4
table comparing the power of the union test and the kernel machine-based test for markers with opposing effects.
see the legend for figure 1d for details.
click here for file
additional file 5
table comparing the power of the union test and the kernel machine-based test for dilute signal.
see the legend for figure 1e for details.
click here for file
additional file 6
table comparing the power of the union test and the kernel machine-based test for correlated predictors.see the legend for figure 1f for details.
click here for file
additional file 7
table comparing the power of the union test and the kernel machine-based test for epistasis.
see the legend for figure 1g for details.
click here for file
additional file 8
table comparing the power of the union test and the kernel machine-based test for differential correlation between markers.
see the legend for figure 1h for details.
click here for file
additional file 9
quantile-quantile plots for genomic control.
a) the discovery sample.
b) the replication sample.
the red lines indicate the expected distribution.
the inflation factors are shown, indicated by λgc.
click here for file
acknowledgements
we dedicate this manuscript to allan morris and fern vaughan.
we thank david t. redden and charles rotimi for discussion and critical review of this manuscript.
funding support was provided by the national institute on aging, the national institute of neurological disorders and stroke, the robert p. & judith n. goldberg foundation, the bumpus foundation, and the harvard neurodiscovery center.
samples were provided by the national institute of neurological disorders and stroke human genetics resource center dna and cell line repository.
genotyping services were provided by the center for inherited disease research and the national institute of neurological disorders and stroke .
the discovery sample used for the analyses described in this manuscript was obtained from the ninds database found at http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?id=phs <dig> v <dig> p <dig> through dbgap accession number phs <dig> v <dig> p <dig>  the replication sample used for the analyses described in this manuscript was obtained from the coriell queue, now available from the ninds database found at http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs <dig> v <dig> p <dig>  the contents of this publication are solely the responsibility of the authors and do not necessarily reflect the official view of the national institutes of health.
this research was supported by the intramural research program of the center for research on genomics and global health , which is supported by the national human genome research institute, the national institute of diabetes and digestive and kidney diseases, the center for information technology, and the office of the director at the national institutes of health  and national institutes of health grant k01dk <dig>  the funding bodies had no role in study design; in the collection, analysis, and interpretation of data; in the writing of the manuscript; and in the decision to submit the manuscript.

 mapping disease susceptibility loci for complex diseases has proven to be challenging. the etiology of a complex disease is multifactorial, with both genetic and environmental risk factors and interactions among these risk factors. genome-wide association studies are an alternative to traditional family-based linkage studies for mapping disease susceptibility loci for complex diseases. the most common design to date for gwa studies is to genotype a dense panel of single-nucleotide polymorphisms for a large collection of unrelated cases and controls. snps are then tested one marker at a time for association most often using contingency table analysis, i.e., a x2test or fisher's exact test with one or two degrees of freedom, depending on the hypothesis being tested. the resulting set of p-values is then adjusted for multiple comparisons, usually a bonferroni correction for the number of tests performed. single marker analysis involves the marginalization of effects over all genetic and environmental backgrounds. consequently, effect sizes for single markers tend to be small, necessitating sample sizes of thousands to obtain reasonable levels of power at genome-wide significance levels. multi-locus methods may have increased power if they account for the correlation of loci due to linkage disequilibrium. however, multi-locus methods generally have been observed to have decreased power because they tend to produce test statistics with larger degrees of freedom. therefore, in designing a powerful and efficient multi-locus method for generating testable hypotheses based on genome-wide snp data, we are motivated to reduce the degrees of freedom. alternatively, consider a biochemical pathway consisting of many genes, which may be located on different chromosomes and hence may segregate independently. this motivation raises the general question of how to test a set of markers that may be biologically related but not correlated via linkage disequilibrium. possibilities include multivariate techniques such as the f test or multiple regression. however, these techniques can also lose power due to large degrees of freedom. a third issue for association analysis is low power for rare variants that medical sequencing is expected to uncover. to overcome low power for variants with frequencies < 5%, specialized methods for grouping variants have been described, but such methods were designed only for analysis of rare variants for which linkage disequilibrium is expected to be negligible. furthermore, the two most commonly used estimators of linkage disequilibrium, r and d', are biased upwards when allele frequencies are low or sample sizes are small. in response to all of these motivations, we propose a flexible multi-locus method based on unions of multiple snps that yields a test that always has only one degree of freedom. the method implicitly accounts for linkage disequilibrium and is appropriate for simultaneous analysis of common and rare variants. we demonstrate by simulation that our method has valid control of the false positive error rate while yielding more power than traditional single marker analysis. we also demonstrate by simulation that our method is sensitive to any source of differential correlation among markers. we then analyze two publicly available data sets for parkinson disease for which no significant associations were detected by single marker analysis. using our new method, we found evidence supporting susceptibility to pd at three loci. methods study samples parkinson disease is the second most common neurodegenerative disorder after alzheimer disease. for the discovery sample, we obtained publicly available data from the progeni and genepd studies of familial pd. briefly, cases and controls were genotyped using illumina humancnv version1_c beadchips and the illumina infinium ii assay protocol. for the replication sample, we obtained publicly available data from a study of sporadic pd from the snp resource at the national institute of neurological disorders and stroke human genetics resource center dna and cell line repository http://ccr.coriell.org/ninds/ as well as clinical data. the original genotyping was performed in the laboratory of drs. andrew singleton and john a. hardy, bethesda, md usa. briefly, snps were genotyped using the illumina infinium i and humanhap assays on cases and neurologically normal controls. the data made available had already been processed for quality control. the framework of logical unions to illustrate the construction of a union based on classical logic, consider two events. let p be the probability of event x and p be the probability of event y. the probability of the union of events x and y, p, is p = p + p - p. to estimate the probability of the intersection p, we need a measure of correlation. one commonly used measure of correlation is the pairwise correlation coefficient, given by r=covvarvar=p−ppp)p). after rearrangement and substitution, p=p+p−pp−rp)p). in the genetic context, let a represent the major allele and b represent the minor allele at a snp. with respect to the minor allele, dominant coding of the three genotypes is given by and recessive coding is given by. additive coding of the number of copies of the minor allele, i.e.,, can be achieved by summing over dominant and recessive variables since + = and thus need not be separately considered. r is a commonly used measure of linkage disequilibrium. thus, unions implicitly account for correlation between markers, including coupling or repulsion phase. it is critical to note that this use of the correlation coefficient encompasses gametic phase disequilibrium as well as other sources of correlation such as epistasis and natural selection. also note that, by definition, p = p, demonstrating the equivalence of analysis of a union of one marker with traditional single marker analysis. there are many possible coding schemes one could implement within this framework. in this study, we construct unions over consecutive, non-overlapping sets of k markers. the choice of non-overlapping sets was made to facilitate calculating the genome-wide testing burden. let xij be a binary indicator variable for the ith individual and the jth union. with two markers under dominant coding for the minor allele, xij={0for genotypea1a1a2a21for genotypesa1b1a2a b1b1a2a a1a1a2b a1a1b2b a1b1a2b a1b1b2b b1b1a2b b1b1b2b21for genotypesa1b1?2? b1b1?2??1?1a2b?1?1b2b2?for genotypesa1a1?2??1?1a2a?1?1?2? in which? represents missing data. under dominant coding, a value of indicates the absence of at least one b allele and a value of indicates the presence of at least one b allele. similarly, with two markers under recessive coding for the minor allele, xij={0for genotypesa1a1a2a a1a1a2b a1b1a2a a1b1a2b21for genotypesa1a1b2b a1b1b2b b1b1a2a b1b1a2b b1b1b2b21for genotypesb1b1?2??1?1b2b2?for genotypesa1a1?2? a1b1?2??1?1a2a?1?1a2b?1?1?2? under recessive coding, a value of indicates the absence of at least one bb genotype and a value of indicates the presence of at least one bb genotype. any individual for which xij =? is removed from analysis of the jth union. note that this choice of implementation relies solely on counting individuals with pre-specified multi-locus genotypes, does not require estimation of pi, pj, or r, and does not assume hardy-weinberg equilibrium. for case-control data, we test for association using fisher's exact test. for the jth union, the × contingency table is constructed by counting the numbers of individuals for which xij is or among cases and controls. the test has only one degree of freedom, regardless of k. we control for multiple comparisons using a bonferroni correction for the number of unions tested. the number of unions tested is given by ⌈#snpsk⌉, in which ┌x┐ is the ceiling function and returns the smallest integer not less than x. this implementation of multiplicity control preserves more power as k increases and as the testing burden is consequently reduced. given snps in the discovery sample and accounting for both dominant and recessive coding over unions of size one to five, the genome-wide significance level was ×∑k=15⌈ k⌉= ×10− all union testing was performed in r. other analyses imputation was performed using mach, version. http://www.sph.umich.edu/csg/abecasis/mach/download/. for the reference panel, we retrieved the combined hapmap phase ii+iii raw genotype files for the ceu sample http://hapmap.ncbi.nlm.nih.gov/downloads/genotypes/latest_phaseii+iii_ncbi_b36/forward/non-redundant/. we filtered the, autosomal ceu snps based on the inclusion of unrelated individuals only, a minor allele frequency ≥, a snp missingness rate ≤ 5%, and an individual missingness rate ≤ 5%. we inferred haplotype phases for the reference data using the settings --rounds --states we conditioned imputation on the maximum-likelihood estimates of the crossover map and the error rate map. we retained all imputed genotype calls for which the posterior probability ≥. we performed haplotype analysis using plink version http://pngu.mgh.harvard.edu/purcell/plink/. briefly, we compared each haplotype against all of the other haplotypes.  simulation analysis we first investigated the validity and power of our proposed method, exploring over sizes of unions of snps ranging from one to five. note that analysis of a union of one snp is identical to single marker analysis. under the null hypothesis of no association, analysis of unions has the expected per comparison error rate and is therefore valid. for small effect sizes, analysis of unions is increasingly more powerful than single marker analysis as the frequency of the minor multi-locus genotype decreases. if a union consists of both risk-increasing and risk-decreasing predictors, effects will cancel and power will be lost. also, if a union consists of too many predictors with no effect on the outcome, then power to detect a predictor with an effect within the union will decrease. analysis of unions of correlated predictors can be more powerful than analysis of unions of independent predictors. epistasis can either increase or decrease power, depending on the directions of the effect sizes for epistatic effects vs. marginal effects. for epistasis to be detectable, minor genotype frequencies must be large so that the joint genotype counts are reasonably large. the union test is a valid test of association of genotype frequencies if the correlation between markers is matched in cases vs. controls. the union test is also a powerful test of differential correlation among markers, even in the absence of differences in genotype frequencies. thus, the union test can simultaneously detect differences in genotype frequencies, epistasis, haplotype structure, natural selection, population structure, and any other process affecting correlation among markers. comparison to other methods we compared the union test to a logistic kernel machine-based test, a set-based test for common variant analysis that adaptively estimates the degrees of freedom given the correlational structure of the markers in the set. both tests were valid at the significance level. the union test was more powerful than the kernel machine-based test for multiple independent markers with nonzero effects. the union test was also more powerful than the kernel machine-based test when the minor genotype frequency was below ~. the kernel machine-based test was more powerful than the union test for markers with opposing effects, dilute signal, correlated predictors, and epistasis. the kernel machine-based test was conservative in the presence of differential correlation between markers under the null hypothesis and had no power to detect differential correlation under the alternative hypothesis. we also compared the union test to the collapsing method, a set-based test for rare variant analysis that collapses the data into a x test with one degree of freedom. the collapsing method is based on aggregating genotypes across all loci in the set such that "an individual is coded as if a rare allele is present at any of the variant sites and as otherwise". thus, the collapsing method is in fact identical to the union test by construction and power is equivalent for both tests. the two main differences between the collapsing method and the union test are that the former assumes that minor allele frequencies are low, i.e., ≤, and that all markers are independent, i.e., there is no linkage disequilibrium. the collapsing method is thus a special case of the union test, with the union test applicable across the entire range of frequencies and across the entire range of linkage disequilibrium. real data analysis data processing for quality control for both samples is depicted in figure for the discovery sample, we retained of the cases, of the controls, and of the snps. for the replication sample, we retained of the cases, of the controls, and of the snps. to investigate the possibility of population stratification in the discovery and replication samples, we estimated the variance inflation factor of the genomic control method. we estimated an inflation factor of for the discovery sample and for the replication sample, indicating a negligible inflation of the false positive error rate due to population stratification. under dominant coding, we identified no regions correlated with pd case/control status at a genome-wide significance level of × 10-. we identified seven independent loci at suggestive levels of association. of these, we replicated the region on chromosome 4p including the gene gak. this region was discovered using unions of one, two, or three markers but was replicated only in the single marker analysis. using haplotype analysis, we found that the association in both samples was driven by the same low-frequency haplotype with an odds ratio of. a in order, the snps are rs and rs to directly compare single marker analysis, haplotype analysis, and our multi-locus method, we examined the first multi-locus union in table which consisted of snps rs and rs using the replication sample. by × contingency table analysis of genotypes, rs was significantly associated with pd but rs was not. similarly, by × contingency table analysis under dominant coding, rs was significantly associated with pd but rs was not. haplotype analysis revealed significant association with pd for the haplotype consisting of the minor allele at rs and the major allele at rs i.e., haplotype tt vs. haplotypes cc and ct with haplotype tc being unobserved. under dominant coding, × contingency table analysis of the union did not yield significant association, because the union tested multi-locus genotype cctt vs. the other eight genotypes, effectively attenuating the signal of haplotype tt. in general, for a union over k snps, the number of joint genotypes is 3k and × contingency table analysis is more powerful than × 3k contingency table analysis. more importantly, as k increases, some joint genotypes likely will be unobserved, leading to cells in the contingency table having counts of zero. collapsing a × 3k contingency table into a × contingency table is one way to address this problem of data sparsity. in our implementation of dominant coding over multiple loci, the reference multi-locus genotype consists solely of diplotypes of the most common haplotype. this example illustrates the trade-off of gaining power by reducing the degrees of freedom vs. losing power by attenuation of signal. under recessive coding, we identified no regions correlated with pd case/control status at a genome-wide significance level of × 10-. at suggestive significance levels, we identified two potentially interesting loci at chromosomes 3p and 22q. at chromosome 3p snp rs was associated with decreased risk in the discovery sample but with increased risk in the replication sample. we investigated this directional inconsistency locus using haplotype analysis. we identified a haplotype block of eight snps including rs in the discovery sample, the haplotype with the strongest effect carried the minor allele at rs and was associated with decreased risk. in the replication sample, the haplotype with the strongest effect carried the major allele and was associated with increased risk. thus, the effects were directionally consistent at the haplotypic level with respect to rs but were directionally inconsistent with respect to marginal effects at this snp. at chromosome 22q the marginal effects for snp rs were directionally inconsistent and different haplotypes drove association in the two samples. a in order, the snps are rs rs rs rs rs rs rs and rs a in order, the snps are rs and rs finally, we examined execution time for the union test. algorithmically, the union test consists of one hypothesis test per union, regardless of the number of markers in the union. for a union of k markers, the union test is approximately k- times faster than the corresponding single marker test, because the computational burden of taking the union is less than the computational burden of performing one fisher's exact test. discussion recent attempts at elucidating the underlying genetic structure of common, complex diseases have been hampered by the use of methods that do not account for correlation in the data, are overly stringent, and are underpowered in the presence of allelic or locus heterogeneity and rare variants. compared to other multi-locus methods designed to address these issues, unions offer many advantages, including accounting for linkage disequilibrium, epistatic interactions within unions, and haplotype structure, with no requirement for the estimation of haplotypes or phase. union testing can be viewed as a generalization of grouping schemes for analysis of independent rare variants that allows for arbitrary patterns of correlation among common and/or rare variants. an additional advantage is that the counting algorithm implicitly accounts for correlation without explicit and potentially biased estimation of linkage disequilibrium measures and without the need of permutation testing to establish correct type i error rates. union tests can be more powerful than single marker tests because they make fuller use of the information in a multivariable framework. power may also be gained by collapsing multi-locus genotype information into a single degree of freedom test. the union test is a composite test of association of genotype frequencies and differential correlation among markers. union tests are also haplotype tests and are sensitive to processes such as natural selection and population structure that may lead to differential correlation among markers. previous analyses did not detect substantial population structure in either pd data set. however, given that our multi-locus method is more powerful than single marker analysis under many genotypic configurations, it is possible that our method is more sensitive to cryptic relatedness or residual population stratification than is single marker analysis. to better address this issue, regression with a covariate for parental ancestry provides a framework compatible with unions. our proposed method is very flexible. the choice of which snps to combine into unions belongs to the investigator. we chose to combine consecutive snps in non-overlapping sets in order to take advantage of linkage disequilibrium, so that union testing may be viewed as a generalization of haplotype testing, while minimizing the correlation between tests. with an appropriate grouping scheme, unions may be used to address allelic heterogeneity, which occurs when the phenotype of interest is caused by different mutations within the same gene. the role of allelic heterogeneity in complex disease is exemplified by the several, individually rare mutations in nod associated with crohn's disease. alternatively, one could focus on rare variants and group snps with minor allele frequencies < 5%. rare variants have been predicted to play an important role in complex diseases, as has been illustrated for colorectal cancer and cholesterol metabolism. unions may also be used to address genetic heterogeneity, which occurs when multiple genes cause the same disease. although we tested unions using contingency table analysis, it is straightforward to adapt our coding for use with generalized linear models and thus to incorporate covariates into the analysis. our method is also flexible with respect to genetic models, such that coding for snps may be any combination of dominant, additive, or recessive. depending on the implementation, union testing can also reduce the testing burden compared to single marker analysis. however, instead of non-overlapping windows as we used herein, a sliding window approach could be employed, at the cost of a higher testing burden. there are two main disadvantages of unions. one disadvantage is that the snp and/or haplotype within a union driving association cannot be inferred. if genome-wide analysis is intended as exploratory and hypothesis-generating, then this disadvantage is not a concern. the second disadvantage is that there are scenarios under which statistical power may be lost, depending on how collapsing is achieved. as with other collapsing methods, the power of union testing decreases if the size of the union is too large, whether or not the included variants are causal. this limitation is more pronounced for markers with large minor allele frequencies. similarly, the power of union testing decreases if both risk-increasing and risk-decreasing effects are included within a union. one approach to address both attenuation and cancellation of signal is to use a mixture model with one component for risk-increasing effects, one component for risk-decreasing effects, and one component for null effects. we illustrated union testing by reanalyzing genome-wide genotype data for parkinson disease. as expected, our findings represent a superset of those originally described. snp rs is kb upstream of the gene ngly, which encodes an enzyme thought to participate in proteasomal degradation of misfolded glycoproteins. snp rs is intronic in the gene gak, which regulates clathrin-mediated membrane trafficking. in the original analysis of these data, gak was found to be associated under an additive model, although our analysis revealed that the association at this locus resulted from a dominant effect. linkage disequilibrium levels are moderate to strong over a kb region including gak as well as five other genes, leading to poor resolution at this locus. snps rs and rs are intronic in the gene syn. synapsins are essential for assembly of synaptic vesicles and modulate neurotransmitter release, with expression of syn being neuron-specific. as with other genome-wide association studies, the power of our analysis was limited by small effect sizes and small sample sizes, particularly a small replication sample. we have developed a unified framework using unions of multi-locus genotypes to powerfully, flexibly, and efficiently analyze genome-wide genotype data, candidate gene data, or medical sequencing data. implementation with standard statistical packages is straightforward. we replicated as significantly affecting pd risk a single, low-frequency haplotype in the region of the gak gene on chromosome and identified two other regions with suggestive evidence of multiple low-frequency risk haplotypes.  recent attempts at elucidating the genetic architecture of complex traits have been hampered by the use of methods that do not account for correlation in the data, are overly stringent, and are underpowered in the presence of allelic or locus heterogeneity and rare variants. to address these issues, we developed a powerful and efficient multivariable method based on unions of variants. our method produces multivariable test statistics with reduced degrees of freedom compared to haplotype-based methods, implicitly accounts for linkage disequilibrium, and reduces the testing burden. our method also generalizes collapsing methods as previously described for analysis of rare variants. thus, we provide a unified framework for the analysis of one or more variants with any pattern of linkage disequilibrium and with any minor genotype frequency. our method is suitable for genome-wide genotype data, candidate gene data, exome sequencing data, and whole genome sequencing data. using our new method, we found evidence supporting susceptibility to pd at three loci. authors' contributions ds designed the study, analyzed the data, and drafted the manuscript. lkv conceived of the study and contributed to the interpretation of the results and to the writing of the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 table comparing the validity of the union test and the kernel machine-based test. see the legend for figure 1a for details. click here for file additional file 2 table comparing the power of the union test and the kernel machine-based test for multiple independent markers with moderate effects. see the legend for figure 1b for details. click here for file additional file 3 table comparing the power of the union test and the kernel machine-based test for multiple independent markers with larger effects. see the legend for figure 1c for details. click here for file additional file 4 table comparing the power of the union test and the kernel machine-based test for markers with opposing effects. see the legend for figure 1d for details. click here for file additional file 5 table comparing the power of the union test and the kernel machine-based test for dilute signal. see the legend for figure 1e for details. click here for file additional file 6 table comparing the power of the union test and the kernel machine-based test for correlated predictors.see the legend for figure 1f for details. click here for file additional file 7 table comparing the power of the union test and the kernel machine-based test for epistasis. see the legend for figure 1g for details. click here for file additional file 8 table comparing the power of the union test and the kernel machine-based test for differential correlation between markers. see the legend for figure 1h for details. click here for file additional file 9 quantile-quantile plots for genomic control. a) the discovery sample. b) the replication sample. the red lines indicate the expected distribution. the inflation factors are shown, indicated by λgc. click here for file acknowledgements we dedicate this manuscript to allan morris and fern vaughan. we thank david t. redden and charles rotimi for discussion and critical review of this manuscript. funding support was provided by the national institute on aging, the national institute of neurological disorders and stroke, the robert p. & judith n. goldberg foundation, the bumpus foundation, and the harvard neurodiscovery center. samples were provided by the national institute of neurological disorders and stroke human genetics resource center dna and cell line repository. genotyping services were provided by the center for inherited disease research and the national institute of neurological disorders and stroke. the discovery sample used for the analyses described in this manuscript was obtained from the ninds database found at http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?id=phs v p through dbgap accession number phs v p the replication sample used for the analyses described in this manuscript was obtained from the coriell queue, now available from the ninds database found at http://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs v p the contents of this publication are solely the responsibility of the authors and do not necessarily reflect the official view of the national institutes of health. this research was supported by the intramural research program of the center for research on genomics and global health, which is supported by the national human genome research institute, the national institute of diabetes and digestive and kidney diseases, the center for information technology, and the office of the director at the national institutes of health and national institutes of health grant k01dk the funding bodies had no role in study design; in the collection, analysis, and interpretation of data; in the writing of the manuscript; and in the decision to submit the manuscript.
to illustrate use of the union-intersection test with real data, we analyzed a publically available data set of autosomal snps genotyped for cases of parkinson disease and neurologically normal controls for which no genome-wide significant results were found by traditional single marker analysis. within the framework, we developed a union-intersection test suitable for genome-wide analysis of single nucleotide polymorphisms, candidate gene data, as well as medical sequencing data. the union-intersection test is a composite test of association of genotype frequencies and differential correlation among markers. we also demonstrated scenarios in which the multi-locus test was more powerful than traditional single marker analysis. we also found suggestive evidence that directionally inconsistent marginal effects from single marker analysis appeared to result from risk being driven by different haplotypes in the two samples for the genes syn and ngly which are involved in neurotransmitter release and proteasomal degradation, respectively. these results illustrate the utility of our unified framework for genome-wide association analysis of common, complex diseases. we identified a single risk haplotype with a directionally consistent effect in both samples in the gene gak, which is involved in clathrin-mediated membrane trafficking. common, complex diseases are hypothesized to result from a combination of common and rare genetic variants.
 common, complex diseases are hypothesized to result from a combination of common and rare genetic variants. we developed a unified framework for the joint association testing of both types of variants. within the framework, we developed a union-intersection test suitable for genome-wide analysis of single nucleotide polymorphisms, candidate gene data, as well as medical sequencing data. the union-intersection test is a composite test of association of genotype frequencies and differential correlation among markers.  we demonstrated by computer simulation that the false positive error rate was controlled at the expected level. we also demonstrated scenarios in which the multi-locus test was more powerful than traditional single marker analysis. to illustrate use of the union-intersection test with real data, we analyzed a publically available data set of autosomal snps genotyped for cases of parkinson disease and neurologically normal controls for which no genome-wide significant results were found by traditional single marker analysis. we also analyzed an independent follow-up sample of cases and controls for replication.  we identified a single risk haplotype with a directionally consistent effect in both samples in the gene gak, which is involved in clathrin-mediated membrane trafficking. we also found suggestive evidence that directionally inconsistent marginal effects from single marker analysis appeared to result from risk being driven by different haplotypes in the two samples for the genes syn and ngly which are involved in neurotransmitter release and proteasomal degradation, respectively. these results illustrate the utility of our unified framework for genome-wide association analysis of common, complex diseases. 