library(simplePHENOTYPES)
library(dplyr)
library(tibble)
human <- book_human()
tibble(markers = nrow(human), individuals = ncol(human) - 5)| markers | individuals |
|---|---|
| 9028 | 279 |
6 Dominance and epistasis
Dominance makes heterozygotes deviate from the midpoint of the two homozygotes; epistasis makes the effect of one locus depend on the genotype at another. This chapter adds both to an additive trait on real human genotypes, checks by hand how each component is built, and explains why the realized genetic variance of a trait with additive and dominance layers on the same loci is not the sum of the two requested shares. It ends with the orthogonal genotypic model, which gives a Fisher-orthogonal partition when that is what a study needs.
6.1 Worked example
6.1.1 Why human genotypes
Dominance is defined by heterozygotes. The maize panel of the earlier chapters is inbred: fewer than half of its markers carry any heterozygote, and those that do carry very few. The HapMap individuals used here are outbred, so every common marker has heterozygotes. The genotype-input chapter explains how this dataset is read; the shared loader below repeats those steps.
6.1.2 Additive and dominance on the same loci
dominance() adds a dominance component. By default it reuses the additive layer’s markers (same_as_add = TRUE), so the trait has ten causal loci that act both additively and through dominance.
ad <- simulate_phenotype(human, seed = 21) |>
additive(prop = 0.4, n_qtn = 10) |>
dominance(prop = 0.2)
print(ad)#> <phenotype_sim> (realized · long format)
#> Genotypes: human Traits: 1 Architecture: independent Seed: 21
#> Variance partition (proportions of V_P):
#> additive 0.40 (10 QTNs, geometric)
#> dominance 0.20 (same QTNs as additive)
#> residual 0.40
#> Requested genetic share = 0.60 realized H² = 0.42
#> Additive + dominance share loci (variance-partition coding); realized
#> shares of V_P for the A+D block:
#> Trait_1: requested 0.60 (unit-variance scale), realized 0.42
#> = Var(A) 0.01 + Var(D) 0.41 + 2Cov(A,D) -0.00
#> = Var(c_A) 0.51 + Var(c_D) 0.25 + 2Cov(c_A,c_D) -0.35
#> ! realized - requested/V_P = (Var(c_A) + Var(c_D) - requested/V_P) +
#> 2Cov(c_A,c_D): the bracket is the covariance among same-type layers
#> (0 for one additive and one dominance layer); the cross term depends
#> on allele frequencies and on which allele is coded +1, and is not a
#> finite-sample effect. For a Fisher-orthogonal additive/dominance split
#> use additive(orthogonal = TRUE, a =, d =).The requested genetic share is 0.6 but the realized ratio printed by the object is well below it. The printed note explains why: the additive and dominance components of the same loci covary, and the covariance is negative here. The ad_report element stores the same accounting as a table.
ad$ad_report |>
select(trait, requested, realized, var_cA, var_cD, cov2_comp)| trait | requested | realized | var_cA | var_cD | cov2_comp |
|---|---|---|---|---|---|
| Trait_1 | 0.6 | 0.4160872 | 0.5096913 | 0.2548456 | -0.3484497 |
var_cA and var_cD are the realized variances of the scaled additive and dominance components, as shares of the phenotypic variance; cov2_comp is twice their covariance. The three add up to the realized genetic share.
6.1.3 Rebuild the genetic value by hand
qtn_table() lists each layer’s markers and effects. The additive component multiplies the genotype code by the effect; the dominance component multiplies a heterozygote indicator (one where the code is 0, zero elsewhere) by the effect. Each component is centered and scaled to its prop, then the two are added. The code below repeats those steps for every individual and compares the result with genetic_values().
qt <- qtn_table(ad)
X <- as.matrix(human[, -(1:5)])
rownames(X) <- human$snp
add_rows <- filter(qt, layer == "additive")
dom_rows <- filter(qt, layer == "dominance")
score_a <- colSums(X[add_rows$snp, ] * add_rows$effect)
score_d <- colSums((X[dom_rows$snp, ] == 0) * dom_rows$effect)
scale_to <- function(score, prop) sqrt(prop) * (score - mean(score)) / sd(score)
g_a <- scale_to(score_a, 0.4)
g_d <- scale_to(score_d, 0.2)
truth <- genetic_values(ad)[, "Trait_1"]
tibble(
id = names(truth),
package = truth,
by_hand = (g_a + g_d)[names(truth)]
) |>
head()| id | package | by_hand |
|---|---|---|
| NA19152 | 0.3203771 | 0.3203771 |
| NA19139 | -0.4414581 | -0.4414581 |
| NA18912 | 0.2215043 | 0.2215043 |
| NA19160 | 0.3546636 | 0.3546636 |
| NA07034 | -0.1420707 | -0.1420707 |
| NA07055 | -0.1325024 | -0.1325024 |
tibble(
quantity = c("Var(additive component)", "Var(dominance component)",
"2 Cov(additive, dominance)", "Var(genetic value)"),
parameter = c(0.4, 0.2, NA, 0.6),
simulated = c(var(g_a), var(g_d), 2 * cov(g_a, g_d), var(g_a + g_d))
)| quantity | parameter | simulated |
|---|---|---|
| Var(additive component) | 0.4 | 0.4000000 |
| Var(dominance component) | 0.2 | 0.2000000 |
| 2 Cov(additive, dominance) | NA | -0.2734594 |
| Var(genetic value) | 0.6 | 0.3265406 |
The two columns agree to rounding. The second table shows that each component has exactly its requested variance and that the covariance term is what moves the genetic variance away from 0.6. The sign of that term is explained in the theory section.
6.1.4 Dominance on separate loci, with chosen effects
same_as_add = FALSE draws new loci for the dominance layer, and effect sets the dominance effects. Here the five dominance loci have no additive effect, so the two components covary only by chance.
ad_separate <- simulate_phenotype(human, seed = 21) |>
additive(prop = 0.4, n_qtn = 10) |>
dominance(prop = 0.2, same_as_add = FALSE, n_qtn = 5,
effect = c(0.5, 0.4, 0.3, 0.2, 0.1))
print(ad_separate)#> <phenotype_sim> (realized · long format)
#> Genotypes: human Traits: 1 Architecture: independent Seed: 21
#> Variance partition (proportions of V_P):
#> additive 0.40 (10 QTNs, geometric)
#> dominance 0.20 (5 QTNs)
#> residual 0.40
#> Requested genetic share = 0.60 realized H² = 0.54qtn_table(ad_separate) |>
filter(layer == "dominance") |>
select(snp, maf, effect, var_explained)| snp | maf | effect | var_explained |
|---|---|---|---|
| rs2472695 | 0.1863799 | 0.5 | 0.0672040 |
| rs11967821 | 0.2974910 | 0.4 | 0.0598709 |
| rs1470339 | 0.3297491 | 0.3 | 0.0337934 |
| rs10108962 | 0.1164875 | 0.2 | 0.0088068 |
| rs11927298 | 0.0537634 | 0.1 | 0.0013730 |
6.1.5 Epistasis: pairs of interacting loci
epistasis() adds interacting sets of markers. n_pairs is the number of sets, interaction the number of markers per set, and interaction_type says whether each position enters through its genotype code ("a") or its heterozygote indicator ("d"). The default is pairwise additive-by-additive.
ep <- simulate_phenotype(human, seed = 22) |>
additive(prop = 0.3, n_qtn = 10) |>
epistasis(prop = 0.2, n_pairs = 3)
print(ep)#> <phenotype_sim> (realized · long format)
#> Genotypes: human Traits: 1 Architecture: independent Seed: 22
#> Variance partition (proportions of V_P):
#> additive 0.30 (10 QTNs, geometric)
#> epistasis 0.20 (3 pairs, 2-way)
#> residual 0.50
#> Requested genetic share = 0.50 realized H² = 0.52qtn_table(ep) |>
filter(layer == "epistasis") |>
select(set, snp, maf, effect)| set | snp | maf | effect |
|---|---|---|---|
| 1 | rs7530799 | 0.0734767 | 0.500 |
| 1 | rs16840608 | 0.0627240 | 0.500 |
| 2 | rs2589559 | 0.2813620 | 0.250 |
| 2 | rs12025436 | 0.1648746 | 0.250 |
| 3 | rs11883855 | 0.0215054 | 0.125 |
| 3 | rs303008 | 0.1935484 | 0.125 |
Each set has one effect, shared by its members, and var_explained is not reported for epistatic members because a single marker has no marginal share in a product term. The product is formed from centered codes, as the next chunk verifies by rebuilding the epistatic component.
epi_rows <- qtn_table(ep) |> filter(layer == "epistasis")
X_centered <- X - rowMeans(X)
set_products <- lapply(split(epi_rows, epi_rows$set), function(rows) {
rows$effect[1] * apply(X_centered[rows$snp, , drop = FALSE], 2, prod)
})
score_e <- Reduce(`+`, set_products)
add_rows <- qtn_table(ep) |> filter(layer == "additive")
score_a <- colSums(X[add_rows$snp, ] * add_rows$effect)
rebuilt <- scale_to(score_a, 0.3) + scale_to(score_e, 0.2)
truth <- genetic_values(ep)[, "Trait_1"]
all.equal(unname(truth), unname(rebuilt[names(truth)]))#> [1] TRUEHigher-order interactions and dominance-type positions follow the same construction. The next simulation stacks an additive-by-dominance layer and a dominance-by-dominance layer, and a three-way additive interaction.
ep_mixed <- simulate_phenotype(human, seed = 22) |>
additive(prop = 0.3, n_qtn = 10) |>
epistasis(prop = 0.1, n_pairs = 2, interaction_type = c("a", "d")) |>
epistasis(prop = 0.1, n_pairs = 2, interaction_type = c("d", "d"))
print(ep_mixed)#> <phenotype_sim> (realized · long format)
#> Genotypes: human Traits: 1 Architecture: independent Seed: 22
#> Variance partition (proportions of V_P):
#> additive 0.30 (10 QTNs, geometric)
#> epistasis 0.10 (2 pairs, 2-way)
#> epistasis 0.10 (2 pairs, 2-way)
#> residual 0.50
#> Requested genetic share = 0.50 realized H² = 0.46ep_three <- simulate_phenotype(human, seed = 22) |>
additive(prop = 0.3, n_qtn = 4) |>
epistasis(prop = 0.2, n_pairs = 2, interaction = 3)
qtn_table(ep_three) |>
filter(layer == "epistasis") |>
select(set, snp, effect)| set | snp | effect |
|---|---|---|
| 1 | rs7530799 | 0.50 |
| 1 | rs16840608 | 0.50 |
| 1 | rs2589559 | 0.50 |
| 2 | rs12025436 | 0.25 |
| 2 | rs11883855 | 0.25 |
| 2 | rs303008 | 0.25 |
6.1.6 The orthogonal genotypic model
In the layers above, the additive and dominance shares are set directly and the components are not orthogonal. additive(orthogonal = TRUE, a =, d =) instead assigns each locus a genotypic value of \(-a\), \(+d\), \(+a\) for zero, one, and two copies of the counted allele, scales the whole genotypic value to prop, and lets the additive and dominance variances emerge from \(a\), \(d\), and the allele frequencies.
orth <- simulate_phenotype(human, h2 = 0.6, seed = 23) |>
additive(orthogonal = TRUE, a = rep(0.5, 10), d = rep(0.5, 10), n_qtn = 10)
print(orth)#> <phenotype_sim> (realized · long format)
#> Genotypes: human Traits: 1 Architecture: independent Seed: 23
#> Variance partition (proportions of V_P):
#> additive 0.19 (orthogonal a/d model: 10 QTNs)
#> dominance 0.31 (emergent)
#> add_dom_cov 0.11 (emergent covariance)
#> residual 0.40
#> Requested genetic share = 0.60 realized H² = 0.58qtn_table(orth) |>
select(snp, maf, effect, d, var_explained) |>
head(4)| snp | maf | effect | d | var_explained |
|---|---|---|---|---|
| rs3812719 | 0.3405018 | 0.5 | 0.5 | 0.1526351 |
| rs7284020 | 0.1003584 | 0.5 | 0.5 | 0.0244860 |
| rs12920707 | 0.1738351 | 0.5 | 0.5 | 0.0434313 |
| rs2949626 | 0.2132616 | 0.5 | 0.5 | 0.0839672 |
The budget now shows an additive row, a dominance row, and an add_dom_cov row. Vectors are passed for a and d because a single number for a is the base of a geometric series, which would give each locus a different degree of dominance. The degree of dominance is d / abs(a), one here, which is complete dominance at every locus. The next table compares no dominance, complete dominance, and overdominance on the same loci.
degree_budget <- function(d) {
budget <- simulate_phenotype(human, h2 = 0.6, seed = 23) |>
additive(orthogonal = TRUE, a = rep(0.5, 10), d = rep(d, 10), n_qtn = 10) |>
getElement("var_budget")
tibble(
degree = d / 0.5,
additive = round(budget$prop[budget$component == "additive"], 3),
dominance = round(budget$prop[budget$component == "dominance"], 3),
add_dom_cov = round(budget$prop[budget$component == "add_dom_cov"], 3)
)
}
bind_rows(degree_budget(0), degree_budget(0.5), degree_budget(1))| degree | additive | dominance | add_dom_cov |
|---|---|---|---|
| 0 | 0.600 | 0.000 | 0.000 |
| 1 | 0.187 | 0.305 | 0.108 |
| 2 | 0.075 | 0.564 | -0.038 |
With d = 0 the model is purely additive and the dominance row is zero. With complete dominance a substantial part of the genotypic variance is dominance variance, and with overdominance the dominance row grows further. The three rows always add up to the requested 0.6; the additive row is not simply what is left over, because the covariance row also changes. The add_dom_cov row is the realized covariance between the additive and dominance parts of the genotypic value; it is close to zero in a large random-mating sample and nonzero in a structured sample such as this multi-population panel.
6.1.7 The one-call model
model = "AD" splits h2 equally between an additive and a dominance layer on shared loci. It is a shortcut for the first simulation of this chapter, not for the orthogonal model.
one_call <- simulate_phenotype(human, h2 = 0.6, n_qtn = 10, model = "AD", seed = 21)
one_call$var_budget| trait | component | prop |
|---|---|---|
| Trait_1 | additive | 0.3 |
| Trait_1 | dominance | 0.3 |
| Trait_1 | residual | 0.4 |
6.2 Assumptions and theory
6.2.1 Three design columns
For individual \(i\) and locus \(j\), let \(x_{ij} \in \{-1, 0, 1\}\) be the genotype code and \(h_{ij} = \mathbb{1}[x_{ij} = 0]\) the heterozygote indicator. The raw components are package definitions (.component_raw() and .epi_unit_column()):
\[ c^{A}_i = \sum_j b_j\, x_{ij}, \qquad c^{D}_i = \sum_j b_j\, h_{ij}, \qquad c^{E}_i = \sum_{s} b_s \prod_{k \in s} \big(t_{ik} - \bar t_k\big), \tag{6.1}\]
where \(s\) indexes interacting sets and \(t_{ik}\) is either \(x_{ik}\) (position type "a") or \(h_{ik}\) (type "d"), centered by its sample mean before the product. Each component is centered and scaled to its prop exactly as in Section A.1, and the genetic value is the sum of the scaled components. The worked example rebuilt the first two and the third from these definitions.
6.2.3 The orthogonal genotypic model
Under orthogonal = TRUE the locus value is \(a_j x_{ij} + d_j h_{ij}\), the classical \(-a\), \(+d\), \(+a\) genotypic values (Falconer and Mackay 1996, genotypic values; page unverified). The additive part of the genetic value uses the average effect of allele substitution
\[ \alpha_j = a_j + d_j\,(1 - 2p_j), \tag{6.4}\]
(Falconer 1985, Genet. Res. 46:337-347; Falconer and Mackay 1996, page unverified; Fisher 1918), implemented in .avg_effect(), and the breeding value is \(A_i = \sum_j \alpha_j (n_{ij} - 2p_j)\) with \(n_{ij}\) the number of copies of the counted allele. The dominance deviation is the remainder \(D = g - A\). Under random mating the two are uncorrelated in expectation; the package reports the realized add_dom_cov instead of assuming it is zero. Linkage disequilibrium between loci does not break this orthogonality in a no-epistasis model, but population structure or selection can.
The expected variances under Hardy-Weinberg proportions are \(V_A = \sum_j 2p_jq_j\alpha_j^2\) and \(V_D = \sum_j (2p_jq_jd_j)^2\) (Falconer and Mackay 1996, page unverified). The worked example reports realized rather than expected values, which is why the rows change with the sample.
6.2.4 Which model to use
- Use the variance-partition layers when the study needs a fixed share of dominance or epistatic variance in a simulated sample and treats the components as simulation quantities. Read
ad_reportbefore interpreting the realized ratio. - Use the orthogonal model when the study needs breeding values, a defined degree of dominance, or an additive/dominance split that estimation methods will be compared against. The breeding value then feeds the selection and prediction chapters.
6.3 Implementation and practical limits
dominance()and dominance-type epistatic positions need heterozygotes. On the inbred maize panel, markers without heterozygotes are skipped with a warning and the remaining loci absorb the layer’s share;filter_geno(hets = "include")restricts the candidate markers beforehand, and an F2 or an outbred population avoids the problem.- The orthogonal model cannot be combined with a separate
dominance()layer, withvary_qtn, or with the correlation-controlled"pleiotropy"and"ld"architectures;qtnalone is accepted. var_explainedis reported per marker for additive and dominance layers only.- Under
architecture = "ld"a dominance layer must reuse the linked loci or supply disjoint linked pairs, andepistasis()is not available. - The realized ratio printed by the object is a sample diagnostic. With shared loci it can differ from the requested share by a structural amount that depends on which allele is coded
1; the genotype-input chapter shows how to control that orientation.
NoteAllele orientation matters here
The sign of the additive-dominance covariance follows \(1 - 2p\) for the counted allele. Reading the same data with the minor allele coded 1 flips the sign at every locus and can move the realized ratio to the other side of the requested share. Record the orientation with your simulation.