diff --git a/.gitignore b/.gitignore
index eac94acc..61ab5ed7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -21,6 +21,7 @@ tmp
 x.sh
 WIP
 dryrun.log
+reset_rose_test_outputs.sh
 POOL_CONTROLS_FEATURE.md
 .snakemake
 .vscode
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1028f860..183cad30 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,25 @@
+## development version
+
+### New Features
+
+- **Treatment-level merged peak workflow + cross-treatment HOMER/ROSE aggregation**: Added consensus peak merging across treatment replicates (no peak re-calling) and parallel treatment-level annotation workflows. New outputs include `treatment_merged` peak sets for each caller/mode, treatment-level HOMER/AME outputs (`annotation/homer_treatment`), treatment-level ROSE outputs (`annotation/rose_treatment`), and global aggregate files across all treatments for HOMER annotation summaries and ROSE enhancer-gene/gene-enhancer mappings. Existing replicate-level outputs remain unchanged.
+- **ROSE gene mapping outputs for super-enhancers**: Integrated Younglab `ROSE_geneMapper.py` into the ROSE rule so super-enhancer runs now emit `rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt` and `rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt` in each ROSE output folder. Added output-contract coverage for low-peak fallback by writing placeholders when ROSE is skipped. (#248)
+- **Slurm run-state markers and CLI run commentary**: `runmode=run` now manages explicit workdir marker files (`pipeline.running`, `pipeline.completed`, `pipeline.failed`, `pipeline.canceled`) plus `pipeline.status.json` sidecar metadata (state, timestamp, job id, reason). Starting a new slurm run in the same workdir rotates old markers and writes a fresh `pipeline.running`. Submission/runtime transitions are handled directly in launcher/submit script logic, and CLI output now includes a concise run summary + monitor hints. (#250)
+- **Enriched `pipeline.status.json` metadata**: Extended the run-state sidecar with audit and timing fields: `user`, `submission_timestamp_utc`, `start_timestamp_utc`, `duration_seconds`, `exit_code`, `tasks_done`, `tasks_total`, and `snakemake_log`. Duration and task counts are populated by the SLURM job itself; submission timestamp is persisted via a `pipeline.submit_epoch` helper file so it survives the headnode-to-compute-node handoff. Stale sbatch scripts from earlier releases are automatically detected and regenerated. (#250, @kopardev)
+- **Live progress summary in `pipeline.running` during `run`**: The generated SLURM submit script now runs a background monitor that periodically parses Snakemake progress lines and writes a human-readable status block to `pipeline.running` (completed steps, percent complete, remaining steps, timestamp). It includes delayed-log handling when `snakemake.log` is not immediately available after submission, plus clean monitor shutdown when Snakemake exits. (#254, @kopardev)
+- **Advanced scripts sync mode in wrapper**: Added `-m=syncscripts` as an advanced wrapper mode to re-sync `workflow/scripts` into an initialized work directory without running `init` or `reset`. The legacy `rescript` mode remains available as a backward-compatible alias, and script sync now excludes `__pycache__/` directories. (#251, @kopardev)
+
+### Documentation Improvements
+
+- **User guide: treatment-level merged outputs documented**: Added `treatment_merged` peak outputs and new treatment-level annotation directories (`homer_treatment`, `rose_treatment`) to `output.md`, including global aggregate output paths under `peaks/annotation/`.
+- **User guide: `pipeline.status.json` field reference**: Added a table documenting all 18 fields in `pipeline.status.json`, their availability (headnode vs SLURM job), and what each represents. (#250, @kopardev)
+
+### Bug Fixes
+
+- **HOMER annotation: removed unsupported `annotatePeaks.pl -preparsedDir` option**: `annotatePeaks.pl` does not accept `-preparsedDir`, which caused HOMER annotation failures in replicate-level, DEG-level, and treatment-level motif rules. The option is now removed from `annotatePeaks.pl` calls while retained on `findMotifsGenome.pl`, where writable scratch-backed preparsed caching is still supported. (#255, @kopardev)
+- **Treatment-level motif Slurm resources**: Added explicit `homer_motif_treatment` and `ame_motif_enrichment_treatment` entries to the Biowulf cluster configuration so treatment-level HOMER/AME jobs use the intended 120GB, 32-thread, 2-day resource requests instead of falling back to the 40GB, 2-hour defaults. (#252)
+- **ROSE geneMapper schema robustness**: Some ROSE runs emitted `rose_input_SuperEnhancers.table.txt` without trailing `enhancerRank`/`isSuper` columns, causing `ROSE_geneMapper.py` to fail when parsing rank values. The `rose` rule now detects this schema variant and normalizes the table before gene mapping so jobs do not hard-fail on this format difference. (#248)
+
 ## CARLISLE 2.8.0
 
 ### New Features
diff --git a/VERSION b/VERSION
index 834f2629..9781df78 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-2.8.0
+2.8.0-dev
diff --git a/carlisle b/carlisle
index e2b70920..96626c26 100755
--- a/carlisle
+++ b/carlisle
@@ -25,7 +25,7 @@ SCRIPTBASENAME=$(readlink -f $(basename $0))
 #define cluster, partitions dependent on host
 hostID=`echo $HOSTNAME`
 PARTITIONS="norm"
-cluster_specific_yaml="cluster_biowulf.yaml"
+cluster_specific_yaml="${CARLISLE_CLUSTER_CONFIG:-cluster_biowulf.yaml}"
 tools_specific_yaml="tools_biowulf.yaml"
 
 # essential files
@@ -36,7 +36,6 @@ ESSENTIAL_FOLDERS="workflow/scripts"
 
 # ## setting PIPELINE_HOME
 export PIPELINE_HOME=$(readlink -f $(dirname "$0"))
-echo "Pipeline Dir: $PIPELINE_HOME"
 # set snakefile
 SNAKEFILE="${PIPELINE_HOME}/workflow/Snakefile"
 
@@ -51,41 +50,78 @@ function get_version (){
 echo "Version: $PIPELINE_VERSION"
 }
 
+function log_info() {
+  echo "INFO  $*"
+}
+
+function log_step() {
+  echo "STEP  $*"
+}
+
+function log_ok() {
+  echo "OK    $*"
+}
+
+function log_warn() {
+  echo "WARN  $*"
+}
+
+function log_error() {
+  echo "ERROR $*"
+}
+
+function log_next() {
+  echo "NEXT  $*"
+}
+
+function log_divider() {
+  echo "------------------------------------------------------------------"
+}
+
 function usage() {
 # This function prints generic usage of the wrapper script.
   echo "${SCRIPTBASENAME}
   --> run CARLISLE
-  Cut And Run anaLysIS pipeLinE
-
-  USAGE:
-    bash ${SCRIPTNAME} -m/--runmode=<RUNMODE> -w/--workdir=<WORKDIR>
-  Required Arguments:
-  1.  RUNMODE: [Type: String] Valid options:
-      *) init : initialize workdir
-      *) run : run with slurm
-      *) reset : DELETE workdir dir and re-init it
-      *) dryrun : dry run snakemake to generate DAG
-      *) unlock : unlock workdir if locked by snakemake
-      *) runlocal : run without submitting to sbatch
-      *) runtest: run on cluster with included hg38 test dataset
-  2.  WORKDIR: [Type: String]: Absolute or relative path to the output folder with write permissions.
-
-  Optional Arguments:
-      *) -f / --force: Force flag will re-initialize a previously initialized workdir
+      Cut And Run anaLysIS pipeLinE
+      Version: ${PIPELINE_VERSION}
+
+  USAGE
+    bash ${SCRIPTNAME} -w/--workdir=<WORKDIR> -m/--runmode=<RUNMODE>
+
+  REQUIRED ARGUMENTS
+    1) RUNMODE [String]
+       Standard modes:
+         - init        initialize workdir
+         - run         run with slurm
+         - dryrun      dry run snakemake to generate DAG
+         - unlock      unlock workdir if locked by snakemake
+         - runlocal    run without submitting to sbatch
+         - runtest     run on cluster with included hg38 test dataset
+
+       Advanced modes:
+         - reset       DELETE workdir and re-initialize it
+         - syncscripts sync workflow/scripts into workdir/scripts
+         - synccluster sync resources/<cluster_config.yaml> into workdir/config/cluster.yaml
+
+    2) WORKDIR [String]
+       Absolute or relative path to the output folder with write permissions.
+
+  ADVANCED ARGUMENTS
+    -f, --force              re-initialize an existing workdir
+    -c, --singcache=<DIR>    set singularity/apptainer cache directory (Biowulf: do not use; CCBR manages cache dir)
+    --cluster-config=<FILE>  cluster yaml filename in resources/ (default: cluster_biowulf.yaml)
+
+  OTHER ARGUMENTS
+    -v, --version            print CARLISLE version
+    -h, --help               show this message
   "
 }
 
 function err() {
 # This function print a error message (argument 1), the usage and then exits with non-zero status
-  cat <<< "
-  #
-  #
-  #
-    $@
-  #
-  #
-  #
-  " && usage && exit 1 1>&2;
+  log_error "$@"
+  usage
+  exit 1 1>&2
 }
 
 function init() {
@@ -97,13 +133,15 @@ function init() {
   if [[ -d $WORKDIR ]] && [[ -z $FORCEFLAG ]]; then
     err "Folder $WORKDIR already exists! If you'd like to re-initialize use the -f/--force flag $FORCEFLAG|"
   fi
+  log_step "[init] Preparing workdir"
+  log_info "Workdir: $WORKDIR"
   mkdir -p $WORKDIR
   export WORKDIR
 
   # copy essential files
   if [[ ! -d $WORKDIR/config ]]; then mkdir $WORKDIR/config; fi
   for f in $ESSENTIAL_FILES; do
-    echo "Copying essential file: $f"
+    log_info "Copying essential file: $f"
     fbn=$(basename $f)
     cat ${PIPELINE_HOME}/$f |\
       envsubst '$PIPELINE_HOME $WORKDIR $PIPELINE_VERSION' \
@@ -111,19 +149,24 @@ function init() {
   done
 
   # rename config dependent on partition used
+  if [[ ! -f ${PIPELINE_HOME}/resources/$cluster_specific_yaml ]]; then
+    err "Cluster config ${PIPELINE_HOME}/resources/$cluster_specific_yaml not found. Set --cluster-config=<FILE> or CARLISLE_CLUSTER_CONFIG."
+  fi
   cp ${PIPELINE_HOME}/resources/$cluster_specific_yaml $WORKDIR/config/cluster.yaml
   cp ${PIPELINE_HOME}/resources/$tools_specific_yaml $WORKDIR/config/tools.yaml
 
   # copy essential folders
   for f in $ESSENTIAL_FOLDERS;do
+    log_step "[init] Syncing folder: $f"
     rsync -avz --no-perms --no-owner --no-group --progress $PIPELINE_HOME/$f $WORKDIR/
   done
 
   #create log and stats folders
-  if [ ! -d $WORKDIR/logs ]; then mkdir -p $WORKDIR/logs;echo "Logs Dir: $WORKDIR/logs";fi
-  if [ ! -d $WORKDIR/stats ];then mkdir -p $WORKDIR/stats;echo "Stats Dir: $WORKDIR/stats";fi
+  if [ ! -d $WORKDIR/logs ]; then mkdir -p $WORKDIR/logs;log_ok "Logs Dir: $WORKDIR/logs";fi
+  if [ ! -d $WORKDIR/stats ];then mkdir -p $WORKDIR/stats;log_ok "Stats Dir: $WORKDIR/stats";fi
 
-  echo "Done Initializing $WORKDIR. You can now edit $WORKDIR/config/config.yaml and $WORKDIR/config/samples.tsv"
+  log_ok "Initialization completed for $WORKDIR"
+  log_next "Edit $WORKDIR/config/config.yaml, $WORKDIR/config/samples.tsv, and $WORKDIR/config/contrasts.tsv"
 
 }
 
@@ -151,9 +194,28 @@ function set_singularity_binds(){
 
 function rescript(){
 # recopy the scripts folder to WORKDIR
+  log_step "[syncscripts] Re-syncing workflow/scripts into workdir"
   check_essential_files
-  rsync -avz --no-perms --no-owner --no-group --progress ${PIPELINE_HOME}/workflow/scripts/ $WORKDIR/scripts/
-  echo "$WORKDIR/scripts folder has been updated!"
+  rsync -avz --no-perms --no-owner --no-group --progress --exclude '__pycache__/' ${PIPELINE_HOME}/workflow/scripts/ $WORKDIR/scripts/
+  log_ok "$WORKDIR/scripts has been updated (excluding __pycache__/)."
+  log_next "Run dryrun to validate script updates before run/runlocal."
+}
+
+function syncscripts(){
+# sync workflow/scripts into WORKDIR/scripts so newly added scripts are available
+  rescript
+}
+
+function synccluster(){
+# sync cluster yaml into WORKDIR/config so scheduler changes can be propagated
+  check_essential_files
+  log_step "[synccluster] Syncing cluster config into workdir"
+  if [[ ! -f ${PIPELINE_HOME}/resources/$cluster_specific_yaml ]]; then
+    err "Cluster config ${PIPELINE_HOME}/resources/$cluster_specific_yaml not found. Set --cluster-config=<FILE> or CARLISLE_CLUSTER_CONFIG."
+  fi
+  cp "${PIPELINE_HOME}/resources/${cluster_specific_yaml}" "${WORKDIR}/config/cluster.yaml"
+  log_ok "${WORKDIR}/config/cluster.yaml has been updated from ${PIPELINE_HOME}/resources/${cluster_specific_yaml}."
+  log_next "Run dryrun to validate scheduler config updates before run/runlocal."
 }
 
 function runcheck(){
@@ -170,8 +232,8 @@ function controlcheck(){
 
   for sample_id in ${control_list[@]}; do
     if [[ $check1 == $sample_id || $check2 == $sample_id ]]; then
-      echo "Controls ($sample_id) cannot be listed in contrast.csv - update and re-run"
-      echo "$check1 okkk $check2"
+      log_error "Controls ($sample_id) cannot be listed in contrasts.tsv."
+      log_next "Update $WORKDIR/config/contrasts.tsv and re-run dryrun."
       exit 0
     fi
   done
@@ -179,6 +241,7 @@ function controlcheck(){
 
 function dryrun() {
 # Dry-run
+  log_step "[dryrun] Running preflight checks"
   runcheck
   module load $SINGULARITY_VERSION
   controlcheck
@@ -190,13 +253,16 @@ function dryrun() {
     mv ${WORKDIR}/dryrun.log ${WORKDIR}/logs/dryrun/dryrun.${modtime}.log
   fi
 
-  echo "========================================" | tee ${WORKDIR}/dryrun.log
-  echo "Singularity/Apptainer Cache Directory: $SING_CACHE_DIR" | tee -a ${WORKDIR}/dryrun.log
-  echo "SINGULARITY_CACHEDIR will be set to: $SING_CACHE_DIR" | tee -a ${WORKDIR}/dryrun.log
-  echo "========================================" | tee -a ${WORKDIR}/dryrun.log
+  echo "------------------------------------------------------------------" | tee ${WORKDIR}/dryrun.log
+  echo "STEP  [dryrun] Validating workflow DAG" | tee -a ${WORKDIR}/dryrun.log
+  echo "INFO  Singularity/Apptainer Cache Directory: $SING_CACHE_DIR" | tee -a ${WORKDIR}/dryrun.log
+  echo "INFO  SINGULARITY_CACHEDIR will be set to: $SING_CACHE_DIR" | tee -a ${WORKDIR}/dryrun.log
+  echo "------------------------------------------------------------------" | tee -a ${WORKDIR}/dryrun.log
   print_versions | tee -a ${WORKDIR}/dryrun.log
 
   run "--dry-run -r" | tee -a ${WORKDIR}/dryrun.log
+  log_ok "Dry-run completed."
+  log_next "Submit run with: carlisle --workdir=$WORKDIR --runmode=run"
 }
 
 function unlock() {
@@ -218,13 +284,14 @@ function runlocal() {
   set_singularity_binds
   if [ "$SLURM_JOB_ID" == "" ];then err "runlocal can only be done on an interactive node"; fi
   module load $SINGULARITY_VERSION
-  echo "========================================"
-  echo "Singularity/Apptainer Cache Directory: $SING_CACHE_DIR"
-  echo "SINGULARITY_CACHEDIR will be set to: $SING_CACHE_DIR"
-  echo "========================================"
+  log_step "[runlocal] Starting local interactive execution"
+  log_divider
+  log_info "Singularity/Apptainer Cache Directory: $SING_CACHE_DIR"
+  log_info "SINGULARITY_CACHEDIR will be set to: $SING_CACHE_DIR"
+  log_divider
   print_versions
 
-  run "--dry-run" && echo "Dry-run was successful .... starting local execution" && run "local"
+  run "--dry-run" && log_ok "Dry-run was successful. Starting local execution." && run "local"
 }
 
 function runtest() {
@@ -234,7 +301,7 @@ function runtest() {
   sed -e "s/PIPELINE_HOME/${PIPELINE_HOME//\//\\/}/g" ${PIPELINE_HOME}/.test/samples.test.tsv > $WORKDIR/config/samples.tsv
   cp ${PIPELINE_HOME}/.test/contrasts.test.tsv $WORKDIR/config/contrasts.tsv
   check_essential_files
-  run "--dry-run" && echo "Dry-run was successful .... starting slurm execution"
+  run "--dry-run" && log_ok "Dry-run was successful. Starting SLURM execution."
   runslurm
 }
 
@@ -243,26 +310,97 @@ function runslurm() {
   runcheck
   module load $SINGULARITY_VERSION
   set_singularity_binds
-  echo "========================================"
-  echo "Singularity/Apptainer Cache Directory: $SING_CACHE_DIR"
-  echo "SINGULARITY_CACHEDIR will be set to: $SING_CACHE_DIR"
-  echo "========================================"
+  log_divider
+  log_step "[run] Preparing SLURM submission"
+  log_info "CARLISLE Run Summary"
+  log_info "Mode:       SLURM"
+  log_info "Workdir:    $WORKDIR"
+  log_info "Partition:  $PARTITIONS"
+  log_info "Max jobs:   ${MAX_JOBS}"
+  log_info "Scheduler:  --max-jobs-per-second ${MAX_JOBS_PER_SECOND}, --max-status-checks-per-second ${MAX_STATUS_CHECKS_PER_SECOND}"
+  log_info "Log file:   ${WORKDIR}/logs/snakemake.log"
+  log_divider
+  log_info "Singularity/Apptainer Cache Directory: $SING_CACHE_DIR"
+  log_info "SINGULARITY_CACHEDIR will be set to: $SING_CACHE_DIR"
+  log_divider
   print_versions
-  run "--dry-run" && echo "Dry-run was successful .... submitting jobs to job-scheduler" && run "slurm"
+  run "--dry-run" && log_ok "Dry-run was successful. Submitting jobs to scheduler." && run "slurm"
 }
 
 function print_versions() {
 # Print versions of Snakemake and Singularity
-  echo "========================================"
-  echo "Tool Versions:"
-  snakemake --version 2>/dev/null && echo "Snakemake version checked" || echo "Snakemake version: unable to determine"
-  singularity --version 2>/dev/null && echo "Singularity version checked" || echo "Singularity/Apptainer version: UNAVAILABLE"
-  echo "========================================"
+  log_divider
+  echo "INFO  Tool Versions:"
+  snakemake --version 2>/dev/null && echo "OK    Snakemake version checked" || echo "WARN  Snakemake version: unable to determine"
+  singularity --version 2>/dev/null && echo "OK    Singularity version checked" || echo "WARN  Singularity/Apptainer version: UNAVAILABLE"
+  log_divider
+}
+
+function json_escape() {
+# Escape a string for safe embedding in JSON values.
+  printf '%s' "$1" | sed 's/\\/\\\\/g; s/"/\\"/g'
+}
+
+function write_pipeline_state_marker() {
+# Write a single state marker file and status sidecar in WORKDIR.
+# Args: state reason [job_id]
+  local state="$1"
+  local reason="${2:-}"
+  local job_id="${3:-NA}"
+  local now_utc
+  local host
+  local sidecar
+  local sidecar_tmp
+  local submission_ts=""
+
+  now_utc=$(date -u +"%Y-%m-%dT%H:%M:%SZ")
+  host=$(hostname 2>/dev/null || echo "unknown")
+  sidecar="${WORKDIR}/pipeline.status.json"
+  sidecar_tmp="${sidecar}.tmp.$$"
+
+  # Persist submission timestamp on first running write; read back on subsequent writes.
+  if [[ "${state}" == "running" && "${reason}" == "submission_started" ]]; then
+    submission_ts="${now_utc}"
+    date +%s > "${WORKDIR}/pipeline.submit_epoch"
+  elif [[ -f "${WORKDIR}/pipeline.submit_epoch" ]]; then
+    local submit_epoch
+    submit_epoch=$(cat "${WORKDIR}/pipeline.submit_epoch" 2>/dev/null || true)
+    if [[ -n "${submit_epoch}" ]]; then
+      submission_ts=$(date -u -d "@${submit_epoch}" +"%Y-%m-%dT%H:%M:%SZ" 2>/dev/null || true)
+    fi
+  fi
+
+  rm -f "${WORKDIR}/pipeline.running" "${WORKDIR}/pipeline.completed" "${WORKDIR}/pipeline.failed" "${WORKDIR}/pipeline.canceled"
+  : > "${WORKDIR}/pipeline.${state}"
+
+  cat > "$sidecar_tmp" << EOF
+{
+  "pipeline": "CARLISLE",
+  "version": "${PIPELINE_VERSION}",
+  "state": "$(json_escape "$state")",
+  "reason": "$(json_escape "$reason")",
+  "runmode": "$(json_escape "${RUNMODE:-unknown}")",
+  "workdir": "$(json_escape "${WORKDIR}")",
+  "user": "$(json_escape "${USER:-unknown}")",
+  "slurm_job_id": "$(json_escape "$job_id")",
+  "host": "$(json_escape "$host")",
+  "submission_timestamp_utc": "$(json_escape "${submission_ts}")",
+  "start_timestamp_utc": null,
+  "duration_seconds": null,
+  "exit_code": null,
+  "tasks_done": null,
+  "tasks_total": null,
+  "snakemake_log": "$(json_escape "${WORKDIR}/logs/snakemake.log")",
+  "timestamp_utc": "${now_utc}"
+}
+EOF
+  mv "$sidecar_tmp" "$sidecar"
+  log_info "State marker updated: pipeline.${state} (reason=${reason}, slurm_job_id=${job_id})"
 }
 
 function preruncleanup() {
 # Cleanup function to rename/move files related to older runs to prevent overwriting them.
-  echo "Running..."
+  log_step "[run] Preparing workdir for a new execution"
 
   # check initialization
   check_essential_files
@@ -288,6 +426,24 @@ function preruncleanup() {
 
 }
 
+function run_jobby_best_effort() {
+# Best-effort jobby generation outside Snakemake callbacks.
+# Ensures a deterministic artifact in logs/ even if submission fails before Snakemake starts.
+  mkdir -p "${WORKDIR}/logs"
+
+  if ! command -v jobby >/dev/null 2>&1; then
+    module load ccbrpipeliner >/dev/null 2>&1 || true
+  fi
+
+  if [[ -f "${WORKDIR}/logs/snakemake.log" ]] && command -v jobby >/dev/null 2>&1; then
+    if ! jobby --tsv "${WORKDIR}/logs/snakemake.log" | tee "${WORKDIR}/logs/snakemake.log.jobby" >/dev/null; then
+      printf '%s\n' "jobby failed while parsing ${WORKDIR}/logs/snakemake.log" > "${WORKDIR}/logs/snakemake.log.jobby"
+    fi
+  else
+    printf '%s\n' "jobby unavailable or snakemake.log missing; submission failed before Snakemake callbacks" > "${WORKDIR}/logs/snakemake.log.jobby"
+  fi
+}
+
 function run() {
 # RUN function
 # argument1 can be:
@@ -301,13 +457,12 @@ function run() {
   if [ "$1" == "local" ]; then
     preruncleanup
 
-    echo ""
-    echo "========================================"
-    echo "Singularity/Apptainer Cache Directory: $SING_CACHE_DIR"
-    echo "SINGULARITY_CACHEDIR will be set to: $SING_CACHE_DIR"
-    echo "========================================"
+    log_step "[runlocal] Launching local Snakemake execution"
+    log_divider
+    log_info "Singularity/Apptainer Cache Directory: $SING_CACHE_DIR"
+    log_info "SINGULARITY_CACHEDIR will be set to: $SING_CACHE_DIR"
+    log_divider
     print_versions
-    echo ""
 
     $EXPORT_SING_CACHE_DIR_CMD
 
@@ -339,13 +494,12 @@ function run() {
 
     preruncleanup
 
-    echo ""
-    echo "========================================"
-    echo "Singularity/Apptainer Cache Directory: $SING_CACHE_DIR"
-    echo "SINGULARITY_CACHEDIR will be set to: $SING_CACHE_DIR"
-    echo "========================================"
+    log_step "[run] Building submit script and submitting to scheduler"
+    log_divider
+    log_info "Singularity/Apptainer Cache Directory: $SING_CACHE_DIR"
+    log_info "SINGULARITY_CACHEDIR will be set to: $SING_CACHE_DIR"
+    log_divider
     print_versions
-    echo ""
 
     REWRITE_SUBMIT_SCRIPT="false"
     if [[ ! -f ${WORKDIR}/submit_script.sbatch ]]; then
@@ -353,7 +507,22 @@ function run() {
     elif ! grep -q -- "--max-jobs-per-second" "${WORKDIR}/submit_script.sbatch"; then
       backup_file="${WORKDIR}/submit_script.sbatch.$(date +%Y%m%d%H%M%S).bak"
       cp "${WORKDIR}/submit_script.sbatch" "${backup_file}"
-      echo "Existing submit script lacked scheduler throttling flags; backed up to ${backup_file}"
+      log_warn "Existing submit script lacked scheduler throttling flags; backed up to ${backup_file}"
+      REWRITE_SUBMIT_SCRIPT="true"
+    elif ! grep -q -- "_write_pipeline_state_marker" "${WORKDIR}/submit_script.sbatch"; then
+      backup_file="${WORKDIR}/submit_script.sbatch.$(date +%Y%m%d%H%M%S).bak"
+      cp "${WORKDIR}/submit_script.sbatch" "${backup_file}"
+      log_warn "Existing submit script lacked pipeline state marker logic; backed up to ${backup_file}"
+      REWRITE_SUBMIT_SCRIPT="true"
+    elif ! grep -q -- "_progress_monitor" "${WORKDIR}/submit_script.sbatch"; then
+      backup_file="${WORKDIR}/submit_script.sbatch.$(date +%Y%m%d%H%M%S).bak"
+      cp "${WORKDIR}/submit_script.sbatch" "${backup_file}"
+      log_warn "Existing submit script lacked progress monitor; backed up to ${backup_file}"
+      REWRITE_SUBMIT_SCRIPT="true"
+    elif ! grep -q -- "submission_timestamp_utc" "${WORKDIR}/submit_script.sbatch"; then
+      backup_file="${WORKDIR}/submit_script.sbatch.$(date +%Y%m%d%H%M%S).bak"
+      cp "${WORKDIR}/submit_script.sbatch" "${backup_file}"
+      log_warn "Existing submit script lacked enriched status fields; backed up to ${backup_file}"
       REWRITE_SUBMIT_SCRIPT="true"
     fi
 
@@ -366,30 +535,152 @@ function run() {
 #SBATCH --time=96:00:00
 #SBATCH --cpus-per-task=2
 
+    set -euo pipefail
+
+    function _json_escape() {
+      printf '%s' "\$1" | sed 's/\\/\\\\/g; s/"/\\"/g'
+    }
+
+    function _write_pipeline_state_marker() {
+      local state="\$1"
+      local reason="\${2:-}"
+      local job_id="\${3:-\${SLURM_JOB_ID:-NA}}"
+      local exit_code="\${4:-null}"
+      local now_utc
+      local host
+      local sidecar="${WORKDIR}/pipeline.status.json"
+      local sidecar_tmp="\${sidecar}.tmp.\$$"
+      local submission_ts=""
+      local duration_seconds=null
+      local tasks_done=null
+      local tasks_total=null
+      local _raw_tasks=""
+      local _submit_epoch=""
+
+      now_utc=\$(date -u +"%Y-%m-%dT%H:%M:%SZ")
+      host=\$(hostname 2>/dev/null || echo "unknown")
+
+      # Read submission timestamp from helper file written at sbatch submission time.
+      if [[ -f "${WORKDIR}/pipeline.submit_epoch" ]]; then
+        _submit_epoch=\$(cat "${WORKDIR}/pipeline.submit_epoch" 2>/dev/null || true)
+        if [[ -n "\${_submit_epoch}" ]]; then
+          submission_ts=\$(date -u -d "@\${_submit_epoch}" +"%Y-%m-%dT%H:%M:%SZ" 2>/dev/null || true)
+        fi
+      fi
+
+      # Compute wall-clock duration from job start.
+      if [[ -n "\${_START_EPOCH:-}" ]]; then
+        duration_seconds=\$(( \$(date +%s) - _START_EPOCH ))
+      fi
+
+      # Parse final task counts from snakemake.log.
+      _raw_tasks=\$(grep -oP '[0-9]+ of [0-9]+ steps \([0-9]+%\) done' "${WORKDIR}/logs/snakemake.log" 2>/dev/null | tail -1 || true)
+      if [[ -n "\${_raw_tasks}" ]]; then
+        tasks_done=\$(printf '%s' "\${_raw_tasks}" | grep -oP '^\d+' || echo null)
+        tasks_total=\$(printf '%s' "\${_raw_tasks}" | grep -oP '(?<=of )\d+' || echo null)
+      fi
+
+      rm -f "${WORKDIR}/pipeline.running" "${WORKDIR}/pipeline.completed" "${WORKDIR}/pipeline.failed" "${WORKDIR}/pipeline.canceled"
+      : > "${WORKDIR}/pipeline.\${state}"
+
+      cat > "\${sidecar_tmp}" << JSON
+{
+  "pipeline": "CARLISLE",
+  "version": "${PIPELINE_VERSION}",
+  "state": "\$(_json_escape "\${state}")",
+  "reason": "\$(_json_escape "\${reason}")",
+  "runmode": "run",
+  "workdir": "\$(_json_escape "${WORKDIR}")",
+  "user": "\$(_json_escape "\${USER:-unknown}")",
+  "slurm_job_id": "\$(_json_escape "\${job_id}")",
+  "host": "\$(_json_escape "\${host}")",
+  "submission_timestamp_utc": "\$(_json_escape "\${submission_ts}")",
+  "start_timestamp_utc": "\$(_json_escape "\${_START_TS:-}")",
+  "duration_seconds": \${duration_seconds},
+  "exit_code": \${exit_code},
+  "tasks_done": \${tasks_done},
+  "tasks_total": \${tasks_total},
+  "snakemake_log": "\$(_json_escape "${WORKDIR}/logs/snakemake.log")",
+  "timestamp_utc": "\${now_utc}"
+}
+JSON
+      mv "\${sidecar_tmp}" "\${sidecar}"
+    }
+
+    trap '_write_pipeline_state_marker canceled "signal_SIGTERM" "\${SLURM_JOB_ID:-NA}" 130; exit 130' SIGTERM
+    trap '_write_pipeline_state_marker canceled "signal_SIGINT" "\${SLURM_JOB_ID:-NA}" 130; exit 130' SIGINT
+
     module load $SNAKEMAKE_VERSION
     module load $SINGULARITY_VERSION
 
     cd \$SLURM_SUBMIT_DIR
 
     # Diagnostic output: show which modules loaded and check tool availability
-    echo "========================================"
-    echo "Loaded modules:"
+    echo "------------------------------------------------------------------"
+    echo "STEP  [submit_script] Environment diagnostics"
+    echo "INFO  Loaded modules:"
     module list 2>&1 | head -20
-    echo ""
-    echo "Tool availability check:"
-    echo "Snakemake: \$(which snakemake)"
-    echo "Singularity: \$(which singularity)"
-    echo "========================================"
-    echo ""
+    echo "INFO  Tool availability check:"
+    echo "INFO  Snakemake: \$(which snakemake)"
+    echo "INFO  Singularity: \$(which singularity)"
+    echo "------------------------------------------------------------------"
 
     $EXPORT_SING_CACHE_DIR_CMD
 
-    echo ""
-    echo "========================================"
-    echo "Singularity/Apptainer Cache Directory: $SING_CACHE_DIR"
-    echo "SINGULARITY_CACHEDIR has been set to: $SING_CACHE_DIR"
-    echo "========================================"
-    echo ""
+    echo "------------------------------------------------------------------"
+    echo "INFO  Singularity/Apptainer Cache Directory: $SING_CACHE_DIR"
+    echo "INFO  SINGULARITY_CACHEDIR has been set to: $SING_CACHE_DIR"
+    echo "------------------------------------------------------------------"
+
+    # Record job start time for duration tracking in pipeline.status.json.
+    _START_EPOCH=\$(date +%s)
+    _START_TS=\$(date -u +"%Y-%m-%dT%H:%M:%SZ")
+
+    # ── Progress monitor ─────────────────────────────────────────────────
+    # Background loop: waits for snakemake.log, then periodically writes
+    # a human-readable progress summary into pipeline.running.
+    _progress_monitor() {
+      local logfile="${WORKDIR}/logs/snakemake.log"
+      local marker="${WORKDIR}/pipeline.running"
+      local waited=0
+      local max_wait=300
+      local raw done_n total_n pct remaining updated
+
+      while [[ ! -f "\$logfile" && "\$waited" -lt "\$max_wait" ]]; do
+        [[ -f "\$marker" ]] || return 0
+        sleep 10
+        waited=\$(( waited + 10 ))
+      done
+
+      if [[ ! -f "\$logfile" ]]; then
+        printf 'Status   : Waiting for snakemake.log\nUpdated  : %s\n' "\$(date '+%Y-%m-%d %H:%M:%S')" > "\$marker"
+      fi
+
+      while true; do
+        [[ -f "\$marker" ]] || break
+        if [[ -f "\$logfile" ]]; then
+          raw=\$(grep -oP '[0-9]+ of [0-9]+ steps \([0-9]+%\) done' "\$logfile" 2>/dev/null | tail -1 || true)
+          if [[ -n "\$raw" ]]; then
+            done_n=\$(printf '%s' "\$raw" | grep -oP '^\d+' || true)
+            total_n=\$(printf '%s' "\$raw" | grep -oP '(?<=of )\d+' || true)
+            pct=\$(printf '%s' "\$raw" | grep -oP '\d+(?=%)' || true)
+            if [[ -n "\$done_n" && -n "\$total_n" && -n "\$pct" ]]; then
+              remaining=\$(( total_n - done_n ))
+              (( remaining < 0 )) && remaining=0
+              updated=\$(date '+%Y-%m-%d %H:%M:%S')
+              printf 'Progress : %s / %s steps complete (%s%%)\nRemaining: %s steps\nUpdated  : %s\n' \
+                "\$done_n" "\$total_n" "\$pct" "\$remaining" "\$updated" > "\$marker"
+            fi
+          elif [[ ! -s "\$marker" ]]; then
+            printf 'Status   : Submitted, waiting for first progress update\nUpdated  : %s\n' "\$(date '+%Y-%m-%d %H:%M:%S')" > "\$marker"
+          fi
+        fi
+        sleep 60
+      done
+    }
+    _progress_monitor &
+    _monitor_pid=\$!
+    # ─────────────────────────────────────────────────────────────────────
 
     snakemake -s $SNAKEFILE \
     --directory $WORKDIR \
@@ -410,20 +701,59 @@ function run() {
     --stats ${WORKDIR}/logs/snakemake.stats \
     2>&1|tee ${WORKDIR}/logs/snakemake.log
 
-    if [ "\$?" -eq "0" ];then
+    snakemake_rc="\$?"
+    kill "\$_monitor_pid" 2>/dev/null || true
+    wait "\$_monitor_pid" 2>/dev/null || true
+    if [ "\${snakemake_rc}" -eq "0" ];then
       snakemake -s $SNAKEFILE \
       --directory $WORKDIR \
       --report ${WORKDIR}/logs/runslurm_snakemake_report.html \
       --use-singularity \
       --singularity-prefix "$SING_CACHE_DIR" \
       --configfile ${WORKDIR}/config/config.yaml
+      report_rc="\$?"
+      if [ "\${report_rc}" -eq "0" ]; then
+        _write_pipeline_state_marker completed "snakemake_and_report_succeeded" "\${SLURM_JOB_ID:-NA}" 0
+      else
+        _write_pipeline_state_marker failed "report_generation_failed" "\${SLURM_JOB_ID:-NA}" "\${report_rc}"
+        exit "\${report_rc}"
+      fi
+    else
+      _write_pipeline_state_marker failed "snakemake_failed" "\${SLURM_JOB_ID:-NA}" "\${snakemake_rc}"
+      exit "\${snakemake_rc}"
     fi
 
 EOF
-      echo "Created ${WORKDIR}/submit_script.sbatch"
+  log_ok "Created ${WORKDIR}/submit_script.sbatch"
+    fi
+
+    # Fresh slurm run in this workdir: rotate prior marker to running.
+    write_pipeline_state_marker running "submission_started"
+
+    if ! sbatch_output=$(sbatch --parsable ${WORKDIR}/submit_script.sbatch 2>&1); then
+      write_pipeline_state_marker failed "sbatch_submission_failed"
+      run_jobby_best_effort
+      log_error "Failed to submit ${WORKDIR}/submit_script.sbatch"
+      log_error "sbatch output: ${sbatch_output}"
+      log_next "Check ${WORKDIR}/submit_script.sbatch and re-run: carlisle --workdir=$WORKDIR --runmode=run"
+      exit 1
+    fi
+
+    # --parsable outputs "JOBID" or "JOBID;CLUSTER" — take the first field
+    sbatch_job_id=$(echo "${sbatch_output}" | cut -d';' -f1 | tr -d '[:space:]')
+    if [[ -z "${sbatch_job_id}" ]]; then
+      sbatch_job_id="NA"
     fi
+    write_pipeline_state_marker running "sbatch_submitted" "${sbatch_job_id}"
 
-    sbatch ${WORKDIR}/submit_script.sbatch
+    log_divider
+    log_ok "Job submitted successfully (SLURM job ID: ${sbatch_job_id})"
+    log_info "Submission output: ${sbatch_output}"
+    log_next "Monitor:    squeue -u \$USER"
+    log_next "Progress:   tail -f ${WORKDIR}/logs/snakemake.log"
+    log_next "Status:     ls -1 ${WORKDIR}/pipeline.*"
+    log_next "Sidecar:    ${WORKDIR}/pipeline.status.json"
+    log_divider
 
   else # for unlock and dryrun
     snakemake $1 -s $SNAKEFILE \
@@ -447,13 +777,13 @@ EOF
 
 function reset() {
 # Delete the workdir and re-initialize it
-  echo "Working Dir: $WORKDIR"
+  log_info "Working Dir: $WORKDIR"
 
   if [ ! -d $WORKDIR ];then err "Folder $WORKDIR does not exist!";fi
 
-  echo "Deleting $WORKDIR"
+  log_warn "Deleting $WORKDIR"
   rm -rf $WORKDIR
-  echo "Re-Initializing $WORKDIR"
+  log_step "Re-initializing $WORKDIR"
   init
 }
 
@@ -481,6 +811,9 @@ function main(){
         -c=*|--singcache=*)
           SING_CACHE_DIR="${i#*=}"
         ;;
+        --cluster-config=*)
+          cluster_specific_yaml="${i#*=}"
+        ;;
         -v|--version)
           get_version && exit 0
         ;;
@@ -493,31 +826,35 @@ function main(){
   done
 
   WORKDIR=$(readlink -m "$WORKDIR")
-  echo "Working Dir: $WORKDIR"
+  log_info "Pipeline Dir: $PIPELINE_HOME"
+  log_info "Version: $PIPELINE_VERSION"
+  log_info "Working Dir: $WORKDIR"
 
   if [[ -z "$SING_CACHE_DIR" ]]; then
     if [[ -n "$SIFCACHE" ]]; then
       SING_CACHE_DIR="$SIFCACHE"
-    elif [[ -d "/data/$USER" ]]; then
-      SING_CACHE_DIR="/data/$USER/.singularity"
+    elif [[ -d "/data/CCBR_Pipeliner/SIFs" ]] && [[ -w "/data/CCBR_Pipeliner/SIFs" ]]; then
+      SING_CACHE_DIR="/data/CCBR_Pipeliner/SIFs"
     else
-      SING_CACHE_DIR="${WORKDIR}/.singularity"
+      SING_CACHE_DIR="/data/$USER/.singularity"
     fi
-    echo "singularity cache dir (--singcache) is not set, using ${SING_CACHE_DIR}"
+    log_info "singularity cache dir (--singcache) is not set, using ${SING_CACHE_DIR}"
   fi
   mkdir -p $SING_CACHE_DIR
   EXPORT_SING_CACHE_DIR_CMD="export SINGULARITY_CACHEDIR=\"${SING_CACHE_DIR}\""
-  echo "Snakemake scheduler limits: -j ${MAX_JOBS}, --max-jobs-per-second ${MAX_JOBS_PER_SECOND}, --max-status-checks-per-second ${MAX_STATUS_CHECKS_PER_SECOND}"
+  log_info "Snakemake scheduler limits: -j ${MAX_JOBS}, --max-jobs-per-second ${MAX_JOBS_PER_SECOND}, --max-status-checks-per-second ${MAX_STATUS_CHECKS_PER_SECOND}"
 
   case $RUNMODE in
     init) init && exit 0;;
+    synccluster) synccluster && exit 0;;
+    syncscripts) syncscripts && exit 0;;
     dryrun) dryrun && exit 0;;
     unlock) unlock && exit 0;;
     run) runslurm && exit 0;;
     runlocal) runlocal && exit 0;;
     reset) reset && exit 0;;
     runtest) runtest && exit 0;;
-    rescript) rescript && exit 0;;               # hidden option for debugging
+    rescript) rescript && exit 0;;               # backward-compatible alias
     printbinds) printbinds && exit 0;;           # hidden option
     *) err "Unknown RUNMODE \"$RUNMODE\"";;
   esac
diff --git a/docs/css/custom.css b/docs/css/custom.css
index c4c8782c..5b78b9e5 100644
--- a/docs/css/custom.css
+++ b/docs/css/custom.css
@@ -44,3 +44,9 @@ div.admonition.warning {
     margin: 10px 0;
     border-radius: 5px;
 }
+
+.docs-last-updated {
+    margin-top: 18px;
+    font-size: 0.92rem;
+    color: #5f6f7a;
+}
diff --git a/docs/overrides/main.html b/docs/overrides/main.html
index c59aff39..dfe09b34 100644
--- a/docs/overrides/main.html
+++ b/docs/overrides/main.html
@@ -1,2 +1,11 @@
 {% extends "base.html" %}
 {% block breadcrumbs %}{% endblock %}
+
+{% block content %}
+	{{ super() }}
+	{% if page.title == "Home" %}
+		<p class="docs-last-updated">
+			Last updated: {{ build_date_utc.strftime("%Y-%m-%d") }}
+		</p>
+	{% endif %}
+{% endblock %}
diff --git a/docs/user-guide/output.md b/docs/user-guide/output.md
index 771b5b63..9fe657c7 100644
--- a/docs/user-guide/output.md
+++ b/docs/user-guide/output.md
@@ -37,10 +37,17 @@ Upon successful completion, CARLISLE generates a comprehensive directory structu
       - **`peak_output/`** – Raw peak calls organized by control mode:
         - **`individual/`** – Peaks called using individual replicate controls (present for all analyses)
         - **`pooled/`** – Peaks called using merged high-depth controls (present when `pool_controls: true`)
+        - **`treatment_merged/`** – Consensus peak sets merged across replicates for each treatment sample (no re-calling), requiring support from at least 2 replicates.
       - **`annotation/`** – Contains enriched feature and pathway analyses, organized by control mode:
         - **`go_enrichment/`** – Results from **[ChIP-Enrich](https://chipenrich.med.umich.edu/)** gene set enrichment, generated when `run_go_enrichment: true` is enabled.
         - **`homer/`** – Output from **[HOMER](http://homer.ucsd.edu/homer/)** motif discovery and annotation.
         - **`rose/`** – Output from **[ROSE](https://github.com/younglab/ROSE)** super-enhancer analysis, generated when `run_rose: true` is specified.
+        - **`homer_treatment/`** – HOMER/AME outputs run on treatment-level merged peak sets.
+        - **`rose_treatment/`** – ROSE outputs run on treatment-level merged peak sets and merged treatment BAMs.
+
+  - **`annotation/homer_treatment/all_treatments.annotation_summary.tsv`** – Global HOMER annotation summary aggregated across all treatment-level runs.
+  - **`annotation/rose_treatment/all_treatments.SuperEnhancers_ENHANCER_TO_GENE.txt`** – Global ROSE enhancer-to-gene mapping aggregated across all treatment-level runs.
+  - **`annotation/rose_treatment/all_treatments.SuperEnhancers_GENE_TO_ENHANCER.txt`** – Global ROSE gene-to-enhancer mapping aggregated across all treatment-level runs.
 
 - **`deeptools/`** – Genome-wide sample correlation and PCA outputs, generated from 10 kb bin read counts across all samples. Two sets of files are produced:
 
@@ -99,44 +106,59 @@ results/
 │   │   ├── gopeaks/
 │   │   │   ├── peak_output/
 │   │   │   │   ├── individual/    # Peaks with individual controls
-│   │   │   │   └── pooled/        # Peaks with pooled controls
+│   │   │   │   ├── pooled/        # Peaks with pooled controls
+│   │   │   │   └── treatment_merged/  # Consensus peaks merged across treatment replicates
 │   │   │   └── annotation/
 │   │   │       ├── individual/
 │   │   │       │   ├── go_enrichment/
 │   │   │       │   ├── homer/
 │   │   │       │   └── rose/
-│   │   │       └── pooled/
-│   │   │           ├── go_enrichment/
-│   │   │           ├── homer/
-│   │   │           └── rose/
+│   │   │       ├── pooled/
+│   │   │       │   ├── go_enrichment/
+│   │   │       │   ├── homer/
+│   │   │       │   └── rose/
+│   │   │       ├── homer_treatment/
+│   │   │       └── rose_treatment/
 │   │   ├── macs2/
 │   │   │   ├── peak_output/
 │   │   │   │   ├── individual/
-│   │   │   │   └── pooled/
+│   │   │   │   ├── pooled/
+│   │   │   │   └── treatment_merged/
 │   │   │   └── annotation/
 │   │   │       ├── individual/
 │   │   │       │   ├── go_enrichment/
 │   │   │       │   ├── homer/
 │   │   │       │   └── rose/
-│   │   │       └── pooled/
-│   │   │           ├── go_enrichment/
-│   │   │           ├── homer/
-│   │   │           └── rose/
+│   │   │       ├── pooled/
+│   │   │       │   ├── go_enrichment/
+│   │   │       │   ├── homer/
+│   │   │       │   └── rose/
+│   │   │       ├── homer_treatment/
+│   │   │       └── rose_treatment/
 │   │   └── seacr/
 │   │       ├── peak_output/
 │   │       │   ├── individual/
-│   │       │   └── pooled/
+│   │       │   ├── pooled/
+│   │       │   └── treatment_merged/
 │   │       └── annotation/
 │   │           ├── individual/
 │   │           │   ├── go_enrichment/
 │   │           │   ├── homer/
 │   │           │   └── rose/
-│   │           └── pooled/
-│   │               ├── go_enrichment/
-│   │               ├── homer/
-│   │               └── rose/
+│   │           ├── pooled/
+│   │           │   ├── go_enrichment/
+│   │           │   ├── homer/
+│   │           │   └── rose/
+│   │           ├── homer_treatment/
+│   │           └── rose_treatment/
 │   └── 0.01/
 │       └── ...
+│   └── annotation/
+│       ├── homer_treatment/
+│       │   └── all_treatments.annotation_summary.tsv
+│       └── rose_treatment/
+│           ├── all_treatments.SuperEnhancers_ENHANCER_TO_GENE.txt
+│           └── all_treatments.SuperEnhancers_GENE_TO_ENHANCER.txt
 └── qc/
     ├── fastqc_raw/
     ├── fqscreen_raw/
@@ -242,6 +264,23 @@ Key columns to interpret:
 
 `rose/` identifies super-enhancers by stitching H3K27ac (or other active mark) peaks within `stitch_distance` (default 12,500 bp), then applying an inflection-point cutoff on the ranked signal curve. The inflection point separates typical enhancers from super-enhancers.
 
+- **`rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt`**: enhancer-centric table generated by Younglab `ROSE_geneMapper.py` with gene linkage columns (`OVERLAP_GENES`, `PROXIMAL_GENES`, `CLOSEST_GENE`). Use this file to map each super-enhancer to candidate target genes.
+- **`rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt`**: gene-centric view of the same mapping, listing genes and their associated super-enhancer regions.
 - **`*_SuperEnhancers.bed`**: the final super-enhancer calls. Overlap with known oncogenes or lineage TFs to nominate driver regulatory elements.
 - **`*_Enhancers_withSuper.bed`**: all stitched enhancers ranked by signal; use this to inspect the inflection point and evaluate whether the cutoff is biologically reasonable.
 - ROSE is only meaningful for **active chromatin marks** (H3K27ac, H3K4me1). Running it on H3K27me3 or TF datasets is not informative.
+
+### Treatment-Level Merged Analysis
+
+CARLISLE now supports treatment-level analysis in parallel with replicate-level outputs:
+
+- **Consensus merged peak sets** (`peak_output/*/treatment_merged/`): generated from replicate peak calls without re-calling. A merged interval is kept when supported by at least 2 replicates.
+- **Coordinate policy**: merged consensus intervals span the leftest-left and rightest-right coordinates of overlapping replicate peaks.
+- **HOMER treatment outputs** (`annotation/homer_treatment/`): motif and annotation analyses run on merged treatment peak sets.
+- **ROSE treatment outputs** (`annotation/rose_treatment/`): super-enhancer analyses run on merged treatment peak sets and merged treatment BAMs.
+
+Global treatment-level aggregates are written to:
+
+- `peaks/annotation/homer_treatment/all_treatments.annotation_summary.tsv`
+- `peaks/annotation/rose_treatment/all_treatments.SuperEnhancers_ENHANCER_TO_GENE.txt`
+- `peaks/annotation/rose_treatment/all_treatments.SuperEnhancers_GENE_TO_ENHANCER.txt`
diff --git a/docs/user-guide/run.md b/docs/user-guide/run.md
index fef706dc..8d732c29 100644
--- a/docs/user-guide/run.md
+++ b/docs/user-guide/run.md
@@ -57,7 +57,7 @@ For cluster execution (`run`, `runtest`), CARLISLE uses scheduler-safe Snakemake
 
 - **`runlocal`** – Executes the workflow on a local interactive node. This mode is suitable for quick testing or smaller datasets but should only be used within a Biowulf interactive session (`sinteractive`). A Snakemake HTML report (`report.html`) is generated in the working directory after a successful run.
 
-- **`run`** – Submits the workflow to the **[Biowulf HPC cluster](https://hpc.nih.gov/)** via SLURM. CARLISLE manages job scheduling, dependencies, and notifications. Email alerts are automatically sent for job start, errors, and completion. The Singularity module is loaded automatically before submission — no manual `module load singularity` is required. A Snakemake HTML report (`report.html`) is generated in the working directory after a successful run.
+- **`run`** – Submits the workflow to the **[Biowulf HPC cluster](https://hpc.nih.gov/)** via SLURM. CARLISLE manages job scheduling and dependencies. The Singularity module is loaded automatically before submission — no manual `module load singularity` is required. A Snakemake HTML report (`report.html`) is generated in the working directory after a successful run.
 
 ### Maintenance Commands
 
@@ -168,6 +168,71 @@ carlisle --runmode=run --workdir=/path/to/output/dir
 
 ---
 
+## Understanding Log Output
+
+CARLISLE prints structured log lines to the terminal during `init`, `dryrun`, `run`, `runlocal`, and other modes. Each line is prefixed with a fixed-width tag that indicates its severity or purpose:
+
+| Prefix | Meaning | When you see it |
+|--------|---------|-----------------|
+| `STEP` | A major pipeline phase is starting | Beginning of `init`, `dryrun`, `run`, Snakemake submission, etc. |
+| `INFO` | Informational detail about the current step | Paths, settings, versions, tool availability |
+| `OK` | A step completed successfully | Files copied, modules loaded, job submitted without error |
+| `WARN` | A non-fatal issue that may need attention | Outdated submit script backed up, optional tool unavailable |
+| `ERROR` | A fatal problem — execution will stop | Missing files, failed module load, bad arguments |
+| `NEXT` | Suggested next action for the user | What to edit or run after the current step completes |
+
+### Example terminal output
+
+A typical `carlisle --runmode=run` invocation looks like this:
+
+```
+------------------------------------------------------------------
+STEP  [run] Preparing SLURM submission
+INFO  CARLISLE Run Summary
+INFO  Mode:       SLURM
+INFO  Workdir:    /data/$USER/project
+INFO  Partition:  norm
+INFO  Max jobs:   100
+INFO  Scheduler:  --max-jobs-per-second 1, --max-status-checks-per-second 0.1
+INFO  Log file:   /data/$USER/project/logs/snakemake.log
+------------------------------------------------------------------
+INFO  Tool Versions:
+7.32.4
+OK    Snakemake version checked
+apptainer version 1.3.6-1.bionic
+OK    Singularity version checked
+------------------------------------------------------------------
+OK    Dry-run was successful. Submitting jobs to scheduler.
+------------------------------------------------------------------
+OK    Job submitted successfully (SLURM job ID: 12345678)
+INFO  Submission output: 12345678
+NEXT  Monitor:    squeue -u $USER
+NEXT  Progress:   tail -f /data/$USER/project/logs/snakemake.log
+NEXT  Status:     ls -1 /data/$USER/project/pipeline.*
+NEXT  Sidecar:    /data/$USER/project/pipeline.status.json
+------------------------------------------------------------------
+```
+
+`NEXT` lines tell you exactly what to run after each stage — copy them directly into your terminal.
+
+### Live progress in `pipeline.running`
+
+While a `run` job is active, CARLISLE updates the `pipeline.running` marker file every 60 seconds with a human-readable progress summary parsed from `snakemake.log`:
+
+```
+Progress : 47 / 312 steps complete (15%)
+Remaining: 265 steps
+Updated  : 2026-07-17 14:32:10
+```
+
+Read it at any time with:
+
+```bash
+cat /path/to/output/dir/pipeline.running
+```
+
+---
+
 ## Monitoring a Running Job
 
 After submitting with `run`, CARLISLE itself exits immediately — the pipeline runs as a background SLURM job. To monitor progress:
@@ -181,12 +246,64 @@ watch -n 30 squeue -u $USER
 
 # View the Snakemake master job log (replace JOBID with the number from squeue)
 cat /path/to/output/dir/logs/snakemake.log
+
+# Check CARLISLE run-state marker (exactly one should exist)
+ls -1 /path/to/output/dir/pipeline.*
+
+# Inspect state metadata
+cat /path/to/output/dir/pipeline.status.json
+```
+
+`pipeline.status.json` contains the following fields:
+
+| Field | Description |
+|---|---|
+| `pipeline` | Pipeline name (`CARLISLE`) |
+| `version` | CARLISLE version at the time of submission |
+| `state` | Current state: `running`, `completed`, `failed`, or `canceled` |
+| `reason` | Machine-readable reason for the state (e.g. `snakemake_and_report_succeeded`, `snakemake_failed`) |
+| `runmode` | Run mode used (e.g. `run`) |
+| `workdir` | Absolute path to the working directory |
+| `user` | Username that submitted the pipeline |
+| `slurm_job_id` | SLURM job ID of the pipeline coordinator job |
+| `host` | Hostname where the state was last written |
+| `submission_timestamp_utc` | UTC timestamp when `carlisle --runmode=run` was invoked |
+| `start_timestamp_utc` | UTC timestamp when the SLURM job began executing |
+| `duration_seconds` | Wall-clock seconds from job start to final state (null while running or on headnode-only writes) |
+| `exit_code` | Exit code of Snakemake (0 = success; null for headnode-only writes) |
+| `tasks_done` | Number of Snakemake steps completed at final state (null if log is unavailable) |
+| `tasks_total` | Total number of Snakemake steps at final state (null if log is unavailable) |
+| `snakemake_log` | Absolute path to the Snakemake master log file |
+| `timestamp_utc` | UTC timestamp of the most recent state write |
+
+For `runmode=run`, CARLISLE writes exactly one state marker file in the workdir:
+
+- `pipeline.running` — run submitted and in progress
+- `pipeline.completed` — run finished successfully
+- `pipeline.failed` — submission or runtime failure
+- `pipeline.canceled` — run interrupted (for example `scancel`/signal)
+
+While the pipeline is running, `pipeline.running` is updated **every 60 seconds** with a live progress summary derived from `logs/snakemake.log`. You can check it at any time:
+
+```bash
+cat <workdir>/pipeline.running
+```
+
+Example output:
+
 ```
+Progress : 42 / 120 steps complete (35%)
+Remaining: 78 steps
+Updated  : 2026-07-17 14:23:01
+```
+
+If the job has just been submitted and Snakemake has not yet produced output, it will show:
 
-Email notifications are automatically sent to your NIH HPC account email (`$USER@nih.gov`) for:
+```
+Status   : Submitted, waiting for first progress update
+Updated  : 2026-07-17 14:05:00
+```
 
-- **Job start** — confirms the pipeline was accepted by SLURM
-- **Job error** — sent if any rule fails; check the log file above for details
-- **Job completion** — confirms all rules finished successfully
+When you start `runmode=run` again in the same workdir, CARLISLE removes any existing `pipeline.*` marker and replaces it with `pipeline.running` for the new run.
 
 After a successful run, a `report.html` file is generated in your working directory — open it in a browser for an interactive summary of all pipeline steps and outputs.
diff --git a/resources/cluster_biowulf.yaml b/resources/cluster_biowulf.yaml
index b743fc1e..2ef7385d 100644
--- a/resources/cluster_biowulf.yaml
+++ b/resources/cluster_biowulf.yaml
@@ -83,6 +83,11 @@ homer_motif:
   time: 02-00:00:00
   threads: 32
 ###################################################################
+homer_motif_treatment:
+  mem: 120g
+  time: 02-00:00:00
+  threads: 32
+###################################################################
 homer_motif_deg:
   mem: 120g
   time: 02-00:00:00
@@ -93,6 +98,11 @@ ame_motif_enrichment:
   time: 02-00:00:00
   threads: 32
 ###################################################################
+ame_motif_enrichment_treatment:
+  mem: 120g
+  time: 02-00:00:00
+  threads: 32
+###################################################################
 ame_motif_enrichment_deg:
   mem: 120g
   time: 02-00:00:00
diff --git a/workflow/Snakefile b/workflow/Snakefile
index be9a5dd2..787fe7c1 100644
--- a/workflow/Snakefile
+++ b/workflow/Snakefile
@@ -393,6 +393,8 @@ def get_rose(wildcards):
                 files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Gateway_Enhancers.bed"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M,s_dist=S_DISTANCE,control_mode=[control_mode]))
                 files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Gateway_SuperEnhancers.bed"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M,s_dist=S_DISTANCE,control_mode=[control_mode]))
                 files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_AllEnhancers.table.txt"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M,s_dist=S_DISTANCE,control_mode=[control_mode]))
+                files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M,s_dist=S_DISTANCE,control_mode=[control_mode]))
+                files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M,s_dist=S_DISTANCE,control_mode=[control_mode]))
 
             base_list_sg = TREATMENT_CONTROL_LIST_POOLED if control_mode == "pooled" else TREATMENT_LIST_SG
             tc_list_sg = _filter_tc_list(base_list_sg, control_mode)
@@ -402,11 +404,90 @@ def get_rose(wildcards):
                 files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Gateway_Enhancers.bed"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list_sg,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G,s_dist=S_DISTANCE,control_mode=[control_mode]))
                 files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Gateway_SuperEnhancers.bed"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list_sg,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G,s_dist=S_DISTANCE,control_mode=[control_mode]))
                 files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_AllEnhancers.table.txt"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list_sg,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G,s_dist=S_DISTANCE,control_mode=[control_mode]))
+                files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list_sg,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G,s_dist=S_DISTANCE,control_mode=[control_mode]))
+                files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list_sg,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G,s_dist=S_DISTANCE,control_mode=[control_mode]))
             if ("seacr_stringent" in PEAKTYPE) or ("seacr_relaxed" in PEAKTYPE):
                 files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Enhancers_withSuper.bed"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list_sg,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S,s_dist=S_DISTANCE,control_mode=[control_mode]))
                 files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Gateway_Enhancers.bed"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list_sg,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S,s_dist=S_DISTANCE,control_mode=[control_mode]))
                 files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Gateway_SuperEnhancers.bed"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list_sg,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S,s_dist=S_DISTANCE,control_mode=[control_mode]))
                 files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_AllEnhancers.table.txt"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list_sg,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S,s_dist=S_DISTANCE,control_mode=[control_mode]))
+                files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list_sg,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S,s_dist=S_DISTANCE,control_mode=[control_mode]))
+                files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt"),qthresholds=QTRESHOLDS,treatment_control_list=tc_list_sg,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S,s_dist=S_DISTANCE,control_mode=[control_mode]))
+    return files
+
+
+def get_treatment_merged_peak_sets(wildcards):
+    files = []
+    for peak_type in PEAKTYPE:
+        caller = peak_type.split("_", 1)[0]
+        files.extend(
+            expand(
+                join(
+                    RESULTSDIR,
+                    "peaks",
+                    "{qthresholds}",
+                    caller,
+                    "peak_output",
+                    "{control_mode}",
+                    "treatment_merged",
+                    "{treatment_sample}.{dupstatus}.{peak_type}.peaks.bed",
+                ),
+                qthresholds=QTRESHOLDS,
+                control_mode=CONTROL_MODES,
+                treatment_sample=TREATMENT_SAMPLES,
+                dupstatus=DUPSTATUS,
+                peak_type=[peak_type],
+            )
+        )
+    return files
+
+
+def get_treatment_motif_enrichment(wildcards):
+    files = []
+    if not config.get("run_motif_enrichment_called_peaks", False):
+        return files
+
+    if ("macs2_narrow" in PEAKTYPE) or ("macs2_broad" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.annotation.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","knownResults.html"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","ame_results.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M))
+    if ("gopeaks_narrow" in PEAKTYPE) or ("gopeaks_broad" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.annotation.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","knownResults.html"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","ame_results.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G))
+    if ("seacr_stringent" in PEAKTYPE) or ("seacr_relaxed" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.annotation.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","knownResults.html"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","ame_results.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S))
+
+    files.append(join(RESULTSDIR, "peaks", "annotation", "homer_treatment", "all_treatments.annotation_summary.tsv"))
+    return files
+
+
+def get_treatment_rose(wildcards):
+    files = []
+    if not config.get("run_rose", False):
+        return files
+    if GENOME not in ["hg19", "hg38", "mm10"]:
+        return files
+
+    if ("macs2_narrow" in PEAKTYPE) or ("macs2_broad" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Enhancers_withSuper.bed"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M,s_dist=S_DISTANCE))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M,s_dist=S_DISTANCE))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M,s_dist=S_DISTANCE))
+    if ("gopeaks_narrow" in PEAKTYPE) or ("gopeaks_broad" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Enhancers_withSuper.bed"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G,s_dist=S_DISTANCE))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G,s_dist=S_DISTANCE))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G,s_dist=S_DISTANCE))
+    if ("seacr_stringent" in PEAKTYPE) or ("seacr_relaxed" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Enhancers_withSuper.bed"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S,s_dist=S_DISTANCE))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S,s_dist=S_DISTANCE))
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S,s_dist=S_DISTANCE))
+
+    files.extend([
+        join(RESULTSDIR, "peaks", "annotation", "rose_treatment", "all_treatments.SuperEnhancers_ENHANCER_TO_GENE.txt"),
+        join(RESULTSDIR, "peaks", "annotation", "rose_treatment", "all_treatments.SuperEnhancers_GENE_TO_ENHANCER.txt"),
+    ])
     return files
 
 def get_go_enrichment(wildcards):
@@ -513,6 +594,9 @@ rule all:
         unpack(get_combined) if config.get("run_motif_enrichment_called_peaks", False) else [],
         unpack(get_deg_motif_enrichment) if config.get("run_motif_enrichment_deg_peaks", False) else [],
         unpack(get_rose),
+        unpack(get_treatment_merged_peak_sets),
+        unpack(get_treatment_motif_enrichment),
+        unpack(get_treatment_rose),
         unpack(get_go_enrichment)
 
 onsuccess:
@@ -553,6 +637,12 @@ if [ -f logs/spooker.log ]; then
 else
     echo "❌ logs/spooker.log NOT created."
 fi
+
+# Fallback state reconciliation if launcher-level transition did not run.
+if [ ! -f pipeline.canceled ]; then
+    rm -f pipeline.running pipeline.failed
+    : > pipeline.completed
+fi
 """)
 
 onerror:
@@ -593,6 +683,12 @@ if [ -f logs/spooker.log ]; then
 else
     echo "❌ logs/spooker.log NOT created."
 fi
+
+# Fallback state reconciliation if launcher-level transition did not run.
+if [ ! -f pipeline.canceled ]; then
+    rm -f pipeline.running pipeline.completed
+    : > pipeline.failed
+fi
 """)
 
 # """
diff --git a/workflow/rules/annotations.smk b/workflow/rules/annotations.smk
index 8f115c56..317d6552 100644
--- a/workflow/rules/annotations.smk
+++ b/workflow/rules/annotations.smk
@@ -1,4 +1,5 @@
 import re
+import os
 
 def get_peak_file(wildcards):
     # MACS2 OPTIONS
@@ -19,6 +20,41 @@ def get_peak_file(wildcards):
     if wildcards.peak_caller_type =="gopeaks_broad":
         bed=join(RESULTSDIR,"peaks",wildcards.qthresholds,"gopeaks","peak_output",wildcards.control_mode,wildcards.treatment_control_list + "." + wildcards.dupstatus + ".broad.peaks.bed")
     return bed
+
+
+def get_treatment_peak_file(wildcards):
+    if wildcards.peak_caller_type.startswith("macs2_"):
+        caller = "macs2"
+    elif wildcards.peak_caller_type.startswith("gopeaks_"):
+        caller = "gopeaks"
+    elif wildcards.peak_caller_type.startswith("seacr_"):
+        caller = "seacr"
+    else:
+        raise ValueError("Unsupported peak_caller_type for treatment peak file: %s" % wildcards.peak_caller_type)
+
+    bed = join(
+        RESULTSDIR,
+        "peaks",
+        wildcards.qthresholds,
+        caller,
+        "peak_output",
+        wildcards.control_mode,
+        "treatment_merged",
+        wildcards.treatment_sample + "." + wildcards.dupstatus + "." + wildcards.peak_caller_type + ".peaks.bed",
+    )
+    return bed
+
+
+def get_treatment_control_bam(wildcards):
+    control_sample = TREATMENT_SAMPLE_TO_CONTROL_SAMPLE.get(wildcards.treatment_sample, "nocontrol")
+    if control_sample == "nocontrol":
+        return []
+    return join(
+        RESULTSDIR,
+        "bam",
+        "pooled_controls",
+        control_sample + "." + wildcards.dupstatus + ".merged.bam",
+    )
 def get_deg_bed(wildcards):
     # DEG-based peak sets produced by diffbb
     # method: AUCbased | fragmentsbased
@@ -683,6 +719,152 @@ rule combine_homer:
         Rscript {params.rscript} {input.peaks_file} {input.annotation} {output.combined_tsv} {output.combined_xlsx}
         """
 
+
+rule homer_motif_treatment:
+    """
+    Run HOMER annotation and motif discovery on treatment-level merged peak sets.
+    """
+    input:
+        peak_file=get_treatment_peak_file,
+    output:
+        annotation=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.annotation.txt"),
+        annotation_summary=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.annotation.summary"),
+        known_html=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","knownResults.html"),
+        target_fasta=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","target.fa"),
+        background_fasta=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","background.fa"),
+    threads: getthreads("homer_motif")
+    envmodules:
+        TOOLS["homer"],
+    params:
+        genome=config["genome"],
+        fa=config["reference"][config["genome"]]["fa"],
+        gtf=config["reference"][config["genome"]]["gtf"],
+        outDir=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs"),
+        hocomoco_motif=config["hocomoco_motifs"],
+    shell:
+        """
+        set -euo pipefail
+
+        if [[ -d "/lscratch/$SLURM_JOB_ID" ]]; then
+            TMPDIR=$(mktemp -d "/lscratch/$SLURM_JOB_ID/tmp.XXXXXX")
+        else
+            TMPDIR=$(mktemp -d "/dev/shm/tmp.XXXXXX")
+        fi
+        preparsedDir="$TMPDIR/preparsedDir"
+        mkdir -p "$preparsedDir"
+
+        mkdir -p {params.outDir}
+
+        num_peaks=$(wc -l < {input.peak_file} || echo 0)
+        if [[ $num_peaks -lt 5 ]]; then
+            echo "# No peaks found for HOMER annotation" > {output.annotation}
+            echo -e "Annotation\tDistance to TSS\tNumber of Peaks\t% of Peaks\tTotal size (bp)\tLog10 p-value\tLog2 Ratio (vs. Genome)\tLogP enrichment (+values depleted)" > {output.annotation_summary}
+            echo "<html><body><h1>No peaks available for motif analysis</h1></body></html>" > {output.known_html}
+            : > {output.target_fasta}
+            : > {output.background_fasta}
+            exit 0
+        fi
+
+        if [[ {params.genome} == "hs1" ]]; then
+            annotatePeaks.pl {input.peak_file} {params.fa} -annStats {output.annotation_summary} -gtf {params.gtf} > {output.annotation}
+            findMotifsGenome.pl {input.peak_file} {params.fa} {params.outDir} \
+                -nomotif -size given -mknown {params.hocomoco_motif} -p {threads} \
+                -dumpFasta -cpg -maxN 0.1 -len 10 \
+                -preparsedDir "$preparsedDir"
+        else
+            annotatePeaks.pl {input.peak_file} {params.genome} -annStats {output.annotation_summary} > {output.annotation}
+            findMotifsGenome.pl {input.peak_file} {params.genome} {params.outDir} \
+                -nomotif -size given -mknown {params.hocomoco_motif} -p {threads} \
+                -dumpFasta -cpg -maxN 0.1 -len 10 \
+                -preparsedDir "$preparsedDir"
+        fi
+        """
+
+
+rule ame_motif_enrichment_treatment:
+    """
+    Run AME motif enrichment on treatment-level HOMER motif FASTAs.
+    """
+    input:
+        target_fasta=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","target.fa"),
+        background_fasta=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","background.fa"),
+    output:
+        ame_results=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs","ame_results.txt"),
+    threads: getthreads("ame_motif_enrichment")
+    envmodules:
+        TOOLS["parallel"],
+        TOOLS["meme"],
+    params:
+        outDir=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.motifs"),
+        hocomoco_memes_tar=config["hocomoco_memes_targz"],
+        python_script=join(SCRIPTSDIR,"_parse_ame_output.py"),
+    shell:
+        """
+        set -euo pipefail
+        cd {params.outDir}
+        if [[ ! -s target.fa || ! -s background.fa ]]; then
+            echo "# No sequences for AME analysis" > {output.ame_results}
+            exit 0
+        fi
+
+        rm -rf tmpdir
+        mkdir -p tmpdir
+        cp target.fa tmpdir/target.fa
+        cp background.fa tmpdir/background.fa
+        cd tmpdir
+
+        printf '%b\n' 'rank\tmotif_DB\tmotif_ID\tmotif_ALT_ID\tconsensus\tp-value\tadjusted-p-value\tE-value\ttests\tFAMP\tn_sequences\tTP\t%TP\tFP\t%FP' > {output.ame_results}
+        cp {params.hocomoco_memes_tar} .
+        tar xzf $(basename {params.hocomoco_memes_tar})
+        ls *.meme 2>/dev/null | sort > memes || true
+        if [[ -s memes ]]; then
+            while read a; do
+                echo "ame --o ${{a}}_ame_out --noseq --control background.fa --seed 12345 --verbose 3 target.fa ${{a}}"
+            done < memes > do_memes
+            parallel -j {threads} < do_memes
+            find . -name 'ame.tsv' -exec cat {{}} \; | \
+                grep -A1 ^rank | grep -v '^--$' | grep -v ^rank | sort | uniq | sort -k7,7g | \
+                python {params.python_script} >> {output.ame_results}
+        fi
+        cd {params.outDir}
+        rm -rf tmpdir
+        """
+
+
+def get_treatment_homer_annotation_summaries(wildcards):
+    files = []
+    if ("macs2_narrow" in PEAKTYPE) or ("macs2_broad" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.annotation.summary"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M))
+    if ("gopeaks_narrow" in PEAKTYPE) or ("gopeaks_broad" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.annotation.summary"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G))
+    if ("seacr_stringent" in PEAKTYPE) or ("seacr_relaxed" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","homer_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.annotation.summary"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S))
+    return files
+
+
+rule aggregate_homer_treatment:
+    """
+    Create one global HOMER annotation summary across all treatment-level runs.
+    """
+    input:
+        summaries=get_treatment_homer_annotation_summaries
+    output:
+        aggregate=join(RESULTSDIR,"peaks","annotation","homer_treatment","all_treatments.annotation_summary.tsv")
+    params:
+        script=join(SCRIPTSDIR, "_aggregate_homer_treatment.py"),
+        input_args=lambda w, input: " ".join([
+            "--input " + os.path.basename(p).replace(".annotation.summary", "") + "::" + p
+            for p in input.summaries
+        ])
+    envmodules:
+        TOOLS["python3"]
+    shell:
+        """
+        set -euo pipefail
+        mkdir -p "$(dirname "{output.aggregate}")"
+        python {params.script} {params.input_args} --output {output.aggregate}
+        """
+
 rule rose:
     """
     Run ROSE with a containerized two-step flow:
@@ -700,6 +882,8 @@ rule rose:
     - `rose_input_Gateway_Enhancers.bed`
     - `rose_input_Gateway_SuperEnhancers.bed`
     - `rose_input_AllEnhancers.table.txt`
+    - `rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt`
+    - `rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt`
     - Cleans up large ROSE intermediates (`gff/`, `mappedGFF/`) after success.
     """
     input:
@@ -726,6 +910,8 @@ rule rose:
         gateway_enhancers=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Gateway_Enhancers.bed"),
         gateway_super=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Gateway_SuperEnhancers.bed"),
         all_enhancers_table=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_AllEnhancers.table.txt"),
+        super_enhancer_to_gene=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt"),
+        super_gene_to_enhancer=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose","{control_mode}","{treatment_control_list}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt"),
     container: config["containers"].get("rose", "docker://nciccbr/ccbr_rose:v1")
     shell:
         """
@@ -799,6 +985,94 @@ rule rose:
                     -o {params.file_base}
             fi
 
+            # Run Younglab ROSE gene mapper on the super-enhancer table.
+            super_table={params.file_base}/rose_input_SuperEnhancers.table.txt
+            if [[ ! -f ROSE_geneMapper.py ]]; then
+                echo "ERROR: ROSE_geneMapper.py not found in {params.rose_root}."
+                exit 1
+            fi
+            if [[ ! -f "$super_table" ]]; then
+                echo "ERROR: Missing ROSE super-enhancer table required for gene mapping: $super_table"
+                exit 1
+            fi
+
+            # Normalize super-enhancer table schema for ROSE_geneMapper compatibility.
+            # Some ROSE outputs can omit enhancerRank/isSuper columns; geneMapper expects them.
+            normalized_super_table={params.file_base}/rose_input_SuperEnhancers.for_geneMapper.table.txt
+            {params.rose_python} - "$super_table" "$normalized_super_table" <<'PY'
+import sys
+
+
+def normalize_super_table(in_path, out_path):
+    newline = chr(10)
+    tab = chr(9)
+
+    with open(in_path, "r") as fin:
+        lines = fin.readlines()
+
+    header_idx = None
+    for i, raw in enumerate(lines):
+        if not raw.startswith("#") and raw.strip():
+            header_idx = i
+            break
+
+    if header_idx is None:
+        raise ValueError("No tabular header found in super-enhancer table")
+
+    header = lines[header_idx].rstrip(newline).split(tab)
+    has_rank = "enhancerRank" in header
+    has_super = "isSuper" in header
+
+    with open(out_path, "w") as fout:
+        for i in range(header_idx):
+            fout.write(lines[i])
+
+        out_header = list(header)
+        if not has_rank:
+            out_header.append("enhancerRank")
+        if not has_super:
+            out_header.append("isSuper")
+        fout.write(tab.join(out_header) + newline)
+
+        rank_counter = 0
+        for raw in lines[header_idx + 1 :]:
+            stripped = raw.rstrip(newline)
+            if not stripped:
+                continue
+            if stripped.startswith("#"):
+                fout.write(raw)
+                continue
+
+            rank_counter += 1
+            cols = stripped.split(tab)
+            if not has_rank:
+                cols.append(str(rank_counter))
+            if not has_super:
+                cols.append("1")
+            fout.write(tab.join(cols) + newline)
+
+    added_cols = []
+    if not has_rank:
+        added_cols.append("enhancerRank")
+    if not has_super:
+        added_cols.append("isSuper")
+
+    if added_cols:
+        print("INFO: Added missing ROSE SuperEnhancers columns for geneMapper: " + ", ".join(added_cols))
+    else:
+        print("INFO: ROSE SuperEnhancers schema already includes enhancerRank/isSuper")
+
+
+normalize_super_table(sys.argv[1], sys.argv[2])
+PY
+
+            {params.rose_python} ROSE_geneMapper.py \
+                -g {params.genome} \
+                -i "$normalized_super_table" \
+                -o {params.file_base}
+
+            rm -f "$normalized_super_table"
+
             # Cleanup large ROSE intermediates not needed downstream.
             rm -rf {params.file_base}/gff {params.file_base}/mappedGFF
         else
@@ -807,8 +1081,240 @@ rule rose:
             echo "Less than 5 usable peaks detected (N=${{num_of_peaks}})" > {output.gateway_enhancers}
             echo "Less than 5 usable peaks detected (N=${{num_of_peaks}})" > {output.gateway_super}
             echo "Less than 5 usable peaks detected (N=${{num_of_peaks}})" > {output.all_enhancers_table}
+            echo "Less than 5 usable peaks detected (N=${{num_of_peaks}})" > {output.super_enhancer_to_gene}
+            echo "Less than 5 usable peaks detected (N=${{num_of_peaks}})" > {output.super_gene_to_enhancer}
         fi
     """
+
+
+rule rose_treatment:
+    """
+    Run ROSE on treatment-level merged peak sets and merged treatment BAMs.
+    """
+    input:
+        peak_file=get_treatment_peak_file,
+        treatment_bam=join(RESULTSDIR,"bam","treatment_merged","{treatment_sample}.{dupstatus}.merged.bam"),
+        control_bam=get_treatment_control_bam,
+    threads: getthreads("rose")
+    params:
+        genome=config["genome"],
+        tss_bed=config["reference"][config["genome"]]["tss_bed"],
+        stitch_distance=config["stitch_distance"],
+        tss_distance=config["tss_distance"],
+        peak_caller_type="{peak_caller_type}",
+        treatment_sample="{treatment_sample}",
+        control_sample=lambda w: TREATMENT_SAMPLE_TO_CONTROL_SAMPLE.get(w.treatment_sample, "nocontrol"),
+        dupstatus="{dupstatus}",
+        file_base=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}"),
+        control_flag=config["macs2_control"],
+        rose_root="/opt/ROSE",
+        rose_python="/opt/conda/envs/rose/bin/python",
+        prep_bed_name="rose_input.prepared.stitched.bed",
+        prep_gff_name="rose_input.prepared.stitched.gff",
+    output:
+        enh_with_super=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Enhancers_withSuper.bed"),
+        gateway_enhancers=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Gateway_Enhancers.bed"),
+        gateway_super=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_Gateway_SuperEnhancers.bed"),
+        all_enhancers_table=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_AllEnhancers.table.txt"),
+        super_enhancer_to_gene=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt"),
+        super_gene_to_enhancer=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt"),
+    container: config["containers"].get("rose", "docker://nciccbr/ccbr_rose:v1")
+    shell:
+        """
+        set -euo pipefail
+        if [[ "{params.genome}" != "hg19" && "{params.genome}" != "hg38" && "{params.genome}" != "mm10" ]]; then
+            echo "ERROR: rule rose_treatment supports only hg19, hg38, and mm10. Found genome={params.genome}"
+            exit 1
+        fi
+
+        mkdir -p {params.file_base}
+
+        treat_bam={input.treatment_bam}
+        control="{params.control_sample}"
+        cntrl_bam={input.control_bam}
+
+        control_arg=""
+        if [[ "$control" == "nocontrol" ]]; then
+            echo "ROSE control BAM omitted for control-free run (nocontrol sentinel)"
+        elif [[ "{params.control_flag}" == "N" ]] && [[ "{params.peak_caller_type}" == "macs2_narrow" || "{params.peak_caller_type}" == "macs2_broad" ]]; then
+            echo "ROSE control BAM omitted for MACS2 with macs2_control=N"
+        else
+            control_arg="--control-bam ${{cntrl_bam}}"
+        fi
+
+        run-prep-rose \
+            --peak-file {input.peak_file} \
+            --peak-format auto \
+            --sample-id {params.treatment_sample} \
+            --treatment-bam ${{treat_bam}} \
+            ${{control_arg}} \
+            --genome {params.genome} \
+            --tss-bed {params.tss_bed} \
+            --stitch-distance {params.stitch_distance} \
+            --tss-distance {params.tss_distance} \
+            --output-dir {params.file_base} \
+            --prepared-bed-name {params.prep_bed_name} \
+            --prepared-gff-name {params.prep_gff_name}
+
+        prep_bed={params.file_base}/{params.prep_bed_name}
+        prep_gff={params.file_base}/{params.prep_gff_name}
+        num_of_peaks=`cat "$prep_bed" | wc -l`
+        if [[ ${{num_of_peaks}} -gt 4 ]]; then
+            cd {params.rose_root}
+            if [[ -n "$control_arg" ]]; then
+                {params.rose_python} ROSE_main.py \
+                    -g {params.genome} \
+                    -i "$prep_gff" \
+                    -r ${{treat_bam}} \
+                    -c ${{cntrl_bam}} \
+                    -s {params.stitch_distance} \
+                    -t {params.tss_distance} \
+                    -o {params.file_base}
+            else
+                {params.rose_python} ROSE_main.py \
+                    -g {params.genome} \
+                    -i "$prep_gff" \
+                    -r ${{treat_bam}} \
+                    -s {params.stitch_distance} \
+                    -t {params.tss_distance} \
+                    -o {params.file_base}
+            fi
+
+            super_table={params.file_base}/rose_input_SuperEnhancers.table.txt
+            normalized_super_table={params.file_base}/rose_input_SuperEnhancers.for_geneMapper.table.txt
+            {params.rose_python} - "$super_table" "$normalized_super_table" <<'PY'
+import sys
+
+
+def normalize_super_table(in_path, out_path):
+    newline = chr(10)
+    tab = chr(9)
+
+    with open(in_path, "r") as fin:
+        lines = fin.readlines()
+
+    header_idx = None
+    for i, raw in enumerate(lines):
+        if not raw.startswith("#") and raw.strip():
+            header_idx = i
+            break
+
+    if header_idx is None:
+        raise ValueError("No tabular header found in super-enhancer table")
+
+    header = lines[header_idx].rstrip(newline).split(tab)
+    has_rank = "enhancerRank" in header
+    has_super = "isSuper" in header
+
+    with open(out_path, "w") as fout:
+        for i in range(header_idx):
+            fout.write(lines[i])
+        out_header = list(header)
+        if not has_rank:
+            out_header.append("enhancerRank")
+        if not has_super:
+            out_header.append("isSuper")
+        fout.write(tab.join(out_header) + newline)
+
+        rank_counter = 0
+        for raw in lines[header_idx + 1 :]:
+            stripped = raw.rstrip(newline)
+            if not stripped:
+                continue
+            if stripped.startswith("#"):
+                fout.write(raw)
+                continue
+            rank_counter += 1
+            cols = stripped.split(tab)
+            if not has_rank:
+                cols.append(str(rank_counter))
+            if not has_super:
+                cols.append("1")
+            fout.write(tab.join(cols) + newline)
+
+
+normalize_super_table(sys.argv[1], sys.argv[2])
+PY
+
+            {params.rose_python} ROSE_geneMapper.py \
+                -g {params.genome} \
+                -i "$normalized_super_table" \
+                -o {params.file_base}
+
+            rm -f "$normalized_super_table"
+            rm -rf {params.file_base}/gff {params.file_base}/mappedGFF
+        else
+            echo "Less than 5 usable peaks detected (N=${{num_of_peaks}})" > {output.enh_with_super}
+            echo "Less than 5 usable peaks detected (N=${{num_of_peaks}})" > {output.gateway_enhancers}
+            echo "Less than 5 usable peaks detected (N=${{num_of_peaks}})" > {output.gateway_super}
+            echo "Less than 5 usable peaks detected (N=${{num_of_peaks}})" > {output.all_enhancers_table}
+            echo "Less than 5 usable peaks detected (N=${{num_of_peaks}})" > {output.super_enhancer_to_gene}
+            echo "Less than 5 usable peaks detected (N=${{num_of_peaks}})" > {output.super_gene_to_enhancer}
+        fi
+        """
+
+
+def get_treatment_rose_super_to_gene_inputs(wildcards):
+    files = []
+    if ("macs2_narrow" in PEAKTYPE) or ("macs2_broad" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M,s_dist=S_DISTANCE))
+    if ("gopeaks_narrow" in PEAKTYPE) or ("gopeaks_broad" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G,s_dist=S_DISTANCE))
+    if ("seacr_stringent" in PEAKTYPE) or ("seacr_relaxed" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_ENHANCER_TO_GENE.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S,s_dist=S_DISTANCE))
+    return files
+
+
+def get_treatment_rose_gene_to_enh_inputs(wildcards):
+    files = []
+    if ("macs2_narrow" in PEAKTYPE) or ("macs2_broad" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","macs2","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_M,s_dist=S_DISTANCE))
+    if ("gopeaks_narrow" in PEAKTYPE) or ("gopeaks_broad" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","gopeaks","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_G,s_dist=S_DISTANCE))
+    if ("seacr_stringent" in PEAKTYPE) or ("seacr_relaxed" in PEAKTYPE):
+        files.extend(expand(join(RESULTSDIR,"peaks","{qthresholds}","seacr","annotation","rose_treatment","{control_mode}","{treatment_sample}.{dupstatus}.{peak_caller_type}.{s_dist}","rose_input_SuperEnhancers_GENE_TO_ENHANCER.txt"),qthresholds=QTRESHOLDS,control_mode=CONTROL_MODES,treatment_sample=TREATMENT_SAMPLES,dupstatus=DUPSTATUS,peak_caller_type=PEAKTYPE_S,s_dist=S_DISTANCE))
+    return files
+
+
+rule aggregate_rose_treatment:
+    """
+    Create global treatment-level ROSE enhancer<->gene mapping files.
+    """
+    input:
+        enhancer_to_gene=get_treatment_rose_super_to_gene_inputs,
+        gene_to_enhancer=get_treatment_rose_gene_to_enh_inputs,
+    output:
+        enhancer_to_gene=join(RESULTSDIR,"peaks","annotation","rose_treatment","all_treatments.SuperEnhancers_ENHANCER_TO_GENE.txt"),
+        gene_to_enhancer=join(RESULTSDIR,"peaks","annotation","rose_treatment","all_treatments.SuperEnhancers_GENE_TO_ENHANCER.txt"),
+    params:
+        script=join(SCRIPTSDIR, "_aggregate_rose_treatment.py"),
+        enhancer_args=lambda w, input: " ".join([
+            "--enhancer-to-gene "
+            + f"{p.split(os.sep)[-6]}|{p.split(os.sep)[-3]}|{p.split(os.sep)[-2].split('.', 1)[0]}|{p.split(os.sep)[-2].split('.', 1)[1].split('.', 1)[0]}|{p.split(os.sep)[-2].rsplit('.', 1)[-1]}"
+            + "::"
+            + p
+            for p in input.enhancer_to_gene
+        ]),
+        gene_args=lambda w, input: " ".join([
+            "--gene-to-enhancer "
+            + f"{p.split(os.sep)[-6]}|{p.split(os.sep)[-3]}|{p.split(os.sep)[-2].split('.', 1)[0]}|{p.split(os.sep)[-2].split('.', 1)[1].split('.', 1)[0]}|{p.split(os.sep)[-2].rsplit('.', 1)[-1]}"
+            + "::"
+            + p
+            for p in input.gene_to_enhancer
+        ])
+    envmodules:
+        TOOLS["python3"]
+    shell:
+        """
+        set -euo pipefail
+        mkdir -p "$(dirname "{output.enhancer_to_gene}")"
+        python {params.script} \
+            {params.enhancer_args} \
+            {params.gene_args} \
+            --output-enhancer-to-gene {output.enhancer_to_gene} \
+            --output-gene-to-enhancer {output.gene_to_enhancer}
+        """
+
 if config["run_go_enrichment"]:
     rule go_enrichment_peaks:
         """
diff --git a/workflow/rules/init.smk b/workflow/rules/init.smk
index 19452d12..d6be2447 100644
--- a/workflow/rules/init.smk
+++ b/workflow/rules/init.smk
@@ -282,6 +282,30 @@ if RUN_WITHOUT_CONTROLS:
 else:
     TREATMENT_LIST_SG=TREATMENT_CONTROL_LIST
 
+# Treatment-sample helper structures used by treatment-level merged-peak analyses.
+TREATMENT_SAMPLES = []
+for t in TREATMENTS:
+    t_sample = re.sub(r'_[0-9]+$', '', t)
+    if t_sample not in TREATMENT_SAMPLES:
+        TREATMENT_SAMPLES.append(t_sample)
+
+TREATMENT_SAMPLE_TO_REPLICATES = {}
+for t_sample in TREATMENT_SAMPLES:
+    TREATMENT_SAMPLE_TO_REPLICATES[t_sample] = [
+        r for r in TREATMENTS if re.sub(r'_[0-9]+$', '', r) == t_sample
+    ]
+
+TREATMENT_SAMPLE_TO_CONTROL_SAMPLE = {}
+for t in TREATMENTS:
+    t_sample = re.sub(r'_[0-9]+$', '', t)
+    c_rep = TREAT_to_CONTRL_DICT.get(t, "nocontrol")
+    if c_rep == "nocontrol":
+        c_sample = "nocontrol"
+    else:
+        c_sample = re.sub(r'_[0-9]+$', '', c_rep)
+    if t_sample not in TREATMENT_SAMPLE_TO_CONTROL_SAMPLE:
+        TREATMENT_SAMPLE_TO_CONTROL_SAMPLE[t_sample] = c_sample
+
 # create duplication and peaktype list
 DUPSTATUS=config["dupstatus"]
 PEAKTYPE=config["peaktype"]
diff --git a/workflow/rules/peakcalls.smk b/workflow/rules/peakcalls.smk
index 36423773..dfcbba88 100644
--- a/workflow/rules/peakcalls.smk
+++ b/workflow/rules/peakcalls.smk
@@ -11,6 +11,7 @@
 localrules: count_peaks
 
 import re
+import os
 
 rule merge_control_bams:
     """
@@ -43,6 +44,37 @@ rule merge_control_bams:
         samtools index {output.merged_bam}
         """
 
+
+rule merge_treatment_bams:
+    """
+    Merge all replicates for a treatment sample into one BAM for treatment-level ROSE.
+    """
+    input:
+        bams=lambda w: expand(
+            join(RESULTSDIR, "bam", "{replicate}.{dupstatus}.bam"),
+            replicate=TREATMENT_SAMPLE_TO_REPLICATES.get(w.treatment_sample, []),
+            dupstatus=w.dupstatus,
+        )
+    output:
+        merged_bam=join(RESULTSDIR, "bam", "treatment_merged", "{treatment_sample}.{dupstatus}.merged.bam"),
+        merged_bai=join(RESULTSDIR, "bam", "treatment_merged", "{treatment_sample}.{dupstatus}.merged.bam.bai"),
+    params:
+        bam_list=lambda w, input: " ".join(input.bams)
+    threads: getthreads("merge_control_bams")
+    envmodules:
+        TOOLS["samtools"]
+    shell:
+        """
+        set -exo pipefail
+        mkdir -p "$(dirname "{output.merged_bam}")"
+        if [[ {threads} -gt 1 ]]; then
+            samtools merge -@ {threads} {output.merged_bam} {params.bam_list}
+        else
+            samtools merge {output.merged_bam} {params.bam_list}
+        fi
+        samtools index {output.merged_bam}
+        """
+
 rule create_pooled_control_fragments:
     """
     Create fragments.bed file from pooled control BAM
@@ -192,6 +224,95 @@ def get_all_peak_files(wildcards):
 
     return files
 
+
+def _caller_from_peak_type(peak_type):
+    if peak_type.startswith("macs2_"):
+        return "macs2"
+    if peak_type.startswith("gopeaks_"):
+        return "gopeaks"
+    if peak_type.startswith("seacr_"):
+        return "seacr"
+    raise ValueError(f"Unsupported peak_type: {peak_type}")
+
+
+def _suffix_from_peak_type(peak_type):
+    return peak_type.split("_", 1)[1]
+
+
+def get_treatment_replicate_peak_files(wildcards):
+    caller = _caller_from_peak_type(wildcards.peak_type)
+    suffix = _suffix_from_peak_type(wildcards.peak_type)
+
+    if caller == "macs2":
+        base_list = TREATMENT_CONTROL_LIST_POOLED if wildcards.control_mode == "pooled" else TREATMENT_LIST_M
+    else:
+        base_list = TREATMENT_CONTROL_LIST_POOLED if wildcards.control_mode == "pooled" else TREATMENT_LIST_SG
+
+    tc_list = _filter_tc_list(base_list, wildcards.control_mode)
+    treatment_reps = set(TREATMENT_SAMPLE_TO_REPLICATES.get(wildcards.treatment_sample, []))
+
+    peaks = []
+    for tc_pair in tc_list:
+        treatment_rep = tc_pair.split("_vs_", 1)[0]
+        if treatment_rep not in treatment_reps:
+            continue
+        peaks.append(
+            join(
+                RESULTSDIR,
+                "peaks",
+                wildcards.qthresholds,
+                caller,
+                "peak_output",
+                wildcards.control_mode,
+                f"{tc_pair}.{wildcards.dupstatus}.{suffix}.peaks.bed",
+            )
+        )
+    return sorted(set(peaks))
+
+
+rule merge_treatment_peaks:
+    """
+    Merge replicate peak sets into treatment-level consensus peaks without re-calling.
+    Output columns: chrom, start, end, unique_peak_id, replicate_ids, source_peak_count, width_guard_flag
+    """
+    input:
+        peaks=get_treatment_replicate_peak_files
+    output:
+        merged=join(RESULTSDIR,"peaks","{qthresholds}","{peak_caller}","peak_output","{control_mode}","treatment_merged","{treatment_sample}.{dupstatus}.{peak_type}.peaks.bed")
+    params:
+        merge_script=join(SCRIPTSDIR, "_merge_treatment_peaks.py"),
+        input_args=lambda w, input: " ".join(
+            [
+                "--input "
+                + os.path.basename(p).split("_vs_", 1)[0]
+                + "::"
+                + p
+                for p in input.peaks
+            ]
+        ),
+        overlap_bp_min=1,
+        min_replicate_support=2,
+        max_merged_width_bp=10000,
+        strict_overlap_bp=50,
+    envmodules:
+        TOOLS["python3"]
+    shell:
+        """
+        set -euo pipefail
+        mkdir -p "$(dirname "{output.merged}")"
+        if [[ -z "{params.input_args}" ]]; then
+            : > {output.merged}
+            exit 0
+        fi
+        python {params.merge_script} \
+            {params.input_args} \
+            --output {output.merged} \
+            --overlap-bp-min {params.overlap_bp_min} \
+            --min-replicate-support {params.min_replicate_support} \
+            --max-merged-width-bp {params.max_merged_width_bp} \
+            --strict-overlap-bp {params.strict_overlap_bp}
+        """
+
 rule macs2_narrow:
     '''
     MACS2 can be run with and without a control. This featured is controlled in the config.yaml file
diff --git a/workflow/scripts/_aggregate_homer_treatment.py b/workflow/scripts/_aggregate_homer_treatment.py
new file mode 100644
index 00000000..98bef31a
--- /dev/null
+++ b/workflow/scripts/_aggregate_homer_treatment.py
@@ -0,0 +1,41 @@
+#!/usr/bin/env python3
+"""Aggregate treatment-level HOMER annotation summaries across treatments."""
+
+from __future__ import annotations
+
+import argparse
+import os
+
+
+def parse_input_spec(spec: str):
+    if "::" not in spec:
+        raise ValueError(f"Invalid --input '{spec}', expected label::path")
+    label, path = spec.split("::", 1)
+    return label, path
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description="Aggregate treatment-level HOMER summaries")
+    parser.add_argument("--input", action="append", required=True, help="label::path")
+    parser.add_argument("--output", required=True)
+    args = parser.parse_args()
+
+    os.makedirs(os.path.dirname(args.output), exist_ok=True)
+    with open(args.output, "w", encoding="utf-8") as out:
+        out.write("source\tannotation\tdistance_to_tss\tnumber_of_peaks\tpercent_of_peaks\ttotal_size_bp\tlog10_p_value\tlog2_ratio\tlogP_enrichment\n")
+        for spec in args.input:
+            source, path = parse_input_spec(spec)
+            if not os.path.exists(path):
+                continue
+            with open(path, "r", encoding="utf-8") as handle:
+                for line in handle:
+                    line = line.rstrip("\n")
+                    if not line or line.startswith("#"):
+                        continue
+                    if line.startswith("Annotation\t"):
+                        continue
+                    out.write(f"{source}\t{line}\n")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/workflow/scripts/_aggregate_rose_treatment.py b/workflow/scripts/_aggregate_rose_treatment.py
new file mode 100644
index 00000000..01d5a70f
--- /dev/null
+++ b/workflow/scripts/_aggregate_rose_treatment.py
@@ -0,0 +1,60 @@
+#!/usr/bin/env python3
+"""Aggregate treatment-level ROSE outputs across treatments."""
+
+from __future__ import annotations
+
+import argparse
+import os
+
+
+def parse_input_spec(spec: str):
+    if "::" not in spec:
+        raise ValueError(f"Invalid --input '{spec}', expected label::path")
+    label, path = spec.split("::", 1)
+    return label, path
+
+
+def parse_label(label: str):
+    parts = label.split("|")
+    if len(parts) != 5:
+        raise ValueError(
+            f"Invalid aggregate label '{label}', expected peak_caller|control_mode|treatment_sample|dup_status|stitch_distance"
+        )
+    return parts
+
+
+def aggregate(input_specs, output_path):
+    os.makedirs(os.path.dirname(output_path), exist_ok=True)
+    with open(output_path, "w", encoding="utf-8") as out:
+        out.write("peak_caller\ttreatment_sample\tcontrol_mode\tdup_status\tstitch_distance\trow\n")
+        for spec in input_specs:
+            source, path = parse_input_spec(spec)
+            peak_caller, control_mode, treatment_sample, dup_status, stitch_distance = parse_label(source)
+            if not os.path.exists(path):
+                continue
+            with open(path, "r", encoding="utf-8") as handle:
+                for line in handle:
+                    line = line.rstrip("\n")
+                    if not line:
+                        continue
+                    if line.startswith("Less than 5 usable peaks detected"):
+                        continue
+                    out.write(
+                        f"{peak_caller}\t{treatment_sample}\t{control_mode}\t{dup_status}\t{stitch_distance}\t{line}\n"
+                    )
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description="Aggregate treatment-level ROSE mapping outputs")
+    parser.add_argument("--enhancer-to-gene", action="append", default=[])
+    parser.add_argument("--gene-to-enhancer", action="append", default=[])
+    parser.add_argument("--output-enhancer-to-gene", required=True)
+    parser.add_argument("--output-gene-to-enhancer", required=True)
+    args = parser.parse_args()
+
+    aggregate(args.enhancer_to_gene, args.output_enhancer_to_gene)
+    aggregate(args.gene_to_enhancer, args.output_gene_to_enhancer)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/workflow/scripts/_merge_treatment_peaks.py b/workflow/scripts/_merge_treatment_peaks.py
new file mode 100644
index 00000000..9f8577b6
--- /dev/null
+++ b/workflow/scripts/_merge_treatment_peaks.py
@@ -0,0 +1,194 @@
+#!/usr/bin/env python3
+"""Merge replicate peak BEDs into treatment-level consensus peak sets.
+
+Consensus policy:
+- overlap >= overlap_bp_min counts as support
+- keep clusters supported by at least min_replicate_support distinct replicates
+- merged interval spans leftest-left/rightest-right
+- width guard: if merged width exceeds max_merged_width_bp, re-cluster with strict_overlap_bp
+  and if still too wide keep interval but mark width_guard_flag=1
+"""
+
+from __future__ import annotations
+
+import argparse
+from collections import defaultdict
+from dataclasses import dataclass
+from typing import Dict, List, Sequence, Tuple
+
+
+@dataclass
+class Peak:
+    chrom: str
+    start: int
+    end: int
+    replicate: str
+
+
+def parse_input_spec(spec: str) -> Tuple[str, str]:
+    if "::" not in spec:
+        raise ValueError(f"Invalid --input value '{spec}'. Expected replicate::path format.")
+    replicate, path = spec.split("::", 1)
+    if not replicate:
+        raise ValueError(f"Invalid --input value '{spec}'. Missing replicate id.")
+    if not path:
+        raise ValueError(f"Invalid --input value '{spec}'. Missing path.")
+    return replicate, path
+
+
+def load_peaks(input_specs: Sequence[str]) -> Dict[str, List[Peak]]:
+    by_chrom: Dict[str, List[Peak]] = defaultdict(list)
+    for spec in input_specs:
+        replicate, path = parse_input_spec(spec)
+        with open(path, "r", encoding="utf-8") as handle:
+            for line in handle:
+                if not line.strip() or line.startswith("#"):
+                    continue
+                parts = line.rstrip("\n").split("\t")
+                if len(parts) < 3:
+                    continue
+                chrom = parts[0]
+                try:
+                    start = int(parts[1])
+                    end = int(parts[2])
+                except ValueError:
+                    continue
+                if end <= start:
+                    continue
+                by_chrom[chrom].append(Peak(chrom=chrom, start=start, end=end, replicate=replicate))
+    return by_chrom
+
+
+def cluster_intervals(peaks: Sequence[Peak], overlap_bp: int) -> List[List[Peak]]:
+    if not peaks:
+        return []
+    sorted_peaks = sorted(peaks, key=lambda p: (p.start, p.end))
+    clusters: List[List[Peak]] = []
+    current = [sorted_peaks[0]]
+    current_end = sorted_peaks[0].end
+
+    for peak in sorted_peaks[1:]:
+        # overlap condition with configurable minimum overlap length.
+        if peak.start <= current_end - overlap_bp:
+            current.append(peak)
+            if peak.end > current_end:
+                current_end = peak.end
+        else:
+            clusters.append(current)
+            current = [peak]
+            current_end = peak.end
+
+    clusters.append(current)
+    return clusters
+
+
+def cluster_to_record(
+    cluster: Sequence[Peak],
+    min_support: int,
+    max_width: int,
+    strict_overlap: int,
+) -> List[Tuple[str, int, int, str, str, int, int]]:
+    reps = sorted({p.replicate for p in cluster})
+    if len(reps) < min_support:
+        return []
+
+    start = min(p.start for p in cluster)
+    end = max(p.end for p in cluster)
+    width = end - start
+
+    if width <= max_width:
+        return [
+            (
+                cluster[0].chrom,
+                start,
+                end,
+                "",
+                ",".join(reps),
+                len(cluster),
+                0,
+            )
+        ]
+
+    # Too wide: re-cluster with stricter overlap and evaluate each subcluster.
+    refined_records: List[Tuple[str, int, int, str, str, int, int]] = []
+    refined_clusters = cluster_intervals(cluster, strict_overlap)
+    for sub in refined_clusters:
+        sub_reps = sorted({p.replicate for p in sub})
+        if len(sub_reps) < min_support:
+            continue
+        sub_start = min(p.start for p in sub)
+        sub_end = max(p.end for p in sub)
+        sub_width = sub_end - sub_start
+        refined_records.append(
+            (
+                sub[0].chrom,
+                sub_start,
+                sub_end,
+                "",
+                ",".join(sub_reps),
+                len(sub),
+                1 if sub_width > max_width else 0,
+            )
+        )
+
+    return refined_records
+
+
+def assign_peak_ids(records: Sequence[Tuple[str, int, int, str, str, int, int]]) -> List[Tuple[str, int, int, str, str, int, int]]:
+    assigned: List[Tuple[str, int, int, str, str, int, int]] = []
+    for index, rec in enumerate(records, start=1):
+        peak_id = f"merged_peak_{index:06d}"
+        assigned.append((rec[0], rec[1], rec[2], peak_id, rec[4], rec[5], rec[6]))
+    return assigned
+
+
+def write_records(records: Sequence[Tuple[str, int, int, str, str, int, int]], output_path: str) -> None:
+    with open(output_path, "w", encoding="utf-8") as handle:
+        for rec in records:
+            handle.write(
+                "\t".join(
+                    [
+                        rec[0],
+                        str(rec[1]),
+                        str(rec[2]),
+                        str(rec[3]),
+                        rec[4],
+                        str(rec[5]),
+                        str(rec[6]),
+                    ]
+                )
+                + "\n"
+            )
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description="Merge replicate peaks into treatment-level consensus sets")
+    parser.add_argument("--input", action="append", required=True, help="replicate::path to peak BED")
+    parser.add_argument("--output", required=True, help="Output BED path")
+    parser.add_argument("--overlap-bp-min", type=int, default=1)
+    parser.add_argument("--min-replicate-support", type=int, default=2)
+    parser.add_argument("--max-merged-width-bp", type=int, default=10000)
+    parser.add_argument("--strict-overlap-bp", type=int, default=50)
+    args = parser.parse_args()
+
+    by_chrom = load_peaks(args.input)
+    all_records: List[Tuple[str, int, int, str, str, int, int]] = []
+
+    for chrom in sorted(by_chrom):
+        clusters = cluster_intervals(by_chrom[chrom], args.overlap_bp_min)
+        for cluster in clusters:
+            all_records.extend(
+                cluster_to_record(
+                    cluster,
+                    min_support=args.min_replicate_support,
+                    max_width=args.max_merged_width_bp,
+                    strict_overlap=args.strict_overlap_bp,
+                )
+            )
+
+    all_records.sort(key=lambda r: (r[0], r[1], r[2]))
+    write_records(assign_peak_ids(all_records), args.output)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/workflow/scripts/_prep_rose_input.py b/workflow/scripts/_prep_rose_input.py
index 48dea01c..d1f80dbb 100644
--- a/workflow/scripts/_prep_rose_input.py
+++ b/workflow/scripts/_prep_rose_input.py
@@ -272,7 +272,8 @@ def expected_rose_outputs(output_dir, sample_id):
         os.path.join(output_dir, "%s_AllStitched.table.txt" % prefix),
         os.path.join(output_dir, "%s_AllEnhancers.table.txt" % prefix),
         os.path.join(output_dir, "%s_SuperEnhancers.table.txt" % prefix),
-        os.path.join(output_dir, "%s_Enhancer_to_Gene.table.txt" % prefix),
+        os.path.join(output_dir, "%s_SuperEnhancers_ENHANCER_TO_GENE.txt" % prefix),
+        os.path.join(output_dir, "%s_SuperEnhancers_GENE_TO_ENHANCER.txt" % prefix),
     ]