diff --git a/CHANGELOG.md b/CHANGELOG.md
index 401b6f4..7ec9444 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,7 @@
## ASPEN development version
+- Add the missing `deeptools_plot` rule to generate the final `plotHeatmap`/`plotProfile` TSS and metagene PDFs from the `deeptools_mat` matrices, and wire these (plus the upstream `deeptools_prep`/`deeptools_mat` outputs) into `rule all` so `dryrun`/`run` actually schedule them; previously nothing requested these outputs so Snakemake silently never ran deepTools past matrix generation. Also fix `deeptools_mat` to use the workdir's configured genome's TSS BED archive (`config[GENOME]["tssBed"]`) instead of an always-hardcoded `hg38_tssbeds.tar.gz`, and drop the unused `{bedtype}` wildcard. Update `docs/outputs.md` to describe the actual output filenames/location (`results/deeptools/{group}.{bamtype}.{TSS,metagene}.{heatmap,profile}.pdf`) instead of the previously documented, never-implemented `protein_coding`/`all_samples` naming scheme. (#7, @kopardev)
+- Improve failed-run diagnostics by recording the failing Snakemake rule and command in `pipeline.status.json` and the human-readable `pipeline.failed` marker, so users can immediately see which step (for example, `deeptools_mat`) triggered the failure without digging through `snakemake.log` manually. (#000, @kopardev)
- Document successful `dryrun` and `run` console output in `docs/deployment.md`, including example screenshots, an explanation of `STEP`/`OK`/`INFO`/`NEXT`, and explicit links between wrapper messages and `WORKDIR` status files. (#147, @kopardev)
- Add attempt-aware resource scaling for the main Snakemake rules, document how baseline `cluster.json` resources interact with retry-based `resources.gres` scaling, and fix dry-run failure handling so failed dry-runs now exit non-zero instead of reporting success. (#152, @kopardev)
- Detect SLURM-killed Snakemake child jobs in classic cluster mode by adding a `--cluster-status` hook and `scancel` integration, so timed out, OOM-killed, cancelled, or otherwise failed child jobs now retry/fail cleanly instead of leaving the master workflow hanging indefinitely. (#148, @kopardev)
diff --git a/aspen b/aspen
index e64d38a..2a89cc3 100755
--- a/aspen
+++ b/aspen
@@ -408,7 +408,11 @@ function dryrun() {
echo "Git Commit/Tag: $GIT_COMMIT_TAG" > ${WORKDIR}/dryrun_git_commit.txt
log_ok "Dry-run completed."
log_info "Full dry-run output captured in ${WORKDIR}/dryrun.log"
- log_next "Submit run with: aspen --workdir=$WORKDIR --runmode=run"
+ if grep -qi "Nothing to be done" "${WORKDIR}/dryrun.log"; then
+ log_info "Dry-run reports no pending jobs for this workdir; skipping NEXT prompt."
+ else
+ log_next "Submit run with: aspen --workdir=$WORKDIR --runmode=run"
+ fi
}
##########################################################################################
diff --git a/config/config.yaml b/config/config.yaml
index 2b474ef..82ee3bc 100644
--- a/config/config.yaml
+++ b/config/config.yaml
@@ -205,6 +205,17 @@ jaccard_min_peaks: 1000
# fastq screen configuration file
fastqscreen_config: "$PIPELINE_HOME/config/fastq_screen_config.txt"
+# deepTools matrix/heatmap parameters
+# these control the TSS-centered and metagene windows used by computeMatrix
+# and can be overridden per project without editing the Snakemake rules
+# Default values match the issue #7 design and are expressed in base pairs.
+deeptools:
+ tss_upstream: 3000
+ tss_downstream: 3000
+ metagene_body_length: 2000
+ metagene_upstream: 2000
+ metagene_downstream: 2000
+
# Containers
masterdocker: "docker://nciccbr/ccbr_atacseq:v13-feat"
cutadaptdocker: "docker://nciccbr/ccbr_atacseq:v11-feat"
diff --git a/docs/css/custom.css b/docs/css/custom.css
index c4c8782..180dd2e 100644
--- a/docs/css/custom.css
+++ b/docs/css/custom.css
@@ -44,3 +44,44 @@ div.admonition.warning {
margin: 10px 0;
border-radius: 5px;
}
+
+/* Increase content width for large monitors */
+@media (min-width: 1400px) {
+ .wy-nav-content-wrap {
+ margin-left: 300px;
+ }
+ .wy-nav-content {
+ max-width: 1100px;
+ }
+}
+
+@media (min-width: 1600px) {
+ .wy-nav-content-wrap {
+ margin-left: 300px;
+ }
+ .wy-nav-content {
+ max-width: 1300px;
+ }
+}
+
+@media (min-width: 1920px) {
+ .wy-nav-content-wrap {
+ margin-left: 300px;
+ }
+ .wy-nav-content {
+ max-width: 1500px;
+ }
+}
+
+/* Keep header and body columns in one shared table layout. */
+.rst-content table {
+ border-collapse: collapse;
+ width: 100%;
+ table-layout: auto;
+}
+
+.rst-content table td,
+.rst-content table th {
+ overflow-wrap: break-word;
+ padding: 8px 12px;
+}
diff --git a/docs/outputs.md b/docs/outputs.md
index 82af38c..f5fbb2f 100644
--- a/docs/outputs.md
+++ b/docs/outputs.md
@@ -143,6 +143,49 @@ Content details:
| spikein | 1 folder per sample | - Per sample spike-in counts.
- Overall scaling factors table. |
| tmp | various | - Can be deleted.
- Blacklist index.
- Intermediate FASTQs.
- Genrich output reads. |
+### deepTools TSS and metagene outputs
+
+ASPEN generates deepTools heatmaps and profile plots to summarize ATAC-seq accessibility around transcription start sites (TSS) and across gene bodies. These exploratory visualizations reveal whether signal concentrates near promoters, distributes across gene bodies, or shows group-specific patterns.
+
+#### Purpose and when generated
+
+These outputs summarize enrichment patterns across many genomic regions, generated after bigWig preparation as part of the pipeline's final outputs (`rule all`). They are derived from normalized bigWig coverage to show accessibility patterns clearly and enable visual quality assessment alongside numerical QC metrics.
+
+#### Output files
+
+| Output file | Region type | Visualization | Use case |
+| --- | --- | --- | --- |
+| `{group}.{bamtype}.TSS.heatmap.pdf` | ±3000 bp around TSS | Heatmap (regions × position) | View each gene's promoter signal strength and distribution |
+| `{group}.{bamtype}.TSS.profile.pdf` | ±3000 bp around TSS | Aggregate line plot (mean ± SE) | Compare average promoter accessibility across samples |
+| `{group}.{bamtype}.metagene.heatmap.pdf` | Gene body ± flanks | Heatmap (regions × normalized length) | Assess whether signal spreads across genes or clusters at promoters |
+| `{group}.{bamtype}.metagene.profile.pdf` | Gene body ± flanks | Aggregate line plot | Understand average accessibility distribution within genes |
+
+Where `{group}` = sample name from manifest (groups can span multiple replicates) and `{bamtype}` = `dedup` or `nondedup`.
+
+#### Data source and normalization
+
+- **Input**: Pooled bigWig coverage per group per BAM type from `results/visualization/{bamtype}/reads_bigwig/`
+- **Normalization**: Spike-in-scaled (or 1.0× if disabled)
+- **Annotation**: Ensembl protein-coding genes (genome-specific BED archive)
+- **Intermediate files**: Saved under `results/deeptools/temp/` for re-runs and debugging
+
+#### Interpretation
+
+- **TSS heatmaps**: Each row = one gene. Red/hot colors near position 0 = strong promoter accessibility.
+- **TSS profiles**: Peak at position 0 = promoter-centered open chromatin (expected for active genes). Flat or shifted peak = different accessibility pattern.
+- **Metagene heatmaps**: Red near left edge = signal at promoters. Uniform red across = open throughout gene bodies. Red at right edge = 3′ bias.
+- **Metagene profiles**: Shape indicates accessibility distribution: peaked = localized; flat = distributed; sloped = progressive.
+
+#### Configuration
+
+Tune window sizes via `config.deeptools.*` in `config.yaml`:
+
+- `tss_upstream/downstream`: Flanking distance around TSS (default: 3000 bp)
+- `metagene_body_length`: Normalized gene body length (default: 2000 bp)
+- `metagene_upstream/downstream`: Flanking distance around gene start/end (default: 2000 bp)
+
+Edit these before running to adjust the window sizes for your analysis needs.
+
!!! note
BAM files from `dedupBam` can be used for downstream footprinting analysis using [CCBR_TOBIAS](https://github.com/CCBR/CCBR_Tobias) pipeline
@@ -532,4 +575,5 @@ This directory contains all .err and .out log files generated by SLURM for jobs
This structure is particularly useful for troubleshooting and debugging, especially when the SLURM job IDs of failed jobs are known. By examining the corresponding .err or .out files, users can efficiently identify the source of errors within specific Snakemake rules and wildcards.
-> DISCLAIMER: This folder hierarchy is significantly different than v1.2.0 and is subject to change with subsequent versions.
+!!! warning "DISCLAIMER"
+ This folder hierarchy is significantly different than v1.2.0 and is subject to change with subsequent versions.
diff --git a/resources/cluster.json b/resources/cluster.json
index 34f6344..259b3d3 100755
--- a/resources/cluster.json
+++ b/resources/cluster.json
@@ -30,6 +30,21 @@
"threads": "16",
"time": "12:00:00"
},
+ "deeptools_prep": {
+ "mem": "32g",
+ "threads": "2",
+ "time": "4:00:00"
+ },
+ "deeptools_mat": {
+ "mem": "256g",
+ "threads": "32",
+ "time": "12:00:00"
+ },
+ "deeptools_plot": {
+ "mem": "16g",
+ "threads": "2",
+ "time": "2:00:00"
+ },
"align": {
"mem": "240g",
"threads": "16",
diff --git a/workflow/Snakefile b/workflow/Snakefile
index 3972f74..1ef5431 100644
--- a/workflow/Snakefile
+++ b/workflow/Snakefile
@@ -71,6 +71,11 @@ rule all:
expand(join(PEAKSDIR, "{peakcaller}", "DiffATAC", "{bamtype}", "{method}", "degs.done"), method=COUNTING_METHODS, peakcaller=PEAKCALLERS, bamtype=BAM_TYPES) if CONTRASTS.shape[0] > 0 else [],
# aggregated diffatac results
expand(join(PEAKSDIR, "{peakcaller}", "DiffATAC", "{bamtype}", "{method}", "all_diff_atacs.tsv"), method=COUNTING_METHODS, peakcaller=PEAKCALLERS, bamtype=BAM_TYPES) if CONTRASTS.shape[0] > 0 else [],
+ # deepTools TSS/metagene heatmaps and profile plots (issue #7)
+ expand(join(RESULTSDIR, "deeptools", "{group}.{bamtype}.TSS.heatmap.pdf"), group=SAMPLES, bamtype=BAM_TYPES),
+ expand(join(RESULTSDIR, "deeptools", "{group}.{bamtype}.TSS.profile.pdf"), group=SAMPLES, bamtype=BAM_TYPES),
+ expand(join(RESULTSDIR, "deeptools", "{group}.{bamtype}.metagene.heatmap.pdf"), group=SAMPLES, bamtype=BAM_TYPES),
+ expand(join(RESULTSDIR, "deeptools", "{group}.{bamtype}.metagene.profile.pdf"), group=SAMPLES, bamtype=BAM_TYPES),
# NOTE (CCBR/ASPEN#144): do NOT rely on `command -v jobby`/`command -v spooker`
# to decide whether the pinned ccbr_tools install is needed -- SLURM jobs
diff --git a/workflow/rules/align.smk b/workflow/rules/align.smk
index c8a7b2c..74c88ce 100644
--- a/workflow/rules/align.smk
+++ b/workflow/rules/align.smk
@@ -481,3 +481,144 @@ bamCoverage \
--scaleFactor $sf \
--numberOfProcessors {threads}
"""
+
+
+#########################################################
+
+rule deeptools_prep:
+# """
+# Prepare per-group deepTools inputs for TSS and metagene plots.
+# Input: bigWig files for a group or pooled sample set.
+# Output: a manifest with sample labels on the first line and paths on the second line.
+# This matches the format consumed by computeMatrix in the deepTools workflow.
+# """
+ input:
+ bw=lambda wildcards: expand(
+ join(RESULTSDIR, "visualization", "{bamtype}", "reads_bigwig", "{replicate}.reads.bw"),
+ bamtype=wildcards.bamtype,
+ replicate=SAMPLE2REPLICATES[wildcards.group],
+ )
+ output:
+ prep=join(RESULTSDIR, "deeptools", "temp", "{group}.{bamtype}.reads.deeptools_prep"),
+ params:
+ group="{group}",
+ bamtype="{bamtype}",
+ container: config["deeptoolsdocker"]
+ threads: 1
+ shell:"""
+set -euo pipefail
+unset PYTHONPATH
+labels=$(printf '%s\n' {input.bw} | xargs -n1 basename | sed 's/\.reads\.bw$//' | paste -sd' ' -)
+bw_list=$(printf '%s\n' {input.bw} | paste -sd' ' -)
+printf '%s\n%s\n' "$labels" "$bw_list" > {output.prep}
+"""
+
+
+rule deeptools_mat:
+# """
+# Compute deepTools matrices for TSS and metagene summaries.
+# Input: the prepared bigWig manifest plus the genome's bundled TSS BED archive.
+# Output: reference-point and scale-region matrices for each group.
+# """
+ input:
+ prep=join(RESULTSDIR, "deeptools", "temp", "{group}.{bamtype}.reads.deeptools_prep"),
+ bedtar=config[GENOME]["tssBed"],
+ output:
+ tss=join(RESULTSDIR, "deeptools", "temp", "{group}.{bamtype}.TSS.mat.gz"),
+ metagene=join(RESULTSDIR, "deeptools", "temp", "{group}.{bamtype}.metagene.mat.gz"),
+ params:
+ group="{group}",
+ bamtype="{bamtype}",
+ tss_upstream=config["deeptools"]["tss_upstream"],
+ tss_downstream=config["deeptools"]["tss_downstream"],
+ metagene_body_length=config["deeptools"]["metagene_body_length"],
+ metagene_upstream=config["deeptools"]["metagene_upstream"],
+ metagene_downstream=config["deeptools"]["metagene_downstream"],
+ container: config["deeptoolsdocker"]
+ threads: getthreads("deeptools_mat")
+ resources:
+ gres=lambda wildcards, attempt: scale_gres("deeptools_mat", attempt)
+ shell:"""
+set -euo pipefail
+unset PYTHONPATH
+
+tmpdir=$(mktemp -d)
+tar -xzf {input.bedtar} -C "$tmpdir"
+cat "$tmpdir"/*.bed > "$tmpdir"/combined.bed
+
+bw_list=$(awk 'NR==2 {{print $0}}' {input.prep})
+labels=$(awk 'NR==1 {{print $0}}' {input.prep})
+
+computeMatrix reference-point \
+ -S $bw_list \
+ -R "$tmpdir/combined.bed" \
+ -p {threads} \
+ --referencePoint TSS \
+ --upstream {params.tss_upstream} \
+ --downstream {params.tss_downstream} \
+ --skipZeros \
+ --samplesLabel $labels \
+ -o {output.tss}
+
+computeMatrix scale-regions \
+ -S $bw_list \
+ -R "$tmpdir/combined.bed" \
+ -p {threads} \
+ --regionBodyLength {params.metagene_body_length} \
+ --upstream {params.metagene_upstream} \
+ --downstream {params.metagene_downstream} \
+ --skipZeros \
+ --samplesLabel $labels \
+ -o {output.metagene}
+"""
+
+
+#########################################################
+
+rule deeptools_plot:
+# """
+# Render final deepTools heatmaps and profile plots from the TSS and metagene matrices.
+# Input: the TSS and metagene matrices produced by deeptools_mat.
+# Output: per-group/bamtype heatmap and profile PDFs (the final, user-facing deepTools products).
+# """
+ input:
+ tss=rules.deeptools_mat.output.tss,
+ metagene=rules.deeptools_mat.output.metagene,
+ output:
+ tss_heatmap=join(RESULTSDIR, "deeptools", "{group}.{bamtype}.TSS.heatmap.pdf"),
+ tss_profile=join(RESULTSDIR, "deeptools", "{group}.{bamtype}.TSS.profile.pdf"),
+ metagene_heatmap=join(RESULTSDIR, "deeptools", "{group}.{bamtype}.metagene.heatmap.pdf"),
+ metagene_profile=join(RESULTSDIR, "deeptools", "{group}.{bamtype}.metagene.profile.pdf"),
+ params:
+ group="{group}",
+ bamtype="{bamtype}",
+ container: config["deeptoolsdocker"]
+ threads: getthreads("deeptools_plot")
+ resources:
+ gres=lambda wildcards, attempt: scale_gres("deeptools_plot", attempt)
+ shell:"""
+set -euo pipefail
+unset PYTHONPATH
+
+plotHeatmap \
+ -m {input.tss} \
+ -out {output.tss_heatmap} \
+ --refPointLabel TSS
+
+plotProfile \
+ -m {input.tss} \
+ -out {output.tss_profile} \
+ --refPointLabel TSS
+
+plotHeatmap \
+ -m {input.metagene} \
+ -out {output.metagene_heatmap} \
+ --startLabel "start" \
+ --endLabel "end"
+
+plotProfile \
+ -m {input.metagene} \
+ -out {output.metagene_profile} \
+ --startLabel "start" \
+ --endLabel "end"
+"""
diff --git a/workflow/scripts/ccbr_pipeline_logging.sh b/workflow/scripts/ccbr_pipeline_logging.sh
index f6ad7a3..70164bd 100644
--- a/workflow/scripts/ccbr_pipeline_logging.sh
+++ b/workflow/scripts/ccbr_pipeline_logging.sh
@@ -90,9 +90,12 @@ function get_git_commitid_tag() {
# place. Positional args:
# 1=state 2=reason 3=job_id 4=host 5=submission_ts 6=start_ts
# 7=duration_seconds 8=exit_code 9=tasks_done 10=tasks_total
+# 11=failed_rule 12=failed_command
function _pipeline_write_status_json() {
local state="$1" reason="$2" job_id="$3" host="$4" submission_ts="$5"
local start_ts="$6" duration_seconds="$7" exit_code="$8" tasks_done="$9" tasks_total="${10}"
+ local failed_rule="${11:-}"
+ local failed_command="${12:-}"
local now_utc sidecar sidecar_tmp start_ts_json
now_utc=$(date -u +"%Y-%m-%dT%H:%M:%SZ")
@@ -119,6 +122,8 @@ function _pipeline_write_status_json() {
[[ -n "${job_id}" ]] && printf 'Job ID : %s\n' "${job_id}"
[[ -n "${exit_code}" ]] && printf 'Exit code: %s\n' "${exit_code}"
[[ -n "${tasks_done}" && -n "${tasks_total}" ]] && printf 'Progress : %s / %s steps done\n' "${tasks_done}" "${tasks_total}"
+ [[ -n "${failed_rule}" ]] && printf 'Rule : %s\n' "${failed_rule}"
+ [[ -n "${failed_command}" ]] && printf 'Command : %s\n' "${failed_command}"
printf 'Log : %s\n' "${SNAKEMAKE_LOG_PATH}"
printf 'Updated : %s\n' "${now_utc}"
} > "${WORKDIR}/pipeline.${state}"
@@ -140,6 +145,8 @@ function _pipeline_write_status_json() {
"exit_code": ${exit_code:-null},
"tasks_done": ${tasks_done:-null},
"tasks_total": ${tasks_total:-null},
+ "failed_rule": $(json_escape_string_or_null "${failed_rule}"),
+ "failed_command": $(json_escape_string_or_null "${failed_command}"),
"snakemake_log": "$(json_escape "${SNAKEMAKE_LOG_PATH}")",
"timestamp_utc": "${now_utc}"
}
@@ -177,13 +184,22 @@ function write_pipeline_state_marker() {
# script (after `module load ccbrpipeliner; source ccbr_pipeline_logging.sh`).
# Requires the caller to export _START_EPOCH/_START_TS right after module load.
# Args: state reason [job_id] [exit_code]
+function json_escape_string_or_null() {
+ local value="${1:-}"
+ if [[ -z "${value}" ]]; then
+ printf 'null'
+ else
+ printf '"%s"' "$(json_escape "${value}")"
+ fi
+}
+
function write_pipeline_state_marker_job() {
local state="$1"
local reason="${2:-}"
local job_id="${3:-${SLURM_JOB_ID:-NA}}"
local exit_code="${4:-}"
local host submission_ts="" submit_epoch duration_seconds=""
- local tasks_done="" tasks_total="" raw_tasks
+ local tasks_done="" tasks_total="" raw_tasks failed_rule="" failed_command=""
host=$(hostname 2>/dev/null || echo "unknown")
@@ -204,9 +220,18 @@ function write_pipeline_state_marker_job() {
tasks_total=$(printf '%s' "${raw_tasks}" | grep -oP '(?<=of )\d+' || true)
fi
+ if [[ -f "${SNAKEMAKE_LOG_PATH}" ]]; then
+ failed_rule=$(grep -E 'Error in rule ' "${SNAKEMAKE_LOG_PATH}" 2>/dev/null | tail -1 | sed -E 's/.*Error in rule ([^:]+):.*/\1/' || true)
+ if [[ -n "${failed_rule}" ]]; then
+ failed_command=$(grep -n 'Error in rule ' "${SNAKEMAKE_LOG_PATH}" 2>/dev/null | tail -1 | cut -d: -f1 | while read -r line_no; do
+ sed -n "1,${line_no}p" "${SNAKEMAKE_LOG_PATH}" 2>/dev/null | grep -E 'computeMatrix|plotHeatmap|plotProfile|bamCoverage|samtools|python .*|bedToBam|bedSort|bwa|fastqc|multiqc|Rscript' | tail -1 || true
+ done | tail -1 || true)
+ fi
+ fi
+
_pipeline_write_status_json "${state}" "${reason}" "${job_id}" "${host}" \
"${submission_ts}" "${_START_TS:-}" "${duration_seconds}" "${exit_code}" \
- "${tasks_done}" "${tasks_total}"
+ "${tasks_done}" "${tasks_total}" "${failed_rule}" "${failed_command}"
}
# Best-effort jobby summary of the snakemake log. Path is parameterized via