Finalize peak calling

[adamklie]

Context

For igvf-data/igvf_sc-islet_10X-Multiome, finalize peak calling approaches.

Tasks

• Review v2 code (bin/5_peak_calling) and results (results/5_peak_calling)
• Call peaks without replicates (for sequence model training)
• Implement peak annotation and post-processing in bin/5_peak_calling for all peak sets (pseudobulk, single-cell, with/without replicates). Annotation for consensus peaks only, post-processing for all
• Comparative analysis: new results vs each other, then new vs old (scratch/2026_01_30/5_peak_analysis)

[adamklie]

Dependencies

• Feeds into Retrain ChromBPNet models at different resolutions #1 (ChromBPNet retraining), Recall peak-to-gene links using Gaulton lab pipeline #2 (peak-to-gene links), Upload v2 data to GCP #6 (upload v2 to GCP)

[adamklie]

Implementation Plan

Step 1: Call peaks without replicates (for ChromBPNet training)

Same cell_type-condition grouping but without replicate reproducibility filter → more peaks for model training.

Status: SLURM job 7792372 submitted
Output: results/5_peak_calling/cell_type-condition_no_replicates/

Files modified:

• tools/single_cell_utilities/snapatac2/call_peaks.sh — conditional replicate_key handling
• bin/5_peak_calling/1_call_peaks.sh — added no-replicate run command

Step 2: Annotate consensus peaks with HOMER

Run annotatePeaks.pl on consensus peak BED files for genomic feature annotation.

Status: SLURM job 7792373 submitted (2 consensus files)
Output: results/5_peak_calling/annotation/

Files modified:

• bin/5_peak_calling/2_annotate_peaks.sh — dynamic annotation_num from input TSV, updated script paths

Step 3: Post-process / filter peaks

Three operations per narrowPeak file:

1. Summit deduplication — cluster overlapping peaks, keep best summit per cluster
2. Blacklist removal — remove ENCODE blacklist regions
3. Top-N filtering — keep top 250k peaks if > 250k remain

Status: SLURM job 7792375 submitted (124 array tasks for existing narrowPeak files)
Output: results/5_peak_calling/postprocessed/

Files modified/created:

• tools/single_cell_utilities/peaks/granges.py — added deduplicate_summits() function
• tools/single_cell_utilities/peaks/filter_peaks.py — integrated summit dedup via --deduplicate_summits flag
• tools/single_cell_utilities/peaks/filter_peaks.sh — new SLURM wrapper for filter_peaks.py
• bin/5_peak_calling/3_postprocess_peaks.sh — implemented (was empty)

Step 4: Comparative analysis (scratch)

Sanity-check notebook comparing peak sets across strategies and versions.

Status: Pending (depends on steps 1-3 completing)
Location: scratch/2026_03_23/5_peak_calling/comparative_analysis.ipynb

Execution order

Step 1 (call peaks no-rep, job 7792372) ──┐
                                           ├─→ Step 3 (filter all, job 7792375 + no-rep) → Step 4 (compare)
Step 2 (annotate, job 7792373) ───────────┘

Path fixes

Updated stale /cellar/users/aklie/projects/igvf/single_cell_utilities/ references to /carter/users/aklie/projects/stimulated_sc-islets/tools/single_cell_utilities/ in:

• call_peaks.sh, annotate_peaks.sh, annotate_peaks.py, filter_peaks.py

[adamklie]

Notes from project reorganization (2026-03-23)

Filed in igvf-data/igvf_sc-islet_10X-Multiome/README.md under Outstanding TODOs:

Plan

• Call peaks within celltype-conditions using pseudoreplicate approach
• Call peaks using condition-timepoint with true replicates
• Compare reproducible vs non-reproducible in scratch

Pseudobulk ATAC status

• 4_cell_annotation/7_pseudobulk_ATAC cell_type pseudobulk almost done (tagAlign creation failed with sort OOM, needs rerun)
• cell_type-condition pseudobulk not started

Current peak calling setup (from README)

• MACS3 v3.0.0 via SnapATAC2 v2.5.3
• Fragments converted to insertion sites; q=0.05; shift=-100; extsize=200; min_len=200
• Consensus peaks via SnapATAC2 v2.5.3

[adamklie]

Status Update (2026-03-23)

Completed

• Peak calling without replicates — SLURM job 7792372 running (~2.5h in, expected 3-5h)
• HOMER annotation — 2 consensus peak sets annotated (job 7792499 completed)
• Post-processing implemented — summit deduplication + blacklist + top-250k filtering
  • 124 narrowPeak files processed (job 7792375, 0 failures)
  • Example: SC.beta_control 313k → 204k peaks after dedup+blacklist
• ATAC pseudobulk generation — notebook created to aggregate peak matrices into sample-level counts for DESeq2, runs for all peak sets (job 7797645)
• Comparative analysis notebook — created at scratch/2026_03_23/5_peak_calling/comparative_analysis.ipynb

Files created/modified

• tools/single_cell_utilities/snapatac2/call_peaks.sh — conditional replicate_key
• tools/single_cell_utilities/peaks/granges.py — deduplicate_summits() function
• tools/single_cell_utilities/peaks/filter_peaks.py — summit dedup flag
• tools/single_cell_utilities/peaks/filter_peaks.sh — new SLURM wrapper
• tools/single_cell_utilities/peaks/annotate_peaks.py — fixed stale import
• bin/5_peak_calling/1_call_peaks.sh — no-rep run command
• bin/5_peak_calling/2_annotate_peaks.sh — dynamic annotation_num
• bin/5_peak_calling/3_postprocess_peaks.sh — implemented (was empty)
• bin/5_peak_calling/4_pseudobulk_peaks.ipynb — new, loops all peak sets

Pending

• Peak calling no-rep job to complete → then run annotation + post-processing on those peaks
• Re-run pseudobulk notebook once no-rep peaks are available
• Run comparative analysis notebook

[adamklie]

Status Update (2026-03-24 11:30 PDT)

No-replicate peaks — FULLY PROCESSED

• Peak calling completed (job 7792372, 3h33m) — 24 narrowPeak files, 521,395 consensus peaks
• HOMER annotation completed (job 7818997) — consensus peaks annotated
• Post-processing completed (job 7818998, 24/24 tasks) — summit dedup + blacklist + 250k cap

All 3 peak sets now fully processed:

Peak set	Consensus peaks	narrowPeak files	Status
cell_type-condition_pseudoreplicates	368,130	24	Done
cell_type-condition-timepoint	266,606	100	Done
cell_type-condition_no_replicates	521,395	24	Done

Pending

• ATAC pseudobulk job 7819044 (generates sample-level counts for DESeq2)
• Comparative analysis notebook
• ChromBPNet training can now begin (Retrain ChromBPNet models at different resolutions #1)