Metadata-Version: 2.4
Name: pygapit-ng
Version: 1.2.4
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.10
Classifier: Programming Language :: Python :: 3.11
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.10
Requires-Dist: matplotlib<4,>=3.10.9
Requires-Dist: numpy<3,>=2.2.6
Requires-Dist: pandas<4,>=2.3.3
Requires-Dist: plotly<8,>=7.0.0
Requires-Dist: scipy<2,>=1.15.3
Requires-Dist: typing-extensions<5,>=4.12; python_version < '3.11'
Provides-Extra: bigdata
Requires-Dist: h5py<4,>=3.16.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,>=2.3.3.260113; extra == 'dev'
Requires-Dist: pyrefly<2,>=1.2.0; 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.6; extra == 'dev'
Requires-Dist: scipy-stubs<2,>=1.15.3.0; extra == 'dev'
Requires-Dist: typing-extensions<5,>=4.12; extra == 'dev'
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.

### 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 | PCA, covariates, supplied kinship alignment, shuffled labels, and missing phenotypes |
| MLM | `GAPIT(model="MLM")` | R-validated top-level workflow and EMMA/P3D statistics | Full MAF-filtered marker set | Variance components, monomorphic markers, missing genotypes, and data alignment |
| CMLM | `GAPIT(model="CMLM")` | R-validated top-level workflow | Full MAF-filtered marker set with fixed 40-group compression | Fixed and automatic compression, native incidence matrices, 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 | Forward/backward selection, final marker statistics, corrected extended BIC, and indefinite-kinship rejection |
| FarmCPU | `GAPIT(model="FarmCPU")` | R-validated complete iterative workflow | Full MAF-filtered marker set | Static-bin selection, 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, 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 | Not yet | 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 | Not yet | 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 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`](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.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.

```bash
pip install pygapit-ng

pip install "pygapit-ng[bigdata]"  # include HDF5 for larger datasets
```

The default installation supports both in-memory analysis and disk-backed
NumPy memory maps. Optional backends are imported only when selected: the
unified storage API silently uses the NumPy backend when `h5py` is unavailable,
while an explicit HDF5 request emits an actionable warning to install
`pygapit-ng[bigdata]` or select `backend="numpy"`.

**Runtime dependencies** are installed automatically: `numpy`, `scipy`,
`pandas`, `matplotlib`, and `plotly`.

### 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 and big-data 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,
    maf_filter,
    open_genotype_store,
    vanraden_kinship,
    write_genotype_store,
)

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

The top-level `GAPIT()` pipeline and the direct `glm_gwas()` and `mlm_gwas()`
APIs accept an open labeled store and read its source data in bounded marker
blocks. Top-level stores supply their own taxa and marker metadata, must contain
finite pre-imputed values, omit `GM`, and currently use VanRaden kinship without
genomic-prediction output. Iterative CMLM, MLMM, FarmCPU, and BLINK scans still
require an in-memory genotype matrix.

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

The automatic writer uses HDF5 when `h5py` is installed and otherwise silently
falls back to a dependency-free NumPy memory-mapped store. Both preserve the
sample-by-marker layout; VanRaden kinship and PCA read both 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 `h5py` is unavailable, the resulting warning 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 NumPy/HDF5.
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_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`](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.
