Skip to content

Snakemake Pipeline

For production workloads, leech includes a Snakemake pipeline that handles data preparation, grid search, training, inference, and model comparison on HPC clusters. The pipeline is task-agnostic -- samples, labels, and comparisons all come from your own config -- but the examples on this page use the tRNA charging/amino acid assay that drove its development, since that is the config it has been run against in production.

Pipeline structure

pipeline/
├── config/
│   ├── config.yaml              # Main configuration
│   ├── samples-alpine.yaml      # Sample definitions (passed via --configfile)
│   └── comparisons_*.tsv        # Comparison specs (pairwise, chemical, etc.)
├── workflow/
│   ├── Snakefile                # Main workflow
│   └── rules/                   # Modular rule files
│       ├── common.smk
│       ├── prepare.smk
│       ├── grid_search.smk
│       ├── train.smk
│       ├── inference.smk
│       ├── evaluate.smk
│       ├── compare_models.smk
│       └── ablation.smk         # Channel ablation (opt-in, see below)
└── cluster/                     # Cluster execution profiles
    ├── slurm/
    └── lsf/

Configuration

Edit pipeline/config/config.yaml to define your samples and parameters.

Sample definitions

pod5, bam, and label are the only keys the pipeline rules read; label is an arbitrary string passed straight to leech data prepare --label. Any other key (like amino_acid below) is just metadata for your own bookkeeping -- shown here because it's how the tRNA charging example names its samples.

samples:
  sample_charged_ala_rep1:
    pod5: "data/charged/ala/rep1.pod5"
    bam: "data/charged/ala/rep1.bam"
    label: "charged"
    amino_acid: "Ala"

  sample_uncharged_ala_rep1:
    pod5: "data/uncharged/ala/rep1.pod5"
    bam: "data/uncharged/ala/rep1.bam"
    label: "uncharged"
    amino_acid: "Ala"

Model and training settings

model: "ConvLSTMDwell"
epochs: 50
batch_size: 128
learning_rate: 0.001
early_stopping_patience: 5

Model comparison mode

Train and compare multiple architectures on the same data:

compare_models: true
models_to_compare:
  - "ConvLSTMBase"
  - "ConvLSTMDwell"
  - "TransformerDwell"
  - "ConvOnly"
  - "TCNDwell"
  - "ResNetDwell"

ConvOnly is deprecated (leech.models.model_tier): it has zero references in escapepod-models and a known ONNX-export defect with no fix planned, so get_model("ConvOnly", ...) raises a DeprecationWarning. Drop it from models_to_compare rather than silencing the warning.

leech model optimize searches signal context windows and dwell offset. It does not search learning rate, batch size, or layer sizes — set those directly under the training hyperparameters above.

use_grid_search: true
grid_search:
  context_grid: "200:1000:200"   # symmetric fallback when left/right unset
  left_contexts: [200, 500, 1000]
  right_contexts: [200, 500]
  dwell_offsets: "0"
grid_search_epochs: 20
grid_search_parallel: 1

Each value may be a list or a start:stop:step range string.

Pairwise comparisons

Comparisons are driven by a TSV spec (comparison_spec_file in config.yaml, pointing at one of pipeline/config/comparisons_*.tsv; see Batch comparisons with a TSV spec). The bundled specs compare amino acid identity, since that's the tRNA charging example the pipeline ships configured for -- comparisons_charged_uncharged.tsv, comparisons_all_pairwise.tsv, and comparisons_chemical_properties.tsv. A different task supplies its own TSV of meta_label label_set pairs.

Feature channels

feature_channels in config.yaml is passed as leech data prepare --feature-channels by prepare_chunks and the grid-search re-prepare rules (reprepare_chunks_optimized_pairwise, reprepare_chunks_optimized_architecture). Unset (the default) omits the flag, so every corpus is prepared with leech's own default -- the historical twelve channels (leech.constants.DEFAULT_FEATURE_CHANNELS) -- and existing runs are unaffected. Set it to a YAML list, or the string "all" for every known channel (the default twelve, then level_skew, level_kurtosis):

feature_channels: [dwell, dwell_log, level_mean, level_skew]

This key is independent of the channel-ablation workflow below, which always prepares its own corpus with every channel and selects a subset per arm at train time instead.

Channel ablation

Does dropping (or adding) a feature channel change the model's held-out performance? pipeline/workflow/rules/ablation.smk runs that sweep: it prepares one train/val/test corpus per pairwise comparison with every known feature channel, then trains and evaluates several arms -- each a named subset of channels -- on that same corpus, repeated over several training seeds so a channel's effect can be judged against ordinary run-to-run noise.

It is off by default (use_channel_ablation: false) and not part of rule all/all_compare_models -- it trains len(arms) * len(seeds) models per pairwise comparison on top of the normal run, so enable it deliberately:

use_channel_ablation: true

ablation:
  leave_one_out:   # default: each of the default twelve
    - dwell
    - level_mean
    # ...
  add_one:         # default: shown here
    - [level_skew]
    - [level_kurtosis]
    - [level_skew, level_kurtosis]
  custom: {}       # free-form name -> channel list, e.g.:
    # motif_only: [dwell, dwell_log, level_mean, level_std]
  seeds: [1, 2, 3]

snakemake --profile ../cluster/slurm all_channel_ablation

Arms. baseline (the default twelve) is always included. leave_one_out generates one arm per listed channel, named loo_<channel>, training on the default twelve minus that channel. add_one generates one arm per listed group, named add_<channel>[_<channel>...], training on the default twelve plus that group. custom adds any named channel list verbatim. A name collision between custom and a generated arm, or a leave_one_out channel outside the default twelve, fails the workflow at parse time rather than silently overwriting or ablating nothing.

Shared corpus. The corpus is prepared and merged once per pairwise comparison, with --feature-channels all, under {chunks_dir}/ablation/all_channels/{pair}/... and {chunks_dir}/merged/pairwise/{pair}/ablation/all_channels/{train,val,test}.npz -- every arm and seed reads the same three files and selects its own channel subset at leech model train --feature-channels time, from that corpus's recorded columns. This is deliberate: a corpus per arm would repeat the same large extraction len(arms) times over data that is already the pipeline's biggest storage cost.

Output. {metrics_dir}/ablation/aggregate/channel_ablation_comparison.tsv.gz has one row per (pair, arm, seed), with the channel list, leech eval test's metrics (accuracy, precision, recall, f1, auroc, auprc), and, per (pair, arm), that metric's mean and sd across seeds and its delta against the pair's baseline arm (<metric>_mean, <metric>_sd, <metric>_delta_vs_baseline). channel_ablation_summary.txt is the human-readable version, arms ranked by f1 (or accuracy, if f1 is absent) per pair -- the same primary-metric convention compare_architectures.py uses for architecture comparison.

Reading the table: a <metric>_delta_vs_baseline smaller in magnitude than that arm's own <metric>_sd is not evidence the channel changed the outcome -- it is within the noise the seeds list exists to measure. Look for deltas that are large relative to both arms' sd, and prefer more seeds over trusting a single-seed delta.

Running the pipeline

SLURM clusters

cd pipeline/workflow

# Dry run
snakemake --profile ../cluster/slurm -n

# Execute
snakemake --profile ../cluster/slurm

LSF clusters

# Using the modern executor plugin (recommended)
snakemake --executor lsf --configfile ../config/samples-alpine.yaml \
  --default-resources lsf_queue=gpuqueue --jobs 100

# Or using the profile
snakemake --profile ../cluster/lsf

Local execution

snakemake --cores 8

Pipeline targets

Single model mode (when compare_models: false):

snakemake --profile ../cluster/slurm all_prepare        # Prepare data only
snakemake --profile ../cluster/slurm all_grid_search    # Grid search only
snakemake --profile ../cluster/slurm all_train          # Train models
snakemake --profile ../cluster/slurm all_infer          # Run inference
snakemake --profile ../cluster/slurm all_single_model   # Full single-model analysis

Model comparison mode (when compare_models: true):

snakemake --profile ../cluster/slurm all_prepare                # Prepare data
snakemake --profile ../cluster/slurm all_grid_search_comparison # Grid search all architectures
snakemake --profile ../cluster/slurm all_train_comparison       # Train all architectures
snakemake --profile ../cluster/slurm all_compare_models         # Full comparison

Channel ablation (when use_channel_ablation: true; see Channel ablation above):

snakemake --profile ../cluster/slurm all_channel_ablation

Output structure

results/
├── chunks/{sample}/              # Prepared training data
│   ├── train.json
│   ├── val.json
│   └── test.json
├── models/
│   ├── grid_search/              # Grid search results
│   ├── {comparison}/              # Trained models, one dir per top-level
│   │   ├── model_best.pt         # comparison (e.g. charged_vs_uncharged)
│   │   ├── config.json
│   │   └── metrics.json
│   └── pairwise/{pair}/
├── inference/                    # Prediction BAMs
│   └── {sample}_predictions.bam
└── metrics/                      # Evaluation results
    ├── {comparison}_summary.tsv
    └── pairwise_summary.tsv