library(simplePHENOTYPES)
library(dplyr)
library(tibble)9 Complex phenotype architectures
Two kinds of complexity are often needed beyond additive, dominance, and epistatic mean effects. The first is a locus that changes how variable a trait is rather than its mean: a variance QTL. The second is a model that no single architecture describes, for example two traits that share some causal loci but not others. This chapter builds both, checks what the simulation produces against what was requested, and explains why a vQTL is counted as residual rather than genetic variance.
9.1 Worked example
9.1.1 Load the package and the data
A variance QTL is easiest to see when all three genotype classes are present, so the first part uses the human HapMap genotypes introduced in the binary chapter and the input chapter. The second part combines two architectures on the bundled maize panel.
source("_data.R")
human <- book_human()
data("SNP55K_maize282_maf04", package = "simplePHENOTYPES")
maize <- as_numeric(SNP55K_maize282_maf04,
to_r = TRUE, to_file = FALSE, verbose = FALSE)9.1.2 Choose a locus whose variance effect can be seen
A vQTL drawn at random may land on a marker with few individuals in one genotype class, which makes its effect hard to display. We therefore pick a marker on which every class has at least 60 individuals and pass it to vqtl() through qtn. The additive layer is still drawn at random.
class_counts <- human |>
select(-(1:5)) |>
apply(1, function(x) table(factor(x, levels = c(-1, 0, 1)))) |>
t() |>
as_tibble() |>
mutate(snp = human$snp, .before = 1)
candidate <- class_counts |>
filter(`-1` >= 60, `0` >= 60, `1` >= 60) |>
slice(1)
candidate| snp | -1 | 0 | 1 |
|---|---|---|---|
| rs4908481 | 60 | 143 | 76 |
9.1.3 A variance QTL on its own locus
We request an additive share of 0.3 on ten markers and a vQTL share of 0.3 on the chosen marker, leaving 0.4 for the ordinary residual.
sim <- simulate_phenotype(human, seed = 31) |>
additive(prop = 0.3, n_qtn = 10) |>
vqtl(prop = 0.3, qtn = candidate$snp)
print(sim)#> <phenotype_sim> (realized · long format)
#> Genotypes: human Traits: 1 Architecture: independent Seed: 31
#> Variance partition (proportions of V_P):
#> additive 0.30 (10 QTNs, geometric)
#> vqtl 0.30 (1 QTNs)
#> residual 0.40
#> Requested genetic share = 0.30 realized H² = 0.29
#> (vqtl is a residual-heterogeneity component and is not counted in
#> broad-sense heritability)The printout lists the vQTL as a separate row of the variance partition and states that it is not counted in the broad-sense heritability. The requested genetic share is 0.3, the additive layer alone. qtn_table() lists the vQTL locus with its layer effect and leaves its var_explained empty: a vQTL explains no variance of the mean.
qtn_table(sim) |>
filter(layer == "vqtl") |>
select(trait, layer, snp, chr, maf, effect, var_explained)| trait | layer | snp | chr | maf | effect | var_explained |
|---|---|---|---|---|---|---|
| Trait_1 | vqtl | rs4908481 | 1 | 0.4713262 | 0.5 | NA |
9.1.4 Residual variance by genotype class
Because the phenotype is genetic value plus residual, subtracting the genetic value recovers the residual. We attach each individual’s genotype at the vQTL and compare the residual variance within the three classes.
genotype_at <- function(geno, marker) {
row <- geno[geno$snp == marker, ]
tibble(id = names(row)[-(1:5)], genotype = as.numeric(row[-(1:5)]))
}
residuals_of <- function(sim, rep = 1) {
genetic <- genetic_values(sim, rep = rep) |>
as_tibble(rownames = "id") |>
rename(genetic_value = Trait_1)
phenotypes_long(sim) |>
filter(rep == !!rep) |>
left_join(genetic, by = "id") |>
mutate(residual = value - genetic_value)
}
by_class <- residuals_of(sim) |>
left_join(genotype_at(human, candidate$snp), by = "id") |>
group_by(genotype) |>
summarise(
n = n(),
mean_residual = mean(residual),
residual_variance = var(residual),
mean_genetic_value = mean(genetic_value)
)
by_class| genotype | n | mean_residual | residual_variance | mean_genetic_value |
|---|---|---|---|---|
| -1 | 60 | -0.0215405 | 0.5045203 | 0.1766729 |
| 0 | 143 | -0.0051763 | 0.6927378 | 0.0131050 |
| 1 | 76 | 0.0267453 | 1.0746859 | -0.1641367 |
The residual variance increases from one homozygote through the heterozygote to the other homozygote, while the mean residual stays near zero in every class. The mean genetic value also differs somewhat between classes. That difference does not come from the vQTL, which adds no mean effect; it comes from the sample association between this marker and the ten additive loci. In a panel drawn from four populations with different allele frequencies such associations are common. We have created a locus that affects the spread of the residual without affecting its average.
9.1.5 Compare with a model without the vQTL
The same seed and additive layer without the variance component give a reference whose residual variance does not depend on the genotype at that marker. Any differences left between classes are sampling variation of a variance computed from 60 to 143 individuals; the replicated comparison in the next section separates such noise from the modelled pattern.
plain <- simulate_phenotype(human, seed = 31) |>
additive(prop = 0.3, n_qtn = 10)
residuals_of(plain) |>
left_join(genotype_at(human, candidate$snp), by = "id") |>
group_by(genotype) |>
summarise(n = n(), residual_variance = var(residual))| genotype | n | residual_variance |
|---|---|---|
| -1 | 60 | 0.5933948 |
| 0 | 143 | 0.6720811 |
| 1 | 76 | 0.8545010 |
9.1.6 Replications and the expected pattern
One realization can show a class difference by chance. With n_reps the residual is redrawn many times while the loci and effects stay the same. The log-linear link described in the theory section predicts that the conditional variances of the three classes are in geometric progression, with the ratio between adjacent classes set by the vQTL share and by the marker’s standardized dosage.
replicated <- simulate_phenotype(human, seed = 31, n_reps = 200) |>
additive(prop = 0.3, n_qtn = 10) |>
vqtl(prop = 0.3, qtn = candidate$snp)
dosage <- genotype_at(human, candidate$snp)
loading <- dosage |>
mutate(loading = sqrt(0.3) * (genotype - mean(genotype)) / sd(genotype)) |>
distinct(genotype, loading)
class_variances <- lapply(seq_len(200), function(r) {
residuals_of(replicated, rep = r) |>
left_join(dosage, by = "id") |>
group_by(genotype) |>
summarise(variance = var(residual), .groups = "drop") |>
mutate(rep = r)
}) |>
bind_rows()
class_variances |>
group_by(genotype) |>
summarise(mean_variance = mean(variance), sd_variance = sd(variance)) |>
left_join(loading, by = "genotype") |>
mutate(relative_factor = exp(loading) / exp(loading[genotype == 0]))| genotype | mean_variance | sd_variance | loading | relative_factor |
|---|---|---|---|---|
| -1 | 0.5138244 | 0.0899570 | -0.8308094 | 0.4557784 |
| 0 | 0.6522261 | 0.0654666 | -0.0450608 | 1.0000000 |
| 1 | 0.9395092 | 0.1181963 | 0.7406877 | 2.1940487 |
The relative_factor column is the conditional variance of the heterogeneous component relative to the heterozygote class, before the component is rescaled and added to the homogeneous residual. The replicated class variances follow the same ordering, but they are not proportional to the factor, because each class’s total residual variance is the sum of the homogeneous part, which is the same in every class, and the heterogeneous part. The sign of the pattern depends on the sign of the vQTL effect, which is drawn with the effect series.
9.1.7 Fixed and random residual mode
In fixed mode the heterogeneous component is centered and scaled so that its sample variance equals the requested share. In random mode the per-individual standard deviations are used directly, with their mean squared value equal to the share, and nothing is rescaled afterwards.
random_sim <- simulate_phenotype(human, seed = 31, residual_mode = "random") |>
additive(prop = 0.3, n_qtn = 10) |>
vqtl(prop = 0.3, qtn = candidate$snp)
fixed_res <- residuals_of(sim)
random_res <- residuals_of(random_sim)
tibble(
quantity = c("Residual variance", "Phenotypic variance",
"Genetic / phenotypic variance"),
parameter = c(0.7, 1, 0.3),
fixed = c(var(fixed_res$residual), var(fixed_res$value),
var(fixed_res$genetic_value) / var(fixed_res$value)),
random = c(var(random_res$residual), var(random_res$value),
var(random_res$genetic_value) / var(random_res$value))
)| quantity | parameter | fixed | random |
|---|---|---|---|
| Residual variance | 0.7 | 0.7511617 | 0.7598039 |
| Phenotypic variance | 1.0 | 1.0491916 | 1.0586924 |
| Genetic / phenotypic variance | 0.3 | 0.2859344 | 0.2833684 |
In neither mode does the total residual variance equal 0.7 exactly. In fixed mode the homogeneous and heterogeneous parts each have the requested sample variance, but their sample covariance is free. In random mode their sample variances fluctuate as well. The realized genetic ratio is therefore a diagnostic, as in the first chapter.
9.1.8 A vQTL on the additive loci
By default vqtl() reuses the additive loci, so the markers that shift the mean also change the spread. The number of vQTL loci then follows the additive layer and n_qtn is ignored.
shared <- simulate_phenotype(human, seed = 31) |>
additive(prop = 0.3, n_qtn = 10) |>
vqtl(prop = 0.3)
print(shared)#> <phenotype_sim> (realized · long format)
#> Genotypes: human Traits: 1 Architecture: independent Seed: 31
#> Variance partition (proportions of V_P):
#> additive 0.30 (10 QTNs, geometric)
#> vqtl 0.30 (same QTNs as additive)
#> residual 0.40
#> Requested genetic share = 0.30 realized H² = 0.27
#> (vqtl is a residual-heterogeneity component and is not counted in
#> broad-sense heritability)qtn_table(shared) |>
count(layer)| layer | n |
|---|---|
| additive | 10 |
| vqtl | 10 |
Independent variance-only loci are requested with same_as_add = FALSE and a locus count, as in vqtl(prop = 0.3, same_as_add = FALSE, n_qtn = 2).
9.1.9 Several components in one model
Mean-effect layers and the variance component can be stacked. Here an additive layer, an epistatic layer, and a vQTL on the additive loci are combined on the maize panel. The maize lines are inbred, so a dominance layer is not used here; the nonadditive chapter shows how to prepare a panel for dominance.
stacked <- simulate_phenotype(maize, seed = 5) |>
additive(prop = 0.3, n_qtn = 10) |>
epistasis(prop = 0.2, n_pairs = 3) |>
vqtl(prop = 0.2)
print(stacked)#> <phenotype_sim> (realized · long format)
#> Genotypes: maize Traits: 1 Architecture: independent Seed: 5
#> Variance partition (proportions of V_P):
#> additive 0.30 (10 QTNs, geometric)
#> epistasis 0.20 (3 pairs, 2-way)
#> vqtl 0.20 (same QTNs as additive)
#> residual 0.30
#> Requested genetic share = 0.50 realized H² = 0.50
#> (vqtl is a residual-heterogeneity component and is not counted in
#> broad-sense heritability)stacked$var_budget| trait | component | prop |
|---|---|---|
| Trait_1 | additive | 0.3 |
| Trait_1 | epistasis | 0.2 |
| Trait_1 | vqtl | 0.2 |
| Trait_1 | residual | 0.3 |
The variance budget lists every component with its share. The genetic share is the sum of the mean-effect layers, 0.5; the vQTL share belongs to the residual side of the budget, and the homogeneous residual takes what remains. The one-call form simulate_phenotype(maize, h2 = 0.5, model = "AE", n_qtn = 10) builds the two mean-effect layers with equal shares; a vQTL has no one-call shortcut and is always added with vqtl().
9.1.10 Partial pleiotropy by combining architectures
Suppose two maize traits share some causal loci and also have loci of their own. No single architecture expresses this, but two models can be combined. The first uses the pleiotropy architecture with a genetic correlation target of 0.8; the second uses independent loci for each trait.
shared_loci <- simulate_phenotype(maize, architecture = "pleiotropy",
n_traits = 2, seed = 10, cor = 0.8) |>
additive(prop = 0.4, n_qtn = 10)
own_loci <- simulate_phenotype(maize, n_traits = 2, seed = 11) |>
additive(prop = 0.4, n_qtn = 10)
tibble(
model = c("shared loci", "own loci"),
genetic_correlation = c(cor(genetic_values(shared_loci))[1, 2],
cor(genetic_values(own_loci))[1, 2])
)| model | genetic_correlation |
|---|---|
| shared loci | 0.8223777 |
| own loci | 0.0626918 |
The two inputs were given different seeds so that they draw different loci. With one seed, both models would draw the same markers for the first trait; the pleiotropy model reuses them for the second trait while the independent model draws a separate set, so the inputs would no longer be partially pleiotropic. complex_phenotypes() warns about the differing seeds and records the first.
combined <- complex_phenotypes(shared_loci, own_loci, h2 = 0.5)#> Warning: Inputs were built with different seeds (10, 11); using the first input's seed
#> (10).print(combined)#> <phenotype_sim> (realized · long format)
#> Genotypes: maize Traits: 2 Architecture: complex Seed: 10
#> Variance partition (proportions of V_P):
#> combined from: pleiotropy + independent
#> genetic [0.50, 0.50]
#> residual [0.50, 0.50]
#> Requested genetic share = [0.50, 0.50] realized H² = [0.47, 0.46]The combined model has architecture complex. Its genetic value for each trait is the sum of the two inputs’ genetic values, rescaled to the requested share, and a new common residual is added; the inputs’ own residuals are discarded. The genetic correlation of the combined traits lies between the two inputs’ correlations.
combined_genetic <- genetic_values(combined)
summed <- genetic_values(shared_loci) + genetic_values(own_loci)
tibble(
quantity = c("Genetic correlation", "Genetic variance, Trait_1",
"Correlation with the summed inputs, Trait_1"),
parameter = c(NA, 0.5, 1),
simulated = c(cor(combined_genetic)[1, 2], var(combined_genetic[, 1]),
cor(combined_genetic[, 1], summed[, 1]))
)| quantity | parameter | simulated |
|---|---|---|
| Genetic correlation | NA | 0.4130415 |
| Genetic variance, Trait_1 | 0.5 | 0.5000000 |
| Correlation with the summed inputs, Trait_1 | 1.0 | 1.0000000 |
coheritability(combined)#> Trait_1 Trait_2
#> Trait_1 0.4690321 0.1928167
#> Trait_2 0.1928167 0.4646217The combined genetic value is an exact linear rescaling of the summed inputs, so the correlation with the sum is one. Each input contributes in proportion to its genetic variance; with equal shares of 0.4, the two contribute equally and the correlation is close to the average of 0.8 and the near-zero value of the independent model. Unequal shares weight the inputs differently.
9.1.12 What a combined model cannot do
The combined object is terminal. Layers cannot be added to it, its qtn_table() is empty, and it cannot be thresholded: liability_threshold() refuses the combined object, and a threshold applied to an input is discarded with that input’s phenotype, because the combined phenotype is rebuilt from the inputs’ genetic values plus a new residual. The inputs must also be built on the same genotypes with the same individuals and the same number of traits.
additive(combined, prop = 0.1)#> Error:
#> ! additive(): a complex_phenotypes() result is terminal -- its genetic value is the rescaled sum of its inputs and it has no layers of its own, so a layer added afterwards would be silently ignored. Add the layer to one of the input models and combine again.nrow(qtn_table(combined))#> [1] 09.2 Assumptions and theory
The vQTL examples use one additive layer, one vQTL layer, and a normal residual. The combined-architecture examples use additive layers only. No selection or crossing is involved.
9.2.1 The variance-QTL model
Let \(s_{iv}\) be the score of individual \(i\) at vQTL layer \(v\), the sum of its coded genotypes at the layer’s loci multiplied by the layer’s effect series, exactly as for an additive layer. The package standardizes each layer’s score and weights it by the square root of the layer’s share \(\pi_v\) to form a loading
\[ L_i = \sum_v \sqrt{\pi_v}\; \frac{s_{iv} - \bar s_v}{s_{s_v}} , \tag{9.1}\]
where the bar and \(s_{s_v}\) denote the sample mean and sample standard deviation of the scores. The heterogeneous residual component is then
\[ \log \operatorname{Var}(E_{v,i} \mid \text{genotype}) = \text{constant} + L_i, \qquad E_{v,i} = z_i \, e^{L_i / 2}, \tag{9.2}\]
with \(z_i\) a standard normal draw. In fixed mode the vector \(E_v\) is centered and scaled to sample variance \(\pi_v\); in random mode the factors \(e^{L_i/2}\) are divided by the root mean square of the factors and multiplied by \(\sqrt{\pi_v}\), so that the mean of the per-individual variance parameters equals \(\pi_v\), and nothing is rescaled afterwards. The component is added to the homogeneous residual. These are package definitions, implemented by .apply_vqtl(); the equation map records them as entries V2-grammar-19 and V2-grammar-20.
The exponential link guarantees a positive conditional variance for any loading. For a single vQTL on one marker, the loading takes three values, one per genotype class, so the conditional variances form a geometric progression across classes. This is why the replication table showed a monotone pattern. The link is a package choice; the documentation of vqtl() cites Rönnegård and Valdar (2011) and Murphy et al. (2022) for the vQTL concept and its detection and states explicitly that neither is the source of this formula.
9.2.2 Why a vQTL is residual variance
The heterogeneous component has mean zero within every genotype class in expectation, so it does not shift the trait mean and contributes nothing to the covariance between relatives through the mean. It changes the residual distribution, not the genetic value. The package therefore keeps it in the residual share of the budget: the residual share is \(1 - \sum_\ell \pi_\ell - \sum_v \pi_v\), with the sum over mean-effect layers \(\ell\) and vQTL layers \(v\) (equation map entry V2-grammar-9, implemented by .realize_phenotype()), and the realized genetic ratio uses the vQTL share only in the denominator (entry V2-grammar-14, .realized_h2()). A detection method that models the residual variance, such as a double generalized linear model, is the appropriate way to recover such a locus; a model for the mean will not find it unless the locus also carries a mean effect.
9.2.3 How the combined model is built
For trait \(t\), let \(G^{(m)}_{it}\) be the genetic value of individual \(i\) in input model \(m\). complex_phenotypes() forms the sum \(\sum_m G^{(m)}_{it}\), centers it, divides by its sample standard deviation, and multiplies by \(\sqrt{h^2_t}\), so the combined genetic value has sample variance \(h^2_t\). A fresh residual with share \(1 - h^2_t\) is drawn under the chosen residual mode and divided by \(\sqrt{\text{reps}}\) when entry means are requested. This is a package definition, implemented in complex_phenotypes() itself.
Two algebraic consequences follow. First, the combined genetic value is a linear function of the summed inputs, which the correlation of one above confirmed. Second, if the inputs have genetic variances \(v_1\) and \(v_2\) and genetic correlations \(r_1\) and \(r_2\) between the two traits, each input gives both traits the same genetic variance, and the inputs are uncorrelated with each other, the genetic correlation of the sum is the variance-weighted average
\[ r = \frac{v_1 r_1 + v_2 r_2}{v_1 + v_2} . \tag{9.3}\]
The inputs are not exactly uncorrelated in a finite sample, and the pleiotropy input’s realized correlation is itself a sample quantity, so the combined correlation is near this value rather than equal to it. With unequal per-trait shares the general form is \(r = \sum_m r_m \sqrt{v_{m1} v_{m2}} \big/ \sqrt{\sum_m v_{m1} \sum_m v_{m2}}\), where \(v_{mt}\) is input \(m\)’s genetic variance for trait \(t\). Both forms assume the inputs’ genetic values are uncorrelated, which linkage disequilibrium between their causal loci would violate.
9.3 Implementation and practical limits
All of the computations in this chapter are R functions: the vQTL loading and heterogeneous residual are built in .apply_vqtl(), the residual budget in .realize_phenotype(), and the combination in complex_phenotypes().
- A vQTL on a marker without variation, or with a constant combined loading, is an error; on a nearly fixed marker it is legal but uninformative, which is why this chapter chose a balanced locus.
- Under the pleiotropy and linkage architectures,
vqtl(qtn = ...)must respect the architecture’s rules for shared or paired loci. - The vQTL share is a marginal phenotypic-variance share. It is not a heritability and does not add to one; requesting
propfor it is mandatory even when the simulation has anh2budget. - A combined model keeps no per-locus state: breeding-value selection and marker-based methods that need effects refuse it, and
mediation_split()returns nothing for it. Apply such methods to the inputs, or build a single model instead. - Inputs built with different seeds are combined under the first seed with a warning; the warning is a reminder that the combined object reproduces only if all inputs are rebuilt identically.