Metadata-Version: 2.4
Name: orthoglmm
Version: 0.1.0
Summary: Orthogonal Generalized Linear Mixed Models
Author-email: OrthoGLMM Team <maintainers@orthoglmm.org>
License-Expression: MIT
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.12
Requires-Python: >=3.12
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: dendropy>=4.5
Requires-Dist: numpy<2.3,>=2.1
Requires-Dist: jax<0.5,>=0.4.35
Requires-Dist: jaxlib<0.5,>=0.4.35
Requires-Dist: numpyro<0.16,>=0.15
Requires-Dist: tqdm>=4.67.1
Requires-Dist: pyparsing>=3.2.3
Requires-Dist: matplotlib>=3.10.3
Requires-Dist: scikit-learn>=1.6.1
Requires-Dist: arviz>=0.21.0
Requires-Dist: pyarrow<23,>=22.0.0
Requires-Dist: PyYAML<7,>=6.0
Requires-Dist: scipy<1.15,>=1.14.1
Requires-Dist: statsmodels<0.15,>=0.14.4
Provides-Extra: test
Requires-Dist: pytest<9,>=8.4.2; extra == "test"
Requires-Dist: scipy<1.15,>=1.14.1; extra == "test"
Requires-Dist: statsmodels<0.15,>=0.14.4; extra == "test"
Dynamic: license-file

# OrthoGLMM

Phylogenetic genome-wide association toolkit for orthogroup-level traits. Tests gene-phenotype associations across species while controlling for shared evolutionary history using phylogenetic GLMMs.

## Quick Start

```bash
pixi install          # set up environment
pixi run pytest       # run tests
pixi run smoke-v1     # run the tiny end-to-end workflow check
```

For production analyses, start with `docs/production-guide.md`, then keep
run-specific TOML configs and runbooks with the active analysis project.

## Test Dataset And Reproducible Example

The repository includes a tiny synthetic test dataset at
[`examples/v1_tiny_orthofinder/`](examples/v1_tiny_orthofinder/). It is not
biological data; it is a self-contained fixture for checking that installation,
input parsing, model fitting, solver-null calibration, and output writing work
from start to finish.

Run the method on the test dataset:

```bash
git clone https://github.com/jguhlin/orthoglmm.git
cd orthoglmm
pixi install
pixi run orthoglmm smoke \
  --example-dir examples/v1_tiny_orthofinder \
  --output-dir runs/v1_tiny_smoke \
  --overwrite
```

Equivalent shortcut:

```bash
pixi run smoke-v1
```

Test-dataset inputs:

| File | Description |
|------|-------------|
| `species.txt` | Species order used to align all matrices and phenotype rows. |
| `traits.csv` | Synthetic binary trait `Trait` plus `BUSCO` and `GC` covariates. |
| `covariance.tsv` | Phylogenetic relatedness matrix aligned to `species.txt`. |
| `species_tree.nwk` | Newick species tree used by fixed-trait clade calibration. |
| `OrthoFinder/Results/` | Minimal OrthoFinder-like orthogroup tables used as feature input. |
| `count_matrix.tsv` | Equivalent orthogroup-by-species count matrix for matrix-mode examples. |

Expected outputs under `runs/v1_tiny_smoke/`:

| Output | Description |
|--------|-------------|
| `prep.npz` and `prep_meta.json` | Aligned model arrays and input metadata. |
| `solver/genome_wide_stats.csv` and `.parquet` | Observed orthogroup scan statistics before empirical calibration. |
| `solver/solver_meta.json` and `solver/solver_summary.json` | Solver settings, provenance, and row counts. |
| `null_chunks/chunk0.parquet`, `null_chunks/chunk1.parquet` | Packed solver-null calibration shards. |
| `calibrated_genome_wide_stats.csv` | Final calibrated result table. Use `empirical_p` and `empirical_q` for inference. |
| `calibrated_top_candidates.csv` | Top calibrated candidates for inspection. |
| `calibrated_calibrator.json` and `calibrated_summary.json` | Calibration metadata and compact run summary. |
| `smoke_summary.json` | Commands executed and validation counts for the example run. |

For the expanded step-by-step version of the same run, see
[`docs/v1-quickstart.md`](docs/v1-quickstart.md).

## Production Backend

| Backend | Command | Best for |
|---------|--------|----------|
| **Deterministic solver** | `pixi run solver` | Genome-wide screening (fast, reproducible) |

The deterministic solver is the supported production path for genome-wide scans. The old Bayesian/SVI `orthoglmm.pipeline` runner has been retired; reusable IO and predictor helpers now live in smaller modules.

## Deterministic Pipeline (Recommended)

### One command (easiest)

```bash
pixi run pipeline \
  --species-list species_list.txt \
  --cov-matrix phylo_relatedness_matrix.tsv \
  --phenotype-file phenotypes.tsv \
  --phenotype-columns Eusocial Sociality \
  --species-tree species_tree.nwk \
  --gene-tree-glob 'data/OG*_tree.txt' \
  --preset standard \
  --outdir runs/my_analysis
```

This runs prepare → solve → calibrate → QQ for each phenotype with recommended production defaults (`--auto-rho`, `--calibration-backend solver_permutation`, automatic null scheme selection, and automatic null binning). With `--species-tree`, the automatic calibration path resolves to `fixed_trait_clade` plus hybrid binning. Results in `runs/my_analysis/{trait}/`.

For repeatable runs, put the same options in TOML and override only what changes:

```toml
[inputs]
species_list = "species_list.txt"
cov_matrix = "phylo_relatedness_matrix.tsv"
phenotype_file = "phenotypes.tsv"
phenotype_columns = ["Eusocial", "Sociality"]
species_tree = "species_tree.nwk"

[gene_data]
gene_tree_glob = "data/OG*_tree.txt"

[calibration]
preset = "deep"
solver_null_scheme = "auto"
solver_null_binning = "auto"

[output]
outdir = "runs/my_analysis"
out_format = "both"
```

```bash
pixi run pipeline --config run.toml --outdir runs/my_analysis_retry
```

You can generate a starter config instead of writing it from scratch:

```bash
pixi run pipeline --write-config-template run.toml
```

Use `--dry-run` to print the exact commands without executing them:

```bash
pixi run pipeline --config run.toml --dry-run
```

Presets are `smoke` for quick checks, `standard` for routine runs, and `deep` for production validation depth. Any explicit CLI flag overrides the config file and preset-derived defaults.

### Scaling across SLURM or Nextflow

The one-command pipeline is the easiest local path, but very large analyses
should split the observed scan and solver-null calibration into scheduler
chunks. The current scheduler-friendly path is manifest based:

```bash
# observed scan chunks are ordinary solver chunks
pixi run orthoglmm-solver \
  --prep-npz RUN/work/prep.npz \
  --count-matrix counts.tsv \
  --chunk-index "$SLURM_ARRAY_TASK_ID" \
  --n-chunks 200 \
  --output-dir "RUN/og/chunks/chunk_${SLURM_ARRAY_TASK_ID}" \
  --auto-rho \
  --resume

# merge observed chunks
pixi run orthoglmm merge-parquet \
  RUN/og/chunks/chunk_*/genome_wide_stats.parquet \
  --output RUN/og/genome_wide_stats.parquet \
  --source-column shard

# build null task manifest from the merged observed scan
pixi run orthoglmm make-tasks \
  --task-output-dir RUN/og/tasks \
  -- \
  --genome-wide RUN/og/genome_wide_stats.parquet \
  --prep-npz RUN/work/prep.npz \
  --count-matrix counts.tsv \
  --calibration-backend solver_permutation \
  --species-tree species_tree.nwk \
  --trait-name Eusocial \
  --solver-null-scheme auto \
  --solver-null-binning auto \
  --solver-null-n-chunks 400 \
  --solver-null-ogs 400 \
  --solver-reps-per-bin 10 \
  --solver-null-perms 1000 \
  --solver-null-min-pool 1000 \
  --out-prefix RUN/og/calibrated

# run one null chunk per array task, then merge and apply
pixi run orthoglmm run-tasks \
  --task-manifest RUN/og/tasks/manifest.json \
  --chunk-indices "$SLURM_ARRAY_TASK_ID" \
  --workers 1 \
  --resume

pixi run orthoglmm merge-null \
  --task-manifest RUN/og/tasks/manifest.json \
  --out RUN/og/merged_nulls.parquet

pixi run orthoglmm apply-null \
  --task-manifest RUN/og/tasks/manifest.json \
  --null-input RUN/og/merged_nulls.parquet \
  --out-prefix RUN/og/calibrated
```

`auto` here means the calibration logic chooses the null scheme and binning:
with a species tree it resolves to `fixed_trait_clade` plus `hybrid`; without a
species tree it resolves to `phylo_matched` plus `matched`. It does not yet
submit SLURM jobs or generate a Nextflow DAG by itself. For complete scheduler
templates, see `docs/tutorials/slurm.md` and `docs/tutorials/nextflow.md`.

### Step-by-step (for HPC sharding or fine control)

```bash
# 1. Prepare
pixi run python -m orthoglmm.prepare_data_cli \
  --species-list species_list.txt \
  --cov-matrix phylo_relatedness_matrix.tsv \
  --phenotype-file phenotypes.tsv \
  --phenotype-column Eusocial \
  --out-prefix prep

# 2. Solve
pixi run solver \
  --prep-npz prep.npz \
  --gene-tree-glob 'data/OG*_tree.txt' \
  --auto-rho \
  --ridge-lambda 0.05 \
  --out-csv stats.csv \
  --out-format both

# 3. Calibrate (mandatory for valid p-values)
pixi run python -m orthoglmm.calibrate_results_cli \
  --genome-wide stats.csv \
  --prep-npz prep.npz \
  --calibration-backend solver_permutation \
  --gene-tree-glob 'data/OG*_tree.txt' \
  --species-tree species_tree.nwk \
  --auto-rho \
  --solver-null-ogs 400 \
  --solver-null-scheme auto \
  --solver-null-binning auto \
  --solver-reps-per-bin 10 \
  --solver-null-perms 1000 \
  --solver-null-min-pool 1000 \
  --out-prefix calibrated

# 4. QQ plots
pixi run python -m orthoglmm.qq_plots \
  --summary_csv calibrated_genome_wide_stats.csv \
  --outdir qq_plots
```

### Model orientation: feature response, trait predictor

The exact production solver fits each orthogroup or feature as the species-level
response and tests the phenotype as the focal predictor:

```text
orthogroup presence/count ~ trait + phylogeny + covariates
```

This orientation is intentional. OrthoGLMM asks whether a gene-content feature
is enriched, depleted, expanded, or reduced in trait-positive lineages after
accounting for shared ancestry and nuisance covariates. A positive effect means
the feature is more present or has higher copy burden in trait-positive species;
a negative effect means it is depleted or lower copy in those species.

The opposite orientation,

```text
trait ~ orthogroup presence/count + phylogeny + covariates
```

is also a valid association question, and for a simple binary trait by binary
orthogroup table with no phylogeny or covariates it detects nearly the same
co-occurrence. In real comparative genome scans the two orientations are not
equivalent: the response likelihood changes, phylogenetic random effects are
attached to a different stochastic process, sparse orthogroups can separate a
binary trait model, and calibration must match the fitted statistic.

We use the feature-as-response model because genome-content presence and copy
burden have mode-specific error structures that the solver can model directly,
and because the production null is designed to preserve the observed trait
geometry while clade-shuffling feature response vectors. This gives a defensible
phylogenetic enrichment scan: candidate features are associated with
trait-positive lineages under calibrated empirical nulls. As with GWAS, these
hits are causal candidates, not causal proof; causality requires follow-up
evidence such as timing, dosage/specificity, independent transitions, functional
annotation, or sequence-level validation.

### Calibration is mandatory

Raw p-values from the solver show genomic inflation (lambda_GC >> 1) and are not valid for inference. The solver-permutation calibration corrects this by re-running the solver on matched null inputs and fitting an empirical map from null z-scores to calibrated p-values. Always pass `--calibration-backend solver_permutation` explicitly.

For highly structured binary traits, keep `--species-tree` in the run and leave
`--solver-null-scheme auto --solver-null-binning auto` unless you have a
specific reason to override it. With a species tree, auto resolves to
`fixed_trait_clade` plus `hybrid`: it keeps the phenotype fixed and
clade-shuffles orthogroup response vectors, uses matched representative bins
for presence/absence, and pools copy-number nulls globally because matched
copy-number bins were over-conservative in v1 validation. Final production runs
for whole-genome transitions should target about 10,000 nulls per populated
presence bin; lower-depth runs are useful diagnostics but can saturate tail
metrics at the empirical p-value floor.

### Auto-rho for strong phylogeny

When phylogenetic signal is strong, the default `rho=0.0` (full phylogenetic covariance) can produce such low effective sample size that every OG falls back to ridge GLM without phylogenetic correction. Use `--auto-rho` on both `run_ridge_solver.py` and `calibrate.py` to automatically find the minimum covariance blending that achieves the target n_eff.

## CLI Entry Points

| Command | Purpose |
|---------|---------|
| `pixi run pipeline` / `pixi run orthoglmm pipeline` | **Full pipeline in one command** (prepare → solve → calibrate → QQ) |
| `pixi run solver` / `pixi run orthoglmm-solver` | Deterministic solver |
| `pixi run orthoglmm-filter-verify` | Candidate verification heuristics (specificity, dosage, LOCO) |
| `pixi run orthoglmm-tools` | Utilities, including `doctor`, tree prep, and follow-up commands |
| `pixi run sequence-triage` | Sequence-evolution follow-up prioritization |
| `pixi run node-scan` | Internal-node content association scan |

Compatibility scripts remain under `scripts/`, but new commands should prefer
the package-backed Pixi tasks or `python -m orthoglmm.*` entry points.

## Testing & Benchmarking

```bash
pixi run pytest              # unit tests
pixi run bench-solver        # solver micro-benchmark
```

## Repository Layout

- `orthoglmm/solver_core.py`, `solver_robust.py` — deterministic PQL/ridge/score solver
- `orthoglmm/pipeline.py` — retired compatibility shim for old imports
- `orthoglmm/predictors.py`, `pipeline_io.py` — predictor and input preparation helpers
- `orthoglmm/nulls.py` — Pagel's lambda estimation, null simulation, isotonic calibration
- `orthoglmm/permutation.py` — clade block derivation and within-block permutation
- `orthoglmm/concat_results.py`, `qq_plots.py` — aggregation and QC utilities
- `scripts/` — production CLI scripts for the deterministic DAG
- `tests/` — unit and property-based tests
- `docs/production-guide.md` — current production workflow and validation guide
- `docs/publication-methods.md` — manuscript and reviewer methods notes
