blupf90 — a Python interface for the BLUPF90 suite

Disclaimer. This package is an independent, unofficial wrapper. It is not affiliated with, endorsed by, or maintained by the University of Georgia Animal Breeding and Genetics group (the authors of the BLUPF90 family of programs). "BLUPF90" refers to the underlying Fortran programs, which are distributed separately by UGA under their own terms — see https://nce.ads.uga.edu/wiki/doku.php. This wrapper is distributed on PyPI as blupf90-wrapper to make the distinction explicit; the Python import name remains blupf90 for brevity, in the same style as PyYAML (import yaml) or beautifulsoup4 (import bs4).

blupf90 is a Python package that lets quantitative geneticists fit animal models with the BLUPF90 family of programs (renumf90, blupf90+) using a familiar R / statsmodels-style formula syntax. You write

import blupf90 as bf

model = bf.univariate(
    "body_weight ~ contemporary_group + sex + cov(age_days)"
    " + ped(animal) + pe(animal)",
    data="phenotypes.txt",
    pedigree="pedigree.txt",
    columns="phenotypes.cols.txt",
)
result = model.fit()
print(result.h2)              # +0.40 +/- 0.02

…and the package handles renum_trait.par rendering, folder layout, local execution or SLURM submission, log parsing, and a tidy result object — so you spend your time on the science, not on shell plumbing.

It is meant for everyday lab work and for reproducible publications: every run leaves a self-contained folder on disk with the par file, data symlinks, BLUPF90 logs, and a parsed summary.

---

Table of contents

1. Features
2. Why this package
3. Installation
4. Package layout
5. Formula grammar
6. Quick start — univariate
7. Quick start — bivariate
8. Batches — many traits, one SLURM array
9. The columns mapping
10. Configuration knobs
11. What gets written on disk
12. Command-line interface
13. Random regression models
14. Testing
15. Compatibility
16. Contributing
17. Citation
18. Licence

---

Features

• R / statsmodels-style formula syntax for univariate, bivariate, and general multi-trait animal models.
• Automatic renum_trait.par rendering with correct column indexing, effect ordering, and multi-trait (CO)VARIANCES blocks.
• OPTION se_covar_function blocks for heritability, repeatability, and genetic / permanent-environment / residual / phenotypic correlations emitted automatically.
• Batch API (univariate_set, bivariate_set, model_set) that submits many models as a single SLURM array job.
• Log parser that turns BLUPF90+'s blup_vc.log into a typed Result object with .h2, .rg, .Va, .Vpe, .Ve, ….
• Random regression (Legendre polynomial) support for longitudinal traits — growth curves, lactation curves, and other time-varying phenotypes. See README_RRM.md.
• Zero runtime dependencies for the parsing / rendering core (NumPy required only for random regression).

---

Why this package

BLUPF90 (Misztal et al., 2002–present) is the de facto standard for variance-component estimation and genetic evaluation in dairy and beef cattle, pigs, poultry, aquaculture, and many other livestock species. Its file-based interface — the renum_*.par parameter file with its many OPTION keywords — is flexible but unforgiving:

• column numbers must match the data file exactly,
• random-effect blocks must be ordered correctly,
• OPTION se_covar_function formulas are written by hand,
• batch submission across many traits or many trait pairs typically ends up as ad-hoc bash + sed.

blupf90 removes that friction by lifting the same workflow into a small, well-typed Python API:

Manual BLUPF90	blupf90 Python
Write renum_trait.par	bf.univariate(formula, ...)
sbatch run.slurm per model	models.submit() (one array job)
grep/awk the log	result.h2, result.rg, .Va, …
Per-trait folder + symlinks	model.setup()

The API is data-agnostic — there is no built-in column map for any particular schema. Bring your own data file, pedigree, and a tiny column-position sidecar, and every formula name is resolved at build time.

---

Installation

The package is pure Python (no compiled extensions) but requires the BLUPF90 executables (renumf90, blupf90+) on PATH for the run / submit / fit actions to do real work. Building the model, inspecting it, and parsing existing logs needs only Python.

# From PyPI (once released):
pip install blupf90-wrapper

# From a checkout of this repository:
pip install .

# Or, for a development install:
pip install -e .[dev]

Note the distribution name is blupf90-wrapper but the import name is blupf90 — you pip install blupf90-wrapper once, then write import blupf90 as bf in every script.

BLUPF90 itself is distributed by the UGA Animal Breeding and Genetics group under its own (free for non-commercial research) licence. Install or module load it separately.

Requirements

• Python ≥ 3.9
• numpy ≥ 1.20 (for the random-regression module only)
• BLUPF90 executables on PATH (for run, submit, fit)

The parsing / rendering core has no third-party runtime dependencies; polars or pandas are only needed if you want to post-process the TSV that ResultSet.write_tsv() produces.

---

Package layout

blupf90/
├── __init__.py    public re-exports + version
├── config.py      Config (paths, columns, SLURM defaults), read_columns()
├── terms.py       Term ABC + FixedClass, Covariate, PedigreeRandom,
│                  DiagonalRandom, LegendreFixed, RRMPedigree,
│                  RRMPermanent, parse_formula()
├── legendre.py    orthonormal Legendre polynomial helpers
├── model.py       Model (formula -> par-file rendering -> run / collect)
├── runner.py      Runner (local), SlurmJob, SlurmArrayJob
├── result.py      VarComp, Result, ResultSet, ModelSet
├── api.py         statsmodels-style functions:
│                    univariate(), bivariate(), from_formula(), fit(),
│                    univariate_set(), bivariate_set(), model_set()
└── __main__.py    `python -m blupf90 <subcommand>` CLI

The split mirrors the conceptual layers of an analysis: configure → specify → render → run → parse. Every layer is independently usable (you can Result.from_log(model, path) on its own to re-parse a finished run, for example).

---

Formula grammar

<response>[ + <response>]  ~  <term>[ + <term>]*

where each <term> is one of

token	semantics	par-file effect
name	cross-classified fixed effect	cross class
cov(name)	linear covariate	cov (regression)
ped(name)	additive genetic via pedigree, Va	random, NUMTYPE 2
pe(name)	diagonal random, permanent environment, Vpe	random, NUMTYPE 1
leg(cov, order)	fixed Legendre-polynomial regression	order + 1 covariates
rrm_ped(cov, order, animal)	additive-genetic random regression	coupled random group
rrm_pe(cov, order, animal)	permanent-environment random regression	coupled random group

• Bivariate models just put two responses on the LHS: "weaning_weight + yearling_weight ~ ...".
• A cov(X) term whose name is also a response is silently dropped (perfect collinearity), and the drop is noted in model.summary().
• At most one ped(...) and one pe(...) per model. The package emits OPTION se_covar_function blocks for heritability, repeatability, and (for bivariate models) the genetic / permanent-environment / residual correlations automatically.
• Random regression (leg, rrm_ped, rrm_pe) is documented separately in README_RRM.md.

---

Quick start — univariate

Example 1: Beef body weight (repeatability model)

A single-trait animal model for body weight measured at multiple ages on the same animal:

import blupf90 as bf

model = bf.univariate(
    "body_weight ~ contemporary_group + sex + age_group"
    " + cov(age_days) + ped(animal) + pe(animal)",
    data="phenotypes.txt",
    pedigree="pedigree.txt",
    columns="phenotypes.cols.txt",     # <-- column-position sidecar
)

print(model.summary())   # effects, function definitions, what will be reported
print(model.par())       # rendered renum_trait.par (no files written)

result = model.fit()     # setup folder + renumf90 + blupf90+ + parse log
print(result.summary())
print("h2:", result.h2)           # VarComp(h2 = +0.40 +/- 0.022)
print("Va:", result.Va)           # VarComp(Va = +95.21 +/- 4.10)
print("Vpe:", result.Vpe)
print("Ve:", result.Ve)

Example 2: Dairy milk yield

Same syntax, different domain — a lactation-mean milk yield model:

model = bf.univariate(
    "milk_yield ~ parity + herd_year_season + dim_class"
    " + cov(dim) + ped(animal) + pe(animal)",
    data="phenotypes.txt",
    pedigree="pedigree.txt",
    columns="phenotypes.cols.txt",
)
result = model.fit()
print(result.h2)                  # VarComp(h2 = +0.30 +/- 0.02)

What's in a Result

Attribute	Type	Meaning
.converged	bool	parsed from "Final Estimates" block
.rounds	int	AI-REML iterations
.n_records	int	observations used
.Va	VarComp	additive genetic variance + SE
.Vpe	VarComp	permanent-environment variance + SE
.Ve	VarComp	residual variance + SE
.h2	VarComp	heritability (univariate) from se_covar_function
.repeatability	VarComp	repeatability (univariate)
.rg	VarComp	genetic correlation (bivariate)
.rpe	VarComp	permanent-environment correlation (bivariate)
.re	VarComp	residual correlation (bivariate)
.h2_t1, .h2_t2	VarComp	heritability of each trait (bivariate)

VarComp is a frozen dataclass with three fields — mean, smean (sample mean), se — and a readable repr like +0.40 +/- 0.022. Convenience attributes such as result.Va_SE are also provided.

---

Quick start — bivariate

A two-trait genetic-correlation model between weaning weight and yearling weight (beef selection index inputs):

model = bf.bivariate(
    "weaning_weight + yearling_weight ~ contemporary_group + sex"
    " + cov(age_days) + ped(animal) + pe(animal)",
    data="phenotypes.txt", pedigree="pedigree.txt", columns="phenotypes.cols.txt",
)

result = model.fit()
print("rg :", result.rg)            # genetic correlation +/- SE
print("rpe:", result.rpe)
print("re :", result.re)
print("h2 :", result.h2_t1, result.h2_t2)

The same syntax applies to any two traits — e.g. milk yield / somatic-cell score in dairy, or backfat / ribeye area in beef carcass evaluation.

---

Batches — many traits, one SLURM array

For a publication-grade analysis you typically fit many traits and many trait pairs in parallel. The batch helpers below build a ModelSet that shares a single Config (one data file, one pedigree, one column map) so the whole batch becomes one sbatch --array submission.

univariate_set — same RHS, one model per trait

Heritability for a panel of production traits:

TRAITS = [
    "birth_weight",     # early growth
    "weaning_weight",   # pre-weaning growth
    "yearling_weight",  # post-weaning growth
    "backfat",          # carcass quality
    "ribeye_area",      # carcass quality
]
RHS = ("contemporary_group + sex"
       " + cov(age_days) + ped(animal) + pe(animal)")

ms = bf.univariate_set(
    traits=TRAITS,
    rhs=RHS,
    data="phenotypes.txt", pedigree="pedigree.txt",
    columns="phenotypes.cols.txt",
    modules=["intel-oneapi", "mkl", "blupf90"],   # for `module load` on SLURM
    root="runs/",
)

job_id = ms.submit()         # one sbatch --array=1-N, one node per trait
# ...wait until `squeue -j <job_id>` is empty...
results = ms.collect()        # parses every blup_vc.log
print(results.summary())
results.write_tsv("heritabilities.tsv")

heritabilities.tsv is a long-format file with one row per (model, derived function) pair plus Va, Va_SE, Vpe, Vpe_SE, Ve, Ve_SE, converged, rounds, n_records repeated on every row for easy pivoting.

bivariate_set — every unordered trait pair

Genetic correlations across the whole trait panel:

ms = bf.bivariate_set(
    traits=TRAITS,                # 5 traits -> C(5, 2) = 10 bivariate models
    rhs=RHS,
    data="phenotypes.txt", pedigree="pedigree.txt",
    columns="phenotypes.cols.txt",
    modules=["intel-oneapi", "mkl", "blupf90"],
    root="runs/",
    # pairs=[("weaning_weight", "yearling_weight"), ...]   # optional subset
)

ms.submit()
results = ms.collect()
results.write_tsv("genetic_correlations.tsv")

Each bivariate model lands in <root>/pair_<a>__<b>/, so univariate and bivariate batches over the same traits never collide.

model_set — arbitrary list of formulas

When models in the batch don't share an RHS, or when mixing univariate and bivariate runs:

ms = bf.model_set([
    "weaning_weight ~ contemporary_group + ped(animal) + pe(animal)",
    "weaning_weight + yearling_weight ~ contemporary_group + cov(age_days)"
    " + ped(animal) + pe(animal)",
    "calving_interval ~ parity + herd_year_season + ped(animal)",
], data="phenotypes.txt", pedigree="pedigree.txt", columns="phenotypes.cols.txt")

---

The columns mapping

Every name used in a formula — responses, fixed effects, covariates, the argument of ped(...) and pe(...) — must appear in the columns mapping with its 1-based position in the data file. Two forms are accepted everywhere:

1. A Python dict (explicit):

columns = {
    "animal":              1,
    "contemporary_group":  2,
    "sex":                 3,
    "age_days":            4,
    "age_group":           5,
    "birth_weight":        6,
    "weaning_weight":      7,
    "yearling_weight":     8,
    "backfat":             9,
    "ribeye_area":        10,
}

2. A path to a sidecar file — parsed by bf.read_columns(...):

• *.cols.txt — two-column name<TAB>index file (recommended; mirrors OPTION fields_passed_to_output output)
• *.csv / *.tsv — header row of the data file itself (column index inferred from header position)
• *.json — {name: int} object

If a formula references a name not in columns, Model.__post_init__ raises a clear error at construction time — before any par file or folder is created.

---

Configuration knobs

Every API entry point accepts the same set of optional keyword arguments, all forwarded to Config.from_paths(...):

keyword	default	meaning
modules	[]	module load line in the SLURM script
slurm	{time=24:00:00, mem=64G, ...}	per-job SLURM resource overrides
missing	-9999	OPTION missing N in the par file
root	current working directory	parent of every per-model folder
em_rounds	50	EM-REML warm-up rounds before AI-REML
maxrounds	500	maximum AI-REML iterations
conv_crit	"1d-10"	convergence criterion (Fortran notation)
tol	"1d-14"	positive-definite tolerance

---

What gets written on disk

For each model in <root>/<model_name>/:

renum_trait.par     # generated by Python
data_clean.txt      # symlink to the data file
pedigree.txt        # symlink to the pedigree file
renum.log           # renumf90 stdout
blup_vc.log         # blupf90+ stdout (parsed by Result.from_log)
renf90.par          # produced by renumf90
renf90.dat          # renumbered data
renaddXX.ped        # renumbered pedigree
solutions           # EBVs / fixed-effect solutions

In <root>/ itself:

logs/               # SLURM .out / .err
_job_<name>.slurm   # single-model SLURM script
_job_models.slurm   # batch / array SLURM script

For univariate batches the model folder is named after the trait (<root>/body_weight/); for bivariate batches it is named <root>/pair_<a>__<b>/.

---

Command-line interface

python -m blupf90 show     <formula> --data ... --pedigree ... --columns ...
python -m blupf90 par      <formula> --data ... --pedigree ... --columns ...
python -m blupf90 setup    <formula> --data ... --pedigree ... --columns ...
python -m blupf90 run      <formula> --data ... --pedigree ... --columns ...
python -m blupf90 submit   <formula> --data ... --pedigree ... --columns ...
python -m blupf90 collect  <formula> --data ... --pedigree ... --columns ...

Every subcommand takes the formula as its single positional argument. --columns accepts any of the sidecar formats listed above. Add --kind univariate (or bivariate) to enforce a runtime check.

For batches use the Python ModelSet API directly — the CLI handles one model at a time on purpose.

---

Random regression models

For longitudinal phenotypes — growth curves in beef, pigs, poultry, and aquaculture; lactation curves in dairy; length-at-age in fish — the classical repeatability model with pe(animal) can inflate heritability estimates by absorbing within-animal autocorrelation into the additive genetic variance. The package provides three random-regression tokens (leg, rrm_ped, rrm_pe) that fit smooth per-animal curves using orthonormal Legendre polynomials:

model = bf.univariate(
    "body_weight ~ contemporary_group + leg(age_days, 2)"
    " + rrm_ped(age_days, 1, animal) + rrm_pe(age_days, 1, animal)",
    data="phenotypes.txt", pedigree="pedigree.txt", columns="phenotypes.cols.txt",
)
result = model.fit()
# h2 is now a function of age_days; reconstruct h2(t) from the raw
# variance components in blup_vc.log (see README_RRM.md section 6).

Full documentation — mathematical background, model choice guidance, worked example, and troubleshooting — is in README_RRM.md.

---

How submit() works

• Model.submit() writes a one-shot SLURM script that cds into <root>/<model>/ and runs renumf90 < renum_trait.par followed by blupf90+ < renf90.par. No Python re-entry is required because the par file is already on disk.
• ModelSet.submit() writes one array script for every model in the set; $SLURM_ARRAY_TASK_ID indexes into a bash array of model names baked into the script. One node per task, all tasks run in parallel, results are gathered with ModelSet.collect() once squeue clears.

Both scripts honour the modules and slurm arguments passed in at construction time.

---

Hot-reload in Jupyter

If you edit blupf90/*.py from a long-lived kernel, hot-reload the package without restarting:

import importlib, sys
for name in list(sys.modules):
    if name == "blupf90" or name.startswith("blupf90."):
        importlib.reload(sys.modules[name])

The result parser is robust to this pattern — it never relies on class-identity comparisons that would silently break after a reload.

---

Testing

Smoke tests live under tests/ and exercise the formula parser, Legendre polynomial engine, par-file renderer, and data-file augmentation logic — no BLUPF90 executables required. Run with:

pip install -e .[dev]
pytest tests/ -v

---

Compatibility

• Tested against RENUMF90 / BLUPF90+ version 1.169 (2024–).
• The generated par files use only the widely-supported OPTION set: sol se, EM-REML, maxrounds, conv_crit, tol, method VCE, se_covar_function, missing, remove_all_missing, alpha_size, and (for random regression) RANDOM_REGRESSION data / RR_POSITION / (CO)VARIANCES_PE.
• The wrapper does not currently emit ssGBLUP (OPTION SNP_file), maternal-genetic (RANDOM animal maternal), or OPTION heterogeneous variance blocks. These can be added by subclassing Term and registering the subclass in terms.TERM_FACTORIES.

---

Contributing

Bug reports, feature requests, and pull requests are welcome. When opening a PR please:

1. add a smoke test under tests/ covering the new behaviour;
2. keep the parsing / rendering core free of runtime dependencies (NumPy is only imported by legendre.py and the RRM setup path);
3. update README.md / README_RRM.md if the public API changes.

The codebase follows PEP 8 with a 100-column line limit and uses ruff for linting.

---

Citation

If you use this package in a publication, please cite the underlying BLUPF90 family of programs. The standard reference is:

Misztal, I., S. Tsuruta, D. A. L. Lourenco, Y. Masuda, I. Aguilar, A. Legarra, and Z. Vitezica (2002–present). Manual for BLUPF90 family of programs. University of Georgia. https://nce.ads.uga.edu/wiki/doku.php

For specific algorithms (e.g. AI-REML, ssGBLUP, single-step) cite the appropriate primary literature as listed in the BLUPF90 manual.

---

Licence

The Python wrapper (this package) is released under the MIT licence (see LICENSE), Copyright (c) 2024-2026 Rajesh Neupane, Juan P. Nani, and Umit Bilginer.

BLUPF90 itself is not included in this distribution and is maintained separately by the University of Georgia Animal Breeding and Genetics group under its own terms. This wrapper is not affiliated with, endorsed by, or maintained by that group. All trademarks and program names belong to their respective owners.