library(simplePHENOTYPES)
# **********
# Thank you for using the simplePHENOTYPES
# For the reference publication, please run: citation("simplePHENOTYPES")
# A Developmental version may be found at: https://github.com/samuelbfernandes/simplePHENOTYPES
# **********
data("SNP55K_maize282_maf04")
geno <- SNP55K_maize282_maf04simplePHENOTYPES can simulate meiosis directly, which means you can build a pedigree from real founder genotypes: cross two lines, self the offspring, produce doubled haploids, and then simulate phenotypes on the result.
These functions implement the algorithms from the isqg
package, reimplemented in Rust. A gamete carries on average one
crossover per Morgan, and chromosomes assort independently. By default
crossovers are spaced by the gamma model of crossover interference
(nu = 2.6, p = 0, AlphaSimR’s default; see
?cross). With interference = "poisson"
recombination follows isqg’s Karlin & Liberman count-location model:
the number of crossovers on a chromosome is Poisson with mean equal to
its length in Morgans and their positions are uniform along it.
There are four functions to learn:
| Function | What it does |
|---|---|
as_population() |
turns genotype data into phased founders |
cross() |
mates two individuals |
selfcross() |
self-pollinates one individual |
double_haploid() |
makes fully homozygous progeny in one step |
Meiosis is driven by genetic distance, in centiMorgans, not physical distance in base pairs. Two markers 1 Mb apart may recombine freely or almost never, depending on where they sit relative to the centromere.
The genetic map lives in the cm column of your genotype
data:
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.5956If yours is empty, as_population() will stop and tell
you so. Build one with synthetic_map():
The map bundled with this data set is synthetic — a model, not a
measured maize linkage map. See vignette("genetic-maps")
for how it is built, how to tune it, and how to substitute a real
map.
as_population() converts genotype data into a
Population: the phased representation meiosis needs. Pick
your founders with individuals.
pop <- as_population(geno, individuals = c("33-16", "38-11"))
pop
# <Population>
# Individuals: 2 Markers: 10650 Chromosomes: 10
# Genetic map: 1497 cM total span (109-219 cM per chromosome)
# Origin: founder
# Pedigree: 2 recorded individuals; generation 0
# IDs: 33-16, 38-11A dosage matrix says an individual is heterozygous but not
which homologue carries which allele, so heterozygotes are
phased arbitrarily. For an inbred panel like this one that is nearly
lossless, because heterozygotes are rare. For outbred material it is a
real assumption — see ?as_population.
n_individuals(pop)
# [1] 2
pop[1] # subset individuals
# <Population>
# Individuals: 1 Markers: 10650 Chromosomes: 10
# Genetic map: 1497 cM total span (109-219 cM per chromosome)
# Origin: founder
# Pedigree: 1 recorded individuals; generation 0
# IDs: 33-16
dim(dosages(pop)) # markers x individuals, coded -1/0/1
# [1] 10650 2
dosages(pop)[1:4, ]
# 33-16 38-11
# ss196422159 -1 -1
# ss196422171 -1 -1
# ss196422173 -1 -1
# ss196422186 -1 -1dosages() collapses the phased haplotypes back to the
familiar -1 / 0 / 1 coding: -1 and 1 are the
two homozygotes, 0 is a heterozygote.
cross() mates two individuals. Each parent must be a
single individual, so subset with [ if your
population holds more.
f1 <- cross(pop[1], pop[2], n = 1, seed = 11)
# 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:
# Toledo, F.H., Perez-Rodriguez, P., Crossa, J. and Burgueno, J. (2019). isqg: A Binary Framework for in Silico Quantitative Genetics. G3 9(8):2425-2428, doi:10.1534/g3.119.400373, when using the crossing pipelines (cross, selfcross, double_haploid).
# This message is displayed once per session.
f1
# <Population>
# Individuals: 1 Markers: 10650 Chromosomes: 10
# Genetic map: 1497 cM total span (109-219 cM per chromosome)
# Origin: cross(33-16 x 38-11)
# Pedigree: 3 recorded individuals; generation 1
# IDs: prog_1The markers worth watching are the ones that can tell the founders apart: both founders homozygous, and for different alleles. At those the F1 must be heterozygous, because it received one homologue from each parent.
d_pop <- dosages(pop)
informative <- which(d_pop[, 1] != d_pop[, 2] & d_pop[, 1] != 0 & d_pop[, 2] != 0)
length(informative)
# [1] 5157
all(dosages(f1)[informative, 1] == 0)
# [1] TRUERequiring both founders to be homozygous matters. A locus where one founder is heterozygous also “differs” between them, but the F1 there is heterozygous only half the time, and including such loci quietly blurs every ratio that follows.
Ask for more progeny with n. Each is an independent
draw, so full sibs differ:
sibs <- cross(pop[1], pop[2], n = 5, seed = 12)
sibs
# <Population>
# Individuals: 5 Markers: 10650 Chromosomes: 10
# Genetic map: 1497 cM total span (109-219 cM per chromosome)
# Origin: cross(33-16 x 38-11)
# Pedigree: 7 recorded individuals; generation 1
# IDs: prog_1, prog_2, prog_3, prog_4, prog_5selfcross() pollinates an individual with itself. Both
gametes come from independent meioses of the same plant, so an F2
segregates:
f2 <- selfcross(f1, n = 200, seed = 13)
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)At the informative markers the classic 1:2:1 ratio appears:
double_haploid() takes a single gamete and doubles it,
producing a completely homozygous individual in one generation. This is
the fast route to a fixed line:
A backcross is just a cross back to one of the parents — there is no separate function:
Every BC1 gets one gamete from the recurrent parent, so at an informative locus it is either homozygous for that parent or still heterozygous — never homozygous for the donor. Counting how many of the two alleles came from the recurrent parent:
recurrent <- dosages(pop[1])[informative, 1]
d_bc <- dosages(bc1)[informative, ]
n_recurrent <- ifelse(d_bc == recurrent, 2, ifelse(d_bc == 0, 1, 0))
round(prop.table(table(n_recurrent)), 3)
# n_recurrent
# 1 2
# 0.497 0.503
mean(n_recurrent) / 2 # fraction of the genome from the recurrent parent
# [1] 0.7512643Half the loci are fixed for the recurrent parent and half are heterozygous, which works out to 75% of alleles — the familiar backcross arithmetic (50% from its own gamete, plus half of the F1’s 50%).
Repeated selfing halves heterozygosity each generation. Single-seed descent is a loop — advance each line by selfing one plant per generation:
het_of <- function(x) mean(dosages(x)[informative, 1] == 0)
set.seed(7)
lines <- vector("list", 20) # 20 independent F2 lines
for (i in seq_along(lines)) lines[[i]] <- selfcross(f1, n = 1)
het <- numeric(6)
het_all <- function(pop_list) {
out <- numeric(length(pop_list))
for (i in seq_along(pop_list)) out[i] <- het_of(pop_list[[i]])
out
}
het[1] <- mean(het_all(lines))
for (g in 2:6) {
for (i in seq_along(lines)) lines[[i]] <- selfcross(lines[[i]], n = 1)
het[g] <- mean(het_all(lines))
}
data.frame(
generation = paste0("F", 2:7),
heterozygosity = round(het, 3),
expected = round(0.5^(1:6), 3)
)
# generation heterozygosity expected
# 1 F2 0.489 0.500
# 2 F3 0.255 0.250
# 3 F4 0.117 0.125
# 4 F5 0.060 0.062
# 5 F6 0.029 0.031
# 6 F7 0.013 0.016Averaged across lines, heterozygosity halves each generation, as
theory says it should. After six generations the lines are essentially
fixed — the same endpoint double_haploid() reaches in one
step, which is why doubled haploids are attractive when the technology
exists for your crop.
Note the averaging. Any single line is a poor estimate of the expectation, because heterozygosity is not scattered locus by locus — it survives in long chromosomal blocks, so a genome carries far fewer independent pieces of information than it has markers:
At F7 individual lines range from essentially fixed to several percent heterozygous. That spread is real biology, not simulation noise, and it is why breeding programmes carry many lines rather than one.
Segregation ratios confirm the alleles are being transmitted correctly, but they say nothing about linkage. To check recombination you have to ask which parent each locus came from.
Doubled haploids make this easy: they are homozygous, so their allele names the founder it came from directly.
on_chr1 <- geno$chr == 1
p1 <- dosages(pop[1])[on_chr1, 1]
p2 <- dosages(pop[2])[on_chr1, 1]
# Informative loci: both founders homozygous, for different alleles.
inf1 <- which(p1 != p2 & p1 != 0 & p2 != 0)
cm1 <- geno$cm[on_chr1][inf1]
from_p1 <- dosages(dh)[on_chr1, ][inf1, ] == p1[inf1]
pick <- c(5, 20, 100, 400)
data.frame(
distance_cM = round(cm1[pick] - cm1[1], 1),
recombination = round(rowMeans(from_p1[pick, ] != from_p1[rep(1, 4), ]), 3),
haldane = round((1 - exp(-2 * (cm1[pick] - cm1[1]) / 100)) / 2, 3)
)
# distance_cM recombination haldane
# ss196422258 0.9 0.020 0.009
# ss196526252 4.0 0.085 0.039
# ss196418562 27.3 0.250 0.210
# ss196426535 90.2 0.445 0.418The observed recombination fractions track Haldane’s map function, \((1 - e^{-2d})/2\) at distance \(d\) Morgans. Nearby loci travel together, distant ones approach one half — and never exceed it, because a recombination fraction above 0.5 is impossible.
This check is worth remembering: it is sensitive to bugs that segregation ratios cannot see. An earlier version of this package reconstructed progeny phase from dosages, which silently produced haplotypes the founders never had. F2 and DH ratios still looked perfect; the recombination fractions came out at 0.99, which is what exposed it.
All randomness is drawn by R, so a pedigree is exactly reproducible.
Either pass seed, or set one beforehand:
A Population is accepted anywhere
simulate_phenotype() takes genotype data, so you can go
straight from pedigree to phenotype:
ph <- simulate_phenotype(f2, seed = 21) |>
additive(prop = 0.6, n_qtn = 5)
ph
# <phenotype_sim> (realized · long format)
# Genotypes: <Population: selfcross(prog_1)> Traits: 1 Architecture: independent Seed: 21
# Variance partition (proportions of V_P):
# additive 0.60 (5 QTNs, geometric)
# residual 0.40
# Requested genetic share = 0.60 realized H² = 0.58
head(phenotypes_long(ph))
# id trait rep value
# 1 self_1 Trait_1 1 -0.800873783
# 2 self_2 Trait_1 1 -1.507204033
# 3 self_3 Trait_1 1 0.296709284
# 4 self_4 Trait_1 1 -2.189896227
# 5 self_5 Trait_1 1 0.473239995
# 6 self_6 Trait_1 1 -0.005319188This is the natural way to simulate a mapping population: build the
pedigree, then lay a genetic architecture over it. See
vignette("simplePHENOTYPES-v2") for the full grammar.
Putting it together — a biparental doubled-haploid mapping population with a simulated trait:
founders <- as_population(geno, individuals = c("33-16", "38-11"))
f1_line <- cross(founders[1], founders[2], n = 1, seed = 1)
dh_pop <- double_haploid(f1_line, n = 150, seed = 2)
trait <- simulate_phenotype(dh_pop, seed = 3) |>
additive(prop = 0.5, n_qtn = 10)
n_individuals(dh_pop)
# [1] 150
head(phenotypes_wide(trait))
# id rep Trait_1
# 1 dh_1 1 -1.0933637
# 2 dh_2 1 0.6667039
# 3 dh_3 1 0.2931373
# 4 dh_4 1 -0.8436566
# 5 dh_5 1 -1.3086055
# 6 dh_6 1 0.1968329“mother must contain exactly one
individual” — the mating functions take one individual at a
time. Subset with [: cross(pop[1], pop[2]),
not cross(pop, pop).
“The genetic map (cm) is all NA” —
meiosis needs recombination distances. Build a map with
synthetic_map() first.
“cm must be non-decreasing within each
chromosome” — sort your markers by chromosome and position
before converting.
Progeny look unexpectedly similar to one parent — check how many markers are actually informative. Two founders that share an allele at a locus produce progeny that are uninformative there no matter how meiosis went.
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 htmltools_0.5.9 buildtools_1.0.0
# [10] lifecycle_1.0.5 cli_3.6.6 sass_0.4.10
# [13] data.table_1.18.6.1 jquerylib_0.1.4 compiler_4.6.1
# [16] sys_3.4.3 tools_4.6.1 evaluate_1.0.5
# [19] bslib_0.12.0 yaml_2.3.12 otel_0.2.0
# [22] jsonlite_2.0.0 rlang_1.3.0