What tidybreed is
tidybreed simulates animal and plant breeding programs.
You build a population step by step — a genome, founder animals, traits,
phenotypes, matings, selection — and the package tracks every
individual, haplotype, and record for you.
What makes it different from other simulators is where the data lives. Everything is stored in a DuckDB database, not in R memory. A simulation with a million animals does not need a million animals’ worth of RAM. The database is a plain file on disk, so runs are resumable, shareable, and queryable with ordinary SQL long after the R session that made them has ended.
The API is pipe-friendly and built on dplyr verbs, so if
you know the tidyverse you already know most of the syntax.
This vignette walks through a complete single-generation simulation. Every code block below actually runs when this page is built, so what you see is real output.
Installation
tidybreed is not on CRAN yet. Install from GitHub —
pak handles the compiled code and dependencies cleanly:
install.packages("pak")
pak::pak("austin-putz/tidybreed")You will need a C++ compiler (Xcode command line tools on macOS, Rtools on Windows). Then:
How the package works
Three ideas explain almost all of the API. They are worth reading before the code.
1. R is the control layer, DuckDB is the storage layer
Your R session holds a small object with a database connection in it.
It does not hold your genotypes. When you ask for a
table you get a lazy reference — nothing is pulled into R until
you call collect(). This is why you should filter first and
collect last.
2. There are two object types
| Object | Created by | What it is |
|---|---|---|
tidybreed_pop |
open_pop() |
The population: a database connection plus a table registry |
tidybreed_table |
get_table() |
A lazy reference to one table, plus any pending filter |
Most functions take a tidybreed_pop and return a
tidybreed_pop, so they chain:
pop <- pop |> define_genome(...) |> define_founder_haplotypes(...)But action functions — the ones that operate on a
subset of rows — take a tidybreed_table as their
first argument instead. That is how you tell them which rows to act
on:
pop <- pop |>
get_table("ind_meta") |> # which table
filter(sex == "F") |> # which rows
add_phenotype("ADG") # what to doadd_founders(), add_phenotype(),
add_tgv(), add_index(),
define_additive_effects(), define_chip(),
add_genotypes(), extract_genotypes(), and
mutate_table() all work this way. Most of them return a
tidybreed_pop, so the chain continues; the
extract_* family is the exception — it hands you data back
as a tibble, since pulling results into R is the whole point of it.
There is no positional row selection anywhere in
tidybreed — no TRUE/FALSE masks,
no integer indices. Database row order is not guaranteed, so selection
always goes through names (id_ind, locus_name)
or a SQL filter.
3. Genetics and observations are separate layers
This one trips people up, so it is worth being explicit. A trait’s genetics and the phenotype someone writes on a clipboard are defined by two different functions:
| You want to specify | Function | Examples |
|---|---|---|
| Genetic component | define_trait() |
description, units
|
| Genetic target and QTL effects | define_additive_effects() |
G, anchor, method,
base_tbl
|
| What gets observed | define_phenotype() |
mean, type, residual_var,
expressed_sex, repeatable
|
For a simple trait the two share a name and you call both. The payoff comes later: one observed phenotype (weaning weight) can be built from several genetic components (direct + maternal), and the split makes that natural instead of a special case.
Open a population
open_pop() creates the database. Here we use
":memory:" so this vignette writes nothing to disk; in a
real run you would omit db_name and get a
sim.duckdb file inside an organised output folder.
pop <- open_pop(pop_name = "demo", db_name = ":memory:")
#> duckdb keeps downloaded extensions and secrets in a temporary directory:
#> ℹ /tmp/RtmpcLi7Ca/duckdb
#> This is removed when the R session ends.
#> • Extensions are re-downloaded each session.
#> • Secrets are lost.
#> ℹ Run duckdb(shared_home = TRUE) (or create ~/.duckdb) to keep them (suitable for most users).
#> ℹ Run duckdb(shared_home = FALSE) to accept the temporary directory (and silence this message).
#> ℹ See ?duckdb_storage for details and alternatives.
#> Opened in-memory population 'demo'
pop
#> ── tidybreed population: demo ──────────────────────────────────────────────────
#> Database in-memory [connected]
#>
#> schema(pop) · describe_table(pop, "name")
#> ────────────────────────────────────────────────────────────────────────────────Printing a population gives you a live summary. Right now there is almost nothing in it — sections appear as you add data.
Define the genome
pop <- pop |>
define_genome(
n_loci = 500, # total loci (SNP and QTL both come from this pool)
n_chr = 5, # chromosomes
chr_len_Mb = 100, # physical length of each chromosome in megabases
cM_per_Mb = 1.0 # genetic-map rate (the default)
)
#> Defined genome: 500 loci across 5 chromosomes | chr lengths (Mb): all equal to 100 Mb
#> Tables written: genome_meta, genome_map, ind_haplotype, ind_genotype, ind_crossover, chr_inheritance, chr_recombination, genome_effects, genome_effect_members, genome_effect_member_origins
#> Views created: genome_effect_terms, genome_effect_loci
pop |> get_table("genome_meta")
#> <tidybreed_table: genome_meta> [500 rows × 5 fields]
#> # A tibble: 10 × 5
#> locus_id locus_name chr chr_name pos_bp
#> <int> <chr> <int> <chr> <dbl>
#> 1 1 Locus_1 1 1 990099
#> 2 2 Locus_2 1 1 1980198
#> 3 3 Locus_3 1 1 2970297
#> 4 4 Locus_4 1 1 3960396
#> 5 5 Locus_5 1 1 4950495
#> 6 6 Locus_6 1 1 5940594
#> 7 7 Locus_7 1 1 6930693
#> 8 8 Locus_8 1 1 7920792
#> 9 9 Locus_9 1 1 8910891
#> 10 10 Locus_10 1 1 9900990genome_meta holds one row per locus with its
physical position (pos_bp). The
genetic map lives in its own table, because a genome
can have several maps — male vs female, or line-specific:
pop |> get_table("genome_map")
#> <tidybreed_table: genome_map> [500 rows × 7 fields]
#> # A tibble: 10 × 7
#> id_genome_map locus_id locus_name sex line_name map_name pos_cM
#> <int> <int> <chr> <chr> <chr> <chr> <dbl>
#> 1 1 1 Locus_1 NA NA default 0.990
#> 2 2 2 Locus_2 NA NA default 1.98
#> 3 3 3 Locus_3 NA NA default 2.97
#> 4 4 4 Locus_4 NA NA default 3.96
#> 5 5 5 Locus_5 NA NA default 4.95
#> 6 6 6 Locus_6 NA NA default 5.94
#> 7 7 7 Locus_7 NA NA default 6.93
#> 8 8 8 Locus_8 NA NA default 7.92
#> 9 9 9 Locus_9 NA NA default 8.91
#> 10 10 10 Locus_10 NA NA default 9.90cM_per_Mb is what turns physical distance into genetic
distance (pos_cM = pos_bp / 1e6 * cM_per_Mb), and it is
written here as the default map — sex and
line_name both NULL,
map_name = "default". Everything distance-driven reads that
map, so this one number sets the crossover rate: a 100 Mb chromosome at
cM_per_Mb = 1.0 is 100 cM, or ~1 crossover per meiosis.
define_genome() is a one-shot call. It writes seven
tables in a single transaction and rolls all of them back if anything
fails, so a bad call leaves the database untouched and you can simply
call it again with corrected arguments. Calling it on a population that
already has a genome is an error.
Per-chromosome rules live in two explicit tables:
chr_inheritance (how many copies of a chromosome an
offspring of each sex inherits from each parent) and
chr_recombination (whether a parent of each sex recombines
it). Every chromosome starts as a normal diploid autosome;
define_chromosome() is how you would declare an X, Y, or
mitochondrial chromosome.
pop |> get_table("chr_inheritance") |> collect()
#> # A tibble: 5 × 5
#> chr_name offspring_sex line_name from_parent_1 from_parent_2
#> <chr> <chr> <chr> <int> <int>
#> 1 1 NA NA 1 1
#> 2 2 NA NA 1 1
#> 3 3 NA NA 1 1
#> 4 4 NA NA 1 1
#> 5 5 NA NA 1 1
pop |> get_table("chr_recombination") |> collect()
#> # A tibble: 5 × 4
#> chr_name parent_sex line_name recombines
#> <chr> <chr> <chr> <lgl>
#> 1 1 NA NA TRUE
#> 2 2 NA NA TRUE
#> 3 3 NA NA TRUE
#> 4 4 NA NA TRUE
#> 5 5 NA NA TRUEExplore the database
Two helpers exist so you never have to guess what is in the database.
schema(pop) lists every table:
schema(pop)
#> ── Schema: demo ─────────────────────────────── 27 tables · 3.3 MiB in memory ──
#> Use describe_table(pop, "name") for column-level details.
#>
#> Genome
#> genome_meta 500 5 Locus-level metadata. One row p...
#> genome_map 500 7 Genetic map in long format. One...
#> chr_inheritance 5 5 Per-chromosome copy counts, key...
#> chr_recombination 5 4 Per-chromosome recombination, k...
#> + 5 empty: genome_effects, genome_effect_members,
#> genome_effect_member_origins, genome_effect_terms,
#> genome_effect_loci
#>
#> Individuals
#> + 4 empty: ind_meta, ind_haplotype, ind_genotype, ind_crossover
#>
#> Genetic model
#> + 2 empty: trait_meta, trait_var_comp
#>
#> Observation model
#> + 5 empty: phenotype_meta, phenotype_components, phenotype_var_comp,
#> phenotype_effects, phenotype_random_effects
#>
#> Selection
#> + 1 empty: index_meta
#>
#> Results
#> + 6 empty: ind_tgv, ind_tgv_total, ind_phenotype, ind_ebv, ind_index,
#> ind_true_indexdescribe_table() explains one table’s columns:
describe_table(pop, "ind_meta")
#> ── ind_meta ─────────────────────────────────────────────────── 0 rows · 6 cols
#> Individual-level metadata. One row per individual. Core columns are managed by the system; user-defined columns can be added via mutate_table() or the ... arguments of add_founders() and add_offspring().
#>
#> Column Type Description
#> ──────────────────────────────────────────────────────────────────────────────
#> id_ind VARCHAR Primary key; format '{line_name}_{n}' (e.g. 'A_1',...
#> id_parent_1 VARCHAR Paternal parent id_ind; NA for founders
#> id_parent_2 VARCHAR Maternal parent id_ind; NA for founders
#> line_name VARCHAR Genetic line name (e.g. 'A', 'Holstein')
#> sex VARCHAR Sex of the individual: 'M' for male, 'F' for female
#> ploidy UTINYINT Genome ploidy; declared at add_founders() time (mu...Remember that get_table() is lazy. The filter below runs
inside DuckDB and only the answer comes back to R:
Create founder haplotypes
Founders are sampled from a pool of haplotypes. Tagging the pool with
a line_name is good practice — it is what lets you run
crossbreeding later.
pop <- pop |>
define_founder_haplotypes(
line_name = "A",
n_haplotypes = 100,
method = "uniform",
min_allele_freq = 0.05,
max_allele_freq = 0.95
)
#> Generating 100 founder haplotypes (method = "uniform")...
#> Created founder_haplotypes for line 'A' (100 haplotypes x 500 loci)Several methods control the allele-frequency spectrum and linkage disequilibrium:
method |
What it gives you |
|---|---|
"uniform" |
Frequencies drawn uniformly between min_allele_freq and
max_allele_freq
|
"fixed" |
Every locus at the same allele_freq
|
"beta" |
Beta-distributed frequencies (beta_shape1,
beta_shape2) — U-shaped is realistic |
"balding_nichols" |
FST-based divergence between lines (fst,
mean_allele_freq) |
"mosaic" |
Haplotypes built from templates, producing simple LD
(n_templates, template_switch_rate) |
"gaussian_copula" |
Correlated frequencies with distance-based decay
(ld_decay_rate) |
Each argument belongs to exactly one method — passing one to the
wrong method is an error that names where it belongs. The two LD methods
walk the genetic map you defined above, resolved for this pool’s own
line_name.
The first four methods also take
exact_freq. By default they draw a
frequency p per locus and then sample each allele as an
independent Bernoulli(p), so the realized pool
frequency scatters around p. With
exact_freq = TRUE each locus instead gets exactly
round(p * n_haplotypes) copies of the 1-allele, on an
independently drawn subset of haplotypes — so the realized frequency
is p, and still no LD is induced. It
defaults to TRUE for "fixed" (a fixed
frequency that drifts is not fixed) and FALSE elsewhere;
turn it on when you want a drift-free base frequency.
One caveat worth knowing early: in "mosaic", templates
are the only source of variation, so n_templates sets the
MAF spectrum, not just the block length. Roughly
2 / (n_templates + 1) of loci come out monomorphic no
matter how many haplotypes you ask for, and you get a warning past 10% —
QTL landing there contribute no genetic variance. Raise
n_templates, or use "gaussian_copula", which
gives an unquantized spectrum.
Add founder individuals
This is the first action function, and it shows why the
tidybreed_table pattern exists: you pipe in the haplotype
pool you want to sample from.
Any extra named argument — like gen below — becomes a
new column on ind_meta, written in the same transaction as
the animals themselves.
pop <- pop |>
get_table("founder_haplotypes") |>
filter(line_name == "A") |>
add_founders(
n_males = 250,
n_females = 250,
line_name = "A",
gen = 0L # custom column
)
#> Added new column 'gen' (INTEGER) to `ind_meta`
#> Added 500 founders (250 males, 250 females) to line 'A'
pop |> get_table("ind_meta") |> collect() |> head()
#> # A tibble: 6 × 7
#> id_ind id_parent_1 id_parent_2 line_name sex ploidy gen
#> <chr> <chr> <chr> <chr> <chr> <int> <int>
#> 1 A_1 NA NA A M 2 0
#> 2 A_2 NA NA A M 2 0
#> 3 A_3 NA NA A M 2 0
#> 4 A_4 NA NA A M 2 0
#> 5 A_5 NA NA A M 2 0
#> 6 A_6 NA NA A M 2 0Watch the type suffix.
gen = 0Lcreates anINTEGERcolumn;gen = 0would create aDOUBLE. The same applies toNA_integer_,NA_character_, and friends.
Haplotypes are stored in long format — one row per individual × parent × locus:
pop |> get_table("ind_haplotype")
#> <tidybreed_table: ind_haplotype> [5e+05 rows × 7 fields]
#> # A tibble: 10 × 7
#> id_ind parent_origin strand line_origin locus_id locus_name allele
#> <chr> <int> <int> <chr> <int> <chr> <int>
#> 1 A_1 1 1 A 1 Locus_1 1
#> 2 A_2 1 1 A 1 Locus_1 0
#> 3 A_3 1 1 A 1 Locus_1 0
#> 4 A_4 1 1 A 1 Locus_1 0
#> 5 A_5 1 1 A 1 Locus_1 1
#> 6 A_6 1 1 A 1 Locus_1 0
#> 7 A_7 1 1 A 1 Locus_1 0
#> 8 A_8 1 1 A 1 Locus_1 0
#> 9 A_9 1 1 A 1 Locus_1 1
#> 10 A_10 1 1 A 1 Locus_1 1Add your own columns
mutate_table() adds or updates columns on any
table. This is the mechanism for keeping simulation state in the
database instead of in parallel R objects.
# Broadcast a value to every row
pop <- pop |> get_table("ind_meta") |> mutate_table(farm = "Iowa")
#> Added new column 'farm' (VARCHAR) to `ind_meta`; 500 rows set
# Update only some rows
pop <- pop |>
get_table("ind_meta") |>
filter(sex == "M") |>
mutate_table(farm = "AI_Stud")
#> Warning: 'farm': replaced 250 existing values in `ind_meta` [250 of 500 rows]
pop |> get_table("ind_meta") |> count(sex, farm) |> collect()
#> # A tibble: 2 × 3
#> sex farm n
#> <chr> <chr> <dbl>
#> 1 F Iowa 250
#> 2 M AI_Stud 250The warning above is deliberate: mutate_table() tells
you whenever it overwrites values that were already there, so you notice
an accidental clobber instead of discovering it three generations
later.
Calling mutate_table() on an empty table still
creates the column, which lets you declare a typed schema before any
data arrives:
pop |> get_table("ind_ebv") |> mutate_table(model_version = NA_character_)Define a trait
The genetic layer. A trait is just a name here; its additive genetic variance target comes with its QTL effects, below.
pop <- pop |>
define_trait("ADG", units = "kg/day")
#> Added trait 'ADG'.
pop |> get_table("trait_meta") |> collect()
#> # A tibble: 1 × 4
#> id_trait trait_name description units
#> <int> <chr> <chr> <chr>
#> 1 1 ADG NA kg/dayNote what is not here: no mean, no residual variance, no trait type. Those are observation-layer properties and come later.
Assign QTL effects
Which loci are QTL is decided by the filter you pipe in. Here
chromosomes 4 and 5 carry the QTL, leaving 1–3 free for a SNP chip.
G = 0.25 is the additive genetic variance target: it is
stored in trait_var_comp, and the sampled effects are
calibrated so they deliver it exactly at the founder pool’s allele
frequencies.
pop <- pop |>
get_table("genome_meta") |>
filter(chr %in% c(4L, 5L)) |>
define_additive_effects("ADG", G = 0.25)
#> The founder pool (pool expectation under random pairing) sees relative spectrum [0.939, 0.939] of the target (sampling LD of 100 haplotypes).
#> Set additive effects for 200 QTL on trait 'ADG' (base: founder_haplotypes; scope: all lines, both parents' copies): exact under the genic anchor; delivered variance 0.25.
pop |> get_table("genome_effect_terms") |> collect() |> head()
#> # A tibble: 6 × 9
#> trait_name effect_owner effect_name id_genome_effect effect_order
#> <chr> <chr> <chr> <int> <dbl>
#> 1 ADG generated NA 1 1
#> 2 ADG generated NA 2 1
#> 3 ADG generated NA 3 1
#> 4 ADG generated NA 4 1
#> 5 ADG generated NA 5 1
#> 6 ADG generated NA 6 1
#> # ℹ 4 more variables: contrast_signature <chr>, family_key <chr>,
#> # scope_description <chr>, genome_value <dbl>Causal-locus membership is implicit: a locus is
causal for a trait if it appears as a member of any term for it —
genome_effect_loci is the view that answers that. There is
no is_QTL flag to keep in sync. The centring frequency
lives on the member as center_value, because a
line-specific variant legitimately centres against its own line’s base
frequency.
Define the phenotype
The observation layer. With an additive target of 0.25 and
residual_var = 0.75, this trait has a heritability of
0.25.
pop <- pop |>
define_phenotype(
"ADG",
type = "continuous",
mean = 1.0,
residual_var = 0.75
)
#> Added phenotype 'ADG' (type: continuous).
pop |> get_table("phenotype_meta") |> collect() |> select(1:6)
#> # A tibble: 1 × 6
#> id_phenotype_meta phenotype_name type mean expressed_sex repeatable
#> <int> <chr> <chr> <dbl> <chr> <lgl>
#> 1 1 ADG continuous 1 both FALSEtype also accepts "count",
"categorical" (with prevalence or
thresholds), and "derived_formula".
Add a fixed effect
Non-genetic terms attach to the phenotype, not the trait. Here males grow 0.30 kg/day faster than females, with females as the reference level.
pop <- pop |>
define_effect_fixed_class(
"ADG",
effect_name = "sex",
source_column = "sex",
levels = c(M = 0.30, F = 0)
)
#> Added fixed-class effect 'sex' to phenotype 'ADG' (2 levels).define_effect_fixed_cov() adds a regression on a
covariate, and define_effect_random() adds a named random
effect such as pen or litter.
Compute breeding values and phenotypes
add_tgv() computes true genetic values — the
simulation’s ground truth, which you would never know in a real
population. It writes ind_tgv, one row per individual,
trait and component; the breeding value is the
additive component (the only one in this purely additive
model). The total over every component is the ind_tgv_total
view.
pop <- pop |> get_table("ind_meta") |> add_tgv("ADG")
#> Computed TGV for 500 individuals on trait 'ADG' (additive).
pop |> get_table("ind_tgv") |> collect() |> head()
#> # A tibble: 6 × 5
#> id_tgv id_ind trait_name component_name tgv_value
#> <int> <chr> <chr> <chr> <dbl>
#> 1 6 A_103 ADG additive -0.457
#> 2 8 A_105 ADG additive -0.206
#> 3 12 A_109 ADG additive -0.178
#> 4 13 A_11 ADG additive 0.130
#> 5 16 A_112 ADG additive -0.153
#> 6 20 A_116 ADG additive 0.0310add_phenotype() builds the observed record: intercept,
plus fixed and random effects, plus the total genetic value, plus a
sampled residual. It calls add_tgv() internally, so you can
go straight to it.
pop <- pop |> get_table("ind_meta") |> add_phenotype("ADG")
#> Computed TGV for 500 individuals on trait 'ADG' (additive).
#> Wrote 500 phenotype records for 'ADG'.
pop |> get_table("ind_phenotype") |> collect() |> head()
#> # A tibble: 6 × 9
#> id_phenotype id_ind phenotype_name pheno_value pheno_number liability_value
#> <int> <chr> <chr> <dbl> <int> <dbl>
#> 1 1 A_1 ADG 1.71 1 NA
#> 2 2 A_10 ADG 1.00 1 NA
#> 3 3 A_100 ADG 1.83 1 NA
#> 4 4 A_101 ADG 2.07 1 NA
#> 5 5 A_102 ADG 1.90 1 NA
#> 6 6 A_103 ADG 2.45 1 NA
#> # ℹ 3 more variables: cat_name <chr>, residual_value <dbl>,
#> # residual_condition_level <chr>The sex effect we defined should show up as a gap of roughly 0.30 between the group means:
pop |>
get_table("ind_phenotype") |>
collect() |>
left_join(
pop |> get_table("ind_meta") |> collect() |> select(id_ind, sex),
by = "id_ind"
) |>
group_by(sex) |>
summarise(n = n(), mean_ADG = mean(pheno_value))
#> # A tibble: 2 × 3
#> sex n mean_ADG
#> <chr> <int> <dbl>
#> 1 F 250 1.02
#> 2 M 250 1.34A second, genetically correlated trait
Real programs select on several traits at once. Define the second trait, then supply a genetic covariance matrix for both.
ADG already has a stored 1 x 1 target. A stored target is never overwritten (not even by the same number), so the ADG-only block is removed first, then the 2 x 2 block is written together with the new effects:
pop <- pop |> define_trait("BF", units = "mm")
#> Added trait 'BF'.
G <- matrix(
c(0.25, 0.08,
0.08, 0.30),
nrow = 2,
dimnames = list(c("ADG", "BF"), c("ADG", "BF"))
)
pop <- pop |>
get_table("trait_var_comp") |>
filter(effect_name == "additive", trait_name_1 == "ADG") |>
remove_rows()
#> Deleted 1 row from `trait_var_comp`
pop <- pop |>
get_table("genome_meta") |>
filter(chr %in% c(4L, 5L)) |>
define_additive_effects(c("ADG", "BF"), G = G)
#> The founder pool (pool expectation under random pairing) sees relative spectrum [0.836, 0.894] of the target (sampling LD of 100 haplotypes).
#> Set correlated additive effects for traits: ADG, BF (method: shared; base: founder_haplotypes; scope: all lines, both parents' copies): exact under the genic anchor; delivered [ADG,ADG = 0.25; ADG,BF = 0.08; BF,BF = 0.3].
pop |> get_table("trait_var_comp") |> collect()
#> # A tibble: 4 × 6
#> id_trait_var_comp effect_name line_name trait_name_1 trait_name_2 cov_value
#> <int> <chr> <chr> <chr> <chr> <dbl>
#> 1 1 additive NA ADG ADG 0.25
#> 2 2 additive NA ADG BF 0.08
#> 3 3 additive NA BF ADG 0.08
#> 4 4 additive NA BF BF 0.3Covariance matrices are stored in the database, not passed around in
R. The one call drew correlated effects for both traits from
MVN(0, G) and then calibrated them, so variances
and the genetic correlation are exactly those of
G at the founder pool’s allele frequencies. BF then gets
its observation layer:
pop <- pop |>
define_phenotype("BF", type = "continuous", mean = 12, residual_var = 0.70)
#> Added phenotype 'BF' (type: continuous).This replaces the ADG effects assigned earlier with a fresh correlated draw. Re-running
define_additive_effects()replaces the variant at the same scope — here the common, all-lines one — and leaves variants at other scopes (line_name,parent_origin) standing. Recompute anything derived from the replaced terms.
pop <- pop |> get_table("ind_meta") |> add_tgv() # no trait_name = all traits
#> Computed TGV for 500 individuals on trait 'ADG' (additive).
#> Computed TGV for 500 individuals on trait 'BF' (additive).
pop |>
get_table("ind_tgv") |>
filter(component_name == "additive") |>
collect() |>
group_by(trait_name) |>
summarise(n = n(), mean_tbv = mean(tgv_value), var_tbv = var(tgv_value))
#> # A tibble: 2 × 4
#> trait_name n mean_tbv var_tbv
#> <chr> <int> <dbl> <dbl>
#> 1 ADG 500 -0.0280 0.211
#> 2 BF 500 -0.00194 0.268The realised variances are in the neighbourhood of the 0.25 and 0.30
we asked for, but they are not exact — and that is worth understanding
rather than glossing over. The calibration is exact for the
genic anchor: the random-mating limit at the founder pool’s
allele frequencies, which a population approaches over generations of
random mating. The animals you actually created are a finite sample from
that pool, with their own linkage disequilibrium, so their realised
variance wobbles around the target, and
define_additive_effects() warns when the pool expectation
is far off. Larger founder populations wobble less. For a
single-generation study, anchor = "realised" makes the
covariance exact in the named individuals themselves.
Make some offspring
Matings are described by an ordinary tibble: one row per offspring,
with id_parent_1, id_parent_2,
sex, and line_name required. Any custom
ind_meta column can ride along.
sires <- pop |> get_table("ind_meta") |> filter(sex == "M") |> collect() |> pull(id_ind)
dams <- pop |> get_table("ind_meta") |> filter(sex == "F") |> collect() |> pull(id_ind)
matings <- tibble::tibble(
id_parent_1 = rep(sires[1:5], each = 4),
id_parent_2 = dams[1:20],
sex = rep(c("M", "F"), 10),
line_name = "A",
gen = 1L
)
head(matings)
#> # A tibble: 6 × 5
#> id_parent_1 id_parent_2 sex line_name gen
#> <chr> <chr> <chr> <chr> <int>
#> 1 A_1 A_251 M A 1
#> 2 A_1 A_252 F A 1
#> 3 A_1 A_253 M A 1
#> 4 A_1 A_254 F A 1
#> 5 A_2 A_255 M A 1
#> 6 A_2 A_256 F A 1add_offspring() simulates meiosis — crossovers are drawn
from the genetic map in genome_map, respecting the
per-chromosome rules in chr_inheritance and
chr_recombination.
pop <- add_offspring(pop, matings)
#> Added 20 offspring (base_seed = 1026959721)
pop |> get_table("ind_meta") |> filter(gen == 1L) |> collect() |> head()
#> # A tibble: 6 × 8
#> id_ind id_parent_1 id_parent_2 line_name sex ploidy gen farm
#> <chr> <chr> <chr> <chr> <chr> <int> <int> <chr>
#> 1 A_501 A_1 A_251 A M 2 1 NA
#> 2 A_502 A_1 A_252 A F 2 1 NA
#> 3 A_503 A_1 A_253 A M 2 1 NA
#> 4 A_504 A_1 A_254 A F 2 1 NA
#> 5 A_505 A_2 A_255 A M 2 1 NA
#> 6 A_506 A_2 A_256 A F 2 1 NABecause you build the mating plan yourself, any design is possible — nested, factorial, reciprocal crosses, assortative mating — without the package needing a special argument for each.
pop
#> ── tidybreed population: demo ──────────────────────────────────────────────────
#> Database in-memory [connected]
#>
#> Genome 500 loci · 5 chr · 495 Mb · 490 cM
#> founder pool: 100 haplotypes
#> Model 2 traits · 2 phenotypes · 200 causal loci
#>
#> Individuals 520
#> by sex 260 F · 260 M
#>
#> Records phenotypes 500 · TGV 1,000
#>
#> schema(pop) · describe_table(pop, "name")
#> ────────────────────────────────────────────────────────────────────────────────SNP chips and genotypes
A chip is a named set of loci. Mark them with a filter:
pop <- pop |>
get_table("genome_meta") |>
filter(chr %in% c(1L, 2L, 3L)) |>
define_chip("50K")
#> Defined chip '50K' with 300 SNPs in column 'is_50K'
pop |> get_table("genome_meta") |> count(is_50K) |> collect()
#> # A tibble: 2 × 2
#> is_50K n
#> <lgl> <dbl>
#> 1 FALSE 200
#> 2 TRUE 300add_genotypes() records which animals were
genotyped — here, only the offspring:
pop <- pop |>
get_table("ind_meta") |>
filter(gen == 1L) |>
add_genotypes("50K")
#> Chip '50K': 20 animal(s) now marked as genotyped.
pop |> get_table("ind_meta") |> count(has_50K) |> collect()
#> # A tibble: 2 × 2
#> has_50K n
#> <lgl> <dbl>
#> 1 FALSE 500
#> 2 TRUE 20extract_genotypes() returns the familiar wide 0/1/2
dosage matrix, ready for GBLUP or export:
geno <- pop |> get_table("ind_meta") |> extract_genotypes("50K")
dim(geno)
#> [1] 20 301
geno[1:5, 1:6]
#> # A tibble: 5 × 6
#> id_ind locus_1 locus_2 locus_3 locus_4 locus_5
#> <chr> <int> <int> <int> <int> <int>
#> 1 A_501 0 0 1 1 1
#> 2 A_502 0 0 2 0 2
#> 3 A_503 0 1 2 0 2
#> 4 A_504 2 0 2 1 2
#> 5 A_505 1 2 2 0 2Haplotypes are the source of truth and dosages are derived on demand,
so nothing is stored twice. If you need dosages repeatedly,
add_dosage() caches them in ind_genotype.
Selection index
Register the weights:
pop <- pop |>
define_index("terminal",
trait_names = c("ADG", "BF"),
index_wts = c(1.2, -0.8))
#> Defined index 'terminal': ADG (1.2), BF (-0.8)
pop |> get_table("index_meta") |> collect()
#> # A tibble: 4 × 5
#> id_index_name index_name trait_name index_weight economic_weight
#> <int> <chr> <chr> <dbl> <dbl>
#> 1 1 NA ADG NA 0
#> 2 2 NA BF NA 0
#> 3 3 terminal ADG 1.2 NA
#> 4 4 terminal BF -0.8 NAThe rows with a missing index_name are the global
economic-weight entries that define_trait() writes for
every trait; the terminal rows are our index.
Because we know the true breeding values, we can compute the true index — useful for measuring how well selection actually worked:
pop <- pop |> get_table("ind_meta") |> add_tgv(index_names = "terminal")
#> Computed TGV for 520 individuals on trait 'ADG' (additive).
#> Computed TGV for 520 individuals on trait 'BF' (additive).
#> Computed true index 'terminal' (index, additive) for 520 individuals.
pop |> get_table("ind_true_index") |> collect() |> head()
#> # A tibble: 6 × 6
#> id_true_index id_ind index_name weight_type component_name true_index_value
#> <int> <chr> <chr> <chr> <chr> <dbl>
#> 1 1 A_1 terminal index additive -0.587
#> 2 2 A_10 terminal index additive -0.639
#> 3 3 A_100 terminal index additive 0.449
#> 4 4 A_101 terminal index additive -0.438
#> 5 5 A_102 terminal index additive -0.401
#> 6 6 A_103 terminal index additive 0.179add_index() does the same arithmetic on any table of
values. In a real program you would run it on ind_ebv after
a BLUP evaluation; here we use the true breeding values. Filter
first so there is exactly one value per individual per trait
(ind_tgv has a row per component):
pop <- pop |>
get_table("ind_tgv") |>
filter(component_name == "additive", trait_name %in% c("ADG", "BF")) |>
add_index("terminal")
#> Computing index 'terminal': ADG (wt=1.2), BF (wt=-0.8)
#> Added index 'terminal' (run #1) for 520 individuals
pop |>
get_table("ind_index") |>
collect() |>
arrange(desc(index_value)) |>
head(5)
#> # A tibble: 5 × 5
#> id_index id_ind index_name index_number index_value
#> <int> <chr> <chr> <int> <dbl>
#> 1 252 A_325 terminal 1 1.66
#> 2 443 A_498 terminal 1 1.57
#> 3 50 A_143 terminal 1 1.53
#> 4 300 A_369 terminal 1 1.38
#> 5 327 A_393 terminal 1 1.37Those top animals are your selection candidates. Pull their IDs,
build a new mating plan, call add_offspring(), and you have
closed the loop on a breeding cycle.
Close the population
close_pop(pop)For an on-disk run, close_pop() flushes everything to
the .duckdb file and you can pick the simulation back up
later:
pop <- restore_pop("tidybreed_output/sim.duckdb")Nothing needs to be re-specified — the genome, traits, variance components, and every animal are already in the file. This is the payoff of keeping configuration in the database rather than in R.
Where to go next
This vignette covered the main path. The package has more, grouped by prefix:
| Prefix | Meaning |
|---|---|
open_ / restore_
|
Create or reopen a population |
define_ |
Write model configuration — how the simulation should behave |
add_ |
Write simulation output — data the model produced |
mutate_ |
Add or update columns on an existing table |
get_ / extract_
|
Read data out, lazily or in analysis format |
remove_ / archive_
|
Delete rows, or stamp a finished replicate |
Functions not shown above, with a pointer to their help page:
| Function | What it does |
|---|---|
?add_ebv |
Run BLUPF90 (or parent average) and store estimated breeding values |
?define_chromosome |
Sex chromosomes, organelles, and achiasmatic meiosis |
?define_phenotype |
Composite traits: maternal effects and social genetic effects, via
components
|
?define_residual_cov |
Correlated or heterogeneous residual (co)variances |
?define_effect_random |
Named random effects such as pen, litter, or herd-year-season |
?define_effect_fixed_cov |
Fixed regression on a covariate |
?archive_replicate |
Collect many replicates into one archive database |
?remove_rows |
Delete rows safely across related tables |
?mutate_group_seq |
Litter and group utilities (mutate_group_*) |
?define_table |
Create your own tables inside the population |
The project README covers installation, global options, the output directory layout, and scenario YAML files in more detail.