Practical 2: A single-locus model in the mouse population
Learning objectives
Calculate genotype and allele frequencies with explicit called denominators; fit and interpret additive and dominance genotype contrasts; distinguish a homozygote contrast from an allele-substitution effect; and explain why variation associated with one marker is not total heritability.
Does body weight differ among genotypes at marker M227, and how can we describe these differences using a single-locus model?
The mouse file now includes markers M227 and M1139 alongside the phenotype and family fields used in Practical 1. We focus on M227. AA, AB and BB are teaching allele labels; they do not establish nucleotide identity or identify a causal variant. Trait units are not supplied in the file metadata, so figures use recorded units.
Load the marker data
The URL fixes a public dataset revision and the checksum verifies its contents. Downloaded data stay in R’s temporary directory.
called missing p_B q_A
1174.0000000 3.0000000 0.3517888 0.6482112
Suggested answer
There are 3 missing calls and 1174 called individuals. B has frequency 0.3518. Missing genotypes do not contribute observed chromosomes. The HWE column is an expectation under random union of gametes at the supplied frequencies, not an HWE test. Family structure and the cross design require consideration before treating these mice as an unrelated equilibrium sample.
2. Inspect body weight by genotype
Use the same complete records for both models below. Explicit exclusions make model comparisons meaningful.
boxplot(BW ~ M227, data = analysis, col ="#86B7A3",xlab ="M227 genotype", ylab ="Body weight (recorded units)")
Figure 1: Observed body weight distributions by M227 genotype. These are unadjusted phenotype comparisons.
Question 2. Describe the direction and overlap of the genotype differences. Could unequal sex or replicate composition contribute?
Suggested answer
Compare both means and within-genotype spread. Unadjusted means combine marker association with other group differences. These plots do not establish that M227 is causal, nor do they separate family and environmental effects.
3. Define additive and dominance contrasts
Write the genotype contributions relative to the midpoint of the two homozygotes as
\[g_{AA}=-a,\qquad g_{AB}=d,\qquad g_{BB}=a.\]
Thus \(a=(g_{BB}-g_{AA})/2\) is half the homozygote difference, and \(d=g_{AB}-(g_{AA}+g_{BB})/2\) is the heterozygote departure from their midpoint. A positive \(d\) does not require AB to exceed both homozygotes.
It assumes common genotype contrasts across sex and replicate groups. The heterozygote indicator must flag AB, not BB. Convert factor labels to character before named-vector lookup to avoid relying on factor level numbers.
fit_ad <-lm(BW ~ sex + reps + x + z, data = analysis)round(coef(fit_ad), 4)
(Intercept) sexMale reps2 reps3 x z
34.1144 8.4130 -0.3706 2.9062 1.6073 0.4405
Show R code
a <-unname(coef(fit_ad)["x"])d <-unname(coef(fit_ad)["z"])# A three-level factor describes the same genotype contrasts.fit_factor <-lm(BW ~ sex + reps + M227, data = analysis)stopifnot(isTRUE(all.equal(unname(fitted(fit_ad)),unname(fitted(fit_factor)), tolerance =1e-10)))
Question 3. Interpret \(a\) and \(d\). Are they the same as the raw genotype means?
Suggested answer
The fitted \(a\) is 1.6073 and \(d\) is 0.4405 in recorded body-weight units. They are genotype contrasts adjusted for sex and replicate. Raw means need not yield the same contrasts. Family dependence is not modelled here, so we use the fit descriptively rather than interpreting ordinary regression p-values as definitive genetic evidence.
Show R code
g <-c(AA =-a, AB = d, BB = a)plot(0:2, g, type ="b", pch =16, xaxt ="n", col ="#247BA0",xlab ="Genotype (B count)", ylab ="Fitted genotype contrast")axis(1, at =0:2, labels =names(g))lines(c(0, 2), c(-a, a), lty =2, col ="#B65C39")legend("topleft", c("Fitted contrasts", "Homozygote midpoint line"),col =c("#247BA0", "#B65C39"), lty =c(1, 2), bty ="n")
Figure 2: Fitted genotype contrasts relative to the homozygote midpoint, adjusted for sex and replicate.
4. Distinguish an allele-substitution effect
The average effect of substituting B for A is the slope of genotypic value on B count in a specified population. Under HWE, with the coding above,
\[\alpha=a+d(q-p).\]
It depends on allele frequency when dominance is present. The coefficient \(a\) of the additive contrast in the joint model is therefore not automatically \(\alpha\). This distinction follows the single-locus average-effect model.
Show R code
alpha_hwe <- a + d*(q-p)c(a = a, d = d, alpha_HWE_model = alpha_hwe)
a d alpha_HWE_model
1.6073358 0.4405312 1.7379192
Question 4. When is \(\alpha=a\)? What changes if we reverse the counted allele?
Suggested answer
They coincide when \(d=0\) or \(p=q=0.5\). Reversing the counted allele changes the signs of \(a\) and \(\alpha\), swaps \(p\) and \(q\), and leaves \(d\) unchanged. The predicted genotype values and variance are unchanged. These labels do not mean B is inherently favorable.
Here \(V_A\) is variance of the linear allele-count prediction, and \(V_D\) is the remaining dominance variance. These are model-based quantities for the supplied frequencies and fitted contrasts; they are not independently estimated genome-wide variance components. Single-locus variance decomposition
Our observed frequencies need not equal HWE proportions. We can compute the variance of the same fitted contrasts under the observed distribution directly:
This is a frequency-weighted population variance with weights summing to one; it differs from the \(n-1\) sample-variance convention. It is variance of fitted genotype contrasts, not the observed within-genotype phenotype variation.
Question 5. Is \(V_A/\operatorname{Var}(BW)\) total narrow-sense heritability?
Suggested answer
No. The numerator describes a single-marker HWE model using fitted effects, while the denominator describes sample phenotypes. Other loci, family/environmental effects, estimation uncertainty and population assumptions remain. An associated marker can tag variation through linkage disequilibrium rather than be causal. Estimating total heritability is the subject of the next practical.
Exercises
Reverse A and B in the model. Check that the substitution effect changes sign but \(V_A\) and \(V_D\) do not.
Hold the fitted \(a,d\) fixed and vary B frequency. Plot \(\alpha\) and \(V_A\). Why can the same genotype effects yield different variances?
Set \(d=0\). Which formulas simplify? Does additive gene action imply equal allele frequencies?
Suggested solutions
Open exercise solutions
Show R code
alpha_reversed <--a + d*(p-q)stopifnot(abs(alpha_reversed + alpha_hwe) <1e-10,abs(2*q*p*alpha_reversed^2- VA) <1e-10)p_grid <-seq(0, 1, length.out =101)alpha_grid <- a + d*(1-2*p_grid)local({ old <-par(mfrow =c(1, 2))on.exit(par(old))plot(p_grid, alpha_grid, type ="l", xlab ="B frequency",ylab ="HWE model substitution effect", col ="#247BA0")plot(p_grid, 2*p_grid*(1-p_grid)*alpha_grid^2, type ="l",xlab ="B frequency", ylab ="HWE model additive variance",col ="#247BA0")})
Show R code
# Hand-solvable example: p=q=0.5, a=2, d=1.# Genotype values -2,1,2 have mean 0.5 and variance 2.25.toy_g <-c(-2, 1, 2)toy_f <-c(0.25, 0.5, 0.25)stopifnot(abs(sum(toy_f*(toy_g-0.5)^2) -2.25) <1e-12,abs(2*0.5*0.5*2^2+ (2*0.5*0.5*1)^2-2.25) <1e-12)
Allele reversal preserves the biological predictions. Frequency changes both the amount of segregating variation and, under dominance, the average substitution effect. With \(d=0\), \(\alpha=a\), \(V_D=0\) and \(V_A=2pqa^2\); this does not require \(p=q\). At fixation genotype variance is zero. The substitution formula at fixation is interpreted as a limiting expression, not an estimable slope in a monomorphic population.
Connection to the next practical
A marker model describes only one part of quantitative variation. Next we use family and pedigree information to estimate genetic parameters, distinguishing additive genetic variance, environmental sources and total phenotypic variance. Population summaries can connect to gpop as its APIs become available; known-truth biological simulations belong in gsim; fitting belongs in estimation tools.
