From bcfe49fa50a3e110a5ffebb8e422442fe73cb0ed Mon Sep 17 00:00:00 2001 From: kopardev Date: Mon, 28 Sep 2026 12:08:50 -0400 Subject: [PATCH 1/7] docs: add deeptools output explanations --- docs/outputs.md | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/docs/outputs.md b/docs/outputs.md index 82af38c..82ad850 100644 --- a/docs/outputs.md +++ b/docs/outputs.md @@ -143,6 +143,25 @@ 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 can generate deepTools heatmaps and profile plots to summarize ATAC-seq accessibility around transcription start sites (TSS) and across gene bodies. These are useful QC and exploratory summaries: they make it easy to see whether signal is concentrated near promoters, spread across gene bodies, or strongly driven by a subset of samples or groups. + +These outputs are built from the pipeline's normalized coverage signal, using the default `deduplicated` read coverage unless a spike-in-normalized branch is enabled. The default annotation set is `protein_coding`, and the resulting filenames make that choice explicit so users can tell which gene model was used. + +| Output type | Meaning | Typical interpretation | +| --- | --- | --- | +| `*_protein_coding.TSS_heatmap.pdf` | Heatmap of signal aligned at each gene TSS | Strong promoter-centered enrichment usually indicates open chromatin near active promoters | +| `*_protein_coding.TSS_profile.pdf` | Average profile of signal around TSS | Helps compare the shape of promoter accessibility across samples and groups | +| `*_protein_coding.metagene_heatmap.pdf` | Heatmap across scaled gene bodies and flanks | Shows whether accessibility is distributed over gene bodies or concentrated near promoters | +| `*_protein_coding.metagene_profile.pdf` | Average gene-body profile | Summarizes whether signal rises near promoters, remains across the body, or falls toward the 3' end | +| `all_samples.*` | Pooled experiment-level summary | Best for the overall pattern across the full dataset | +| `group-level` or sample-level files | Per-group or per-sample summaries | Useful for identifying whether a pooled effect is driven by one group or by a shared biological pattern | + +A TSS heatmap is a matrix of signal centered on each gene's transcription start site; each row is a region and each column is a position relative to that TSS. A metagene plot instead summarizes signal across many genes after aligning them by their start and end coordinates, so the user can see the average signal across promoter, gene body, and downstream regions. In practice, a strong TSS-centered signal often reflects promoter accessibility, while a broader gene-body signal suggests a more distributed accessibility pattern. + +The documentation should be interpreted in context: deepTools outputs are exploratory summaries, not a replacement for peak calling or differential accessibility results. They are most useful alongside FRiP, TSS enrichment, library complexity, and peak annotation QC metrics. If spike-in normalization is enabled, the spike-in-normalized plots are especially helpful for making sample-to-sample comparisons when total sequencing depth differs substantially. + !!! note BAM files from `dedupBam` can be used for downstream footprinting analysis using [CCBR_TOBIAS](https://github.com/CCBR/CCBR_Tobias) pipeline From ff079b16f55993e2d204ce5d54635b50a49ac1ad Mon Sep 17 00:00:00 2001 From: kopardev Date: Mon, 28 Sep 2026 12:19:45 -0400 Subject: [PATCH 2/7] fix(deeptools): make matrix windows configurable --- config/config.yaml | 11 +++++ resources/cluster.json | 10 +++++ workflow/rules/align.smk | 92 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 113 insertions(+) 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/resources/cluster.json b/resources/cluster.json index 34f6344..824a8b9 100755 --- a/resources/cluster.json +++ b/resources/cluster.json @@ -30,6 +30,16 @@ "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" + }, "align": { "mem": "240g", "threads": "16", diff --git a/workflow/rules/align.smk b/workflow/rules/align.smk index c8a7b2c..409a833 100644 --- a/workflow/rules/align.smk +++ b/workflow/rules/align.smk @@ -481,3 +481,95 @@ 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 bundled TSS BED archive. +# Output: reference-point and scale-region matrices for each group and annotation class. +# """ + input: + prep=join(RESULTSDIR, "deeptools", "temp", "{group}.{bamtype}.reads.deeptools_prep"), + bedtar=join(RESOURCESDIR, "tssBed", "hg38_tssbeds.tar.gz"), + output: + tss=join(RESULTSDIR, "deeptools", "temp", "{group}.{bamtype}.{bedtype}.TSS.mat.gz"), + metagene=join(RESULTSDIR, "deeptools", "temp", "{group}.{bamtype}.{bedtype}.metagene.mat.gz"), + params: + bedtype="{bedtype}", + 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} +""" + From 743addfd8f07fd5aaacc34abb8459523034ce754 Mon Sep 17 00:00:00 2001 From: kopardev Date: Mon, 28 Sep 2026 12:24:17 -0400 Subject: [PATCH 3/7] fix(logging): record failed rule in pipeline status MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Include the failing Snakemake rule and command in pipeline.status.json and pipeline.failed so failed runs point directly to the offending step.\n\n⚡ Generated using AI ⚡ --- CHANGELOG.md | 1 + workflow/scripts/ccbr_pipeline_logging.sh | 29 +++++++++++++++++++++-- 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 401b6f4..3ee6f2a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,6 @@ ## ASPEN development version +- 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/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 From 754e42ac57880d9901744f44fac616079fff7179 Mon Sep 17 00:00:00 2001 From: kopardev Date: Mon, 28 Sep 2026 12:48:30 -0400 Subject: [PATCH 4/7] feat(dryrun): suppress NEXT prompt when nothing to be done MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Skip the "Submit run with..." NEXT hint after a dry-run that reports "Nothing to be done", since there is no pending work to submit. ⚡ Generated using AI ⚡ --- aspen | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) 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 } ########################################################################################## From 9cccf04dedc6144fff80aa687f649055ee1844ef Mon Sep 17 00:00:00 2001 From: kopardev Date: Mon, 28 Sep 2026 12:48:43 -0400 Subject: [PATCH 5/7] feat(deeptools): add plotHeatmap/plotProfile rule and wire outputs into rule all MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit deeptools_prep/deeptools_mat outputs were never requested by rule all (or any other rule), so Snakemake never scheduled them regardless of config/container fixes -- dry-run always reported "Nothing to be done" even when the .mat.gz files had never been created. - Add deeptools_plot rule: runs plotHeatmap/plotProfile on the TSS and metagene matrices to produce the final, user-facing heatmap/profile PDFs (results/deeptools/{group}.{bamtype}.{TSS,metagene}.{heatmap,profile}.pdf). - Wire the new PDF outputs into rule all so they're actually scheduled. - 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, since ASPEN supports multiple genomes. - Drop the unused {bedtype} wildcard from deeptools_mat/deeptools_plot; it was never populated by any config/rule and unused in the shell logic. - Add a deeptools_plot entry to resources/cluster.json. - Update docs/outputs.md to describe the actual output filenames and location instead of the previously documented, never-implemented protein_coding/all_samples naming scheme. Validated with a real dry-run against a live workdir: previously reported "Nothing to be done", now correctly schedules deeptools_prep, deeptools_mat, and deeptools_plot jobs. ⚡ Generated using AI ⚡ --- CHANGELOG.md | 1 + docs/outputs.md | 16 +++++------ resources/cluster.json | 5 ++++ workflow/Snakefile | 5 ++++ workflow/rules/align.smk | 62 ++++++++++++++++++++++++++++++++++++---- 5 files changed, 75 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3ee6f2a..7ec9444 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,6 @@ ## 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) diff --git a/docs/outputs.md b/docs/outputs.md index 82ad850..7b5a870 100644 --- a/docs/outputs.md +++ b/docs/outputs.md @@ -145,18 +145,18 @@ Content details: ### deepTools TSS and metagene outputs -ASPEN can generate deepTools heatmaps and profile plots to summarize ATAC-seq accessibility around transcription start sites (TSS) and across gene bodies. These are useful QC and exploratory summaries: they make it easy to see whether signal is concentrated near promoters, spread across gene bodies, or strongly driven by a subset of samples or groups. +ASPEN generates deepTools heatmaps and profile plots to summarize ATAC-seq accessibility around transcription start sites (TSS) and across gene bodies. These are useful QC and exploratory summaries: they make it easy to see whether signal is concentrated near promoters, spread across gene bodies, or strongly driven by a subset of samples or groups. -These outputs are built from the pipeline's normalized coverage signal, using the default `deduplicated` read coverage unless a spike-in-normalized branch is enabled. The default annotation set is `protein_coding`, and the resulting filenames make that choice explicit so users can tell which gene model was used. +These outputs are built per sample group (as defined in the sample manifest) and per BAM type (`dedup`/`nondedup`), from that group's pooled, spike-in-scaled `reads` bigWig coverage (`results/visualization/{bamtype}/reads_bigwig/`); if spike-in normalization is disabled, the scaling factor defaults to 1.0. All outputs live under `results/deeptools/`, with intermediate `computeMatrix` matrices (`*.mat.gz`) kept under `results/deeptools/temp/` for debugging or re-runs. | Output type | Meaning | Typical interpretation | | --- | --- | --- | -| `*_protein_coding.TSS_heatmap.pdf` | Heatmap of signal aligned at each gene TSS | Strong promoter-centered enrichment usually indicates open chromatin near active promoters | -| `*_protein_coding.TSS_profile.pdf` | Average profile of signal around TSS | Helps compare the shape of promoter accessibility across samples and groups | -| `*_protein_coding.metagene_heatmap.pdf` | Heatmap across scaled gene bodies and flanks | Shows whether accessibility is distributed over gene bodies or concentrated near promoters | -| `*_protein_coding.metagene_profile.pdf` | Average gene-body profile | Summarizes whether signal rises near promoters, remains across the body, or falls toward the 3' end | -| `all_samples.*` | Pooled experiment-level summary | Best for the overall pattern across the full dataset | -| `group-level` or sample-level files | Per-group or per-sample summaries | Useful for identifying whether a pooled effect is driven by one group or by a shared biological pattern | +| `{group}.{bamtype}.TSS.heatmap.pdf` | Heatmap of signal aligned at each gene TSS | Strong promoter-centered enrichment usually indicates open chromatin near active promoters | +| `{group}.{bamtype}.TSS.profile.pdf` | Average profile of signal around TSS | Helps compare the shape of promoter accessibility across samples and groups | +| `{group}.{bamtype}.metagene.heatmap.pdf` | Heatmap across scaled gene bodies and flanks | Shows whether accessibility is distributed over gene bodies or concentrated near promoters | +| `{group}.{bamtype}.metagene.profile.pdf` | Average gene-body profile | Summarizes whether signal rises near promoters, remains across the body, or falls toward the 3' end | + +`{group}` is the sample name from the manifest (so a pooled/experiment-wide summary can be produced by defining a group whose replicates span the whole dataset); `{bamtype}` is `dedup` or `nondedup`. A TSS heatmap is a matrix of signal centered on each gene's transcription start site; each row is a region and each column is a position relative to that TSS. A metagene plot instead summarizes signal across many genes after aligning them by their start and end coordinates, so the user can see the average signal across promoter, gene body, and downstream regions. In practice, a strong TSS-centered signal often reflects promoter accessibility, while a broader gene-body signal suggests a more distributed accessibility pattern. diff --git a/resources/cluster.json b/resources/cluster.json index 824a8b9..259b3d3 100755 --- a/resources/cluster.json +++ b/resources/cluster.json @@ -40,6 +40,11 @@ "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 409a833..c24228d 100644 --- a/workflow/rules/align.smk +++ b/workflow/rules/align.smk @@ -517,17 +517,16 @@ 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 bundled TSS BED archive. -# Output: reference-point and scale-region matrices for each group and annotation class. +# 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=join(RESOURCESDIR, "tssBed", "hg38_tssbeds.tar.gz"), + bedtar=config[GENOME]["tssBed"], output: - tss=join(RESULTSDIR, "deeptools", "temp", "{group}.{bamtype}.{bedtype}.TSS.mat.gz"), - metagene=join(RESULTSDIR, "deeptools", "temp", "{group}.{bamtype}.{bedtype}.metagene.mat.gz"), + tss=join(RESULTSDIR, "deeptools", "temp", "{group}.{bamtype}.TSS.mat.gz"), + metagene=join(RESULTSDIR, "deeptools", "temp", "{group}.{bamtype}.metagene.mat.gz"), params: - bedtype="{bedtype}", group="{group}", bamtype="{bamtype}", tss_upstream=config["deeptools"]["tss_upstream"], @@ -573,3 +572,54 @@ computeMatrix scale-regions \ -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" +""" + From 32bfdfdd260c1e4abdd57a4cee8320a49aff0cd6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 29 Sep 2026 15:41:16 +0000 Subject: [PATCH 6/7] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- workflow/rules/align.smk | 1 - 1 file changed, 1 deletion(-) diff --git a/workflow/rules/align.smk b/workflow/rules/align.smk index c24228d..74c88ce 100644 --- a/workflow/rules/align.smk +++ b/workflow/rules/align.smk @@ -622,4 +622,3 @@ plotProfile \ --startLabel "start" \ --endLabel "end" """ - From 5adbde99cf40972a85fffaf689be527aedaefa17 Mon Sep 17 00:00:00 2001 From: kopardev Date: Tue, 29 Sep 2026 17:39:43 -0400 Subject: [PATCH 7/7] docs(outputs): clarify deepTools summaries and widen tables MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Document when deepTools outputs are generated, how to interpret them, and which normalization and annotation defaults apply. Expand the documentation content area and keep table columns aligned on large displays. _commit message is ai-generated_ ⚡ Generated using AI ⚡ --- docs/css/custom.css | 41 +++++++++++++++++++++++++++++++++++++ docs/outputs.md | 49 ++++++++++++++++++++++++++++++++++----------- 2 files changed, 78 insertions(+), 12 deletions(-) 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 7b5a870..f5fbb2f 100644 --- a/docs/outputs.md +++ b/docs/outputs.md @@ -145,22 +145,46 @@ Content details: ### 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 are useful QC and exploratory summaries: they make it easy to see whether signal is concentrated near promoters, spread across gene bodies, or strongly driven by a subset of samples or groups. +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. -These outputs are built per sample group (as defined in the sample manifest) and per BAM type (`dedup`/`nondedup`), from that group's pooled, spike-in-scaled `reads` bigWig coverage (`results/visualization/{bamtype}/reads_bigwig/`); if spike-in normalization is disabled, the scaling factor defaults to 1.0. All outputs live under `results/deeptools/`, with intermediate `computeMatrix` matrices (`*.mat.gz`) kept under `results/deeptools/temp/` for debugging or re-runs. +#### Purpose and when generated -| Output type | Meaning | Typical interpretation | -| --- | --- | --- | -| `{group}.{bamtype}.TSS.heatmap.pdf` | Heatmap of signal aligned at each gene TSS | Strong promoter-centered enrichment usually indicates open chromatin near active promoters | -| `{group}.{bamtype}.TSS.profile.pdf` | Average profile of signal around TSS | Helps compare the shape of promoter accessibility across samples and groups | -| `{group}.{bamtype}.metagene.heatmap.pdf` | Heatmap across scaled gene bodies and flanks | Shows whether accessibility is distributed over gene bodies or concentrated near promoters | -| `{group}.{bamtype}.metagene.profile.pdf` | Average gene-body profile | Summarizes whether signal rises near promoters, remains across the body, or falls toward the 3' end | +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. -`{group}` is the sample name from the manifest (so a pooled/experiment-wide summary can be produced by defining a group whose replicates span the whole dataset); `{bamtype}` is `dedup` or `nondedup`. +#### Output files -A TSS heatmap is a matrix of signal centered on each gene's transcription start site; each row is a region and each column is a position relative to that TSS. A metagene plot instead summarizes signal across many genes after aligning them by their start and end coordinates, so the user can see the average signal across promoter, gene body, and downstream regions. In practice, a strong TSS-centered signal often reflects promoter accessibility, while a broader gene-body signal suggests a more distributed accessibility pattern. +| 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 | -The documentation should be interpreted in context: deepTools outputs are exploratory summaries, not a replacement for peak calling or differential accessibility results. They are most useful alongside FRiP, TSS enrichment, library complexity, and peak annotation QC metrics. If spike-in normalization is enabled, the spike-in-normalized plots are especially helpful for making sample-to-sample comparisons when total sequencing depth differs substantially. +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 @@ -551,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.