Agricultural Statistics & Data Science · Statistical genetics and biometrics · also Statistical Genomics, Genomic Data Analysis
Statistics for molecular-breeding data: genetic map construction, interval mapping of QTL, genome-wide association with mixed models, multiple-testing control, genomic prediction (GBLUP, ridge and Bayesian alphabet), and the analysis of high-throughput expression data.
In the exams: AICE-JRF Bioinformatics Unit 5 'Statistical Genomics' (heavy); AICE-JRF Ag Stat Unit 3 (present); AICE-JRF Genetics & Plant Breeding molecular units (present)
Statistical genomics is the statistics of DNA-level data in breeding and genetics: turning raw marker calls into clean genotype matrices, ordering markers into linkage maps, locating the genes behind quantitative traits by QTL mapping and genome-wide association, controlling the flood of false positives that comes from testing hundreds of thousands of markers, predicting breeding values from genome-wide markers, and finding the genes whose expression differs between treatments. It grew from classical quantitative genetics, which partitions trait variance into genetic and environmental parts using pedigrees, and from the arrival of cheap high-density genotyping and sequencing, which made the genome itself the data.
The subject sits at the centre of modern plant and animal breeding. Dairy cattle breeding programmes in many countries moved to genomic selection within a few years of its demonstration, halving generation intervals; maize, wheat, rice and other crop programmes, public and private, use marker-assisted selection, association panels and genomic prediction; aquaculture and forestry follow. A statistician in a breeding programme or genomics laboratory works daily with the methods in this manual.
The manual begins with markers and data quality, then linkage maps and map functions, QTL mapping from single-marker tests to composite interval mapping, linkage disequilibrium and population structure, genome-wide association with mixed models, multiple testing, genomic relationship matrices and GBLUP, genomic selection with ridge regression and Bayesian methods, and gene expression and RNA-seq analysis. A final chapter gathers the descriptive and inferential tools that run through all of them: allele frequencies, Hardy-Weinberg tests, likelihood ratios and LOD scores, permutation, shrinkage and cross-validation. It is written to the level of agricultural PhD fellowship examinations in bioinformatics and statistical genetics and for practitioners who must choose and defend an analysis. Worked examples show each computation; questions in examinations are of the same kind.
A molecular marker is a detectable DNA variant at a known location, used to track the inheritance of a chromosome segment. Markers differ in how many alleles they show, whether heterozygotes can be distinguished, how many can be assayed at once and their cost.
Restriction fragment length polymorphisms (RFLPs), the first DNA markers, detect variation in fragment lengths after restriction digestion; they are codominant but laborious. Randomly amplified polymorphic DNA (RAPD) and amplified fragment length polymorphism (AFLP) markers are PCR-based and scored as presence or absence of a band, so they are dominant: the heterozygote cannot be distinguished from the band-present homozygote, which loses information in F₂ mapping. Simple sequence repeats (SSRs, microsatellites) are tandem repeats of 1 to 6 base pairs whose length varies; they are codominant, often multi-allelic and highly informative, and were the workhorse of mapping in the 1990s and 2000s. Single nucleotide polymorphisms (SNPs) are single-base differences; each is usually biallelic and codominant and less informative alone, but they are by far the most abundant variants and can be genotyped by the thousands or millions on arrays or by sequencing. Genotyping-by-sequencing (GBS) and related reduced-representation methods sequence a fraction of the genome in many individuals at once, producing many SNPs cheaply but with more missing data.
The informativeness of a marker is measured by its expected heterozygosity (gene diversity), Hₑ = 1 − Σpᵢ², where pᵢ are allele frequencies, and by the polymorphism information content, PIC = 1 − Σpᵢ² − ΣΣ_{i<j} 2pᵢ²pⱼ². For an SSR with three alleles at frequencies 0.5, 0.3 and 0.2, Hₑ = 1 − (0.25 + 0.09 + 0.04) = 0.62 and PIC ≈ 0.548. A biallelic SNP has maximum Hₑ of 0.5 and maximum PIC of 0.375, at allele frequency 0.5. Highly informative multi-allelic markers are preferred for diversity and parentage studies; dense SNPs are preferred for association and genomic prediction.
SNP genotypes are usually coded as the number of copies of a reference (or minor) allele: 0, 1 or 2 for diploids. Some methods centre the codes (−1, 0, 1) or subtract twice the allele frequency. Polyploid crops such as potato, wheat and sugarcane need dosage-aware coding, and inbred lines are nearly all 0 or 2.
Raw calls contain errors and gaps, and quality control comes before every analysis. Standard filters remove markers with a low call rate (a high proportion of missing genotypes), markers with a low minor allele frequency (MAF), which carry little information and give unstable test statistics, and individuals with high missingness. Thresholds are chosen per study; values such as a call rate below 90 to 95% or a MAF below 1 to 5% are common, and the choice should be reported.
In outbred populations, strong departure from Hardy-Weinberg proportions at a SNP can signal genotyping error (for example heterozygotes miscalled as homozygotes), and extreme departures are filtered; in breeding populations under selection, inbreeding or structure, departures are expected and the filter must be used cautiously. Excess heterozygosity of an individual suggests contamination or mixed samples; unexpected relationships between supposedly unrelated samples, or sex inconsistencies in animals, reveal sample mix-ups. In family data, Mendelian inconsistencies (a calf genotype impossible given its parents) flag errors in genotypes or in recorded parentage, and parentage verification and discovery by SNP is now standard in livestock.
Missing genotypes are imputed from the patterns of surrounding markers, using linkage disequilibrium in a reference panel or haplotype-based hidden Markov models. Imputation also raises low-density panels to higher density: many animals are genotyped cheaply on a small panel and imputed to the density of a reference panel. Imputation accuracy, usually measured as the correlation between imputed and true genotypes in masked data, is lower for rare alleles and for individuals poorly related to the reference panel.
Genes on the same chromosome are linked: they tend to be inherited together unless a crossover between them during meiosis separates them. The recombination fraction r is the proportion of gametes that are recombinant. It ranges from 0 (complete linkage) to 0.5 (independent assortment, as for genes on different chromosomes or far apart on one chromosome).
Linkage maps are built from mapping populations derived from two parents differing at many markers. In a backcross (F₁ × one parent) or a doubled-haploid population, each individual shows directly whether a meiosis was recombinant between two markers, and r̂ = number of recombinants/total, with standard error √(r̂(1 − r̂)/n). In an F₂, both parental meioses contribute and r is estimated by maximum likelihood. Recombinant inbred lines (RILs), made by selfing F₂ plants for several generations, accumulate recombination over generations; for RILs by selfing the proportion of lines recombinant between two loci is R = 2r/(1 + 2r), which must be converted back to r. Multi-parent populations (nested association mapping, MAGIC) combine many founders to increase diversity and recombination.
Before mapping, each marker is tested for segregation distortion, a departure from the expected Mendelian ratio (1:1 in a backcross or doubled-haploid population, 1:2:1 for codominant markers in an F₂), by a chi-square goodness-of-fit test. Distortion arises from selection against gametes or zygotes, or from genotyping error, and affects estimates of r.
Evidence for linkage is expressed as a LOD score, the base-10 logarithm of the likelihood ratio of the data at the estimated r to the data under independence (r = 0.5): LOD = log₁₀[L(r̂)/L(0.5)]. For a backcross with k recombinants among n, LOD = k log₁₀ r̂ + (n − k) log₁₀(1 − r̂) + n log₁₀ 2. With 30 recombinants among 200, r̂ = 0.15 and LOD ≈ 23.5. A LOD of 3, meaning the data are 1,000 times more likely under linkage than under independence, is the traditional threshold for declaring linkage between two markers in experimental populations.
Map construction has three steps. Markers are first assigned to linkage groups, sets of markers connected by pairwise LOD scores above a threshold (and recombination fractions below a maximum), which ideally correspond to chromosomes. Within each group, markers are ordered: for a few markers all orders can be compared by likelihood, but for many, heuristic algorithms seek the order that minimises the total number of recombination events or the sum of adjacent recombination fractions, or maximises the likelihood. Finally, distances between adjacent markers are estimated and converted to map units.
Recombination fractions are not additive: for three loci in order A-B-C, r_AC is less than r_AB + r_BC, because double crossovers between A and C restore the parental combination. Map distance d, measured in Morgans (M) or centiMorgans (cM; 1 cM = 0.01 M), is defined as the expected number of crossovers per meiosis in the interval and is additive. A map function converts r to d.
Haldane's map function assumes crossovers occur independently (no interference): d = −½ ln(1 − 2r), with inverse r = ½(1 − e^(−2d)). Kosambi's map function allows for interference, the tendency of one crossover to reduce the chance of another nearby: d = ¼ ln[(1 + 2r)/(1 − 2r)]. For small r both give d ≈ r. For r = 0.2, Haldane gives 25.5 cM and Kosambi 21.2 cM; Kosambi distances are always shorter for the same r. Kosambi's function is widely used for plant maps; Haldane's suits organisms with little interference. Whichever is used must be reported, because distances differ.
Map length depends on population size, marker density and genotyping errors: each error at a marker adds apparent double crossovers and inflates distances, so error-checking of singletons and suspicious double recombinants is part of map construction. With genome sequences now available for most crop and livestock species, physical maps (in base pairs) and genetic maps (in cM) are compared; recombination per megabase varies along chromosomes, low near centromeres and high towards the ends in many species.
A quantitative trait locus (QTL) is a chromosome region containing one or more genes that affect a quantitative trait. QTL mapping in a biparental population tests for association between marker genotypes and trait values: if a marker is linked to a QTL, individuals with different marker genotypes differ in mean trait value.
Single-marker analysis tests each marker separately, by a t-test (backcross, DH or RIL: two genotype classes) or one-way ANOVA or regression (F₂: three classes). It is simple, needs no map, and is easily extended to covariates. Its weaknesses are fundamental. The difference between marker classes confounds QTL effect and distance: in a backcross, if the QTL has effect a (difference between QTL genotypes) and lies at recombination fraction r from the marker, the expected difference between marker classes is a(1 − 2r). A large QTL far from the marker and a small QTL close to it can give the same difference. And individuals missing the marker genotype are lost.
In an F₂ the three marker classes allow estimation of additive and dominance effects at a QTL (if the marker is very tightly linked). With genotype means QQ, Qq and qq, the additive effect is a = (QQ − qq)/2 and the dominance effect d = Qq − (QQ + qq)/2. With means 12, 11 and 6, a = 3 and d = 2, and the degree of dominance d/a = 0.67 indicates partial dominance. The proportion of phenotypic variance explained by a marker or QTL is its R² from the regression.
Interval mapping, introduced by Lander and Botstein in 1989, uses a pair of flanking markers to infer the probability of each QTL genotype at every position along the interval, and at each position fits a mixture model: the trait distribution is a mixture of normal distributions, one per QTL genotype, weighted by those genotype probabilities. Maximum likelihood (by the EM algorithm) estimates the QTL effects at each position, and the LOD score compares the likelihood with a QTL there against no QTL. Plotting LOD against position gives a LOD profile, whose peaks above a threshold indicate QTLs. Interval mapping separates position from effect, estimates effects less biased by distance, and handles missing marker data. Haley-Knott regression approximates the mixture by regressing the trait on the expected genotype codes at each position; it is fast and close to maximum likelihood except where markers are sparse or effects large.
Interval mapping still assumes one QTL at a time, so a second QTL on the same chromosome can create ghost peaks between the two, and QTLs elsewhere inflate the residual variance and reduce power. Composite interval mapping (CIM), developed by Zeng in 1994, adds selected markers elsewhere in the genome as cofactors in the model at each test position, absorbing the effects of other QTLs; markers within a window around the test position are excluded so they do not absorb the QTL being tested. Multiple interval mapping (MIM) fits several QTLs simultaneously, with their epistatic interactions, and selects the model by an information criterion.
Significance thresholds are set by permutation tests (Churchill and Doerge, 1994): the trait values are shuffled among individuals many times (often 1,000), the genome scan is repeated, and the 95th percentile of the maximum LOD over the genome becomes the 5% genome-wide threshold. This accounts for the number of tests, the map and the trait distribution. QTL position is given with a support interval, often the 1.5-LOD drop interval, which is commonly wide (10 to 30 cM in typical populations of a few hundred), far too wide to identify a gene without fine-mapping.
A known bias, the Beavis effect, is that in small populations the effects of detected QTLs are overestimated, because only QTLs whose estimated effects happened to be large pass the threshold. Validation in independent populations is needed before marker-assisted selection.
Linkage disequilibrium (LD) is the non-random association of alleles at different loci in a population. For two biallelic loci with alleles A/a and B/b, the coefficient of LD is D = p_AB − p_A p_B, where p_AB is the frequency of the AB haplotype and p_A, p_B the allele frequencies. D = 0 means linkage equilibrium. Because the range of D depends on allele frequencies, two standardised measures are used.
D′ = D/D_max, where D_max = min(p_A p_b, p_a p_B) if D > 0 and min(p_A p_B, p_a p_b) if D < 0. |D′| = 1 means at least one of the four haplotypes is absent, often indicating no historical recombination between the loci; D′ is inflated in small samples and for rare alleles.
r² = D²/(p_A p_a p_B p_b), the squared correlation between the allele indicators at the two loci. It is the measure most relevant to association mapping, because the power to detect a causal variant through a marker depends on the r² between them: roughly, a study using the marker needs 1/r² times as many individuals as one using the causal variant directly. For p_AB = 0.4, p_A = 0.5 and p_B = 0.6: D = 0.4 − 0.30 = 0.10; D_max = min(0.5 × 0.4, 0.5 × 0.6) = 0.20, so D′ = 0.5; r² = 0.01/(0.5 × 0.5 × 0.6 × 0.4) = 0.167.
LD decays through recombination: in a large random-mating population, D_t = (1 − c)ᵗ D₀ after t generations, where c is the recombination fraction between the loci. For c = 0.01 and D₀ = 0.2, after 50 generations D ≈ 0.121. LD is created by mutation, genetic drift in small populations, selection, admixture of populations with different allele frequencies, and population bottlenecks such as domestication. The extent of LD, often summarised as the distance at which average r² falls below a chosen value, differs greatly: it is short in outcrossing species with large historical populations (maize landraces, many forest trees) and long in self-pollinating crops (rice, wheat, soybean) and in livestock breeds with small effective population size. LD extent sets the marker density needed for association mapping and genomic prediction and the resolution with which a signal can be localised.
A sample of varieties or animals is rarely a single random-mating population. It contains subpopulations (for example indica and japonica rice, or flint and dent maize, or breeds), and within them families of related individuals. Structure and relatedness create LD between unlinked loci and correlations between genotype and phenotype that have nothing to do with causation: if one subpopulation is taller for any reason and differs in allele frequency at thousands of markers, all those markers appear associated with height.
Structure is described by principal component analysis (PCA) of the centred and scaled genotype matrix, whose leading components separate the main groups; by model-based clustering, which assigns each individual a proportion of ancestry from each of K inferred groups, with K chosen by a likelihood-based or cross-validation criterion; and by differentiation statistics such as F_ST, the proportion of total genetic variance due to differences between subpopulations. F_ST = (H_T − H_S)/H_T, where H_T is expected heterozygosity in the pooled population and H_S the average within subpopulations.
Relatedness among individuals is described by a kinship or genomic relationship matrix computed from markers (Chapter 7), which captures both family structure and finer cryptic relatedness that pedigrees miss. In association analysis, structure enters as fixed covariates (principal components or ancestry proportions, the Q matrix) and relatedness as the covariance of a random polygenic effect (the K matrix).
Structure is also a resource. The same diversity that confounds association provides alleles for breeding; knowledge of subpopulations guides core-collection sampling, heterotic grouping in hybrid breeding, and the design of training sets for genomic prediction. Designed populations (biparental, nested association mapping, MAGIC) avoid structure by construction, at the price of fewer alleles than a diversity panel.
A genome-wide association study (GWAS) tests each of many markers, spread over the genome, for association with a trait in a sample of individuals that are not a designed cross: a diversity panel of varieties or landraces, a breeding population, or animals in a breed. Because LD in such samples extends over shorter distances than in a biparental cross, GWAS can localise signals more finely, and it samples many alleles at once.
The simplest analysis regresses the phenotype on each marker in turn (a general linear model, GLM). In structured samples this gives many false positives. The Q + K mixed model of Yu and colleagues (2006) became the standard correction: y = Xβ + Sα + Qv + Zu + e, where y is the phenotype vector, β fixed effects such as environment or trial, α the marker effect being tested (S its genotype codes), v the effects of population structure (Q, ancestry proportions or principal components), u a random polygenic effect with u ~ N(0, Kσ²_g) where K is the kinship matrix, and e ~ N(0, Iσ²_e) the residual. The random effect absorbs the background genetic similarity among related individuals, so that a marker is judged by whether it explains more than relatedness predicts.
Fitting a mixed model for every marker is computationally heavy. Efficient approaches estimate the variance components σ²_g and σ²_e once, under the null model without marker effects, and reuse the ratio for all markers, or transform the model by the eigen-decomposition of K so that each test becomes a weighted regression. Multi-locus methods add significant markers as cofactors iteratively, which recovers power lost when the polygenic term absorbs part of a large effect; leaving the tested chromosome out of K (leave-one-chromosome-out) serves the same purpose.
For multi-environment data, the phenotype is often first summarised as best linear unbiased estimates or predictions (BLUEs or BLUPs) of each genotype across trials; for binary traits and counts, generalised linear mixed models are used; and for traits measured over time, random-regression or function-valued models. The marker effect estimate, its standard error, the test statistic and p-value are recorded for every marker.
Results are displayed in a Manhattan plot, −log₁₀(p) for each marker against its genomic position, with chromosomes in alternating colours; peaks of adjacent markers rising above a threshold line indicate associated regions. The quantile-quantile (Q-Q) plot of observed against expected −log₁₀(p) under the null shows whether the test statistics are calibrated: if most markers have no effect, the points should follow the diagonal except for a tail of true associations at the top. Early, general departure from the diagonal indicates inflation from uncorrected structure or relatedness.
The genomic inflation factor λ summarises calibration: the median of the observed chi-square (1 df) statistics divided by its expected median under the null, 0.455. λ close to 1 indicates good calibration; λ = 1.2 suggests 20% inflation of the test statistics. A λ well below 1 can indicate over-correction. Inflation can also arise from true polygenicity in very large samples, so λ is a diagnostic, not proof of confounding.
Power in GWAS depends on sample size, the proportion of variance explained by the causal variant, its allele frequency, the LD (r²) between causal variant and markers, and the significance threshold. Rare alleles and small effects need very large samples; most agricultural GWAS panels of a few hundred lines detect only moderately large effects. Significant marker effects in small samples are overestimated (the winner's curse, the GWAS analogue of the Beavis effect), so effect sizes should be validated.
A significant marker is not the causal gene: it is in LD with the causal variant, which may lie anywhere in the LD block. Follow-up uses fine-mapping, candidate genes annotated in the region, biparental validation, expression data and, increasingly, gene editing. Results from association panels should be reported with the population, the model (covariates and kinship), the multiple-testing threshold, λ and the Q-Q plot, so that readers can judge them.
Testing hundreds of thousands of markers, or tens of thousands of genes, at the conventional 5% level would yield thousands of false positives by chance alone: with 500,000 independent tests and no true effects, about 25,000 would be significant at p < 0.05. Multiple-testing procedures control an error rate over the whole family of tests.
The family-wise error rate (FWER) is the probability of at least one false positive among all tests. The Bonferroni correction controls it by testing each hypothesis at α/m, where m is the number of tests: for 500,000 SNPs at α = 0.05 the per-test threshold is 1 × 10⁻⁷, a −log₁₀(p) threshold of 7. Bonferroni holds whatever the dependence between tests, but is conservative when tests are correlated, as SNPs in LD are. The Šidák correction, 1 − (1 − α)^(1/m), is exact for independent tests and only slightly less strict.
Holm's step-down procedure is uniformly more powerful than Bonferroni while still controlling the FWER: order the p-values p₍₁₎ ≤ … ≤ p₍ₘ₎, and reject H₍₁₎, H₍₂₎, … as long as p₍ᵢ₎ ≤ α/(m − i + 1), stopping at the first failure.
Because SNPs in LD are not independent, an effective number of independent tests, M_eff, smaller than m, can be estimated from the eigenvalues of the marker correlation matrix or from the LD structure, and used in place of m. Permutation of phenotypes, recording the minimum p-value over the genome each time, gives an empirical genome-wide threshold that accounts for the correlation exactly, as in QTL mapping. In human genetics a threshold of 5 × 10⁻⁸ became conventional for genome-wide significance; in crops and livestock, thresholds depend on marker number and LD and should be justified for each study.
When many true effects are expected, as in gene-expression studies with thousands of differentially expressed genes, controlling the FWER is too strict: it buys near-certainty of no false positives at the cost of missing most true effects. The false discovery rate (FDR), introduced by Benjamini and Hochberg in 1995, is the expected proportion of false positives among the rejected hypotheses. Controlling the FDR at 5% means that, on average, no more than 5% of the discoveries are false.
The Benjamini-Hochberg (BH) procedure: order the m p-values p₍₁₎ ≤ … ≤ p₍ₘ₎; find the largest i such that p₍ᵢ₎ ≤ (i/m)q, where q is the target FDR; reject the hypotheses with the i smallest p-values. It controls the FDR at q for independent tests and under positive dependence, which covers most genomic applications; the Benjamini-Yekutieli variant holds under any dependence at the cost of a stricter threshold.
Example: ten genes have ordered p-values 0.001, 0.008, 0.012, 0.020, 0.041, 0.060, 0.074, 0.205, 0.512 and 0.900, with q = 0.05. The BH critical values (i/10) × 0.05 are 0.005, 0.010, 0.015, 0.020, 0.025, …. The fourth p-value, 0.020, equals its critical value 0.020; the fifth, 0.041, exceeds 0.025, and so do all later ones. So the largest i is 4, and four genes are declared significant. Bonferroni at 0.05/10 = 0.005 would declare only one.
Storey's q-value is the FDR analogue of the p-value: the minimum FDR at which a test would be declared significant. It also estimates π₀, the proportion of true nulls, from the distribution of p-values (which is flat for nulls), making FDR control less conservative when many hypotheses are false. Local FDR estimates the probability that a particular test is null given its statistic. Reports should state which error rate was controlled, at what level and by what procedure.
Classical animal and plant breeding predicts breeding values with best linear unbiased prediction (BLUP), using a numerator relationship matrix A computed from the pedigree: the expected proportion of genes shared identical by descent. The pedigree gives expectations only; full sibs share on average half their genes, but actual sharing varies around that value because of Mendelian sampling. Dense markers measure the realised sharing.
The most widely used genomic relationship matrix (GRM) is VanRaden's first method (2008). Let M be the n × m matrix of genotypes coded 0, 1, 2 (copies of one allele), pⱼ the frequency of the counted allele at marker j, and Z = M − 2P, where each column of P holds pⱼ. Then G = ZZ′/[2Σⱼ pⱼ(1 − pⱼ)]. Scaling by 2Σpq makes G analogous to A: diagonal elements average about 1 + F (F the inbreeding coefficient) and off-diagonals estimate twice the kinship relative to the base population defined by the allele frequencies used. A second version standardises each marker separately before averaging. The allele frequencies should ideally be those of the base population; using current-sample frequencies shifts the scale, which matters when G is combined with A.
Because each element of G is estimated from finite markers, G is noisy for small m and can be singular; a small value is added to the diagonal or G is blended with A (for example 0.95G + 0.05A) before inversion. G also captures relationships unknown from pedigree, detects pedigree errors, and estimates relationships in crops where pedigrees are incomplete or complex.
Genomic BLUP (GBLUP) replaces A with G in the standard mixed model: y = Xb + Zg + e, with g ~ N(0, Gσ²_g) and e ~ N(0, Iσ²_e), where b are fixed effects, g the genomic breeding values, and Z the incidence matrix linking records to individuals. Henderson's mixed model equations give the BLUP of g: [X′X X′Z; Z′X Z′Z + G⁻¹λ] [b; g] = [X′y; Z′y], with λ = σ²_e/σ²_g = (1 − h²)/h². The variance components are estimated by restricted maximum likelihood (REML). Individuals without phenotypes but with genotypes receive predictions through their genomic relationships with phenotyped individuals: this is how young animals or new lines are evaluated before they have records.
GBLUP is equivalent to ridge-regression BLUP of marker effects (SNP-BLUP): fitting y = Xb + Wa + e with each marker effect aⱼ ~ N(0, σ²_a) independently, where W is the centred genotype matrix, gives breeding values g = Wâ identical to GBLUP when σ²_g = 2Σpⱼ(1 − pⱼ)σ²_a. GBLUP works with an n × n matrix and SNP-BLUP with an m × m system, so the cheaper form is chosen according to whether there are more individuals or more markers.
Single-step GBLUP (ssGBLUP) combines genotyped and non-genotyped animals in one evaluation by using a blended relationship matrix H, whose inverse is A⁻¹ plus a correction for the genotyped animals, G⁻¹ − A₂₂⁻¹, where A₂₂ is the pedigree relationship among genotyped animals. It is now the standard in many national dairy, beef, pig and poultry evaluations, because only part of each population is genotyped.
The accuracy of genomic predictions depends on the size of the reference population, heritability, the effective number of chromosome segments (related to effective population size and genome length), and the relationship between reference and candidates. Daetwyler and colleagues gave an approximation for the expected accuracy, r = √[Nh²/(Nh² + M_e)], where N is the reference size and M_e the number of independent segments. With N = 1,000, h² = 0.5 and M_e = 1,000, r ≈ √(500/1,500) = 0.58.
Genomic selection (GS), proposed by Meuwissen, Hayes and Goddard in 2001, predicts the breeding value of an individual from all its marker genotypes simultaneously, without first testing which markers are significant. A training (reference) population of individuals with both genotypes and phenotypes (or progeny-test results) is used to estimate the effects of all markers together; the genomic estimated breeding values (GEBVs) of selection candidates, which may have genotypes only, are then the sums of their marker genotypes times the estimated effects.
The statistical problem is that there are far more markers than individuals (p ≫ n), so ordinary least squares cannot be used: the equations have no unique solution and would overfit. All GS methods shrink marker effects towards zero; they differ in their assumptions about how effects are distributed.
The breeding value of GS comes from the breeder's equation in its genetic-gain form, ΔG per year = i r σ_A/L, where i is the selection intensity, r the accuracy of selection, σ_A the additive genetic standard deviation, and L the generation interval. GS raises r for young candidates without records and, above all, shortens L: dairy bulls can be selected as calves rather than after a progeny test of about five years, and crop lines can be selected at early generations or in off-season nurseries without field phenotyping. Dairy cattle programmes in many countries roughly doubled their rate of genetic gain after adopting GS, largely through shorter generation intervals.
Accuracy is evaluated by cross-validation: the training data are split into k folds, each fold is predicted from a model trained on the others, and the correlation between predicted values and observed phenotypes (predictive ability) is computed; dividing by √h² estimates accuracy with respect to true breeding values. Validation should mimic the intended use: random cross-validation within closely related families overstates the accuracy achievable on new families or in new environments, so forward validation (predicting a later generation) or leave-one-family-out schemes are more honest.
Ridge-regression BLUP (RR-BLUP) assumes that all marker effects come from one normal distribution with a common variance, aⱼ ~ N(0, σ²_a). The estimates minimise the residual sum of squares plus λΣaⱼ², shrinking all effects equally towards zero. It is equivalent to GBLUP, works well when the trait is controlled by many genes of small effect (an 'infinitesimal' architecture), and is robust and fast.
The Bayesian alphabet, introduced by Meuwissen and colleagues, relaxes the common-variance assumption. In BayesA each marker has its own variance σ²_aⱼ, drawn from a scaled inverse chi-square prior; marginally, effects follow a t-distribution with heavier tails than the normal, so large effects are shrunk less than small ones. In BayesB, a proportion π of markers have zero effect and the remaining 1 − π have their own variances as in BayesA; this mixture suits traits with a few loci of large effect among many null markers. BayesC uses a common variance for the non-zero effects, and BayesCπ estimates π from the data rather than fixing it. The Bayesian LASSO uses a double-exponential (Laplace) prior, giving strong shrinkage of small effects. These models are fitted by Markov chain Monte Carlo (Gibbs sampling, with Metropolis-Hastings steps for BayesB), which is slower than BLUP and needs convergence checks and sensible priors.
Penalised regressions from machine learning offer non-Bayesian analogues: the LASSO (L1 penalty) sets many effects exactly to zero, and the elastic net combines L1 and L2 penalties. Non-parametric methods, reproducing kernel Hilbert space regression and machine learning models, capture epistatic and non-additive effects but rarely outperform GBLUP for additive breeding values.
In practice, for most polygenic traits (yield, milk yield, growth), GBLUP and the Bayesian methods give similar accuracies; BayesB-type methods gain where major genes exist (for example some disease resistances and quality traits). The size and relatedness of the training population and the quality of phenotypes matter far more than the choice of method, and in crops genotype-by-environment interaction must be modelled, by multi-environment genomic models with environment-specific effects or environmental covariates, for predictions to transfer to new environments.
Gene expression studies ask which genes are transcribed at different levels between conditions: drought-stressed and well-watered plants, resistant and susceptible varieties, infected and healthy tissue, animals on different diets. Microarrays measured expression by hybridising labelled RNA to probes on a chip and reading fluorescence intensities; in two-colour arrays the log-ratio M = log₂(R/G) of the two dyes and the average log intensity A = ½log₂(RG) were analysed, with MA-plots and loess normalisation to remove intensity-dependent dye bias. RNA sequencing (RNA-seq) has largely replaced arrays: RNA is converted to cDNA, sequenced, and the reads mapped to a reference genome or transcriptome and counted per gene, giving a count matrix of genes by samples.
Design principles from field experiments apply fully. Biological replicates (independent plants, animals or plots) are essential, because the inference is about the population of plants or animals, not about one sample; technical replicates (re-sequencing the same library) measure only the technical noise, which for RNA-seq is small. Three biological replicates per condition is a common minimum, and more are needed to detect small fold changes. Samples processed in different batches (dates of extraction, sequencing lanes or runs) differ systematically, so treatments must be balanced across batches, and batch included in the model; a design that confounds treatment with batch cannot be rescued by analysis. Randomisation, blocking and recording all processing variables are as important as sequencing depth, and for a fixed budget more replicates usually buy more power than deeper sequencing.
Raw counts depend on the library size (total reads per sample) and the composition of the library: if a few genes are very highly expressed in one condition, they take a larger share of reads and make all other genes look lower. Normalisation corrects for both. Counts per million (CPM) divide by library size; RPKM/FPKM and TPM also divide by gene length, which is needed to compare genes within a sample but not to compare the same gene between samples. Composition-robust methods compute a scaling factor per sample from genes assumed mostly not differentially expressed: the trimmed mean of M-values (TMM) and the median-of-ratios method are the standard choices.
RNA-seq counts are over-dispersed relative to the Poisson: between biological replicates, the variance exceeds the mean. The negative binomial (NB) distribution models this with Var(Y) = μ + φμ², where φ is the dispersion; for a gene with mean count 100 and φ = 0.1, the variance is 100 + 1,000 = 1,100, eleven times the Poisson variance. Each gene is fitted by a NB generalised linear model with log link, log μ = log(size factor) + Xβ, so that treatment effects are log fold changes. With few replicates, gene-wise dispersions are unstable, so empirical Bayes methods shrink each gene's dispersion towards a trend fitted across all genes, borrowing strength. Tests are Wald or likelihood-ratio tests of coefficients, and the p-values are adjusted by the Benjamini-Hochberg FDR.
An alternative transforms counts to log-CPM and estimates a mean-variance relation to give each observation a precision weight, then applies weighted linear models with empirical Bayes moderated t-statistics, the approach developed for microarrays.
Results are summarised by the log₂ fold change (a gene with mean normalised counts of 400 under drought and 100 under control has log₂FC = 2, a fourfold increase), its FDR-adjusted p-value, volcano plots (−log₁₀ p against log₂FC) and MA-plots. Fold-change estimates for low-count genes are noisy and are often shrunk towards zero for ranking. Lists of differentially expressed genes are interpreted with gene-set enrichment against functional annotations (gene ontology terms, pathways), which raises its own multiple-testing and background-choice issues. Expression levels can also be mapped as traits: expression QTLs (eQTLs) link genetic variants to gene expression.
Genomic data are summarised first by allele and genotype frequencies. For a biallelic locus with genotype counts n_AA, n_Aa and n_aa in N individuals, the allele frequency is p = (2n_AA + n_Aa)/2N. Under Hardy-Weinberg equilibrium (HWE), reached after one generation of random mating in a large population without selection, migration or mutation, genotype frequencies are p², 2pq and q². The chi-square test compares observed counts with Np², 2Npq and Nq², with 1 degree of freedom (three classes, minus one, minus one estimated allele frequency). For 360 AA, 480 Aa and 160 aa among 1,000 individuals, p = 0.6, expected counts are 360, 480 and 160, and χ² = 0: perfect agreement. Exact tests are preferred for rare alleles and small samples.
The inbreeding coefficient estimated from heterozygosity is F = 1 − H_obs/H_exp: a population of selfing lines has F near 1, a random-mating population near 0. In crops, observed heterozygosity is a check on line purity; in livestock, genomic inbreeding is now measured from runs of homozygosity (long homozygous stretches indicating recent common ancestry).
Diversity across a germplasm collection is summarised by the number of alleles, expected heterozygosity, allelic richness corrected for sample size, and genetic distances between accessions (for example from the proportion of shared alleles), from which dendrograms and principal coordinate plots are drawn. Effective population size Ne, the size of an ideal population that would show the same rate of drift or inbreeding, is estimated from LD or from changes in inbreeding over time; it is far smaller than census size in most livestock breeds, because few sires contribute to each generation.
Most genomic tests are likelihood-based. The likelihood ratio statistic, 2 ln(L₁/L₀), has an approximately chi-square distribution under the null with degrees of freedom equal to the number of extra parameters. The LOD score is the same comparison on the base-10 scale: LOD = log₁₀(L₁/L₀), so the likelihood ratio statistic equals 2 ln 10 × LOD ≈ 4.605 × LOD. A LOD of 3 corresponds to a chi-square of 13.8. Testing a variance component on the boundary of its parameter space (σ² = 0) gives a mixture of chi-square distributions, so the naive chi-square p-value is conservative, by a factor of about two for a single variance component.
When theoretical null distributions are doubtful, resampling helps. Permutation tests shuffle phenotypes relative to genotypes, destroying any true association while keeping the correlation among markers and the phenotype distribution, and so give exact genome-wide thresholds. The bootstrap resamples individuals with replacement to give confidence intervals for QTL positions, heritabilities or LD measures. Cross-validation estimates predictive performance on new data, guarding against overfitting in prediction models.
Genomic data are high-dimensional: often p ≫ n. Then unbiased estimators are unstable or undefined, and shrinkage (ridge, LASSO, empirical Bayes, hierarchical Bayesian priors) trades a little bias for much less variance. The same idea runs through GBLUP, the Bayesian alphabet, moderated t-statistics and dispersion shrinkage in RNA-seq.
Variance components link genomics to quantitative genetics. The SNP-based (genomic) heritability, estimated by REML with a GRM, is the proportion of phenotypic variance explained by all markers jointly; it is usually lower than pedigree-based heritability because markers do not capture all causal variation in perfect LD (the 'missing heritability'). Finally, the multiplicity of choices in genomic analysis (filters, models, thresholds) makes transparent reporting and reproducible pipelines essential: a result should be reproducible from the raw genotype and phenotype files with the stated software versions and parameters.
A doctoral student in genetics, breeding or agricultural statistics uses this subject to analyse thesis data: cleaning SNP data from a mapping population or diversity panel, building a linkage map, running interval or composite interval mapping with permutation thresholds, or a mixed-model GWAS with Q-Q plots and FDR control, and fitting GBLUP and a Bayesian model with cross-validation. Examinations ask for LOD scores, map distances, LD measures, multiple-testing thresholds and the assumptions of each method, so the student must be able to compute and to explain the choice of model.
A quantitative geneticist or biometrician five years into work at a breeding company, research institute or genomics laboratory runs routine genomic evaluations: quality control and imputation of new genotypes, parentage checks, GBLUP or single-step evaluations, cross-validation of training sets, GWAS of new traits, and RNA-seq analyses for colleagues. Such a person advises on mapping population size, number of replicates, batch balancing and which markers to put on a custom panel for marker-assisted selection, and documents every pipeline step so that colleagues can rerun it.
A senior statistician with twenty-five years of experience leads the genomic evaluation of a national livestock population or the genomic-selection strategy of a crop breeding programme: designs and maintains reference populations, chooses models and validation schemes, sets rules for publishing genomic breeding values, manages the transition from pedigree BLUP to single-step evaluation, and explains accuracy and its limits to breeders and farmers. Such a person also decides on investments in genotyping density and phenotyping, and represents the programme in international evaluation collaborations.
A professor or research leader with a PhD and long experience develops methods: models for genotype-by-environment interaction in genomic prediction, multi-omics integration, polyploid and multi-parent analyses, scalable algorithms for millions of animals, and causal fine-mapping. Such a person trains statisticians for breeding programmes, advises international and national breeding networks on genomic strategies, reviews grant and journal work, and leads efforts on open data, shared reference populations and reproducible pipelines that let smaller programmes in many countries use genomic methods well.
The manual gives the main models and formulas but not full derivations of the EM algorithm for interval mapping, the MCMC samplers of the Bayesian alphabet or the REML algorithms. Polyploid genetics, multi-parent population analysis, haplotype-based and epistatic models, multi-trait and multi-environment genomic prediction, single-cell and long-read transcriptomics, and causal fine-mapping are introduced only briefly or not at all. Quality-control thresholds, significance thresholds and typical accuracies vary between species, populations and studies; the manual gives common ranges and says they must be justified per study. Software is described by method, not by product name.
Original practice questions written to the syllabus and pattern of the exams named in each question's examTags. They are not past papers and do not reproduce any board's questions.
These notes are free to read here, with no account.
What the Laisha app adds:
Free for ten days, then a subscription. We would rather write that here than have you find it later.