SAFARI-IGHJ-Heavy

Reviewer-Hardened Validation Suite for SAFARI-IGHJ

An 11-subcommand validation and analysis extension for the SAFARI-IGHJ pipeline. Provides negative controls, benchmarks, parameter robustness analysis, and publication-ready statistics — without modifying the base pipeline or its outputs.

---

Why This Exists

SAFARI-IGHJ discovers IGHJ gene segments from mammalian genome assemblies. SAFARI-IGHJ-Heavy wraps the base pipeline with reviewer-class controls and validations:

• Negative controls — false-positive rate estimation at two pipeline stages
• Benchmarking — precision/recall on synthetic ground-truth genomes
• Parameter robustness — 36-combination grid search proving result stability
• Pseudogene sensitivity — two-dimension lenient scan showing what strict criteria exclude
• Assembly quality — fragmentation stress testing across N50 levels
• Method comparison — SAFARI vs naive BLAST baseline (+ optional IgBLAST)
• Phylogenetic statistics — bootstrap CI and FDR-corrected correlations (computed, never hardcoded)
• Expression support — RNA-seq coverage parser for candidate loci
• Scope documentation — formal IGHD/IGHV limitation report with citations
• Dependency verification — full environment health check

All original pipeline outputs remain bit-identical — this suite adds validation layers only.

Architecture

safari-ighj (base pipeline)          safari-ighj-heavy (this package)
─────────────────────────────        ────────────────────────────────
 tBLASTn → locus → RSS → classify    11 subcommands wrapping the base:

 Genome ──► run_single() ──► TSVs     safari run          (baseline wrapper)
                  │                   safari benchmark     (precision/recall)
                  │                   safari lenient-scan  (pseudo detection)
                  ▼                   safari negative-control (FP estimation)
          Unchanged outputs           safari method-compare   (vs naive BLAST)
                                      safari robustness-sweep (36-param grid)
                                      safari expression-support (RNA-seq)
                                      safari phylo-stats      (CI + FDR)
                                      safari frag-stress-test  (N50 curve)
                                      safari ighd-lite-scan    (scope report)
                                      safari doctor            (dep checker)

Installation

# 1. Install the base pipeline first
pip install safari-ighj
# Or: git clone https://github.com/jpierrevd/safari-ighj.git && pip install -e safari-ighj/

# 2. Install safari-ighj-heavy
git clone https://github.com/jpierrevd/safari-ighj-heavy.git
cd safari-ighj-heavy
pip install -e .

# 3. Verify installation
safari doctor

External Dependencies

Tool	Required	Used by
BLAST+	Yes	run, benchmark, lenient-scan, negative-control, method-compare, robustness-sweep
Python >= 3.10	Yes	All subcommands
pandas	Yes	All subcommands (auto-installed)
scipy	Yes	phylo-stats (auto-installed)
matplotlib	Yes	frag-stress-test, robustness-sweep (auto-installed)
MAFFT	Optional	Phylogenetic alignment
IQ-TREE	Optional	Phylogenetic tree inference
IgBLAST	Optional	method-compare --igblast
rpy2 + R	Optional	phylo-stats --pgls (PGLS correction)

Run safari doctor to check which dependencies are available.

Quick Start

# Run all subcommands on the bundled toy genome
bash examples/run_example.sh

# Or run individual subcommands:
safari benchmark --outdir results/benchmark
safari doctor

Subcommands

safari run — Baseline Pipeline Wrapper

Invokes the base SAFARI-IGHJ pipeline. Produces identical outputs plus a SHA256 manifest.

safari run \
    --genome /path/to/genome.fasta \
    --species "Bos taurus" \
    --outdir results/bos_taurus

safari benchmark — Precision / Recall

Runs the pipeline on a synthetic toy genome with 5 planted IGHJ genes (3 Functional, 1 ORF, 1 Pseudo). Computes TP, FP, FN, precision, recall, and F1.

safari benchmark --outdir results/benchmark

Output: benchmark_results.tsv, benchmark_summary.json

safari lenient-scan — Two-Dimension Pseudogene Detection

Explores what the strict pipeline excludes by relaxing parameters along two independent dimensions:

1. RSS relaxation — lower BLAST thresholds to find more candidates
2. Coding integrity — classify pseudo-like sequences (stop codons, truncations, missing WG.G motif)

safari lenient-scan \
    --genome /path/to/genome.fasta \
    --species "Bos taurus" \
    --outdir results/lenient_scan

Output: lenient_scan_results.tsv, lenient_scan_summary.md

safari negative-control — False-Positive Estimation

Two-level negative control design:

Level	What it tests	Method
NC1	BLAST query specificity	Full pipeline on dinucleotide-shuffled genome
NC2 (primary)	RSS + classification FP rate	Random genomic windows through downstream pipeline

safari negative-control \
    --genome /path/to/genome.fasta \
    --species "Bos taurus" \
    --outdir results/negative_control

Output: negative_control_results.tsv

safari method-compare — SAFARI vs Naive BLAST Baseline

Compares locus-filtered SAFARI results against all BLAST hits (no locus clustering). Quantifies the candidates excluded by SAFARI's dominant-cluster selection. Optional IgBLAST comparison.

safari method-compare \
    --genome /path/to/genome.fasta \
    --species "Bos taurus" \
    --outdir results/method_compare \
    --igblast  # optional, requires IgBLAST installation

Output: method_compare.tsv, method_compare_summary.md

safari robustness-sweep — Parameter Grid Search

Tests 36 parameter combinations (4 IC thresholds x 3 locus windows x 3 min_pident values) to demonstrate result stability across the parameter space.

safari robustness-sweep \
    --genome /path/to/genome.fasta \
    --species "Bos taurus" \
    --outdir results/robustness_sweep

Output: robustness_sweep.tsv, robustness_heatmap.png

safari expression-support — RNA-seq Coverage Parser

Parses user-provided BEDtools genomecov -bg output and classifies IGHJ candidate loci as expressed, low_coverage, or no_coverage.

safari expression-support \
    --candidates results/run/candidates.tsv \
    --coverage-bed coverage.bedgraph \
    --outdir results/expression_support

Note: SAFARI does not run read alignment. Users provide their own BEDtools coverage file. Non-zero coverage indicates transcriptional activity but is not proof of V(D)J recombination.

Output: expression_support.tsv

safari phylo-stats — Correlation Statistics

Computes Pearson r, Spearman rho, bootstrap 95% CI, Bonferroni correction, and Benjamini-Hochberg FDR from any SAFARI master summary TSV. All values are computed — nothing is hardcoded.

safari phylo-stats \
    --summary-tsv safari_master_summary.tsv \
    --outdir results/phylo_stats \
    --pgls  # optional, requires rpy2 + R

Circular metrics (enrichment_ratio, RSS_integrity) are explicitly excluded with documented reasons.

Output: phylo_stats.tsv, phylo_stats_report.md

safari frag-stress-test — Assembly Fragmentation Sensitivity

Simulates assembly fragmentation at multiple N50 levels (500, 200, 100, 50, 10 kb) and measures recall degradation.

safari frag-stress-test --outdir results/frag_stress_test

Output: frag_stress_results.tsv, frag_stress_plot.png

safari ighd-lite-scan — IGHD/IGHV Scope Report

Generates a formal scope limitation document explaining why IGHD and IGHV gene segments require dedicated pipelines, with citations to published tools (IgDetective, IMPre, IgDiscover).

safari ighd-lite-scan --outdir results/ighd_scan

Output: IGHD_IGHV_SCOPE_LIMITATION.md

safari doctor — Dependency Checker

Reports the status of all 11 runtime dependencies with clear pass/warn/fail indicators.

safari doctor

Running Tests

# Run all 33 tests
pytest tests/ -v

# With canonical dataset (optional — for phylo-stats integration tests)
SAFARI_CANONICAL_TSV=/path/to/safari_master_summary.tsv pytest tests/ -v

All tests use the bundled toy genome (~600 KB) and complete in under 30 seconds.

Project Structure

safari-ighj-heavy/
├── safari/                         # Main package
│   ├── __init__.py
│   ├── cli.py                      # Entry point (argparse, 11 subcommands)
│   ├── toy_data.py                 # Synthetic genome generator (600 KB, 5 planted genes)
│   └── subcommands/
│       ├── run.py                  # Baseline pipeline wrapper
│       ├── benchmark.py            # Precision/recall on ground truth
│       ├── lenient_scan.py         # Two-dimension pseudogene detection
│       ├── negative_control.py     # NC1 (BLAST specificity) + NC2 (pipeline FP)
│       ├── method_compare.py       # SAFARI vs naive BLAST + optional IgBLAST
│       ├── robustness_sweep.py     # 36-parameter grid search
│       ├── expression_support.py   # RNA-seq coverage parser
│       ├── phylo_stats.py          # Bootstrap CI + BH-FDR correlations
│       ├── frag_stress_test.py     # Assembly N50 sensitivity
│       ├── ighd_limitation.py      # IGHD/IGHV scope report
│       └── doctor.py               # Dependency checker
├── tests/                          # 33 unit tests
│   ├── conftest.py                 # Shared fixtures (toy genome, canonical TSV)
│   ├── test_toy_data.py            # Toy genome generation + pseudo detection
│   ├── test_run.py                 # Baseline wrapper
│   ├── test_benchmark.py           # Benchmark metrics
│   ├── test_negative_control.py    # FP estimation
│   ├── test_phylo_stats.py         # Correlation statistics
│   ├── test_doctor.py              # Dependency checker
│   └── test_ighd.py                # Scope report
├── examples/
│   ├── toy_genome.fasta            # Bundled synthetic genome (~600 KB)
│   ├── toy_truth.tsv               # Ground-truth gene positions
│   └── run_example.sh              # Runs all subcommands on toy data
├── pyproject.toml                  # Package configuration
├── CHANGELOG.md
├── LICENSE                         # MIT
└── README.md

Relationship to SAFARI-IGHJ

	safari-ighj	safari-ighj-heavy (this repo)
Purpose	IGHJ gene discovery pipeline	Validation and analysis suite
Modifies base pipeline?	—	No (bit-identical outputs)
CLI	safari-ighj	safari (11 subcommands)
Dependency	Standalone	Requires safari-ighj
Tests	37 unit tests	33 unit tests

Citation

If you use SAFARI-IGHJ-Heavy in your research, please cite:

Viana de Correia, J.P. (2026). SAFARI-IGHJ-Heavy: Reviewer-Hardened Validation Suite for SAFARI-IGHJ. GitHub. https://github.com/jpierrevd/safari-ighj-heavy

Also cite the base pipeline:

Viana de Correia, J.P. (2026). SAFARI-IGHJ: Systematic Analysis of Framework Antibody RSS Immunogenomics. GitHub. https://github.com/jpierrevd/safari-ighj

License

MIT License. See LICENSE for details.

Contributing

Contributions are welcome. Please open an issue to discuss proposed changes before submitting a pull request.