Tidy Breeding Program Simulation • tidybreed

Breeding program simulation, backed by DuckDB, driven by dplyr.

📖 Documentation: https://austin-putz.github.io/tidybreed/

[!IMPORTANT] tidybreed is very much in the “alpha testing” phase, I do not recommend building on it just yet. I am however looking for valuable feedback to finish the API prior to version 1.0.0

A pipe-friendly (%>% or |>) R package for breeding program simulation backed by DuckDB. Design large-scale genomic simulations without running out of memory — all data lives on disk in a DuckDB database and is queried lazily via dplyr.

It is easier to show than explain. The API often resembles the following steps:

pop |>                       # stores 'db_conn' connection to your .duckdb file
  get_table("ind_meta") |>   # stores 1 record per individual to track pedigree, sex, line, any user defined fields
  filter(
    sex == "M",                  # subset to males only
    line_name == "Angus",        # subset to Angus genetic line only
    birth_date == Sys.Date(),    # user added column to track birth dates
    status == "calf"           # user defined column to track animal 'status'
  ) |>
  add_phenotype(
    c("birth_weight", "stillborn"),       # list of phenotypes to 'collect' on filtered individuals
    phenotype_date = Sys.Date()           # user added column, added to ind_phenotype table
  )

Steps:

1️⃣ User identifies a table to pass via get_table()
2️⃣ User filters/subsets the table to identify rows/records/animals
- filter() is from dplyr, users often already know this package well
3️⃣ User passes that table to either a define_*() function or add_*() function
- define_*() ➡️ functions often add records as meta data
- add_*() ➡️ functions often add data on individuals (e.g. TBV, EBV, Phenotype, Genotype, etc)
4️⃣ New records are stored on disk to your .duckdb file (database)
- database allows for efficient storage and IO
- users allowed to insert rows themselves without me having to provide ‘helper functions’

[!NOTE] tidybreed uses a lightweight S3 object as a handle to a DuckDB-backed breeding program. Rather than storing simulation state in R objects, the breeding program is represented as a persistent relational database. Functions operate by querying subsets of the database, performing calculations in R or C++, and writing results back to the database. The S3 object primarily manages access to the underlying database rather than storing biological state.

Motivation

I’ve both struggled to learn other simulation software and also found them extremely limiting for several reasons. tidybreed started as a fresh conceptual idea on simulation software to allow users to easy store custom data (tables + columns) and easily manipulate any data they want (metadata or otherwise).

I also try to make it as intuitive and easy to learn as possible. Users can utilize many functions they already know from tidyverse such as filter() or slice_*() functions (for selection for instance).

License

MIT License

tidybreed is released under the MIT License. You are free to use, modify, and distribute it, including in commercial and proprietary projects, provided the copyright notice is retained.

[!CAUTION] NO WARRANTY — This software is provided as-is, without warranty of any kind, express or implied. The authors accept no liability for any damages or losses arising from its use.

Design

Main		Reason or Description
R	→	Standard for most scientists; flexible design of custom breeding programs
DuckDB	→	Columnar, embedded, no server needed; handles datasets larger than RAM
Pipe everything	→	Filter individuals, add phenotypes/genotypes/EBVs/index with `tidyverse` verbs
Customizable	→	Add your own table, any column within any table with `mutate_table()`; query with standard `dplyr` and SQL will function in the background. `DBI` allows users to interact with the database directly in any way they want.

Databases ✅ are efficient enough for our large simulations, both in speed and memory savings. Users can pull from different tables using the get_table() function. We no longer need to fill up our RAM with data from generation 1 in a long term simulation pipeline.
SQL ✅ can add any custom table or field/column of any type, allowing users to define a DATE that would allow almost near perfect “digital twins”. My understanding is DuckDB was written very efficiently with C++ mostly. Users can DBI their way to modifying any table at any time and often do not have to worry about “breaking” tidybreed if I build the functions correctly and truly modular.
R ✅ is great for many things and has become a standard for us in research, however memory is a massive issue compared to python/julia/etc and it’s simply not efficient enough for massive datasets of this size. This will be moderated by the use of SQL and C++ with Rcpp and many operations will (hopefully) never be pulled into R, but can be orchestrated through R.
Pipes ✅ allow users to insert filter steps that is critical for most operations to “point” to certain individuals or groups to calculate TBVs, EBVs, phenotypes, or to extract data such as genotypes/QTL.

What I try to AVOID:

❌ Storing metadata
- such a n_loci or such, a truly composable/modular system using SQL can easily calculate this from a table as needed with very little “cost” to timing of the simulation
- metadata is stored, just in tables such as the variance components but other fixed metadata is silly most of the time when it can be taken from a table at any moment in time cheaply
❌ helper/wrapper functions
- any function that exists, but the user could easily do it themselves
- e.g. a function to design a specific mating scheme
- users love them, developers love them, but they are almost always a bad idea leading to spaghetti code long term.
- Why? 👉 There is a nearly infinite number of ways users could mate a list of dams to a list of sires, I cannot program every possible option for users, therefore I designed a simple tibble as input and users just list exactly what offspring they want per row (sire, dam, sex, etc). No helper function needed to create the mating design, you can do this yourself easily without a function using sample(), rep(), etc. and building a tibble yourself. This also avoids me having to break your code later changing functions and arguments.
❌ generation 💀 🆘 🚫
- generation was a leftover archaic artifact from QMSim and the original AlphaSim I think
- generation has no useful application besides making some simulations easier to follow or design in a generic/fixed way
- generation only exists in simple simulations and selection experiments, not real breeding programs (for the most part)
- users can add a field/column called “gen” or “generation” within tidybreed in every table they want, however I do not enforce you to use such a silly idea unless doing some simple simulations and need to create something quick/easy

Installation

Install pak then install tidybreed from GitHub:

WINDOWS People ➡️ YOU NEED RTools SETUP, please read below

install.packages("pak", repos = "https://packagemanager.posit.co/cran/latest")
library(pak)
pak::pak("austin-putz/tidybreed")
library(tidybreed)

System requirements: a C++ compiler

As of v0.53.0, tidybreed ships compiled C++17 code (the recombination kernel used by add_offspring()). Installing from GitHub is a source install, so each machine needs a working C++ toolchain. tidybreed’s header-only build dependencies (Rcpp, dqrng, BH, sitmo) are pulled in automatically by pak — no system libraries are required — but the compiler itself must be present. If it is missing, the install stops with a platform-specific message telling you exactly what to install.

Platform	What to install	Command / link
Windows 11 / 10	Rtools matching your R version (Rtools45 for R 4.5.x, Rtools44 for R 4.4.x)	Download & run the installer from https://cran.r-project.org/bin/windows/Rtools/, then restart R. Rtools puts `g++`/`make` on `PATH` for you.
Ubuntu / Debian	`r-base-dev` + `build-essential`	`sudo apt-get update && sudo apt-get install -y r-base-dev build-essential`
Fedora / RHEL / Rocky	`gcc-c++` + `make`	`sudo dnf install -y gcc-c++ make`
macOS	Xcode Command Line Tools (provides `clang++`)	`xcode-select --install`

[!NOTE] macOS: the Command Line Tools alone are enough — you do not need the full Xcode app. If clang++ was already installed with a previous R/dev setup, no action is needed. (Apple Silicon and Intel are both supported; the kernel is portable C++17 with no architecture-specific code.)

Once the toolchain is in place, re-run the pak::pak(...) command above. To force a fresh recompile after adding a compiler:

pak::pak("austin-putz/tidybreed", upgrade = TRUE)

To force the pure-R reference kernel (e.g. to cross-check results, or on a machine with no compiler where you install a binary build), set an environment variable before loading the package:

Sys.setenv(TIDYBREED_KERNEL = "r")   # "auto" (default) uses the compiled C++ kernel

[!WARNING] Pre-v1.0.0 packages are considered beta and subject to breaking changes. Pin your version to avoid surprises.

packageVersion("tidybreed")   # check installed version, e.g. '0.46.1'
pak::pak("austin-putz/tidybreed@v0.46.1")   # pin to a specific release

Browse all releases on the GitHub Releases page.

Global Options

Set package-wide options at the top of your script (or in a startup / config script) so simulations can be parameterized without changing function calls throughout the codebase. Every option has a built-in default and is entirely optional — set only the ones you need. Options prefixed tidybreed. are read by open_pop() (folder layout) and archive_replicate() (multi-replicate archiving).

options(
  tidybreed.pop_name        = "my_project",
  tidybreed.base_dir        = "~/path/to/project/",
  tidybreed.output          = "tidybreed_output",
  tidybreed.scenario        = "baseline_scenario",
  tidybreed.tools           = c("blupf90", "plink", "JWAS"),
  tidybreed.db_name         = "sim.duckdb",
  tidybreed.replicate       = 1L,
  tidybreed.archive_path    = "~/path/to/project/archive/",
  tidybreed.db_name_archive = "baseline_scenario_all_reps.duckdb",
  tidybreed.quiet           = FALSE
)

tidybreed.pop_name — Label stored on the pop object (pop$pop_name); shown in print(pop). Default: "sim". Purely descriptive — does not affect file paths.
tidybreed.base_dir — Root folder (layer 1); default: getwd()
tidybreed.output — Output subfolder name (layer 2); default: "tidybreed_output"
tidybreed.scenario — Scenario subfolder (layer 3). Default NULL auto-generates a YYYYMMDD_HHMMSS folder so runs never overwrite each other. Set explicitly (e.g. "baseline") to reuse the same folder, such as in an HPC array job.
tidybreed.tools — Character vector of tool subfolders created at layer 4 (e.g. c("blupf90", "plink")). Default NULL skips tool folder creation.
tidybreed.db_name — Working DuckDB file name; default: "sim.duckdb". Use ":memory:" for an in-memory database (skips all folder creation — useful for tests).
tidybreed.replicate — Integer replicate number stamped on archive rows; default: 1L. Auto-increments after each successful archive_replicate() call.
tidybreed.archive_path — Directory for the archive DuckDB file. Default NULL places the archive next to the working database.
tidybreed.db_name_archive — Archive DuckDB file name (e.g. "all_reps.duckdb"). Default NULL skips archiving entirely.
tidybreed.quiet — Suppress the startup banner on library(tidybreed); default: FALSE.

[!NOTE] Archive path resolution (archive_replicate()) — first non-NULL wins:

Explicit archive_path argument passed to archive_replicate().

file.path(tidybreed.archive_path, tidybreed.db_name_archive) — when both options are set.

file.path(dirname(pop$db_path), tidybreed.db_name_archive) — archive lands next to the working database.

tidybreed.db_name_archive is NULL — no archive written; only reset phases run.

Directory Layout

The recommended project structure separates scenario configs, per-scenario simulation databases, and the merged archive:

~/projects/swine/                        ← tidybreed.base_dir
│
├── scenarios/                           ← one YAML per scenario
│   ├── baseline.yaml
│   └── scenario_b.yaml                  ← next scenario to run
│
├── results/                                ← merged multi-replicate archive
│   └── baseline_scenario_all_reps.duckdb   ← written by archive_replicate()
│
└── tidybreed_output/                    ← tidybreed.output
    ├── baseline/                        ← tidybreed.scenario = "baseline"
    │   ├── sim.duckdb                   ← tidybreed.db_name
    │   ├── blupf90/                     ← subfolder per tool in tidybreed.tools
    |   |-- JWAS/
    │   └── plink/
    │
    └── scenario_b/                      ← tidybreed.scenario = "scenario_b"
        ├── sim.duckdb
        ├── blupf90/
    |   |-- JWAS/
        └── plink/

Scenario YAML format

Save scenario parameters in scenarios/baseline.yaml and read them at the top of your script with yaml::read_yaml():

# scenarios/baseline.yaml
pop_name:  swine
base_dir:  ~/projects/swine/
output:    tidybreed_output
scenario:  baseline
db_name:   sim.duckdb
tools:
  - blupf90
  - plink

cfg <- yaml::read_yaml("scenarios/baseline.yaml")

options(
  tidybreed.pop_name  = cfg$pop_name,
  tidybreed.base_dir  = cfg$base_dir,
  tidybreed.output    = cfg$output,
  tidybreed.scenario  = cfg$scenario,
  tidybreed.db_name   = cfg$db_name,
  tidybreed.tools     = cfg$tools,
  tidybreed.replicate = as.integer(Sys.getenv("REP", 1))  # set per run
)

Swap scenarios by pointing to a different YAML file; the rest of the script stays identical.

[!TIP] I HIGHLY encourage the use of purrr::chuck() to pull from the stored yaml structure (cfg above). Unlike $, it throws an error when the element is missing instead of returning NULL and silently breaking your pipeline without a clear warning.

Core Concept

Every function accepts and returns a tidybreed_pop object — a thin wrapper around a DuckDB connection. Chain operations with |> (or %>%):

pop <- pop |>
  get_table("ind_meta") |>       # identify which DB table to work with
  filter(                        # filter to the individuals you want
    sex == "M",                  # filter to males
    gen == 1L                    # filter to generation 1 (this was "user defined") 
  ) |>  
  add_phenotype("ADG")           # add records to ind_phenotype

All tables are queryable at any time. Use collect() to pull results into R as a tibble.

Simulation Workflow

1. Open a population

open_pop() creates or re-opens a DuckDB-backed population object. Use clean = TRUE to start fresh.

pop <- open_pop(
  clean = TRUE   # wipe and recreate the database
)

print(pop)

2. Define the genome

pop <- pop |>
  define_genome(
    n_loci      = 10000,  # total number of loci (SNP + QTL)
    n_chr       = 18,     # number of chromosomes
    chr_len_Mb  = 50,     # physical length per chromosome, Mb (scalar or length n_chr)
    cM_per_Mb   = 1.0     # genetic-map rate (scalar or length n_chr)
  )

pop |> get_table("genome_meta")   # 1 row per locus — physical map (pos_bp)
pop |> get_table("genome_map")    # 1 row per locus — genetic map (pos_cM)

define_genome() is the single call that builds the genome. It writes seven tables in one transaction — genome_meta, genome_map, ind_haplotype, ind_genotype, ind_crossover, chr_inheritance, and chr_recombination — and rolls all of them back if anything fails, so a bad call leaves the database exactly as it found it and the same pop can be reused for a corrected call.

Physical and genetic coordinates live in separate tables. genome_meta.pos_bp is the single physical source of truth (BIGINT, 1-based, VCF/PLINK convention); pos_cM = pos_bp / 1e6 * cM_per_Mb is written to genome_map as the default map (sex = NULL, line_name = NULL, map_name = "default"). Sex- and line-specific maps are added later as extra genome_map rows — rows, never a schema change. Everything distance-driven (founder LD, recombination) reads the resolved map, so cM_per_Mb sets the crossover rate: a 50 Mb chromosome at cM_per_Mb = 1.0 is 50 cM, i.e. ~0.5 crossovers per meiosis.

Optional arguments:

Argument	Purpose
`locus_names`	Custom locus names (length `n_loci`). Must be unique, non-`NA`, non-empty. Default `Locus_1 … Locus_n`.
`chr_names`	Custom chromosome names (length `n_chr`), same rules. Default `"1" … "n"`.
`recombines_M` / `recombines_F`	Genome-wide per-parent-sex recombination defaults (both `TRUE`). Set one to `FALSE` for a whole-genome achiasmatic sex (e.g. Drosophila males). Seeded into `chr_recombination`.

[!IMPORTANT] define_genome() may be called once per population. It errors if any of the seven genome tables already exists — there is no partial re-definition. n_loci and n_chr must be whole numbers with n_loci >= n_chr, and chr_len_Mb / cM_per_Mb must be finite and strictly positive. All of these are checked before anything is written.

Per-chromosome exceptions (sex chromosomes, organelles) are set afterwards with define_chromosome(), which writes chr_inheritance rows (keyed by offspring sex) or chr_recombination rows (keyed by producing-parent sex) — one concern per call.

3. Define founder haplotypes

Generate haplotype pools for each genetic line. Six methods are available: four that draw a per-locus allele frequency and sample alleles independently (no LD), and two that build LD along the genetic map.

# Uniform allele frequencies between min and max (no LD)
pop <- pop |>
  define_founder_haplotypes(
    line_name       = "A",
    n_haplotypes    = 1000,
    method          = "uniform",
    min_allele_freq = 0.01,
    max_allele_freq = 0.99
  )

# Fixed frequency — realized exactly at every locus (no LD)
pop <- pop |>
  define_founder_haplotypes(line_name = "B", n_haplotypes = 1000,
    method = "fixed", allele_freq = 0.5)

# Beta(0.5, 0.5) — U-shaped MAF spectrum, biologically realistic (no LD)
pop <- pop |>
  define_founder_haplotypes(line_name = "C", n_haplotypes = 1000,
    method = "beta", beta_shape1 = 0.5, beta_shape2 = 0.5)

# Balding-Nichols (FST-based drift around an ancestral mean; no LD)
pop <- pop |>
  define_founder_haplotypes(line_name = "D", n_haplotypes = 1000,
    method = "balding_nichols", fst = 0.1, mean_allele_freq = 0.5)

# Mosaic — LD via Li-Stephens haplotype-block copying
pop <- pop |>
  define_founder_haplotypes(line_name = "E", n_haplotypes = 1000,
    method = "mosaic", n_templates = 32, template_switch_rate = 1.0)

# Gaussian copula — LD via AR(1) latent normal (fast, unquantized MAF)
pop <- pop |>
  define_founder_haplotypes(line_name = "F", n_haplotypes = 1000,
    method = "gaussian_copula", ld_decay_rate = 0.25)

Each method owns its own arguments, and passing one to the wrong method is a hard error naming the method it belongs to:

Method	Arguments	LD
`"uniform"` (default)	`min_allele_freq` (0.01), `max_allele_freq` (0.99)	none
`"fixed"`	`allele_freq` (0.5)	none
`"beta"`	`beta_shape1` (0.5), `beta_shape2` (0.5)	none
`"balding_nichols"`	`fst` (0.1), `mean_allele_freq` (0.5)	none
`"mosaic"`	`n_templates`, `template_switch_rate` (1.0 per cM)	blocks
`"gaussian_copula"`	`ld_decay_rate` (1.0; ρ = exp(−λ·d_cM))	AR(1) decay

All four frequency-based methods also accept exact_freq. When TRUE, each locus gets exactly round(p × n_haplotypes) copies of the 1-allele on an independently drawn random subset of haplotypes, so the realized pool frequency equals the target with no binomial scatter (and still no LD — each locus draws its own subset). When FALSE, alleles are independent Bernoulli(p) draws and realized frequencies scatter around p with sd sqrt(p(1-p)/n_haplotypes).

exact_freq defaults to TRUE for method = "fixed" (a “fixed” frequency that drifts is not fixed) and FALSE for the distribution-based methods, where drawing a frequency and then sampling binomially is the correct generative model. Set it TRUE on those methods when you want a drift-free base frequency. Frequencies live on a 1 / n_haplotypes grid; for "fixed" an off-grid request warns and names the frequency actually used.

# Beta frequencies, realized exactly (no binomial scatter around the draw)
pop <- pop |>
  define_founder_haplotypes(line_name = "G", n_haplotypes = 1000,
    method = "beta", exact_freq = TRUE)

Both LD methods read the genetic map (genome_map) resolved for that pool’s own line_name, so founder LD is built on the same map that later drives recombination for that line.

[!TIP] In "mosaic", n_templates controls the MAF spectrum, not just block length. Templates are the only source of allelic variation, so a locus where all templates agree is monomorphic regardless of n_haplotypes — roughly 2 / (n_templates + 1) of loci, with MAF quantized to multiples of 1 / n_templates. A warning fires when more than 10% of loci come out monomorphic, because QTL placed there contribute nothing to sum(2pq a²). Raise n_templates, or use "gaussian_copula" for a dense, unquantized MAF spectrum. Note also that the template re-draw is uniform over all templates including the current one (the standard Li-Stephens kernel — it makes realized LD invariant to marker density), so observable template changes occur at template_switch_rate × (n_templates − 1) / n_templates.

Calling with line_name = NULL (the default) stores one shared pool that add_founders() falls back to when no named pool exists for the line it is building. Calling again with a line_name that already has a pool — or again with NULL once a NULL-line pool exists — is an error.

Every call also writes a founder_allele_freq column to genome_meta holding the per-locus frequency of the pool written most recently.

[!NOTE] founder_allele_freq is informational only — no other tidybreed function reads it, and each call rewrites it for every locus, so in a multi-line setup it describes only the last pool written. For per-line Falconer centering use base = "current_pop" with a line-filtered base_tbl in define_additive_effects(); base = "founder_haplotypes" recomputes the base frequency by pooling all lines together (which overstates within-line heterozygosity — Wahlund — and under-scales target_add_var).

4. Add founder individuals

Sample haplotypes for founders. Pass any custom ind_meta column as ... arguments — they are written atomically with the new rows.

pop <- pop |>
  get_table("founder_haplotypes") |>  # pass this table only
    filter(line_name == "A") |>       # neat way to filter and pass only haplotypes related to founder line "A"
  add_founders(
    n_males    = 400,
    n_females  = 1600,
    line_name  = "A",                   # required
    birth_date = sampled_birth_dates,   # user sampled birth dates earlier
    alive      = TRUE,                  # user defined - show as alive to filter by later
    active     = FALSE                  # user defined - or set 'status' with character/VARCHAR
  )

5. Add custom columns with `mutate_table()`

The real power of tidybreed is freely adding columns to any table so your simulation state is stored in the database — no parallel R objects to maintain.

# Add user-defined fields to ind_meta (declare schema before data arrives)
pop |>
  get_table("ind_meta") |>     # stores 1 row / individual to store pedigree (sire/dam) and sex
  mutate_table(
    status        = NA_character_,   # VARCHAR: production status (update with `mutate_table()`)
    birth_date    = as.Date(NA),     # DATE, initialize column with missing
    puberty_date  = as.Date(NA),
    mate_date     = as.Date(NA),
    farrow_date   = as.Date(NA),
    wean_date     = as.Date(NA),
    cull_date     = as.Date(NA),
    off_test_date = as.Date(NA),
    alive         = TRUE,            # BOOLEAN with default value
    active        = FALSE,
    .set_default  = TRUE             # if TRUE -> write this value as the column default when creating new rows
  )

Types are inferred from R values:

R value	DuckDB type
`0L`, `NA_integer_`	`INTEGER`
`0.0`, `NA_real_`	`DOUBLE`
`TRUE`/`FALSE`	`BOOLEAN`
`"text"`, `NA_character_`	`VARCHAR`
`as.Date(...)`	`DATE`
`Sys.time()`	`TIMESTAMP`

Users can then update the ‘status’ of an animal to closely mimic a real program

pop |> 
  get_table("ind_meta") |>     # stores 1 row / individual to store pedigree (sire/dam) and sex
    filter(
      sex == "M",
      off_test_date == current_date
    ) |> 
  mutate_table(
    status = "after-test-boar"
  )

Which now replaces the old ‘status’ of “on-test” perhaps. The next filter you set can filter out “after-test-boar” animals only to perform selection in any way you want.

Add descriptions to columns for documentation:

pop |>
  get_table("ind_meta") |>
  define_schema_description("status",     "Production status (e.g. gestation, lactation)") |>
  define_schema_description("birth_date", "Date of birth") |>
  define_schema_description("alive",      "Is the animal alive?")

[!TIP] If you get lost with tables (and you will…) please use schema(pop), summary(pop), and pop |> describe_table("table_name") to get your feet back.

schema(pop)                                # grouped by pipeline stage, empty tables collapsed
schema(pop, show_empty = TRUE)             # one row per table
schema(pop, order = "rows")                # flat, biggest first
schema(pop, order = "size", sizes = TRUE)  # on-disk bytes (issues a CHECKPOINT)

schema() prints tables in the order a population is built — Genome, Founders, Individuals, Genetic model, Observation model, Selection, Results — with the database size in the header. The grouping is also returned as a table_group column, so it can be filtered like any other data.

Pick out one table to focus on:

pop |> describe_table("ind_meta")    # view all column descriptions

6. Define a SNP chip

Filter genome_meta to the loci you want, then call define_chip(). Adds an is_<chip_name> boolean column to genome_meta.

# Random 9,000-locus chip
pop |>
  get_table("genome_meta") |>
    slice_sample(n = 9000) |>         # simple dplyr function you already know
  define_chip(chip_name = "9k")       # just name the chip, now `is_9k` field in `genome_meta`

# QTL are defined as loci NOT on the chip
pop |>
  get_table("genome_meta") |>
    filter(is_9k != TRUE) |>          # filter out SNP chip loci to become QTL
  define_additive_effects(...)

This power architecture allows users to basically pass whatever they want and since you can define your own field to genome_meta and filter yourself, the possibilities are quite advanced and you know exactly what is going on.

7. Define variance components

Use a single entry point for all variance/covariance matrices. effect_name = "gen_add" routes to trait_var_comp; "residual" and named random effects (e.g. "pen") route to phenotype_var_comp.

# Additive genetic (co)variance — 3 traits (ADG, WWD, WWM)
vars.mat.add <- matrix(c(
  0.0045, 0.00, 0.00,
  0.00,   0.04, 0.00,
  0.00,   0.00, 0.13),
  nrow = 3, byrow = TRUE,
  dimnames = list(c("ADG", "WWD", "WWM"),
                  c("ADG", "WWD", "WWM")))

# store matrix in a table called 'trait_var_comp'
pop <- pop |>
  define_effect_cov_matrix(effect_name = "gen_add", cov_matrix = vars.mat.add)

# Residual (co)variance — 2 phenotypes (ADG, WW)
vars.mat.res <- matrix(c(
  0.0067, 0.00,
  0.00,   0.45),
  nrow = 2, byrow = TRUE,
  dimnames = list(c("ADG", "WW"),
                  c("ADG", "WW")))

# store matrix in 'phenotype_var_comp' table
pop <- pop |>
  define_effect_cov_matrix(effect_name = "residual", cov_matrix = vars.mat.res)

# Named random effect — pen effect on ADG only (1x1)
vars.mat.pen <- matrix(0.0005, nrow = 1, byrow = TRUE,
  dimnames = list("ADG", "ADG"))

# add pen covariances to the 'phenotype_var_comp' table
pop <- pop |>
  define_effect_cov_matrix(effect_name = "pen", cov_matrix = vars.mat.pen)

8. Define traits and phenotypes (two-layer design)

Genetic layer (define_trait) — one row per underlying genetic “trait”:

pop <- pop |>
  define_trait(
    trait_name      = "ADG",
    description     = "Average Daily Gain",
    units           = "kg/d",
    target_add_mean = 0,      # TBV mean in base population
    overwrite       = TRUE
  )

Observation layer (define_phenotype) — what animals actually receive records for:

pop <- pop |>
  define_phenotype(
    phenotype_name = "ADG",
    type           = "continuous",   # "continuous", "count", "categorical", "derived_formula"
    mean           = 0.92,           # phenotypic population mean
    expressed_sex  = "both",         # "both", "M", or "F"
    min_value      = 0,
    overwrite      = TRUE
  )

Why the split? - composite traits exist that can be made up of multiple underlying “traits”, while the phenotype is physically observed. Some phenotypes are actually a combination of 3 “traits” such as the direct/maternal/nurse dam model from Egbert Knols research. Or social genetic effects where the phenotype is a combination of individual genetics but also the effects from different pen mates genetics.

For composite phenotypes (e.g. weaning weight = direct + maternal):

pop <- pop |>
  define_trait(
    trait_name      = "WWD",
    description     = "Weaning Weight - Direct",
    units           = "kg",
    target_add_mean = 0,      # TBV mean in base population
    overwrite       = TRUE
  ) |>
  define_trait(
    trait_name      = "WWM",
    description     = "Weaning Weight - Maternal",
    units           = "kg",
    target_add_mean = 0,      # TBV mean in base population
    overwrite       = TRUE
  )

pop <- pop |>
  define_phenotype(
    phenotype_name           = "WW",
    type                     = "continuous",         # simple continuous phenotype, 
    formula_tbv              = "WWD + dam(WWM)",     # DSL: self + dam TBV for 'MILK'
    mean                     = 6.0,
    expressed_sex            = "both",
    missing_component_action = "skip",   # skip founders with no dam record
    overwrite                = TRUE
  )

Derived phenotypes computed from other phenotypes:

pop <- pop |>
  define_phenotype(
    phenotype_name = "FCR",
    type           = "derived_formula",  # specify to make sure it builds this phenotype relative to other phenotypes observed
    formula        = "ADFI / ADG",       # DSL allows for more than just division, but this is most common
    expressed_sex  = "both",
    overwrite      = TRUE
  )

Add QTL effects for a trait (or multiple correlated traits at once):

# Single trait
pop |>
  get_table("genome_meta") |>
    filter(is_9k != TRUE) |>           # QTL are non-SNP-chip loci
  define_additive_effects(
    trait_name      = "ADG",
    distribution    = "normal",      # mostly set to normal for now until I can figure out other ways to sample multivariate
    scale_to_target = TRUE,          # scale to target additive variance
    base            = "current_pop"  # standardize to current animals
  )

# Multiple correlated traits in one call (draws from MVN with G matrix)
pop |>
  get_table("genome_meta") |>
    filter(is_9k != TRUE) |>
  define_additive_effects(
    trait_name      = c("WWD", "WWM"),   # need to add together since they are genetically correlated
    distribution    = "normal",
    scale_to_target = TRUE,
    base            = "current_pop"
  )

Add fixed effects:

pop |>
  define_effect_fixed_class(
    phenotype_name = "ADG",
    effect_name    = "sex",
    source_column  = "sex",
    levels         = c(M = 0.08, F = 0),     # male advantage in ADG
    source_table   = "ind_meta",
    overwrite      = TRUE
  )

9. Calculate True Breeding Values

# All animals in ind_meta
pop <- pop |>
    get_table("ind_meta") |>      # pass table with no filter, so "select all" 
  add_tbv(trait_name = "ADG")     # calculate TBV for all individuals (notice we use "traitn_name" here)

Then the real great part of the package is to simply extract any table, and utilize the functions you already know to verify it worked (or didn’t…).

# Check means by trait
pop |>
  get_table("ind_tbv") |>
  collect() |>                              # pull into R memory
  group_by(trait_name) |>
  summarise(MeanTBV = mean(tbv_value))

Define a selection index first, then compute true index values from TBVs (ground truth for monitoring genetic trend):

# Define the index weights (must exist before add_tbv uses index_names)
pop |>
  define_index(
    index_name   = "maternal",
    trait_names  = c("ADG", "WWD", "WWM"),
    index_wts    = c(ADG = 0.5, WWD = 0.3, WWM = 0.2),
    economic_wts = c(ADG = 0.5, WWD = 0.3, WWM = 0.2)
  )

# Compute true index values from TBVs and write to ind_true_index
pop |>
  get_table("ind_meta") |>
  add_tbv(index_names = "maternal")       # will save the true breeding value index in 'ind_true_index' table)

pop |> 
  get_table("ind_true_index") |> 
  collect()

10. Add Genotypes

# Add 9k genotypes to ONLY CURRENT MALES
pop |>
  get_table("ind_meta") |>
    filter(
      sex == "M"
    ) |>
  add_genotypes(chip_name = "9k")

# Extract genotypes for downstream analysis
pop |>
  get_table("ind_meta") |>
    filter(
      sex == "M"          # extract only for males
    ) |> 
  extract_genotypes(chip_name = "9k")

Of course, you can all use pull(id_ind) to extract a character vector and use that to filter rows as well to get a more precise set of animals for genotyping, phenotyping, or whatever you want.

11. Add Phenotypes

[!IMPORTANT] tidybreed is explicit in the separation of “trait” vs “phenotype”.

trait ➡️ is linked to the genome via QTL effects

phenotype ➡️ is linked to the observed phenotype

This allows us to separate something like weaning weight into 2️⃣ components weaning weight direct and weaning weight maternal.

# Phenotype male animals that reached off-test today
pop |>
  get_table("ind_meta") |>
    filter(
      sex == "M",
      off_test_date == cur_date         # power to loop over dates allows real breeding program dynamics
    ) |>
  add_phenotype(
    phenotype_name = c("ADG", "BF", "IMF"),   # add ADG, Backfat, and IMF for these filtered animals
    pheno_date     = cur_date          # this is a user defined column in the 'ind_phenotype' table, pass it here
  )

# Phenotype all females for age at puberty (sex-limited trait)
pop |>
  get_table("ind_meta") |>
    filter(sex == "F") |>              # only females get age at puberty
  add_phenotype(phenotype_name = "AP")

Check counts:

pop |> 
  get_table("ind_phenotype") |> 
  filter(phenotype_name == "ADG") |> 
  collect() |> 
  count()

12. Run Evaluations (EBVs)

[!WARNING] This is probably the sketchiest part of any simulation software as it’s nearly impossible to accurately solve all models you can simulate. Simulation is far easier than solving equations. If users define a very complex structure to simulation, we cannot guarantee there is a solver out there for it. BLUPF90 is nice and easy, however extremely limited and we get what we get from it. Safest always to run user defined parameter files and programs. PLEASE BE CAREFUL!

Run BLUPF90 to estimate breeding values:

pop <- pop |>
  get_table("ind_meta") |>
  add_ebv(
    "ADG",
    software  = "blupf90",
    model     = "blup",
    phenotype = pop |>                      # this was critical to implement as some records will need to be sampled "early"
      get_table("ind_phenotype") |>
      filter(pheno_date < cur_date | is.na(pheno_date)),  # exclude future records for evaluations
    eval_date = cur_date
  )

pop |> get_table("ind_ebv") |> filter(trait_name == "ADG")

13. Define and Calculate a Selection Index

# Define a multi-trait selection index
pop |>
  define_index(
    index_name   = "maternal",
    trait_names  = c("ADG", "WWD", "WWM"),
    index_wts    = c(1,2,3),    # simple vector of index weights
    economic_wts = c(4,5,6)     # we can use later to establish profitability
  )

# Calculate index values from latest EBVs (run after add_ebv)
pop |>
  get_table("ind_ebv") |>
    filter(
      eval_date == cur_date         # only extract the latest EBVs to calc INDEX
    ) |>
  add_index("maternal", index_date = cur_date)

pop |> get_table("ind_index")

Of course in crossbreeding, each line may have their own index, just name it and pass only those line animals to the add_index() function and boom, you have your custom indexes per line.

14. Select Parents

Pull candidate IDs and then select from any table:

pull() may be new to you, it pulls a vector back into the R object to use within filter() downstream

# Step 1: identify candidates
male_candidates <- pop |>
  get_table("ind_meta") |>
    filter(
      sex == "M", 
      status %in% c("after-test-boar", "breeding-boar")
    ) |>
  pull(id_ind)

# Step 2: select top animals by index
selected_males <- pop |>
  get_table("ind_index") |>
    filter(
      id_ind %in% male_candidates, 
      index_date == latest_index_date
    ) |>
  slice_max(index_value, n = 10) |>
  pull(id_ind)

# Update status
pop |>
  get_table("ind_meta") |>
  filter(id_ind %in% selected_males) |>
  mutate_table(status = "breeding-boar")

15. Add Offspring

Build a mating plan as a tibble (one row per offspring), then call add_offspring():

(again I hate helper functions as I cannot imagine the infinite ways all these combinations could arise, for now just use the full tibble and learn how to use rep and samle yourself)

data.matings <- tibble(
  id_parent_1   = rep(sampled_sires, times = nw_per_dam),  # sire IDs
  id_parent_2   = rep(selected_dams, times = nw_per_dam),  # dam IDs
  line_name     = "A",
  sex           = sample(c("M", "F"), size = n_offspring, replace = TRUE),
  conc_date     = cur_date,
  birth_date    = cur_date + 116L,
  on_test_date  = cur_date + 116L + 70L,
  off_test_date = cur_date + 116L + 160L
)

pop |> add_offspring(data.matings)

116L refers to the average gestation length of sows these days, and 70/160 refers to roughly the on-test and off-test ages of animals tested (weights, feed intake, ultrasounds, etc).

16. Utility: mutate_group helpers

I made a few exceptions to my “I hate helper functions” mentality, you can use some mutate_group_*() functions to create groups.

Add sequence numbers, named labels, or concatenated strings within groups:

# Sequential number within each litter (group by dam ID)
pop |>
  get_table("ind_meta") |>
    filter(birth_date == cur_date) |>
  mutate_group_seq(
    group_col  = "id_parent_2",   # group by dam
    new_col    = "piglet_number"
  )

# Named labels within a group
pop |>
  get_table("ind_meta") |>
  mutate_group_named(
    group_col  = "id_parent_2",
    name_col   = "id_ind",
    new_col    = "litter_members"
  )

17. Archive and Restore

Save a replicate’s DuckDB to an archive database for multi-replicate analysis:

archive_replicate(pop, rep = 1L)

[!IMPORTANT] options(tidybreed.replicate) will increase by 1 automatically after each successful archive_replicate() call.

Restore a population from an existing DuckDB file (e.g. to resume a run):

pop <- restore_pop(db_path = "~/path/to/project/tidybreed_output/sim.duckdb")

Database Tables

Table	Rows	Description
`genome_meta`	1 per locus	Locus positions; user chip columns added via `define_chip()`
`founder_haplotypes`	1 per (haplotype × locus)	Haplotype pool sampled by `add_founders()`
`ind_haplotype`	2 per (individual × locus)	Phased haplotypes, long format (paternal / maternal)
`ind_genotype`	1 per (individual × locus)	0/1/2 dosages, long format; on-demand cache filled by `add_dosage()`
`genome_effects`	1 per (locus × trait × effect type)	Additive QTL effect sizes
`ind_meta`	1 per individual	Pedigree, sex, line; user date/status columns added via `mutate_table()`
`ind_phenotype`	1 per (individual × phenotype record)	Long-format phenotype records
`ind_tbv`	1 per (individual × trait)	True breeding values (simulation ground truth)
`ind_true_index`	1 per (individual × index × weight type)	True index values from TBVs
`ind_ebv`	1 per (individual × trait × evaluation)	Estimated breeding values from BLUP/GBLUP
`ind_index`	1 per (individual × index × run)	Computed selection index values
`trait_meta`	1 per genetic trait	Genetic-layer configuration (variance target, units)
`trait_var_comp`	1 per (effect × trait pair)	Additive genetic (co)variance matrix entries
`phenotype_effects`	1 per (phenotype × effect)	Fixed and random model effects
`phenotype_random_effects`	1 per (phenotype × effect × level)	Sampled random effect deviations
`phenotype_meta`	1 per observed phenotype	Observation-layer config (mean, type, sex expression)
`phenotype_components`	1 per (phenotype × component)	Composite phenotype wiring (self/dam/sire/group)
`phenotype_var_comp`	1 per (effect × phenotype pair)	Residual and random effect (co)variance entries
`index_meta`	1 per (index × trait)	Selection index weight definitions

Function Overview

Population & Genome

Function	Purpose
`open_pop()`	Open (or create) a DuckDB-backed population object
`define_genome()`	Build the genome in one transaction: physical map (`genome_meta`), genetic map (`genome_map`), and default chromosome rules (`chr_inheritance`, `chr_recombination`)
`define_chromosome()`	Override inheritance (by offspring sex) or recombination (by parent sex) for one chromosome — sex chromosomes, organelles
`define_founder_haplotypes()`	Generate haplotype pools per line (`uniform`, `fixed`, `beta`, `balding_nichols`, `mosaic`, `gaussian_copula`; `exact_freq` for drift-free frequencies)
`restore_pop()`	Reopen an existing population from a DuckDB file
`close_pop()`	Safely close the DuckDB connection
`print.tidybreed_pop()`	Print population summary

Individuals

Function	Purpose
`add_founders()`	Sample haplotypes and create founder rows in `ind_meta`
`add_offspring()`	Add progeny given a mating-plan `tibble` (1 row per offspring)

Tables & Queries

Function	Purpose
`get_table()`	Return a lazy `dplyr` tbl from any database table
`mutate_table()`	Add or update columns in any table (scalar or vector)
`remove_rows()`	Delete rows from a table by filter (with safety confirmation)
`schema()`	Print all table schemas
`describe_table()`	Print column descriptions for a table
`define_schema_description()`	Register a description string for a user-defined column

Genome & Chips

Function	Purpose
`define_chip()`	Mark filtered loci as members of a named SNP chip
`add_genotypes()`	adds a TRUE within `has_<chip_name>` field of ‘ind_meta’
`extract_genotypes()`	Pull genotypes into R for a chip or set of QTL effects

Traits & Model Configuration

Function	Purpose
`define_trait()`	Register a genetic-layer trait in `trait_meta`
`define_trait_simple()`	Convenience wrapper: `define_trait()` + `define_additive_effects()`
`define_phenotype()`	Register an observed phenotype in `phenotype_meta`
`define_additive_effects()`	Assign QTL effects to filtered loci (single or correlated multi-trait)
`define_effect_cov_matrix()`	Load a (co)variance matrix into `trait_var_comp` or `phenotype_var_comp`
`define_effect_fixed_class()`	Add a discrete fixed-effect level-to-shift mapping
`define_effect_fixed_cov()`	Add a linear regression fixed covariate
`define_effect_random()`	Add a named random effect
`define_effect_intercept()`	Set the overall phenotypic mean (intercept)
`define_residual_cov()`	Write residual (co)variance entries to `phenotype_var_comp`

Simulation Output

Function	Purpose
`add_tbv()`	Compute and store true breeding values in `ind_tbv`
`add_phenotype()`	Sample and store phenotype records in `ind_phenotype`
`add_ebv()`	Run BLUPF90 or parent average; store results in `ind_ebv`
`add_index()`	Compute weighted index from `ind_ebv` (or any table); store in `ind_index`

Selection Index

Function	Purpose
`define_index()`	Register index weights in `index_meta`
`add_index()`	Calculate and store index values from EBVs or TBVs

Group / Litter Utilities

Function	Purpose
`mutate_group_seq()`	Add a within-group sequence number column
`mutate_group_named()`	Add a column with named labels within groups
`mutate_group_concatenate()`	Concatenate values within groups into a string column

Replication & Archiving

Function	Purpose
`archive_replicate()`	Save a replicate’s database to an archive DuckDB
`summary_pop()`	Print a structured summary of the population state

tidybreed