pyGAPIT — Genome Association and Prediction Integrated Tool (Python)

Compatibility target: pyGAPIT currently targets R GAPIT 3.5, using the official GAPIT3.5 by Jiabo Wang & Zhiwu Zhang as its upstream reference. That tag currently resolves to commit 8d6651c. This identifies the intended upstream baseline; exact numerical and interface parity has not yet been established. GAPIT v4 compatibility is not currently targeted.

It provides GAPIT-style inputs for the GWAS models GLM, MLM, CMLM, MLMM, FarmCPU, and BLINK, together with gBLUP, cBLUP, and sBLUP genomic-selection functions. It is not a drop-in replacement for every R GAPIT parameter.

Compatibility policy

GAPIT 3.5 is a pinned behavioral reference, not a requirement to reproduce every historical implementation detail. pyGAPIT follows its numerical results when the upstream behavior is statistically valid and well defined. A safer Python behavior is used when GAPIT produces invalid statistics, silently changes the requested model, discards data labels, or relies on an obsolete runtime convention. Each such intentional divergence must be documented and covered by a characterization or regression test.

The GAPIT() entry point currently retains several R-derived parameter names for migration compatibility. New public interfaces and replacement parameters will use descriptive snake_case names, explicit option semantics, and typed result objects. Legacy spellings will be deprecated gradually rather than removed without a transition period.

Compatibility evidence

The table below reports tested behavior, not an estimate of feature completeness. R-validated means that an automated cross-language test compares pyGAPIT with the pinned GAPIT 3.5 source. Python-only means that the interface has ordinary regression tests but does not yet have a direct R comparison. Not yet identifies an evidence gap rather than a known model failure.

The official-data column currently refers to GAPIT's bundled maize diversity panel and the EarHT trait. GLM, MLM, CMLM, MLMM, FarmCPU, and BLINK comparisons cover every marker retained by the shared MAF filter; gBLUP covers the complete set of phenotyped taxa. The CMLM regression uses a fixed 40-group compression; automatic compression selection remains covered by the smaller cross-language workflow test. The MLMM regression omits added covariates so GAPIT's extended-BIC penalty is valid; its broken all-NA seqQTN output for this null-model optimum is characterized, while pyGAPIT returns an empty QTN array. The BLINK reference forwards the already supplied PCA covariates into GAPIT's two BIC calls; the unmodified upstream path is also executed to lock its different QTN selection as an intentional divergence. Other rows must not be interpreted as official-data parity until a corresponding regression is added.

Intentional divergences are tested rather than hidden. They include replacing GAPIT's invalid BLINK NaN/infinite statistics with documented valid outputs, retaining PCA covariates during BLINK BIC selection, rejecting CMLM requests that silently change the model or return a singular-fit sentinel, correcting MLMM's extended-BIC marker penalty and invalid null-model QTN sentinel, and rejecting materially indefinite supplied kinship matrices. The detailed test inventory and divergence rationale are maintained in tests/cross_language/README.md; the required R-backed suite runs in .github/workflows/ci.yml.

SUPER selection is another intentional modernization. GAPIT 3.5 ranks pseudo-QTN sets with the legacy GAPIT.get.LL approximation but fits the final sBLUP with a kinship mixed model. pyGAPIT uses the same EMMA REML objective for selection and final prediction, groups bins by chromosome labels rather than a numeric chromosome-offset encoding, and preserves one-marker matrices that R otherwise drops to vectors.

Installation

Regular use

pyGAPIT requires Python 3.10 or newer. Install the pygapit-ng distribution from PyPI; the Python import package and command-line entry point remain pygapit. R and the GAPIT reference repository are not required.

pip install pygapit-ng

pip install "pygapit-ng[bigdata]"  # include HDF5, Zarr, and Dask for larger datasets

Runtime dependencies are installed automatically: numpy, scipy, pandas, matplotlib, seaborn, plotly, joblib, biopython, and jinja2.

Development

Use an editable installation when changing pyGAPIT itself. The development dependencies provide pytest, Ruff, Pyrefly, BasedPyright, rpy2, and the type stubs. The GAPIT reference repository is still not needed for ordinary Python development.

pip install -e ".[dev]"

pip install -e ".[bigdata]"      # editable install with big-data support
pip install -e ".[dev,bigdata]"  # development tools plus big-data support

The equivalent Pixi environment is:

pixi install -e dev
pixi shell -e dev

Run the ordinary Python checks without the R reference repository:

pixi run -e dev ruff check .
pixi run -e dev ruff format --check .
pixi run -e dev pyrefly check -p all
pixi run -e dev basedpyright
pixi run -e dev pytest tests/test_pygapit.py

GAPIT 3.5 alignment validation

This level is only for maintainers checking numerical behavior against the pinned R GAPIT 3.5 source. It uses the same development environment, but also initializes the GAPIT/ Git submodule. The submodule is not used by the installed pyGAPIT package at runtime.

# Fetch the pinned GAPIT 3.5 reference only when running alignment tests
git submodule update --init --recursive

pixi install -e dev
pixi run -e dev pytest tests/cross_language -q

For one environment containing development and big-data dependencies, use full:

pixi install -e full
pixi shell -e full

If Pixi is unavailable, development and validation require an existing R installation and its MASS package before installing the Python development dependencies:

pip install -e ".[dev]"
pytest tests/cross_language -q

Performance baseline

A deterministic manual benchmark covers PCA, VanRaden kinship, GLM, MLM, FarmCPU, and BLINK without adding machine-dependent timing thresholds to CI:

pixi run -e full python benchmarks/run_baseline.py --output benchmarks/results/baseline.json

See benchmarks/README.md for workload and memory measurement details. Benchmark reports are evidence for comparing commits on the same machine; they are not universal performance guarantees.

Quick start

import pandas as pd
from pygapit import GAPIT

# Load GAPIT-style tabular data
Y = pd.read_csv("mdp_traits.txt", sep="\t")  # phenotype
GD = pd.read_csv("mdp_numeric.txt", sep="\t")  # numeric genotype
GM = pd.read_csv("mdp_SNP_information.txt", sep="\t")  # SNP map

# Select one trait so GAPIT returns one GAPITResult instead of a result dict
result = GAPIT(
    Y=Y,
    GD=GD,
    GM=GM,
    model="BLINK",
    trait="EarHT",
    PCA_total=3,
)

print(result.GWAS.head())  # full GWAS results table
print(f"h²    = {result.h2:.3f}")  # heritability
print(f"λ     = {result.lambda_gc:.3f}")  # genomic inflation factor
print(f"QTNs  = {0 if result.QTNs is None else len(result.QTNs)}")

Comparable R GAPIT call when Y contains the selected trait:

myGAPIT <- GAPIT(Y=myY, GD=myGD, GM=myGM, model="Blink", PCA.total=3)

Input data formats

pyGAPIT accepts GAPIT-style phenotype, numeric-genotype/map, and HapMap inputs:

Phenotype file (Y)

Tab-delimited. First column = Taxa names, remaining columns = trait values.

Taxa    EarHT   dpoll
33-16   64.75   64.5
38-11   69.12   61.0
4226    65.5    59.5

Numeric genotype (GD) + map (GM)

GD: First column = taxa names, remaining = SNP dosages (0/1/2).

taxa        PZB00859.1  PZA01271.1  ...
33-16       2           0           ...
38-11       2           2           ...

GM: Three columns: SNP name, Chromosome, Position (bp).

SNP         Chromosome  Position
PZB00859.1  1           157104
PZA01271.1  1           1947984

HapMap genotype (G)

Standard HapMap format with IUPAC allele codes.

result = GAPIT(Y=Y, G=hapmap_df, model="BLINK")

GWAS models

# Run multiple models simultaneously
result = GAPIT(
    Y=Y, GD=GD, GM=GM, model=["GLM", "MLM", "FarmCPU", "BLINK"], trait="EarHT"
)
# Returns a dict keyed by "EarHT_GLM", "EarHT_MLM", etc.

Genomic selection

# gBLUP — best for polygenic traits
result = GAPIT(Y=Y, GD=GD, GM=GM, model="gBLUP", trait="EarHT")

# cBLUP — compressed-kinship prediction
result = GAPIT(Y=Y, GD=GD, GM=GM, model="cBLUP", trait="EarHT")

# sBLUP — MLM scan followed by SUPER pseudo-QTN selection and prediction
result = GAPIT(
    Y=Y,
    GD=GD,
    GM=GM,
    model="sBLUP",
    trait="EarHT",
    super_bin_size=10_000,
    super_qtn_counts=[10, 20, 40, 60, 80, 100],
)

# Run gBLUP prediction after BLINK
result = GAPIT(Y=Y, GD=GD, GM=GM, model="BLINK", trait="EarHT", buspred=True)

# FarmCPU + buspred uses sBLUP when FarmCPU identifies QTNs;
# otherwise prediction falls back to gBLUP
result = GAPIT(Y=Y, GD=GD, GM=GM, model="FarmCPU", trait="EarHT", buspred=True)

# Direct sBLUP requires pseudo-QTN column indices from a prior GWAS
from pygapit import sblup

prediction = sblup(y, X0, GD_array, qtn_indices=selected_qtns)

# Access prediction results
print(result.Pred)
#      Taxa    BLUE    BLUP     PEV   gBreedingValue  Prediction
# 0  33-16   67.4   -2.65   89.3      -2.65          64.75

Genomic prediction interfaces

The prediction API has three distinct layers:

• GAPIT(...) is the aligned end-to-end GWAS and genomic-selection workflow.
• rrblup, gblup, cblup, and sblup are canonical full-data fits; they are exported from both pygapit and pygapit.gs and return typed result objects.
• cross_validate_rrblup and cross_validate_gblup perform fold-local validation and return every out-of-fold prediction and fold metric.

The uppercase pygapit.models.genomic_prediction.RR_BLUP and GBLUP functions remain compatibility wrappers with their historical (training_gebv, cv_correlation) return shape. New code should use the typed fit and validation interfaces above.

from pygapit import rrblup

fit = rrblup(y, GD_array, lambda_=2.5)
print(fit.intercept, fit.regularization)
print(fit.effects, fit.marker_means)
print(fit.gebv)  # centered genomic component
print(fit.prediction)  # intercept + GEBV

Prediction cross-validation

Evaluate prediction on held-out samples independently of the full-data GAPIT(...).Pred fit:

from pygapit import cross_validate_rrblup, cross_validate_gblup

# y: finite phenotype vector; GD_array: numeric sample-by-marker array.
# All arrays must already share the same sample order.
cv = cross_validate_rrblup(y, GD_array, n_folds=5, seed=42)
print(cv.pearson_r, cv.rmse)
print(cv.predictions, cv.fold_ids)  # one held-out prediction per input sample
print(cv.fold_pearson_r, cv.fold_rmse, cv.regularization)

# Keep each family entirely within one fold.
grouped = cross_validate_rrblup(y, GD_array, n_folds=5, groups=family_ids, seed=42)

# Reuse the same split with a phenotype-independent, precomputed kinship.
comparison = cross_validate_gblup(y, K, fold_ids=cv.fold_ids)

These initial CV APIs fit an intercept only. Remove missing phenotypes and align all arrays before calling. RR-BLUP learns marker imputation means, centering, and REML regularization separately in every training fold; markers entirely missing in training contribute zero. A supplied lambda_ must be positive and is treated as externally fixed, so tuning it from these same held-out phenotypes would invalidate the assessment. gBLUP estimates a GLS intercept and variance components in each fold, treating the supplied kinship as fixed. Do not supply a kinship constructed using held-out phenotype-based marker selection. This API does not evaluate cBLUP/sBLUP selection pipelines.

For gBLUP, "fold-local" describes model fitting over the fixed supplied kinship. If its allele frequencies or centering were estimated from all genotypes, held-out genotype information participated in preprocessing (a transductive evaluation), even without phenotype leakage. Strict inductive evaluation requires training-only preprocessing and construction of both training kinship and test-to-training kinship in each fold. The precomputed-K API cannot reconstruct that pipeline; RR-BLUP provides training-only marker preprocessing directly.

gBLUP checks the full kinship for positive semidefiniteness once before fitting. It uses a symmetric copy and accepts negative eigenvalues only within eps * n * max(spectral_radius, 1), allowing numerical roundoff in rank-deficient kinships. It rejects materially indefinite matrices and preserves caller data.

With seed=None, ungrouped folds follow input order. An integer seed shuffles samples; grouped splits balance group sizes without splitting groups. Groups can be NumPy arrays or Pandas Series, including object/string columns containing non-missing strings or integers. Mixed object ID types and missing IDs are rejected. Explicit integer fold_ids override n_folds and cannot be combined with groups or seed. Every sample receives one prediction; every training fold needs at least two samples. Constant vectors and singleton test folds have undefined Pearson correlation (NaN), while RMSE remains available. Results own read-only copies of their arrays; save the fold IDs to reproduce an exact split.

The older pygapit.models.genomic_prediction.RR_BLUP and GBLUP functions retain their (training_gebv, cv_correlation) return shape and now use this validation machinery. RR-BLUP fits an unpenalized intercept and returns GEBV without it, correcting its previous uncentered behavior. The canonical pygapit.gblup result and top-level GAPIT prediction contract are unchanged.

Full-fit tests independently check centering, the intercept, and marker effects against augmented mixed-model equations in R with GAPIT's REML estimator, covering fixed/estimated penalties and both solver dimensions. The returned centered GEBV also matches the BLUP component of bundled GAPIT.EMMAxP3D for the same centered kernel; adding the intercept matches its BLUE + BLUP output. This is not a claim of direct rrBLUP::mixed.solve package parity.

RR-BLUP uses K = Z_centered @ Z_centered.T / m and lambda = m * delta, consistent with its marker-effect parameterization; see Endelman (2011). When markers outnumber training samples, the ridge solve uses sample space. Inputs and the REML kernel still reside in memory; this is not streaming CV.

Output files

When file_output=True (default), pyGAPIT writes to output_dir:

The returned GAPITResult.output_files records every file written for that analysis. With file_output=False, output_files is None and output_dir is not created.

Parameter reference

The main supported GAPIT-style parameters are:

When Z is supplied, KI represents covariance among the random-effect levels (the columns of Z). NumPy inputs are aligned positionally. DataFrame inputs are aligned by phenotype taxa and random-effect labels, and mismatched labels are rejected. An explicit p_threshold takes precedence over FDRcut during BLINK candidate selection. Multiple-analysis plots join models by SNP, chromosome, and position before drawing them on a shared genomic axis; they are written only when file_output=True.

marker_workspace_mib controls batches used by wide-matrix PCA, VanRaden kinship construction, and direct and top-level GLM/MLM scans, including the MLM scan reused by sBLUP. PCA and VanRaden accumulate sample-space cross-products without retaining a second full centered genotype matrix, while MLM whitens markers one batch at a time. The setting bounds one principal sample-by-marker workspace, not total process memory: the input genotype, sample-space matrices, result arrays, and native BLAS allocations remain outside it. Very small budgets still process at least one marker.

Command-line interface

# Basic GWAS
pygapit --Y traits.txt --GD geno.txt --GM map.txt --model BLINK

# Multiple models, custom output directory
pygapit --Y traits.txt --GD geno.txt --GM map.txt \
        --model GLM MLM BLINK FarmCPU \
        --PCA_total 5 --output_dir results/

# Genomic prediction
pygapit --Y traits.txt --GD geno.txt --GM map.txt --model gBLUP

# Phenotype simulation
pygapit --Y traits.txt --GD geno.txt --GM map.txt \
        --model BLINK --h2 0.7 --NQTN 20

Using individual functions

import numpy as np

from pygapit import (
    blink_gwas,
    bonferroni_threshold,
    build_covariate_matrix,
    compute_pca,
    emma_remle,
    farmcpu_gwas,
    gblup,
    genomic_inflation_factor,
    glm_gwas,
    manhattan_plot,
    mlm_gwas,
    qq_plot,
    vanraden_kinship,
)

# Compute kinship
K = vanraden_kinship(
    GD_array,
    marker_workspace_mib=32.0,
)  # (n, n) VanRaden matrix

# PCA for structure control
pca = compute_pca(
    GD_array,
    n_components=3,
    marker_workspace_mib=32.0,
)
X0 = build_covariate_matrix(pca, n_pcs=3)

# REML variance components
remle = emma_remle(y, X0, K)
print(f"h² = {remle.h2:.3f}")

# Run BLINK GWAS
result = blink_gwas(y, X0, GD_array, max_iterations=10, ld_threshold=0.7)
lam = genomic_inflation_factor(result.p_values)
thresh = bonferroni_threshold(len(result.p_values))
sig = (result.p_values <= thresh).sum()
print(f"λ = {lam:.3f},  {sig} significant SNPs")

# Genomic prediction
gs = gblup(y, X0, K)
print(f"Prediction accuracy (r): {np.corrcoef(y, gs.prediction)[0, 1]:.3f}")

# Plots
manhattan_plot(
    snp_names, chromosomes, positions, result.p_values, save_path="manhattan.pdf"
)
qq_plot(result.p_values, save_path="qq.pdf")

Mathematical models

Mixed Linear Model (MLM)

y = X·β + u + e
u ~ N(0, K·σ²g),   e ~ N(0, I·σ²e)

Variance components estimated by REML via EMMA (Kang et al. 2008): spectral decomposition of K → grid search + Brent's method for optimal δ = σ²e/σ²g. P3D approximation: δ estimated once from null model, fixed for all m SNP tests.

VanRaden Kinship (2009)

K = ZZ' / [2 · Σⱼ pⱼ(1-pⱼ)]
Z = GD - 2p          (column-centered 0/1/2 coding)
p = alternate-allele frequencies

BLINK iteration

Loop until convergence:
  1. GLM-1: sort markers by p-value
             LD-prune candidates (r² > threshold)
             select cofactors by BIC minimization
  2. GLM-2: test all m markers with cofactor set as fixed effects
             → updated p-values

BIC = -2·logL + k·log(n) — replaces expensive REML from FarmCPU.

Henderson's MME (gBLUP)

[X'X        X'Z        ] [β]   [X'y]
[Z'X   Z'Z + δ·K⁻¹     ] [u] = [Z'y]

BLUP = û,   BLUE = X·β̂
PEV  = diag(C⁻¹)ᵤᵤ · σ²g

Citation

Citation metadata for pyGAPIT is available in CITATION.cff. If you use pyGAPIT, please also cite the original GAPIT papers:

• Wang J., Zhang Z. (2021) GAPIT Version 3. Genomics, Proteomics & Bioinformatics https://doi.org/10.1016/j.gpb.2021.08.005
• Huang M. et al. (2019) BLINK. GigaScience https://doi.org/10.1093/gigascience/giy154
• Liu X. et al. (2016) FarmCPU. PLOS Genetics https://doi.org/10.1371/journal.pgen.1005767
• Kang H.M. et al. (2008) EMMA. Genetics 178:1709–1723
• VanRaden P.M. (2009) Kinship. J. Dairy Sci. 91:4414–4423

License

pyGAPIT is distributed under the GNU General Public License v3.0 or later, consistent with the original R GAPIT license.