Skip to content

ROADMAP.md

Latest commit

 

History

History
386 lines (304 loc) · 14.1 KB

ROADMAP.md

File metadata and controls

386 lines (304 loc) · 14.1 KB
Raw

crispr-sieve Architecture Roadmap

Overview

crispr-sieve replaces the heuristic ML classifier in find_edited_reads.py with a richer biological feature set and a better-trained model. It slots into nf-core/scge as a drop-in Nextflow process downstream of the existing indel detection step.

Integration Context: find_edited_reads.py

find_edited_reads.py is the current indel detection script inside nf-core/scge. For WGS DRAGEN tumor-normal data it:

  1. Fetches reads from an edited BAM over target intervals (clustered from a target VCF)
  2. Detects indels via CIGAR ops, SA tags, and soft-clip realignment
  3. Compares against a matched control BAM
  4. Optionally applies a pre-trained .pkl sklearn model (predict_reads_at_position) that extracts 17 read-level features and predicts a CRISPR probability

Limitations driving crispr-sieve:

  • Heuristic filters (--max-in-control, --max-mutation-distance) are too blunt
  • 17-feature model lacks biological priors (repair thermodynamics, off-target cleavage)
  • No strand bias or Phred-weighted evidence
  • Result: false-positive rate too high for clinical scientist trust

Existing 17-feature set (baseline to beat)

Feature Type
read_pair_gap template length
read_insertion CIGAR ins count
read_deletion CIGAR del count
read_mismatch NM tag count
read_softclip CIGAR soft-clip bp
read_del_vs_control delta from control fraction
read_ins_vs_control delta from control fraction
read_mismatch_vs_control delta from control fraction
deletion_exclusive_to_edited binary exclusivity flag
control_has_same_variant binary exclusivity flag
distance_to_closest_pam bp from PAM
is_on_target_site binary
is_at_any_target_site binary
total_indel_size sum ins+del bp
indel_size_category ordinal (0–3)
insertion_to_deletion_ratio float
indel_complexity_score float (capped at 10)

Existing output TSV schema (preserved exactly)

Cluster  Chromosome  Start  End  Ontarget  Gene  indel_type  indel_fraction
indel_allele_fraction  indel_size  indel_bases  num_edited  num_control
total_edited  total_control  significance  prediction  probability  model_info

crispr-sieve appends columns (sieve_score, sieve_label, optional shap_report_path) without touching existing columns, so the downstream scge report requires no changes.

Integration Point in nf-core/scge

nf-core/scge pipeline DAG

DRAGEN tumor-normal
      │
      ▼
find_edited_reads.py          ← existing; outputs indel TSV + flagged BAM regions
      │  (indel_tsv, edited_bam, control_bam, reference_fasta, target_vcf)
      ▼
[crispr-sieve]                ← new Nextflow process
      │  (enriched TSV: original cols + sieve_score + sieve_label)
      │  (optional: shap_report.tsv)
      ▼
Annotate SVs / VEP
      │
      ▼
scge report

Mode B — Full Replacement: crispr-sieve receives the same inputs as find_edited_reads.py and produces an enriched version of the same TSV. The existing .pkl model and predict_reads_at_position function are retired. find_edited_reads.py remains in the repository as the read-only reference implementation and performance baseline.

Tool Stack (WGS-Corrected)

Concern Tool Notes
Read parsing pysam Reuse existing CIGAR/SA/softclip logic verbatim
VCF I/O cyvcf2 Faster than pysam.VariantFile for target parsing
Repair prediction InDelphi + FORECasT Pure Python; sequence-context inputs; WGS-compatible
Off-target scoring CRISPR-Net Subprocess wrapper; returns cleavage probability
ML classifier XGBoost + scikit-learn Binary-compatible with existing .pkl format
Explainability SHAP Native XGBoost support; per-site feature attributions
Pipeline Nextflow DSL2 nf-core module structure; stub mode for dry-run testing
Environments pixi pixi.toml lockfile; multi-feature environments
Testing pytest + synthetic BAM fixtures pysam builds test BAMs with crafted CIGAR strings

Excluded for WGS (amplicon/ECS-only tools — not applicable):

  • CRISPResso2 spatial window quantification
  • CRISPECTOR 2.0 Bayesian posterior
  • ampliCan positional noise map / UMI collapse

Feature Set (WGS-Appropriate)

Tier 1 — Baseline (from find_edited_reads.py, refactored into sieve/features/cigar.py)

All 17 existing features are preserved to maintain backward-compatible model inputs and to allow direct A/B comparison against the existing .pkl model.

Tier 2 — New CIGAR/Phred features (sieve/features/cigar.py)

Feature Description
strand_bias_score Forward vs. reverse read asymmetry at the site; strong asymmetry → PCR artifact
phred_confidence Mean base quality of reads carrying the indel at the indel position
homopolymer_length Length of homopolymer run at/near the indel; ≥4 identical bases = high artifact risk

Tier 3 — Repair priors (sieve/features/repair.py)

Feature Description
predicted_repair_match Top-probability InDelphi/FORECasT outcome vs. observed indel sequence
repair_entropy Shannon entropy of predicted outcome distribution; low entropy = deterministic MMEJ = likely real

Requires sequence fetch from reference FASTA at ±50 bp of the cut site (3–4 bp upstream of PAM).

Tier 4 — Off-target prior (sieve/features/offtarget.py)

Feature Description
off_target_cleavage_prob CRISPR-Net score; biophysically un-cleavable sites provide negative prior

Repository Layout

crispr-sieve/
├── workflow/
│   ├── main.nf                  # Nextflow DSL2 entry point
│   ├── modules/
│   │   └── sieve/
│   │       └── main.nf          # nf-core-compatible process (stub mode supported)
│   └── conf/
│       └── params.config        # default params (mirrors find_edited_reads.py defaults)
├── sieve/                       # installable Python package
│   ├── io/
│   │   ├── tsv.py               # ingest find_edited_reads.py TSV output
│   │   ├── cram.py              # pysam BAM/CRAM reader (CIGAR logic from find_edited_reads.py)
│   │   └── vcf.py               # cyvcf2 wrapper for target VCF
│   ├── features/
│   │   ├── base.py              # FeatureExtractor ABC
│   │   ├── cigar.py             # strand_bias, phred_confidence, homopolymer (+ baseline 17)
│   │   ├── spatial.py           # distance_to_cut_site (precise PAM-relative)
│   │   ├── repair.py            # InDelphi + FORECasT wrappers
│   │   └── offtarget.py         # CRISPR-Net subprocess wrapper
│   ├── labels.py                # 5%/1% threshold labeling + cross-sample hard negatives
│   ├── matrix.py                # assemble feature dicts → DataFrame
│   └── classify.py              # XGBoost + leave-one-donor-out CV + SHAP
├── tests/
│   ├── conftest.py              # synthetic BAM fixtures via pysam
│   ├── fixtures/                # static VCF + guide sequences
│   ├── test_cigar.py
│   ├── test_repair.py
│   ├── test_labels.py
│   └── test_classify.py
├── models/                      # trained model artifacts (.pkl)
├── find_edited_reads.py         # reference implementation (read-only baseline)
├── pixi.toml                    # all environments + task shortcuts
├── pyproject.toml               # sieve package definition
├── ROADMAP.md                   # this file
└── CLAUDE.md

Pipeline Orchestration: Nextflow DSL2

  • Each tool is a module in workflow/modules/
  • stub mode stubs all processes for DAG validation without running tools
  • workflow/modules/sieve/main.nf is structured for eventual nf-core module submission
  • CLI interface mirrors find_edited_reads.py exactly:
process CRISPR_SIEVE {
    tag "$meta.id"
    label 'process_medium'

    input:
    tuple val(meta), path(edited_bam), path(edited_bai)
    tuple val(meta), path(control_bam), path(control_bai)
    path target_vcf
    path reference_fasta
    path indel_tsv          // output of find_edited_reads.py

    output:
    tuple val(meta), path("*.sieve.tsv"), emit: scored_tsv
    tuple val(meta), path("*.shap.tsv"),  emit: shap_report, optional: true

    stub:
    """
    touch ${meta.id}.sieve.tsv
    """
}

Environment Management: pixi

pixi.toml replaces all .yaml environment files. Each feature block is independently installable; environments compose from features.

[project]
name = "crispr-sieve"
channels = ["conda-forge", "bioconda"]
platforms = ["linux-64", "osx-arm64"]

[dependencies]
python = ">=3.11"
pysam = "*"
cyvcf2 = "*"
xgboost = "*"
scikit-learn = "*"
shap = "*"
pandas = "*"
numpy = "*"
pytest = "*"

[feature.repair.dependencies]
# InDelphi + FORECasT; pinned for reproducibility
indelphi-model = "*"

[feature.offtarget.dependencies]
# CRISPR-Net subprocess binary
crispr-net = "*"

[environments]
default = ["repair"]
full = ["repair", "offtarget"]
ci = []   # core only for fast CI

[tasks]
test   = "pytest tests/"
train  = "python -m sieve.classify --train"
score  = "python -m sieve.classify --score"
stub   = "nextflow run workflow/main.nf -stub"

Plugin Architecture

Each feature extractor is a self-contained class implementing the FeatureExtractor ABC:

# sieve/features/base.py
from abc import ABC, abstractmethod
from dataclasses import dataclass

@dataclass
class SiteRecord:
    chrom: str
    start: int
    end: int
    indel_type: str          # from find_edited_reads.py TSV
    indel_fraction: float
    indel_allele_fraction: float
    pam_position: int | None
    guide_sequence: str | None
    reference_context: str | None   # ±50 bp fetched by cram.py

class FeatureExtractor(ABC):
    @abstractmethod
    def extract(self, site: SiteRecord) -> dict[str, float]:
        """Return feature name → value for this site."""
        ...

Extractors are enabled/disabled via pixi.toml feature flags. The matrix.py assembler collects all enabled extractor outputs into a single DataFrame row per site.

Diffusion Levels

Resolution 0 — skeleton
  SiteRecord dataclass
  FeatureExtractor ABC
  Indel detection refactored from find_edited_reads.py into sieve/io/cram.py (verbatim logic)
  Nextflow module skeleton; stub mode works; CLI interface matches find_edited_reads.py exactly
  pixi.toml with all dependencies
  Passthrough classifier (sieve_score = indel_fraction heuristic)
  E2E: nextflow run workflow/main.nf -stub passes; output TSV schema extends find_edited_reads.py
  ─────────────────────────────────────────── add enriched CIGAR features

Resolution 1 — parity + improvement
  strand_bias.py, phred.py, homopolymer.py in sieve/features/cigar.py
  Full 20-feature XGBoost replacing existing .pkl model
  Labeled dataset from ECS cohort (5%/1% threshold via labels.py)
  Leave-one-donor-out CV; ROC-AUC > existing .pkl model = drop-in replacement criterion
  ─────────────────────────────────────────── add biological priors

Resolution 2 — repair priors
  sieve/features/repair.py: InDelphi + FORECasT
  Features: predicted_repair_match, repair_entropy
  Requires reference FASTA fetch at ±50 bp of cut site
  ─────────────────────────────────────────── add cleavage priors

Resolution 3 — off-target scoring
  sieve/features/offtarget.py: CRISPR-Net subprocess wrapper
  Feature: off_target_cleavage_prob
  Negative prior: biophysically un-cleavable sites → artifact
  ─────────────────────────────────────────── optimize + deploy

Resolution 4 — production
  Optuna hyperparameter sweep
  SHAP feature importance report (per-site + aggregate)
  Confidence interval on sieve_score
  Nextflow module ready for nf-core/scge PR

Testing Strategy

Synthetic BAM Fixtures

tests/conftest.py uses pysam to build tiny BAMs with surgically crafted CIGAR strings:

  • Known deletions at known positions relative to mock PAM sites
  • Allows exact assertions on strand_bias_score, homopolymer_length, etc.
  • No real sequencing data required in CI

Mocked Subprocesses

CRISPR-Net is called as a subprocess. monkeypatch returns known probability values so off-target tests are deterministic and fast.

Contract Tests for Labels

tests/test_labels.py asserts specific (edited_frac, control_frac) → label mappings:

  • (0.08, 0.005) → Label 1 (true edit)
  • (0.08, 0.02) → Label 0 (above control threshold)
  • (0.03, 0.005) → Label 0 (below edit threshold)

Stub Nextflow Run

pixi run stub
# validates full DAG without running any tools

Integration Test

One small cohort: 2 synthetic samples, 5 sites each — runs full pipeline end to end using pixi run test.

Baseline Comparison

tests/test_classify.py loads the existing models/site14_site5_combined_model.pkl and asserts that the new XGBoost model achieves higher ROC-AUC on a held-out synthetic dataset.

Domain Knowledge Reference

  • Cas9 cleavage: 3–4 bp upstream of the PAM; indels distant from this coordinate are suspect
  • True edits appear on both strands; strong forward/reverse asymmetry → PCR artifact
  • Homopolymer runs (≥4 identical bases): polymerase slippage mimics real edits
  • MMEJ repair outcomes are thermodynamically predictable from flanking microhomology
  • Multiplexed editing (NS0011: 21 simultaneous targets): chimeric reads and soft-clipping in CIGAR signal chromosomal translocations
  • Low repair_entropy = deterministic MMEJ outcome = strong positive prior for true edit
  • High off_target_cleavage_prob required to call a site real; un-cleavable sites → artifact