This vignette re-creates every example from the original (v1) README
/ vignette using the v2 grammar, so a v1 user can see the direct
translation. The legacy create_phenotypes() is
still shipped (frozen, bug-fix only) if you need
bit-for-bit continuity with old scripts; everything below shows the
new way.
Two structural differences to keep in mind:
h2 and effect sizes; in v2 each layer takes a
prop (share of phenotypic variance) and the trait’s h² is
the sum of those proportions. You can still pass a target
h2 to simulate_phenotype() for the one-call
shortcut.big_add_QTN_effect is gone from the
grammar (it lives only in the frozen
create_phenotypes()). Give a large-effect QTN an explicit
effect instead.v1
create_phenotypes(
geno_obj = SNP55K_maize282_maf04,
add_QTN_num = 3,
add_effect = 0.2,
big_add_QTN_effect = 0.9,
rep = 10,
h2 = 0.7,
model = "A",
home_dir = tempdir())v2 — an additive layer taking 70% of the variance
over three QTNs. The one-call form (h2 +
n_qtn) is the closest analogue:
ph <- simulate_phenotype(geno, h2 = 0.7, n_qtn = 3, n_reps = 10, seed = 1)
ph
# <phenotype_sim> (realized · long format)
# Genotypes: geno Traits: 1 Architecture: independent Seed: 1
# Variance partition (proportions of V_P):
# additive 0.70 (3 QTNs, geometric)
# residual 0.30
# Requested genetic share = 0.70 realized H² = 0.68To reproduce the v1 “one big QTN plus a geometric tail”, give the
effects explicitly instead of big_add_QTN_effect:
simulate_phenotype(geno, seed = 1) |>
additive(prop = 0.7, n_qtn = 3, effect = c(0.9, 0.2, 0.2^2))
# <phenotype_sim> (realized · long format)
# Genotypes: geno Traits: 1 Architecture: independent Seed: 1
# Variance partition (proportions of V_P):
# additive 0.70 (3 QTNs, geometric)
# residual 0.30
# Requested genetic share = 0.70 realized H² = 0.61In v1, three pleiotropic traits shared QTNs and the genetic correlation was an uncontrolled by-product of the per-trait effect sizes. v2 lets you set the genetic correlation directly (the PleioArch engine), which is the main reason to prefer the grammar here.
v1
create_phenotypes(
geno_obj = SNP55K_maize282_maf04,
add_QTN_num = 3, dom_QTN_num = 4,
big_add_QTN_effect = c(0.3, 0.3, 0.3),
h2 = c(0.2, 0.4, 0.4),
add_effect = c(0.04, 0.2, 0.1),
dom_effect = c(0.04, 0.2, 0.1),
ntraits = 3, rep = 10, architecture = "pleiotropic",
model = "AD", seed = 10, to_r = TRUE)v2 — additive + dominance on the same shared QTNs,
three traits. Without cor the correlation is emergent (as
in v1); with cor it is realized on target:
ph <- simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 3,
cor = 0.5, seed = 11) |>
additive(prop = 0.3, n_qtn = 3) |>
dominance(prop = 0.1, same_as_add = TRUE)
# Warning: dominance(): some (but not all) selected loci have no heterozygous
# individuals (or are all heterozygous), so those loci contribute nothing and the
# remaining het-bearing loci absorb the layer's `prop` -- the realized dominance
# rests on fewer loci than requested. Choose het-bearing loci (qtn =) or
# pre-filter with filter_geno(hets = "include") if that is not intended.
round(cor(genetic_values(ph)), 2)
# Trait_1 Trait_2 Trait_3
# Trait_1 1.00 0.22 0.54
# Trait_2 0.22 1.00 0.18
# Trait_3 0.54 0.18 1.00cor also accepts a full matrix for asymmetric
targets:
R <- matrix(c(1, 0.3, -0.5,
0.3, 1, -0.2,
-0.5, -0.2, 1), 3)
simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 3,
cor = R, seed = 10) |>
additive(prop = 0.4, n_qtn = 100) |>
genetic_values() |>
cor() |>
round(2)
# Trait_1 Trait_2 Trait_3
# Trait_1 1.00 0.16 -0.54
# Trait_2 0.16 1.00 -0.17
# Trait_3 -0.54 -0.17 1.00v1 had a dedicated architecture = "partially".
v2 builds it by combining a pleiotropic model (shared loci) with
an independent model (trait-specific loci) via
complex_phenotypes() — which is exactly what
“partial pleiotropy” means: some loci are shared, some are
trait-specific.
v1
cor_matrix <- matrix(c(1, 0.3, -0.9,
0.3, 1, -0.5,
-0.9, -0.5, 1), 3)
create_phenotypes(
geno_obj = SNP55K_maize282_maf04,
ntraits = 3, pleio_a = 3,
trait_spec_a_QTN_num = c(4, 10, 1),
h2 = c(0.2, 0.4, 0.8),
cor = cor_matrix,
architecture = "partially", model = "A", rep = 20, to_r = TRUE)v2 — the shared (pleiotropic) component carries the
correlation; the independent component adds trait-specific loci;
complex_phenotypes() merges them under a common
heritability:
cor_matrix <- matrix(c(1, 0.3, -0.9,
0.3, 1, -0.5,
-0.9, -0.5, 1), 3)
shared <- simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 3,
cor = cor_matrix, seed = 10) |>
additive(prop = 0.5, n_qtn = 3)
specific <- simulate_phenotype(geno, architecture = "independent", n_traits = 3,
seed = 10) |>
additive(prop = 0.5, n_qtn = 5)
partial <- complex_phenotypes(shared, specific, h2 = c(0.2, 0.4, 0.8))
partial
# <phenotype_sim> (realized · long format)
# Genotypes: geno Traits: 3 Architecture: complex Seed: 10
# Variance partition (proportions of V_P):
# combined from: pleiotropy + independent
# genetic [0.20, 0.40, 0.80]
# residual [0.80, 0.60, 0.20]
# Requested genetic share = [0.20, 0.40, 0.80] realized H² = [0.20, 0.39, 0.80]The realized correlation is now driven by the shared loci while each
trait still carries its own private QTNs — the definition of partial
pleiotropy. (A grammar layer uses one QTN count per call, so the v1
per-trait trait_spec_a_QTN_num = c(4, 10, 1) becomes a
single n_qtn here; if you need different trait-specific
counts, build one independent model per trait and combine them.)
v1’s architecture = "LD" with type_of_ld
maps directly onto v2’s architecture = "ld". The two traits
get distinct causal loci that sit in linkage
disequilibrium; ld_type = "direct" (the v2 default) links
the two traits’ causal SNPs to each other, "indirect" links
each to a shared hidden cause-of-LD locus.
ld_max/ld_min become
r2_max/r2_min.
v1
create_phenotypes(
geno_obj = SNP55K_maize282_maf04,
add_QTN_num = 3,
h2 = c(0.2, 0.4),
add_effect = c(0.02, 0.05),
rep = 5, seed = 200, architecture = "LD",
ld_max = 0.8, ld_min = 0.2, type_of_ld = "indirect", model = "A")v2
ph_ld <- simulate_phenotype(geno, architecture = "ld", n_traits = 2,
ld_type = "indirect", r2_max = 0.8, r2_min = 0.2,
seed = 200) |>
additive(prop = c(0.2, 0.4), n_qtn = 3)
# distinct causal loci per trait: QTN_t1 / QTN_t2 are the linked pair
qtn_table(ph_ld)[, c("trait", "snp", "QTN_t1", "QTN_t2", "ld_r2")]
# trait snp QTN_t1 QTN_t2 ld_r2
# 1 Trait_1 ss196439432 ss196439432 ss196439453 0.2529085
# 2 Trait_1 ss196445359 ss196445359 ss196445290 0.2033483
# 3 Trait_1 ss196493090 ss196493090 ss196493119 0.3019851
# 4 Trait_2 ss196439453 ss196439432 ss196439453 0.2529085
# 5 Trait_2 ss196445290 ss196445359 ss196445290 0.2033483
# 6 Trait_2 ss196493119 ss196493090 ss196493119 0.3019851v1 took a QTN_list; in v2 every layer accepts a
qtn argument (marker names or indices). The other layers
stay random.
v1
QTN_list <- list()
QTN_list$add[[1]] <- c("ss196523212")
QTN_list$dom[[1]] <- c("ss196510214", "ss196472187")
QTN_list$epi[[1]] <- c("ss196530605", "ss196475446")
create_phenotypes(
geno_obj = SNP55K_maize282_maf04,
add_QTN_num = 1, dom_QTN_num = 2, epi_QTN_num = 1, epi_interaction = 2,
h2 = 0.9, QTN_list = QTN_list, ntraits = 1, model = "ADE")The frozen v1 engine now rejects this exact call with an explicit
error: a single-trait run cannot take a QTN_list, and a
single-trait "ADE"/"DE" model is not supported
(both used to fail with cryptic messages). v1 accepts
QTN_list for ntraits >= 2; the v2 form
below has no such restriction.
v2 — pass the loci to each layer;
epistasis(qtn=) takes an n_pairs × interaction
matrix:
simulate_phenotype(geno, seed = 1) |>
additive(prop = 0.5, qtn = "ss196523212") |>
dominance(prop = 0.2, qtn = c("ss196510214", "ss196472187")) |>
epistasis(prop = 0.2,
qtn = matrix(c("ss196530605", "ss196475446"), nrow = 1))
# Warning: dominance(): some (but not all) selected loci have no heterozygous
# individuals (or are all heterozygous), so those loci contribute nothing and the
# remaining het-bearing loci absorb the layer's `prop` -- the realized dominance
# rests on fewer loci than requested. Choose het-bearing loci (qtn =) or
# pre-filter with filter_geno(hets = "include") if that is not intended.
# <phenotype_sim> (realized · long format)
# Genotypes: geno Traits: 1 Architecture: independent Seed: 1
# Variance partition (proportions of V_P):
# additive 0.50 (1 QTNs, geometric)
# dominance 0.20 (2 QTNs)
# epistasis 0.20 (1 pairs, 2-way)
# residual 0.10
# Requested genetic share = 0.90 realized H² = 0.89v1 used QTN_variance = TRUE. In v2 the per-QTN variance
is always available from qtn_table():
ph <- simulate_phenotype(geno, h2 = 0.6, seed = 1) |>
additive(prop = 0.4, n_qtn = 5) |>
dominance(prop = 0.2, same_as_add = TRUE)
# Warning: dominance(): some (but not all) selected loci have no heterozygous
# individuals (or are all heterozygous), so those loci contribute nothing and the
# remaining het-bearing loci absorb the layer's `prop` -- the realized dominance
# rests on fewer loci than requested. Choose het-bearing loci (qtn =) or
# pre-filter with filter_geno(hets = "include") if that is not intended.
qtn_table(ph)[, c("trait", "layer", "snp", "effect", "var_explained")]
# trait layer snp effect var_explained
# 1 Trait_1 additive ss196522658 0.50000 0.2684611843
# 2 Trait_1 additive ss196507928 0.25000 0.0653346862
# 3 Trait_1 additive ss196483086 0.12500 0.0167240360
# 4 Trait_1 additive ss196481140 0.06250 0.0041347281
# 5 Trait_1 additive ss196474743 0.03125 0.0010306992
# 6 Trait_1 dominance ss196522658 0.50000 0.0000000000
# 7 Trait_1 dominance ss196507928 0.25000 0.1784766028
# 8 Trait_1 dominance ss196483086 0.12500 0.0000000000
# 9 Trait_1 dominance ss196481140 0.06250 0.0018930570
# 10 Trait_1 dominance ss196474743 0.03125 0.0009431359v1 controlled output with output_format /
out_geno / to_r. In v2 the object is always
returned, and the exporters are separate:
ph <- simulate_phenotype(geno, h2 = 0.5, n_qtn = 3, n_traits = 2, seed = 1)
phenotypes_long(ph) # long (default), one row per id x trait x rep
phenotypes_wide(ph) # wide, one column per trait
write_phenotypes(ph, "pheno.txt") # to disk (long)
write_phenotypes(ph, "pheno_wide.csv", format = "wide", sep = ",")v1 read files through geno_obj / geno_file
/ geno_path. In v2, convert to the numeric object once with
as_numeric() (HapMap / VCF / GDS / PLINK / numeric) and
pass it to simulate_phenotype():
v1’s cor_res becomes resid_cor.
cor is the genetic correlation;
resid_cor is the correlation between the traits’
residuals (NULL, the default, keeps them
independent, exactly as before). It takes one value for every pair or a
full n_traits x n_traits correlation matrix, and also
exists on complex_phenotypes():
Each trait’s residual variance still equals its h2
target (the draws are mixed through the Cholesky factor of the matrix
and re-standardized per trait), so the realized heritability and the
printed variance budget do not change. The correlation is a target: the
realized sample correlation agrees with it up to sampling error of order
1/sqrt(n). A vqtl() component is drawn
independently per trait and dilutes it. Like all grammar draws it owes
v1 no bit-parity.
A few v1 knobs do not have a grammar equivalent yet and remain only
in the frozen create_phenotypes():
big_add_QTN_effect — intentionally
dropped (use an explicit effect).architecture = "partially" as a single
call — recreated with complex_phenotypes(), as shown
above.