Generate DiffATAC counts matrices/results from both filtered.bam and dedup.bam #138

Background

filtered.bam (quality-filtered) is produced before Picard MarkDuplicates runs, so it retains PCR duplicates; dedup.bam is derived from filtered.bam by marking + removing duplicates. Currently, both the Tn5 nicking-sites and read-counts matrices used for DiffATAC (DESeq2/edgeR) are generated from tn5sites_bam/reads_bam, which are built from filtered.bam (see create_tn5bams in workflow/rules/align.smk) — i.e. duplicates are not removed before counting today.

Standard ATAC-seq practice (ENCODE ATAC pipeline, ArchR, Signac, DiffBind/csaw) counts from deduplicated fragments, since PCR duplication rate varies by sample (input amount, PCR cycles, library complexity) and isn't corrected by DESeq2/edgeR size factors, which can confound differential accessibility calls.

Counter-consideration: Tn5 insertion is biased toward accessible chromatin, so independent molecules can genuinely converge on the same insertion coordinate at highly accessible sites — Picard's paired-end duplicate definition (requires both fragment ends to match) partially mitigates this, and ROI windows (fixed_width: 500 by default) further dilute the effect, but it's not eliminated.

Proposal

Rather than picking one BAM source, generate both counts matrices and DiffATAC results — one from filtered.bam, one from dedup.bam — in separate output folders, so users can compare/choose.

Plan

Add a BAM_TYPES = ["dedup", "filtered"] axis and thread it through existing rules as a wildcard (not duplicated rules):

1. align.smk — create_tn5bams: wildcard on {bamtype} so it runs once per BAM type, writing to visualization/{bamtype}/tn5sites_bam/... and visualization/{bamtype}/reads_bam/....
2. create_bigwigs: same wildcard threading (open question: needed for both types, or just counts/DESeq2?).
3. diffatac.smk — get_counts_table / scale_counts_table / diffatac: add {bamtype} wildcard so outputs land in peaks/{peakcaller}/fixed_width/counts/{bamtype}/... and peaks/{peakcaller}/DiffATAC/{bamtype}/{method}/....
4. init.smk: define BAM_TYPES and extend get_diffatac_input()/final target list to expand over it.
5. Docs: update docs/outputs.md and docs/extra.md to explain both folders and give guidance (recommend dedup for standard differential testing; filtered retained for comparison, with the Tn5-convergence caveat noted).
6. CHANGELOG: user-facing entry for the new dedup-based counts/DiffATAC outputs alongside the existing filtered-based ones.

Open questions

• Should dedup be the default recommended path (highlighted in docs) while filtered is kept for backward compatibility, or present both as equal options?
• Duplicate bigwigs per bamtype, or just counts matrices/DiffATAC results (to limit extra compute/storage)?
• Naming: dedup/filtered, or something like dedup/nondedup for folder-name clarity?

_{⚡ Generated using AI ⚡}