pipeline-dnaseseq

ENCODE DNase-seq Pipeline: FASTQ to Hotspots and Footprints

When to Use

• User wants to run a DNase-seq processing pipeline from FASTQ to hotspots and footprints
• User asks about "DNase-seq pipeline", "DNase hypersensitive sites", "Hotspot2", "footprinting", or "DHS"
• User needs to process DNase-seq data for chromatin accessibility and TF footprint analysis
• Example queries: "process my DNase-seq FASTQs", "call DNase hypersensitive sites", "run footprinting analysis on DNase-seq"

Execute the ENCODE DNase-seq pipeline for chromatin accessibility profiling, producing DNase hypersensitive sites (DHSs) via Hotspot2 and transcription factor footprints.

Pipeline Overview

FASTQ -> Trim -> BWA-MEM align -> Filter/dedup -> Hotspot2 -> DHS peaks
                                       |                        |
                                    Signal track         Footprinting (HINT)

ENCODE Repository

• GitHub: ENCODE-DCC/dnase-seq-pipeline
• Container: encodedcc/dnase-seq-pipeline
• WDL: Available for Cromwell execution
• This skill: Nextflow DSL2 reimplementation for portability

Core Tools and Versions

| Tool | Version | Purpose | Citation | |------|---------|---------|----------| | BWA-MEM | 0.7.17 | Alignment | Li & Durbin 2009 | | samtools | 1.19 | BAM operations | Li et al. 2009 | | Picard | 3.1.1 | Duplicate marking | Broad Institute | | Hotspot2 | 2.3.1 | DHS calling (ENCODE standard) | John et al. 2011 | | bedtools | 2.31.0 | Genomic arithmetic | Quinlan & Hall 2010 | | HINT-ATAC | 0.13.2 | TF footprinting | Li et al. 2019 | | FastQC | 0.12.1 | Read quality | Andrews (Babraham) | | MultiQC | 1.21 | Aggregated QC | Ewels et al. 2016 |

Key Literature

1. John et al. 2011 - "Chromatin accessibility pre-determines glucocorticoid

receptor binding patterns" (Nature Genetics, ~600 citations) DOI: 10.1038/ng.759

1. Thurman et al. 2012 - "The accessible chromatin landscape of the human

genome" (Nature, ~3,000 citations) DOI: 10.1038/nature11232

1. Vierstra et al. 2020 - "Global reference mapping of human transcription

factor footprints" (Nature, ~600 citations) DOI: 10.1038/s41586-020-2528-x

1. Amemiya et al. 2019 - "The ENCODE Blacklist" (Scientific Reports, ~1,372 citations)

DOI: 10.1038/s41598-019-45839-z

1. Li et al. 2019 - "Identification of transcription factor binding sites using

ATAC-seq" (Genome Biology) -- HINT-ATAC footprinting DOI: 10.1186/s13059-019-1642-2

Execution

Quick Start (Local)

nextflow run main.nf \
    -profile local \
    --reads '/data/fastq/*_R{1,2}.fastq.gz' \
    --bwa_index '/ref/bwa_index/genome.fa' \
    --chrom_sizes '/ref/hg38.chrom.sizes' \
    --hotspot_index '/ref/hotspot2_index/' \
    --blacklist '/ref/hg38-blacklist.v2.bed' \
    --outdir results/ \
    -resume

SLURM HPC

nextflow run main.nf \
    -profile slurm \
    --reads '/data/fastq/*_R{1,2}.fastq.gz' \
    --bwa_index '/ref/bwa_index/genome.fa' \
    --chrom_sizes '/ref/hg38.chrom.sizes' \
    --hotspot_index '/ref/hotspot2_index/' \
    --blacklist '/ref/hg38-blacklist.v2.bed' \
    --outdir results/ \
    -resume

Cloud (GCP / AWS)

nextflow run main.nf \
    -profile gcp \
    --reads 'gs://bucket/fastq/*_R{1,2}.fastq.gz' \
    --bwa_index 'gs://bucket/ref/genome.fa' \
    --chrom_sizes 'gs://bucket/ref/hg38.chrom.sizes' \
    --hotspot_index 'gs://bucket/ref/hotspot2_index/' \
    --blacklist 'gs://bucket/ref/hg38-blacklist.v2.bed' \
    --outdir 'gs://bucket/results/' \
    -resume

Resource Requirements

| Step | CPUs | RAM | Time (per sample) | |------|------|-----|-------------------| | BWA-MEM align | 8 | 16 GB | 1-2 hours | | Filter/dedup | 4 | 8 GB | 30-60 min | | Hotspot2 | 4 | 8 GB | 30-60 min | | Signal generation | 2 | 4 GB | 15-30 min | | Footprinting | 4 | 8 GB | 1-2 hours | | Total | 8 | 16 GB | 3-6 hours |

Pipeline Parameters

| Parameter | Default | Description | |-----------|---------|-------------| | --reads | required | Glob pattern to paired FASTQ files | | --bwa_index | required | Path to BWA genome index (.fa) | | --chrom_sizes | required | Chromosome sizes file | | --hotspot_index | required | Hotspot2 mappability index directory | | --blacklist | required | ENCODE blacklist BED file | | --outdir | ./results | Output directory | | --fdr | 0.05 | Hotspot2 FDR threshold | | --skip_footprint | false | Skip footprinting analysis | | --motif_db | null | JASPAR motif database for footprinting |

Output Files

results/
  fastqc/                          # Raw read quality
  alignment/
    {sample}.filtered.bam          # Filtered, deduplicated BAM
    {sample}.filtered.bam.bai
  hotspots/
    {sample}.hotspots.fdr0.05.bed  # DHS peaks (primary output)
    {sample}.peaks.narrowPeak      # narrowPeak format
    {sample}.density.bw            # Signal track (bigWig)
    {sample}.allcalls.bed          # All hotspot calls (unfiltered)
    {sample}.SPOT.txt              # SPOT score
  footprints/
    {sample}.footprints.bed        # TF footprints
    {sample}.footprint_scores.txt  # Per-motif footprint scores
  qc/
    {sample}.flagstat.txt
    {sample}.insert_sizes.txt
  multiqc/
    multiqc_report.html

QC Thresholds (ENCODE Standards)

| Metric | Pass | Warning | Fail | |--------|------|---------|------| | SPOT score (Signal Portion of Tags) | >0.4 | 0.2-0.4 | <0.2 | | Hotspot count | >50,000 | 20,000-50,000 | <20,000 | | Mapping rate | >80% | 60-80% | <60% | | Duplication rate | <30% | 30-50% | >50% | | NRF (Non-Redundant Fraction) | >0.8 | 0.7-0.8 | <0.7 | | PBC1 (PCR Bottleneck Coefficient 1) | >0.9 | 0.7-0.9 | <0.7 | | Insert size peak | 50-150 bp | Variable | Abnormal |

SPOT Score

The SPOT score (Signal Portion of Tags) is the fraction of reads falling within hotspots. It is the DNase-seq equivalent of FRiP for ChIP-seq.

Higher SPOT = more enrichment in accessible regions = better library quality.

Hotspot2 vs MACS2

IMPORTANT: ENCODE uses Hotspot2 for DNase-seq, NOT MACS2.

| Feature | Hotspot2 | MACS2 | |---------|----------|-------| | Designed for | DNase-seq | ChIP-seq | | Background model | Local tag density + mappability | Dynamic Poisson | | ENCODE standard | Yes (DNase-seq) | Yes (ChIP-seq/ATAC-seq) | | Mappability correction | Built-in | Not available | | Output | Hotspots + peaks | Peaks only |

Hotspot2 accounts for mappability variation across the genome, which is critical for DNase-seq because DNase I cuts accessible chromatin regardless of whether it is uniquely mappable.

Critical Pitfalls

DNase-seq vs ATAC-seq

These are different assays measuring the same biology (chromatin accessibility):

• DNase-seq: Uses DNase I enzyme, requires more input material
• ATAC-seq: Uses Tn5 transposase, works on fewer cells
• Analysis pipelines differ: Hotspot2 for DNase-seq, MACS2 for ATAC-seq
• Data are largely concordant but not identical

Fragment Size Distribution

DNase-seq produces a characteristic fragment size distribution:

• Peak at ~50-100 bp (sub-nucleosomal fragments at DHS)
• Secondary peak at ~150-200 bp (mononucleosomal fragments)
• Long tail of larger fragments
• If distribution is abnormal, check library preparation protocol

Mappability Index

Hotspot2 requires a pre-computed mappability index. These are read-length and genome-build specific:

• hg38 / 36 bp: Use ENCODE-provided index
• hg38 / 76 bp: Use ENCODE-provided index
• hg38 / 150 bp: May need to generate custom index
• Wrong mappability index = incorrect peak calls

Blacklist Filtering

Always filter peaks against the ENCODE blacklist (Amemiya et al. 2019):

bedtools intersect -a hotspots.bed -b hg38-blacklist.v2.bed -v > hotspots_filtered.bed

Blacklist regions produce artifactual signal in accessibility assays.

Footprinting Analysis

Transcription factor footprinting detects bound TFs f