This vignette exercises every option added in version 2, end to end, in one runnable document. It is deliberately exhaustive rather than gentle: if you want a guided introduction, start with vignette("simplePHENOTYPES-v2"), then vignette("breeding-populations") and vignette("genetic-maps").

Everything below actually runs when this vignette is built, so the printed output is real, not illustrative.

library(simplePHENOTYPES)
data("SNP55K_maize282_maf04")
geno <- SNP55K_maize282_maf04
dim(geno)
#> [1] 10650   285

The first five columns are marker metadata (snp, allele, chr, pos, cm); the rest are individuals.

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

1. Genotype input

Converting from other formats

as_numeric() reads HapMap, VCF, GDS and PLINK bed/ped and returns the numeric format used throughout.

hmp <- data.frame(
  `rs#`   = c("snp1", "snp2", "snp3"),
  alleles = c("C/T", "A/G", "G/T"),
  chrom   = c(1, 1, 2),
  pos     = c(100, 500, 200),
  strand  = "+", `assembly#` = NA, center = NA,
  protLSID = NA, assayLSID = NA, panelLSID = NA, QCcode = NA,
  L1 = c("C", "A", "G"), L2 = c("T", "G", "T"),
  L3 = c("C", "A", "T"), L4 = c("T", "G", "G"),
  check.names = FALSE, stringsAsFactors = FALSE
)

as_numeric(hmp, to_r = TRUE, verbose = FALSE)
#>    snp allele chr pos cm L1 L2 L3 L4
#> 1 snp1    C/T   1 100 NA  1 -1  1 -1
#> 2 snp2    A/G   1 500 NA  1 -1  1 -1
#> 3 snp3    G/T   2 200 NA  1 -1 -1  1

Coding schemes and allele orientation:

# 0 / 1 / 2 instead of -1 / 0 / 1
as_numeric(hmp, to_r = TRUE, code_as = "012", verbose = FALSE)[, 1:8]
#>    snp allele chr pos cm L1 L2 L3
#> 1 snp1    C/T   1 100 NA  2  0  2
#> 2 snp2    A/G   1 500 NA  2  0  2
#> 3 snp3    G/T   2 200 NA  2  0  0

# Orient against a supplied reference allele rather than by frequency
as_numeric(hmp, to_r = TRUE, verbose = FALSE,
           method = "reference", ref_allele = c("C", "A", "G"))[, 1:8]
#>    snp allele chr pos cm L1 L2 L3
#> 1 snp1    C/T   1 100 NA  1 -1  1
#> 2 snp2    A/G   1 500 NA  1 -1  1
#> 3 snp3    G/T   2 200 NA  1 -1 -1

Reading from files (not run here, since it needs your data):

as_numeric("my_genotypes.hmp.txt", to_r = TRUE)
as_numeric("my_genotypes.vcf",     to_r = TRUE)
as_numeric("my_genotypes.bed",     to_r = TRUE)   # needs .bim and .fam alongside

Filtering markers with `filter_geno()`

filter_geno() is one place to trim a genotype before simulating: it drops monomorphic markers, filters by minor-allele frequency and by heterozygosity, and prunes by linkage disequilibrium. A marker must pass every requested filter.

data("SNP55K_maize282_maf04")
geno <- SNP55K_maize282_maf04

# MAF window + only markers that carry heterozygotes (e.g. before dominance)
g <- filter_geno(geno, maf_below = 0.48, hets = "include", verbose = FALSE)
nrow(g)
#> [1] 3884

Heterozygosity matters for dominance: a dominance deviation acts on heterozygotes, so on a near-inbred panel dominance() errors unless the loci have heterozygotes. Pre-filtering fixes that:

simulate_phenotype(g, seed = 60) |>
  additive(prop = 0.4, n_qtn = 5) |>
  dominance(prop = 0.1)
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: g   Traits: 1   Architecture: independent   Seed: 60
#>   Variance partition (proportions of V_P):
#>     additive    0.40   (5 QTNs, geometric)
#>     dominance   0.10   (same QTNs as additive)
#>     residual    0.50
#>   Requested genetic share = 0.50   realized H² = 0.49
#>   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.49
#>       = Var(A) 0.59 + Var(D) 0.17 + 2Cov(A,D) -0.26
#>       = Var(c_A) 0.37 + Var(c_D) 0.09 + 2Cov(c_A,c_D) 0.03
#>   ! 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 =).

LD pruning and haplotype blocks. The three sliding-window pruners take a c(window, step, threshold) specification; blocks groups markers into Gabriel et al. (2002) haplotype blocks and keeps one tag per block. (Shown on a single chromosome to keep the vignette quick.)

chr1 <- geno[geno$chr == geno$chr[1], ][1:400, ]

# pairwise composite-r2 pruning (r2 > 0.2 in 50-marker windows)
nrow(filter_geno(chr1, indep_pairwise = c(50, 5, 0.2), verbose = FALSE))
#> [1] 272

# phased (EM haplotype) r2 pruning
nrow(filter_geno(chr1, indep_pairphase = c(50, 5, 0.2), verbose = FALSE))
#> [1] 272

# variance-inflation-factor pruning (PLINK's --indep)
nrow(filter_geno(chr1, indep = c(50, 5, 2), verbose = FALSE))
#> [1] 262

# Gabriel haplotype blocks: keep one tag marker per block
nrow(filter_geno(chr1, blocks = TRUE, block_max_kb = 2000, verbose = FALSE))
#> [1] 336

Using the VIF pruner or the block definition prints a one-per-session citation for the method it implements (PLINK’s VIF procedure, or Gabriel et al. for the blocks); the pairwise pruners are generic r^2 operations and print nothing.

2. The genetic map

Meiosis needs distances in centiMorgans. SNP55K_maize282_maf04 ships with a synthetic map in cm — modelled from the physical positions, not a published maize linkage map, and not to be quoted as measured recombination distance.

chrom_length <- aggregate(cm ~ chr, data = geno, FUN = max)
chrom_length
#>    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

round(sum(chrom_length$cm))             # total genetic length, cM
#> [1] 1497

synthetic_map() builds such a map for your own data. Every argument:

# Default: pericentromeric suppression, length from cm_per_mb
cm_default <- synthetic_map(geno$chr, geno$pos)

# Uniform map: cM strictly proportional to bp
cm_uniform <- synthetic_map(geno$chr, geno$pos, suppression = 0)

# Fixed total length per chromosome
cm_fixed <- synthetic_map(geno$chr, geno$pos,
                          total_cm = c(200, 180, 175, 170, 160,
                                       130, 130, 125, 115, 110))

# Explicit centromere positions, and a wider/stronger suppressed region
chr1 <- geno$pos[geno$chr == 1]
mid <- aggregate(pos ~ chr, data = geno, FUN = median)
cm_custom <- synthetic_map(geno$chr, geno$pos,
                           centromere = mid$pos,
                           suppression = 0.95, width = 0.25)

round(c(default = max(cm_default[geno$chr == 1]),
        uniform = max(cm_uniform[geno$chr == 1]),
        fixed   = max(cm_fixed[geno$chr == 1]),
        custom  = max(cm_custom[geno$chr == 1])), 1)
#> default uniform   fixed  custom 
#>   219.4   219.4   200.0   219.4

The suppression model concentrates recombination on the chromosome arms:

along_chr1 <- data.frame(
  pos = chr1,
  cm  = cm_default[geno$chr == 1],
  bin = cut(chr1, breaks = 8, labels = FALSE)
)

span <- aggregate(cbind(cm, pos) ~ bin, data = along_chr1,
                  FUN = function(x) max(x) - min(x))
span$cM_per_Mb <- round(span$cm / (span$pos / 1e6), 2)

span[, c("bin", "cM_per_Mb")]     # low near the centromere, high on the arms
#>   bin cM_per_Mb
#> 1   1      1.06
#> 2   2      0.96
#> 3   3      0.65
#> 4   4      0.26
#> 5   5      0.26
#> 6   6      0.65
#> 7   7      0.96
#> 8   8      1.06

To use your own real map, just assign it:

geno$cm <- my_real_map_in_cM

3. The simulation grammar

A simulation starts with simulate_phenotype() and gains variance components as layers. Each prop is a proportion of total phenotypic variance, so heritability is the sum of the genetic proportions.

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

Every layer

full <- simulate_phenotype(geno, seed = 2) |>
  additive(prop = 0.3, n_qtn = 5) |>
  dominance(prop = 0.1) |>
  epistasis(prop = 0.1, n_pairs = 2, interaction = 2) |>
  vqtl(prop = 0.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.
#> In addition to citing:
#> Fernandes, S.B., Lipka, A.E. simplePHENOTYPES: SIMulation of pleiotropic, linked and epistatic phenotypes. BMC Bioinformatics 21, 491 (2020). https://doi.org/10.1186/s12859-020-03804-y
#> Please also cite:
#> Murphy et al. (2022), Heredity 129:93-102, doi:10.1038/s41437-022-00541-1, for the plant variance-QTL (vGWAS) context (the log-variance link used here is this package's own simulation choice, not Murphy et al.'s generator).
#> This message is displayed once per session.
full
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: geno   Traits: 1   Architecture: independent   Seed: 2
#>   Variance partition (proportions of V_P):
#>     additive    0.30   (5 QTNs, geometric)
#>     dominance   0.10   (same QTNs as additive)
#>     epistasis   0.10   (2 pairs, 2-way)
#>     vqtl        0.10   (same QTNs as additive)
#>     residual    0.40
#>   Requested genetic share = 0.50   realized H² = 0.47
#>   Additive + dominance share loci (variance-partition coding); realized
#>   shares of V_P for the A+D block:
#>     Trait_1: requested 0.40 (unit-variance scale), realized 0.36
#>       = Var(A) 0.26 + Var(D) 0.11 + 2Cov(A,D) -0.01
#>       = Var(c_A) 0.29 + Var(c_D) 0.10 + 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 =).
#>   (vqtl is a residual-heterogeneity component and is not counted in
#>    broad-sense heritability)

Layer options

# Geometric effects with an explicit base (0.5, 0.25, 0.125, ...)
simulate_phenotype(geno, seed = 3) |>
  additive(prop = 0.5, n_qtn = 3, effect = 0.5)
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: geno   Traits: 1   Architecture: independent   Seed: 3
#>   Variance partition (proportions of V_P):
#>     additive    0.50   (3 QTNs, geometric)
#>     residual    0.50
#>   Requested genetic share = 0.50   realized H² = 0.51

# Explicit effect series, one per QTN
simulate_phenotype(geno, seed = 3) |>
  additive(prop = 0.5, n_qtn = 3, effect = c(0.6, 0.3, 0.1))
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: geno   Traits: 1   Architecture: independent   Seed: 3
#>   Variance partition (proportions of V_P):
#>     additive    0.50   (3 QTNs, geometric)
#>     residual    0.50
#>   Requested genetic share = 0.50   realized H² = 0.51

# Dominance and vQTL on their own loci rather than the additive ones
simulate_phenotype(geno, seed = 4) |>
  additive(prop = 0.3, n_qtn = 5) |>
  dominance(prop = 0.1, same_as_add = FALSE, n_qtn = 3) |>
  vqtl(prop = 0.1, same_as_add = FALSE, n_qtn = 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.
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: geno   Traits: 1   Architecture: independent   Seed: 4
#>   Variance partition (proportions of V_P):
#>     additive    0.30   (5 QTNs, geometric)
#>     dominance   0.10   (3 QTNs)
#>     vqtl        0.10   (2 QTNs)
#>     residual    0.50
#>   Requested genetic share = 0.40   realized H² = 0.43
#>   (vqtl is a residual-heterogeneity component and is not counted in
#>    broad-sense heritability)

# Three-way epistasis
simulate_phenotype(geno, seed = 5) |>
  additive(prop = 0.3, n_qtn = 4) |>
  epistasis(prop = 0.2, n_pairs = 2, interaction = 3)
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: geno   Traits: 1   Architecture: independent   Seed: 5
#>   Variance partition (proportions of V_P):
#>     additive    0.30   (4 QTNs, geometric)
#>     epistasis   0.20   (2 pairs, 3-way)
#>     residual    0.50
#>   Requested genetic share = 0.50   realized H² = 0.56

Replications and per-trait proportions

reps <- simulate_phenotype(geno, n_traits = 2, n_reps = 3, seed = 6) |>
  additive(prop = c(0.6, 0.3), n_qtn = 4)     # different h2 per trait
table(phenotypes_long(reps)$rep)
#> 
#>   1   2   3 
#> 560 560 560
head(phenotypes_wide(reps))
#>      id rep    Trait_1    Trait_2
#> 1  4226   1  0.4445307  0.9735806
#> 2  4722   1  0.6251630  0.9778744
#> 3 33-16   1 -1.3646453  0.3652453
#> 4 38-11   1 -1.0963849  0.8503272
#> 5  A188   1 -1.4675454 -0.9355720
#> 6  A239   1  0.6525548  0.2324499

One call instead of a pipe

simulate_phenotype(geno, h2 = 0.5, n_qtn = 3, seed = 7)              # model "A"
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: geno   Traits: 1   Architecture: independent   Seed: 7
#>   Variance partition (proportions of V_P):
#>     additive    0.50   (3 QTNs, geometric)
#>     residual    0.50
#>   Requested genetic share = 0.50   realized H² = 0.53
simulate_phenotype(geno, h2 = 0.6, n_qtn = 4, 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   (4 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.46 + Var(D) 15.61 + 2Cov(A,D) -34.45
#>       = Var(c_A) 0.33 + Var(c_D) 0.33 + 2Cov(c_A,c_D) -0.03
#>   ! 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 =).

Fine control: fixed loci, phase, means, subsets

Set the causal loci for one layer while the rest are drawn at random. Give marker names or column indices:

simulate_phenotype(geno, h2 = 0.5, seed = 1) |>
  additive(qtn = c("ss196442916", "ss196439337", "ss196480535"))
#> <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.52

phase = "repulsion" alternates the effect signs so linked increasing and decreasing alleles oppose one another (the realized variance is still prop; only the sign structure changes):

ph_rep <- simulate_phenotype(geno, h2 = 0.5, seed = 1) |>
  additive(n_qtn = 6, phase = "repulsion")
qtn_table(ph_rep)$effect
#> [1]  0.500000 -0.250000  0.125000 -0.062500  0.031250 -0.015625

A per-trait intercept with mean, and a subset of individuals with individuals:

ph_m <- simulate_phenotype(geno, n_traits = 2, h2 = 0.5, mean = c(10, 20),
                           seed = 1) |>
  additive(n_qtn = 5)
colMeans(phenotypes_wide(ph_m)[, c("Trait_1", "Trait_2")])
#> Trait_1 Trait_2 
#>      10      20

ph_s <- simulate_phenotype(geno, h2 = 0.5,
                           individuals = c("33-16", "38-11", "A188", "B73"),
                           seed = 1) |>
  additive(n_qtn = 3)
n_individuals_sim <- ph_s$n_ind
n_individuals_sim
#> [1] 4

vary_qtn = TRUE gives each replication its own QTNs, and distinct_chr = TRUE puts each trait’s QTNs on separate chromosomes:

ph_v <- simulate_phenotype(geno, h2 = 0.5, n_reps = 3, vary_qtn = TRUE,
                           seed = 1) |>
  additive(n_qtn = 4)
# a different QTN set per replication
lengths(ph_v$layers[[1]]$qtn_reps)
#> [1] 1 1 1

ph_d <- simulate_phenotype(geno, n_traits = 2, h2 = 0.4, distinct_chr = TRUE,
                           seed = 1) |>
  additive(n_qtn = 3)
lapply(ph_d$layers[[1]]$qtn, function(i) unique(geno$chr[i]))
#> [[1]]
#> [1] 1 3
#> 
#> [[2]]
#> [1]  4 10  6

4. Genetic architectures

Independent

Trait-specific QTNs, no enforced sharing.

ind <- simulate_phenotype(geno, n_traits = 2, h2 = 0.4, seed = 8) |>
  additive(n_qtn = 5)

qtn <- ind$layers[[1]]$qtn
head(qtn[[1]], 3)      # trait 1
#> [1] 7872 6268 9136
head(qtn[[2]], 3)      # trait 2: different loci
#> [1] 7671 1033 2806

Pleiotropy with a target genetic correlation

Shared QTNs whose effects are drawn from a multivariate normal with cross-trait covariance cor * sqrt(V1 * V2), so the realized genetic correlation targets cor – scattering around it, and converging as the shared QTNs and the individuals grow for loci in approximate linkage equilibrium. Works for any number of traits.

ph2 <- simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 2,
                          cor = 0.6, seed = 9) |>
  additive(prop = 0.5, n_qtn = 200)
#> In addition to citing:
#> Fernandes, S.B., Lipka, A.E. simplePHENOTYPES: SIMulation of pleiotropic, linked and epistatic phenotypes. BMC Bioinformatics 21, 491 (2020). https://doi.org/10.1186/s12859-020-03804-y
#> Please also cite:
#> Prado et al. (in preparation), when controlling the correlation in the pleiotropic architecture.
#> This message is displayed once per session.

g <- genetic_values(ph2)
round(cor(g[, 1], g[, 2]), 3)
#> [1] 0.517

Any number of traits, with a full target matrix including negative correlations:

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

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

round(total / 20, 2)      # compare with R above
#>         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

Partial pleiotropy — only some of each trait’s genetic variance is shared — via pi:

part <- simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 3,
                           cor = 0.5, pi = c(1, 0.8, 0.6), seed = 10) |>
  additive(prop = 0.5, n_qtn = 100)

# Every trait keeps n_qtn loci; a shared core carries the covariance
qtn_pi <- part$layers[[1]]$qtn
lengths(qtn_pi)
#> [1]  80 100 100

shared <- intersect(intersect(qtn_pi[[1]], qtn_pi[[2]]), qtn_pi[[3]])
length(shared)
#> [1] 80

Concentrating variance in major pleiotropic loci. This is available but costs precision in the realized correlation, so both arguments default to zero:

simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 2,
                   cor = 0.6, n_pleio_major = 5, prop_var_major = 0.5,
                   seed = 11) |>
  additive(prop = 0.5, n_qtn = 100)
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: geno   Traits: 2   Architecture: pleiotropy   Seed: 11
#>   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.49, 0.49]

Impossible requests are an error rather than a silent approximation:

simulate_phenotype(geno, architecture = "pleiotropy", n_traits = 2,
                   cor = 0.9, pi = 0.5, seed = 12) |>
  additive(prop = 0.5, n_qtn = 20)
#> Error:
#> ! Biological constraint violated: cor^2 (0.81) cannot exceed pi_target * pi_secondary (0.25).

Linkage disequilibrium

The causal variant sits near, but not on, the marker a study would test.

simulate_phenotype(geno, architecture = "ld", n_traits = 2,
                   ld_type = "indirect", seed = 13) |>
  additive(prop = 0.4, n_qtn = 3)
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: geno   Traits: 2   Architecture: ld   Seed: 13
#>   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.37, 0.41]

simulate_phenotype(geno, architecture = "ld", n_traits = 2,
                   ld_type = "direct", r2_max = 0.9, r2_min = 0.4, seed = 14) |>
  additive(prop = 0.4, n_qtn = 3)
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: geno   Traits: 2   Architecture: ld   Seed: 14
#>   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.39, 0.43]

Combining architectures

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

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

5. Breeding populations

Meiosis is simulated directly, so a pedigree can be built from real founders and phenotyped without leaving the package.

pop <- as_population(geno, individuals = c("33-16", "38-11", "A188", "B73"))
pop
#> <Population>
#>   Individuals: 4   Markers: 10650   Chromosomes: 10
#>   Genetic map: 1497 cM total span (109-219 cM per chromosome)
#>   Origin: founder
#>   Pedigree: 4 recorded individuals; generation 0
#>   IDs: 33-16, 38-11, A188, B73
n_individuals(pop)
#> [1] 4
dim(dosages(pop))
#> [1] 10650     4

Crossing, selfing, doubled haploids

f1  <- cross(pop[1], pop[2], n = 1, seed = 21)
f2  <- selfcross(f1, n = 200, seed = 22)
dh  <- double_haploid(f1, n = 200, seed = 23)
bc1 <- cross(f1, pop[1], n = 100, seed = 24)      # backcross to founder 1

c(F1 = n_individuals(f1), F2 = n_individuals(f2),
  DH = n_individuals(dh), BC1 = n_individuals(bc1))
#>  F1  F2  DH BC1 
#>   1 200 200 100

Checking the genetics came out right:

d <- dosages(pop)
inf <- which(d[, 1] != d[, 2] & d[, 1] != 0 & d[, 2] != 0)  # informative loci

# F1 is heterozygous everywhere the founders differ
all(dosages(f1)[inf, 1] == 0)
#> [1] TRUE

# F2 segregates 1:2:1
round(prop.table(table(dosages(f2)[inf, ])), 3)
#> 
#>    -1     0     1 
#> 0.248 0.502 0.249

# Doubled haploids are fully homozygous, ~1:1
any(dosages(dh) == 0)
#> [1] FALSE
round(prop.table(table(dosages(dh)[inf, ])), 3)
#> 
#>    -1     1 
#> 0.501 0.499

# Backcross to founder 1: half heterozygous, half like the recurrent parent
round(prop.table(table(dosages(bc1)[inf, ])), 3)
#> 
#>    -1     0     1 
#> 0.224 0.510 0.266

Recurrent selfing

Single-seed descent: each generation one plant is selfed, and the next generation descends from one of its progeny.

gen <- f1
het <- numeric(5)
for (i in 1:5) {
  gen <- selfcross(gen[1], n = 100, seed = 30 + i)
  het[i] <- mean(dosages(gen)[inf, ] == 0)
}
round(het, 3)     # heterozygosity halves each generation
#> [1] 0.500 0.248 0.110 0.083 0.042

Reproducibility

identical(dosages(selfcross(f1, n = 10, seed = 99)),
          dosages(selfcross(f1, n = 10, seed = 99)))
#> [1] TRUE

Phenotyping a pedigree

A Population is accepted anywhere simulate_phenotype() takes genotypes.

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

6. Inspecting a simulation

Nothing is written to disk unless you ask for it. Everything the original create_phenotypes() wrote to Log_Sim.txt and its companion files lives in the returned object instead.

print() shows the requested variance budget alongside what was actually realized:

# dominance needs heterozygous loci; on this near-inbred panel filter to them
# first (otherwise dominance() errors -- see filter_geno() above)
gh <- filter_geno(geno, hets = "include", verbose = FALSE)
ph_i <- simulate_phenotype(gh, h2 = 0.5, seed = 60) |>
  additive(prop = 0.3, n_qtn = 5) |>
  dominance(prop = 0.2)
ph_i
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: gh   Traits: 1   Architecture: independent   Seed: 60
#>   Variance partition (proportions of V_P):
#>     additive    0.30   (5 QTNs, geometric)
#>     dominance   0.20   (same QTNs as additive)
#>     residual    0.50
#>   Requested genetic share = 0.50   realized H² = 0.52
#>   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.52
#>       = Var(A) 0.11 + Var(D) 0.33 + 2Cov(A,D) 0.08
#>       = Var(c_A) 0.32 + Var(c_D) 0.21 + 2Cov(c_A,c_D) -0.01
#>   ! 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 =).

genetic_values() returns the genetic component, before the residual — the equivalent of the old Genetic_values.txt:

gv <- genetic_values(ph_i)
head(gv)
#>           Trait_1
#> 4226   0.19793736
#> 4722   0.68309547
#> 33-16  0.25858212
#> 38-11 -0.22657598
#> A188   0.19793736
#> A239   0.01600307

# Realized heritability, computed by hand from the same pieces
round(var(gv[, 1]) / var(phenotypes_wide(ph_i)$Trait_1), 3)
#> [1] 0.52

The var_explained column of qtn_table() gives each additive/dominance QTN’s marginal share of phenotypic variance (NA for epistatic sets, where variance is not per-locus). qtn_table() names the causal loci and their effects — the equivalent of Additive_QTNs.txt and QTN_effects_summary.txt. Epistatic sets are expanded one row per member, sharing a set number and effect:

qtn_table(ph_i)
#>      trait     layer set         snp chr       pos       maf  effect  d
#> 1  Trait_1  additive  NA ss196506235  10   3989555 0.4178571 0.50000 NA
#> 2  Trait_1  additive  NA ss196490729   9  58451982 0.4839286 0.25000 NA
#> 3  Trait_1  additive  NA ss196478754   7 103537211 0.4303571 0.12500 NA
#> 4  Trait_1  additive  NA ss196416215   5 176219732 0.4982143 0.06250 NA
#> 5  Trait_1  additive  NA ss196473344   6 114188703 0.4660714 0.03125 NA
#> 6  Trait_1 dominance  NA ss196506235  10   3989555 0.4178571 0.50000 NA
#> 7  Trait_1 dominance  NA ss196490729   9  58451982 0.4839286 0.25000 NA
#> 8  Trait_1 dominance  NA ss196478754   7 103537211 0.4303571 0.12500 NA
#> 9  Trait_1 dominance  NA ss196416215   5 176219732 0.4982143 0.06250 NA
#> 10 Trait_1 dominance  NA ss196473344   6 114188703 0.4660714 0.03125 NA
#>    var_explained QTN_t1 QTN_t2 ld_r2
#> 1   0.2339089856   <NA>   <NA>    NA
#> 2   0.0621020335   <NA>   <NA>    NA
#> 3   0.0146819794   <NA>   <NA>    NA
#> 4   0.0038853560   <NA>   <NA>    NA
#> 5   0.0009529365   <NA>   <NA>    NA
#> 6   0.1690325144   <NA>   <NA>    NA
#> 7   0.0043666733   <NA>   <NA>    NA
#> 8   0.0115779447   <NA>   <NA>    NA
#> 9   0.0002729171   <NA>   <NA>    NA
#> 10  0.0003362554   <NA>   <NA>    NA

Everything else the log recorded is a field on the object:

ph_i$seed
#> [1] 60
ph_i$n_traits
#> [1] 1
ph_i$var_budget
#>     trait component prop
#> 1 Trait_1  additive  0.3
#> 2 Trait_1 dominance  0.2
#> 3 Trait_1  residual  0.5

Note that a vqtl() layer is an explicit residual-heterogeneity variance share, not part of the broad-sense h2 budget. It varies the conditional residual variance rather than adding genetic value; print() says so when one is present:

simulate_phenotype(geno, h2 = 0.5, seed = 61) |>
  additive(prop = 0.3, n_qtn = 4) |>
  vqtl(prop = 0.2)
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: geno   Traits: 1   Architecture: independent   Seed: 61
#>   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
#>   ! Incomplete h2 allocation: genetic layers sum to 0.30 of h2 = 0.50; extracting phenotypes will error until the budget is filled (SPEC 4.1).
#>   (vqtl is a residual-heterogeneity component and is not counted in
#>    broad-sense heritability)

plot() gives a four-panel diagnostic summary – variance partition, phenotype distribution, QTN effects, and (for two or more traits) the genetic-value scatter:

ph_plot <- simulate_phenotype(geno, n_traits = 2, h2 = 0.5, seed = 62) |>
  additive(prop = 0.4, n_qtn = 8) |>
  dominance(prop = 0.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.
plot(ph_plot)

7. Output

In memory

ph3 <- simulate_phenotype(geno, n_traits = 2, seed = 50) |>
  additive(prop = 0.5, n_qtn = 3)

head(phenotypes_long(ph3))
#>      id   trait rep      value
#> 1  4226 Trait_1   1 -1.4368624
#> 2  4722 Trait_1   1 -0.3112517
#> 3 33-16 Trait_1   1  0.3905008
#> 4 38-11 Trait_1   1  1.5292108
#> 5  A188 Trait_1   1 -1.0562369
#> 6  A239 Trait_1   1  0.2306071
head(phenotypes_wide(ph3))
#>      id rep    Trait_1     Trait_2
#> 1  4226   1 -1.4368624  1.11941705
#> 2  4722   1 -0.3112517 -0.06187057
#> 3 33-16   1  0.3905008 -0.82000948
#> 4 38-11   1  1.5292108 -1.34667914
#> 5  A188   1 -1.0562369  1.03062968
#> 6  A239   1  0.2306071  1.93143009

Written to disk

Each call writes exactly the file you name — nothing else is created alongside it.

out_dir <- file.path(tempdir(), "sp_demo")
dir.create(out_dir, showWarnings = FALSE)

write_phenotypes(ph3, file = file.path(out_dir, "phenotypes_long.txt"))
write_phenotypes(ph3, file = file.path(out_dir, "phenotypes_wide.csv"),
                 format = "wide", sep = ",")

list.files(out_dir)
#> [1] "phenotypes_long.txt" "phenotypes_wide.csv"
head(read.delim(file.path(out_dir, "phenotypes_long.txt")))
#>      id   trait rep      value
#> 1  4226 Trait_1   1 -1.4368624
#> 2  4722 Trait_1   1 -0.3112517
#> 3 33-16 Trait_1   1  0.3905008
#> 4 38-11 Trait_1   1  1.5292108
#> 5  A188 Trait_1   1 -1.0562369
#> 6  A239 Trait_1   1  0.2306071

8. Transcriptome and expression-mediated phenotypes

A transcriptome layer lets a phenotype depend on gene expression instead of (or alongside) the markers. There are two pieces: a generator that simulates a genetically controlled genes x individuals expression matrix, and a transcriptome() layer that turns expression into a phenotype under the same variance budget as every other layer.

Simulating expression

simulate_transcriptome() needs no reference expression and no gene annotation. Each gene gets cis-eQTL near its (synthetic) position, trans effects through a few co-expression modules, and a per-gene heritability h2:

tx <- simulate_transcriptome(geno, n_genes = 200, seed = 70)
tx
#> <transcriptome_sim>
#>   Genes: 200   Individuals: 280   Factors: 5
#>   Coordinates: synthetic   Profile: generic_bulk
#>   Realized expression h2 Var(G)/Var(P): median 0.16  [0.01, 0.63]
#>   cis-eQTL: 301 over 200 genes; trans hubs: 7
dim(tx$expression)               # genes x individuals
#> [1] 200 280
head(tx$genes[, c("gene_id", "chr", "h2_target", "h2_realized", "n_cis")])
#>    gene_id chr h2_target h2_realized n_cis
#> 1 gene0001   1 0.5242829   0.5367475     2
#> 2 gene0002  10 0.3370640   0.3658942     3
#> 3 gene0003   6 0.1026378   0.1000263     1
#> 4 gene0004   1 0.3476490   0.3403233     1
#> 5 gene0005   9 0.1899271   0.1981428     2
#> 6 gene0006   9 0.4289051   0.4395164     1
head(tx$cis_eqtl)                # the cis truth: gene, marker, effect
#>    gene_id         snp chr       pos     effect
#> 1 gene0001 ss196427424   1 146441165  0.6041505
#> 2 gene0001 ss196427410   1 146203350  0.4799341
#> 3 gene0002 ss196511886  10 123574157  0.1103348
#> 4 gene0002 ss196504582  10 123956232  0.2402038
#> 5 gene0002 ss196499533  10 122808698 -0.1422082
#> 6 gene0003 ss196472103   6  85206564 -0.1539420

Useful arguments: h2 (a single value, one per gene, or "beta" for a realistic spread), cis_fraction (share of genetic variance that is cis), cis_window (bp), n_factors and residual_module_fraction (co-expression), and annotation (real gene coordinates). With epistasis > 0, genes also get additive-by-additive eQTL pairs, returned in epi_eqtl:

tx_epi <- simulate_transcriptome(geno, n_genes = 100, epistasis = 0.3, seed = 72)
head(tx_epi$epi_eqtl)
#>              gene_id        snp1        snp2 chr1 chr2      effect   prod_mean
#> ss196490285 gene0001 ss196490285 ss196494524    9    9  0.01483969  0.14667092
#> ss196515066 gene0001 ss196515066 ss196421592    2    1  0.21294751  0.22690051
#> ss196490571 gene0002 ss196490571 ss196470080    9    6  0.15622813  0.17396684
#> ss196463366 gene0003 ss196463366 ss196441917    5    2 -0.35730691  0.12364796
#> ss196426936 gene0004 ss196426936 ss196485784    1    8 -0.36233064 -0.02584184
#> ss196445264 gene0004 ss196445264 ss196443472    3    3  0.07602859 -0.01285714

A phenotype from expression

Hand the transcriptome to simulate_phenotype() with transcriptome = and add a transcriptome() layer. Its prop is required and is a separate variance category: it is not drawn from the h2 budget. Pick the causal genes with n_genes (random) or genes (explicit), and optionally fix their slopes:

ph_tx <- simulate_phenotype(geno, h2 = 0.5, seed = 71, transcriptome = tx) |>
  transcriptome(prop = 0.3, n_genes = 20) |>
  additive(prop = 0.2, n_qtn = 10)
ph_tx
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: geno   Traits: 1   Architecture: independent   Seed: 71
#>   Variance partition (proportions of V_P):
#>     transcriptome 0.30   (20 genes)
#>     additive    0.20   (10 QTNs, geometric)
#>     residual    0.50
#>   Requested genetic share = 0.20   realized H² = 0.21
#>   Expression-mediated (derived): genetic 0.06 + environmental 0.23 + cov 0.00 of V_P
#>   (the genetic-mediated share is included in realized H² above)

Passing transcriptome = TRUE derives the transcriptome from the same genome instead of supplying one.

Genetic and environmental mediation

For a derived transcriptome the component splits into a genetically mediated part (expression traced to the genome), which counts toward realized H², and an environmental part, which does not. mediation_split() reports both and their covariance as shares of the phenotypic variance:

mediation_split(ph_tx)
#>     trait genetic_mediated env_mediated  covariance
#> 1 Trait_1       0.05837516    0.2304993 0.002054562

qtn_table() lists the causal genes alongside any marker QTNs (layer tells them apart):

head(qtn_table(ph_tx), 4)
#>     trait         layer set      snp chr pos maf     effect  d var_explained
#> 1 Trait_1 transcriptome  NA gene0003  NA  NA  NA  0.1016006 NA  0.0002094589
#> 2 Trait_1 transcriptome  NA gene0020  NA  NA  NA  0.1562553 NA  0.0004954224
#> 3 Trait_1 transcriptome  NA gene0025  NA  NA  NA -0.8194778 NA  0.0136263812
#> 4 Trait_1 transcriptome  NA gene0036  NA  NA  NA -0.7116215 NA  0.0102755317
#>   QTN_t1 QTN_t2 ld_r2
#> 1   <NA>   <NA>    NA
#> 2   <NA>   <NA>    NA
#> 3   <NA>   <NA>    NA
#> 4   <NA>   <NA>    NA

A supplied expression matrix (expression =) is different: its genetic content is unknown, so the whole component is treated as environmental and is left out of the genetic value and of H².

Without genotypes

With an observed expression matrix and no genome, the phenotype is built from transcriptome() layers alone (there is no h2 or marker layer to set):

set.seed(1)
expr <- matrix(rnorm(50 * 100), 50, 100,
               dimnames = list(paste0("g", 1:50), paste0("i", 1:100)))
simulate_phenotype(expression = expr, seed = 75) |>
  transcriptome(prop = 0.4, n_genes = 5)
#> <phenotype_sim>  (realized · long format)
#>   Genotypes: <expression>   Traits: 1   Architecture: independent   Seed: 75
#>   Variance partition (proportions of V_P):
#>     transcriptome 0.40   (5 genes)
#>     residual    0.60
#>   Requested genetic share = 0.00   realized H² = 0.00

Reusing an architecture in a new population

predict() applies a fixed transcriptome architecture (the same eQTL effects) to new genotypes, so a model learned in one population can be evaluated in another:

ids   <- seq_len(ncol(geno) - 5L)
train <- geno[, c(1:5, 5L + ids[ids %% 2 == 1])]
test  <- geno[, c(1:5, 5L + ids[ids %% 2 == 0])]
tx_train <- simulate_transcriptome(train, n_genes = 100, seed = 73)
tx_test  <- predict(tx_train, test, seed = 74)
dim(tx_test$expression)
#> [1] 100 140

Calibrating to observed expression

Give mimic = a genes x individuals matrix (columns named as the genotype individuals) to calibrate the generator to it: per-gene mean and variance, a GREML heritability per gene, and the co-expression strength. The estimates are returned in $calibration. Effects and loadings are still drawn de novo, so this matches the distribution of the data, not individual eQTL.

# here the "observed" matrix is the simulated one above; use your own in practice
tx_m <- simulate_transcriptome(geno, mimic = tx$expression, seed = 77)
names(tx_m$calibration)
#> [1] "source"    "n_factors" "kappa"     "h2"

RNA-seq counts

observe_counts() adds a negative-binomial count layer on top of a transcriptome_sim (library size, baseline, coupling to expression and dispersion per gene):

tx_c <- observe_counts(tx, seed = 76)
tx_c$counts[1:4, 1:4]
#>          4226 4722 33-16 38-11
#> gene0001   35   35    18    14
#> gene0002   12   36    22    23
#> gene0003   26   39    19    12
#> gene0004   12   18    17    25

Directed regulatory networks and tissue specificity are outside this version’s scope. See docs/SPEC-transcriptome.md for the full design.

9. Original implementation of create_phenotypes

create_phenotypes() is the original interface, kept frozen so published results stay reproducible. It receives bug fixes only; new work should use the grammar above.

Its output always lands in a folder of its own, so a directory is never littered with loose files:

legacy_home <- file.path(tempdir(), "sp_legacy")
dir.create(legacy_home, showWarnings = FALSE)

invisible(create_phenotypes(
  geno_obj = geno, add_QTN_num = 3, add_effect = 0.2,
  rep = 2, h2 = 0.5, model = "A",
  home_dir = legacy_home, verbose = FALSE
))
#> Simulation completed!
#> Results are saved at:/tmp/RtmpPuyO0l/sp_legacy/simplePHENOTYPES_output

list.files(legacy_home)                                  # one folder
#> [1] "simplePHENOTYPES_output"
list.files(file.path(legacy_home, "simplePHENOTYPES_output"))
#> [1] "Additive_QTNs.txt"                   "Genetic_values.txt"                 
#> [3] "Log_Sim.txt"                         "Simulated_Data_2_Reps_Herit_0.5.txt"

Naming the folder yourself, and re-running without overwriting:

invisible(create_phenotypes(
  geno_obj = geno, add_QTN_num = 3, add_effect = 0.2,
  rep = 1, h2 = 0.5, model = "A",
  home_dir = legacy_home, output_dir = "my_run", verbose = FALSE
))
#> Simulation completed!
#> Results are saved at:/tmp/RtmpPuyO0l/sp_legacy/my_run
invisible(create_phenotypes(
  geno_obj = geno, add_QTN_num = 3, add_effect = 0.2,
  rep = 1, h2 = 0.5, model = "A",
  home_dir = legacy_home, output_dir = "my_run", verbose = FALSE
))
#> Directory name provided by 'output_dir' alredy exists! 
#> Creating: /tmp/RtmpPuyO0l/sp_legacy/my_run(1)
#> Simulation completed!
#> Results are saved at:/tmp/RtmpPuyO0l/sp_legacy/my_run(1)

list.files(legacy_home)      # my_run and my_run(1): nothing is overwritten
#> [1] "my_run"                  "my_run(1)"              
#> [3] "simplePHENOTYPES_output"

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

Complete reference: every option in one place