# ATAC-seq vs DNase-seq: Concordance and Differences

Reference guide for understanding when ATAC-seq and DNase-seq are interchangeable and when platform-specific considerations matter.

## The Corces et al. 2017 Concordance Study

Corces et al. (2017, Nature Methods, 733 citations) performed the definitive comparison of ATAC-seq and DNase-seq for chromatin accessibility profiling. Key findings:

- ATAC-seq captures approximately **75% of DNase-seq hypersensitive sites**
- Sites detected by both assays show **highly concordant signal intensities** (Pearson r > 0.8)
- ATAC-specific sites tend to be **weaker** DNase sites (present but below DNase calling threshold)
- DNase-specific sites tend to be in **GC-rich regions** where Tn5 insertion is biased
- For regulatory element discovery, the two assays are **largely interchangeable**

**Implication for aggregation**: Combining ATAC-seq and DNase-seq peaks in a union set is scientifically justified. Both assays detect the same biological signal (open chromatin) through different enzymatic mechanisms.

## Enzymatic Mechanism Differences

| Property | ATAC-seq | DNase-seq |
|----------|----------|-----------|
| Enzyme | Tn5 transposase | DNase I endonuclease |
| Mechanism | Inserts sequencing adapters into open chromatin | Cleaves exposed DNA |
| Fragment sizes | Nucleosome-free (<150bp) + mono-nucleosome (~200bp) | Continuous size range |
| Resolution | High (single-nucleotide Tn5 insertion sites) | High (DNase I cut sites) |
| Input cells | 500 - 50,000 cells | 100,000 - 1,000,000 cells |
| Protocol complexity | Simple (1-2 hour protocol) | Complex (multi-day) |

## Known Biases

### Tn5 Sequence Preference (ATAC-seq)

Tn5 transposase has a mild sequence preference for insertion:
- Slight bias toward **10bp periodicity** matching nucleosome wrapping
- Moderate **GC content preference** at insertion sites
- These biases affect peak **intensity** (signal strength), not peak **presence**
- For union-based aggregation, this is acceptable since we care about detection, not quantification

### DNase I Cleavage Bias (DNase-seq)

DNase I also has sequence preferences:
- Slight preference for **WW dinucleotides** (W = A or T)
- Less GC-biased than Tn5
- Can create false hotspots at highly accessible repetitive elements

### Mitochondrial DNA Contamination (ATAC-seq)

ATAC-seq libraries typically contain **30-80% mitochondrial reads** because:
- Mitochondrial DNA is highly accessible (no histones)
- Tn5 readily inserts into mtDNA
- ENCODE pipeline removes these reads post-alignment
- Always verify mtDNA fraction in QC metrics

DNase-seq does **not** have this issue because the protocol uses nuclei (mitochondria excluded).

## Fragment Size Information

ATAC-seq uniquely provides chromatin structure information through fragment sizes:

| Fragment Class | Size Range | Represents |
|---------------|------------|------------|
| Sub-nucleosomal | < 100 bp | Open chromatin / Tn5 artifacts |
| Nucleosome-free | 100-150 bp | True open chromatin regions |
| Mono-nucleosome | 180-250 bp | Single nucleosome wrapping |
| Di-nucleosome | 315-475 bp | Two nucleosomes |

This information is **not available** in DNase-seq. For peak-based aggregation, fragment size filtering is typically already done by the ENCODE pipeline.

## When to Combine ATAC + DNase

### Recommended (same tissue, different assays)

Combine when building a comprehensive accessibility catalog:
- Maximizes genomic coverage
- Captures sites each assay alone might miss
- Union approach handles platform-specific sensitivity differences

### Caution Needed

- **Quantitative comparisons**: Do not directly compare signal intensities between ATAC and DNase peaks. Normalize separately.
- **Footprinting analysis**: Tn5 and DNase I have different cleavage profiles. Do not combine for transcription factor footprinting.
- **Nucleosome positioning**: Only ATAC-seq provides fragment size-based nucleosome information.

### Not Recommended

- **Differential accessibility**: Do not treat ATAC and DNase as replicates in differential analysis. Platform effects dominate.
- **Single-cell comparisons**: scATAC-seq and scDNase-seq (if available) have very different sparsity patterns.

## ENCODE Data Availability

| Assay | ENCODE Experiments | Period | Peak Caller |
|-------|-------------------|--------|-------------|
| DNase-seq | ~800+ experiments | 2007-present | Hotspot2, MACS2 |
| ATAC-seq | ~400+ experiments | 2015-present | MACS2 |

For tissues with both assay types available, combining yields the most comprehensive map. For tissues with only one assay, either is sufficient for accessibility cataloging.

## Practical Workflow for Combined Aggregation

```bash
# 1. Download ATAC and DNase peaks separately
# 2. Filter both through ENCODE blacklist
# 3. Apply signalValue filter to each sample independently
# 4. Tag peaks by sample AND assay type

awk -v sid="atac_donor1" 'BEGIN{OFS="\t"} {$4=sid; print}' atac_sample.narrowPeak > tagged.bed
awk -v sid="dnase_donor1" 'BEGIN{OFS="\t"} {$4=sid; print}' dnase_sample.narrowPeak >> tagged.bed

# 5. Union merge (treat as equal sources)
cat *.tagged.bed | bedtools sort | bedtools merge -c 4 -o count_distinct > union.bed

# 6. Optionally annotate whether support comes from ATAC, DNase, or both
```

## References

- Corces et al. 2017, Nature Methods -- definitive ATAC vs DNase concordance study (733 citations)
- Buenrostro et al. 2013, Nature Methods -- original ATAC-seq protocol (6,800+ citations)
- Thurman et al. 2012, Nature -- DNase-seq regulatory landscape (ENCODE Phase 2)
- Yan et al. 2020, Genome Biology -- comprehensive Tn5 sequence bias characterization
- Sung et al. 2014, Nature Methods -- DNase-seq protocol optimization
- Zhao & Boyle 2020, NAR Genomics -- F-Seq2 peak caller for both ATAC and DNase

# Accessibility Aggregation — Literature References

**Last updated:** 2026-03-07
**Purpose:** Reference catalog for the accessibility-aggregation skill — papers supporting the union-based approach for aggregating ATAC-seq and DNase-seq chromatin accessibility peaks across experiments, and the comparability of these two accessibility assay platforms.

---

## ATAC-seq / DNase-seq Comparability

---

### Corces et al. 2017 — Omni-ATAC and DNase-seq comparison

- **Citation:** Corces MR, Trevino AE, Hamilton EG, Greenside PG, Sinnott-Armstrong NA, Vesuna S, Satpathy AT, Rubin AJ, Montine KS, Wu B, Kathiria A, Cho SW, Mumbach MR, Carter AC, Kasowski M, Orloff LA, Risca VI, Kundaje A, Khavari PA, Montine TJ, Greenleaf WJ, Chang HY. An improved ATAC-seq protocol reduces background and enables interrogation of frozen tissues. Nature Methods, 14(10):959-962, 2017.
- **DOI:** [10.1038/nmeth.4396](https://doi.org/10.1038/nmeth.4396)
- **PMID:** 28846090 | **PMC:** PMC5623106
- **Citations:** ~1,925
- **Key findings:** Established that ATAC-seq and DNase-seq identify largely overlapping sets of accessible chromatin regions, with ATAC-seq capturing ~75% of DNase-seq hypersensitive sites. The concordance supports combining ATAC-seq and DNase-seq peaks in union aggregation maps for comprehensive accessibility catalogs. Platform-specific differences (Tn5 insertion bias in ATAC-seq, requirement for higher cell input for DNase-seq) affect sensitivity at individual sites but not the validity of the union approach.

---

### Thurman et al. 2012 — The accessible chromatin landscape

- **Citation:** Thurman RE, Rynes E, Humbert R, et al. The accessible chromatin landscape of the human genome. Nature, 489(7414):75-82, 2012.
- **DOI:** [10.1038/nature11232](https://doi.org/10.1038/nature11232)
- **PMID:** 22955617 | **PMC:** PMC3721348
- **Citations:** ~3,000
- **Key findings:** Definitive DNase-seq accessibility atlas across 125 human cell and tissue types. Identified ~2.9 million unique DNase I hypersensitive sites covering ~40% of the genome, with individual cell types having 100,000-200,000 DHSs. Demonstrated that distal accessible sites (enhancers) are highly cell-type-specific while promoter-proximal sites are shared. This atlas provides the conceptual framework for the accessibility aggregation approach — combining sites across samples to build comprehensive tissue-level maps.

---

### Buenrostro et al. 2013 — ATAC-seq method

- **Citation:** Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nature Methods, 10(12):1213-1218, 2013.
- **DOI:** [10.1038/nmeth.2688](https://doi.org/10.1038/nmeth.2688)
- **PMID:** 24097267 | **PMC:** PMC3959825
- **Citations:** ~5,600
- **Key findings:** Original ATAC-seq method paper establishing Tn5 transposase-based chromatin accessibility profiling. ATAC-seq requires far fewer cells than DNase-seq (500-50,000 vs >100,000) and uses a simpler protocol, making it the dominant accessibility assay in modern experiments. The nucleosomal ladder pattern in ATAC-seq fragment sizes (sub-nucleosomal <150 bp, mono-nucleosome 150-300 bp) provides additional biological information not available from DNase-seq. Both assays fundamentally measure the same biology — open chromatin regions where regulatory factors bind.

---

## Peak Calling and Quality

---

### Zhao & Boyle 2021 — F-Seq2: improved peak calling for accessibility data

- **Citation:** Zhao Z, Boyle AP. F-Seq2: improving the feature density based peak caller with dynamic statistics. NAR Genomics and Bioinformatics, 3(1):lqab012, 2021.
- **DOI:** [10.1093/nargab/lqab012](https://doi.org/10.1093/nargab/lqab012)
- **PMID:** 33655203 | **PMC:** PMC7899645
- **Citations:** ~30
- **Key findings:** Introduced F-Seq2, an improved feature density peak caller for DNase-seq and ATAC-seq that provides proper test statistics (p-values) enabling IDR analysis across replicates. Unlike the original F-Seq which used kernel density estimation without statistical testing, F-Seq2 models background with a dynamic Poisson distribution and supports both narrow and broad peak modes. Provides an alternative to MACS2 (for ATAC-seq) and Hotspot2 (for DNase-seq) that works uniformly across both accessibility platforms, which is advantageous for mixed-platform aggregation.

---

### Amemiya et al. 2019 — ENCODE Blacklist

- **Citation:** Amemiya HM, Kundaje A, Boyle AP. The ENCODE Blacklist: Identification of Problematic Regions of the Genome. Scientific Reports, 9:9354, 2019.
- **DOI:** [10.1038/s41598-019-45839-z](https://doi.org/10.1038/s41598-019-45839-z)
- **PMID:** 31249361
- **Citations:** ~1,372
- **Key findings:** Blacklist v2 filtering is essential before accessibility peak aggregation. Both Tn5 (ATAC-seq) and DNase I preferentially access certain repetitive regions in the blacklist, creating artifact peaks. Without filtering, these artifacts would appear as high-confidence accessible regions in aggregated maps because they are consistently detected across experiments — but they represent technical artifacts, not genuine regulatory elements.

---

### Orchard et al. 2020 — ataqv: ATAC-seq quality control

- **Citation:** Orchard P, Kyono Y, Hensley J, Kitzman JO, Parker SCJ. Quantification, Dynamic Visualization, and Validation of Bias in ATAC-Seq Data with ataqv. Cell Systems, 10(3):298-306.e4, 2020.
- **DOI:** [10.1016/j.cels.2020.02.009](https://doi.org/10.1016/j.cels.2020.02.009)
- **PMID:** 32213349 | **PMC:** PMC7138743
- **Citations:** ~62
- **Key findings:** Identified TSS enrichment as the primary quality indicator for ATAC-seq experiments. Analysis of 2,009 public ATAC-seq datasets revealed a 10-fold range in quality metrics, establishing that quality-gating is essential before aggregation. ATAC-seq experiments with TSS enrichment <6 should be excluded from aggregation to prevent dilution of the accessibility signal with background noise.

---

## Data Integration Framework

---

### ENCODE Project Consortium 2020 — Expanded encyclopaedias of DNA elements

- **Citation:** ENCODE Project Consortium et al. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature, 583(7818):699-710, 2020.
- **DOI:** [10.1038/s41586-020-2493-4](https://doi.org/10.1038/s41586-020-2493-4)
- **PMID:** 32728249 | **PMC:** PMC7410828
- **Citations:** ~2,500
- **Key findings:** ENCODE Phase 3 paper establishing the cCRE (candidate cis-Regulatory Element) classification system that uses integrated accessibility data. Both ATAC-seq and DNase-seq DHSs serve as core inputs for identifying cCREs, supporting the equivalence of these platforms for accessibility cataloging. The cCRE classification combines accessibility peaks with histone mark data to distinguish promoter-like, enhancer-like, and CTCF-only elements.

---

### Gorkin et al. 2020 — Atlas of chromatin landscapes in mouse development

- **Citation:** Gorkin DU, Barozzi I, Zhao Y, et al. An atlas of dynamic chromatin landscapes in mouse fetal development. Nature, 583(7818):744-751, 2020.
- **DOI:** [10.1038/s41586-020-2093-3](https://doi.org/10.1038/s41586-020-2093-3)
- **PMID:** 32728240 | **PMC:** PMC7402670
- **Citations:** ~301
- **Key findings:** Integrated accessibility data across tissues and developmental stages using union peak sets combined with histone marks for chromatin state annotation. Demonstrated that comprehensive accessibility catalogs from union merging enable discovery of tissue-specific enhancers and developmental regulatory switches. Supports the aggregation-then-annotate workflow used in this skill.

---

## Computational Tools

---

### Quinlan & Hall 2010 — BEDTools

- **Citation:** Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics, 26(6):841-842, 2010.
- **DOI:** [10.1093/bioinformatics/btq033](https://doi.org/10.1093/bioinformatics/btq033)
- **PMID:** 20110278
- **Citations:** ~12,000
- **Key findings:** Core computational tool for accessibility peak aggregation. `bedtools merge` performs the union merge; `bedtools intersect -v` removes blacklisted regions; `bedtools multiIntersect` tracks which samples support each accessible region. The same BEDTools workflow is used for both ATAC-seq and DNase-seq peaks, enabling seamless cross-platform aggregation.

#!/usr/bin/env python3
"""Validate narrowPeak files from ENCODE ATAC-seq/DNase-seq accessibility aggregation.

Checks narrowPeak format, coordinate validity, typical accessibility peak sizes,
Tn5 bias artifacts, and reports summary statistics.

Usage:
    python validate_peaks.py input.narrowPeak [--blacklist hg38-blacklist.v2.bed] [--assay atac|dnase]
    python validate_peaks.py input.narrowPeak --assay atac --blacklist hg38-blacklist.v2.bed

Plain and gzipped (.gz) inputs and blacklists are both accepted.
"""

import argparse
import gzip
import statistics
import sys
from collections import Counter, defaultdict
from pathlib import Path

VALID_CHROMS = {f"chr{i}" for i in range(1, 23)} | {"chrX", "chrY", "chrM"}
NARROWPEAK_COLS = 10

# Typical ATAC-seq peak widths: 100-500bp (nucleosome-free regions)
# Peaks >2kb are suspicious for accessibility data
ATAC_LARGE_THRESHOLD = 2000
DNASE_LARGE_THRESHOLD = 3000
ATAC_SMALL_THRESHOLD = 50  # Sub-nucleosomal Tn5 artifacts

MAX_COLUMN_ERRORS = 5
MAX_WARNINGS = 20


def parse_args():
    parser = argparse.ArgumentParser(
        description="Validate narrowPeak files from ATAC-seq/DNase-seq accessibility aggregation.",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog=(
            "Examples:\n"
            "  python validate_peaks.py sample.narrowPeak\n"
            "  python validate_peaks.py sample.narrowPeak --assay atac\n"
            "  python validate_peaks.py sample.narrowPeak --blacklist hg38-blacklist.v2.bed\n"
        ),
    )
    parser.add_argument("input", type=Path, help="Input narrowPeak file")
    parser.add_argument(
        "--blacklist",
        type=Path,
        default=None,
        help="ENCODE blacklist BED file (e.g., hg38-blacklist.v2.bed)",
    )
    parser.add_argument(
        "--assay",
        choices=["atac", "dnase", "unknown"],
        default="unknown",
        help="Assay type for assay-specific checks. Default: unknown",
    )
    return parser.parse_args()


def open_text(path):
    """Open a plain or gzipped text file for reading."""
    if str(path).endswith(".gz"):
        return gzip.open(path, "rt")
    return open(path)


def quartiles(values):
    """25th percentile, median and 75th percentile (interpolated)."""
    median = statistics.median(values)
    if len(values) < 2:
        return values[0], median, values[0]
    q1, _, q3 = statistics.quantiles(values, n=4, method="inclusive")
    return q1, median, q3


def load_blacklist(path):
    """Load blacklist regions as a dict of chrom -> list of (start, end)."""
    regions = defaultdict(list)
    with open_text(path) as f:
        for line in f:
            if line.startswith("#") or line.strip() == "":
                continue
            parts = line.strip().split("\t")
            if len(parts) >= 3:
                regions[parts[0]].append((int(parts[1]), int(parts[2])))
    for chrom in regions:
        regions[chrom].sort()
    return regions


def overlaps_blacklist(chrom, start, end, blacklist):
    """Check if a region overlaps any blacklist interval."""
    if chrom not in blacklist:
        return False
    for bl_start, bl_end in blacklist[chrom]:
        if bl_start >= end:
            break
        if bl_end > start:
            return True
    return False


def validate_accessibility_peaks(input_path, blacklist_path, assay):
    errors = []
    warnings = []
    dropped_warnings = 0
    chrom_counts = Counter()
    peak_sizes = []
    signal_values = []
    start_positions = Counter()  # Track exact start positions for Tn5 pileup detection
    total_lines = 0
    # Every line excluded from the statistics below counts as malformed, so
    # total_lines == valid_peaks + bad_lines always holds.
    bad_lines = 0
    column_errors = 0
    blacklist_overlaps = 0
    chrm_peaks = 0
    large_peaks = 0
    tiny_peaks = 0

    large_threshold = ATAC_LARGE_THRESHOLD if assay == "atac" else DNASE_LARGE_THRESHOLD

    blacklist = None
    if blacklist_path:
        if not blacklist_path.exists():
            print(f"ERROR: Blacklist file not found: {blacklist_path}", file=sys.stderr)
            sys.exit(1)
        blacklist = load_blacklist(blacklist_path)

    if not input_path.exists():
        print(f"ERROR: Input file not found: {input_path}", file=sys.stderr)
        sys.exit(1)

    with open_text(input_path) as f:
        for line_num, line in enumerate(f, 1):
            if line.startswith("#") or line.startswith("track") or line.startswith("browser"):
                continue
            line = line.strip()
            if not line:
                continue

            total_lines += 1
            fields = line.split("\t")

            # Accessibility peaks should always be narrowPeak
            if len(fields) != NARROWPEAK_COLS:
                bad_lines += 1
                column_errors += 1
                if column_errors <= MAX_COLUMN_ERRORS:
                    errors.append(
                        f"Line {line_num}: expected {NARROWPEAK_COLS} columns (narrowPeak), got {len(fields)}"
                    )
                elif column_errors == MAX_COLUMN_ERRORS + 1:
                    errors.append("... suppressing further column-count errors")
                continue

            chrom = fields[0]
            if chrom not in VALID_CHROMS:
                if not chrom.startswith("chr"):
                    errors.append(f"Line {line_num}: invalid chromosome '{chrom}'")
                    bad_lines += 1
                    continue
                if len(warnings) < MAX_WARNINGS:
                    warnings.append(
                        f"Line {line_num}: non-standard chromosome '{chrom}' (not in chr1-22, chrX, chrY, chrM)"
                    )
                else:
                    dropped_warnings += 1

            try:
                start = int(fields[1])
                end = int(fields[2])
            except ValueError:
                errors.append(f"Line {line_num}: non-integer coordinates")
                bad_lines += 1
                continue

            if start < 0:
                errors.append(f"Line {line_num}: negative start coordinate ({start})")
            if end < 0:
                errors.append(f"Line {line_num}: negative end coordinate ({end})")
            if start >= end:
                errors.append(f"Line {line_num}: start ({start}) >= end ({end})")
            # an impossible interval is malformed: count it once and keep it out of the statistics
            if start < 0 or end < 0 or start >= end:
                bad_lines += 1
                continue

            # signalValue, pValue and qValue are checked before anything is counted: a row with
            # an unusable value is malformed and stays out of every statistic
            scores = {}
            for col_idx, col_name in [(6, "signalValue"), (7, "pValue"), (8, "qValue")]:
                try:
                    scores[col_name] = float(fields[col_idx])
                except ValueError:
                    errors.append(f"Line {line_num}: invalid {col_name} in column {col_idx + 1}")
            if scores.get("signalValue", 0) < 0:
                errors.append(f"Line {line_num}: negative signalValue ({scores['signalValue']})")
            if len(scores) < 3 or scores["signalValue"] < 0:
                bad_lines += 1
                continue

            peak_size = end - start
            peak_sizes.append(peak_size)
            chrom_counts[chrom] += 1
            signal_values.append(scores["signalValue"])

            # Track start positions for Tn5 pileup detection (ATAC-specific)
            if assay == "atac":
                start_positions[(chrom, start)] += 1

            if chrom == "chrM":
                chrm_peaks += 1

            if peak_size > large_threshold:
                large_peaks += 1

            if peak_size < ATAC_SMALL_THRESHOLD:
                tiny_peaks += 1

            # Blacklist overlap check
            if blacklist and overlaps_blacklist(chrom, start, end, blacklist):
                blacklist_overlaps += 1

    if total_lines == 0:
        print(f"ERROR: no data rows in {input_path} (only comments, headers or blank lines)", file=sys.stderr)
        sys.exit(1)

    valid_peaks = total_lines - bad_lines

    # --- Tn5 Pileup Detection (ATAC-specific) ---
    tn5_pileup_count = 0
    if assay == "atac" and start_positions:
        tn5_pileup_count = sum(1 for count in start_positions.values() if count >= 5)

    # --- Report Statistics ---
    assay_label = assay.upper() if assay != "unknown" else "Accessibility"
    print(f"=== {assay_label} NarrowPeak Validation Report ===")
    print(f"File: {input_path}")
    print(f"Assay: {assay}")
    print()

    print("--- Summary ---")
    print(f"Data lines: {total_lines:,}")
    print(f"Valid peaks: {valid_peaks:,}")
    print(f"Malformed lines: {bad_lines}")
    if blacklist_path:
        print(f"Blacklist overlaps: {blacklist_overlaps:,} ({100 * blacklist_overlaps / max(valid_peaks, 1):.1f}%)")
    print()

    if peak_sizes:
        sorted_sizes = sorted(peak_sizes)
        n = len(sorted_sizes)
        size_q1, size_median, size_q3 = quartiles(sorted_sizes)

        # Count peaks in expected accessibility range (100-500bp)
        in_range = sum(1 for s in peak_sizes if 100 <= s <= 500)
        in_range_pct = 100 * in_range / max(n, 1)

        print("--- Peak Size Distribution ---")
        print(f"Min:    {sorted_sizes[0]:,} bp")
        print(f"25th:   {size_q1:,.1f} bp")
        print(f"Median: {size_median:,.1f} bp")
        print(f"75th:   {size_q3:,.1f} bp")
        print(f"Max:    {sorted_sizes[-1]:,} bp")
        print(f"In typical range (100-500bp): {in_range:,} ({in_range_pct:.1f}%)")
        print()

    if signal_values:
        sorted_sig = sorted(signal_values)
        sig_q1, sig_median, sig_q3 = quartiles(sorted_sig)
        print("--- SignalValue Distribution ---")
        print(f"Min:    {sorted_sig[0]:.2f}")
        print(f"25th:   {sig_q1:.2f}")
        print(f"Median: {sig_median:.2f}")
        print(f"75th:   {sig_q3:.2f}")
        print(f"Max:    {sorted_sig[-1]:.2f}")
        print()

    print("--- Chromosome Distribution ---")
    for chrom in sorted(chrom_counts.keys(), key=lambda c: (len(c), c)):
        count = chrom_counts[chrom]
        pct = 100 * count / max(valid_peaks, 1)
        print(f"  {chrom:<6} {count:>8,}  ({pct:5.1f}%)")
    print()

    # --- Warnings ---
    if chrm_peaks > 0:
        msg = (
            f"WARNING: {chrm_peaks:,} peaks on chrM. "
            f"Mitochondrial peaks are common ATAC-seq artifacts (high mito read fraction)."
        )
        print(msg, file=sys.stderr)

    if large_peaks > 0:
        msg = (
            f"WARNING: {large_peaks:,} peaks exceed {large_threshold}bp. "
            f"Accessibility peaks are typically 100-500bp. "
            f"Large peaks may indicate artifacts or broad-mark contamination."
        )
        print(msg, file=sys.stderr)

    if tiny_peaks > 0 and assay == "atac":
        msg = (
            f"WARNING: {tiny_peaks:,} peaks are <{ATAC_SMALL_THRESHOLD}bp. "
            f"Very narrow peaks in ATAC-seq can be Tn5 insertion artifacts."
        )
        print(msg, file=sys.stderr)

    if tn5_pileup_count > 0:
        msg = (
            f"WARNING: {tn5_pileup_count:,} genomic positions have 5+ peaks "
            f"sharing the exact same start coordinate. This may indicate Tn5 "
            f"insertion bias (positional pileup artifact)."
        )
        print(msg, file=sys.stderr)

    if blacklist_overlaps > 0:
        msg = (
            f"WARNING: {blacklist_overlaps:,} peaks overlap ENCODE blacklist regions. Remove these before aggregation."
        )
        print(msg, file=sys.stderr)

    for w in warnings:
        print(w, file=sys.stderr)
    if dropped_warnings:
        print(f"  ... and {dropped_warnings} more warning(s) suppressed", file=sys.stderr)

    # --- Errors ---
    if errors:
        print(f"\n--- Errors ({len(errors)}) ---", file=sys.stderr)
        for e in errors[:50]:
            print(f"  {e}", file=sys.stderr)
        if len(errors) > 50:
            print(f"  ... and {len(errors) - 50} more errors", file=sys.stderr)

    has_errors = len(errors) > 0
    if has_errors:
        print(f"\nRESULT: FAIL -- {len(errors)} error(s) found", file=sys.stderr)
    else:
        print("\nRESULT: PASS -- file is valid narrowPeak")

    return 1 if has_errors else 0


if __name__ == "__main__":
    args = parse_args()
    exit_code = validate_accessibility_peaks(args.input, args.blacklist, args.assay)
    sys.exit(exit_code)