orthoglmm 0.1.0

Orthogonal Generalized Linear Mixed Models

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

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/. 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:

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:

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.

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)

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:

[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"

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

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

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

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

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:

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

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

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,

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

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