Workflow Overview¶

This page provides a high-level view of the aa-tRNA-seq pipeline architecture and data flow.

Pipeline Architecture¶

The pipeline is organized into modular Snakemake rule files:

flowchart TB
    subgraph rules[Rule Files]
        REF[aatrnaseq-reference.smk<br/>Reference validation/build]
        A[aatrnaseq-process.smk<br/>Core processing]
        B[aatrnaseq-charging.smk<br/>Charging analysis]
        C[aatrnaseq-qc.smk<br/>Quality control]
        D[aatrnaseq-modifications.smk<br/>Modification calling]
        OR[aatrnaseq-odds-ratios.smk<br/>Odds ratio analysis]
        R[aatrnaseq-report.smk<br/>QC report]
        CL[clean.smk<br/>On-demand cleanup]
        E[demux.smk<br/>WDX + LDX/FDX + EDX<br/><i>conditional</i>]
    end

    subgraph common[common.smk]
        F[Sample parsing]
        G[Helper functions]
        H[Output definitions]
    end

    I[Snakefile<br/>Main entry] --> rules
    I --> common

Complete Pipeline Flow¶

Standard Pipeline (No Demultiplexing)¶

flowchart TB
    subgraph Input
        A[POD5 files<br/>per sample]
    end

    subgraph Processing[aatrnaseq-process.smk]
        B[stage_pod5<br/>Symlink raw POD5s]
        C[rebasecall<br/>Dorado basecalling]
        E[bwa_align<br/>Align, dorado tags carried through]
        CM[calmd<br/>MD/NM tags for the ref]
        F[classify_charging<br/>escpod classify]
        G2[add_adapter_tags<br/>PT tags]
        G3[finalize_bam<br/>Hardlink final BAM]
    end

    subgraph Charging[aatrnaseq-charging.smk]
        H[get_cca_trna<br/>Extract probabilities]
        I[get_cca_trna_cpm<br/>Calculate CPM]
    end

    subgraph QC[aatrnaseq-qc.smk]
        J[base_calling_error<br/>Error metrics]
        K[align_stats<br/>Read statistics]
        L[read_attrition<br/>Where reads were lost]
    end

    subgraph Mods[aatrnaseq-modifications.smk]
        M[bam_to_coverage<br/>Coverage tracks]
        N[modkit_pileup<br/>Modification consensus]
        O[modkit_extract_calls<br/>Per-read mods]
        P[modkit_extract_full<br/>Full export]
    end

    subgraph OddsRatios[aatrnaseq-odds-ratios.smk]
        Q[compute_odds_ratios<br/>Pairwise OR]
    end

    subgraph Report[aatrnaseq-report.smk]
        R[render_combined_qc_report<br/>QC report]
    end

    A --> B --> C --> E --> CM --> F --> G2 --> G3

    G3 --> H --> I
    G3 --> J
    G3 --> K
    G3 --> L
    G3 --> M
    G3 --> N
    G3 --> O
    G3 --> P
    O --> Q
    H --> Q
    J --> R
    H --> R
    I --> R
    K --> R

With Demultiplexing (WarpDemuX + EDX)¶

flowchart TB
    subgraph Input
        A[Pooled POD5 files<br/>per run]
    end

    subgraph WDX[WDX Demultiplexing]
        B[warpdemux<br/>Barcode prediction]
        C[parse_warpdemux<br/>Create mapping]
        D[extract_sample_reads<br/>Per-sample IDs]
        E[split_pod5<br/>Split by WDX barcode]
    end

    subgraph Standard[Standard Pipeline]
        F[rebasecall]
    end

    subgraph EDX[EDX Early Splitting]
        G[detect_edx_adapters<br/>3' adapter ID per read]
        H[extract_edx_read_ids<br/>the read list bwa_align aligns]
    end

    subgraph Downstream[Downstream Processing]
        K[bwa_align → classify_charging → ...]
    end

    A --> B --> C --> D --> E --> F --> G --> H
    H --> K

With Demultiplexing (escapepod LDX/FDX)¶

The escapepod backend writes no split POD5: a per-read .p5s sidecar records each read's barcode call, the run is basecalled once as a whole, and the per-sample uBAM is cut from that afterwards — rejoining the standard pipeline at bwa_align.

flowchart TB
    subgraph Input
        A[Raw POD5 files<br/>per run]
    end

    subgraph LDX[LDX/FDX Demultiplexing]
        B[escapepod_demux<br/>--annotate, LDX axis, writes .p5s]
        C[escapepod_demux_fdx<br/>FDX axis, dual-index samples only]
        D[ldx_run_read_ids<br/>select_demux_reads.py joins axes]
        E[rebasecall_ldx_run<br/>whole run, one dorado pass]
        F[ldx_split_parent_map<br/>dorado split-read parents]
        G[extract_ldx_sample_reads<br/>per-sample read IDs]
        H[split_ldx_ubam<br/>per-sample uBAM]
    end

    subgraph Downstream[Downstream Processing]
        I[bwa_align → calmd → classify_charging → ...]
    end

    A --> B
    A -.-> C
    B --> D
    C -.-> D
    D --> E --> F --> G --> H --> I
    B -.-> G
    C -.-> G

Rule Categories¶

Processing Rules¶

Core data processing from raw signal to classified reads:

Charging Analysis Rules¶

Extract and summarize charging classification:

Quality Control Rules¶

Generate QC metrics and statistics:

Modification Rules¶

RNA modification calling with Modkit:

Odds Ratio Rules¶

Per-tRNA pairwise modification odds ratios:

Report Rules¶

QC report generation:

Reference Rules¶

Validate, build, or trim the alignment reference (aatrnaseq-reference.smk):

Demultiplexing Rules¶

Loaded from demux.smk when either backend is enabled (warpdemux.enabled or ldx.enabled); a ruleorder picks the active backend's version of the outputs they share. See Demultiplexing for the full flow.

WarpDemuX (WDX), signal-based:

escapepod (LDX/FDX), basecall-based:

EDX (3' adapter), either backend:

Maintenance Rules¶

Key Processing Steps¶

1. POD5 Staging¶

stage_pod5 lays a directory of symlinks over a sample's raw POD5 files — nothing is copied. Both consumers of a sample's signal (dorado, escpod classify) take a directory and read it recursively, so a sample pooled from several runs, or a run that keeps pod5_pass/pod5_fail apart, is covered without merging:

1
2
3

run1/pod5_pass/*.pod5  ─┐
run1/pod5_fail/*.pod5  ─┼──► pod5/{sample}/<run>/<pod5_pass|pod5_fail|pod5>/*.pod5 (symlinks)
run2/pod5/*.pod5       ─┘

On the LDX path there is no per-sample staging step at all: rebasecall_ldx_run basecalls the raw run directly, and the per-sample uBAM is cut out afterwards.

2. Basecalling¶

Dorado re-basecalls with:

• Move tables (--emit-moves) required by the charging model
• Modification calling (--modified-bases pseU m5C inosine_m6A)
• High-accuracy model (rna004_sup@v6.0.0)

3. Alignment¶

samtools fastq -T '*' streams the uBAM into bwa mem -C so dorado's tags (move table, MM/ML modbase calls, RG) ride the FASTQ comment onto the aligned records — no FASTQ file is written. BWA MEM runs with RNA-optimized parameters:

• -x ont2d preset for ONT reads
• -K 100000000 pins the input batch regardless of thread count, since the tag-bearing comment is ~13x the read
• -F 2324: unmapped, reverse-strand, secondary and supplementary records removed, leaving primary forward alignments only

4. MD/NM Recomputation¶

calmd recomputes MD/NM against the reference (bwa mem does not emit MD on its own). This runs as its own pass right after alignment.

5. Charging Classification¶

escpod classify analyzes signal at the CCA 3' end:

• Input: POD5 (signal) + aligned BAM (with move tables) + reference + model bundle
• Output: the same BAM records with a cl tag (0-255 score) added
• Threshold: ≥200 = charged (the bundle's own recommended operating point)

6. Adapter Position Tagging¶

The add_adapter_tags rule adds PT tags with adapter boundaries:

• Uses parasail Smith-Waterman alignment
• Detects 5' and 3' adapter positions
• Can infer 5' adapter from alignment position when truncated

Resource Requirements¶

GPU Rules¶

These rules always request a GPU; classify_charging optionally does (see GPU Configuration):

CPU-Intensive Rules¶

Memory-Intensive Rules¶

Configuration Points¶

Key parameters that affect pipeline behavior:

Next Steps¶

• Rules Reference - Detailed rule documentation
• Scripts Reference - Python scripts documentation
• Demultiplexing - WarpDemuX setup guide