Lecture 6 of the Molecular Genomics series. Testing every marker against a phenotype, and the two problems that creates.
We have dense markers, we know they are correlated with whatever lies near them, and we have a measure of how strong that correlation is. A genome-wide association study is the most direct use of those facts: test every marker against a phenotype and see which ones respond.
Before dense genotyping was affordable, association work was necessarily restricted to candidate genes: pick a gene with a plausible biological connection to the trait, genotype variants in it, and test. The approach has an obvious weakness. It can only find what was already suspected, and for quantitative traits influenced by many loci of small effect, biological intuition is a poor guide to where those loci are. The hit rate was low.
Genome-wide association removes the requirement to guess. Every marker on the chip is tested, so a signal can be found anywhere, including in regions with no previously known connection to the trait. The trade is that you are now performing tens or hundreds of thousands of tests at once, and that creates its own problem, dealt with below.
The candidate gene approach is looking for a needle in a haystack, having first guessed which haystack. A genome-wide scan searches all the haystacks, but with many markers and many of them correlated, it is closer to searching for a needle in a stack of needles: the difficulty is no longer finding something, it is telling the real signal from the many that arise by chance.
For each marker separately, we ask whether the genotype at that marker is associated with the phenotype. The simplest formulation regresses the phenotype on the genotype coded as an allele count.
Here \(y_i\) is the phenotype of animal \(i\), \(x_{ij} \in \{0,1,2\}\) is the number of copies of the reference allele that animal \(i\) carries at marker \(j\), and \(\beta_j\) is the effect of substituting one allele for the other. The null hypothesis is \(\beta_j = 0\), the marker has no effect; the alternative is that it does.
The output is one estimate and one \(p\)-value per marker, conventionally plotted against genomic position as \(-\log_{10}(p)\), which is the Manhattan plot. Note carefully what \(\beta_j\) represents. Except by coincidence, the marker is not causal. The estimate is the effect of a causal variant somewhere nearby, attenuated by the imperfect correlation between the marker and that variant. Its magnitude is therefore not a direct estimate of the biological effect, and it is diluted by exactly the factor \(r^{2}\) discussed in the previous lecture.
The single-marker regression above is the starting point, not the endpoint. Three broad families are used.
| Approach | Idea | Main limitation |
|---|---|---|
| Single-marker linear models | Fit each marker separately, as above, usually with a random polygenic effect added to absorb relatedness | Ignores that markers are correlated; the same signal appears at many neighbouring markers |
| Penalised multiple regression, such as lasso and elastic net | Fit all markers simultaneously with a penalty that shrinks most effects to zero | Selection among correlated markers is unstable; which marker is retained can be close to arbitrary |
| Bayesian methods | Fit all markers with an explicit prior on the distribution of effects, allowing most to be near zero and a few to be large | Computationally heavy; results depend on the prior, which is a modelling choice |
In simulation, methods that model all markers jointly generally outperform single-marker regression. On real data the differences are usually much smaller than the simulation results suggest, a pattern that recurs throughout this field and that we will meet again in Lecture 10.
The null hypothesis at each marker is that it has no effect, and we reject when the \(p\)-value falls below a threshold, conventionally \(p<0.05\). The definition of that threshold is the probability of rejecting the null when it is in fact true.
That is acceptable for one test. Perform 50,000 tests on markers that all have no effect, and by construction
of them will be significant at the 5% level. Every one is a false positive. Publishing a list of 2,500 loci associated with a trait, none of which are real, is a failure mode that GWAS had to solve before it could be trusted.
A significance threshold controls the error rate per test. When many tests are performed, the number of errors accumulates. Any threshold used genome-wide must be adjusted for how many tests are being made.
Two errors are possible, and the correction for multiple testing trades one against the other.
Making the threshold more stringent reduces false positives and increases false negatives. There is no setting that avoids both, and the right balance depends on what the result is for. A study intended to nominate a small number of loci for expensive follow-up work should be stringent, since a false positive wastes a laboratory year. A study intended to build a prediction model can tolerate false positives more readily, since a marker with no real effect contributes little.
Bonferroni correction. The simplest approach divides the desired overall error rate by the number of tests.
This controls the probability of making even one false positive across the whole study. It is where the widely quoted genome-wide significance threshold of \(5\times10^{-8}\) in human studies comes from: roughly \(0.05\) divided by a million independent tests.
Bonferroni is conservative, and knowing why matters. It assumes the tests are independent. In a GWAS they emphatically are not: markers in linkage disequilibrium test overlapping information, so the number of effective independent tests is smaller than the number of markers, often by a factor of several. Dividing by the full marker count therefore sets the bar higher than necessary and produces avoidable false negatives.
False discovery rate. An alternative is to control not the probability of any false positive but the expected proportion of false positives among the markers declared significant. Accepting that 5% of a list of 100 hits are spurious is often a more sensible target than insisting on a 95% probability that none of them are, particularly when the purpose is to identify regions for further study.
Permutation. The most defensible option is to shuffle the phenotypes among animals, breaking any genotype-phenotype relationship while preserving the LD structure exactly, rerun the whole analysis, and record the best \(p\)-value obtained. Repeating this many times builds the empirical distribution of the best \(p\)-value under the null, from which a threshold follows directly. It is correct by construction and expensive to compute.
The most serious source of false positives in GWAS is not multiple testing, for which the corrections above work. It is confounding by population structure.
The reasoning behind a GWAS assumes that animals differ at a marker for reasons unrelated to why they differ in phenotype, except through the causal variant. In real populations that assumption fails. Animals are related, belong to families, breeds and geographic groups, and those groups differ both in allele frequency at very many markers and, usually, in phenotype.
Suppose a sample contains two subpopulations. One is a European taurine breed, the other an African indigenous population, and the taurine animals have higher milk yield for reasons that include a century of selection and better management.
The two groups also differ in allele frequency at thousands of markers scattered across the genome, simply because they are different populations. Every one of those markers will therefore be associated with milk yield in the combined sample, because marker genotype predicts group membership and group membership predicts yield.
None of those markers has any causal relationship with milk production. The association is entirely a consequence of having sampled two populations and not accounted for it.
The classic demonstration of this effect used ancestral origin and body weight in laboratory mice, where the confounding could be constructed deliberately and its magnitude measured. The same mechanism operates in any structured sample, and livestock samples are always structured.
Note that these corrections are not free. Population structure is correlated with genuine adaptive differences between populations, so correcting for structure also removes some real signal. A variant that is genuinely at high frequency in one breed because it was selected there is difficult to distinguish from a variant that differs in frequency for no reason.
A GWAS identifies genomic regions whose variation is statistically associated with phenotypic variation in the sample studied. It is worth being precise about the limits of that statement.
It finds regions, not genes. The resolution is set by LD, which in livestock can extend over hundreds of kilobases containing many genes. Moving from an associated region to a causal variant requires denser data, sequence, and functional evidence, and is a separate research programme.
It works well for common variants and poorly for rare ones. Power depends on the number of copies of the variant in the sample. A variant at 1% frequency in 1,000 animals is present in about 20 copies, which is not enough to estimate an effect unless that effect is very large. Rare variants of large effect, which include most Mendelian disorders, are better found by other designs.
It explains less variance than heritability implies. For most quantitative traits, the significant hits together account for a modest fraction of the genetic variance. This is not a paradox but a consequence of the architecture: if a trait is affected by thousands of loci each contributing a fraction of a percent, then almost all of them fall below any defensible significance threshold. Recognising this is precisely what motivates genomic selection, which abandons significance testing altogether.
Setting a threshold. A GWAS uses 42,000 markers after quality control. Compute the Bonferroni threshold for an overall error rate of 5%. If linkage disequilibrium means there are effectively only 12,000 independent tests, what threshold would be appropriate, and what is the practical consequence of using the wrong one?
Bonferroni on the raw marker count: \(0.05/42{,}000 = 1.19\times10^{-6}\).
On the effective number of independent tests: \(0.05/12{,}000 = 4.17\times10^{-6}\), which is about 3.5 times less stringent.
The practical consequence of using \(1.19\times10^{-6}\) is a loss of power. Markers with \(p\)-values between \(1.19\times10^{-6}\) and \(4.17\times10^{-6}\) are declared non-significant although the correct family-wise threshold would have retained them. In a study with a handful of true signals of moderate effect, this can be the difference between a positive and a negative result.
The effective number of tests is usually estimated from the eigenvalues of the marker correlation matrix, or empirically by permutation, which is why permutation remains the reference method despite its cost.
Diagnosing structure. A GWAS for body weight in a sample combining three breeds produces 3,000 significant markers spread evenly across every chromosome, with no clear peaks. Before concluding that body weight is highly polygenic, what would you check, and what pattern would confirm your suspicion?
Check for confounding by population structure. Genuine polygenic architecture produces many weak signals, but they are not usually distributed with uniform density across all chromosomes; real GWAS signals cluster into peaks because markers near a causal variant are in LD with it and with each other.
The specific diagnostic is the quantile-quantile plot of observed against expected \(p\)-values, summarised by the genomic inflation factor \(\lambda\), the ratio of the median observed test statistic to its expectation under the null. Values near 1 are healthy; a \(\lambda\) of, say, 1.8 indicates that the whole distribution of test statistics is inflated, which is the signature of structure rather than of many true effects.
The confirming pattern is that the inflation disappears when principal components of the genomic relationship matrix, or a genomic relationship mixed model, are fitted. If 3,000 signals collapse to four peaks after correction, the original result was confounding.
An honest caveat: correcting for structure can also remove genuine between-breed differences, so a large drop in signal is not proof that the original signals were all spurious. It is evidence that they cannot be distinguished from structure with this sample.
Effect size attenuation. A causal variant has a true allele substitution effect of 10 kg of milk. The nearest genotyped marker has \(r^{2}=0.36\) with it. What effect size do you expect the GWAS to estimate at that marker, and why does this matter for interpreting published effect sizes?
The expected estimated effect at a marker in LD with a causal variant is attenuated by a factor of \(r\), not \(r^{2}\). With \(r^{2}=0.36\), \(r=0.6\), so the expected estimate is \(0.6 \times 10 = 6\) kg, assuming similar allele frequencies at the two loci.
The proportion of variance explained is attenuated by \(r^{2}\), so the marker captures 36% of the variance the causal variant explains. That is why sample size requirements scale as \(1/r^{2}\).
This matters for interpretation in two ways. Reported effect sizes from GWAS are systematically downward-biased estimates of the underlying causal effects, so quoting them as biological effect sizes overstates the precision of what is known. And the degree of attenuation differs from marker to marker, so the relative magnitudes of effects at different loci are distorted, not merely scaled. Comparing effect sizes between two associated regions without accounting for their differing LD to the underlying causal variants is not meaningful.
Missing heritability. A trait has heritability 0.35. A GWAS in 8,000 animals finds twelve genome-wide significant markers which together explain 4% of the phenotypic variance. Is the heritability estimate wrong, is the GWAS wrong, or neither?
Neither. The two numbers measure different things and the gap between them is expected.
Heritability is the proportion of phenotypic variance attributable to additive genetic effects in total, estimated from resemblance between relatives without reference to any particular locus. The GWAS figure is the variance explained by the small subset of loci whose effects were large enough to survive a genome-wide significance threshold.
If the trait is affected by thousands of loci each explaining a hundredth of a percent, then with 8,000 animals almost none of them reach \(p<10^{-6}\), while collectively they account for the full 0.35. The twelve that were detected are the extreme tail of the effect size distribution, plus possibly some that were lucky.
Three things reduce the gap: much larger samples, which move more loci above the threshold; estimating variance explained by all markers jointly rather than only the significant ones, which typically recovers a large share of the heritability; and sequence data, which captures rare and structural variants that chips miss.
The important conclusion is methodological. If most of the genetic variance sits in loci that will never be individually significant, then selecting on significant loci is the wrong strategy. Using all markers regardless of significance is the right one, and that is genomic selection.
