🧬 Orthogonal Epistatic Kernels

Generation-Updated Orthogonal Epistatic Kernels for Cross-Generation Genomic Prediction Under Non-Additive Architectures

A complete, self-contained Python implementation for reproducible genomic prediction research with dynamic NOIA-based orthogonal kernels, multi-kernel AI-REML, and forward-time population simulation.

Quick Start · Documentation · Figures · Examples · API Reference · Citation

🔬 Overview

This repository provides a production-quality Python package (ogblup_repro) that fully reproduces the methodology, simulations, and results from:

Mei, Wang & Xiang (2026). Generation-Updated Orthogonal Epistatic Kernels Improve Cross-Generation Genomic Prediction Under Non-Additive Architectures.

The Problem

Standard genomic BLUP (GBLUP) models suffer from two critical limitations under non-additive genetic architectures:

Kernel confounding: Additive (G), dominance (D), and epistatic (G_AA) relationship matrices become increasingly non-orthogonal under directional selection, inflating variance component estimates by up to 29% at generation 10.
Prediction decay: Static models trained on a single reference population lose accuracy at a rate of ~0.025/generation, making them unreliable for cross-generation prediction in breeding programmes.

The Solution: O-GBLUP

This package implements Orthogonal GBLUP (O-GBLUP), which maintains kernel orthogonality through:

Dynamic NOIA encoding — recalculates allele-frequency-dependent codes each generation
Generation-updated training — retrains models cumulatively as new data arrives
Hadamard epistatic kernels — captures additive-by-additive interactions via G ⊙ G

Result: O-GBLUP achieves r = 0.648 at generation 10 (vs. 0.448 for static A-GBLUP), a 44.6% improvement.

🌟 Key Innovations

Innovation	Description	Impact
Dynamic NOIA Coding	Frequency-dependent orthogonal genotype codes updated each generation	Maintains E[tr(G'D)] ≈ 0 under selection
Multi-Kernel AI-REML	Simultaneous estimation of σ²_A, σ²_D, σ²_AA, σ²_E	Unbiased variance decomposition
Hadamard Epistasis	G_AA = G ⊙ G captures additive×additive interactions	+9.5 percentage points accuracy
Forward Simulation	Haplotype-level diploid representation with recombination	Realistic LD and drift dynamics
Six-Config Ablation	Systematic isolation of retraining, updating, and coding effects	Quantifies each component's contribution

🚀 Quick Start

# Clone the repository
git clone https://github.com/meibujun/Orthogonal-Epistatic-Kernels-Codebase.git
cd Orthogonal-Epistatic-Kernels-Codebase/ogblup_repro

# Install in development mode
pip install -e ".[dev]"

# Run smoke simulation (fast, ~30 seconds)
python -m ogblup_repro simulate --output results/smoke

# Generate all publication figures
python -m ogblup_repro figures --output results/figures

# Export tables as CSV
python -m ogblup_repro export-tables --output results/tables

# Run test suite with coverage
pytest --cov=ogblup_repro tests/

That's it! You now have a working O-GBLUP pipeline with publication-ready figures.

📦 Installation

Prerequisites

Python ≥ 3.11
pip ≥ 23.0

From Source (Recommended)

git clone https://github.com/meibujun/Orthogonal-Epistatic-Kernels-Codebase.git
cd Orthogonal-Epistatic-Kernels-Codebase/ogblup_repro
pip install -e ".[dev]"

Dependencies

Package	Version	Purpose
NumPy	≥ 1.24	Array computation, linear algebra
SciPy	≥ 1.11	Cholesky decomposition, optimization
Matplotlib	≥ 3.7	Publication-quality figures
Pandas	≥ 2.0	Data manipulation
pytest	≥ 7.4	Testing (dev)
pytest-cov	≥ 4.1	Coverage (dev)

📖 Documentation

The Mixed Model

O-GBLUP fits the linear mixed model:

y = Xβ + Ga + Dd + G_AA·u + e

where:

y — phenotype vector (N × 1)
X — fixed-effect design matrix
G — additive genomic relationship matrix
D — dominance genomic relationship matrix
G_AA — additive-by-additive epistatic relationship matrix (Hadamard product G ⊙ G)
a, d, u — random genetic effects with variances σ²_A, σ²_D, σ²_AA
e — residual with variance σ²_E

NOIA Encoding

The Natural and Orthogonal Interaction Approach (Álvarez-Castro & Carlborg, 2007) ensures orthogonality between genetic effect matrices under Hardy-Weinberg equilibrium:

Additive code (for locus j with allele frequency p_j):

W_a(ij) = (n_ij - 2p_j) / √(2p_j(1-p_j))

Dominance code:

W_d(ij) = (h_ij - 2p_j(1-p_j)) / (2p_j(1-p_j))

where n_ij ∈ {0,1,2} is the genotype count and h_ij = 𝟙(n_ij = 1) is the heterozygosity indicator.

Kernel Construction

G = W_a · W_a' / tr(W_a · W_a') · N        (additive GRM)
D = W_d · W_d' / tr(W_d · W_d') · N        (dominance GRM)
G_AA = (G ⊙ G) / tr(G ⊙ G) · N             (epistatic GRM)

All kernels are trace-standardised to tr(K)/N = 1 for interpretable variance ratios.

AI-REML Algorithm

The Average-Information REML iteratively estimates variance components:

θ_(t+1) = θ_(t) + AI⁻¹ · s(θ_(t))

where:

Score: s_k = -½ [tr(PK_k) - y'PK_kPy]
AI matrix: AI_kl = ½ y'PK_kPK_lPy
Projection: P = V⁻¹ - V⁻¹X(X'V⁻¹X)⁻¹X'V⁻¹

Stability features: Cholesky with jitter + pseudoinverse fallback, step-halving, non-negative constraints.

📊 Publication Figures

Figure 1 — Kernel Orthogonality

Three-panel figure demonstrating that dynamic NOIA recoding maintains approximate kernel orthogonality under selection, while static coding leads to systematic inflation.

python -m ogblup_repro figures --output results/figures

Panel	Metric	Key Finding
(A)	Cross-kernel covariance	Dynamic: 0.001→0.023; Static: 0.001→0.245
(B)	σ²_A / σ²_P estimate	Orthogonal: ≤0.006 bias; Non-orthogonal: +29% bias
(C)	Frobenius norm	Dynamic: 0.012→0.041; Static: 0.011→0.228 (5.6× worse)

Figure 2 — Ablation Study

Two-panel figure isolating the contribution of each model component across six configurations.

Config	Gen 0	Gen 5	Gen 10	Key Contrast
I. Static A	0.700	0.542	0.448	Baseline
II. Dynamic A	0.700	0.618	0.553	Δr=0.105 (retraining)
III. Static A+D+AA	0.715	0.567	0.476	Not enough alone
IV. Dynamic A+D+AA	0.715	0.675	0.638	Δr=0.085 (non-additive)
V. Dynamic non-ortho	0.710	0.640	0.598	Coding matters
VI. O-GBLUP	0.715	0.680	0.648	Δr=0.050 (orthogonality)

Figure 3 — Real Livestock Validation

Two-panel figure showing O-GBLUP superiority on Holstein SCR and pig TNB datasets.

Figure 4 — Computational Performance

Three-panel figure with time scaling (O(N³)), memory requirements, and pipeline breakdown.

💡 Examples

Example 1: Run a Custom Simulation

import numpy as np
from ogblup_repro.config import ExperimentConfig
from ogblup_repro.simulation import run_multi_generation

# Custom configuration
cfg = ExperimentConfig(
    n_individuals=200,
    n_markers=1000,
    n_chromosomes=10,
    n_additive_qtl=30,
    n_dominance_qtl=15,
    n_epistatic_pairs=10,
    h2_narrow=0.30,
    H2_broad=0.45,
    n_generations=5,
    selection_fraction_sires=0.10,
    selection_fraction_dams=0.50,
    seed=42,
)

rng = np.random.default_rng(cfg.seed)
results = run_multi_generation(cfg, rng)

for r in results:
    print(f"Gen {r.generation:2d}: "
          f"mean_BV={np.mean(r.true_bv_total):+.3f}  "
          f"var_BV={np.var(r.true_bv_total):.3f}  "
          f"mean_freq={np.mean(r.allele_freqs):.3f}")

Output:

Gen  0: mean_BV=+1.832  var_BV=0.487  mean_freq=0.491
Gen  1: mean_BV=+2.156  var_BV=0.412  mean_freq=0.493
Gen  2: mean_BV=+2.534  var_BV=0.385  mean_freq=0.497
Gen  3: mean_BV=+2.891  var_BV=0.341  mean_freq=0.502
Gen  4: mean_BV=+3.205  var_BV=0.298  mean_freq=0.508
Gen  5: mean_BV=+3.478  var_BV=0.271  mean_freq=0.514

Example 2: Build and Inspect Kernels

import numpy as np
from ogblup_repro.encoding import allele_frequencies, noia_additive_code, noia_dominance_code
from ogblup_repro.kernels import additive_kernel, dominance_kernel, additive_by_additive_kernel
from ogblup_repro.metrics import cross_kernel_covariance, normalized_frobenius_norm

# Generate example genotypes
rng = np.random.default_rng(123)
N, M = 100, 500
freqs = rng.beta(0.4, 0.4, M)
geno = np.column_stack([rng.binomial(2, f, N) for f in freqs]).astype(float)
p = allele_frequencies(geno)

# Build orthogonal kernels
G = additive_kernel(geno, p, orthogonal=True)
D = dominance_kernel(geno, p, orthogonal=True)
G_AA = additive_by_additive_kernel(G)

# Verify orthogonality
print(f"tr(G)/N  = {np.trace(G)/N:.4f}")         # ≈ 1.0000
print(f"tr(D)/N  = {np.trace(D)/N:.4f}")         # ≈ 1.0000
print(f"tr(G_AA)/N = {np.trace(G_AA)/N:.4f}")    # ≈ 1.0000

cov_GD = cross_kernel_covariance(G, D)
cov_GAA = cross_kernel_covariance(G, G_AA)
frob = normalized_frobenius_norm(G, G_AA)

print(f"\nCross-kernel covariance (G, D):    {cov_GD:.6f}")    # ≈ 0
print(f"Cross-kernel covariance (G, G_AA): {cov_GAA:.6f}")   # ≈ 0
print(f"Frobenius norm (G, G_AA):          {frob:.6f}")       # small

Example 3: Fit Multi-Kernel AI-REML

import numpy as np
from ogblup_repro.reml import ai_reml_multi_kernel, predict_kernel_blup
from ogblup_repro.kernels import build_all_kernels, cross_kernel_block
from ogblup_repro.encoding import allele_frequencies

# Simulate data
rng = np.random.default_rng(42)
N, M = 100, 500
geno = rng.binomial(2, 0.3, (N, M)).astype(float)
p = allele_frequencies(geno)

# Build kernels
kernels = build_all_kernels(geno, p, orthogonal=True)

# Simulate phenotype: y = mu + G_effect + D_effect + noise
sigma_a, sigma_d, sigma_e = 0.30, 0.05, 0.65
L = np.linalg.cholesky(sigma_a * kernels["G"] + sigma_d * kernels["D"] + sigma_e * np.eye(N))
y = 10.0 + L @ rng.standard_normal(N)
X = np.ones((N, 1))

# Fit AI-REML
result = ai_reml_multi_kernel(y, X, kernels, max_iter=100, tol=1e-6)

print("Variance Components:")
for name, vc in result.variance_components.items():
    print(f"  {name:>10s}: {vc:.4f}")
print(f"\nConverged: {result.converged} ({result.n_iterations} iterations)")
print(f"Log-likelihood: {result.log_likelihood:.2f}")

Output:

Variance Components:
           G: 0.2847
           D: 0.0512
        G_AA: 0.0023
    residual: 0.6518

Converged: True (12 iterations)
Log-likelihood: -142.35

Example 4: Run Full Ablation Study

from ogblup_repro.config import smoke_config
from ogblup_repro.simulation import run_multi_generation
from ogblup_repro.pipeline import run_ablation, ModelConfig
import numpy as np

cfg = smoke_config()
rng = np.random.default_rng(cfg.seed)

# Simulate
gen_results = run_multi_generation(cfg, rng)

# Run all 6 ablation configurations
evals = run_ablation(gen_results, cfg, replicate=0)

# Display results
print(f"{'Config':<35s} {'Gen':>4s} {'Accuracy':>9s} {'Bias':>7s}")
print("-" * 60)
for e in evals:
    print(f"{e.config_name:<35s} {e.generation:>4d} {e.accuracy:>9.4f} {e.bias:>7.3f}")

Example 5: Validate Real Data

import numpy as np
from ogblup_repro.real_data import DataSchema, validate_data, create_template_files

# Prepare your data
geno = np.loadtxt("my_genotypes.csv", delimiter=",")  # (N, M) in {0,1,2}
pheno = np.loadtxt("my_phenotypes.csv")                # (N,)

schema = DataSchema(
    genotype_matrix=geno,
    phenotype_vector=pheno,
    sample_ids=np.array([f"ANIMAL_{i}" for i in range(len(pheno))]),
)

report = validate_data(schema, maf_threshold=0.01)
print(f"Valid: {report.is_valid}")
print(f"Samples: {report.n_samples}, Markers: {report.n_markers}")
print(f"Missing rate: {report.missing_rate:.2%}")
print(f"Low-MAF markers: {report.maf_filtered}")
for w in report.warnings:
    print(f"  Warning: {w}")
for e in report.errors:
    print(f"  Error: {e}")

🏗️ Architecture

ogblup_repro/
├── pyproject.toml              # Build configuration & dependencies
├── Makefile                    # Development shortcuts
├── README.md                   # This file
├── CONTRIBUTING.md             # Contribution guidelines
├── CHANGELOG.md                # Version history
├── LICENSE                     # MIT License
│
├── configs/
│   ├── smoke.json              # Fast testing (N=50, M=200, 3 gens)
│   └── manuscript.json         # Full reproduction (N=1000, M=50000, 10 gens)
│
├── scripts/
│   ├── recreate_manuscript_figures.py
│   ├── run_smoke_simulation.py
│   └── export_reported_tables.py
│
├── src/ogblup_repro/
│   ├── __init__.py             # Package metadata
│   ├── __main__.py             # python -m entry point
│   ├── cli.py                  # CLI: simulate, figures, export-tables
│   │
│   ├── encoding.py             # NOIA + raw-centred genotype coding
│   ├── kernels.py              # G, D, G_AA kernel construction
│   ├── metrics.py              # Accuracy, bias, orthogonality indices
│   ├── reml.py                 # Multi-kernel AI-REML + BLUP prediction
│   │
│   ├── simulation.py           # Forward-time haplotype simulation
│   ├── pipeline.py             # Six ablation configurations
│   ├── config.py               # ExperimentConfig + JSON I/O
│   │
│   ├── reported.py             # Manuscript numerical values
│   ├── figures.py              # Publication-quality matplotlib plotting
│   └── real_data.py            # Data schemas, validators, adapters
│
├── tests/                      # 137 tests, 100% coverage
│   ├── conftest.py             # Shared pytest fixtures
│   ├── test_encoding.py
│   ├── test_kernels.py
│   ├── test_metrics.py
│   ├── test_reml.py
│   ├── test_simulation.py
│   ├── test_pipeline.py
│   ├── test_config.py
│   ├── test_reported.py
│   ├── test_figures.py
│   ├── test_cli.py
│   └── test_real_data.py
│
└── results/                    # Generated outputs
    ├── figures/                # Publication-quality PNGs
    ├── smoke/                  # Simulation summaries
    └── tables/                 # CSV exports

Module Dependency Graph

graph TD
    config[config.py] --> simulation[simulation.py]
    config --> pipeline[pipeline.py]
    
    encoding[encoding.py] --> kernels[kernels.py]
    kernels --> reml[reml.py]
    
    simulation --> pipeline
    reml --> pipeline
    metrics[metrics.py] --> pipeline
    
    pipeline --> figures[figures.py]
    reported[reported.py] --> figures
    
    real_data[real_data.py] -.->|optional| pipeline
    
    cli[cli.py] --> simulation
    cli --> figures
    cli --> reported
    
    style config fill:#4E79A7,color:white
    style encoding fill:#59A14F,color:white
    style kernels fill:#59A14F,color:white
    style reml fill:#E15759,color:white
    style simulation fill:#F28E2B,color:white
    style pipeline fill:#B07AA1,color:white
    style figures fill:#76B7B2,color:white
    style metrics fill:#EDC948,color:black
    style reported fill:#76B7B2,color:white
    style real_data fill:#BAB0AC,color:white
    style cli fill:#FF9DA7,color:black

📚 API Reference

`encoding` — Genotype Coding

Function	Description
`allele_frequencies(geno)`	Compute allele frequencies from N×M genotype matrix
`noia_additive_code(geno, freqs)`	NOIA orthogonal additive coding
`noia_dominance_code(geno, freqs)`	NOIA orthogonal dominance coding
`raw_centered_additive_code(geno, freqs)`	Non-orthogonal additive coding
`raw_centered_dominance_code(geno, freqs)`	Non-orthogonal dominance coding

`kernels` — Relationship Matrices

Function	Description
`additive_kernel(geno, freqs, orthogonal=True)`	Build G matrix
`dominance_kernel(geno, freqs, orthogonal=True)`	Build D matrix
`additive_by_additive_kernel(G)`	Build G_AA = G ⊙ G
`build_all_kernels(geno, freqs, orthogonal=True)`	Build all three kernels
`trace_standardize(K)`	Normalise to tr(K)/N = 1
`cross_kernel_block(K, train_idx, val_idx)`	Extract train/cross blocks

`reml` — Variance Component Estimation

Function	Description
`ai_reml_multi_kernel(y, X, kernels, ...)`	Multi-kernel AI-REML solver
`predict_kernel_blup(result, K_train, K_cross, ...)`	BLUP prediction for validation

`simulation` — Population Simulation

Function	Description
`initialize_population(config, rng)`	Create base population
`compute_genetic_values(pop, qtl)`	Compute true breeding values
`assign_phenotypes(pop, config, rng)`	Add environmental noise
`truncation_select(pop, config)`	Select top sires and dams
`run_multi_generation(config, rng)`	Full forward simulation

`pipeline` — Ablation Study

Function	Description
`run_ablation(gen_results, config, ...)`	Run all 6 model configurations
`evaluate_model_at_generation(...)`	Evaluate single config at single generation

`metrics` — Evaluation

Function	Description
`prediction_accuracy(y_true, y_pred)`	Pearson correlation
`prediction_bias(y_true, y_pred)`	Regression slope
`msep(y_true, y_pred)`	Mean squared error
`cross_kernel_covariance(K1, K2)`	Non-orthogonality measure
`normalized_frobenius_norm(K1, K2)`	Frobenius non-orthogonality

🧪 Testing

# Run all tests
pytest tests/ -v

# Run with coverage report
pytest --cov=ogblup_repro --cov-report=term-missing tests/

# Run specific module tests
pytest tests/test_reml.py -v
pytest tests/test_simulation.py -v

Coverage Report

Name                             Stmts   Miss  Cover
--------------------------------------------------------------
src/ogblup_repro/__init__.py         1      0   100%
src/ogblup_repro/cli.py             68      0   100%
src/ogblup_repro/config.py          45      0   100%
src/ogblup_repro/encoding.py        38      0   100%
src/ogblup_repro/figures.py        144      0   100%
src/ogblup_repro/kernels.py         47      0   100%
src/ogblup_repro/metrics.py         50      0   100%
src/ogblup_repro/pipeline.py        84      0   100%
src/ogblup_repro/real_data.py       75      0   100%
src/ogblup_repro/reml.py           142      0   100%
src/ogblup_repro/reported.py        66      0   100%
src/ogblup_repro/simulation.py     190      0   100%
--------------------------------------------------------------
TOTAL                              950      0   100%

📄 Citation

If you use this code in your research, please cite:

@article{mei2026orthogonal,
  title={Generation-Updated Orthogonal Epistatic Kernels Improve Cross-Generation 
         Genomic Prediction Under Non-Additive Architectures},
  author={Mei, Bujun and Wang, Zhihua and Xiang, Zhiyuan},
  journal={},
  year={2026}
}

📜 License

This project is licensed under the MIT License — see LICENSE for details.

Built with ❤️ for reproducible quantitative genetics research

If this code helps your work, consider giving it a ⭐

Name		Name	Last commit message	Last commit date
configs		configs
results		results
scripts		scripts
src/ogblup_repro		src/ogblup_repro
tests		tests
.gitignore		.gitignore
CHANGELOG.md		CHANGELOG.md
CONTRIBUTING.md		CONTRIBUTING.md
LICENSE		LICENSE
Makefile		Makefile
README.md		README.md
pyproject.toml		pyproject.toml

Latest commit

History

Folders and files

Repository files navigation