offtarget: optional per-read XC tags on a target-window BAM for IGV r…

conf/modules.config

@@ -70,9 +70,13 @@ process {
     // Python scripts are single-threaded, so they get 1 CPU (the process_medium label gave 4).
 
     // ECS truth arm — long-running (~15 min–2.5 h, no checkpoint) but light on RAM/CPU.
+    // --offtarget_tagged_bam is the exception: it holds a per-read tag map for the whole
+    // target panel until the second (write) pass, which measured ~4 GB on top of the 1.4 GB
+    // baseline for the 1149-interval AAVS1 site14 panel at ~11,000x. Give it room rather
+    // than let an opt-in review artifact OOM the truth arm.
     withName: 'ECS_INDELS' {
         cpus   = { 1 }
-        memory = { check_max( 8.GB  * task.attempt, 'memory' ) }   // was process_highmem = 64 GB
+        memory = { check_max( (params.offtarget_tagged_bam ? 24.GB : 8.GB) * task.attempt, 'memory' ) }
         time   = { check_max( 8.h   * task.attempt, 'time'   ) }   // longest observed run ~2.5 h
     }
 

docs/OFFTARGET_WORKFLOW.md

@@ -55,6 +55,63 @@ no IGV needed. Off by default (`offtarget_snapshots = false`) because it renders
 pile up at the cut site; right (unedited normal) — clean. Pass `--snapshots` to `run_offtarget.sh`
 to turn this on, so the confirmed off-targets (PLCB2, CNNM3) get the same tumor/normal packet.*
 
+### Read-level tags for IGV (optional)
+
+`<id>.offtarget_analysis.tsv` gives you `indel_fraction = 0.42` and leaves you to work out by eye
+which reads made up the 0.42. Add `--offtarget_tagged_bam true` and the ECS caller also writes
+`<id>.tagged.bam` (+ `.bai`), in which **every read carries an `XC` string tag naming how the caller
+classified it** — so the pileup shows you its reasoning directly.
+
+In IGV: load the BAM → right-click the track → **Color alignments by** → **tag** → type `XC`.
+
+| tag | meaning |
+|---|---|
+| `Edited_Deletion_<N>bp` | deletion of N bp called from the CIGAR or a split alignment |
+| `Edited_Insertion_<N>bp` | insertion of N bp |
+| `Edited_Duplication_<N>bp` | tandem duplication from a split alignment |
+| `Edited_BND_<chrom>` | breakend — the read's other end maps to `<chrom>` (e.g. the transgene contig) |
+| `Edited_SoftClip` | event recovered by realigning a soft clip |
+| `Unedited_WT` | spans the target cleanly, no event — this is the denominator |
+| `Skipped_Duplicate` / `_LowMapQ` / `_Mismatches` / `_Secondary` / `_Supplementary` / `_Unmapped` | excluded by a read filter, shown so you can see *why* it isn't counted |
+| `Skipped_Unevaluable` | a classifier ran but found nothing near enough to a PAM to call |
+| `Skipped_NoSpan` | inside the padded window but doesn't span the target and carries no event |
+
+Each alignment appears **exactly once**, so the depth IGV shows is real. Where a read falls in two
+overlapping target windows and gets classified differently, the most specific call wins — the table
+above is in precedence order, top to bottom.
+
+**Don't expect the tag counts to equal the TSV columns** — they count different things, and the BAM
+is the more literal of the two:
+
+- tags are per **alignment record**, while `indel_reads` is per **fragment** (the caller collapses
+  R1/R2 by read name) and is taken *after* the site-level filters. So `Edited_*` records run higher
+  than `indel_reads` — roughly 2× where both mates cover the cut site. On the AAVS1 site14
+  on-target that is 17,396 `Edited_*` records against `indel_reads` = 10,821.
+- the BAM covers the whole padded window, so it also holds reads that never entered `total_reads`.
+  Across the AAVS1 site14 panel that is 65% `Unedited_WT`, 35% `Skipped_*` and 0.7% `Edited_*`; at
+  a *1 bp* target specifically, the ±150 bp pad means most records are `Skipped_NoSpan` (69% at the
+  on-target). They are kept deliberately — a pileup cropped to reads spanning a single base is
+  unreadable in IGV.
+
+Use the tags to see *which* reads drove a call and why; use the TSV for the number.
+
+**Off by default, and worth keeping that way for routine runs.** Output is restricted to the target
+windows (target ±150 bp) rather than the whole CRAM, but it still scales with target count and
+depth. Measured on one AAVS1 site14 ECS sample (1,149 target intervals → 1,145 merged windows,
+~11,000× at the on-target):
+
+| | without | with `--offtarget_tagged_bam` |
+|---|---|---|
+| output | — | **800 MB** `.bam` + 1.5 MB `.bai`, 21.8 M reads |
+| peak RSS | 1.4 GB | **6.3 GB** (the tag map is held until the write pass) |
+| wall clock | ~30 min | ~46 min |
+
+`ECS_INDELS` is given 24 GB instead of 8 GB when the flag is set, so enabling it does not OOM the
+truth arm. Still check the size on one sample before turning it on across a cohort — and note the
+debug log this pipeline learned that lesson from, `offtarget_ecs_unevaluable_log`, reached
+0.1–1 TB/sample and filled the work directory. Restricting the run with `--regions` keeps both
+numbers small when you only want to review a handful of loci.
+
 ## What runs depends on the samplesheet
 
 The `datatype` column decides:
@@ -115,6 +172,7 @@ You hand it one samplesheet with these columns:
 | `offtarget_germline_max_ctrl_if` | 0.05 | matched-normal indel frac above this = germline/artifact, not a somatic edit (`label` 0) |
 | `offtarget_hotspot_pad` | 25 | bp window to match a worklist hit to a predicted hotspot |
 | `offtarget_snapshots` | false | render IGV-style pileup PNGs for LIKELY EDITs |
+| `offtarget_tagged_bam` | false | emit `<id>.tagged.bam` with per-read `XC` tags for IGV review (above) |
 | `offtarget_rescue` | true | high-evidence rescue (below); `false` = call on model score alone |
 | `offtarget_rescue_min_ifrac` | 0.15 | rescue: minimum indel fraction in the edited sample |
 | `offtarget_rescue_min_conc` | 0.5 | rescue: minimum positional concordance (clonality) |

modules/local/ecs_indels.nf

@@ -13,20 +13,26 @@ process ECS_INDELS {
     output:
     tuple val(meta), path("${meta.id}.offtarget_analysis.tsv"), emit: indels_file
     tuple val(meta), path("${meta.id}.offtarget_edits.vcf"),    emit: indels_vcf
+    tuple val(meta), path("${meta.id}.tagged.bam*"), optional: true, emit: tagged_bam
     path "versions.yml", emit: versions
 
     script:
     // The per-read "unevaluable reads" log is a debug artifact that is NOT an emitted
     // output and is not consumed downstream, yet it grows to ~0.1-1 TB per sample and
     // was the sole cause of multi-TB work-dir bloat / ENOSPC. Off unless explicitly asked.
     def unevaluable = params.offtarget_ecs_unevaluable_log ? "-u ${meta.id}.unevaluable_reads.txt" : ""
+    // Review aid, not a pipeline input: a BAM of the target windows in which every read
+    // carries an XC tag naming the per-read call, for colouring the pileup in IGV. Window-
+    // restricted, but it still scales with target count -- off unless explicitly asked.
+    def tagged_bam = params.offtarget_tagged_bam ? "--tagged-bam-out ${meta.id}.tagged.bam" : ""
     """
     python ${projectDir}/bin/find_edited_reads.py \\
         --fasta ${reference} \\
         --edited-bam ${edited_cram} \\
         --control-bam ${control_cram} \\
         --target-file ${target_file} \\
         ${unevaluable} \\
+        ${tagged_bam} \\
         --vcf-out ${meta.id}.offtarget_edits.vcf \\
         -o ${meta.id}.offtarget_analysis.tsv
 
@@ -37,5 +43,6 @@ process ECS_INDELS {
     """
 
     stub:
-    "touch ${meta.id}.offtarget_analysis.tsv ${meta.id}.offtarget_edits.vcf versions.yml"
+    def tagged_bam = params.offtarget_tagged_bam ? "${meta.id}.tagged.bam ${meta.id}.tagged.bam.bai" : ""
+    "touch ${meta.id}.offtarget_analysis.tsv ${meta.id}.offtarget_edits.vcf ${tagged_bam} versions.yml"
 }

modules/local/get_indels.nf

@@ -12,6 +12,7 @@ process GET_INDELS {
     output:
     tuple val(meta), path("${meta.id}.offtarget_analysis.tsv"), emit: indels_file
     tuple val(meta), path("${meta.id}.offtarget_edits.vcf"), emit: indels_vcf
+    tuple val(meta), path("${meta.id}.tagged.bam*"), optional: true, emit: tagged_bam
     //tuple val(meta), path("${meta.id}.ml_results.txt"), emit: ml_results
     //tuple val(meta), path("${meta.id}.fp_filtered.txt"), emit: fp_log
     path "versions.yml",    emit: versions
@@ -25,8 +26,12 @@ process GET_INDELS {
         hotspot_file ? "--target-file ${hotspot_file}" : ""
     ].join(' ').trim()
 
-    """    
-    find_edited_reads.py ${inputs} -u ${meta.id}.unevaluable_reads.txt --vcf-out ${meta.id}.offtarget_edits.vcf -o ${meta.id}.offtarget_analysis.tsv
+    // Optional review aid: a window-restricted BAM whose reads carry an XC tag naming
+    // the per-read call, for colouring the pileup in IGV. Off unless explicitly asked.
+    def tagged_bam = params.offtarget_tagged_bam ? "--tagged-bam-out ${meta.id}.tagged.bam" : ""
+
+    """
+    find_edited_reads.py ${inputs} ${tagged_bam} -u ${meta.id}.unevaluable_reads.txt --vcf-out ${meta.id}.offtarget_edits.vcf -o ${meta.id}.offtarget_analysis.tsv
 
     cat <<-END_VERSIONS > versions.yml
     ${task.process}:
@@ -35,9 +40,11 @@ process GET_INDELS {
     """
 
     stub:
+    def tagged_bam = params.offtarget_tagged_bam ? "touch ${meta.id}.tagged.bam ${meta.id}.tagged.bam.bai" : ""
     """
     touch ${meta.id}.offtarget_analysis.tsv
     touch ${meta.id}.offtarget_edits.vcf
+    ${tagged_bam}
 
     cat <<-END_VERSIONS > versions.yml
     ${task.process}:

nextflow.config

@@ -59,6 +59,7 @@ params {
     offtarget_top                 = 60      // top-N to print / snapshot
     offtarget_snapshots           = false   // render IGV-style pileup PNGs for LIKELY EDITs
     offtarget_ecs_unevaluable_log = false   // write per-read unevaluable-reads debug log (very large: ~0.1-1 TB/sample); off by default
+    offtarget_tagged_bam          = false   // emit a target-window BAM whose reads carry an XC tag with the per-read call, for IGV review; off by default
     offtarget_ecs_edit_threshold  = 0.0     // ECS indel_fraction strictly above this = an ECS edit (ecs_is_edit)
     offtarget_germline_max_ctrl_if = 0.05   // matched-normal (WGS) indel frac above this = germline/artifact, not a somatic edit (label 0)
     offtarget_hi_score            = 0.60    // WGS shape score >= this = detected (recall curve)

nextflow_schema.json

@@ -166,6 +166,11 @@
                     "help_text": "Off by default. This log reaches roughly 0.1-1 TB per sample and was the sole cause of work-directory ENOSPC failures; enable only when debugging read evaluation.",
                     "hidden": true
                 },
+                "offtarget_tagged_bam": {
+                    "type": "boolean",
+                    "description": "Emit a per-read tagged BAM from the ECS caller for review in IGV",
+                    "help_text": "Off by default. Writes <id>.tagged.bam(.bai) in which every read carries an XC string tag naming how the caller classified it (Edited_Deletion_5bp, Unedited_WT, Skipped_LowMapQ, ...), so a reviewer can load it in IGV and use Color alignments by -> tag -> XC. Output is restricted to the target windows (target +/- 150 bp), not the whole genome, but it still scales with target count and depth: one AAVS1 site14 ECS sample (1149 targets, ~11000x) produced a 800 MB BAM of 21.8 M reads and pushed ECS_INDELS from 1.4 GB to 6.3 GB peak RSS, so the process is given 24 GB instead of 8 GB while this is set. Check the size on one sample before enabling it across a cohort. See docs/OFFTARGET_WORKFLOW.md."
+                },
                 "offtarget_ecs_edit_threshold": {
                     "type": "number",
                     "default": 0.0,

tests/conftest.py

@@ -80,3 +80,174 @@ def workspace(tmp_path):
     ], columns=["rank", "sample", "chrom", "start", "alt", "dragen_af", "score", "verdict_pon"]
     ).to_csv(d / "worklist_pon.csv", index=False)
     return d
+
+
+# ---------------------------------------------------------------------------
+# Synthetic aligned-read workspace for find_edited_reads.py
+#
+# The glue fixtures above are pure tables; the read-level tagging path needs
+# real alignments. This builds the smallest workspace that reproduces the two
+# things that make tagging non-trivial:
+#
+#   * two targets 200 bp apart -- far enough that pyranges clusters them
+#     separately (slack = --target-window = 150), close enough that their
+#     +/-150 fetch windows overlap, so some reads are visited twice;
+#   * a read pair sitting in that overlap which one target calls Unedited_WT
+#     and the other cannot evaluate, i.e. the tag-conflict case.
+# ---------------------------------------------------------------------------
+
+ECS_CHROM = "chr1"
+ECS_CONTIG_LEN = 3000
+ECS_TARGET_A = 1001        # 1-based; the edited site
+ECS_TARGET_B = 1201        # 1-based; a quiet site
+ECS_DEL_LEN = 5
+ECS_N_EDIT = 12            # read pairs carrying the deletion
+ECS_N_WT = 8               # read pairs spanning target A cleanly
+ECS_N_DUP = 2              # duplicate-flagged pairs
+ECS_N_LOWMAPQ = 2          # MAPQ below the default floor of 20
+ECS_N_MISMATCH = 2         # NM above the default ceiling of 4
+ECS_N_OVERLAP = 6          # pairs inside the window overlap, spanning neither target
+ECS_N_SPAN_B = 3           # pairs in the overlap that DO span target B (conflict case)
+
+
+def _write_ecs_fasta(path):
+    import random
+    random.seed(7)
+    seq = "".join(random.choice("ACGT") for _ in range(ECS_CONTIG_LEN))
+    with open(path, "w") as fh:
+        fh.write(f">{ECS_CHROM}\n")
+        for i in range(0, ECS_CONTIG_LEN, 60):
+            fh.write(seq[i:i + 60] + "\n")
+    import pysam
+    pysam.faidx(str(path))
+    return pysam.FastaFile(str(path))
+
+
+def _ecs_pair(fasta, name, r1_start, r1_cigar, r2_start, mapq=60, dup=False,
+              r1_nm=0, r2_nm=0):
+    """One properly-paired FR pair (R1 forward, R2 reverse) with real sequence."""
+    import pysam
+
+    def qseq(start, cigar):
+        out, ref = [], start
+        for op, ln in cigar:
+            if op == 0:                       # M
+                out.append(fasta.fetch(ECS_CHROM, ref, ref + ln)); ref += ln
+            elif op == 2:                     # D
+                ref += ln
+            elif op == 1:                     # I
+                out.append("A" * ln)
+            elif op == 4:                     # S
+                out.append("T" * ln)
+        return "".join(out)
+
+    span = r2_start + 100 - r1_start
+    reads = []
+    for is_read1, start, cigar, mate_start, nm in (
+            (True, r1_start, r1_cigar, r2_start, r1_nm),
+            (False, r2_start, [(0, 100)], r1_start, r2_nm)):
+        a = pysam.AlignedSegment()
+        a.query_name = name
+        a.query_sequence = qseq(start, cigar)
+        a.reference_id = 0
+        a.reference_start = start
+        a.mapping_quality = mapq
+        a.cigar = cigar
+        a.next_reference_id = 0
+        a.next_reference_start = mate_start
+        a.template_length = span if is_read1 else -span
+        a.query_qualities = pysam.qualitystring_to_array("I" * len(a.query_sequence))
+        a.is_paired = True
+        a.is_proper_pair = True
+        a.is_read1 = is_read1
+        a.is_read2 = not is_read1
+        a.is_reverse = not is_read1
+        a.mate_is_reverse = is_read1
+        a.is_duplicate = dup
+        a.set_tag("NM", nm, value_type="i")
+        a.set_tag("MC", "100M", value_type="Z")
+        reads.append(a)
+    return reads
+
+
+def _build_ecs_cram(fasta, fasta_path, out_bam, reads):
+    import pysam
+    header = {"HD": {"VN": "1.6", "SO": "coordinate"},
+              "SQ": [{"SN": ECS_CHROM, "LN": ECS_CONTIG_LEN}]}
+    reads.sort(key=lambda r: r.reference_start)
+    tmp = str(out_bam) + ".tmp.bam"
+    with pysam.AlignmentFile(tmp, "wb", header=header) as out:
+        for r in reads:
+            out.write(r)
+    pysam.sort("-o", str(out_bam), tmp)
+    Path(tmp).unlink()
+    pysam.index(str(out_bam))
+    # production opens its inputs as CRAM, so hand the caller a CRAM
+    cram = str(out_bam)[:-4] + ".cram"
+    with pysam.AlignmentFile(str(out_bam)) as src, \
+            pysam.AlignmentFile(cram, "wc", template=src,
+                                reference_filename=str(fasta_path)) as out:
+        for r in src:
+            out.write(r)
+    pysam.index(cram)
+    return cram
+
+
+@pytest.fixture
+def ecs_reads_workspace(tmp_path):
+    """Synthetic reference + targets VCF + edited/control CRAMs for the ECS caller.
+
+    Returns a dict of paths plus the read counts the tags should reconcile with.
+    """
+    pysam = pytest.importorskip("pysam")
+    d = tmp_path
+    fasta_path = d / "ref.fa"
+    fasta = _write_ecs_fasta(fasta_path)
+
+    header = pysam.VariantHeader()
+    header.contigs.add(ECS_CHROM, length=ECS_CONTIG_LEN)
+    targets = d / "targets.vcf"
+    with pysam.VariantFile(str(targets), "w", header=header) as vout:
+        for pos in (ECS_TARGET_A, ECS_TARGET_B):
+            rec = vout.new_record()
+            rec.chrom = ECS_CHROM
+            rec.pos = pos
+            rec.id = "."
+            rec.ref = fasta.fetch(ECS_CHROM, pos - 1, pos)
+            rec.alts = ("N",)
+            rec.filter.add("PASS")
+            vout.write(rec)
+
+    del_cigar = [(0, 60), (2, ECS_DEL_LEN), (0, 40)]
+    edited = []
+    for i in range(ECS_N_EDIT):
+        edited += _ecs_pair(fasta, f"edit{i}", 940, del_cigar, 1180, r1_nm=ECS_DEL_LEN)
+    for i in range(ECS_N_WT):
+        edited += _ecs_pair(fasta, f"wt{i}", 940, [(0, 100)], 1180)
+    for i in range(ECS_N_DUP):
+        edited += _ecs_pair(fasta, f"dup{i}", 940, [(0, 100)], 1180, dup=True)
+    for i in range(ECS_N_LOWMAPQ):
+        edited += _ecs_pair(fasta, f"lowq{i}", 940, [(0, 100)], 1180, mapq=3)
+    for i in range(ECS_N_MISMATCH):
+        edited += _ecs_pair(fasta, f"mm{i}", 940, [(0, 100)], 1180, r1_nm=6, r2_nm=6)
+    for i in range(ECS_N_OVERLAP):
+        edited += _ecs_pair(fasta, f"both{i}", 1080, [(0, 100)], 1220)
+    # R1 lands in both padded windows: target A cannot evaluate it, target B
+    # calls it reference. Unedited_WT must win, and it must be written once.
+    for i in range(ECS_N_SPAN_B):
+        edited += _ecs_pair(fasta, f"spanb{i}", 1150, [(0, 100)], 1220)
+
+    control = []
+    for i in range(ECS_N_EDIT + ECS_N_WT):
+        control += _ecs_pair(fasta, f"ctl{i}", 940, [(0, 100)], 1180)
+
+    return {
+        "dir": d,
+        "fasta": fasta_path,
+        "targets": targets,
+        "edited": _build_ecs_cram(fasta, fasta_path, d / "edited.bam", edited),
+        "control": _build_ecs_cram(fasta, fasta_path, d / "control.bam", control),
+        "n_edit_reads": ECS_N_EDIT,
+        "del_len": ECS_DEL_LEN,
+        "total_records": len(edited),
+    }