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.

The next major release has deliberate visualization and Python-runtime breaks; see Migrating from pyGAPIT 1.2.4 for direct replacements.

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 and cBLUP and sBLUP cover 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. GAPIT 3.5's FarmCPU and Blink functions return only their terminal seqQTN; iteration.output does not add state history to that return value. Their R evidence therefore freezes terminal iterative states, not every loop state. GAPIT's static-bin FarmCPU path also returns no variance components; pyGAPIT preserves its fixed result fields as vg=ve=h2=0.0 rather than deriving a pseudo-kinship REML fit. Other rows must not be interpreted as official-data parity until a corresponding regression is added.

Storage evidence is reported separately from model evidence. Every listed model is exercised by the backend parity matrix on ndarray, NumPy mmap, HDF5, and Zarr inputs. Those tests prove that storage paths preserve the established in-memory statistics. In addition, the CMLM and cBLUP cross-language workflows run a NumPy mmap store directly against the pinned GAPIT 3.5 result. HDF5 and Zarr currently inherit their R-alignment evidence transitively from the strict backend matrix; they are not presented as independent R executions.

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.12 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 and Zarr for larger datasets
pip install "pygapit-ng[styles]"   # optional Seaborn and SciencePlots styles

The default installation supports both in-memory analysis and disk-backed NumPy memory maps. For automatic writes, the unified storage API prefers HDF5, then Zarr, and finally the dependency-free NumPy backend. Explicit requests for an unavailable optional backend emit an actionable warning to install pygapit-ng[bigdata] or select backend="numpy".

Runtime dependencies are installed automatically: numpy, scipy, pandas, matplotlib, holoviews, plotly, bokeh, and datashader. Install pygapit-ng[styles] to add the optional Seaborn and SciencePlots Matplotlib styles.

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 ".[lint,dev]"

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

The default Pixi environment uses the highest supported Python version and includes development, big-data, and optional plotting-style dependencies:

pixi install
pixi shell

Run the ordinary Python checks without the R reference repository:

pixi run ruff check .
pixi run ruff format --check .
pixi run pyrefly check -p all
pixi run basedpyright
pixi run 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
pixi run pytest tests/cross_language -q

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 ".[lint,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 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 source-read batches used by PCA, VanRaden kinship construction, and direct and top-level GLM/MLM scans, including the MLM scan reused by sBLUP. Wide-matrix PCA and VanRaden accumulate sample-space cross-products without retaining a complete centered genotype matrix, while MLM whitens markers one batch at a time. When PCA has fewer retained markers than samples, it retains the complete centered matrix while that matrix remains within twice the workspace budget; larger tall matrices use an exact two-pass sample-batched accumulation. The setting bounds source blocks, not total process memory: the input genotype, PCA Gram and result arrays, sample-space matrices, and native BLAS allocations remain outside it. Very small budgets still process at least one row or marker.

Disk-backed genotypes use the same chunk-readable interface as in-memory arrays:

from pygapit import (
    compute_pca,
    import_hapmap_genotype_store,
    import_numeric_genotype_store,
    maf_filter,
    open_genotype_store,
    vanraden_kinship,
    write_genotype_store,
)

# Convert numeric GD/GM files directly into a disk-backed store. The input GD
# matrix is read once in bounded sample-row blocks instead of being retained in
# full. The same workspace budget controls each parsed float64 block.
import_numeric_genotype_store(
    "genotype.zarr",
    "genotype.txt",
    "marker-map.txt",
    marker_chunk_size=2048,
    marker_workspace_mib=64,
)

# HapMap stores markers in rows, so its genotype calls are converted in one
# sequential marker-block pass.
import_hapmap_genotype_store(
    "hapmap-store.zarr",
    "genotype.hmp.txt",
    marker_chunk_size=2048,
)

write_genotype_store("genotype-store", genotype)
with open_genotype_store("genotype-store") as store:
    filtered, kept = maf_filter(store, marker_workspace_mib=64)
    kinship = vanraden_kinship(filtered, marker_workspace_mib=64)
    pca = compute_pca(filtered, n_components=3, marker_workspace_mib=64)

Numeric GD files store samples in rows and markers in columns. Their importer therefore streams consecutive sample rows across all markers. Ordinary imputation needs one genotype pass; mean imputation uses a statistics pass and a writing pass so that missing values receive whole-dataset marker means. marker_workspace_mib selects the sample-row batch size, while marker_chunk_size controls the destination store's physical marker chunks.

HapMap files already store markers in rows. Their bounded importer reads marker rows sequentially, numericalizes and imputes each block, then writes its sample-by-marker transpose directly to the destination store. Marker metadata is read separately without loading the complete genotype-call matrix.

Existing labeled stores can also be converted between backends without loading the full genotype matrix:

with open_genotype_store("genotype-store") as store:
    write_genotype_store("genotype.zarr", store, marker_chunk_size=2048)

The top-level GAPIT() pipeline supports GLM, MLM, CMLM, MLMM, FarmCPU, BLINK, gBLUP, cBLUP, and sBLUP with an open labeled store. The direct glm_gwas(), mlm_gwas(), cmlm_gwas(), mlmm_gwas(), farmcpu_gwas(), blink_gwas(), cblup(), select_super_qtns(), and sblup() APIs accept the same chunk-readable input; gblup() consumes the already constructed kinship matrix. cBLUP streams marker blocks while constructing VanRaden kinship, but its complete sample-by-sample kinship matrix remains resident in memory during compression selection. Top-level stores supply their own taxa and marker metadata, must contain finite pre-imputed values, omit GM, and currently use VanRaden kinship. SUPER/sBLUP materializes only the largest requested pseudo-QTN candidate pool and reuses its prefixes; MLMM and FarmCPU materialize only bounded marker batches and their small selected-cofactor designs. BLINK bounds parent-store reads for marker scans and LD pruning, while retaining only its sample-count-limited BIC candidate pool.

maf_filter() returns a GenotypeView when its input is a store, so filtering does not copy the complete genotype matrix. GenotypeView also supports sample subsets and marker reordering; sparse marker selections are translated into contiguous reads from the parent store where possible. Stores that expose a physical marker-chunk width, including HDF5 stores written by pyGAPIT, can coalesce nearby sparse selections within a chunk while limiting the read span to four times the requested unique markers.

Without an explicit backend, the writer prefers HDF5, then Zarr, and silently falls back to a dependency-free NumPy memory-mapped store. An explicit backend="zarr" or a .zarr path selects the optional backend implemented with the Zarr Python 3.x API. pyGAPIT currently writes Zarr format 2 stores for a stable, interoperable on-disk layout. All three backends preserve the sample-by-marker layout; VanRaden kinship and PCA read them in contiguous marker blocks. PCA's exact tall-matrix branch switches to sample batches under memory pressure as described above. A .h5/.hdf5 filename or explicit backend="hdf5" requests HDF5. If an explicitly selected optional backend is unavailable, the error includes the pygapit-ng[bigdata] installation command. Explicit backend="numpy" is always available.

Disk stores read the layout produced by the corresponding writer and check the completion flag and schema version. Reading errors propagate from the backend. Returned marker blocks remain usable after the store is closed. Use store.taxa, store.marker_ids, store.chromosomes, and store.positions for typed NumPy metadata arrays; store.marker_map provides a pandas table.

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,
    mlm_gwas,
    qq_plot,
    save_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 carry a typed notebook backend preference.
manhattan_plot = manhattan(snp_names, chromosomes, positions, result.p_values)
qq = qq_plot(result.p_values)
save_plot(manhattan_plot, "manhattan.pdf")
save_plot(qq, "qq.pdf")

HoloViews plots and output backends

Every plotting function returns a thin Visualization with a typed set of supported backends and a preferred display backend. Use backend="bokeh", "matplotlib", or "plotly" at construction to choose how that object displays in Jupyter. The default is Bokeh for 2D plots and Plotly for 3D PCA. Assigning plot.backend changes only that object and never changes the current global HoloViews backend.

The Bokeh renderer does not implement HoloViews Scatter3D; use Matplotlib or Plotly for pca_plot_3d. Its backend argument, mutable backend property, render(), and save_plot() carry this restriction in their types, and unsupported selections raise at runtime.

Seaborn and SciencePlots remain optional Matplotlib output styles. They are deliberately not plot-construction arguments: use save_plot() when a style preset is wanted. A static filename or style selects Matplotlib automatically. The helper applies the style in a local context and falls back with a warning when the optional pygapit-ng[styles] extra is unavailable. The SciencePlots preset includes no-latex, so a TeX installation is not required.

# One HoloViews specification with a per-object notebook backend.
plot = manhattan(
    snp_names,
    chromosomes,
    positions,
    result.p_values,
    effects=result.effects,
)
plot  # Last expression in Jupyter: interactive Bokeh display.
plot.backend = "plotly"
plot  # The same object now displays with Plotly.

save_plot(plot, "manhattan.html")  # Uses plot.backend.
save_plot(plot, "manhattan.pdf")  # Static suffix selects Matplotlib.
matplotlib_figure = plot.render("matplotlib")

# Optional publication style, applied only while saving this figure.
save_plot(plot, "manhattan-science.pdf", style="science")

Advanced HoloViews composition remains available through plot.specification; ordinary notebook, rendering, and saving workflows do not need to unwrap it.

For 250,000 markers or more, large_data="auto" uses HoloViews and Datashader to retain the maximum -log10(p) value in each output pixel. Exact significant markers are overlaid so their locations and interactive hover data remain available. Use large_data="points" to force exact points or large_data="aggregate" to force bounded rasterization at any size.

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.