git:20260920.92639e7 to git:20260920.ce7840d

---

name: pipeline-dnaseseq

description: "Execute ENCODE DNase-seq pipeline from FASTQ to hotspots and footprints. Child of pipeline-guide. Provides Nextflow execution with Docker and cloud deployment. Use when processing DNase-seq data, calling DNase hypersensitive sites, performing footprinting analysis. Trigger on: DNase-seq pipeline, DNase hypersensitive, DHS, Hotspot2, footprinting, DNase I, chromatin accessibility DNase."

---

# 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**: built from `scripts/Dockerfile` in this skill (`docker build -t encode-toolkit/pipeline-dnaseseq:1.0.0 scripts/`); override with `--container`

- **WDL**: Available for Cromwell execution

- **This skill**: Nextflow DSL2 reimplementation for portability

## Core Tools and Versions

| Tool | Version | Purpose | Citation |

|------|---------|---------|----------|

| BWA-MEM | 0.7.18 | Alignment | Li & Durbin 2009 |

| samtools | 1.19 | BAM operations | Li et al. 2009 |

| Picard | 3.1.1 | Duplicate marking | Broad Institute |

| Hotspot2 | 2.1.2 | DHS calling (ENCODE standard) | John et al. 2011 |

| modwt | 1.0 | Wavelet smoothing used by Hotspot2 | Stam Lab |

| bedtools | 2.31.0 | Genomic arithmetic | Quinlan & Hall 2010 |

| HINT (RGT) | 1.0.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

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

genome" (Nature, ~3,000 citations)

DOI: 10.1038/nature11232

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

factor footprints" (Nature, ~600 citations)

DOI: 10.1038/s41586-020-2528-x

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

DOI: 10.1038/s41598-019-45839-z

5. **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)

```bash

-profile local \

--reads '/data/fastq/*_R{1,2}.fastq.gz' \

--bwa_index '/ref/bwa_index/genome.fa' \

--chrom_sizes '/ref/hg38.chrom.sizes' \

--hotspot_center_sites '/ref/hotspot2/hg38.center_sites.n100.starch' \

--hotspot_mappable '/ref/hotspot2/hg38.mappable_only.bed' \

--rgt_data '/ref/rgtdata' \

--blacklist '/ref/hg38-blacklist.v2.bed' \

--outdir results/ \

-resume

```

### SLURM HPC

```bash

-profile slurm \

--reads '/data/fastq/*_R{1,2}.fastq.gz' \

--bwa_index '/ref/bwa_index/genome.fa' \

--chrom_sizes '/ref/hg38.chrom.sizes' \

--hotspot_center_sites '/ref/hotspot2/hg38.center_sites.n100.starch' \

--hotspot_mappable '/ref/hotspot2/hg38.mappable_only.bed' \

--rgt_data '/ref/rgtdata' \

--blacklist '/ref/hg38-blacklist.v2.bed' \

--outdir results/ \

-resume

```

### Cloud (GCP / AWS)

```bash

```

## 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 |

|-----------|---------|-------------|

| `--outdir` | `./results` | Output directory |

| `--skip_footprint` | `false` | Skip footprinting analysis |

| `--organism` | `hg38` | Genome name registered in the RGT data directory, used by HINT footprinting |

| `--rgt_data` | required unless `--skip_footprint` | RGT data directory with the genome for `--organism` set up (see below) |

## Output Files

```

results/

alignment/

{sample}.filtered.bam.bai

hotspots/

{sample}.SPOT.txt # SPOT score

footprints/

qc/

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 needs a center-sites file, which is derived from a mappable-regions BED. Both are

read-length and genome-build specific. Create the center sites once per genome with the script

that ships with Hotspot2 (it is on the PATH inside the image):

```bash

# chrom_sizes.bed is a BED file: chromosome, 0, length

awk 'BEGIN{OFS="\t"} {print $1, 0, $2}' hg38.chrom.sizes | sort-bed - > chrom_sizes.bed

extractCenterSites.sh -c chrom_sizes.bed -M hg38.mappable_only.bed -o hg38.center_sites.n100.starch

```

The workflow ends at footprint calling; motif matching against JASPAR is a separate

downstream step (see the `jaspar-motifs` skill).

### RGT Data Directory (footprinting)

HINT reads genome sequence and annotation from an RGT data directory, which is several GB and

is not part of the image. Create it once, then pass it with `--rgt_data`:

```bash

pip install RGT==1.0.2 # creates ~/rgtdata with setupGenomicData.py

cd ~/rgtdata && python setupGenomicData.py --hg38

```

Then run with `--rgt_data ~/rgtdata`, or use `--skip_footprint` to stop after hotspot calling.

### Blacklist Filtering

```bash

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

```

## Footprinting Analysis

```bash

rgt-hint footprinting \

--paired-end \

--organism hg38 \

--output-location footprints/ \

sample.filtered.bam \

```

### Interpretation

- Footprints are depressions in the DNase signal where a bound TF protects DNA

- Requires deep sequencing (>100M reads) for reliable footprints

- Sensitivity varies by TF: pioneer factors have shallow footprints

- Vierstra et al. 2020 provides a global reference map for comparison

## Provenance Integration

After pipeline completion, log all outputs:

```python

encode_log_derived_file(

file_path="/results/hotspots/sample1.hotspots.fdr0.05.bed",

source_accessions=["ENCSR...", "ENCFF..."],

description="DNase hypersensitive sites from ENCODE DNase-seq pipeline",

file_type="DHS_peaks",

tool_used="BWA 0.7.18 + Hotspot2 2.1.2",

parameters="FDR 0.05, blacklist filtered, ENCODE hg38 mappability index"

)

```

## Reference Files

Detailed step-by-step documentation is provided in the `references/` directory:

1. `01-qc-trimming.md` -- Read QC and adapter trimming

2. `02-alignment.md` -- BWA-MEM alignment for DNase-seq

3. `03-filtering.md` -- BAM filtering, deduplication, blacklist removal

4. `04-hotspot-calling.md` -- Hotspot2 DHS detection and signal generation

## Walkthrough: Processing ENCODE DNase-seq from FASTQ to Hypersensitive Sites

**Goal**: Process raw DNase-seq FASTQ files through the ENCODE pipeline to generate DNase I hypersensitive site (DHS) peak calls.

**Context**: DNase-seq identifies open chromatin via DNase I enzyme digestion. The pipeline uses BWA alignment and Hotspot2 for DHS identification.

### Step 1: Find DNase-seq experiment

```

encode_search_experiments(assay_title="DNase-seq", biosample_term_name="K562", organism="Homo sapiens")

```

Expected output:

```json

{

"total": 8,

"results": [

{"accession": "ENCSR000DNS", "assay_title": "DNase-seq", "biosample_summary": "K562", "status": "released"}

]

}

```

### Step 2: List and download FASTQ files

```

```

### Step 3: Run the DNase-seq pipeline

```bash

```

Key pipeline steps:

1. Quality trimming

2. BWA-MEM alignment

4. Hotspot2 DHS calling

5. Signal track generation

### Step 4: Validate output quality

| Metric | Threshold | Purpose |

|---|---|---|

| SPOT score | > 0.4 | Signal portion of tags |

| Duplicate rate | < 30% | Library complexity |

### Step 5: Compare with ATAC-seq

```

encode_search_experiments(assay_title="ATAC-seq", biosample_term_name="K562", organism="Homo sapiens")

```

**Interpretation**: DNase-seq and ATAC-seq both measure accessibility but with different biases. Compare peaks from both assays -- concordant peaks are high confidence.

### Integration with downstream skills

- DHS peaks feed into -> **accessibility-aggregation** alongside ATAC-seq peaks

- Footprint data feeds into -> **motif-analysis** for TF binding prediction

- Signal tracks feed into -> **visualization-workflow**

- Peaks integrate with -> **regulatory-elements** for cCRE classification

## Code Examples

### 1. Survey DNase-seq availability

```

```

Expected output:

```json

{

"facets": {

"organ": {"blood": 45, "brain": 30, "liver": 20, "heart": 15, "lung": 12}

}

}

```

### 2. Check for existing DHS peaks

```

```

Expected output:

```json

{

"files": [

{"accession": "ENCFF800DHS", "output_type": "peaks", "file_format": "bed narrowPeak", "file_size_mb": 1.5}

]

}

```

### 3. Track DNase-seq experiments

```

encode_track_experiment(accession="ENCSR000DNS", notes="K562 DNase-seq for accessibility comparison with ATAC-seq")

```

Expected output:

```json

{

"status": "tracked",

"accession": "ENCSR000DNS",

"notes": "K562 DNase-seq for accessibility comparison with ATAC-seq"

}

```

## Integration

| This skill produces... | Feed into... | Purpose |

|---|---|---|

| DHS peaks (narrowPeak) | **accessibility-aggregation** | Union merge with ATAC-seq peaks |

| TF footprints | **motif-analysis** | Validate motif predictions with footprint evidence |

| Signal tracks (bigWig) | **visualization-workflow** | Genome browser display |

| Accessible regions | **regulatory-elements** | cCRE classification |

| DHS coordinates | **variant-annotation** | Annotate variants in hypersensitive sites |

| QC metrics | **quality-assessment** | Validate SPOT score and sensitivity |

| Pipeline parameters | **data-provenance** | Record BWA/Hotspot2 versions |

| DHS peak regions | **jaspar-motifs** | Scan accessible sites for known TF motifs |

## Related Skills

- `pipeline-guide` -- Parent skill with compute resource assessment and cloud setup

- `accessibility-aggregation` -- Aggregate DHS data across samples/tissues

- `quality-assessment` -- Evaluate pipeline output quality metrics

- `data-provenance` -- Track all pipeline inputs, outputs, and parameters

- `download-encode` -- Download ENCODE DNase-seq FASTQ files for pipeline input

- `publication-trust` -- Verify literature claims backing analytical decisions

## Presenting Results

When reporting DNase-seq pipeline results:

- **Hotspot counts**: Report total Hotspot2 DHS calls at the specified FDR threshold and the number remaining after blacklist filtering

- **Signal-to-noise (SPOT score)**: Report the SPOT score prominently (>0.4 pass, 0.2-0.4 warning, <0.2 fail). This is the DNase-seq equivalent of FRiP

- **Output paths**: Provide paths to hotspot BED files, narrowPeak files, signal bigWig tracks, and footprint results

- **Mappability note**: Confirm which Hotspot2 mappability index was used and that it matches the read length

- **Next steps**: Suggest `motif-analysis` for TF motif enrichment in DHS peaks, or `accessibility-aggregation` for merging DHS data across samples

## For the request: "$ARGUMENTS"