Optional UMI-aware preprocessing for UMI-ATAC-seq (`--with_umi`)

Description of feature

Summary

We'd like to propose adding an optional --with_umi path to nf-core/atacseq that implements UMI-aware preprocessing and deduplication for the UMI-ATAC-seq protocol (Zhu et al., Communications Biology, 2020: ATAC-seq with unique molecular identifiers improves quantification and footprinting). This issue outlines the motivation, the proposed design, and the implementation plan, and asks a few architectural questions before we start work. We raised this informally in Slack first (https://nfcore.slack.com/archives/CE5EL6326/p1786349372219949) and are opening this issue to continue the discussion and, if there's interest, align on design before submitting a PR.

Motivation

Standard ATAC-seq deduplication removes PCR duplicates by mapping coordinate. In highly accessible chromatin, two independent Tn5 insertions can land on the exact same coordinate ("natural duplicates") and get wrongly discarded. This problem grows with sequencing depth and disproportionately affects transcription-factor footprinting, since footprints concentrate in the most accessible regions.

UMI-ATAC-seq addresses this by carrying a UMI (≥6 bp) on the Tn5 adapter, between the TruSeq adapter and the 19 bp mosaic end (ME), so Read 1 has the structure [UMI][ME][genomic]. Reads sharing a coordinate but carrying different UMIs are independent insertions and are retained, rather than being collapsed as duplicates.

On the authors' rice panicle dataset (NCBI BioProject PRJNA602139), UMI-aware dedup (UD) vs. coordinate dedup (CD):

Finding	Headline result
Quantification concordance with the undeduplicated truth (ND)	R²(UD,ND) = 0.995 vs. R²(CD,ND) = 0.949
Duplication rate vs. chromatin accessibility	CD duplication rate rises with accessibility; UMI duplication rate stays flat
Reads rescued	~20% of coordinate-flagged duplicates, ~6% of total mapped reads, concentrated in the most accessible bins
Footprint sensitivity	UD recovers 5,718 footprints missed by CD, vs. only 76 lost (~+50%)
Biological validity of rescued footprints	UD-only footprints enriched for FIMO TF motif overlap vs. shuffled controls (P = 4.9×10⁻³⁴)

Why not just reuse the authors' code

We reviewed the reference implementation accompanying the paper. It's a proof-of-concept (three short scripts glued by README shell commands, 7 stars, no releases, no tests, no license) with two specific limitations we don't want to inherit:

• Exact-string UMI matching in dedup (rm_umi_dup.py) — a single sequencing error in the 6 bp UMI creates a false-unique read and inflates counts. umi_tools dedup (directional-adjacency method) handles this correctly and should be strictly better.
• In-memory ME-trimming script (remove_me.py) with a hardcoded n_jobs=10 — doesn't scale and isn't suited to a production Nextflow pipeline.

So our intent is to reuse the concept from the paper, not the code, and reimplement it natively on top of existing, tested nf-core modules.

Proposed design

Add a --with_umi flag, off by default, that only affects the read-processing/deduplication layer. Everything downstream (alignment, filtering, Tn5 shifting, peak calling, QC, MultiQC) stays on the current code path.

Step	Paper's tool	Proposed nf-core module	Notes
UMI extraction	umi_tools extract	umitools/extract	Same convention as nf-core/rnaseq's --umitools_bc_pattern; extracts the UMI from R1 and appends it to the read name
Mosaic End (ME) trimming + re-pairing	remove_me.py + bbmap repair.sh	cutadapt, anchored 5′ adapter (-g) on the ME sequence	Single step; keeps R1/R2 in sync natively, replacing two bespoke scripts
Alignment	BWA-mem	existing BWA/Bowtie2	Unchanged
Filtering (MAPQ, mito/chloroplast)	SAMtools	existing atacseq filters	Unchanged, including the +4/−5 Tn5 shift
Deduplication	rm_umi_dup.py (exact match)	umitools/dedup	Only when --with_umi is set; replaces Picard MarkDuplicates for this path; adds UMI error-correction
Peak calling / QC	MACS2	existing atacseq	Unchanged

Expected execution graph when --with_umi is enabled:

FASTQ
  │
  ▼
UMI extraction (umitools/extract)
  │
  ▼
Mosaic End trimming (cutadapt)
  │
  ▼
Alignment (existing)
  │
  ▼
Filtering (existing)
  │
  ▼
UMI deduplication (umitools/dedup)
  │  
  ▼
Tn5 shift → Peak calling → QC + MultiQC

New parameters:

• --with_umi — enables the UMI-aware path (default: false, standard path unchanged)
• --umi_pattern — UMI barcode pattern for umitools/extract (default matches the paper's 6 bp UMI, e.g. NNNNNN, overridable)
• --me_sequence — Mosaic End sequence for the cutadapt trim (default AGATGTGTATAAGAGACAG, the standard Tn5/Nextera ME, overridable since it's shared across other Tn5-based assays)

We'd implement the new steps as a single local subworkflow (e.g. umi_atacseq_preprocess) that wraps umitools/extract → cutadapt ME-trim and hands off a channel compatible with the existing alignment subworkflow, so the main workflow only needs to branch on --with_umi to call it. Dedup routing (umitools/dedup vs. the existing coordinate-based step) would be a conditional in the existing dedup stage rather than a separate subworkflow.

Validation plan

We'd validate against the paper's own dataset before proposing this as production-ready:

• Primary dataset: rice young panicle UMI-ATAC-seq, PRJNA602139 (7 samples, the paper's own data)
• Negative/sanity controls: Arabidopsis SRX2000808 and human HEK293 SRX3511086 (standard ATAC, no UMI) to confirm the accessibility-vs-duplication trend still reproduces without UMIs

For each rice sample we'd process three ways from the same alignment — no dedup (ND), coordinate dedup (CD), UMI dedup (UD) — and check:

Metric	Target
Reads rescued by UMI dedup	Same order of magnitude as paper (~20% of CD-flagged duplicates, ~6% of total mapped reads)
Duplication rate vs. accessibility decile	CD rate rises with accessibility; UMI dup rate stays flat
Quantification concordance	R²(UD,ND) > R²(UD,CD)
Peak count / fold-enrichment (CD vs UD)	Small difference in count, improved fold-enrichment on top peaks
Footprints identified (UD vs CD)	UD yields substantially more (paper: +50%+)
UD-only footprint ↔ FIMO motif overlap	Significant enrichment over shuffled controls
Standard ATAC QC	TSS enrichment, FRiP, insert-size periodicity within expected ranges

Where our numbers diverge from the paper, we'd expect the cause to be umitools/dedup's error-correction vs. the authors' exact-match dedup, and would document that explicitly (with a direct comparison) rather than treating it as a discrepancy to hide.

Questions for the community

1. Would support for this protocol be considered in scope for nf-core/atacseq?
2. If so, would a local subworkflow specific to this pipeline be preferred, or is there interest in new reusable nf-core/modules components? (To our knowledge the UMI-ATAC-seq protocol as published is ATAC-seq-specific, though the ME-trim step generalizes to other Tn5-based assays.)
3. For protocol-specific steps like Mosaic End trimming, is it generally preferable to keep them local to the pipeline, or is there interest in modularizing despite limited reuse?
4. Are there existing architectural patterns in other nf-core pipelines (e.g. nf-core/rnaseq's --with_umi) that we should follow for this kind of optional, protocol-specific preprocessing branch?
5. Is there existing test data we should use, or should we prepare a small subsampled UMI-ATAC-seq test dataset for CI?

We're happy to do the implementation and validation work (on a fork/feature branch, reproducing the paper's footprinting-gain figures on PRJNA602139) and contribute it back as a PR if there's appetite for this.

References

• Zhu, T., Liao, K., Zhou, R. et al. ATAC-seq with unique molecular identifiers improves quantification and footprinting. Commun Biol 3, 675 (2020). https://doi.org/10.1038/s42003-020-01403-4
• Reference implementation (proof-of-concept, not proposed for direct reuse): https://github.com/tzhu-bio/UMI-ATAC-seq
• Data: NCBI BioProject PRJNA602139 (rice), SRX2000808 (Arabidopsis), SRX3511086 (human HEK293)

Start by reviewing the existing alignment, filtering, deduplication, peak-calling, and QC stages, then compare them with the proposed local umi_atacseq_preprocess subworkflow and --with_umi branch. Done means the optional extraction, Mosaic End trimming, and UMI deduplication path integrates without changing the standard path and is validated against PRJNA602139 plus the listed negative controls using the proposed QC metrics.