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_tbv(), 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 architecture | define_trait() |
target_add_var, target_add_mean,
expressed_parent, units
|
| 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/RtmpGwVAPE/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
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 ─────────────────────────────── 22 tables · 1.9 MiB in memory ──
#> Use describe_table(pop, "name") for column-level details.
#>
#> Genome
#> genome_meta 500 5 Locus-level metadata. One row per l...
#> genome_map 500 7 Genetic map in long format. One row...
#> chr_inheritance 5 5 Per-chromosome copy counts, keyed b...
#> chr_recombination 5 4 Per-chromosome recombination, keyed...
#> + 1 empty: genome_effects
#>
#> 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
#> + 5 empty: ind_tbv, 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. target_add_var is the additive
genetic variance the QTL effects will be scaled to hit.
pop <- pop |>
define_trait("ADG", target_add_var = 0.25, units = "kg/day")
#> Added trait 'ADG'.
pop |> get_table("trait_meta") |> collect()
#> # A tibble: 1 × 6
#> id_trait trait_name description units expressed_parent target_add_mean
#> <int> <chr> <chr> <chr> <chr> <dbl>
#> 1 1 ADG NA kg/day both 0Note 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.
pop <- pop |>
get_table("genome_meta") |>
filter(chr %in% c(4L, 5L)) |>
define_additive_effects("ADG")
#> Set additive effects for 200 QTL on trait 'ADG' (base: founder_haplotypes).
pop |> get_table("genome_effects") |> collect() |> head()
#> # A tibble: 6 × 7
#> id_genome_effect locus_name line_name trait_name genome_effect_type
#> <int> <chr> <chr> <chr> <chr>
#> 1 1 Locus_301 NA ADG additive
#> 2 2 Locus_302 NA ADG additive
#> 3 3 Locus_303 NA ADG additive
#> 4 4 Locus_304 NA ADG additive
#> 5 5 Locus_305 NA ADG additive
#> 6 6 Locus_306 NA ADG additive
#> # ℹ 2 more variables: genome_value <dbl>, base_allele_freq <dbl>QTL membership is implicit: a locus is a QTL for a
trait if it has a row in genome_effects. There is no
is_QTL flag to keep in sync. base_allele_freq
is stored alongside each effect because it is what centres breeding
values.
Define the phenotype
The observation layer. With target_add_var = 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_tbv() computes true breeding values — the
simulation’s ground truth, which you would never know in a real
population.
pop <- pop |> get_table("ind_meta") |> add_tbv("ADG")
#> Computed TBV for 500 individuals on trait 'ADG'.
pop |> get_table("ind_tbv") |> collect() |> head()
#> # A tibble: 6 × 4
#> id_tbv id_ind trait_name tbv_value
#> <int> <chr> <chr> <dbl>
#> 1 5 A_5 ADG 0.778
#> 2 13 A_13 ADG -0.929
#> 3 16 A_16 ADG 1.02
#> 4 22 A_22 ADG -0.0430
#> 5 25 A_25 ADG -0.132
#> 6 31 A_31 ADG -0.334add_phenotype() builds the observed record: intercept,
plus fixed and random effects, plus the breeding value, plus a sampled
residual. It calls add_tbv() internally, so you can go
straight to it.
pop <- pop |> get_table("ind_meta") |> add_phenotype("ADG")
#> Computed TBV for 500 individuals on trait 'ADG'.
#> Wrote 500 phenotype records for 'ADG'.
pop |> get_table("ind_phenotype") |> collect() |> head()
#> # A tibble: 6 × 5
#> id_phenotype id_ind phenotype_name pheno_value pheno_number
#> <int> <chr> <chr> <dbl> <int>
#> 1 1 A_1 ADG 1.71 1
#> 2 2 A_2 ADG 1.58 1
#> 3 3 A_3 ADG -0.319 1
#> 4 4 A_4 ADG 0.848 1
#> 5 5 A_5 ADG 1.92 1
#> 6 6 A_6 ADG 2.94 1The 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 0.976
#> 2 M 250 1.38A second, genetically correlated trait
Real programs select on several traits at once. Define the second trait, then supply a genetic covariance matrix.
pop <- pop |> define_trait("BF", target_add_var = 0.30, 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 |> define_effect_cov_matrix("gen_add", G)
#> Stored 'gen_add' covariance matrix for: ADG, BF.
pop |> get_table("trait_var_comp") |> collect()
#> # A tibble: 4 × 5
#> id_trait_var_comp effect_name trait_name_1 trait_name_2 cov_value
#> <int> <chr> <chr> <chr> <dbl>
#> 1 1 gen_add ADG ADG 0.25
#> 2 2 gen_add ADG BF 0.08
#> 3 3 gen_add BF ADG 0.08
#> 4 4 gen_add BF BF 0.3Covariance matrices are stored in the database, not passed around in
R. Now one call draws correlated effects for both traits from
MVN(0, G):
pop <- pop |>
get_table("genome_meta") |>
filter(chr %in% c(4L, 5L)) |>
define_additive_effects(c("ADG", "BF"))
#> Set correlated additive effects for traits: ADG, BF (method: shared)
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()for a trait always overwrites its rows ingenome_effects, so recompute anything derived from them.
pop <- pop |> get_table("ind_meta") |> add_tbv() # no trait_name = all traits
#> Computed TBV for 500 individuals on trait 'ADG'.
#> Computed TBV for 500 individuals on trait 'BF'.
pop |>
get_table("ind_tbv") |>
collect() |>
group_by(trait_name) |>
summarise(n = n(), mean_tbv = mean(tbv_value), var_tbv = var(tbv_value))
#> # A tibble: 2 × 4
#> trait_name n mean_tbv var_tbv
#> <chr> <int> <dbl> <dbl>
#> 1 ADG 500 -0.00520 0.220
#> 2 BF 500 0.00673 0.312The 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. target_add_var scales the QTL
effects so the additive variance comes out right in the founder
haplotype pool. The animals you actually created are a finite
sample from that pool, so their realised variance wobbles around the
target. Larger founder populations wobble less.
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 = 1393264544)
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 QTL
#>
#> Individuals 520
#> by sex 260 F · 260 M
#>
#> Records phenotypes 500 · TBV 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 TRUE 300
#> 2 FALSE 200add_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 1 0 1 1 1
#> 2 A_502 1 0 2 1 2
#> 3 A_503 1 1 2 0 2
#> 4 A_504 1 0 2 0 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_tbv(index_names = "terminal")
#> Computed TBV for 520 individuals on trait 'ADG'.
#> Computed TBV for 520 individuals on trait 'BF'.
#> Computed true index 'terminal' (index) for 520 individuals.
pop |> get_table("ind_true_index") |> collect() |> head()
#> # A tibble: 6 × 5
#> id_true_index id_ind index_name weight_type true_index_value
#> <int> <chr> <chr> <chr> <dbl>
#> 1 1 A_1 terminal index -0.455
#> 2 2 A_10 terminal index 0.338
#> 3 3 A_100 terminal index 0.167
#> 4 4 A_101 terminal index 0.431
#> 5 5 A_102 terminal index 0.522
#> 6 6 A_103 terminal index 0.0660add_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 TBVs. Filter first
so there is exactly one value per individual per trait:
pop <- pop |>
get_table("ind_tbv") |>
filter(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 381 A_441 terminal 1 1.64
#> 2 509 A_89 terminal 1 1.44
#> 3 465 A_517 terminal 1 1.43
#> 4 486 A_68 terminal 1 1.33
#> 5 352 A_415 terminal 1 1.28Those 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 |
?define_trait_simple |
Shortcut wrapper: define_trait() +
define_additive_effects()
|
?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.