Fix create_tn5bams failure from negative BED coordinates in Tn5 conversion

[kopardev]

Summary

create_tn5bams fails intermittently/repeatedly for some replicates when building Tn5 BAM outputs, due to malformed BED coordinates generated during Tn5 site conversion.

Observed Failure

In a recent run, Snakemake repeatedly retried create_tn5bams and failed at the first bedToBam step with:

• Error: malformed BED entry ... Start Coordinate detected that is < 0. Exiting.

Example run context:

• Workdir: /data/Boufraqech_group/analysis/.temp/aspen_run_for_Ying
• Rule: create_tn5bams
• Failing command stage: bedToBam -i ...tn5sites.bed -g ...hs1_chrR.genome

Root Cause

Current logic in workflow/scripts/ccbr_atac_bam2tn5bed.py:

1. Applies a reverse-read shift of -5 without lower-bound guarding.
2. Derives start/end from pair-level min/max boundaries, then applies shifts tied to read orientation in a way that can produce incorrect or out-of-range cut positions.
3. Writes BED rows without enforcing valid BED interval constraints.

This can emit entries with negative starts (especially for alignments near chromosome start), which causes bedToBam to abort.

Why This Matters

• Blocks pipeline progress for affected samples.
• Causes repeated cluster retries and wasted compute.
• Risks downstream inconsistency in tn5/reads counting and visualization paths.

Proposed Fix

Update workflow/scripts/ccbr_atac_bam2tn5bed.py to:

1. Use strand-correct per-read cut-site derivation for ATAC conventions (+4/-5), aligned to ENCODE/deepTools semantics.
2. Guard and sanitize coordinates before BED write:
   • enforce start >= 0
   • enforce end > start
   • clamp to valid bounds as needed
3. Apply the same validation/sanitization policy to both outputs:
   • tn5sites.bed
   • reads.bed
4. Emit counters/warnings for clipped/skipped entries to aid QC.

Validation Plan

1. Re-run failing replicate(s) and confirm no malformed BED / negative-coordinate bedToBam errors.
2. Assert generated BED intermediates have no invalid intervals.
3. Confirm both BAM outputs and indexes are produced:
   • tn5sites.bam(.bai)
   • reads.bam(.bai)
4. Smoke-test downstream counting and bigwig generation for both methods (tn5sites, reads).

Acceptance Criteria

• create_tn5bams completes successfully on previously failing replicates.
• No malformed BED coordinate errors in logs.
• Downstream outputs remain complete and consistent.

References

• Script: workflow/scripts/ccbr_atac_bam2tn5bed.py
• Rule: workflow/rules/align.smk (create_tn5bams)

_{⚡ Generated using AI ⚡}

[kopardev]

Fixed on the issue_117 branch in commit 4841431 (fix(align): guard Tn5/read BED coordinates against negative/out-of-range values):

• ccbr_atac_bam2tn5bed.py: added parse_genome_sizes(), clamp_interval(), and tn5_cutsite() helpers. The script now requires --genomefile/-g, clamps every Tn5 and read interval into [0, chrom_len], skips unmapped/secondary/supplementary pairs and unknown chromosomes, writes per-read (strand-aware) Tn5 cutsite entries, and reports pairs_seen/pairs_written/pairs_skipped/unknown_chrom/tn5_clipped/tn5_skipped/reads_clipped/reads_skipped counters to stderr for QC.
• align.smk: passes --genomefile {params.genomefile} (already defined for the existing bedToBam calls) through to the script.

Validated via a smoke test (0 clipped/skipped bad intervals) and confirmed via dry-run that --genomefile is wired through correctly for all replicates. This also unblocked a real run that was previously failing at bedToBam with Start Coordinate detected that is < 0 for hs1_chrR.

_{⚡ Generated using AI ⚡}