 Methodology article
 Open Access
 Published:
Combining controls can improve power in twostage association studies
BMC Geneticsvolume 19, Article number: 89 (2018)
Abstract
Background
High dimensional case control studies are ubiquitous in the biological sciences, particularly genomics. To maximise power while constraining cost and to minimise type1 error rates, researchers typically seek to replicate findings in a second experiment on independent cohorts before proceeding with further analyses. This can be an expensive procedure, particularly when control samples are difficult to recruit or ascertain; for example in interdisease comparisons, or studies on degenerative diseases.
Results
This paper presents a method in which control (or case) samples from the discovery cohort are reused in a replication study. The theoretical implications of this method are discussed and simulated genomewide association study (GWAS) tests are used to compare performance against the standard approach in a range of circumstances.
Using similar methods, a procedure is proposed for ‘partial replication’ using a new independent cohort consisting of only controls. This methods can be used to provide some validation of findings when a full replication procedure is not possible.
The new method has differing sensitivity to confounding in study cohorts compared to the standard procedure, which must be considered in its application. Type1 error rates in these scenarios are analytically and empirically derived, and an online tool for comparing power and error rates is provided.
Conclusions
In several common study designs, a sharedcontrol method allows a substantial improvement in power while retaining type1 error rate control. Although careful consideration must be made of all necessary assumptions, this method can enable more efficient use of data in GWAS and other applications.
Background
Highdimensional casecontrol studies have become a mainstay of investigation of pathophysiology in complex diseases and traits. An important part of their analysis is the process of replication [1], in which the results of a highdimensional study are used to inform the design of a second study at a subset of the original variables, with a joint analysis used to determine overall association.
Replicating studies in this way has the advantage of increasing the effective study sample sizes without requiring measurement of all variables in all samples. It also serves to protect against falsepositives due to systematic errors in the original datasets, by retesting association in a second nominally independent dataset.
Replication has a significant cost, and can require large numbers of samples, especially when associated variables have small effects (i.e. [2]). There is therefore a need to minimise the number of additional samples which need to be analysed. This paper presents a method to perform replication by combining controls in both the original ‘discovery’ and second ‘replication’ datasets, potentially reducing the number of new samples required. Sharedcontrol approaches can improve study efficiency in many related applications in which studies are compared [3–8].
Results from original and replication datasets for which some or all controls are shared cannot be directly compared due to the correlation between test statistics directly resulting from shared controls even under the null hypothesis [5]; use of the same thresholds in a sharedcontrol design as used in an independentcontrols design will lead to higher type1 error rates. This paper demonstrates a simple adaptation to a standard design to account for the changed correlation structure and retain control of type1 error rate, only requiring a change to one pvalue threshold.
An important purpose of replication is control against falsepositives arising from variables for which confounding causes an apparent casecontrol difference in one of the discovery or replication phase experiments, but not the other. The action of sharing control samples results in a different spectrum of sensitivity to variables of this type. It necessitates a sacrifice of type1 error rate control in variables for which confounding affects the discoveryphase control cohort, but improves type1 error rate control in variables for which confounding affects the replicationphase control cohort. The type1 error rate is largely equivalent to an independentcontrols design in variables affected by confounding in either case cohort.
The new spectrum of false positive rates can be advantageous in circumstances where control samples in the replication cohort are less wellascertained than those in the discovery cohort. This may be the case in studies on degenerative disease, where control ascertainment is generally uncertain, and populationsourced controls may be used for replication. The sharedcontrol design can reduce power losses from misspecified controls in the replication cohort, as well as reducing falsepositive rates caused by confounding in the cohort.
When used with shared cases instead of controls, this method can be adapted to a ‘partial replication’ procedure where only a new control set is used. Although not equivalent to a full replication in an independent dataset, the procedure enables improvement in type1 error rates and control over confounding. This is applicable in studies on rare traits, where all available samples need to be included in the discovery analysis for adequate power.
Throughout this paper, GWAS terminology will be used (singlenucleotide polymorphisms (SNPs), allele frequency, variants etc) although the method is applicable to any highdimensional case control study. ‘Controls’ will be considered to generally be samples unaffected by a disease or trait of interest, although the method can be applied with case/control labels swapped, or applied to comparisons between subgroups of a case group.
Differences in power (at fixed type1 error rate) between standard (independentcontrols) and new (sharedcontrol) methods are established by considering hypothesis tests typical of those found in GWAS. Asymptotic analytical results are established where possible, but all type 1/type 2 error rates are readily tractable empirically to good accuracy given study sizes and proposed pvalue thresholds, and a tool is provided to do this at https://wallacegroupliley.shinyapps.io/replication_shared/.
Results
Overview of method
We assume a GWAS dataset of a set of cases C_{1} and controls C_{0} used in a ‘discovery’ phase of a GWAS or similar study, and corresponding sets of cases and controls \(C_{1}^{\prime }\), \(C_{0}^{\prime }\) in the replication phase. We assume that C_{0} and C_{1} are genotyped at a set of SNPs S and \(C_{0}^{\prime }\), \(C_{1}^{\prime }\) at a set S^{′}⊆S.
For each SNP we designate μ_{1}, μ_{0}, \(\mu _{1}^{\prime }\), \(\mu _{0}^{\prime }\) as the population minor allele frequency in the corresponding group, and m_{1}, m_{0}, \(m_{1}^{\prime }\), \(m_{0}^{\prime }\) as the observed allele frequency (so E(m_{i})=μ_{i}). We designate two null hypotheses; \(H_{0}^{\cup }:(\mu _{1}=\mu _{0}) \cup (\mu _{1}^{\prime }=\mu _{0}^{\prime })\) and \(H_{0}^{=}:(\mu _{1}=\mu _{0} = \mu _{1}^{\prime }=\mu _{0}^{\prime })\), noting that \(H_{0}^{\cup } \supseteq H_{0}^{=}\). In a typical conservative GWAS approach, we seek to test against \(H_{0}^{\cup }\), since μ_{1}≠μ_{0} or \(\mu _{1}^{\prime } \neq \mu _{0}^{\prime }\) may hold at nondisease associated SNPs due to confounding in the original or replication studies respectively. The alternative null hypothesis \((\mu _{1}=\mu _{0} \cap \mu _{1}^{\prime }=\mu _{0}^{\prime })\), which implies \(H_{0}^{=}\) and is implied by \(H_{0}^{\cup }\), is more appropriate than \(H_{0}^{=}\) in cases where replication is performed in a different population than discovery. However, this situation is not adaptable to a sharedcontrol design.
A typical twostage genetic testing procedure [9], which we will refer to as method A, begins by comparing genotypes of C_{1} and C_{0} at SNPs S generating pvalues p_{d} (discovery). A subset S^{′} of SNPs reaching putative significance level p_{d}<α are genotyped in \(C_{0}^{\prime }\) and \(C_{1}^{\prime }\), with genotypes compared to generate pvalues p_{r} (replication stage). Finally, genotypes are compared between \(C_{0} \cup C_{0}^{\prime }\) and \(C_{1} \cup C_{1}^{\prime }\) at SNPs S^{′} to generate pvalues p_{m} (metaanalytic stage). SNPs are designated as ‘hits’ if p_{d}<α,p_{r}<β,p_{m}<γ for some β, γ, and all effects have the same direction. The values α, β, γ may not be explicitly stated in some study designs, although they are usually implicitly present. This is discussed further in the “Choice of thresholds” section below.
The main modification proposed in this paper, denoted as method B, differs at the replication stage in that \(C_{1}^{\prime }\) is compared with \(C_{0} \cup C_{0}^{\prime }\) at S^{′} instead of just C_{0} (Fig. 1). The pvalues resulting from the modified replication stage are termed p_{s}, and the criterion to designate a hit changed to p_{d}<α,p_{s}<β^{∗},p_{m}<γ, with all effects in the same direction. The threshold β^{∗} is chosen to conserve type1 error rate between methods (see “Methods” section and Additional file 1: Appendix 1). This requires estimation of systematic correlation between Z scores, which may be estimated either empirically or (in some cases) analytically.
A second modification, denoted method C, combines C_{0} and \(C_{0}^{\prime }\) at both the discovery and replication phase (see Fig. 1). This is analogous to a situation in which only a single control cohort is available, and a choice must be made to split it between discovery and replication procedures or to use it for both. In this case, \(C_{0} \cup C_{0}^{\prime }\) is compared with C_{1} at SNPs S in the discovery phase to produce pvalues p_{c}, then \(C_{0} \cup C_{0}^{\prime }\) is compared with \(C_{1}^{\prime }\) at SNPs S^{′} at the replication phase and compared with \(C_{1} \cup C_{1}^{\prime }\) at the metaanalytic stage to produce pvalues p_{s} and p_{m} as before. A hit is determined by p_{c}<α,p_{s}<β^{⊥},p_{m}<γ, with all effects in the same direction. Again, β^{⊥} is chosen to maintain the type1 error rate between methods.
General properties
For SNPs in \(H_{0}^{=}\), the overall type1 error rate is conserved between methods by the definition of β^{∗}, β^{⊥} (Eq. 4) at a level P_{0}. It is shown in Additional file 1: Appendix 2.2 that β>β^{∗}>β^{⊥}. For SNPs in \(H_{0}^{\cup } \setminus H_{0}^{=}\) the type1 error rates differ between methods. Such SNPs may be characterised by the group(s) amongst C_{0}, C_{1}, \(C_{0}^{\prime }\), \(C_{1}^{\prime }\) in which their expected minor allele frequency (MAF) is aberrant from the expected MAF in the population which the group ostensibly represents. ‘Aberrance’ is taken to mean an incorrect expected value from systematic measurement error or uncorrected confounding, rather than random deviance around a correct expected value.
Bounds on type1 error rates with aberrance in each group are shown in Table 1. Methods B and C necessitate sacrificing bounds on error rates with aberrance in C_{0} and \(C_{0},C_{0}^{\prime }\) respectively. The bound on error with aberrance in \(C_{1}^{\prime }\) improves through methods AC. In the “Methods” section, it is shown that the type1 error with aberrance in \(C_{0}^{\prime }\) decreases from methods A to B, and the error with aberrance in \(C_{1}^{\prime }\) increases from A through C, although the upper bound is the same for both.
Bias in effect size estimates
If a set of variants in a study are selected based on pvalue (either by ordering all pvalues and selecting some number, or by choosing all with a pvalue below some threshold), the observed casecontrol odds ratios at those variants are upwardsbiased when used as estimates of the true oddsratios of these variants between cases and controls in the population [10]. This bias is highest amongst variants for which the true logoddsratio is 0 (nonassociated).
A standard replication procedure can be considered as enabling an unbiased effect size estimate [11]; for nonassociated variants, this estimate has expectation 0. If controls are reused in the replication procedure, the estimate of effect size for associated variants from the replication procedure is no longer unbiased (since the original control samples are reused), and summary statistics from the replication procedure cannot be used directly as estimates of effect size (although estimates can still be made by considering summary statistics p_{d}, p_{r} if these can be calculated). After sharing controls, the effect size estimate for null variants in the replication procedure is similarly biased, and the adjustment β→β^{∗}/β^{⊥} corresponds to an adjustment for this effect.
Differences in power between methods
The power difference between methods B and A was analysed systematically by considering the behaviour of GWAS data across a range of values of \((n_{0},n_{1},n_{0}^{\prime },n_{1}^{\prime })\). In each calculation, genetic data was considered for a single common SNP with average minor allele frequency across cases and controls equal to 0.1, with a given effect size between cases and controls quantified by logodds ratio. Varying ratios \(n_{1}/n_{1}^{\prime }\), \(n_{0}^{\prime }/n_{1}^{\prime }\) were considered, with \(n_{0} + n_{0}^{\prime } + n_{1} + n_{1}^{\prime }\) held constant at 20,000 samples (Fig. 1).
At large effect sizes (in GWAS terms, large allelic differences between case and control cohorts) both methods have power approaching 1, so the difference is slight. Similarly, at very small effect sizes, both methods have power near zero. Since the only power differences are at moderate effect sizes, the main metric for power difference used in this paper was the average effect size difference (Fig. 2). Considering power of A and power of B as functions Power_{A}(x), Power_{B}(x) of an underlying logodds ratio x, the average power difference was defined as
The maximum power difference:
was also considered.
Figure 3 shows average power difference at various study sizes for typical α, β, γ values (α=5×10^{−6}, β=5×10^{−4}, γ=5×10^{−8}). The difference is typically highest when the ratio of controls to cases is high in the discovery cohort and low or equal in the replication cohort, and the number of cases in the discovery cohort is larger than the number in the replication cohort. Power to detect SNPs in H_{1} is typically highest in method C, secondhighest in method B, and lowest in method A.
Recommended applications
To demonstrate areas where this approach is applicable, several examples are constructed or sourced from the GWAS field in which the procedure of sharing controls or cases will improve power or type1 error profile of the twostage testing procedure or enable some form of orthogonal replication to be performed.
Assumptions
In order to use method B or C, it must be assumed that cohort C_{0} and \(C_{0}^{\prime }\) are sampled from similar enough populations to be comparable to C_{1} and \(C_{1}^{\prime }\). A reasonable check on whether the method is appropriate is whether the cohorts C_{0} and \(C_{0}^{\prime }\) could be interchanged without compromising matching between cases and controls in the discovery or validation studies (possibly with the inclusion of strata or covariates in the genetic risk model). An important caveat of methods B and C is sacrifice of control over errors arising from aberrance in C_{0} (method B) or \(C_{0} \cup C_{0}^{\prime }\) (method C), so an assumption must be made that variables affected by confounding or measurement error in these cohorts are understood to be distinguishable from true associations by qualitycontrol measures only. Variants which are aberrant in the same direction in both discovery and control cohorts  that is, \(\text {sign}(\mu _{1}\mu _{0})=\text {sign}(\mu _{1}^{\prime }\mu _{0}^{\prime }) \neq 0\)  cannot be distinguished from true associations without the use of external data.
Posthoc assessment of all putative hits should be performed to check for genotyping errors [12] and assess whether the hit could have arisen from aberrance in C_{0}.
Conventional GWAS
Method B is applicable in several cases in large conventional GWAS, particularly when then ratio of controls to cases in the discovery cohort is larger than that in the replication cohort. In a relatively recent GWAS on rheumatoid arthritis [13] with comparable sample populations for discovery and replication cohorts, method B could be used to attain greater power than method A for a fixed type1 error rate. Assuming that summary statistics are wellapproximated by binomial tests of allelic differences (so covariates and strata used in computation of summary statistics have only small effects), the improvement in power is around 4% for SNPs with an oddsratio of 1.3, MAF 0.1, and is positive across all odds ratios. More than 2000 additional controls in \(C_{0}^{\prime }\) would be needed to increase power by this amount (Fig. 4, top left).
Small power advantages such as this may make minimal difference in a single study, although since they require no extra cost, are worth attaining if possible. The power of method B is generally considerably higher than method A when n_{0}>n_{1} and \(n_{0}^{\prime } \approx n_{1}^{\prime }\). Power advantages may be more substantial in some cases; for example, a study with \((n_{0},n_{1},n_{0}^{\prime },n_{1}^{\prime })=(15 000,5000,5000,5000)\), method B enables a power increase of up to 8% (Fig. 4, top right panel). To achieve comparable performance with method A, around 2000 additional controls would be necessary in the replication cohort. Method B with \((n_{0},n_{0}^{\prime })=(15000,5000)\) is also more powerful than method A would be if controls were divided equally between C_{0} and \(C_{0}^{\prime }\) (see Fig. 4, top right panel).
Difficult control ascertainment
An important application of the method presented in this paper is in studies for which ‘control’ samples are expensive or difficult to ascertain. This is often the case in comparative studies between disease subtypes. In such studies, sharing controls can improve power substantially, especially if a proportion of samples in the replication cohort are falsely assigned to the control cohort (see “Methods” section).
An international GWAS on frontotemporal dementia in 2014 [14] is an example in which sharing controls may be beneficial. The study had sample sizes \((n_{0},n_{1},n_{0}^{\prime },n_{1}^{\prime })=(4308,2154,5094,1372)\). Control samples in the discovery phase were assessed for current neurological disease, and were used in previous studies on Parkinson’s disease, indicating a high degree of reliability. Control samples in the replication phase were collected from the same geographic distribution as cases, but were not explicitly used in previous neurological studies, suggesting better control ascertainment amongst the discovery cohort.
In this study, sharing controls could allow for a more stronglyascertained control cohort, and reduce the effects of confounders affecting \(C_{1}^{\prime }\) (see Fig. 4, bottom left panel). At typical values α=1×10^{−4}, β=1×10^{−3}, γ=5×10^{−8}, power is nearly equivalent between the two methods assuming all controls are genuine. However, with 10% misascertainment in \(C_{1}^{\prime }\), the power advantage of method B is up to 5%. Given the nearidentical distribution of cases in the discovery and validation cohort, cases could alternatively be shared, leading to a power increase of up to 6%.
Prospective study design
Studies may be planned and powered with the assumption that samples may be shared. For certain restrictions on sample numbers, this can provide the potential for greater power than would be attainable by restricting to an independentcontrols design. For instance, if we seek to validate hits on a GWAS with 10,000 controls and 5000 cases, and can afford to genotype a further 10,000 samples, power is higher after recruiting 4000 additional controls and 6000 additional cases and sharing controls than can be achieved from any independentcontrol study design (Fig. 4, bottom right panel).
This may be a common scenario if controls are sourced from a known database rather than specifically recruited for the study.
Partial replication
In circumstances where case recruitment is difficult, as in studies of rare diseases, an assessment of replicability may be made by retesting results from a discovery phase with a new control set only. This can enable the use of control cohorts which only partially match the case cohort.
In a GWAS on pemphigus vulgaris [15], a rare disease primarily affecting individuals of Ashkenazi Jewish ethnicity, the discovery cohorts were sampled from Jewish populations, with age and population matched controls. Control cohorts were small (\((n_{0},n_{1},n_{0}^{\prime },n_{1}^{\prime })=(100,400,59,285)\)), potentially due to difficulty recruiting both ethnically and geographicallymatched controls.
Method C could be used in this instance to enable a larger control set and greater power. If a control cohort of Ashkenazi individuals could be assembled without requiring geographic matching with the case set, it would be inappropriate to use as a sole control cohort against the existing case cohort, due to the potential for geographic confounding. However, such a cohort could be used as either C_{0} or \(C_{0}^{\prime }\) in method C, with the existing ethnically and geographically matched controls serving as the other cohort. In this way, the power advantage of the larger cohort could be used while maintaining control over potential aberrance in the larger control group.
Method C enables computation of power and type1 error rates, and comparison to alternative designs with cases split into smaller independent discovery and validation cohorts (method A). Testing a case cohort against two separate control cohorts is almost always more powerful for a fixed type1 error rate than splitting the case cohort in two and performing method A (see Additional file 1: Figures S1 and S2).
Choice of thresholds
The designation of explicit thresholds α, β, γ in a twostage study may not appear to reflect many reallife designs, but in general most studies will use it in some form, even if the thresholds are not directly stated. Heuristically, α is used as an initial ‘triage’ step, to reduce data dimensionality, β (which is usually less stringent than α to allow for some regression to the mean in true associations) is used as a check, and γ is used as a definitive test for association amongst candidate variants.
Because studies are usually limited by cost or resources, a given number of variants are selected to pass through to the replication step, rather than following up all variants passing a predetermined threshold, which complicates assessment of summary statistics [11]. However, in practice, researchers will have an implicit or explicit maximum allowable pvalue for a variant to proceed to replication. If, for example, resources were available to followup 100 variants, but the 100th smallest Bonferronicorrected p_{d} value was >1, the variant would not generally be followed up. It is this implicit threshold  representing the maximum allowable p_{d} value which would be be deemed acceptable  which is considered to be α. A similar implicit threshold at the replication stage is the effective value of β. If no thresholds α, β are used (that is, α=β=1), then the procedure can be considered as a standard metaanalysis of the discovery and replication studies, and cannot be improved upon by combining controls at the replication stage.
If the method proposed in this paper is to be considered in a study, the values α, β, γ should be determined by the values which would otherwise have been used in a standard replication procedure. In the context of GWAS analysis, the threshold γ=5×10^{−8} should be retained, and the values α, β should reflect the implicit maximum allowable level above. The corresponding β^{∗}/β^{⊥} values can then be determined. As in any statistical procedure, the overall falsepositive rate should be considered along with the cost of following up falsepositives.
Discussion
This paper proposes a method to improve efficiency of data use in a replication procedure, adding to the body of methods for comparison of highdimensional casecontrol studies. For many common study sizes, the method can reduce the cost of replication, or increase power of discovery. The adapted method is simple to apply, only requiring modification of a single association threshold.
A standard replication procedure (or more general comparison of casecontrol studies) with independent control datasets does not make use of the information that the unconditional expected values of variables in control datasets are, in principle, the same. Conditional on p_{d}≤α, m_{0} is biased away from m_{1} (since the effect size is biased upwards), and this bias is greatest for nonassociated variants. If the observed difference m_{1}−m_{0} is large even accounting for this bias, and the observed difference \(m_{1}^{\prime }  m_{0}^{\prime }\) is small but consistent in direction with m_{1}−m_{0}, we intuitively expect that the variant is disease associated, with the observed m_{1}−m_{0} value being larger than its unconditional expectation, and the \(m_{1}^{\prime } m_{0}^{\prime }\) value being smaller. In a standard replication procedure, the variant would be declared null on the basis of \(m_{1}^{\prime }m_{0}^{\prime }\) being small, but in the shared controls procedure, some information from the first study is allowed to propagate through to the second. A metaanalysis in which observed values of both m_{1} and m_{0} are allowed to propagate information is stronger still, but this cannot in itself detect aberrance in \(C_{1}^{\prime }\).
Correspondingly, a more stringent threshold β^{∗}/β^{⊥} is needed to account for the bias in \(m_{1}^{\prime }\) conditioning on p_{d}<α, and the differential in power between the standard replication procedure and the two proposed here relates to the tradeoff between these two effects. By considering which method has the highest power in a given circumstance, the same dataset can in theory yield more information when controls are shared, while retaining some of the systematic errordetecting properties of the standard replication procedure.
The most important caveat of these methods is the loss of systematic type1 error rate control for null SNPs which are aberrant in C_{0}. Control of such errors must not be sacrificed entirely, but in some circumstances it may be satisfactory to assess such errors on a SNPbySNP basis. Such assessment is important and standard for all proposed GWAS hits under any method [16] in the interests of quality control. In method C, control over aberrance in \(C_{0}^{\prime }\) is additionally lost; however, since this method is largely applicable when \(C_{0} \cup C_{0}^{\prime }\) is a single homogeneous control (or case) cohort, there is no way that aberrance in the cohort can be systematically identified by comparison with other cohorts.
Somewhat better control of the type1 error rate can often be achieved for SNPs with aberrance in C_{1} or \(C_{0}^{\prime }\). This may incentivise the use of this method when confidence in the representativeness of these cohorts is low compared to that of C_{0}. The type 1 error rate is somewhat increased for SNPs with aberrance in \(C_{1}^{\prime }\), although as it remains bounded by α, this increase is not a major problem.
The twostage validation procedure is similar to a metaanalysis of the discovery and validation experiments, for which several adaptations to sharedcontrol designs have been proposed [3, 4]. However, there are several important distinctions which necessitate an alternative approach in this case. Firstly, not all variables are measured in the second (replication) study; we are restricted to analysis of variables reaching a given observed effect size. Secondly, the studies to be ‘metaanalysed’ are not complete, in the sense that there may be residual confounding; a strong effect size in the metaanalysis alone is not adequate evidence for association and some level of association (with consistent direction) is additionally required in both constituent studies.
The method is inapplicable when replication is performed on cohorts from completely distinct geographic groups, although there can be some difference in geographic distribution between control sets if this is controlled for in computing summary statistics. The method is most applicable when control groups are sampled from similar populations and genotyped on similar platforms. The method proposed in this paper is not universally applicable, and may only yield a modest increase in power, at the cost of changing sensitivity to different types of errors. However, it is in the interest of all researchers to use data as efficiently as possible, and methods such as this which may provide improvements without additional cost in resources should be considered as analytical options.
The widespread discoveries of the GWAS field have led to corresponding increases in complexity of phenotypic definitions, with everfiner delineations of disease types of everrarer prevalence. The genetic analysis of such complex phenotypes is necessarily comparative; there is little use understanding the genetics of a rare disease subtype except in the context of the genetics of the disease in general. Such analyses necessitate GWAS and other comparative studies between rare phenotypic types [17], with ‘controls’ meaning the bettercharacterised disease subphenotype in this sense, as well as between cases and controls. Rare disease subtypes are often afflicted with ascertainment difficulties, leading to varying degrees of expected aberrance in disease cohorts. Within this paradigm, the applicability of this method is likely to expand.
Conclusions
This paper details a method in which controls are shared in the replication phase of a twostage association study. Sharing controls can improve the power of the twostage procedure at a fixed type1 error rate. The action of sharing controls changes the spectrum of sensitivity to systematic errors caused by confounders affecting one of the study cohorts, and this should be accounted for if the sharedcontrol design is used. Adaptations of the method can enable a partial replication to be performed with only a new control cohort, or to enable robustness to misascertainment of control samples in the replication cohort.
Methods
Definitions
Denote z_{d}, z_{r}, z_{s}, z_{m}, z_{c} as the signed zscores corresponding to p_{d}, p_{r}, p_{s}, p_{m} and p_{c} respectively (where subscripts d, r, s, m, c are as defined in the “Results” section), so z_{d}=±Φ^{−1}(p_{d}/2) and so on (where Φ, Φ^{−1} are the standard normal CDF and quantile functions). Define \(\phantom {\dot {i}\!}z_{\alpha }, z_{\beta }, z_{\beta ^{*}}, z_{\beta ^{\perp }}, z_{\gamma }\) as the positive corresponding thresholds for α, β, β^{∗}, β^{⊥}, γ respectively, so z_{α}=−Φ^{−1}(x/2) etc. Other than (z_{d},z_{r}), all pairs of zscores are correlated under \(H_{0}^{=}\), with correlation estimable from sample sizes or empirically if covariates are used (Additional file 1: Appendix 1). Denote ρ_{ij} as the correlation between z_{i} and z_{j}, (i,j)∈{d,r,s,m,c}^{2}, and set
For i∈{d,r,s,m,c} define ζ_{i}=E(z_{i}), where the expectation is conditional on the SNP in question. For SNPs in \(H_{0}^{=}\), ζ_{i}≡0 for all i, but this may not hold for SNPs in \(H_{0}^{\cup } \setminus H_{0}^{=}\). In theoretical working, aberrance in groups is characterised by values ζ rather than logodds ratios. Define R_{A}, R_{B}, R_{C} as the falsepositive rates for a SNP of interest in methods A, B and C respectively.
General type 1 error rate
The values β^{∗}, β^{⊥} are chosen to satisfy
thus conserving the type 1 error rate (denoted P_{0}) against \(H_{0}^{=}\) between methods (Fig. 5). If no threshold is used on p_{m} (ie, γ=1), then β^{∗}, β^{⊥} satisfy
since \(z_{d} \perp \perp z_{r}H_{0}^{=}\). Definition Eq. (4) will be considered a generalisation of definition Eq. (5), with results established first for β^{∗} as per definition Eq. (5) and extending where possible to definition Eq. (4).
For β^{∗} defined as per definition Eq. (5) we have (see Additional file 1: Appendix 2)
approaching from above, so \(z_{\beta ^{*}} > \max \left (z_{\beta }, \sqrt {1\rho _{ds}^{2}}z_{\beta }\right. \left. + \rho _{ds} z_{\alpha }{\vphantom {\sqrt {1\rho _{ds}^{2}}z_{\beta }}} \right)\) and \(z_{\beta ^{\perp }} > \max \left (z_{\beta },\sqrt {1\rho _{cs}^{2}}z_{\beta } + \rho _{cs} z_{\alpha } \right)\). As defined by Eq. 5, \(\phantom {\dot {i}\!}z_{\beta ^{*}}, z_{\beta ^{\perp }}\) are also asymptotically linear in z_{α}, z_{γ}, z_{β} as the former two tend to ∞, with some constraints (Additional file 1: Appendix 2.1), although the limit does not necessarily approach from above. For both definitions, β^{⊥}<β^{∗}<β (Additional file 1: Appendix 2.2).
Empirical computations
Define N_{Σ}(z) as the pdf of the multivariate normal with mean 0 and variance Σ at z. Determination of covariance is described in Additional file 1: Appendix 1. Given ζ_{d}, ζ_{r}, ζ_{s}, ζ_{m}, the probability of rejecting the null for a given SNP using method A is
and using method B
If \(\frac {n_{0}}{n_{1}}=\frac {n_{0}^{\prime }}{n_{1}^{\prime }}\), matrix Σ_{A} is singular (Additional file 1: Appendix 1), in which case z_{m}=ρ_{dm}z_{d}+ρ_{vm}z_{v} and the expression above may be reduced to a twodimensional integral over a more complex region (Fig. 5). Matrix Σ_{C} is generally singular, so the formula \(z_{m} = \frac {\rho _{cs}\rho _{sm}\rho _{cm}}{\rho _{cs}^{2}  1} z_{d} + \frac {\rho _{cs}\rho _{cm}  \rho _{sm}}{\rho _{cs}^{2}  1} z_{s}\) is used to reduce the integral in a similar way. A similar formula may be used if Σ_{B} is nearly singular.
In Fig. 3, mean power difference is determined as the integral of the power difference with respect to the logodds ratio over the real line, as discussed in the “Results” section.
Study sizes, odds ratios and allele frequencies
Consider a study with n_{0} controls and n_{1} cases, with underlying allele frequencies μ_{0} and μ_{1} in cases and observed allele frequencies m_{0}, m_{1}. Let Z be a signed Zscore derived from a GWAS pvalue against the null hypothesis μ_{0}=μ_{1}. Considering Z to be proportional to the logoddsratio divided by its standard error, we have:
Setting δ=m_{1}−m_{0}, \(\bar {m} = \frac {n_{0} m_{0} + n_{1} m_{1}}{n_{0}+ n_{1}} m_{0}=\bar {m}  k \delta \), \(m_{1}=\bar {m} + k \delta \) for some k, we have
so
where \(\bar {\mu } = \frac {n_{0}\mu _{0}+n_{1} \mu _{1}}{n_{0}+n_{1}}\). Hence
where \(\bar {\mu }\) varies between definitions (though it is taken to be approximately equal).
Estimation of covariance between Z scores
Correlation between Zscores under H_{0} can be computed analytically with the following formulas (with ρ_{dr}=0):
More general formulae are given in Additional file 1: Appendix 1.
Empirical estimation of covariance and ζ values
The above formulae allow ζ and ρ to be estimated in empirical computations. The estimates may be poor if covariates or strata are used in the computation of z_{i}. Correlation may be estimated in several ways:

1.
If strata alone are used, or covariates are adjusted for in an analogous way to strata, correlations ρ_{ij} between zscores is estimable using analytic formulas (see Additional file 1: Appendix 1).

2.
If a set of variants known to be in \(H_{0}^{=}\) is available, the sample correlation between observed zscores at these variants can be used as an estimator for values ρ_{ij},

3.
A set of genotypes can be simulated for each sample for a set of variants in \(H_{0}^{=}\). Zscores corresponding to these variants can then be computed under the same correlation structure, and the sample correlation between these Zscores.
Estimates of values ζ corresponding to given logodds ratios and minor allele frequencies can be estimated in a similar way to method 1; that is, by simulating variants with given underlying oddsratios between cases and controls, computing z scores using the same method and covariance structure as used in the main study, and setting the relevant value of ζ to the observed mean z score.
False ascertainment
In general, for a true association, \(\mu _{0}=\mu _{0}^{\prime }\) and \(\mu _{1}=\mu _{1}^{\prime }\). If some proportion κ of samples in \(C_{0}^{\prime }\) are incorrectly assigned and come from the case population, then \(\mu _{0}^{\prime } = (1\kappa)\mu _{0} + \kappa \mu _{1}\). This lowers the absolute values of ζ_{r}, ζ_{s} and ζ_{m}, reducing the probability that the z_{r} score for the SNP will reach the requisite threshold β and hence reducing the power to detect the SNP using method A. This loss of power is lowered when using methods B or C.
Type 1 error rates
Aberrance in C _{1}
For SNPs aberrant in only C_{1} we have ζ_{d}≠0, ζ_{c}≠0, ζ_{m}≠0, and ζ_{r}=ζ_{s}=0.
R_{A}, R_{B}, R_{C} can be considered as functions of ζ_{d}. As ζ_{d}→0, R_{A},R_{B},R_{C}→P_{0} (Eq. 4). As ζ_{d}→±∞, \(R_{A} \to \frac {\beta }{2}\), \(R_{B} = \frac {\beta ^{*}}{2}\) and \(R_{C} = \frac {\beta ^{\perp }}{2}\). For positive ζ_{d} both R_{A} and R_{B} are increasing (and both are symmetric in ζ_{d}) so \(R_{A}<\frac {\beta }{2}\), \(R_{B} < \frac {\beta ^{*}}{2}\), \(R_{C} < \frac {\beta ^{\perp }}{2}\) for all ζ_{d}.
Since β^{⊥}<β^{∗}<β (often substantially), methods B and C are generally better at rejecting \(H_{0}^{=}\) for such SNPs. In the simplified case where z_{γ}=1, R_{A}≥R_{B} universally (Additional file 1: Appendix 3.1). This typically holds for all z_{γ}, except for small deviations in pathological cases.
In general, we consider aberrance which is only still present after any strata or covariates have been accounted for in the computation of z scores. If strata or covariates remove the effective aberrance between groups, the type1 error rate is equivalent to that under \(H_{0}^{=}\).
Aberrance in \(C_{1}^{\prime }\)
For SNPs aberrant in \(C_{1}^{\prime }\), we have ζ_{d}=0, ζ_{c}=0, ζ_{r}≠0, ζ_{s}≠0 and ζ_{m}≠0.
Again, R_{A},R_{B},R_{C}→P_{0} as ζ_{r}→0. As ζ_{r}→±∞, \(R_{A},R_{B},R_{C} \to \frac {\alpha }{2}\), and both are bounded by \(\frac {\alpha }{2}\). Although R_{B} and R_{C} are typically higher than R_{A} in this case, since both have the same (typically conservative) upper bound, this is not typically a large sacrifice in type 1 error.
In the simplified case where γ=1, an approximate upper bound on R_{B}−R_{A} is given by (Additional file 1: Appendix 4)
where
In practice, there is typically a similarly small difference between R_{C}, R_{B} and R_{A} in the general case.
Aberrance in \(C_{0}^{\prime }\)
For SNPs aberrant in \(C_{0}^{\prime }\), ζ_{d}=0, ζ_{r}≠0, ζ_{c}≠0, ζ_{s}≠0 and ζ_{m}≠0. As for SNPs with aberrance in \(C_{1}^{\prime }\), R_{A},R_{B},R_{C}→P_{0} as ζ_{r}→0 and as ζ_{r}→±∞, \(R_{A},R_{B} \to \frac {\alpha }{2}\), both bounded above by \(\frac {\alpha }{2}\). R_{C}, however, tends to 1 as ζ_{d}→∞.
In method B the cohort C_{0} has a correcting effect on the replication study, meaning ζ_{s}<ζ_{r} and R_{B}<R_{A}.
For the simplified case where γ=1, a similar bound to 14 holds for the difference R_{A}−R_{B} (note signs are reversed) with
in the place of k. The improvement in type1 error rate for a SNP with aberrance in \(C_{0}^{\prime }\) is generally larger than the loss with the same aberrance in \(C_{1}^{\prime }\) (see methods), meaning that if aberrances are of similar prevalence and size in \(C_{1}^{\prime }\) and \(C_{0}^{\prime }\), method B will typically have a lower type1 error rate than method A.
Aberrance in C _{0}
Aberrance in C_{0} represents a serious problem in casecontrol study comparison. Falsepositive rates are generally worse under method B, and tend to 1 as E(z)→∞. If aberrances of this type are expected to be very frequent, this may preclude use of methods B or C.
However, aberrances of this type may be best detected retrospectively by examining aberrances between control groups at SNPs declared ‘hits’. This procedure is already a necessary qualitycontrol procedure in method A [12, 16], as method A does not provide any control over differences between C_{0} and \(C_{0}^{\prime }\). The number of SNPs reaching significance in the twostage procedure is usually small enough that this examination is readily tractable.
Aberrance in two or more cohorts
If SNPs are aberrant in both C_{1} and \(C_{1}^{\prime }\), or in both C_{0} and \(C_{0}^{\prime }\), the effect on R_{A} and R_{B} is similar. If both cohorts are aberrant in the same direction, there is no way to differentiate the SNP from a genuine association on the basis of the genotype data alone. If cohorts are aberrant in different directions, then in both methods, the type1 error rate is lower than for a null SNP with no aberration or aberration in only one cohort, as effect sizes for the discovery and replication cohorts are biased in opposite directions. The same typically holds if \(C_{0}^{\prime }\) and C_{1}, or C_{0} and \(C_{1}^{\prime }\), are biased in the same direction.
If \(C_{0}^{\prime }\) and \(C_{1}^{\prime }\) or C_{0} and C_{1} are both biased in the same direction, R_{A} is generally lower than R_{B}, as ζ_{s}≠0. Both R_{A} and R_{B} are bounded by \(\frac {\alpha }{2}\) in this case. In addition, a systematic bias in both replication groups (or both discovery groups) is likely to be due to a known confounder, the effect of which can be removed by performing a stratified test (as is typically good practice when confounders are known). Aberrance in opposite directions leads to R_{B}>R_{A} in the first case, and a scenario similar to aberrance in C_{0} in the second case.
Aberrance in three or more cohorts corresponds to a chaotic scenario in which neither methods A,B, or C will reliably provide FPR control. Aberrance of this extent is typically detectable and removable using quality control procedures.
Abbreviations
 GWAS:

Genomewide association study
 MAF:

Minor allele frequency
 SNP:

Singlenucleotide polymorphism
References
 1
Wason JM, Dudbridge F. A general framework for twostage analysis of genomewide association studies and its application to casecontrol studies. Am J Hum Genet. 2012; 90(5):760–73.
 2
Fuchsberger C, Flannick J, Teslovich TM, Mahajan A, Agarwala V, Gaulton KJ, Ma C, Fontanillas P, Moutsianas L, McCarthy DJ, et al.The genetic architecture of type 2 diabetes. Nature. 2016; 536(7614):41–7.
 3
Lin D, Sullivan PF. Metaanalysis of genomewide association studies with overlapping subjects. Am J Hum Genet. 2009; 85(6):862–72.
 4
Han B, Duong D, Sul JH, de Bakker PI, Eskin E, Raychaudhuri S. A general framework for metaanalyzing dependent studies with overlapping subjects in association mapping. Hum Mol Genet. 2016; 90:1857–66.
 5
Bhattacharjee S, Rajaraman P, Jacobs KB, Wheeler WA, Melin BS, Hartge P, Yeager M, Chung CC, Chanock SJ, Chatterjee N, et al. A subsetbased approach improves power and interpretation for the combined analysis of genetic association studies of heterogeneous traits. Am J Hum Genet. 2012; 90(5):821–35.
 6
Zaykin DV, Kozbur DO. Pvalue based analysis for shared controls design in genomewide association studies. Genet Epidemiol. 2010; 34(7):725–38.
 7
Liley J, Wallace C. A method for identifying genetic heterogeneity within phenotypically defined disease subgroups. PLoS Genet. 2015; 49(2):310–6.
 8
Fortune MD, Guo H, Burren O, Schofield E, Walker NM, Ban M, Sawcer SJ, Bowes J, Worthington J, Barton A, Eyre S, Todd JA, Wallace C. Statistical colocalization of genetic risk variants for related autoimmune diseases in the context of common controls. Nat Genet. 2015; 47:839–46.
 9
Chanock SJ, Manolio T, Boehnke M, Boerwinkle E, Hunter DJ, Thomas G, Hirschhorn JN, Abecasis G, Altshuler D, BaileyWilson JE, et al.Replicating genotype–phenotype associations. Nature. 2007; 447(7145):655.
 10
Garner C. Upward bias in odds ratio estimates from genomewide association studies. Genet Epidemiol. 2007; 31(4):288–95.
 11
Bowden J, Dudbridge F. Unbiased estimation of odds ratios: combining genomewide association scans with replication studies. Genet Epidemiol. 2009; 33(5):406–18.
 12
Anderson CA, Pettersson FH, Clarke GM, Cardon LR, Morris AP, Zondervan KT. Data quality control in genetic casecontrol association studies. Nat Protoc. 2010; 5(9):1564–73.
 13
Stahl EA, Raychaudhuri S, Remmers EF, Xie G, Eyre S, Thomson BP, Li Y, Kurreeman FAS, Zhernakova A, Hinks A, Guiducci C, Chen R, Alfredsson L, Amos CI. Genomewide association study metaanalysis identifies seven new rheumatoid arthritis risk loci. Nat Genet. 2010; 42(6):508–16.
 14
Ferrari R, Hernandez DG, Nalls MA, Rohrer JD, Ramasamy A, Kwok JB, DobsonStone C, Brooks WS, Schofield PR, Halliday GM, et al. Frontotemporal dementia and its subtypes: a genomewide association study. Lancet Neurol. 2014; 13(7):686–99.
 15
Sarig O, Bercovici S, Zoller L, Goldberg I, Indelman M, Nahum S, Israeli S, Sagiv N, De Morentin HM, Katz O, et al. Populationspecific association between a polymorphic variant in ST18, encoding a proapoptotic molecule, and pemphigus vulgaris. J Investig Dermatol. 2012; 132(7):1798–805.
 16
The Wellcome Trust Case Control Consortium. Genomewide association study of 14000 cases of seven common diseases and 3000 shared controls. Nature. 2007; 447(7145):661–78.
 17
Liley J, Todd JA, Wallace C. A method for identifying genetic heterogeneity within phenotypically defined disease subgroups. Nat Genet. 2016.
Acknowledgments
JL would like to thank Dr Chris Wallace and Dr Jenn Asimit for helpful comments and review of this work.
Funding
This work was mostly performed while JL was funded by the NIHR Cambridge Biomedical Research Centre and on the Wellcome Trust PhD programme in Mathematical Genomics and Medicine at the University of Cambridge. During its completion, JL was funded by the Wellcome Trust (107881). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Availability of data and materials
This work did not involve the generation of new data other than by simulation. Software to simulate data and and run methods is freely available and usable at https://wallacegroupliley.shinyapps.io/replication_shared/.
Author information
Affiliations
Contributions
The author read and approved the final manuscript.
Corresponding author
Correspondence to James Liley.
Ethics declarations
Ethics approval and consent to participate
This work analyses only simulated data and metadata (study sizes and sample collection methods) from previously published work. All data were handled in accordance with the University of Cambridge policy on the ethics of research involving human participants and personal data, available at https://www.researchintegrity.admin.cam.ac.uk/researchethics.
Consent for publication
Not applicable.
Competing interests
The author declares that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file
Additional file 1
Supplementary figures and appendices. Supplementary figures showing additional power comparisons, and appendices pertaining to the method. (PDF 941 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Casecontrol study
 Replication
 GWAS