From 38bbb7da321a68e9a1ea1c178ffc93c2552e0335 Mon Sep 17 00:00:00 2001 From: kopardev Date: Thu, 10 Sep 2026 11:00:52 -0400 Subject: [PATCH 1/6] chore(gitignore): ignore local plans/ directory MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add plans/ alongside existing plans_reviews/ for local, untracked planning docs (e.g. dedup vs filtered counts matrix plan, tracked in issue #138). ⚡ Generated using AI ⚡ --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 380880a..a09f180 100644 --- a/.gitignore +++ b/.gitignore @@ -6,5 +6,6 @@ aspen.bak site/ tmp_test/ plans_reviews/ +plans/ .snakemake/log/ node_modules From eb9df5d03a001c0fc84d1bced3df1e3a40dc7840 Mon Sep 17 00:00:00 2001 From: kopardev Date: Fri, 11 Sep 2026 13:48:37 -0400 Subject: [PATCH 2/6] feat(diffatac): generate dedup and nondedup counts/DiffATAC outputs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Generate Tn5 nicking-site and read counts matrices, and the corresponding DiffATAC/DESeq2 results, from both dedup.bam (PCR/optical duplicates removed) and filtered.bam (duplicates retained, labeled nondedup), written to separate dedup/nondedup output subfolders. Previously only the duplicate-retaining filtered.bam was used. dedup is recommended for standard differential accessibility testing. Refs #138 ⚡ Generated using AI ⚡ --- VERSION | 2 +- workflow/Snakefile | 16 ++++++++-------- workflow/rules/diffatac.smk | 21 +++++++++++---------- workflow/rules/init.smk | 5 +++++ 4 files changed, 25 insertions(+), 19 deletions(-) diff --git a/VERSION b/VERSION index 26aaba0..0c64b5c 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -1.2.0 +1.2.0-dev diff --git a/workflow/Snakefile b/workflow/Snakefile index faf1515..723aaa3 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -57,20 +57,20 @@ rule all: # spikein counts and scaling factors expand(join(RESULTSDIR, "spikein", "{replicate}", "{replicate}.counts"), replicate=REPLICATES), join(RESULTSDIR, "spikein", "scaling_factors.tsv"), - expand(join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "ROI.{peakcaller}.{method}_scaled_counts.tsv"), method=COUNTING_METHODS, peakcaller=PEAKCALLERS), + expand(join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "{bamtype}", "ROI.{peakcaller}.{method}_scaled_counts.tsv"), method=COUNTING_METHODS, peakcaller=PEAKCALLERS, bamtype=BAM_TYPES), # visualizations - expand(join(RESULTSDIR, "visualization", "{method}_bam", "{replicate}.{method}.bam"), method=COUNTING_METHODS, replicate=REPLICATES), - expand(join(RESULTSDIR, "visualization", "{method}_bigwig", "{replicate}.{method}.bw"), method=COUNTING_METHODS, replicate=REPLICATES), - expand(join(RESULTSDIR, "visualization", "reads_bed", "{replicate}.reads.bed.gz"), replicate=REPLICATES), # these will be used by chromVar + expand(join(RESULTSDIR, "visualization", "{bamtype}", "{method}_bam", "{replicate}.{method}.bam"), method=COUNTING_METHODS, replicate=REPLICATES, bamtype=BAM_TYPES), + expand(join(RESULTSDIR, "visualization", "{bamtype}", "{method}_bigwig", "{replicate}.{method}.bw"), method=COUNTING_METHODS, replicate=REPLICATES, bamtype=BAM_TYPES), + expand(join(RESULTSDIR, "visualization", "{bamtype}", "reads_bed", "{replicate}.reads.bed.gz"), replicate=REPLICATES, bamtype=BAM_TYPES), # these will be used by chromVar # regions of interest (ROI) gtf expand(join(PEAKSDIR, "{peakcaller}", "fixed_width", "ROI.{peakcaller}.gtf"), peakcaller=PEAKCALLERS), expand(join(PEAKSDIR, "{peakcaller}", "fixed_width", "ROI.{peakcaller}.bed.annotated.gz"), peakcaller=PEAKCALLERS), - # ROI counts matrix using both tn5sites and reads for each peakcaller - expand(join(PEAKSDIR, "{peakcaller}", "fixed_width", "counts", "ROI.{peakcaller}.{method}_counts.tsv"), method=COUNTING_METHODS, peakcaller=PEAKCALLERS), + # ROI counts matrix using both tn5sites and reads for each peakcaller, for both dedup and nondedup BAM types + expand(join(PEAKSDIR, "{peakcaller}", "fixed_width", "counts", "{bamtype}", "ROI.{peakcaller}.{method}_counts.tsv"), method=COUNTING_METHODS, peakcaller=PEAKCALLERS, bamtype=BAM_TYPES), # diffatac - expand(join(PEAKSDIR, "{peakcaller}", "DiffATAC", "{method}", "degs.done"), method=COUNTING_METHODS, peakcaller=PEAKCALLERS) if CONTRASTS.shape[0] > 0 else [], + 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", "{method}", "all_diff_atacs.tsv"), method=COUNTING_METHODS, peakcaller=PEAKCALLERS) if CONTRASTS.shape[0] > 0 else [], + 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 [], on_complete = f""" for cmd in spooker jobby; do diff --git a/workflow/rules/diffatac.smk b/workflow/rules/diffatac.smk index 2ecbd53..e0d7e91 100644 --- a/workflow/rules/diffatac.smk +++ b/workflow/rules/diffatac.smk @@ -109,13 +109,14 @@ ls -alrth $(dirname {output.roi_annotated}) rule get_counts_table: input: bam_files = lambda wildcards: expand( - join(RESULTSDIR, "visualization", "{method}_bam", "{replicate}.{method}.bam"), + join(RESULTSDIR, "visualization", "{bamtype}", "{method}_bam", "{replicate}.{method}.bam"), + bamtype=wildcards.bamtype, method=wildcards.method, replicate=REPLICATES ), gtf = join(RESULTSDIR,"peaks","{peakcaller}","fixed_width","ROI.{peakcaller}.gtf"), output: - counts = join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "ROI.{peakcaller}.{method}_counts.tsv"), + counts = join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "{bamtype}", "ROI.{peakcaller}.{method}_counts.tsv"), params: scriptsdir = SCRIPTSDIR, script = "_featureCounts_header_fix.py" @@ -149,9 +150,9 @@ localrules: scale_counts_table rule scale_counts_table: input: scaling_factors = rules.compute_scaling_factors.output.scaling_factors, - counts = join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "ROI.{peakcaller}.{method}_counts.tsv"), + counts = join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "{bamtype}", "ROI.{peakcaller}.{method}_counts.tsv"), output: - scaledcounts = join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "ROI.{peakcaller}.{method}_scaled_counts.tsv") + scaledcounts = join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "{bamtype}", "ROI.{peakcaller}.{method}_scaled_counts.tsv") params: scriptsdir=SCRIPTSDIR, script="_scale_counts.py", @@ -168,9 +169,9 @@ python {params.scriptsdir}/{params.script} --counts {input.counts} --scaling_fac rule diffatac: input: - counts = join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "ROI.{peakcaller}.{method}_counts.tsv") + counts = join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "{bamtype}", "ROI.{peakcaller}.{method}_counts.tsv") output: - degsdone = join(PEAKSDIR, "{peakcaller}", "DiffATAC", "{method}", "degs.done") + degsdone = join(PEAKSDIR, "{peakcaller}", "DiffATAC", "{bamtype}", "{method}", "degs.done") params: contrasts = config['contrasts'], scriptsdir = SCRIPTSDIR, @@ -244,11 +245,11 @@ touch {output.degsdone} rule diffatac_aggregate: input: - counts = join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "ROI.{peakcaller}.{method}_counts.tsv"), - degsdone = join(RESULTSDIR,"peaks","{peakcaller}","DiffATAC","{method}","degs.done"), + counts = join(RESULTSDIR, "peaks", "{peakcaller}", "fixed_width", "counts", "{bamtype}", "ROI.{peakcaller}.{method}_counts.tsv"), + degsdone = join(RESULTSDIR,"peaks","{peakcaller}","DiffATAC","{bamtype}","{method}","degs.done"), output: - alldegshtml = join(RESULTSDIR,"peaks","{peakcaller}","DiffATAC","{method}","all_diff_atacs.html"), - alldegstsv = join(RESULTSDIR,"peaks","{peakcaller}","DiffATAC","{method}","all_diff_atacs.tsv"), + alldegshtml = join(RESULTSDIR,"peaks","{peakcaller}","DiffATAC","{bamtype}","{method}","all_diff_atacs.html"), + alldegstsv = join(RESULTSDIR,"peaks","{peakcaller}","DiffATAC","{bamtype}","{method}","all_diff_atacs.tsv"), params: contrasts = config['contrasts'], scriptsdir = SCRIPTSDIR, diff --git a/workflow/rules/init.smk b/workflow/rules/init.smk index b490fb9..4200e09 100644 --- a/workflow/rules/init.smk +++ b/workflow/rules/init.smk @@ -155,6 +155,11 @@ except KeyError: PEAKCALLERS = ["genrich", "macs2"] COUNTING_METHODS = ["tn5sites", "reads"] +# counts matrices (and DiffATAC results) are generated from both the +# duplicate-removed "dedup" BAM and the duplicate-retaining "nondedup" +# (filtered) BAM, so end users can compare/choose which to use downstream. +# See docs/extra.md and issue #138 for background/guidance. +BAM_TYPES = ["dedup", "nondedup"] if contrastsfileexists: check_readaccess(config["contrasts"]) try: From da77b6d3a4a722a8b718000a2fa3c30863e9ae92 Mon Sep 17 00:00:00 2001 From: kopardev Date: Fri, 11 Sep 2026 13:49:01 -0400 Subject: [PATCH 3/6] fix(spikein): fail fast with clear error on 0 spike-in reads MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add an explicit check in _compute_downsampling_scaling_factors.py that exits with an actionable error message when a replicate has 0 reads aligned to the spike-in genome, instead of crashing with an opaque ZeroDivisionError. The compute_scaling_factors rule now also captures the script's stderr in a dedicated results/spikein/compute_scaling_factors.log so the message isn't buried in the main snakemake log. Fixes #139 ⚡ Generated using AI ⚡ --- workflow/rules/align.smk | 37 ++++++++++++++----- .../_compute_downsampling_scaling_factors.py | 17 +++++++++ 2 files changed, 44 insertions(+), 10 deletions(-) diff --git a/workflow/rules/align.smk b/workflow/rules/align.smk index 12c09d4..8d6900a 100644 --- a/workflow/rules/align.smk +++ b/workflow/rules/align.smk @@ -288,6 +288,8 @@ rule compute_scaling_factors: counts=expand(join(RESULTSDIR,"spikein","{replicate}","{replicate}.counts"),replicate=REPLICATES), output: scaling_factors=join(RESULTSDIR,"spikein","scaling_factors.tsv"), + log: + join(RESULTSDIR,"spikein","compute_scaling_factors.log"), params: replicates=REPLICATES, workdir=RESULTSDIR, @@ -306,7 +308,7 @@ unset PYTHONPATH # Check if spike-in is enabled if [[ "{params.spikein}" == "true" ]];then # Compute scaling factors using the provided script - cat {input.counts} | python {params.scriptsdir}/{params.script} | sort > {output.scaling_factors} + cat {input.counts} | python {params.scriptsdir}/{params.script} 2> {log} | sort > {output.scaling_factors} else # If spike-in is not enabled, set scaling factor to 1.0 for all replicates cat {input.counts} | sort > ${{TMPDIR}}/allcounts.txt @@ -320,25 +322,39 @@ fi ######################################################### +def get_bam_by_bamtype(wildcards): +# """ +# Resolve the source BAM for a given {bamtype} wildcard value: +# dedup -> dedup.bam (PCR/optical duplicates removed by Picard MarkDuplicates) +# nondedup -> filtered.bam (quality-filtered but duplicates retained) +# """ + if wildcards.bamtype == "dedup": + return join(ALIGNDIR,"dedupBam",wildcards.replicate+".dedup.bam") + elif wildcards.bamtype == "nondedup": + return join(ALIGNDIR,"filteredBam",wildcards.replicate+".filtered.bam") + else: + raise ValueError("Unknown bamtype wildcard: "+wildcards.bamtype) + rule create_tn5bams: # """ # Create tn5 and reads bam files for creating bigwigs -# Input: filtered BAM file from alignment step +# Input: dedup or filtered (nondedup) BAM file from alignment step, selected via {bamtype} # Output: tn5 and reads BAM files -# This rule generates tn5 and reads BAM files from the filtered BAM file obtained after alignment. +# This rule generates tn5 and reads BAM files from the {bamtype} BAM file obtained after alignment. # The tn5 BAM file contains the positions of Tn5 insertions, while the reads BAM file contains # the positions of the filtered original reads. These BAM files are used to create bigwig files for # visualization in genome browsers. # """ # group: "TAD" input: - bam=rules.align.output.filteredBam, + bam=get_bam_by_bamtype, output: - tn5bam=join(RESULTSDIR,"visualization","tn5sites_bam","{replicate}.tn5sites.bam"), - readsbed=join(RESULTSDIR,"visualization","reads_bed","{replicate}.reads.bed.gz"), # these will be used by chromVar - readsbam=join(RESULTSDIR,"visualization","reads_bam","{replicate}.reads.bam"), + tn5bam=join(RESULTSDIR,"visualization","{bamtype}","tn5sites_bam","{replicate}.tn5sites.bam"), + readsbed=join(RESULTSDIR,"visualization","{bamtype}","reads_bed","{replicate}.reads.bed.gz"), # these will be used by chromVar + readsbam=join(RESULTSDIR,"visualization","{bamtype}","reads_bam","{replicate}.reads.bam"), params: replicate="{replicate}", + bamtype="{bamtype}", genome=GENOME, genomefile=GENOMEFILE, scriptsdir=SCRIPTSDIR, @@ -396,7 +412,7 @@ samtools index {output.readsbam} rule create_bigwigs: # """ # Create spike-in scaled bigwigs for visualization -# Input: tn5 and reads BAM files, scaling factors +# Input: tn5 and reads BAM files (per {bamtype}: dedup or nondedup), scaling factors # Output: tn5 and reads bigwig files # This rule generates bigwig files from tn5 and reads BAM files, applying the scaling factors # computed from spike-in counts to normalize the coverage. The resulting bigwig files can be @@ -407,11 +423,12 @@ rule create_bigwigs: readsbam=rules.create_tn5bams.output.readsbam, scaling_factors=rules.compute_scaling_factors.output.scaling_factors output: - tn5bw=join(RESULTSDIR,"visualization","tn5sites_bigwig","{replicate}.tn5sites.bw"), - readsbw=join(RESULTSDIR,"visualization","reads_bigwig","{replicate}.reads.bw"), + tn5bw=join(RESULTSDIR,"visualization","{bamtype}","tn5sites_bigwig","{replicate}.tn5sites.bw"), + readsbw=join(RESULTSDIR,"visualization","{bamtype}","reads_bigwig","{replicate}.reads.bw"), params: spikein=str(SPIKEIN).lower(), replicate="{replicate}", + bamtype="{bamtype}", scriptsdir=SCRIPTSDIR, script="_print_replicate_scaling_factor.py" container: config["deeptoolsdocker"] diff --git a/workflow/scripts/_compute_downsampling_scaling_factors.py b/workflow/scripts/_compute_downsampling_scaling_factors.py index 58b8b9f..01bcdb0 100644 --- a/workflow/scripts/_compute_downsampling_scaling_factors.py +++ b/workflow/scripts/_compute_downsampling_scaling_factors.py @@ -7,6 +7,23 @@ for row in data: row[1] = int(row[1]) +# Fail loudly and clearly if any replicate has 0 spike-in reads, rather than +# letting the division below raise an opaque ZeroDivisionError. +zero_count_replicates = [row[0] for row in data if row[1] == 0] +if zero_count_replicates: + sys.exit( + "ERROR: the following replicate(s) have 0 reads aligned to the " + f"spike-in genome: {', '.join(zero_count_replicates)}.\n" + "Spike-in scaling factors cannot be computed. This usually means " + "either (a) the spike-in material was not actually present in this " + "library, (b) the spike-in genome/index is misconfigured, or (c) " + "this sample should not use spike-in normalization at all.\n" + "Fix by verifying the spike-in genome/index path in config.yaml, " + "confirming spike-in was added during library prep, or excluding " + "this replicate / setting `spikein: false` if none of your samples " + "have spike-in material." + ) + # Find the smallest value in column 2 min_value = min(row[1] for row in data) From 54d25f08fe3dc34b5849e3b25619940fd199384e Mon Sep 17 00:00:00 2001 From: kopardev Date: Fri, 11 Sep 2026 13:49:23 -0400 Subject: [PATCH 4/6] docs: overhaul documentation for accuracy and ELI15 clarity MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Audit and fix documentation across docs/overview.md, docs/outputs.md, docs/introduction.md, docs/deployment.md, and docs/index.md: - Fix stale claims left over from the dedup/nondedup change. - Clarify how consensus peaks (pooled/consensus.bed/fixed-width ROI) and regions of interest are generated, with a why-two-rounds explanation and a table of which output to use for which downstream task. - Add ENCODE-referenced QC rule-of-thumb callouts for library complexity (NRF/PBC1/PBC2), TSS enrichment, and FRiP. - Consolidate scattered spike-in "when should I use this" guidance (previously split across overview.md and deployment.md) into a single decision-tree callout in overview.md, cross-linked from deployment.md. - Note that the embedded `aspen --help` output in deployment.md is an illustrative, point-in-time snapshot. - Remove the outdated, unreferenced docs/extra.md page; its content is superseded by docs/outputs.md. Update CHANGELOG.md with entries for this docs pass and the spike-in zero-reads fix. Refs #138, #139 ⚡ Generated using AI ⚡ --- CHANGELOG.md | 6 ++ docs/deployment.md | 6 +- docs/extra.md | 5 -- docs/index.md | 5 ++ docs/introduction.md | 16 ++--- docs/outputs.md | 166 +++++++++++++++++++++++++++++++------------ docs/overview.md | 86 ++++++++++++++++------ 7 files changed, 207 insertions(+), 83 deletions(-) delete mode 100644 docs/extra.md diff --git a/CHANGELOG.md b/CHANGELOG.md index f7a58fe..9a3bd7a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,9 @@ +## ASPEN 1.3.0 + +- Generate Tn5 nicking-site and read counts matrices, and the corresponding DiffATAC/DESeq2 results, from both `dedup.bam` (PCR/optical duplicates removed) and `filtered.bam` (duplicates retained, labeled `nondedup`), written to separate `dedup`/`nondedup` output subfolders under `visualization/` and `peaks/{peakcaller}/{fixed_width/counts,DiffATAC}/`; previously only the duplicate-retaining `filtered.bam` was used. `dedup` is recommended for standard differential accessibility testing. (#138, @kopardev) +- Fail fast with a clear, actionable error message when a replicate has 0 reads aligned to the spike-in genome, instead of crashing with an opaque `ZeroDivisionError` in `_compute_downsampling_scaling_factors.py`; the error is now also captured in a dedicated `results/spikein/compute_scaling_factors.log`. (#139, @kopardev) +- Overhaul documentation for accuracy and clarity across `docs/overview.md`, `docs/outputs.md`, `docs/introduction.md`, `docs/deployment.md`, and `docs/index.md`: fix stale claims left over from the dedup/nondedup change, clarify how consensus peaks and ROIs are generated, add ENCODE-referenced QC rule-of-thumb callouts (library complexity, TSS enrichment, FRiP), consolidate scattered spike-in "when should I use this" guidance into a single decision-tree callout, and remove the outdated, unreferenced `docs/extra.md` page whose content is now covered by `docs/outputs.md`. (#138, @kopardev) + ## ASPEN 1.2.0 - Add `hs1_chrR` as a supported reference genome (T2T-CHM13 + chrR rDNA unit from [Paralkar lab](https://github.com/vikramparalkar/rDNA-Mapping-Genomes)). (#117, @kopardev) diff --git a/docs/deployment.md b/docs/deployment.md index 9a65142..ca96788 100644 --- a/docs/deployment.md +++ b/docs/deployment.md @@ -74,6 +74,8 @@ To explore all possible options of the `aspen` command you can either run it wit Here is what help looks like: +> **Note**: This is illustrative example output captured at doc-writing time — exact values (e.g. `pipeline_home`, `git commit/tag`, `aspen_version`) will differ depending on which ASPEN version/branch is installed at your site. Run `aspen --help` yourself to see the current values for your installation. + ```bash ########################################################################################## @@ -197,6 +199,8 @@ contrasts_fdr_cutoff: 0.05 ASPEN supports spike-in normalization, which is useful for controlling technical variability or comparing global shifts in chromatin accessibility across samples. Spike-in reads (e.g., from _Drosophila melanogaster_ or _E. coli_) are aligned separately and used to compute normalization factors that are applied to host genome accessibility counts. +Not sure whether you need this for your experiment? See ["Should I turn on spike-in normalization?"](overview.md#optional-spike-in-normalization) in the Overview docs for a decision guide. + To enable spike-in normalization, edit the `config.yaml` file that was generated during `init`. You can find it in your output directory (`/config.yaml`). Open the file and locate the following lines: @@ -236,8 +240,6 @@ Once enabled, ASPEN will: - Normalize accessibility counts using spike-in-derived scaling factors. - Report both normalized and raw counts in the output tables and reports. -This step is optional but highly recommended when you expect global changes in chromatin accessibility due to treatments or perturbations. - ### 🛠️ Dry Run the Pipeline Before executing the full analysis and after editing the `config.yaml` as needed, perform a dry run to visualize the workflow and identify potential issues: diff --git a/docs/extra.md b/docs/extra.md deleted file mode 100644 index 48d3b79..0000000 --- a/docs/extra.md +++ /dev/null @@ -1,5 +0,0 @@ -Once the ROIs are established, ASPEN generates two distinct count matrices: - -- **Tn5 Nicking Sites Count Matrix**: This matrix quantifies the frequency of Tn5 transposase insertion events at each ROI. The `filtered.bam` file (with PCR duplicates and not `dedup.bam`) is used for determining Tn5 sites and counting them. The Tn5 transposase preferentially inserts into accessible regions of the chromatin, and the number of insertion events serves as a proxy for chromatin accessibility. By counting these insertion sites, researchers can accurately infer the openness of chromatin regions under different experimental conditions. - -- **Read Counts Matrix**: This matrix records the number of sequencing reads mapped to each ROI. The `filtered.bam` file (with PCR duplicates and not `dedup.bam`) is used for counting. While Tn5 nicking sites provide a more direct measure of chromatin accessibility, read counts are included in ASPEN as they have been widely used in recent publications. Analyzing both matrices together offers a comprehensive view of chromatin accessibility dynamics. diff --git a/docs/index.md b/docs/index.md index 71209cb..978d5e0 100644 --- a/docs/index.md +++ b/docs/index.md @@ -1,3 +1,8 @@ ## Background +!!! tip "New here?" + If you just want to **run the pipeline**, skip ahead to + [Running ASPEN](deployment.md). Come back here later for the + scientific background on ATAC-seq and why ASPEN is built the way it is. + The Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) has revolutionized genomics by providing a rapid and sensitive method to assess chromatin accessibility across the genome. This technique offers profound insights into gene regulation, epigenetic modifications, and the dynamic landscape of the chromatin environment. To facilitate the comprehensive analysis of ATAC-seq data, the Center for Cancer Research (CCR) Collaborative Bioinformatics Resource (CCBR) has developed ASPEN (**A**tac **S**eq **P**ip**E**li**N**e), an automated, robust, reproducible pipeline designed using the Snakemake pipelining framework to streamline the complexities inherent in ATAC-seq data processing. diff --git a/docs/introduction.md b/docs/introduction.md index f5b6b1c..eb003f7 100644 --- a/docs/introduction.md +++ b/docs/introduction.md @@ -10,12 +10,10 @@ ATAC-seq is a powerful technique that enables the identification of open chromat Despite its advantages, ATAC-seq data analysis presents several challenges: -- **📊 Data Complexity**: The technique generates large datasets with varying fragment sizes corresponding to nucleosome-free regions, mono-nucleosomes, and multi-nucleosomes, necessitating sophisticated computational tools for accurate interpretation. - -- **✅ Quality Control**: Ensuring high-quality data requires meticulous assessment at multiple stages, including evaluating sequencing quality, fragment size distribution, and enrichment of reads in regulatory regions. - -- **📈 Peak Calling**: Identifying regions of open chromatin (peaks) demands precise computational methods to distinguish true signals from background noise, especially in datasets with low signal-to-noise ratios. - -- **🔬 Spike-in Normalization Challenges**: While spike-in controls are invaluable for normalization, their implementation requires careful consideration. Variability in spike-in recovery or differences in chromatin accessibility between samples can complicate normalization efforts, necessitating robust experimental design and computational strategies. - -- **🔄 Reproducibility**: Achieving consistent results across different experiments and conditions is crucial for the reliability of biological interpretations. +| Challenge | Why it's hard | How ASPEN helps | +|---|---|---| +| 📊 **Data complexity** | Datasets contain a mix of fragment sizes corresponding to nucleosome-free regions, mono-nucleosomes, and multi-nucleosomes, requiring careful computational interpretation | Reports a [fragment length distribution](overview.md#fragment-length-distribution) QC metric so you can visually confirm the expected nucleosome-free/mono-nucleosome pattern | +| ✅ **Quality control** | Requires assessing sequencing quality, fragment size distribution, and enrichment of reads in regulatory regions all together | Runs FastQC/MultiQC plus ATAC-specific metrics ([FRiP](overview.md#fraction-of-reads-in-peaks-frip), [TSS enrichment](overview.md#transcription-start-site-tss-enrichment), [library complexity](overview.md#library-complexity)) automatically for every sample | +| 📈 **Peak calling** | Distinguishing true accessible regions from background noise is especially hard in low signal-to-noise datasets | Runs both MACS2 and Genrich (see [Peak Detection](overview.md#peak-detection)) so results can be cross-checked against each other | +| 🔬 **Spike-in normalization** | Variability in spike-in recovery or accessibility differences between samples can complicate normalization | Provides optional, transparent [spike-in normalization](overview.md#optional-spike-in-normalization) with a documented failure mode if spike-in reads are missing | +| 🔄 **Reproducibility** | Achieving consistent results across experiments and conditions requires more than just consistent code | Runs entirely in versioned Docker/Singularity containers (see "Containerization" above) so the same environment is reproduced on every run | diff --git a/docs/outputs.md b/docs/outputs.md index 17ccfec..8822ee4 100644 --- a/docs/outputs.md +++ b/docs/outputs.md @@ -68,14 +68,28 @@ WORKDIR ├── peaks │   ├── genrich │   │   ├── DiffATAC - │   │   │   ├── reads - │   │   │   └── tn5sites + │   │   │   ├── dedup + │   │   │   │   ├── reads + │   │   │   │   └── tn5sites + │   │   │   └── nondedup + │   │   │   ├── reads + │   │   │   └── tn5sites │   │   └── fixed_width + │   │   └── counts + │   │   ├── dedup + │   │   └── nondedup │   └── macs2 │   ├── DiffATAC - │   │   ├── reads - │   │   └── tn5sites + │   │   ├── dedup + │   │   │   ├── reads + │   │   │   └── tn5sites + │   │   └── nondedup + │   │   ├── reads + │   │   └── tn5sites │   └── fixed_width + │   └── counts + │   ├── dedup + │   └── nondedup ├── QC │   ├── fastqc │   ├── fld @@ -96,11 +110,18 @@ WORKDIR │   ├── genrichReads │   └── trim └── visualization - ├── reads_bam - ├── reads_bed - ├── reads_bigwig - ├── tn5sites_bam - └── tn5sites_bigwig + ├── dedup + │ ├── reads_bam + │ ├── reads_bed + │ ├── reads_bigwig + │ ├── tn5sites_bam + │ └── tn5sites_bigwig + └── nondedup + ├── reads_bam + ├── reads_bed + ├── reads_bigwig + ├── tn5sites_bam + └── tn5sites_bigwig ``` Content details: @@ -108,19 +129,15 @@ Content details: | Folder | SubFolder | Description | | ------------- | ------------------- | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | | alignment | qsortedBam | - Query sorted Bowtie2 alignments in BAM format.
- Excludes unmapped and platform/vendor quality failing reads.
- Used for Genrich peak calling. | -| alignment | filteredBam | - Filtered BAM files after excluding non-primary, supplementary, and MAPQ <=5 alignments.
- Used for counting reads/tn5 nicks.
- Derived from `qsortedBam`. | -| alignment | dedupBam | - Deduplicated filtered BAM files.
- PCR or optical duplicates marked with PicardTools and excluded.
- Can be used downstream with CCBR_TOBIAS pipeline.
- Derived from `filteredBam`. | +| alignment | filteredBam | - Filtered BAM files after excluding non-primary, supplementary, and MAPQ <=5 alignments.
- This is the `nondedup` BAM (i.e. before duplicate removal) used for the `nondedup` counts/DiffATAC variant.
- Derived from `qsortedBam`. | +| alignment | dedupBam | - Deduplicated filtered BAM files (PCR/optical duplicates marked with PicardTools and excluded).
- Derived from `filteredBam`.
- Used for MACS2 peak calling, FRiP, and the recommended `dedup` counts/DiffATAC variant.
- Can also be used downstream with CCBR_TOBIAS pipeline. | | alignment | tagAlign | - `tagAlign.gz` files used for MACS2 peak calling.
- Derived from `dedupBam`. | -| peaks | genrich & macs | - Genrich/MACS2 peak calls (raw, consensus, fixed-width).
- Contains ROI files with Diff-ATAC results if `contrasts.tsv` is provided.
- Calculated with DESeq2 using both read counts and tn5 nicking sites in ROI. | +| peaks | genrich & macs | - Genrich/MACS2 peak calls (raw, consensus, fixed-width).
- Contains ROI files with Diff-ATAC results if `contrasts.tsv` is provided.
- Calculated with DESeq2 using both read counts and tn5 nicking sites in ROI.
- Counts matrices and DiffATAC results are generated separately for `dedup` (PCR duplicates removed, recommended) and `nondedup` (duplicates retained) BAMs, in like-named subfolders. | +| visualization | dedup / nondedup | - Tn5 nicking-site/read BAMs, BEDs, and bigwigs, split into `dedup` and `nondedup` subfolders depending on whether PCR duplicates were removed from the source BAM before generating them. | | QC | various | - Flagstats.
- Dupmetrics.
- Read counts.
- Motif enrichments.
- FLD stats.
- Fqscreen.
- FRiP.
- ChIPSeeker results.
- TSS enrichments.
- Preseq.
- Homer/AME motif enrichments.
- MultiQC. | | QC | peak_annotation | detailed peak annotations described below | | 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. | -| visualization | reads_bam | - Tn5 nick adjusted reads in BAM format.
- Derived from `filteredBam`. | -| visualization | reads_bed | - Tn5 nick adjusted reads in BED format.
- Derived from `reads_bam`.
- Can be used by ChromVar. | -| visualization | reads_bigwig | - Tn5 nick adjusted reads in BIGWIG format.
- Scaled using spike-in scaling factors if present.
- Derived from `reads_bam`. | -| visualization | tn5sites_bam | - Tn5 nicking sites in BAM format.
- Derived from `filteredBam`. | -| visualization | tn5sites_bigwig | - Tn5 nicking sites in BIGWIG format.
- Scaled using spike-in scaling factors if present.
- Derived from `tn5sites_bam`. | !!! note BAM files from `dedupBam` can be used for downstream footprinting analysis using [CCBR_TOBIAS](https://github.com/CCBR/CCBR_Tobias) pipeline @@ -133,16 +150,40 @@ BAM files from `dedupBam` can also be converted to BED format and processed with #### How consensus peaks are generated +**Why two rounds?** A single replicate's peak calls are noisy — some "peaks" +are just sequencing/technical artifacts that happened to occur in one +replicate. Pooling and filtering peaks *within* a sample (Round 1) gives a +reproducibility-checked, per-sample truth set. But differential accessibility +testing needs to compare counts *across* samples/conditions, which requires +one shared set of regions measured identically everywhere — a per-sample +peak set alone doesn't give you that. Round 2 solves this by merging the +per-sample consensus peaks into a single, fixed-width Region of Interest +(ROI) set shared by every sample, so the same coordinates can be quantified +and compared across the whole experiment. + ASPEN produces peaks at two levels of consensus. Understanding the distinction is important for interpreting output files and configuring the pipeline correctly. +ASPEN generates three related-but-distinct peak representations, not three +competing choices: `pooled.narrowPeak` (diagnostic — how many peaks exist in +the pooled replicate data, before requiring individual-replicate +reproducibility), `consensus.bed` (Round 1 — variable-width, +reproducibility-filtered peaks per sample, useful for peak annotation and for +comparing against other pipelines' conventional peak outputs), and the +`fixed_width` consensus (Round 2 — same width for every ROI, specifically so +that read/Tn5 counts and DESeq2 differential testing aren't confounded by +peak-length differences or broken by daisy-chained merges across many +samples). **Use the fixed-width ROI set for differential accessibility +testing; use `consensus.bed`/`pooled.narrowPeak` for peak annotation, QC, or +cross-pipeline comparison.** + **Round 1 — Per-sample consensus** (`*.macs2.consensus.bed` / `*.genrich.consensus.bed`) > _"Which peaks are reproducible across biological replicates of the same sample?"_ 1. All replicate tagAlign files for a sample are **pooled** and passed to MACS2/Genrich together → produces a large set of candidate peaks from the pooled data. 2. Each candidate peak is checked for overlap with peaks called in each **individual replicate**. -3. A peak is retained in the consensus only if it overlaps peaks in **≥ `consensus_min_replicates`** replicates (default: 1, meaning present in at least 1 replicate). Adjust `consensus_min_replicates` in `config.yaml` to make this filter stricter. +3. A peak is retained in the consensus only if it overlaps peaks in **≥ `consensus_min_replicates`** replicates (default: 2, meaning present in at least 2 replicates). Adjust `consensus_min_replicates` in `config.yaml` to make this filter stricter. ``` replicate 1 peaks ─┐ @@ -155,8 +196,9 @@ replicate 3 peaks ─┘ (≥ min_replicates) > _"Across ALL samples, what is the unified set of fixed-width windows for differential accessibility analysis?"_ 1. Per-sample consensus peaks are converted to **fixed-width windows** (default: 500 bp, centered on the peak summit) using the method from [Corces et al. 2018](https://doi.org/10.1038/nmeth.4396). Window width is controlled by `fixed_width` in `config.yaml`. -2. Fixed-width peaks from all samples are **merged** → `ROI.macs2.bed` — the master set of regions used for DESeq2 read counting and differential accessibility analysis. -3. P-values across the fixed-width consensus are re-normalized (`*.renormalized.fixed_width.consensus.narrowPeak`) to account for the pooling. +2. Each fixed-width peak is additionally required to overlap peaks in **≥ `roi_min_replicates`** samples/replicates (default: 1) and to meet the **`roi_min_spm`** signal-per-million threshold (default: 2) — analogous to, but independent from, the Round 1 `consensus_min_*` filters above. Adjust `roi_min_replicates`/`roi_min_spm` in `config.yaml` to change how strict this cross-sample filter is. +3. Fixed-width peaks passing that filter are **merged** → `ROI.macs2.bed` — the master set of regions used for DESeq2 read counting and differential accessibility analysis. +4. P-values across the fixed-width consensus are re-normalized (`*.renormalized.fixed_width.consensus.narrowPeak`) to account for the pooling. ``` sample1.consensus.bed ─► fixed-width peaks ─┐ @@ -166,11 +208,27 @@ sample3.consensus.bed ─► fixed-width peaks ─┘ !!! tip "Config knobs that control consensus" -| Parameter | Default | Effect | -| -------------------------- | ------- | ---------------------------------------------------------------------------- | -| `consensus_min_replicates` | `1` | Min. replicates a peak must appear in to be retained in per-sample consensus | -| `consensus_min_spm` | `5` | Min. signal-per-million reads threshold for a peak to be included | -| `fixed_width` | `500` | Width (bp) of fixed-width peaks used to build the ROI set | +| Parameter | Default | Round | Effect | +| -------------------------- | ------- | ----- | ---------------------------------------------------------------------------- | +| `consensus_min_replicates` | `2` | 1 | Min. replicates a peak must appear in to be retained in per-sample consensus | +| `consensus_min_spm` | `5` | 1 | Min. signal-per-million reads threshold for a peak to be included | +| `roi_min_replicates` | `1` | 2 | Min. samples/replicates a fixed-width peak must appear in to be kept in the ROI set | +| `roi_min_spm` | `2` | 2 | Min. signal-per-million reads threshold for a fixed-width peak to be kept in the ROI set | +| `fixed_width` | `500` | 1 | Width (bp) of fixed-width peaks used to build the ROI set | + +#### Counts matrices: reads vs Tn5 nicking sites, and `dedup` vs `nondedup` + +Once the ROIs are established, ASPEN generates two distinct count matrices: + +- **Tn5 Nicking Sites Count Matrix**: This matrix quantifies the frequency of Tn5 transposase insertion events at each ROI. The number of insertion events serves as a proxy for chromatin accessibility, since Tn5 transposase preferentially inserts into accessible regions of the chromatin. By counting these insertion sites, researchers can accurately infer the openness of chromatin regions under different experimental conditions. + +- **Read Counts Matrix**: This matrix records the number of sequencing reads mapped to each ROI. While Tn5 nicking sites provide a more direct measure of chromatin accessibility, read counts are included in ASPEN as they have been widely used in recent publications. Analyzing both matrices together offers a comprehensive view of chromatin accessibility dynamics. + +Both count matrices (and the corresponding DiffATAC/DESeq2 results) are generated **twice** — once from `dedup.bam` (PCR/optical duplicates removed) and once from `filtered.bam` (quality-filtered, duplicates retained, labeled `nondedup`) — under separate `dedup`/`nondedup` output folders (see the folder tree above). Previously only the duplicate-retaining `filtered.bam` was used, which was inconsistent with how MACS2 peaks are called (from `dedup.bam`-derived `tagAlign.gz`) and how FRiP is computed (also from `dedup.bam`). + +**Recommendation: use the `dedup` counts/DiffATAC results for standard differential accessibility analysis.** PCR duplication rate varies across samples (input amount, PCR cycles, library complexity), and this variability is not corrected by DESeq2/edgeR size factors — leaving duplicates in can confound differential calls. This also matches standard ATAC-seq practice (e.g. the ENCODE ATAC-seq pipeline, ArchR, Signac). + +The `nondedup` outputs are retained for comparison/legacy reasons, with one caveat worth knowing: because Tn5 preferentially inserts into a limited set of highly accessible positions, independent DNA molecules can genuinely land on the same insertion coordinate at very open/narrow sites, and naive deduplication can occasionally discard some true (non-PCR) signal there. Picard's paired-end duplicate definition (which requires *both* fragment ends to match, not just one nick site) substantially — though not perfectly — mitigates this. If you suspect this is affecting your data, cross-check the per-replicate NRF/PBC (`QC/preseq`) library-complexity metrics: unexpectedly high duplication despite high library complexity points to this effect rather than true PCR over-amplification. #### Peak Annotation folder @@ -237,18 +295,32 @@ WORKDIR ├── sample1.consensus.macs2.peakfiles ├── sample1.replicate.macs2.peakfiles ├── DiffATAC - │   ├── reads - │   │   ├── all_diff_atacs.html - │   │   ├── all_diff_atacs.tsv - │   │   ├── degs.done - │   │   ├── sample2_vs_sample1.html - │   │   └── sample2_vs_sample1.tsv - │   └── tn5sites - │   ├── all_diff_atacs.html - │   ├── all_diff_atacs.tsv - │   ├── degs.done - │   ├── sample2_vs_sample1.html - │   └── sample2_vs_sample1.tsv + │   ├── dedup + │   │   ├── reads + │   │   │   ├── all_diff_atacs.html + │   │   │   ├── all_diff_atacs.tsv + │   │   │   ├── degs.done + │   │   │   ├── sample2_vs_sample1.html + │   │   │   └── sample2_vs_sample1.tsv + │   │   └── tn5sites + │   │   ├── all_diff_atacs.html + │   │   ├── all_diff_atacs.tsv + │   │   ├── degs.done + │   │   ├── sample2_vs_sample1.html + │   │   └── sample2_vs_sample1.tsv + │   └── nondedup + │   ├── reads + │   │   ├── all_diff_atacs.html + │   │   ├── all_diff_atacs.tsv + │   │   ├── degs.done + │   │   ├── sample2_vs_sample1.html + │   │   └── sample2_vs_sample1.tsv + │   └── tn5sites + │   ├── all_diff_atacs.html + │   ├── all_diff_atacs.tsv + │   ├── degs.done + │   ├── sample2_vs_sample1.html + │   └── sample2_vs_sample1.tsv ├── sample2 │   ├── sample2_replicate1.macs2.narrowPeak │   ├── sample2_replicate1.macs2.narrowPeak_motif_enrichment @@ -294,8 +366,12 @@ WORKDIR ├── sample1.renormalized.fixed_width.consensus.narrowPeak ├── sample1.renormalized.fixed_width.consensus.narrowPeak.annotated.gz ├── counts - │   ├── ROI.macs2.reads_counts.tsv - │   └── ROI.macs2.tn5sites_counts.tsv + │   ├── dedup + │   │   ├── ROI.macs2.reads_counts.tsv + │   │   └── ROI.macs2.tn5sites_counts.tsv + │   └── nondedup + │   ├── ROI.macs2.reads_counts.tsv + │   └── ROI.macs2.tn5sites_counts.tsv ├── sample2_replicate1.macs2.fixed_width.narrowPeak ├── sample2_replicate2.macs2.fixed_width.narrowPeak ├── sample2.fixed_width.consensus.narrowPeak @@ -322,13 +398,13 @@ Some of the key output files are: | `*.narrowPeak_motif_enrichment/ame_results.txt` | motif enrichment results from AME tool from MEME suite using HOCOMOCO v11 database | | `*.narrowPeak_motif_enrichment/knownResults.txt` | motif enrichment results using HOMER with HOCOMOCO v11 database | | `*.macs2.consensus.bed` | consensus peak call between multiple replicates of each sample. **Note:** consensus bed annotations are located in `QC/peak_annotations` | -| `DiffATAC/reads` | folder containing differential open chromatin results:
- computated using read counts in MACS2 regions of interest (ROIs)
- `all_diff_atacs.html` HTML report aggregated across all contrasts from `contrasts.tsv`
- `all_diff_atacs.tsv` DESeq2 results in TSV format aggregated across all contrasts from `contrasts.tsv`
- HTML and TSV file each per contrast in `contrasts.tsv` | -| `DiffATAC/tn5sites` | folder containing differential open chromatin results:
- computated using Tn5 nicking site counts in MACS2 regions of interest (ROIs)
- `all_diff_atacs.html` HTML report aggregated across all contrasts from `contrasts.tsv`
- `all_diff_atacs.tsv` DESeq2 results in TSV format aggregated across all contrasts from `contrasts.tsv`
- HTML and TSV file each per contrast in `contrasts.tsv` | +| `DiffATAC/dedup/reads`, `DiffATAC/nondedup/reads` | folder containing differential open chromatin results:
- computated using read counts in MACS2 regions of interest (ROIs)
- `dedup` uses `dedup.bam`-derived counts (PCR/optical duplicates removed, recommended); `nondedup` uses `filtered.bam`-derived counts (duplicates retained)
- `all_diff_atacs.html` HTML report aggregated across all contrasts from `contrasts.tsv`
- `all_diff_atacs.tsv` DESeq2 results in TSV format aggregated across all contrasts from `contrasts.tsv`
- HTML and TSV file each per contrast in `contrasts.tsv` | +| `DiffATAC/dedup/tn5sites`, `DiffATAC/nondedup/tn5sites` | folder containing differential open chromatin results:
- computated using Tn5 nicking site counts in MACS2 regions of interest (ROIs)
- `dedup` uses `dedup.bam`-derived counts (PCR/optical duplicates removed, recommended); `nondedup` uses `filtered.bam`-derived counts (duplicates retained)
- `all_diff_atacs.html` HTML report aggregated across all contrasts from `contrasts.tsv`
- `all_diff_atacs.tsv` DESeq2 results in TSV format aggregated across all contrasts from `contrasts.tsv`
- HTML and TSV file each per contrast in `contrasts.tsv` | | `fixed_width` | `fixed_width` can be set in `config.yaml` to create peaks of a user defined fixed width (default 500bp). This folder contains:
- individual replicate `*.fixed_width.narrowPeak` files
- `*.renormalized.fixed_width.consensus.narrowPeak` per sample; [_Corces et. al._](https://doi.org/10.1038/nmeth.4396) method is used for consensus calling; used to generate MACS2 regions of interest (ROI) peaks which are used to generate a reads or Tn5 sites counts matrix for DESeq2
- ROI related files: `ROI.macs2.bed`, `ROI.macs2.bed.annotated.gz`, `ROI.macs2.annotation_summary`, `ROI.macs2.annotation_distribution` | -| `fixed_width/counts/ROI.macs2.read_counts.tsv` | read counts in MACS2 ROIs using featureCounts | -| `fixed_width/counts/ROI.reads_scaled_counts.tsv` | `ROI.macs2.read_counts.tsv` scaled using spike-in scaling factors | -| `fixed_width/counts/ROI.tn5sites_counts.tsv` | Tn5 nicking site counts in MACS2 ROIs using featureCounts | -| `fixed_width/counts/ROI.tn5sites_scaled_counts.tsv` | `ROI.macs2.tn5sites_counts.tsv` scaled using spike-in scaling factors | +| `fixed_width/counts/{dedup,nondedup}/ROI.macs2.read_counts.tsv` | read counts in MACS2 ROIs using featureCounts; generated separately from `dedup.bam` (duplicates removed, recommended) and `filtered.bam`/`nondedup` (duplicates retained) | +| `fixed_width/counts/{dedup,nondedup}/ROI.reads_scaled_counts.tsv` | `ROI.macs2.read_counts.tsv` scaled using spike-in scaling factors | +| `fixed_width/counts/{dedup,nondedup}/ROI.tn5sites_counts.tsv` | Tn5 nicking site counts in MACS2 ROIs using featureCounts; generated separately from `dedup.bam` (duplicates removed, recommended) and `filtered.bam`/`nondedup` (duplicates retained) | +| `fixed_width/counts/{dedup,nondedup}/ROI.tn5sites_scaled_counts.tsv` | `ROI.macs2.tn5sites_counts.tsv` scaled using spike-in scaling factors | ### Genrich output folder diff --git a/docs/overview.md b/docs/overview.md index 0b02b83..545a630 100644 --- a/docs/overview.md +++ b/docs/overview.md @@ -22,15 +22,22 @@ Employs FastQC for initial quality checks and MultiQC for aggregated reporting, ### 🚫 **Blacklist Filtering** -Reads aligning to known blacklisted regions and mitochondrial DNA are systematically removed during the preprocessing stage. Blacklisted regions, as defined by ENCODE and other genomic consortia, represent genomic loci prone to artifacts due to their repetitive nature or unusually high signal, which can confound downstream analyses. Similarly, mitochondrial reads, which often constitute a significant proportion of sequencing data, are excluded to prevent skewing results and to ensure computational resources are focused on nuclear chromatin. This filtering step enhances the accuracy and reliability of peak detection by concentrating on biologically meaningful and interpretable regions of the genome. +Two categories of reads are systematically removed during the preprocessing stage, since both would otherwise confound peak detection: + +| Filtered out | Why | Source of coordinates | +|---|---|---| +| Blacklisted regions | Repetitive/artifact-prone genomic loci with unusually high, spurious signal | ENCODE (and other genomic consortia) blacklists, bundled per-genome in `resources/blacklistFa/` | +| Mitochondrial reads (`chrM`/`chrR`) | Often make up a large fraction of ATAC-seq reads (mitochondrial DNA is highly accessible); excluding them focuses compute and peak calling on nuclear chromatin | The reference genome's `chrM`/`chrR` contig | + +This filtering step enhances the accuracy and reliability of peak detection by concentrating on biologically meaningful and interpretable regions of the genome. ### 🔗 **Read Alignment** Reads are aligned to the reference genome using Bowtie2, a fast and memory-efficient aligner optimized for short-read sequencing data. This step produces several key output files essential for subsequent analyses: -- **`filtered.bam`**: A BAM file with low-quality alignments, reads mapping to blacklisted regions removed, ensuring high-quality data for differential open chromatin analysis. +- **`filtered.bam`**: A BAM file with low-quality alignments and reads mapping to blacklisted regions removed. This is the **`nondedup`** BAM referenced elsewhere in these docs — i.e. the file *before* PCR/optical duplicate removal — and is used for the `nondedup` variant of the counts matrices and DiffATAC results. - **`qsorted.bam`**: A quality-sorted BAM file containing all aligned reads, serving as the primary input for Genrich peak calling. -- **`dedup.bam`**: A BAM file with PCR duplicates removed, used for advanced downstream analyses such as transcription factor footprinting with TOBIAS. +- **`dedup.bam`**: A BAM file derived from `filtered.bam` with PCR/optical duplicates additionally removed via PicardTools. Used for MACS2 peak calling, FRiP, transcription factor footprinting with TOBIAS, and the **recommended `dedup`** variant of the counts matrices and DiffATAC results. - **`tagAlign.gz`**: A compressed file containing alignment information in a simplified format, specifically prepared for MACS2 peak calling. This alignment process ensures that only high-confidence reads are retained, providing a robust foundation for accurate and reproducible ATAC-seq data analysis. @@ -51,17 +58,36 @@ Genrich is integrated into the pipeline to complement MACS2. This tool is partic By combining the strengths of MACS2 and Genrich, ASPEN delivers a comprehensive and reliable peak detection framework, facilitating downstream analyses and enabling researchers to uncover critical insights into chromatin accessibility and gene regulation. -### 🤝 **Consensus Peaks** - -For datasets with multiple replicates, ASPEN employs several strategies to derive consensus peaks across replicates: +!!! tip "Which peak caller should I use?" + Both MACS2 and Genrich run automatically for every sample — you + don't have to choose upfront. **Genrich is generally + recommended for ATAC-seq** because it was purpose-built for + assays like ATAC-seq/DNase-seq: it models Tn5 transposase cut + sites directly (`-j` ATAC-seq mode) rather than adapting a + ChIP-seq fragment-shift model, and it combines biological + replicates natively via Fisher's method instead of requiring a + separate consensus step. **MACS2** was originally designed for + ChIP-seq and requires ATAC-specific parameter workarounds to + approximate cut-site signal, but remains the field standard — + it's included because many reviewers and downstream tools + expect to see MACS2 peaks, and comparing both gives an extra + sanity check. CCBR's internal benchmarking on ASPEN's ATAC-seq + data has generally found Genrich peaks to be higher quality. If + your MACS2 and Genrich DiffATAC results disagree substantially + for a given region, treat that region as lower-confidence + rather than assuming one caller is unconditionally "right". -- **`consensus.bed`**: This file represents consensus peaks derived using the "Consensus MAX" strategy, as described by [_Yang et al._](https://doi.org/10.5936/csbj.201401002). It identifies overlapping peaks across replicates to define a shared set of regions. +### 🤝 **Consensus Peaks** -- **`pooled.narrowPeak`**: In this approach, reads from all replicates are pooled together, and peak calling is performed on the combined dataset to generate a unified set of "pooled" peaks. +For datasets with multiple replicates, ASPEN generates three related-but-distinct peak representations — not three competing choices, but two sequential stages plus one diagnostic side-output: -- **`fixed_width.consensus.narrowPeak`**: This method calculates consensus peaks by renormalizing p-value scores across replicates, following the approach outlined by [_Corces et al._](https://doi.org/10.1038/nmeth.4396). +| Output | Stage | Width | Best used for | +|---|---|---|---| +| `pooled.narrowPeak` | Diagnostic (pre-filter) — peaks called on all replicates' reads pooled together | Variable | Sanity-checking how strict the reproducibility filter (`consensus_min_replicates`/`consensus_min_spm`) is | +| `consensus.bed` | Round 1 — per-sample, reproducibility-filtered via the ["Consensus MAX" strategy](https://doi.org/10.5936/csbj.201401002) (_Yang et al._) | Variable | Peak annotation (ChIPseeker) and comparison against other pipelines' conventional peak outputs | +| `fixed_width.consensus.narrowPeak` | Round 2 — cross-sample ROI set, p-values renormalized across replicates following [_Corces et al._](https://doi.org/10.1038/nmeth.4396) | Fixed (`fixed_width` config, default 500 bp) | Differential accessibility testing (DESeq2) — fixed width avoids confounding read/Tn5 counts with peak-length differences and avoids daisy-chained merges across many samples | -Among these methods, the fixed-width consensus peak strategy is recommended for its robustness and reproducibility. The other outputs are provided primarily for exploratory purposes. +**Use the fixed-width ROI set for differential accessibility testing; use `consensus.bed`/`pooled.narrowPeak` for peak annotation, QC, or cross-pipeline comparison.** ### 🏷️ **Peak Annotation** @@ -71,7 +97,7 @@ ASPEN integrates ChIPseeker to perform comprehensive annotation of identified pe Regions of Interest (ROIs) are generated through a systematic process to ensure consistency and reliability in downstream analyses: -- 🧩 **Replicate-Level Peaks**: `.narrowPeak` files from individual replicates are converted into fixed-width peaks, centering on their summits and extending ±250 bp to achieve a uniform width of 500 bp. +- 🧩 **Replicate-Level Peaks**: `.narrowPeak` files from individual replicates are converted into fixed-width peaks, centered on their summits and extended symmetrically to the width configured by `fixed_width` in `config.yaml` (default: 500 bp, i.e. ±250 bp from the summit). - 🧬 **Sample-Level Consensus Peaks**: Fixed-width replicate peak files are combined to generate a consensus set of peaks at the sample level, ensuring reproducibility across replicates. - 🔗 **Merged ROIs**: Sample-level consensus peak files are merged using the [Coerses method](https://doi.org/10.1038/nmeth.4396) to create a comprehensive set of ROIs. @@ -87,16 +113,25 @@ ASPEN employs custom scripts to analyze the distribution of fragment lengths wit ### 📚 **Library Complexity** -To evaluate the sufficiency of sequencing depth and detect potential biases introduced during PCR amplification, ASPEN utilizes Preseq to estimate library complexity. This metric helps determine whether the sequencing effort is adequate to capture the diversity of the library, ensuring that the data is representative of the underlying chromatin landscape. By identifying potential saturation or over-representation of certain fragments, researchers can assess the reliability of their sequencing results. +To evaluate the sufficiency of sequencing depth and detect potential biases introduced during PCR amplification, ASPEN utilizes Preseq to estimate library complexity, reporting the Non-Redundant Fraction (`NRF`) and PCR Bottlenecking Coefficients (`PBC1`, `PBC2`). This metric helps determine whether the sequencing effort is adequate to capture the diversity of the library, ensuring that the data is representative of the underlying chromatin landscape. By identifying potential saturation or over-representation of certain fragments, researchers can assess the reliability of their sequencing results. + +!!! tip "Rule of thumb" + Per [ENCODE's ATAC-seq data standards](https://www.encodeproject.org/atac-seq/#standards), the preferred values are **`NRF` > 0.9, `PBC1` > 0.9, and `PBC2` > 3**. Lower values indicate a less complex library (e.g. over-amplified by PCR), which can inflate apparent signal at a subset of loci rather than reflecting true biological accessibility. ### 🧬 **Transcription Start Site (TSS) Enrichment** ASPEN calculates TSS enrichment scores, a widely recognized quality metric for ATAC-seq data. These scores measure the accumulation of sequencing reads around transcription start sites (TSS), which are hallmark regions of open chromatin. High TSS enrichment scores indicate well-prepared libraries with minimal technical artifacts, as they reflect the accessibility of promoter regions and the integrity of the chromatin preparation process. +!!! tip "Rule of thumb" + [ENCODE's ATAC-seq data standards](https://www.encodeproject.org/atac-seq/#standards) define annotation-dependent TSS enrichment cutoffs — for example, using a GRCh38 RefSeq TSS annotation: **< 5 is concerning, 5-7 is acceptable, and > 7 is ideal**. **Caveat:** ASPEN builds its TSS bins from GENCODE (not RefSeq) gene annotations (see `resources/tssBed/`), so ENCODE's exact per-annotation cutoffs may not transfer precisely to ASPEN's TSS enrichment values — treat these numbers as directional guidance (aim for high single digits or higher) rather than an exact pass/fail threshold. + ### 📊 **Fraction of Reads in Peaks (FRiP)** The Fraction of Reads in Peaks (FRiP) score quantifies the proportion of sequencing reads that fall within identified peaks, serving as a measure of the signal-to-noise ratio in the dataset. Higher FRiP scores indicate datasets with strong, biologically meaningful signals and minimal background noise. Additionally, ASPEN computes the fraction of reads localized to specific genomic features, such as promoters, enhancers, and DNase hypersensitive sites (DHS). These feature-specific FRiP scores provide further insights into the quality and biological relevance of the data. +!!! tip "Rule of thumb" + Per [ENCODE's ATAC-seq data standards](https://www.encodeproject.org/atac-seq/#standards), a FRiP score **> 0.3** indicates high-quality data, though values **> 0.2** may still be acceptable. Consistently lower scores suggest poor signal-to-noise and warrant a closer look at library prep or peak-calling parameters. + --- ## 🧬 **Motif Enrichment Analysis** @@ -117,7 +152,7 @@ ASPEN extends its core functionalities by offering advanced capabilities for dif ### 🛠️ **Refinement of Peak Calls to Fixed-Width Peaks** -Following the initial peak calling, ASPEN standardizes the identified peaks by centering them on their summits and extending them symmetrically to a uniform width of 500 base pairs (bp). This refinement ensures consistency across datasets and facilitates downstream analyses. The approach of extending peak summits by ±250 bp is a widely accepted practice in ATAC-seq data analysis, as it focuses on the most accessible regions of chromatin while minimizing variability in peak sizes. This standardization enhances the comparability of results across experimental conditions and replicates. +Following the initial peak calling, ASPEN standardizes the identified peaks by centering them on their summits and extending them symmetrically to the width configured by `fixed_width` in `config.yaml` (default: 500 bp, i.e. ±250 bp from the summit). This refinement ensures consistency across datasets and facilitates downstream analyses. Extending peak summits to a fixed width is a widely accepted practice in ATAC-seq data analysis, as it focuses on the most accessible regions of chromatin while minimizing variability in peak sizes. This standardization enhances the comparability of results across experimental conditions and replicates. ### 📍 **Definition of Regions of Interest (ROIs)** @@ -127,9 +162,11 @@ ASPEN processes the fixed-width peak calls for each replicate across all samples Once the ROIs are established, ASPEN generates two key count matrices to quantify chromatin accessibility: -- **Tn5 Nicking Sites Count Matrix**: This matrix quantifies the frequency of Tn5 transposase insertion events at each ROI. The `filtered.bam` file is used to determine Tn5 insertion sites and count them. Tn5 transposase preferentially inserts into accessible chromatin regions, and the number of insertion events serves as a direct proxy for chromatin accessibility. This metric provides a precise measure of chromatin openness across experimental conditions. +- **Tn5 Nicking Sites Count Matrix**: This matrix quantifies the frequency of Tn5 transposase insertion events at each ROI. Tn5 transposase preferentially inserts into accessible chromatin regions, and the number of insertion events serves as a direct proxy for chromatin accessibility. + +- **Read Counts Matrix**: This matrix records the number of sequencing reads mapped to each ROI. While Tn5 nicking sites offer a more direct measure of chromatin accessibility, read counts are included as they are widely utilized in recent publications. -- **Read Counts Matrix**: This matrix records the number of sequencing reads mapped to each ROI. Similar to the Tn5 nicking sites matrix, the `filtered.bam` file is used for this quantification. While Tn5 nicking sites offer a more direct measure of chromatin accessibility, read counts are included as they are widely utilized in recent publications. Together, these matrices provide a comprehensive view of chromatin accessibility dynamics. +Both matrices (and the corresponding DiffATAC/DESeq2 results) are generated **twice** — once from `dedup.bam` (PCR/optical duplicates removed) and once from `filtered.bam` (duplicates retained, labeled `nondedup`) — under separate `dedup`/`nondedup` output folders. **The `dedup` results are recommended** for standard differential accessibility testing, since PCR duplication rate varies across samples and is not corrected by DESeq2/edgeR size factors. See [Counts matrices: reads vs Tn5 nicking sites, and `dedup` vs `nondedup`](outputs.md#counts-matrices-reads-vs-tn5-nicking-sites-and-dedup-vs-nondedup) in the Outputs page for the full explanation, including when the `nondedup` outputs are still useful. ### 📈 **Differential Accessibility Analysis Using DESeq2** @@ -143,9 +180,7 @@ To enhance the biological interpretation of differential accessibility results, #### 🧪 **Optional: Spike-In Normalization** -To account for potential global shifts in chromatin accessibility—particularly in perturbation studies where widespread chromatin compaction or relaxation may occur—ASPEN optionally supports **spike-in normalization**. - -Spike-in normalization involves the use of exogenous DNA or cells from a different species (e.g., _Drosophila melanogaster_ or _E. coli_) that are added in equal proportions across all experimental samples prior to lysis and tagmentation. These spike-in reads serve as an internal control to correct for technical variation and global accessibility shifts that may not be captured by traditional normalization strategies. +As discussed in [Challenges in ATAC-seq Data Analysis](introduction.md#challenges-in-atac-seq-data-analysis), perturbation studies can cause global shifts in chromatin accessibility that traditional normalization can't correct for. ASPEN addresses this with optional **spike-in normalization**: exogenous DNA or cells from a different species (e.g., _Drosophila melanogaster_ or _E. coli_) are added in equal proportions across all experimental samples prior to lysis and tagmentation, serving as an internal control for technical variation and global accessibility shifts. In ASPEN, if spike-in data is present: @@ -153,14 +188,21 @@ In ASPEN, if spike-in data is present: - The number of spike-in reads is counted for each sample. - A normalization factor is calculated based on spike-in counts and applied to the accessibility read counts from the host genome. -This spike-in-derived scaling factor allows the comparison of chromatin accessibility across conditions even when global chromatin accessibility levels differ (e.g., treatment-induced repression or global decondensation). This method is particularly valuable in experiments involving: +This spike-in-derived scaling factor allows the comparison of chromatin accessibility across conditions even when global chromatin accessibility levels differ (e.g., treatment-induced repression or global decondensation). -- Transcription factor knockdowns/knockouts -- Chromatin remodeler inhibition -- Drug-induced chromatin modulation +!!! tip "Should I turn on spike-in normalization?" + Ask yourself: do I expect a **global, genome-wide shift** in chromatin accessibility between my conditions — rather than just **localized** changes at a handful of specific regulatory elements? + + - **Yes** (e.g. a chromatin remodeler inhibitor, a broad transcription factor knockdown/knockout, drug-induced chromatin modulation) → turn on spike-in normalization. Standard depth-based normalization (DESeq2 size factors) assumes *most* regions are unchanged between conditions — that assumption breaks down under a genome-wide shift, and spike-in gives you an external, biology-independent scale instead. + - **No** (you expect differences to be confined to specific loci/pathways, with most of the genome unchanged) → spike-in is probably unnecessary. It adds experimental complexity (extra reagents, a second alignment step, and a "0 spike-in reads" failure mode to manage — see the warning below) for a scenario DESeq2's built-in normalization already handles well. + + See [Enabling Spike-In Normalization](deployment.md#enabling-spike-in-normalization-optional) for the config steps once you've decided to use it. ASPEN performs spike-in-aware normalization transparently, and reports both raw and normalized counts in the final output matrix for differential analysis. This ensures flexibility in downstream interpretation while preserving the ability to adjust for systemic experimental artifacts. +!!! warning "What if a replicate has 0 spike-in reads?" +If `spikein: true` is set but a replicate has **zero reads** aligned to the spike-in genome (e.g. the spike-in material wasn't actually added to that library, the spike-in genome/index is misconfigured, or a genuinely contamination-free host-only library), ASPEN cannot compute a scaling factor for it and the run will **fail with a clear error message** naming the affected replicate(s) rather than silently producing `Inf`/`NaN` normalized counts. To resolve this: verify the spike-in genome/index path in `config.yaml`, confirm spike-in material was actually included during library prep for that replicate, or set `spikein: false` if none of your samples have spike-in material. + ### 📊 **Reporting** ASPEN enhances differential chromatin accessibility analysis by providing an interactive HTML report generated from DESeq2 results, featuring various visualizations. Additionally, it offers a TSV (tabl-delimited) file with integrated gene annotations, compatible with Microsoft Excel, enabling efficient data manipulation and facilitating the identification of genes near regions with altered accessibility. The Excel file is aggregated across all different contrasts queried in the project. From 58469c9715ea3e5519d447965550a0220a774d03 Mon Sep 17 00:00:00 2001 From: kopardev Date: Fri, 11 Sep 2026 13:55:43 -0400 Subject: [PATCH 5/6] chore(changelog): use 'development version' heading pre-release MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The 1.3.0 heading was premature -- per repo convention the CHANGELOG heading stays as 'development version' while VERSION is a -dev suffix, and only gets renamed to the actual release version (e.g. 1.3.0) by the draft-release automation when preparing the release. ⚡ Generated using AI ⚡ --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9a3bd7a..1961861 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,4 +1,4 @@ -## ASPEN 1.3.0 +## ASPEN development version - Generate Tn5 nicking-site and read counts matrices, and the corresponding DiffATAC/DESeq2 results, from both `dedup.bam` (PCR/optical duplicates removed) and `filtered.bam` (duplicates retained, labeled `nondedup`), written to separate `dedup`/`nondedup` output subfolders under `visualization/` and `peaks/{peakcaller}/{fixed_width/counts,DiffATAC}/`; previously only the duplicate-retaining `filtered.bam` was used. `dedup` is recommended for standard differential accessibility testing. (#138, @kopardev) - Fail fast with a clear, actionable error message when a replicate has 0 reads aligned to the spike-in genome, instead of crashing with an opaque `ZeroDivisionError` in `_compute_downsampling_scaling_factors.py`; the error is now also captured in a dedicated `results/spikein/compute_scaling_factors.log`. (#139, @kopardev) From 1b06ebb1a69f3635b3fba955a892b73b106fb3e8 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 11 Sep 2026 17:58:57 +0000 Subject: [PATCH 6/6] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- docs/index.md | 6 ++-- docs/introduction.md | 14 +++++----- docs/outputs.md | 66 ++++++++++++++++++++++---------------------- docs/overview.md | 62 ++++++++++++++++++++--------------------- 4 files changed, 74 insertions(+), 74 deletions(-) diff --git a/docs/index.md b/docs/index.md index 978d5e0..67714c2 100644 --- a/docs/index.md +++ b/docs/index.md @@ -1,8 +1,8 @@ ## Background !!! tip "New here?" - If you just want to **run the pipeline**, skip ahead to - [Running ASPEN](deployment.md). Come back here later for the - scientific background on ATAC-seq and why ASPEN is built the way it is. +If you just want to **run the pipeline**, skip ahead to +[Running ASPEN](deployment.md). Come back here later for the +scientific background on ATAC-seq and why ASPEN is built the way it is. The Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) has revolutionized genomics by providing a rapid and sensitive method to assess chromatin accessibility across the genome. This technique offers profound insights into gene regulation, epigenetic modifications, and the dynamic landscape of the chromatin environment. To facilitate the comprehensive analysis of ATAC-seq data, the Center for Cancer Research (CCR) Collaborative Bioinformatics Resource (CCBR) has developed ASPEN (**A**tac **S**eq **P**ip**E**li**N**e), an automated, robust, reproducible pipeline designed using the Snakemake pipelining framework to streamline the complexities inherent in ATAC-seq data processing. diff --git a/docs/introduction.md b/docs/introduction.md index eb003f7..47c707e 100644 --- a/docs/introduction.md +++ b/docs/introduction.md @@ -10,10 +10,10 @@ ATAC-seq is a powerful technique that enables the identification of open chromat Despite its advantages, ATAC-seq data analysis presents several challenges: -| Challenge | Why it's hard | How ASPEN helps | -|---|---|---| -| 📊 **Data complexity** | Datasets contain a mix of fragment sizes corresponding to nucleosome-free regions, mono-nucleosomes, and multi-nucleosomes, requiring careful computational interpretation | Reports a [fragment length distribution](overview.md#fragment-length-distribution) QC metric so you can visually confirm the expected nucleosome-free/mono-nucleosome pattern | -| ✅ **Quality control** | Requires assessing sequencing quality, fragment size distribution, and enrichment of reads in regulatory regions all together | Runs FastQC/MultiQC plus ATAC-specific metrics ([FRiP](overview.md#fraction-of-reads-in-peaks-frip), [TSS enrichment](overview.md#transcription-start-site-tss-enrichment), [library complexity](overview.md#library-complexity)) automatically for every sample | -| 📈 **Peak calling** | Distinguishing true accessible regions from background noise is especially hard in low signal-to-noise datasets | Runs both MACS2 and Genrich (see [Peak Detection](overview.md#peak-detection)) so results can be cross-checked against each other | -| 🔬 **Spike-in normalization** | Variability in spike-in recovery or accessibility differences between samples can complicate normalization | Provides optional, transparent [spike-in normalization](overview.md#optional-spike-in-normalization) with a documented failure mode if spike-in reads are missing | -| 🔄 **Reproducibility** | Achieving consistent results across experiments and conditions requires more than just consistent code | Runs entirely in versioned Docker/Singularity containers (see "Containerization" above) so the same environment is reproduced on every run | +| Challenge | Why it's hard | How ASPEN helps | +| ----------------------------- | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| 📊 **Data complexity** | Datasets contain a mix of fragment sizes corresponding to nucleosome-free regions, mono-nucleosomes, and multi-nucleosomes, requiring careful computational interpretation | Reports a [fragment length distribution](overview.md#fragment-length-distribution) QC metric so you can visually confirm the expected nucleosome-free/mono-nucleosome pattern | +| ✅ **Quality control** | Requires assessing sequencing quality, fragment size distribution, and enrichment of reads in regulatory regions all together | Runs FastQC/MultiQC plus ATAC-specific metrics ([FRiP](overview.md#fraction-of-reads-in-peaks-frip), [TSS enrichment](overview.md#transcription-start-site-tss-enrichment), [library complexity](overview.md#library-complexity)) automatically for every sample | +| 📈 **Peak calling** | Distinguishing true accessible regions from background noise is especially hard in low signal-to-noise datasets | Runs both MACS2 and Genrich (see [Peak Detection](overview.md#peak-detection)) so results can be cross-checked against each other | +| 🔬 **Spike-in normalization** | Variability in spike-in recovery or accessibility differences between samples can complicate normalization | Provides optional, transparent [spike-in normalization](overview.md#optional-spike-in-normalization) with a documented failure mode if spike-in reads are missing | +| 🔄 **Reproducibility** | Achieving consistent results across experiments and conditions requires more than just consistent code | Runs entirely in versioned Docker/Singularity containers (see "Containerization" above) so the same environment is reproduced on every run | diff --git a/docs/outputs.md b/docs/outputs.md index 8822ee4..5432529 100644 --- a/docs/outputs.md +++ b/docs/outputs.md @@ -126,18 +126,18 @@ WORKDIR Content details: -| Folder | SubFolder | Description | -| ------------- | ------------------- | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | -| alignment | qsortedBam | - Query sorted Bowtie2 alignments in BAM format.
- Excludes unmapped and platform/vendor quality failing reads.
- Used for Genrich peak calling. | -| alignment | filteredBam | - Filtered BAM files after excluding non-primary, supplementary, and MAPQ <=5 alignments.
- This is the `nondedup` BAM (i.e. before duplicate removal) used for the `nondedup` counts/DiffATAC variant.
- Derived from `qsortedBam`. | -| alignment | dedupBam | - Deduplicated filtered BAM files (PCR/optical duplicates marked with PicardTools and excluded).
- Derived from `filteredBam`.
- Used for MACS2 peak calling, FRiP, and the recommended `dedup` counts/DiffATAC variant.
- Can also be used downstream with CCBR_TOBIAS pipeline. | -| alignment | tagAlign | - `tagAlign.gz` files used for MACS2 peak calling.
- Derived from `dedupBam`. | -| peaks | genrich & macs | - Genrich/MACS2 peak calls (raw, consensus, fixed-width).
- Contains ROI files with Diff-ATAC results if `contrasts.tsv` is provided.
- Calculated with DESeq2 using both read counts and tn5 nicking sites in ROI.
- Counts matrices and DiffATAC results are generated separately for `dedup` (PCR duplicates removed, recommended) and `nondedup` (duplicates retained) BAMs, in like-named subfolders. | -| visualization | dedup / nondedup | - Tn5 nicking-site/read BAMs, BEDs, and bigwigs, split into `dedup` and `nondedup` subfolders depending on whether PCR duplicates were removed from the source BAM before generating them. | -| QC | various | - Flagstats.
- Dupmetrics.
- Read counts.
- Motif enrichments.
- FLD stats.
- Fqscreen.
- FRiP.
- ChIPSeeker results.
- TSS enrichments.
- Preseq.
- Homer/AME motif enrichments.
- MultiQC. | -| QC | peak_annotation | detailed peak annotations described below | -| 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. | +| Folder | SubFolder | Description | +| ------------- | ------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| alignment | qsortedBam | - Query sorted Bowtie2 alignments in BAM format.
- Excludes unmapped and platform/vendor quality failing reads.
- Used for Genrich peak calling. | +| alignment | filteredBam | - Filtered BAM files after excluding non-primary, supplementary, and MAPQ <=5 alignments.
- This is the `nondedup` BAM (i.e. before duplicate removal) used for the `nondedup` counts/DiffATAC variant.
- Derived from `qsortedBam`. | +| alignment | dedupBam | - Deduplicated filtered BAM files (PCR/optical duplicates marked with PicardTools and excluded).
- Derived from `filteredBam`.
- Used for MACS2 peak calling, FRiP, and the recommended `dedup` counts/DiffATAC variant.
- Can also be used downstream with CCBR_TOBIAS pipeline. | +| alignment | tagAlign | - `tagAlign.gz` files used for MACS2 peak calling.
- Derived from `dedupBam`. | +| peaks | genrich & macs | - Genrich/MACS2 peak calls (raw, consensus, fixed-width).
- Contains ROI files with Diff-ATAC results if `contrasts.tsv` is provided.
- Calculated with DESeq2 using both read counts and tn5 nicking sites in ROI.
- Counts matrices and DiffATAC results are generated separately for `dedup` (PCR duplicates removed, recommended) and `nondedup` (duplicates retained) BAMs, in like-named subfolders. | +| visualization | dedup / nondedup | - Tn5 nicking-site/read BAMs, BEDs, and bigwigs, split into `dedup` and `nondedup` subfolders depending on whether PCR duplicates were removed from the source BAM before generating them. | +| QC | various | - Flagstats.
- Dupmetrics.
- Read counts.
- Motif enrichments.
- FLD stats.
- Fqscreen.
- FRiP.
- ChIPSeeker results.
- TSS enrichments.
- Preseq.
- Homer/AME motif enrichments.
- MultiQC. | +| QC | peak_annotation | detailed peak annotations described below | +| 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. | !!! note BAM files from `dedupBam` can be used for downstream footprinting analysis using [CCBR_TOBIAS](https://github.com/CCBR/CCBR_Tobias) pipeline @@ -152,9 +152,9 @@ BAM files from `dedupBam` can also be converted to BED format and processed with **Why two rounds?** A single replicate's peak calls are noisy — some "peaks" are just sequencing/technical artifacts that happened to occur in one -replicate. Pooling and filtering peaks *within* a sample (Round 1) gives a +replicate. Pooling and filtering peaks _within_ a sample (Round 1) gives a reproducibility-checked, per-sample truth set. But differential accessibility -testing needs to compare counts *across* samples/conditions, which requires +testing needs to compare counts _across_ samples/conditions, which requires one shared set of regions measured identically everywhere — a per-sample peak set alone doesn't give you that. Round 2 solves this by merging the per-sample consensus peaks into a single, fixed-width Region of Interest @@ -208,13 +208,13 @@ sample3.consensus.bed ─► fixed-width peaks ─┘ !!! tip "Config knobs that control consensus" -| Parameter | Default | Round | Effect | -| -------------------------- | ------- | ----- | ---------------------------------------------------------------------------- | -| `consensus_min_replicates` | `2` | 1 | Min. replicates a peak must appear in to be retained in per-sample consensus | -| `consensus_min_spm` | `5` | 1 | Min. signal-per-million reads threshold for a peak to be included | -| `roi_min_replicates` | `1` | 2 | Min. samples/replicates a fixed-width peak must appear in to be kept in the ROI set | +| Parameter | Default | Round | Effect | +| -------------------------- | ------- | ----- | ---------------------------------------------------------------------------------------- | +| `consensus_min_replicates` | `2` | 1 | Min. replicates a peak must appear in to be retained in per-sample consensus | +| `consensus_min_spm` | `5` | 1 | Min. signal-per-million reads threshold for a peak to be included | +| `roi_min_replicates` | `1` | 2 | Min. samples/replicates a fixed-width peak must appear in to be kept in the ROI set | | `roi_min_spm` | `2` | 2 | Min. signal-per-million reads threshold for a fixed-width peak to be kept in the ROI set | -| `fixed_width` | `500` | 1 | Width (bp) of fixed-width peaks used to build the ROI set | +| `fixed_width` | `500` | 1 | Width (bp) of fixed-width peaks used to build the ROI set | #### Counts matrices: reads vs Tn5 nicking sites, and `dedup` vs `nondedup` @@ -228,7 +228,7 @@ Both count matrices (and the corresponding DiffATAC/DESeq2 results) are generate **Recommendation: use the `dedup` counts/DiffATAC results for standard differential accessibility analysis.** PCR duplication rate varies across samples (input amount, PCR cycles, library complexity), and this variability is not corrected by DESeq2/edgeR size factors — leaving duplicates in can confound differential calls. This also matches standard ATAC-seq practice (e.g. the ENCODE ATAC-seq pipeline, ArchR, Signac). -The `nondedup` outputs are retained for comparison/legacy reasons, with one caveat worth knowing: because Tn5 preferentially inserts into a limited set of highly accessible positions, independent DNA molecules can genuinely land on the same insertion coordinate at very open/narrow sites, and naive deduplication can occasionally discard some true (non-PCR) signal there. Picard's paired-end duplicate definition (which requires *both* fragment ends to match, not just one nick site) substantially — though not perfectly — mitigates this. If you suspect this is affecting your data, cross-check the per-replicate NRF/PBC (`QC/preseq`) library-complexity metrics: unexpectedly high duplication despite high library complexity points to this effect rather than true PCR over-amplification. +The `nondedup` outputs are retained for comparison/legacy reasons, with one caveat worth knowing: because Tn5 preferentially inserts into a limited set of highly accessible positions, independent DNA molecules can genuinely land on the same insertion coordinate at very open/narrow sites, and naive deduplication can occasionally discard some true (non-PCR) signal there. Picard's paired-end duplicate definition (which requires _both_ fragment ends to match, not just one nick site) substantially — though not perfectly — mitigates this. If you suspect this is affecting your data, cross-check the per-replicate NRF/PBC (`QC/preseq`) library-complexity metrics: unexpectedly high duplication despite high library complexity points to this effect rather than true PCR over-amplification. #### Peak Annotation folder @@ -391,19 +391,19 @@ WORKDIR Some of the key output files are: -| File | Description | -| --------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | -| `*.macs2.narrowPeak` | peak calls from MACS2 filtered by q-value for each samples each replicate | -| `*.macs2.unfiltered.narrowPeak` | peak calls from MACS2 (unfiltered) for each samples each replicate | -| `*.narrowPeak_motif_enrichment/ame_results.txt` | motif enrichment results from AME tool from MEME suite using HOCOMOCO v11 database | -| `*.narrowPeak_motif_enrichment/knownResults.txt` | motif enrichment results using HOMER with HOCOMOCO v11 database | -| `*.macs2.consensus.bed` | consensus peak call between multiple replicates of each sample. **Note:** consensus bed annotations are located in `QC/peak_annotations` | -| `DiffATAC/dedup/reads`, `DiffATAC/nondedup/reads` | folder containing differential open chromatin results:
- computated using read counts in MACS2 regions of interest (ROIs)
- `dedup` uses `dedup.bam`-derived counts (PCR/optical duplicates removed, recommended); `nondedup` uses `filtered.bam`-derived counts (duplicates retained)
- `all_diff_atacs.html` HTML report aggregated across all contrasts from `contrasts.tsv`
- `all_diff_atacs.tsv` DESeq2 results in TSV format aggregated across all contrasts from `contrasts.tsv`
- HTML and TSV file each per contrast in `contrasts.tsv` | -| `DiffATAC/dedup/tn5sites`, `DiffATAC/nondedup/tn5sites` | folder containing differential open chromatin results:
- computated using Tn5 nicking site counts in MACS2 regions of interest (ROIs)
- `dedup` uses `dedup.bam`-derived counts (PCR/optical duplicates removed, recommended); `nondedup` uses `filtered.bam`-derived counts (duplicates retained)
- `all_diff_atacs.html` HTML report aggregated across all contrasts from `contrasts.tsv`
- `all_diff_atacs.tsv` DESeq2 results in TSV format aggregated across all contrasts from `contrasts.tsv`
- HTML and TSV file each per contrast in `contrasts.tsv` | -| `fixed_width` | `fixed_width` can be set in `config.yaml` to create peaks of a user defined fixed width (default 500bp). This folder contains:
- individual replicate `*.fixed_width.narrowPeak` files
- `*.renormalized.fixed_width.consensus.narrowPeak` per sample; [_Corces et. al._](https://doi.org/10.1038/nmeth.4396) method is used for consensus calling; used to generate MACS2 regions of interest (ROI) peaks which are used to generate a reads or Tn5 sites counts matrix for DESeq2
- ROI related files: `ROI.macs2.bed`, `ROI.macs2.bed.annotated.gz`, `ROI.macs2.annotation_summary`, `ROI.macs2.annotation_distribution` | -| `fixed_width/counts/{dedup,nondedup}/ROI.macs2.read_counts.tsv` | read counts in MACS2 ROIs using featureCounts; generated separately from `dedup.bam` (duplicates removed, recommended) and `filtered.bam`/`nondedup` (duplicates retained) | +| File | Description | +| -------------------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | +| `*.macs2.narrowPeak` | peak calls from MACS2 filtered by q-value for each samples each replicate | +| `*.macs2.unfiltered.narrowPeak` | peak calls from MACS2 (unfiltered) for each samples each replicate | +| `*.narrowPeak_motif_enrichment/ame_results.txt` | motif enrichment results from AME tool from MEME suite using HOCOMOCO v11 database | +| `*.narrowPeak_motif_enrichment/knownResults.txt` | motif enrichment results using HOMER with HOCOMOCO v11 database | +| `*.macs2.consensus.bed` | consensus peak call between multiple replicates of each sample. **Note:** consensus bed annotations are located in `QC/peak_annotations` | +| `DiffATAC/dedup/reads`, `DiffATAC/nondedup/reads` | folder containing differential open chromatin results:
- computated using read counts in MACS2 regions of interest (ROIs)
- `dedup` uses `dedup.bam`-derived counts (PCR/optical duplicates removed, recommended); `nondedup` uses `filtered.bam`-derived counts (duplicates retained)
- `all_diff_atacs.html` HTML report aggregated across all contrasts from `contrasts.tsv`
- `all_diff_atacs.tsv` DESeq2 results in TSV format aggregated across all contrasts from `contrasts.tsv`
- HTML and TSV file each per contrast in `contrasts.tsv` | +| `DiffATAC/dedup/tn5sites`, `DiffATAC/nondedup/tn5sites` | folder containing differential open chromatin results:
- computated using Tn5 nicking site counts in MACS2 regions of interest (ROIs)
- `dedup` uses `dedup.bam`-derived counts (PCR/optical duplicates removed, recommended); `nondedup` uses `filtered.bam`-derived counts (duplicates retained)
- `all_diff_atacs.html` HTML report aggregated across all contrasts from `contrasts.tsv`
- `all_diff_atacs.tsv` DESeq2 results in TSV format aggregated across all contrasts from `contrasts.tsv`
- HTML and TSV file each per contrast in `contrasts.tsv` | +| `fixed_width` | `fixed_width` can be set in `config.yaml` to create peaks of a user defined fixed width (default 500bp). This folder contains:
- individual replicate `*.fixed_width.narrowPeak` files
- `*.renormalized.fixed_width.consensus.narrowPeak` per sample; [_Corces et. al._](https://doi.org/10.1038/nmeth.4396) method is used for consensus calling; used to generate MACS2 regions of interest (ROI) peaks which are used to generate a reads or Tn5 sites counts matrix for DESeq2
- ROI related files: `ROI.macs2.bed`, `ROI.macs2.bed.annotated.gz`, `ROI.macs2.annotation_summary`, `ROI.macs2.annotation_distribution` | +| `fixed_width/counts/{dedup,nondedup}/ROI.macs2.read_counts.tsv` | read counts in MACS2 ROIs using featureCounts; generated separately from `dedup.bam` (duplicates removed, recommended) and `filtered.bam`/`nondedup` (duplicates retained) | | `fixed_width/counts/{dedup,nondedup}/ROI.reads_scaled_counts.tsv` | `ROI.macs2.read_counts.tsv` scaled using spike-in scaling factors | -| `fixed_width/counts/{dedup,nondedup}/ROI.tn5sites_counts.tsv` | Tn5 nicking site counts in MACS2 ROIs using featureCounts; generated separately from `dedup.bam` (duplicates removed, recommended) and `filtered.bam`/`nondedup` (duplicates retained) | +| `fixed_width/counts/{dedup,nondedup}/ROI.tn5sites_counts.tsv` | Tn5 nicking site counts in MACS2 ROIs using featureCounts; generated separately from `dedup.bam` (duplicates removed, recommended) and `filtered.bam`/`nondedup` (duplicates retained) | | `fixed_width/counts/{dedup,nondedup}/ROI.tn5sites_scaled_counts.tsv` | `ROI.macs2.tn5sites_counts.tsv` scaled using spike-in scaling factors | ### Genrich output folder diff --git a/docs/overview.md b/docs/overview.md index 545a630..ef5d646 100644 --- a/docs/overview.md +++ b/docs/overview.md @@ -24,10 +24,10 @@ Employs FastQC for initial quality checks and MultiQC for aggregated reporting, Two categories of reads are systematically removed during the preprocessing stage, since both would otherwise confound peak detection: -| Filtered out | Why | Source of coordinates | -|---|---|---| -| Blacklisted regions | Repetitive/artifact-prone genomic loci with unusually high, spurious signal | ENCODE (and other genomic consortia) blacklists, bundled per-genome in `resources/blacklistFa/` | -| Mitochondrial reads (`chrM`/`chrR`) | Often make up a large fraction of ATAC-seq reads (mitochondrial DNA is highly accessible); excluding them focuses compute and peak calling on nuclear chromatin | The reference genome's `chrM`/`chrR` contig | +| Filtered out | Why | Source of coordinates | +| ----------------------------------- | --------------------------------------------------------------------------------------------------------------------------------------------------------------- | ----------------------------------------------------------------------------------------------- | +| Blacklisted regions | Repetitive/artifact-prone genomic loci with unusually high, spurious signal | ENCODE (and other genomic consortia) blacklists, bundled per-genome in `resources/blacklistFa/` | +| Mitochondrial reads (`chrM`/`chrR`) | Often make up a large fraction of ATAC-seq reads (mitochondrial DNA is highly accessible); excluding them focuses compute and peak calling on nuclear chromatin | The reference genome's `chrM`/`chrR` contig | This filtering step enhances the accuracy and reliability of peak detection by concentrating on biologically meaningful and interpretable regions of the genome. @@ -35,7 +35,7 @@ This filtering step enhances the accuracy and reliability of peak detection by c Reads are aligned to the reference genome using Bowtie2, a fast and memory-efficient aligner optimized for short-read sequencing data. This step produces several key output files essential for subsequent analyses: -- **`filtered.bam`**: A BAM file with low-quality alignments and reads mapping to blacklisted regions removed. This is the **`nondedup`** BAM referenced elsewhere in these docs — i.e. the file *before* PCR/optical duplicate removal — and is used for the `nondedup` variant of the counts matrices and DiffATAC results. +- **`filtered.bam`**: A BAM file with low-quality alignments and reads mapping to blacklisted regions removed. This is the **`nondedup`** BAM referenced elsewhere in these docs — i.e. the file _before_ PCR/optical duplicate removal — and is used for the `nondedup` variant of the counts matrices and DiffATAC results. - **`qsorted.bam`**: A quality-sorted BAM file containing all aligned reads, serving as the primary input for Genrich peak calling. - **`dedup.bam`**: A BAM file derived from `filtered.bam` with PCR/optical duplicates additionally removed via PicardTools. Used for MACS2 peak calling, FRiP, transcription factor footprinting with TOBIAS, and the **recommended `dedup`** variant of the counts matrices and DiffATAC results. - **`tagAlign.gz`**: A compressed file containing alignment information in a simplified format, specifically prepared for MACS2 peak calling. @@ -59,33 +59,33 @@ Genrich is integrated into the pipeline to complement MACS2. This tool is partic By combining the strengths of MACS2 and Genrich, ASPEN delivers a comprehensive and reliable peak detection framework, facilitating downstream analyses and enabling researchers to uncover critical insights into chromatin accessibility and gene regulation. !!! tip "Which peak caller should I use?" - Both MACS2 and Genrich run automatically for every sample — you - don't have to choose upfront. **Genrich is generally - recommended for ATAC-seq** because it was purpose-built for - assays like ATAC-seq/DNase-seq: it models Tn5 transposase cut - sites directly (`-j` ATAC-seq mode) rather than adapting a - ChIP-seq fragment-shift model, and it combines biological - replicates natively via Fisher's method instead of requiring a - separate consensus step. **MACS2** was originally designed for - ChIP-seq and requires ATAC-specific parameter workarounds to - approximate cut-site signal, but remains the field standard — - it's included because many reviewers and downstream tools - expect to see MACS2 peaks, and comparing both gives an extra - sanity check. CCBR's internal benchmarking on ASPEN's ATAC-seq - data has generally found Genrich peaks to be higher quality. If - your MACS2 and Genrich DiffATAC results disagree substantially - for a given region, treat that region as lower-confidence - rather than assuming one caller is unconditionally "right". +Both MACS2 and Genrich run automatically for every sample — you +don't have to choose upfront. **Genrich is generally +recommended for ATAC-seq** because it was purpose-built for +assays like ATAC-seq/DNase-seq: it models Tn5 transposase cut +sites directly (`-j` ATAC-seq mode) rather than adapting a +ChIP-seq fragment-shift model, and it combines biological +replicates natively via Fisher's method instead of requiring a +separate consensus step. **MACS2** was originally designed for +ChIP-seq and requires ATAC-specific parameter workarounds to +approximate cut-site signal, but remains the field standard — +it's included because many reviewers and downstream tools +expect to see MACS2 peaks, and comparing both gives an extra +sanity check. CCBR's internal benchmarking on ASPEN's ATAC-seq +data has generally found Genrich peaks to be higher quality. If +your MACS2 and Genrich DiffATAC results disagree substantially +for a given region, treat that region as lower-confidence +rather than assuming one caller is unconditionally "right". ### 🤝 **Consensus Peaks** For datasets with multiple replicates, ASPEN generates three related-but-distinct peak representations — not three competing choices, but two sequential stages plus one diagnostic side-output: -| Output | Stage | Width | Best used for | -|---|---|---|---| -| `pooled.narrowPeak` | Diagnostic (pre-filter) — peaks called on all replicates' reads pooled together | Variable | Sanity-checking how strict the reproducibility filter (`consensus_min_replicates`/`consensus_min_spm`) is | -| `consensus.bed` | Round 1 — per-sample, reproducibility-filtered via the ["Consensus MAX" strategy](https://doi.org/10.5936/csbj.201401002) (_Yang et al._) | Variable | Peak annotation (ChIPseeker) and comparison against other pipelines' conventional peak outputs | -| `fixed_width.consensus.narrowPeak` | Round 2 — cross-sample ROI set, p-values renormalized across replicates following [_Corces et al._](https://doi.org/10.1038/nmeth.4396) | Fixed (`fixed_width` config, default 500 bp) | Differential accessibility testing (DESeq2) — fixed width avoids confounding read/Tn5 counts with peak-length differences and avoids daisy-chained merges across many samples | +| Output | Stage | Width | Best used for | +| ---------------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------- | -------------------------------------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| `pooled.narrowPeak` | Diagnostic (pre-filter) — peaks called on all replicates' reads pooled together | Variable | Sanity-checking how strict the reproducibility filter (`consensus_min_replicates`/`consensus_min_spm`) is | +| `consensus.bed` | Round 1 — per-sample, reproducibility-filtered via the ["Consensus MAX" strategy](https://doi.org/10.5936/csbj.201401002) (_Yang et al._) | Variable | Peak annotation (ChIPseeker) and comparison against other pipelines' conventional peak outputs | +| `fixed_width.consensus.narrowPeak` | Round 2 — cross-sample ROI set, p-values renormalized across replicates following [_Corces et al._](https://doi.org/10.1038/nmeth.4396) | Fixed (`fixed_width` config, default 500 bp) | Differential accessibility testing (DESeq2) — fixed width avoids confounding read/Tn5 counts with peak-length differences and avoids daisy-chained merges across many samples | **Use the fixed-width ROI set for differential accessibility testing; use `consensus.bed`/`pooled.narrowPeak` for peak annotation, QC, or cross-pipeline comparison.** @@ -116,21 +116,21 @@ ASPEN employs custom scripts to analyze the distribution of fragment lengths wit To evaluate the sufficiency of sequencing depth and detect potential biases introduced during PCR amplification, ASPEN utilizes Preseq to estimate library complexity, reporting the Non-Redundant Fraction (`NRF`) and PCR Bottlenecking Coefficients (`PBC1`, `PBC2`). This metric helps determine whether the sequencing effort is adequate to capture the diversity of the library, ensuring that the data is representative of the underlying chromatin landscape. By identifying potential saturation or over-representation of certain fragments, researchers can assess the reliability of their sequencing results. !!! tip "Rule of thumb" - Per [ENCODE's ATAC-seq data standards](https://www.encodeproject.org/atac-seq/#standards), the preferred values are **`NRF` > 0.9, `PBC1` > 0.9, and `PBC2` > 3**. Lower values indicate a less complex library (e.g. over-amplified by PCR), which can inflate apparent signal at a subset of loci rather than reflecting true biological accessibility. +Per [ENCODE's ATAC-seq data standards](https://www.encodeproject.org/atac-seq/#standards), the preferred values are **`NRF` > 0.9, `PBC1` > 0.9, and `PBC2` > 3**. Lower values indicate a less complex library (e.g. over-amplified by PCR), which can inflate apparent signal at a subset of loci rather than reflecting true biological accessibility. ### 🧬 **Transcription Start Site (TSS) Enrichment** ASPEN calculates TSS enrichment scores, a widely recognized quality metric for ATAC-seq data. These scores measure the accumulation of sequencing reads around transcription start sites (TSS), which are hallmark regions of open chromatin. High TSS enrichment scores indicate well-prepared libraries with minimal technical artifacts, as they reflect the accessibility of promoter regions and the integrity of the chromatin preparation process. !!! tip "Rule of thumb" - [ENCODE's ATAC-seq data standards](https://www.encodeproject.org/atac-seq/#standards) define annotation-dependent TSS enrichment cutoffs — for example, using a GRCh38 RefSeq TSS annotation: **< 5 is concerning, 5-7 is acceptable, and > 7 is ideal**. **Caveat:** ASPEN builds its TSS bins from GENCODE (not RefSeq) gene annotations (see `resources/tssBed/`), so ENCODE's exact per-annotation cutoffs may not transfer precisely to ASPEN's TSS enrichment values — treat these numbers as directional guidance (aim for high single digits or higher) rather than an exact pass/fail threshold. +[ENCODE's ATAC-seq data standards](https://www.encodeproject.org/atac-seq/#standards) define annotation-dependent TSS enrichment cutoffs — for example, using a GRCh38 RefSeq TSS annotation: **< 5 is concerning, 5-7 is acceptable, and > 7 is ideal**. **Caveat:** ASPEN builds its TSS bins from GENCODE (not RefSeq) gene annotations (see `resources/tssBed/`), so ENCODE's exact per-annotation cutoffs may not transfer precisely to ASPEN's TSS enrichment values — treat these numbers as directional guidance (aim for high single digits or higher) rather than an exact pass/fail threshold. ### 📊 **Fraction of Reads in Peaks (FRiP)** The Fraction of Reads in Peaks (FRiP) score quantifies the proportion of sequencing reads that fall within identified peaks, serving as a measure of the signal-to-noise ratio in the dataset. Higher FRiP scores indicate datasets with strong, biologically meaningful signals and minimal background noise. Additionally, ASPEN computes the fraction of reads localized to specific genomic features, such as promoters, enhancers, and DNase hypersensitive sites (DHS). These feature-specific FRiP scores provide further insights into the quality and biological relevance of the data. !!! tip "Rule of thumb" - Per [ENCODE's ATAC-seq data standards](https://www.encodeproject.org/atac-seq/#standards), a FRiP score **> 0.3** indicates high-quality data, though values **> 0.2** may still be acceptable. Consistently lower scores suggest poor signal-to-noise and warrant a closer look at library prep or peak-calling parameters. +Per [ENCODE's ATAC-seq data standards](https://www.encodeproject.org/atac-seq/#standards), a FRiP score **> 0.3** indicates high-quality data, though values **> 0.2** may still be acceptable. Consistently lower scores suggest poor signal-to-noise and warrant a closer look at library prep or peak-calling parameters. --- @@ -191,7 +191,7 @@ In ASPEN, if spike-in data is present: This spike-in-derived scaling factor allows the comparison of chromatin accessibility across conditions even when global chromatin accessibility levels differ (e.g., treatment-induced repression or global decondensation). !!! tip "Should I turn on spike-in normalization?" - Ask yourself: do I expect a **global, genome-wide shift** in chromatin accessibility between my conditions — rather than just **localized** changes at a handful of specific regulatory elements? +Ask yourself: do I expect a **global, genome-wide shift** in chromatin accessibility between my conditions — rather than just **localized** changes at a handful of specific regulatory elements? - **Yes** (e.g. a chromatin remodeler inhibitor, a broad transcription factor knockdown/knockout, drug-induced chromatin modulation) → turn on spike-in normalization. Standard depth-based normalization (DESeq2 size factors) assumes *most* regions are unchanged between conditions — that assumption breaks down under a genome-wide shift, and spike-in gives you an external, biology-independent scale instead. - **No** (you expect differences to be confined to specific loci/pathways, with most of the genome unchanged) → spike-in is probably unnecessary. It adds experimental complexity (extra reagents, a second alignment step, and a "0 spike-in reads" failure mode to manage — see the warning below) for a scenario DESeq2's built-in normalization already handles well.