The simplePHENOTYPES v2 grammar

library(simplePHENOTYPES)
data("SNP55K_maize282_maf04")
geno <- SNP55K_maize282_maf04

This vignette walks through simulating phenotypes with the simplePHENOTYPES v2 grammar, using the bundled maize 282 panel throughout. The original create_phenotypes() interface is covered at the end: it is still supported, but it is superseded, and new work should use the grammar described here.

Two companion vignettes go deeper on the topics this one only introduces:

  • vignette("breeding-populations") — crossing, selfing and doubled haploids: building a pedigree from real founders.
  • vignette("genetic-maps") — what cm is for, and how to build or replace it.

The example data

SNP55K_maize282_maf04 is 10,650 markers scored on 280 inbred lines, coded -1 / 0 / 1. The first five columns are marker metadata.

dim(geno)
# [1] 10650   285
geno[1:5, 1:8]
#           snp allele chr     pos        cm 4226 4722 33-16
# 1 ss196422159    A/G   1  379844 0.0000000   -1   -1    -1
# 2 ss196422171    G/A   1  613257 0.2482394   -1   -1    -1
# 3 ss196422173    G/A   1  659354 0.2972609    1    1    -1
# 4 ss196422186    G/C   1  992572 0.6515820    1   -1    -1
# 5 ss196500940    G/A   1 2044264 1.7694445    1    1     1

The genetic map

cm holds a genetic map in centiMorgans. Earlier versions of the package left this column empty, which meant nothing in the package could simulate meiosis.

aggregate(cm ~ chr, data = geno, FUN = max)
#    chr       cm
# 1    1 219.3573
# 2    2 172.2582
# 3    3 168.4068
# 4    4 175.4795
# 5    5 158.9096
# 6    6 123.2703
# 7    7 128.5558
# 8    8 128.1356
# 9    9 114.2056
# 10  10 108.5956

This map is synthetic. It is generated from the physical positions by synthetic_map() and is not taken from any published maize linkage map, so its cM values should not be quoted as measured recombination distances. It exists so the crossing functions have realistic distances to work with. If you have a real map for your population, use that instead — just put it in cm.

To build one for your own data:

geno$cm <- synthetic_map(geno$chr, geno$pos)

The map only matters for the crossing functions; the grammar below does not use it. See vignette("genetic-maps") for the recombination model, how to tune it, and how to validate a real map.

Reading genotype data

as_numeric() converts a genotype file into the numeric format used above, handling HapMap, VCF, GDS and PLINK bed/ped input.

num <- as_numeric("my_genotypes.hmp.txt")

It handles the other formats the same way:

num <- as_numeric("my_genotypes.vcf", to_r = TRUE)

Simulating phenotypes: the grammar

A simulation starts with simulate_phenotype(), which fixes the genetic architecture, and then adds variance components as layers. Each layer’s prop is a requested marginal proportion of total phenotypic variance. Broad-sense heritability is based on the additive, dominance, and epistatic genetic values; the vQTL layer is residual heterogeneity and is not included.

ph <- simulate_phenotype(geno, seed = 1) |>
  additive(prop = 0.5, n_qtn = 3)
ph
# <phenotype_sim>  (realized · long format)
#   Genotypes: geno   Traits: 1   Architecture: independent   Seed: 1
#   Variance partition (proportions of V_P):
#     additive    0.50   (3 QTNs, geometric)
#     residual    0.50
#   Requested genetic share = 0.50   realized H² = 0.43

The pipe runs eagerly — there is no terminal simulate() call, and every object already carries realized phenotypes.

head(phenotypes_long(ph))
#      id   trait rep      value
# 1  4226 Trait_1   1  1.0748090
# 2  4722 Trait_1   1 -0.1587611
# 3 33-16 Trait_1   1  1.1578901
# 4 38-11 Trait_1   1 -2.2299774
# 5  A188 Trait_1   1  0.8716323
# 6  A239 Trait_1   1  0.4134403

Dominance, epistasis and variance QTL

The other three layers work the same way. By default dominance() and vqtl() reuse the additive QTNs; set same_as_add = FALSE to draw fresh ones.

ph_ade <- simulate_phenotype(geno, seed = 2) |>
  additive(prop = 0.4, n_qtn = 5) |>
  dominance(prop = 0.1) |>
  epistasis(prop = 0.1, n_pairs = 2)
# 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.
ph_ade
# <phenotype_sim>  (realized · long format)
#   Genotypes: geno   Traits: 1   Architecture: independent   Seed: 2
#   Variance partition (proportions of V_P):
#     additive    0.40   (5 QTNs, geometric)
#     dominance   0.10   (same QTNs as additive)
#     epistasis   0.10   (2 pairs, 2-way)
#     residual    0.40
#   Requested genetic share = 0.60   realized H² = 0.55
#   Additive + dominance share loci (variance-partition coding); realized
#   shares of V_P for the A+D block:
#     Trait_1: requested 0.50 (unit-variance scale), realized 0.45
#       = Var(A) 0.35 + Var(D) 0.11 + 2Cov(A,D) -0.00
#       = Var(c_A) 0.37 + Var(c_D) 0.09 + 2Cov(c_A,c_D) -0.02
#   ! 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 =).
ph_v <- simulate_phenotype(geno, seed = 3) |>
  additive(prop = 0.3, n_qtn = 4) |>
  vqtl(prop = 0.2)
ph_v
# <phenotype_sim>  (realized · long format)
#   Genotypes: geno   Traits: 1   Architecture: independent   Seed: 3
#   Variance partition (proportions of V_P):
#     additive    0.30   (4 QTNs, geometric)
#     vqtl        0.20   (same QTNs as additive)
#     residual    0.50
#   Requested genetic share = 0.30   realized H² = 0.31
#   (vqtl is a residual-heterogeneity component and is not counted in
#    broad-sense heritability)

Multiple traits and genetic correlation

architecture = "pleiotropy" puts the same QTNs behind several traits. With two traits, cor sets the genetic correlation through the PleioArch algorithm, which draws the pleiotropic QTN effects from a bivariate normal.

ph_pleio <- simulate_phenotype(geno, architecture = "pleiotropy",
                               n_traits = 2, cor = 0.6, seed = 4) |>
  additive(prop = 0.5, n_qtn = 100)
ph_pleio
# <phenotype_sim>  (realized · long format)
#   Genotypes: geno   Traits: 2   Architecture: pleiotropy   Seed: 4
#   Variance partition (proportions of V_P):
#     additive    0.50   (100 QTNs, geometric)
#     residual    [0.50, 0.50]
#   Requested genetic share = [0.50, 0.50]   realized H² = [0.52, 0.50]

wide <- phenotypes_wide(ph_pleio)
round(cor(wide$Trait_1, wide$Trait_2), 2)
# [1] 0.3

Two things are worth knowing before you compare that number with cor.

It is measured on the wrong scale. cor is the genetic correlation, but what you observe above is the phenotypic one, and the residuals are drawn independently. With uncorrelated residuals the two are related by \(r_P \approx r_G \sqrt{h^2_1 h^2_2}\), so a genetic correlation of 0.6 between traits with \(h^2 = 0.5\) shows up as a phenotypic correlation near 0.3.

It is realized, not fixed. cor is the target the effect draws aim at, not a value hit exactly in any one simulation. The QTN effects are random draws, so the realized correlation scatters around the target – with few shared QTNs it is also pulled slightly toward 0 on average – and the scatter narrows as QTNs are added (for loci in approximate linkage equilibrium), down to the sampling spread of a correlation over the individuals you simulated. Measuring the genetic correlation directly, over 30 seeds:

scatter <- data.frame(n_qtn = c(20, 100, 400), mean = NA_real_, sd = NA_real_)

for (i in seq_len(nrow(scatter))) {
  realized <- numeric(30)
  for (s in 1:30) {
    sim <- simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 2,
                              cor = 0.6, h2 = 0.5, seed = s) |>
      additive(n_qtn = scatter$n_qtn[i])
    g <- genetic_values(sim)
    realized[s] <- cor(g[, 1], g[, 2])
  }
  scatter$mean[i] <- round(mean(realized), 3)
  scatter$sd[i]   <- round(sd(realized), 3)
}

scatter
#   n_qtn  mean    sd
# 1    20 0.587 0.135
# 2   100 0.622 0.075
# 3   400 0.604 0.058

The mean stays near the target and the scatter narrows as QTNs are added, but on a real panel not indefinitely, and not as \(1/\sqrt{n_{qtn}}\): the correlation is measured over a fixed number of individuals, whose sampling spread no number of QTNs removes, linkage disequilibrium among the causal loci limits it further, and with few shared QTNs the realized value is also pulled slightly toward 0 on average. If you need a tight correlation in one particular data set, use many QTNs — and check what you actually got rather than assuming.

Any number of traits. The same engine covers three, five or twenty traits — cor is either a scalar applied to every pair or a full n_traits x n_traits matrix. The genetic covariance matrix must be positive semi-definite, which for two traits is the familiar cor^2 <= pi_1 * pi_2 and for more traits also rules out mutually inconsistent requests; impossible targets are an error, not a silent approximation.

R <- matrix(c(1, 0.8, -0.4,
              0.8, 1, -0.2,
             -0.4, -0.2, 1), 3, 3)

# Averaged over 20 seeds, since any single simulation scatters around the target
total <- matrix(0, 3, 3)
for (s in 1:20) {
  ph3 <- simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 3,
                            cor = R, h2 = 0.5, seed = s) |>
    additive(n_qtn = 300)
  total <- total + cor(genetic_values(ph3))
}

round(total / 20, 2)
#         Trait_1 Trait_2 Trait_3
# Trait_1    1.00    0.79   -0.39
# Trait_2    0.79    1.00   -0.19
# Trait_3   -0.39   -0.19    1.00

Major QTNs trade precision for realism. n_pleio_major and prop_var_major concentrate part of the pleiotropic variance in a few large-effect loci. Both default to zero, because a locus carrying a large share of the variance makes the realized correlation hinge on that one effect draw, and the convergence above disappears — the standard deviation stops improving with n_qtn altogether. Set them when a major pleiotropic gene is the object of study, not otherwise.

Linked, not pleiotropic: the “ld” architecture

architecture = "ld" is the spurious-pleiotropy scenario: two traits whose distinct causal loci sit in linkage disequilibrium, so the traits covary through linkage even though no locus is causal for both. It is a two-trait construct (n_traits = 2) — one causal SNP per trait per locus-pair. ld_type picks how the pair is linked: "direct" (the default) puts trait 1’s and trait 2’s causal SNPs directly in LD; "indirect" links each to a shared hidden, non-causal “cause-of-LD” locus. r2_min/r2_max set the window the linked pair’s squared correlation must fall in. qtn_table() reports the linked pair on each row: QTN_t1 / QTN_t2 (the two traits’ causal SNPs) and ld_r2 (the r² between them).

ph_ld <- simulate_phenotype(geno, architecture = "ld", n_traits = 2,
                            seed = 5) |>
  additive(prop = 0.4, n_qtn = 3)
ph_ld
# <phenotype_sim>  (realized · long format)
#   Genotypes: geno   Traits: 2   Architecture: ld   Seed: 5
#   Variance partition (proportions of V_P):
#     additive    0.40   (3 QTNs, geometric)
#     residual    [0.60, 0.60]
#   Requested genetic share = [0.40, 0.40]   realized H² = [0.41, 0.37]

# distinct causal loci per trait, with the r2 between 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 ss196430644 ss196430644 ss196430642 0.6130081
# 2 Trait_1 ss196512868 ss196512868 ss196512866 0.4963766
# 3 Trait_1 ss196475922 ss196475922 ss196475924 0.7901685
# 4 Trait_2 ss196430642 ss196430644 ss196430642 0.6130081
# 5 Trait_2 ss196512866 ss196512868 ss196512866 0.4963766
# 6 Trait_2 ss196475924 ss196475922 ss196475924 0.7901685

Combining architectures

complex_phenotypes() merges complete models under a common heritability, weighting them by their genetic variances. This is how partial pleiotropy is built: combine a pleiotropic model with an independent one.

pleio <- simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 2,
                            cor = 0.5, seed = 6) |>
  additive(prop = 0.4, n_qtn = 5)
indep <- simulate_phenotype(geno, n_traits = 2, seed = 6) |>
  additive(prop = 0.3, n_qtn = 5)

both <- complex_phenotypes(pleio, indep, h2 = 0.5)
both
# <phenotype_sim>  (realized · long format)
#   Genotypes: geno   Traits: 2   Architecture: complex   Seed: 6
#   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.55, 0.50]

One call instead of a pipe

Supplying h2 together with n_qtn realizes a complete model in a single call.

simulate_phenotype(geno, h2 = 0.6, n_qtn = 5, model = "AD", seed = 7)
# 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: 7
#   Variance partition (proportions of V_P):
#     additive    0.30   (5 QTNs, geometric)
#     dominance   0.30   (same QTNs as additive)
#     residual    0.40
#   Requested genetic share = 0.60   realized H² = 0.62
#   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.62
#       = Var(A) 19.11 + Var(D) 15.34 + 2Cov(A,D) -33.83
#       = Var(c_A) 0.33 + Var(c_D) 0.33 + 2Cov(c_A,c_D) -0.04
#   ! 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 =).

Writing results

Phenotypes come back as an object; write_phenotypes() writes them to a file you name, either as a long table (one row per individual and trait) or a wide one (one column per trait). sep sets the field separator.

out_long <- file.path(tempdir(), "phenotypes_long.txt")
out_wide <- file.path(tempdir(), "phenotypes_wide.csv")
write_phenotypes(ph_pleio, file = out_long, format = "long")
write_phenotypes(ph_pleio, file = out_wide, format = "wide", sep = ",")
unlink(c(out_long, out_wide))
head(phenotypes_wide(ph_pleio))
#      id rep    Trait_1    Trait_2
# 1  4226   1 -0.3901323  0.1663806
# 2  4722   1 -1.2747409  0.4291148
# 3 33-16   1  0.9748761 -1.0156686
# 4 38-11   1 -0.4453048 -0.7165219
# 5  A188   1 -0.2157302  0.6824694
# 6  A239   1  0.6451052 -0.8919026

Multi-generation crosses

The crossing functions simulate meiosis directly, so they can build a pedigree from real founder genotypes. Phase the founders into a Population, then cross, self or make doubled haploids:

pop <- as_population(geno, individuals = c("33-16", "38-11"))
f1  <- cross(pop[1], pop[2], n = 1, seed = 11)
f2  <- selfcross(f1, n = 200, seed = 12)
f2
# <Population>
#   Individuals: 200   Markers: 10650   Chromosomes: 10
#   Genetic map: 1497 cM total span (109-219 cM per chromosome)
#   Origin: selfcross(prog_1)
#   Pedigree: 203 recorded individuals; generation 2
#   IDs: self_1, self_2, self_3, self_4, self_5, self_6, ... (194 more)

A Population is accepted anywhere simulate_phenotype() takes geno, so a pedigree can be phenotyped without any conversion:

ph_f2 <- simulate_phenotype(f2, seed = 14) |>
  additive(prop = 0.6, n_qtn = 5)
ph_f2
# <phenotype_sim>  (realized · long format)
#   Genotypes: <Population: selfcross(prog_1)>   Traits: 1   Architecture: independent   Seed: 14
#   Variance partition (proportions of V_P):
#     additive    0.60   (5 QTNs, geometric)
#     residual    0.40
#   Requested genetic share = 0.60   realized H² = 0.54

vignette("breeding-populations") covers this properly: backcrosses, recurrent selfing to inbred lines, doubled haploids, and how to check that segregation and recombination are behaving.

Original implementation of create_phenotypes

create_phenotypes() is the original v1 function. It is frozen: it keeps its original code paths and seed handling so that existing scripts keep reproducing their published results, and it receives bug fixes only. New work should use the grammar above, which is why it is documented as superseded.

create_phenotypes(
  geno_obj = geno,
  add_QTN_num = 3,
  add_effect = 0.2,
  rep = 10,
  h2 = 0.7,
  model = "A",
  to_r = TRUE,
  home_dir = tempdir()
)

Session information

sessionInfo()
# R version 4.6.1 (2026-06-24)
# Platform: x86_64-pc-linux-gnu
# Running under: Ubuntu 26.04.1 LTS
# 
# Matrix products: default
# BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
# LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
# 
# locale:
#  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
# [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
# 
# time zone: Etc/UTC
# tzcode source: system (glibc)
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
# [1] simplePHENOTYPES_2.0.0.9003 rmarkdown_2.32             
# 
# loaded via a namespace (and not attached):
#  [1] digest_0.6.39       R6_2.6.1            fastmap_1.2.0      
#  [4] xfun_0.61           maketools_1.3.2     cachem_1.1.0       
#  [7] knitr_1.52          RhpcBLASctl_0.23-42 htmltools_0.5.9    
# [10] buildtools_1.0.0    lifecycle_1.0.5     cli_3.6.6          
# [13] gdsfmt_1.49.8       SNPRelate_1.47.0    sass_0.4.10        
# [16] data.table_1.18.6.1 jquerylib_0.1.4     compiler_4.6.1     
# [19] sys_3.4.3           tools_4.6.1         evaluate_1.0.5     
# [22] bslib_0.12.0        yaml_2.3.12         otel_0.2.0         
# [25] jsonlite_2.0.0      rlang_1.3.0