Metadata-Version: 2.4
Name: pygapit-ng
Version: 2.0.0rc1
Summary: Python implementation of GAPIT: Genome Association and Prediction Integrated Tool (GLM, MLM, CMLM, MLMM, FarmCPU, BLINK, gBLUP, cBLUP, sBLUP)
Project-URL: Homepage, https://github.com/ou6326/pygapit
Project-URL: Documentation, https://github.com/ou6326/pygapit#readme
Project-URL: Changelog, https://github.com/ou6326/pygapit/blob/main/CHANGELOG.md
Project-URL: Source, https://github.com/ou6326/pygapit
Project-URL: Bug Tracker, https://github.com/ou6326/pygapit/issues
Author: Yiwu Ouyang, PyGAPIT Contributors
License-Expression: GPL-3.0-or-later
License-File: LICENSE
Keywords: BLINK,BLUP,EMMA,FarmCPU,GWAS,MLMM,SNP,association-study,bioinformatics,genomic-selection,genomics,heritability,kinship,mixed-model,plant-breeding,quantitative-genetics
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.12
Requires-Dist: bokeh<3.11,>=3.8.2
Requires-Dist: datashader<0.20,>=0.17.0
Requires-Dist: holoviews<2,>=1.23.0
Requires-Dist: matplotlib<4,>=3.9.0
Requires-Dist: numpy<3,>=2.2.6
Requires-Dist: pandas<4,>=2.3.0
Requires-Dist: plotly<8,>=4.14.3
Requires-Dist: scipy<2,>=1.15.3
Provides-Extra: bigdata
Requires-Dist: h5py<4,>=3.16.0; extra == 'bigdata'
Requires-Dist: zarr<4,>=3.2.0; extra == 'bigdata'
Provides-Extra: dev
Requires-Dist: basedpyright<2,>=1.39.10; extra == 'dev'
Requires-Dist: black<27,>=26.5.1; extra == 'dev'
Requires-Dist: pandas-stubs<4,>=3.0.5; extra == 'dev'
Requires-Dist: pyrefly<2,>=1.3.1; extra == 'dev'
Requires-Dist: pytest-cov<8,>=7.1.0; extra == 'dev'
Requires-Dist: pytest<10,>=9.1.1; extra == 'dev'
Requires-Dist: rpy2<4,>=3.6.7; extra == 'dev'
Requires-Dist: ruff<0.17,>=0.16.8; extra == 'dev'
Requires-Dist: scipy-stubs<2,>=1.17.1.5; extra == 'dev'
Provides-Extra: styles
Requires-Dist: scienceplots<3,>=2.2.2; extra == 'styles'
Requires-Dist: seaborn<0.14,>=0.13.2; extra == 'styles'
Description-Content-Type: text/markdown

# pyGAPIT — Genome Association and Prediction Integrated Tool (Python)

**Compatibility target:** pyGAPIT currently targets R **GAPIT 3.5**, using the official [GAPIT3.5](https://github.com/jiabowang/GAPIT/tree/GAPIT3.5) by Jiabo Wang & Zhiwu Zhang as its upstream reference. That tag currently resolves to commit [`8d6651c`](https://github.com/jiabowang/GAPIT/tree/8d6651c719484c9f6c844144783dca1e4ef85b3e). 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](MIGRATING.md) 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.

| Model | Public interface | GAPIT 3.5 numerical evidence | Official maize regression | Current boundary coverage |
|---|---|---|---|---|
| GLM | `GAPIT(model="GLM")` | R-validated top-level workflow | Full MAF-filtered marker set | Missing-genotype imputation, phenotype subsetting, MAF filtering, PCA/kinship, covariates, supplied kinship alignment, shuffled labels, final-table ordering, and BH adjustment |
| MLM | `GAPIT(model="MLM")` | R-validated complete preprocessing workflow and EMMA/P3D statistics | Full MAF-filtered marker set | `Middle` missing-genotype imputation, phenotype subsetting, MAF filtering, PCA/kinship, P3D/REML variance components and heritability, final-table ordering and BH adjustment, monomorphic markers, covariates, supplied kinship, and data alignment |
| CMLM | `GAPIT(model="CMLM")` | R-validated complete preprocessing and compression workflow | Full MAF-filtered marker set with fixed 40-group compression | `Middle` missing-genotype imputation, phenotype subsetting, MAF filtering, PCA/kinship, fixed and automatic compression, group membership and compressed kinship, incidence P3D/REML variance components and heritability, final-table ordering and BH adjustment, redundant levels, near-collinear covariates, and invalid designs |
| MLMM | `GAPIT(model="MLMM")` | R-validated top-level workflow | Full MAF-filtered marker set without added covariates | Every forward/backward cofactor state, final marker statistics, corrected extended BIC, and indefinite-kinship rejection |
| FarmCPU | `GAPIT(model="FarmCPU")` | R-validated terminal iterative results | Full MAF-filtered marker set | Static-bin selection, terminal pseudo-QTNs, final p-values, and effects |
| BLINK | `GAPIT(model="BLINK")` | R-validated iterative workflow with the upstream missing-CV BIC call characterized | Full MAF-filtered marker set against the corrected-CV reference | PCA-aware BIC, FDR candidates, LD pruning, terminal zero/one/multiple-QTN paths, and invalid-statistic normalization |
| gBLUP | `gblup()` and `GAPIT(..., prediction_model="gBLUP")` | R-validated direct and prediction workflows | Full EarHT prediction set | BLUE, BLUP, PEV, predictions, and variance components |
| cBLUP | `cblup()` and `GAPIT(model="cBLUP")` | R-validated direct and top-level workflows | Full EarHT prediction set | Compression selection, native-incidence BLUE/BLUP/PEV, predictions, and variance components |
| sBLUP | `sblup()`, `GAPIT(model="sBLUP")`, or a prediction override | R-validated corrected SUPER selection and direct/top-level prediction | Full EarHT prediction set | Stable genomic-bin selection, configurable QTN counts, single-QTN support, BLUE/BLUP/PEV, predictions, and variance components |

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`](tests/cross_language/README.md); the required
R-backed suite runs in [`.github/workflows/ci.yml`](.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.

```bash
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.

```bash
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 default Pixi environment uses the highest supported Python version and
includes development, big-data, and optional plotting-style dependencies:

```bash
pixi install
pixi shell
```

Run the ordinary Python checks without the R reference repository:

```bash
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.

```bash
# 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:

```bash
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:

```bash
pixi run python benchmarks/run_baseline.py --output benchmarks/results/baseline.json
```

See [`benchmarks/README.md`](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

```python
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:**
```r
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.
```python
result = GAPIT(Y=Y, G=hapmap_df, model="BLINK")
```

---

## GWAS models

| Model    | Method type  | Uses kinship | Multi-QTN | Power   | Speed    |
|----------|-------------|-------------|-----------|---------|----------|
| `GLM`    | Single-locus | No (PCs)    | No        | Low     | Fastest  |
| `MLM`    | Single-locus | Yes (global) | No        | Medium  | Fast     |
| `CMLM`   | Single-locus | Compressed  | No        | Medium+ | Fast     |
| `MLMM`   | Multi-locus  | Yes (global) | Yes       | High    | Moderate |
| `FarmCPU`| Multi-locus  | Pseudo-QTN  | Yes       | High    | Moderate |
| `BLINK`  | Multi-locus  | No          | Yes       | High    | Fast     |

```python
# 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

```python
# 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.

```python
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:

```python
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)](https://doi.org/10.3835/plantgenome2011.08.0024).
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`:

| File | Content |
|------|---------|
| `GAPIT.BLINK.EarHT.GWAS.Results.csv` | Full GWAS table: SNP, Chr, Pos, P.value, maf, effect, FDR |
| `GAPIT.BLINK.EarHT.Prediction.csv` | BLUE, BLUP, PEV, GEBV per individual; written only when `buspred=True` and prediction succeeds |
| `GAPIT.EarHT.Kinship.csv` | Selected or supplied kinship matrix |
| `GAPIT.EarHT.PCA.csv` | PC scores per individual |
| `GAPIT.BLINK.EarHT.Manhattan.pdf` | Manhattan plot |
| `GAPIT.BLINK.EarHT.QQ.pdf` | QQ plot with λ annotation |
| `GAPIT.EarHT.Kinship.pdf` | Kinship heatmap |
| `GAPIT.EarHT.PCA.pdf` | 2D PCA scatter |

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:

| R parameter | Python parameter | Default | Description |
|-------------|-----------------|---------|-------------|
| `model` | `model` | `"BLINK"` | Model(s) to run |
| `PCA.total` | `PCA_total` | `3` | Number of PCs as covariates |
| `maf.threshold` | `maf_threshold` | `0.05` | Minimum MAF filter |
| `SNP.impute` | `SNP_impute` | `"middle"` | Missing genotype imputation |
| — | `marker_workspace_mib` | `32.0` | Target MiB for one temporary marker-work matrix |
| `file.output` | `file_output` | `True` | Write result files |
| `cutOff` | `cutOff` | Bonferroni | Significance threshold |
| `LD` | `LD` | `0.7` | LD threshold for BLINK pruning |
| `group.from` | `group_from` | `1` | Min groups for CMLM |
| `group.to` | `group_to` | n | Max groups for CMLM |
| `bin.size` | `bin_size` | `5000000` | Bin size (bp) for FarmCPU |
| — | `super_bin_size` | `10000` | Bin size (bp) for SUPER pseudo-QTN selection |
| — | `super_qtn_counts` | `10,20,...,100` | Candidate pseudo-QTN counts evaluated by REML |
| `h2` | `h2` | `None` | Heritability for simulation |
| `NQTN` | `NQTN` | `None` | QTNs for simulation |
| `buspred` | `buspred` | `False` | Run GS after GWAS |
| `FDRcut` | `FDRcut` | `False` | Use the GAPIT 3.5 FDR cutoff for BLINK pseudo-QTNs |
| `kinship.algorithm` | `kinship_algorithm` | `"VanRaden"` | `"VanRaden"` or `"Zhang"` |
| `Z` | `Z` | `None` | Incidence matrix; combines with `KI` as `Z @ KI @ Z.T` |
| — | `prediction_model` | `None` | Override prediction with `gBLUP`, `cBLUP`, or `sBLUP` |
| `Multiple_analysis` | `Multiple_analysis` | `False` | Write combined Manhattan and QQ plots by trait |

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:

```python
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:

```python
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

```bash
# 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

```python
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.

| Function | HoloViews element | Matplotlib | Bokeh | Plotly | Datashader |
|---|---|---:|---:|---:|---:|
| `manhattan` | Overlay / DynamicMap | Yes | Yes | Yes | `large_data="auto"` or `"aggregate"` |
| `qq_plot` | Overlay | Yes | Yes | Yes | No |
| `kinship_heatmap` | Image | Yes | Yes | Yes | No |
| `pca_plot_2d` | Overlay | Yes | Yes | Yes | No |
| `pca_plot_3d` | Scatter3D | Yes | No | Yes | No |
| `gs_scatter` | Overlay | Yes | Yes | Yes | No |
| `phenotype_distribution` | Overlay | Yes | Yes | Yes | No |

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.

```python
# 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`](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](LICENSE), consistent with the
original R GAPIT license.
