# QC and Trimming for DNase-seq Data

DNase-seq reads are standard Illumina paired-end reads. Quality control and
trimming follow conventional practices with attention to fragment size.

## Pre-Trimming QC with FastQC

```bash
fastqc --threads 4 --outdir fastqc_raw/ sample_R1.fastq.gz sample_R2.fastq.gz
```

Key checks:
- Per-base quality (expect Phred >28)
- Adapter content (Illumina universal adapter)
- Sequence duplication (DNase-seq libraries can show moderate duplication)
- GC content (should match genome, ~40% for human)

## Adapter Trimming with Trim Galore

```bash
trim_galore \
    --paired \
    --quality 20 \
    --phred33 \
    --length 20 \
    --cores 4 \
    --fastqc \
    sample_R1.fastq.gz \
    sample_R2.fastq.gz
```

`--fastqc` runs FastQC on the trimmed reads. The workflow publishes the trimmed
FASTQs, the trimming reports and those trimmed-read FastQC reports together in
`results/trim_galore/`, and feeds the reports to MultiQC.

### Parameter Rationale

| Parameter | Value | Reason |
|-----------|-------|--------|
| `--quality 20` | Phred 20 | Standard quality cutoff |
| `--length 20` | 20 bp | DNase-seq can have short fragments; keep short reads |
| `--cores 4` | 4 | Parallel processing |

DNase-seq fragments can be very short (sub-nucleosomal), so the minimum
read length after trimming is set lower than typical (20 bp vs 36 bp).

## Post-Trimming Verification

After trimming, verify:
- >95% of reads pass quality filter
- Adapter contamination removed
- Read length distribution shows reads from 20-150 bp

## Fragment Size Distribution Check

After alignment, verify the insert size distribution. The workflow uses
`samtools stats`, which reports the same distribution and is read by MultiQC:

```bash
samtools stats sample.filtered.bam > sample.insert_sizes.txt

# The insert-size histogram is the IS block of that file
grep ^IS sample.insert_sizes.txt | cut -f 2,3
```

Picard `CollectInsertSizeMetrics` is an alternative, but it always writes a
histogram PDF and therefore requires R, which the pipeline image does not ship.
Run it manually in an environment that has R:

```bash
picard CollectInsertSizeMetrics \
    INPUT=sample.filtered.bam \
    OUTPUT=insert_sizes.txt \
    HISTOGRAM_FILE=insert_size_histogram.pdf \
    MINIMUM_PCT=0.05
```

Expected DNase-seq fragment sizes:
- Primary peak: 50-100 bp (sub-nucleosomal, DHS fragments)
- Secondary peak: ~170 bp (mononucleosomal)
- Tail extending to 500+ bp
- If no sub-nucleosomal peak, library may have size-selection issues

## Single-End vs Paired-End

**Not supported by this workflow -- for reference only.** The Nextflow pipeline
is paired-end only: it always passes `trim_galore --paired`, aligns read 1 and
read 2 together, and filters on properly paired reads. There is no
`--single_end` parameter. Single-end FASTQs have to be processed outside the
pipeline.

ENCODE DNase-seq data exists in both SE and PE formats. For single-end data,
trim with:
```bash
trim_galore \
    --quality 20 \
    --phred33 \
    --length 20 \
    --cores 4 \
    --fastqc \
    sample.fastq.gz
```

Paired-end is preferred for:
- Better duplicate detection
- Insert size QC
- Fragment-level analysis
- Footprinting (requires properly paired reads)

# BWA-MEM Alignment for DNase-seq

DNase-seq uses standard BWA-MEM alignment with paired-end reads.
Unlike Hi-C, DNase-seq reads come from contiguous fragments and can be
aligned as normal paired-end data.

## Genome Index Preparation

```bash
# Build BWA index (one-time, ~1 hour for human genome)
bwa index -a bwtsw genome.fa
```

Requires ~8 GB disk space for the human genome.

`--bwa_index` is the **prefix** the index was built on, i.e. the FASTA path
(`/ref/bwa_index/genome.fa`), not a directory. The workflow stages every file
matching `<prefix>*` -- the FASTA plus `.amb`, `.ann`, `.bwt`, `.pac`, `.sa` --
into the task directory and passes the basename to `bwa mem`.

## Paired-End Alignment

```bash
bwa mem -t 8 -M \
    genome.fa \
    sample_R1_val_1.fq.gz \
    sample_R2_val_2.fq.gz \
    | samtools view -@ 4 -bS - \
    | samtools sort -@ 4 -o sample_sorted.bam

samtools index sample_sorted.bam
```

### BWA-MEM Parameters

| Parameter | Value | Reason |
|-----------|-------|--------|
| `-t 8` | 8 threads | Parallel alignment |
| `-M` | Mark shorter split as secondary | Picard compatibility |

Unlike ChIP-seq or ATAC-seq, DNase-seq does not need special alignment flags.
Standard BWA-MEM settings work well.

## Single-End Alignment

**Not supported by this workflow -- for reference only.** The Nextflow pipeline
is paired-end only and has no `--single_end` parameter. For older ENCODE
DNase-seq datasets with single-end reads, align them outside the pipeline:

```bash
bwa mem -t 8 -M \
    genome.fa \
    sample_trimmed.fq.gz \
    | samtools view -@ 4 -bS - \
    | samtools sort -@ 4 -o sample_sorted.bam

samtools index sample_sorted.bam
```

## Alignment QC

```bash
samtools flagstat sample_sorted.bam > sample_flagstat.txt
```

Expected metrics:
- **Mapping rate**: >80% (typically 85-95%)
- **Properly paired**: >90% of mapped reads
- **Supplementary**: <5%
- **Secondary**: <5%

### MAPQ Distribution

```bash
samtools view sample_sorted.bam | \
    awk '{print $5}' | sort -n | uniq -c | sort -rn | head -20
```

Most reads should have MAPQ >= 30. A large fraction of MAPQ 0 reads suggests
the sample has high repeat content or contamination.

## Mitochondrial Reads

DNase-seq can contain significant mitochondrial DNA (chrM reads):

```bash
# Count mitochondrial reads
samtools idxstats sample_sorted.bam | awk '$1 == "chrM" {print "chrM reads:", $3}'
total=$(samtools view -c sample_sorted.bam)
chrm=$(samtools view -c sample_sorted.bam chrM)
echo "chrM fraction: $(echo "scale=4; $chrm / $total" | bc)"
```

Mitochondrial fraction:
- <5%: Normal
- 5-20%: Elevated but acceptable
- >20%: Poor nuclear enrichment

chrM reads are removed in the filtering step.

## Memory and Time Estimates

| Read Count | Threads | RAM | Time |
|-----------|---------|-----|------|
| 50M reads | 8 | 12 GB | 30-60 min |
| 100M reads | 8 | 12 GB | 1-2 hours |
| 200M reads | 8 | 12 GB | 2-4 hours |

BWA-MEM memory usage is dominated by index size (~8 GB for human genome),
not read count. RAM requirements are relatively stable.

## Alignment Considerations for Footprinting

If footprinting analysis is planned downstream, alignment quality is critical:
- Use paired-end data (SE data has lower footprint resolution)
- Do NOT remove soft-clipped bases (they indicate fragment boundaries)
- Higher depth is better: 100M+ reads for reliable footprints
- MAPQ 30 filter is applied downstream, not during alignment

# Filtering and Deduplication for DNase-seq

After alignment, remove low-quality reads, duplicates, mitochondrial reads,
and blacklisted regions to produce a clean BAM for Hotspot2 analysis.

## Quality Filtering

```bash
samtools view -b -h \
    -q 30 \
    -F 1804 \
    -f 2 \
    sample_sorted.bam \
    | samtools sort -@ 4 -o sample_filtered.bam
```

### Filter Flags

| Flag | Meaning |
|------|---------|
| `-q 30` | MAPQ >= 30 (uniquely mapped) |
| `-F 4` | Remove unmapped |
| `-F 256` | Remove secondary alignments |
| `-F 512` | Remove QC-failed reads |
| `-F 1024` | Remove duplicates (after marking) |
| `-f 2` | Keep only properly paired |

Combined `-F 1804` removes unmapped + secondary + QC-fail + duplicates.

## Remove Mitochondrial Reads

```bash
# Get list of non-chrM chromosomes
samtools idxstats sample_filtered.bam | \
    awk '$1 != "chrM" && $1 != "*" {print $1}' > chroms.txt

# Filter to nuclear chromosomes only
samtools view -b -h \
    sample_filtered.bam \
    $(cat chroms.txt | tr '\n' ' ') \
    > sample_nuclear.bam

samtools index sample_nuclear.bam
```

## Duplicate Marking with Picard

```bash
picard MarkDuplicates \
    INPUT=sample_nuclear.bam \
    OUTPUT=sample_dedup.bam \
    METRICS_FILE=sample_dup_metrics.txt \
    REMOVE_DUPLICATES=true \
    VALIDATION_STRINGENCY=LENIENT \
    ASSUME_SORTED=true

samtools index sample_dedup.bam
```

### Library Complexity Metrics

From the Picard output, extract:

```bash
grep -A 1 'LIBRARY' sample_dup_metrics.txt | tail -1 | \
    awk '{
        print "Total pairs:", $3;
        print "Unique pairs:", $3 - $7;
        print "Duplication rate:", $9;
        print "Estimated library size:", $10
    }'
```

## Blacklist Removal

Remove reads overlapping ENCODE blacklist regions:

```bash
# Download blacklist
wget -q https://github.com/Boyle-Lab/Blacklist/raw/master/lists/hg38-blacklist.v2.bed.gz
gunzip hg38-blacklist.v2.bed.gz

# Remove blacklisted reads
bedtools intersect \
    -a sample_dedup.bam \
    -b hg38-blacklist.v2.bed \
    -v \
    > sample_final.bam

samtools index sample_final.bam
```

## NRF and PBC Metrics

ENCODE reports library complexity metrics:

```bash
# NRF = Non-Redundant Fraction = unique / total
# PBC1 = PCR Bottleneck Coefficient 1 = 1-position / 1-position-or-more

samtools view -F 1804 -f 2 sample_sorted.bam | \
    awk 'BEGIN{OFS="\t"} {print $3, $4, $3, $4}' | \
    sort | uniq -c | \
    awk 'BEGIN{mt=0; m0=0; m1=0; m2=0}
    {mt+=$1; m0++; if($1==1) m1++; if($1==2) m2++}
    END {
        print "NRF:", m0/mt;
        print "PBC1:", m1/m0;
        print "PBC2:", m1/m2
    }'
```

| Metric | Pass | Warning | Fail |
|--------|------|---------|------|
| NRF | >0.8 | 0.7-0.8 | <0.7 |
| PBC1 | >0.9 | 0.7-0.9 | <0.7 |
| PBC2 | >3.0 | 1.0-3.0 | <1.0 |

## Final BAM Statistics

```bash
samtools flagstat sample_final.bam > sample_final_flagstat.txt
echo "Final read count: $(samtools view -c sample_final.bam)"
```

# Hotspot2 DHS Calling and Signal Generation

Hotspot2 is the ENCODE-standard peak caller for DNase-seq data. It identifies
DNase I Hypersensitive Sites (DHSs) using a local tag density model that
accounts for mappability variation across the genome.

## Hotspot2 Execution

`hotspot2.sh` takes options followed by two positional arguments: the input BAM
and the output directory. There is no `-s` (input) or `-o` (output) option.
`-c` expects the chromosome sizes as a sorted BED file (chromosome, 0, length),
not the two-column `.chrom.sizes` file, so convert it first. This is exactly what
`scripts/main.nf` runs:

```bash
# hotspot2.sh wants chromosome sizes as a sorted BED file with column 2 set to 0
awk 'BEGIN {OFS="\t"} {print $1, 0, $2}' hg38.chrom.sizes | sort-bed - > chrom_sizes.bed

hotspot2.sh \
    -c chrom_sizes.bed \
    -C /ref/hotspot2/hg38.center_sites.n100.starch \
    -M /ref/hotspot2/hg38.mappable_only.bed \
    -f 0.05 \
    -F 0.05 \
    sample.filtered.bam \
    hotspot2_out
```

### Key Parameters

| Parameter | Value | Description |
|-----------|-------|-------------|
| `-c` | `chrom_sizes.bed` | Required. Chromosome sizes as a sorted BED file: chromosome, 0, length |
| `-C` | `center_sites.starch` | Required. Center sites, made once per genome by `extractCenterSites.sh` |
| `-M` | `mappable_regions.bed` | Optional. The mappable-regions BED the center sites were derived from (read-length specific) |
| `-f` | 0.05 | Hotspot FDR threshold. Every output file is named after this value |
| `-F` | 0.05 | Site-call FDR threshold. Must be >= `-f`, or `hotspot2.sh` aborts |
| `-p` | default | Peak definition (peak shape), not a protocol name. The workflow leaves it at the default |
| positional 1 | `in.bam` | Input BAM |
| positional 2 | `outdir` | Output directory |

`-F` may not be stricter than `-f`. Because `--fdr` is user-supplied, `main.nf`
passes `max(--fdr, 0.05)` for `-F` so that a loose `--fdr` cannot abort the run.
If you call `hotspot2.sh` by hand with a variable FDR, do the same:

```bash
SITECALL_FDR=$(echo "$FDR 0.05" | awk '{print ($1 > $2) ? $1 : $2}')
hotspot2.sh -c chrom_sizes.bed -C center_sites.starch -f "$FDR" -F "$SITECALL_FDR" in.bam outdir
```

### Center Sites and Mappability

Hotspot2 needs a center-sites archive, which is derived from a mappable-regions
BED. Both are genome-build and read-length specific. There is no
`hotspot2-mappability` program; the script that ships with Hotspot2 is
`extractCenterSites.sh`, and it is on the PATH inside the pipeline image:

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

Mappable-regions BEDs are produced separately (for example with Umap/Bismap or
the ENCODE-provided files) for the read length of the library. Pass the same
file to `hotspot2.sh -M` that was used to build the center sites.

## Output Files

`hotspot2.sh` names every output after the input BAM's basename and the `-f`
threshold, never `-F`. For `sample.filtered.bam` with `-f 0.05`:

| File | Format | Description |
|------|--------|-------------|
| `sample.filtered.hotspots.fdr0.05.starch` | starch | Hotspot regions passing the FDR cutoff |
| `sample.filtered.peaks.fdr0.05.starch` | starch | Individual peaks within hotspots |
| `sample.filtered.peaks.fdr0.05.narrowpeaks.starch` | starch | The same peaks in narrowPeak columns |
| `sample.filtered.SPOT.fdr0.05.txt` | text | SPOT score |
| `sample.filtered.allcalls.starch` | starch | All site calls before FDR filtering |
| `sample.filtered.density.starch` | starch | Per-base cleavage density |
| `sample.filtered.density.bw` | bigWig | The same cleavage density as a bigWig |

Every `.starch` file is a compressed BEDOPS archive, not text. Run `unstarch`
before feeding it to `bedtools`, `awk`, or anything else that expects BED.

## Convert to BED and narrowPeak

```bash
unstarch hotspot2_out/sample.filtered.hotspots.fdr0.05.starch > sample.hotspots.fdr0.05.bed
unstarch hotspot2_out/sample.filtered.peaks.fdr0.05.narrowpeaks.starch > sample.peaks.narrowPeak
unstarch hotspot2_out/sample.filtered.allcalls.starch > sample.allcalls.bed
cp hotspot2_out/sample.filtered.SPOT.fdr0.05.txt sample.SPOT.txt
```

These four files are what the workflow publishes to `results/hotspots/`.

The workflow filters the BAM against the blacklist before Hotspot2 runs, so the
peaks are already blacklist-free. A peak-level intersect is only needed if you
call hotspots on a BAM that was not filtered, or if you want to apply an
additional list:

```bash
bedtools intersect -a sample.peaks.narrowPeak -b /ref/hg38-blacklist.v2.bed -v \
    > sample.DHS.narrowPeak
```

## SPOT Score

The SPOT score (Signal Portion of Tags) is computed by Hotspot2:

```bash
cat sample.SPOT.txt
```

If computing manually from the converted BED:
```bash
total_tags=$(samtools view -c sample.filtered.bam)
tags_in_hotspots=$(bedtools intersect \
    -a sample.filtered.bam \
    -b sample.hotspots.fdr0.05.bed \
    -u -bed | wc -l)
echo "SPOT score: $(echo "scale=4; $tags_in_hotspots / $total_tags" | bc)"
```

## Signal Track Generation

Hotspot2 writes its own cleavage-density bigWig (`*.density.bw`, above). The
workflow does not publish it; it builds a separate read-per-million fragment
coverage track with bedtools and publishes that as
`results/hotspots/<sample>.density.bw`:

```bash
# RPM scale factor from the filtered, properly paired reads
total=$(samtools view -c -F 1804 -f 2 sample.filtered.bam)
scale=$(echo "scale=10; 1000000 / $total" | bc)

bedtools genomecov \
    -ibam sample.filtered.bam \
    -bg \
    -pc \
    -g /ref/hg38.chrom.sizes \
    | awk -v s=$scale 'BEGIN{OFS="\t"} {$4=$4*s; print}' \
    | sort -k1,1 -k2,2n \
    > sample_rpm.bedGraph

bedGraphToBigWig sample_rpm.bedGraph /ref/hg38.chrom.sizes sample.density.bw
```

Use the Hotspot2 density track instead when you want per-base cut counts rather
than fragment coverage.

## Peak Annotation

Annotate DHSs with genomic features (a manual downstream step; the workflow
stops at peak calling):

```bash
# Count peaks by category
total=$(wc -l < sample.peaks.narrowPeak)
echo "Total DHSs: $total"

# Overlap with gene promoters (TSS +/- 2kb)
promoter=$(bedtools intersect -a sample.peaks.narrowPeak -b promoters.bed -u | wc -l)
echo "Promoter DHSs: $promoter ($(echo "scale=1; $promoter*100/$total" | bc)%)"

# Overlap with known enhancers
enhancer=$(bedtools intersect -a sample.peaks.narrowPeak -b enhancers.bed -u | wc -l)
echo "Enhancer DHSs: $enhancer ($(echo "scale=1; $enhancer*100/$total" | bc)%)"
```

# TF Footprinting Analysis for DNase-seq

Transcription factor footprinting detects protein-DNA interactions from
DNase-seq cleavage patterns. Bound TFs protect DNA from cleavage, creating
a "footprint" -- a local depression in the DNase-seq signal.

## Prerequisites

Reliable footprinting requires:
- **Deep sequencing**: >100 million uniquely mapped reads
- **Paired-end data**: Better resolution than single-end
- **High SPOT score**: >0.4 (clean signal)
- **Peak calls**: DHS regions from Hotspot2
- **An RGT data directory**: HINT reads genome sequence and annotation from
  `$RGTDATA`; the workflow stages the directory passed with `--rgt_data`

## HINT Footprinting (DNase-seq mode)

HINT (from the Regulatory Genomics Toolbox) works for both DNase-seq and
ATAC-seq data. This is what the workflow runs:

```bash
rgt-hint footprinting \
    --dnase-seq \
    --paired-end \
    --organism hg38 \
    --output-location footprints/ \
    --output-prefix sample \
    sample.filtered.bam \
    sample.peaks.narrowPeak
```

It writes `footprints/sample.bed`, which the workflow publishes as
`results/footprints/sample.footprints.bed`. No other footprinting output is
produced.

### Key Parameters

| Parameter | Value | Description |
|-----------|-------|-------------|
| `--dnase-seq` | Flag | Use DNase-seq cleavage model (not ATAC-seq) |
| `--paired-end` | Flag | Use paired-end fragment information |
| `--organism` | hg38 | Genome build for bias correction; must be set up in the RGT data directory |
| `--output-prefix` | sample | Base name of the output BED |

**Important**: Use `--dnase-seq` for DNase-seq data and `--atac-seq` for
ATAC-seq data. They have different cleavage bias models, and using the wrong
one silently produces wrong footprints.

## Motif Matching in Footprints (manual, optional)

Motif matching is not part of the workflow. To run it afterwards, note that
`rgt-motifanalysis matching --motif-dbs` takes **directories of `.pwm` files**
in the pre-2016 JASPAR format, not a `.meme` file. The RGT data directory ships
such directories under `motifs/`:

```bash
export RGTDATA=/ref/rgtdata

rgt-motifanalysis matching \
    --organism hg38 \
    --input-files results/footprints/sample.footprints.bed \
    --output-location motif_matches/ \
    --motif-dbs /ref/rgtdata/motifs/jaspar_vertebrates
```

To use a motif collection RGT does not ship, convert it to `.pwm` files in
their own directory first with RGT's `createPwm.py`, then point `--motif-dbs`
at that directory.

## Wellington Footprinting (alternative, manual)

Wellington (pyDNase) uses a different statistical approach. It is **not** in the
pipeline image or in `dnaseseq-env.yml`; install it separately with
`pip install pyDNase`.

```bash
wellington_footprints.py \
    -A \
    -p 20 \
    -fdrlimit 0.01 \
    sample.peaks.narrowPeak \
    sample.filtered.bam \
    wellington_out/
```

### HINT vs Wellington Comparison

| Feature | HINT | Wellington |
|---------|------|-----------|
| Bias correction | Sequence-specific | Position-based |
| Speed | Moderate | Fast |
| Sensitivity | Higher | More conservative |
| DNase + ATAC | Both | Both |
| Active development | Yes | Limited |
| In the pipeline image | Yes (RGT 1.0.2) | No (`pip install pyDNase`) |

## Footprint Quality Assessment

None of the checks below are computed by the workflow; run them on the
published footprint BED.

### Per-Motif Footprint Depth

```bash
# Calculate average footprint score at known CTCF sites
bedtools intersect \
    -a CTCF_motif_sites.bed \
    -b results/footprints/sample.footprints.bed \
    -wa -wb \
    | awk '{print $NF}' \
    | awk '{sum+=$1; n++} END {print "Mean CTCF footprint score:", sum/n}'
```

### Footprint vs Background Signal Ratio

Good footprints show:
- Clear signal depression at the motif center
- Flanking shoulders of higher cleavage
- Depth-to-flank ratio > 1.5

### Expected Footprint Counts

| Sequencing Depth | Expected Footprints |
|------------------|---------------------|
| 50M reads | Unreliable |
| 100M reads | 50,000-100,000 |
| 200M reads | 100,000-200,000 |
| 500M reads | 200,000-400,000 |

## Aggregate Footprint Visualization (manual)

Generate aggregate footprint profiles across all instances of a motif. This
needs the motif matches from the step above:

```bash
rgt-hint differential \
    --organism hg38 \
    --bc \
    --nc 8 \
    --mpbs-files motif_matches/sample_mpbs.bed \
    --reads-files sample.filtered.bam \
    --conditions sample \
    --output-location diff_footprints/
```

This produces per-motif aggregate profiles showing the average cleavage
pattern across all binding sites, which is more robust than individual
footprint calls.

## Comparison with Vierstra 2020 Reference Map

The Vierstra et al. 2020 reference map provides a global catalog of human
TF footprints from 243 DNase-seq datasets. Use it to:

1. Validate your footprint calls
2. Compare tissue-specific footprinting
3. Identify novel TF binding events

```bash
# Download Vierstra reference footprints
# Available at: https://www.vierstra.org/resources/dgf

# Compare overlap
bedtools intersect \
    -a results/footprints/sample.footprints.bed \
    -b vierstra_consensus_footprints.bed \
    -u | wc -l
```

# DNase-seq Pipeline — Literature References

**Last updated:** 2026-03-07
**Purpose:** Reference catalog for the pipeline-dnaseseq skill — papers defining ENCODE DNase-seq processing standards, DHS calling with Hotspot2, transcription factor footprinting, and quality assessment for DNase I hypersensitivity assays.

---

## DNase-seq Method & Standards

---

### Thurman et al. 2012 — The accessible chromatin landscape of the human genome

- **Citation:** Thurman RE, Rynes E, Humbert R, Vierstra J, Maurano MT, Haugen E, Sheffield NC, Stergachis AB, Wang H, Vernot B, Garg K, John S, Sandstrom R, Bates D, Boatman L, Canfield TK, Diegel M, Dunn D, Ebersol AK, Frum T, Giste E, Johnson AK, Johnson EM, Kutyavin T, Laber B, Lee K, Lotakis D, Neph S, Neri F, Nguyen ED, Qu H, Reynolds AP, Roach V, Safi A, Sanchez ME, Sanyal A, Shafer A, Simon JM, Song L, Vong S, Weaver M, Yan Y, Zhang Z, Zhang Z, Lenhard B, Tewari M, Dorschner MO, Hansen RS, Navas PA, Stamatoyannopoulos G, Iyer VR, Lieb JD, Sunyaev SR, Akey JM, Sabo PJ, Kaul R, Furey TS, Dekker J, Crawford GE, Stamatoyannopoulos JA. 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:** Generated the most comprehensive map of DNase I hypersensitive sites (DHSs) across 125 diverse human cell and tissue types as part of ENCODE Phase 2. Identified ~2.9 million unique DHSs covering ~40% of the genome, with individual cell types having 100,000-200,000 DHSs. Demonstrated that distal DHSs (enhancers) are highly cell-type-specific while promoter DHSs are more shared. Established that DHSs mark virtually all classes of cis-regulatory elements and that 95% of the genome lies within 10 kb of a DHS in at least one cell type. This paper validated DNase-seq as a primary assay for the ENCODE encyclopedia.

---

### John et al. 2011 — Hotspot2: chromatin accessibility peak calling

- **Citation:** John S, Sabo PJ, Thurman RE, Sung MH, Biddie SC, Johnson TA, Hager GL, Stamatoyannopoulos JA. Chromatin accessibility pre-determines glucocorticoid receptor binding patterns. Nature Genetics, 43(3):264-268, 2011.
- **DOI:** [10.1038/ng.759](https://doi.org/10.1038/ng.759)
- **PMID:** 21258342 | **PMC:** PMC3049959
- **Citations:** ~600
- **Key findings:** Introduced the Hotspot algorithm (predecessor to Hotspot2) for identifying statistically significant regions of chromatin accessibility from DNase-seq data. The algorithm uses a local background model that accounts for mappability variation across the genome — critical because DNase I cuts accessible chromatin regardless of whether fragments map uniquely. Hotspot2 (the current version used in ENCODE) extends this with improved FDR control and peak refinement. Unlike MACS2 (designed for ChIP-seq), Hotspot2 is specifically calibrated for the diffuse signal pattern of DNase-seq data and is the ENCODE standard for DNase-seq peak calling.

---

### Vierstra et al. 2020 — Global reference mapping of TF footprints

- **Citation:** Vierstra J, Lazar J, Sandstrom R, Halow J, Lee K, Bates D, Diegel M, Dunn D, Neri F, Haugen E, Rynes E, Reynolds A, Nelson J, Johnson A, Frerker M, Buckley M, Kaul R, Meuleman W, Stamatoyannopoulos JA. Global reference mapping of human transcription factor footprints. Nature, 583(7818):729-736, 2020.
- **DOI:** [10.1038/s41586-020-2528-x](https://doi.org/10.1038/s41586-020-2528-x)
- **PMID:** 32728250 | **PMC:** PMC7410830
- **Citations:** ~600
- **Key findings:** Created a global reference map of TF footprints from 243 human cell types using deep DNase-seq (>200M reads per sample). Identified >4.5 million unique TF footprints genome-wide, demonstrating that ~2% of the genome is occupied by TFs at any given time. Established the analytical framework for digital genomic footprinting: bias-corrected DNase cleavage profiles around TF motifs, with bound sites showing a characteristic protection pattern (footprint) flanked by elevated cleavage. Published as part of ENCODE Phase 3, providing the reference standard for interpreting TF footprints from DNase-seq data.

---

### 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 DNase-seq as a primary accessibility assay alongside ATAC-seq for identifying candidate cis-Regulatory Elements (cCREs). DNase-seq DHSs are a core input for the cCRE classification system, with DHS signal strength and overlap with histone marks determining element type (promoter-like, enhancer-like, CTCF-only, K4m3). Established uniform processing standards: BWA-MEM alignment, Hotspot2 peak calling, SPOT score as the primary QC metric (>0.4 for acceptable quality), and GRCh38/mm10 as reference assemblies.

---

### Hitz et al. 2023 — The ENCODE Uniform Analysis Pipelines

- **Citation:** Hitz BC, Lee JW, Jolanki O, et al. The ENCODE Uniform Analysis Pipelines. Research Square (preprint), 2023.
- **DOI:** [10.21203/rs.3.rs-311127/v1](https://doi.org/10.21203/rs.3.rs-311127/v1)
- **Citations:** ~84
- **Key findings:** Documents the official ENCODE DNase-seq pipeline specification: BWA-MEM alignment, Picard MarkDuplicates, Hotspot2 for DHS calling with FDR control, and SPOT score computation. Specifies that DNase-seq uses Hotspot2 (not MACS2) because Hotspot2's mappability-aware background model is critical for the DNase-seq signal profile. Pipeline infrastructure uses Docker containers and WDL for reproducible execution.

---

## TF Footprinting

---

### Li et al. 2019 — HINT-ATAC: TF footprinting from chromatin accessibility

- **Citation:** Li Z, Schulz MH, Look T, Begus M, Zenke M, Costa IG. Identification of transcription factor binding sites using ATAC-seq. Genome Research, 29(11):1850-1858, 2019.
- **DOI:** [10.1186/s13059-019-1642-2](https://doi.org/10.1186/s13059-019-1642-2)
- **PMID:** 31649060
- **Citations:** ~300
- **Key findings:** Introduced HINT-ATAC, a computational framework for TF footprinting from chromatin accessibility data (ATAC-seq and DNase-seq). HINT-ATAC models the enzymatic cleavage bias of both DNase I and Tn5 transposase, providing bias-corrected footprint scores that improve over raw cleavage signal. Supports both DNase-seq and ATAC-seq data with enzyme-specific bias correction models. Used in this pipeline's optional footprinting analysis step (`rgt-hint footprinting --dnase-seq`), producing a footprint BED that can be compared against the Vierstra et al. 2020 reference atlas.

---

### Neph et al. 2012 — An expansive human regulatory lexicon encoded in TF footprints

- **Citation:** Neph S, Vierstra J, Stergachis AB, Reynolds AP, Haugen E, Vernot B, Thurman RE, John S, Sandstrom R, Johnson AK, Maurano MT, Humbert R, Rynes E, Wang H, Vong S, Lee K, Bates D, Diegel M, Roach V, Dunn D, Neri J, Schafer A, Hansen RS, Kutyavin T, Giste E, Weaver M, Canfield T, Sabo P, Zhang M, Balasundaram G, Byron R, MacCoss MJ, Akey JM, Bender MA, Groudine M, Kaul R, Stamatoyannopoulos JA. An expansive human regulatory lexicon encoded in transcription factor footprints. Nature, 489(7414):83-90, 2012.
- **DOI:** [10.1038/nature11212](https://doi.org/10.1038/nature11212)
- **PMID:** 22955618 | **PMC:** PMC3736026
- **Citations:** ~1,200
- **Key findings:** First genome-wide census of TF footprints from DNase-seq data across 41 diverse human cell types. Identified ~8.4 million distinct TF occupancy events involving 475 sequence motifs. Demonstrated that digital genomic footprinting from deep DNase-seq can detect individual TF binding events with nucleotide precision. Revealed regulatory motif co-occurrence patterns and cell-type-specific regulatory lexicons. Established the computational methodology for DNase-seq footprinting that forms the basis for this pipeline's footprinting analysis.

---

## Supplementary Tools (Non-DNase-Specific)

See pipeline-chipseq/references/literature.md for detailed descriptions of shared tools.

---

### Li & Durbin 2009 — BWA

- **DOI:** [10.1093/bioinformatics/btp324](https://doi.org/10.1093/bioinformatics/btp324) | **PMID:** 19451168 | **Citations:** ~25,000
- **DNase-seq role:** Read alignment with BWA-MEM. Same alignment approach as ChIP-seq. MAPQ filtering removes multi-mapped reads that could create false DHS calls.

---

### Li et al. 2009 — SAMtools

- **DOI:** [10.1093/bioinformatics/btp352](https://doi.org/10.1093/bioinformatics/btp352) | **PMID:** 19505943 | **Citations:** ~53,700
- **DNase-seq role:** BAM sorting, filtering, and indexing. Used for MAPQ filtering, coordinate sorting, and computing alignment statistics.

---

### Broad Institute — Picard MarkDuplicates

- **URL:** [https://broadinstitute.github.io/picard/](https://broadinstitute.github.io/picard/)
- **DNase-seq role:** PCR duplicate marking and library complexity estimation. NRF, PBC1, PBC2 metrics computed from Picard output. DNase-seq libraries can have moderate duplication rates due to the enzymatic digestion step.

---

### Quinlan & Hall 2010 — BEDTools

- **DOI:** [10.1093/bioinformatics/btq033](https://doi.org/10.1093/bioinformatics/btq033) | **PMID:** 20110278 | **Citations:** ~12,000
- **DNase-seq role:** Used for blacklist region filtering, SPOT score calculation (reads in hotspots), and peak annotation with gene models.

---

### Amemiya et al. 2019 — ENCODE Blacklist

- **DOI:** [10.1038/s41598-019-45839-z](https://doi.org/10.1038/s41598-019-45839-z) | **PMID:** 31249361 | **Citations:** ~1,372
- **DNase-seq role:** Blacklist v2 filtering is essential for DNase-seq because DNase I can preferentially cut accessible regions in blacklisted repetitive elements, creating artifact peaks.

---

### Ewels et al. 2016 — MultiQC

- **DOI:** [10.1093/bioinformatics/btw354](https://doi.org/10.1093/bioinformatics/btw354) | **PMID:** 27312411 | **Citations:** ~6,800
- **DNase-seq role:** Aggregates QC metrics from FastQC, BWA alignment, Picard deduplication, and Hotspot2 into a unified HTML report for batch-level quality assessment.

---

### Andrews 2010 — FastQC

- **URL:** [https://www.bioinformatics.babraham.ac.uk/projects/fastqc/](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
- **DNase-seq role:** Raw read quality assessment. DNase-seq libraries typically show characteristic fragment size distributions with a peak at sub-nucleosomal sizes (50-100 bp).

---

### Martin 2011 — Cutadapt (basis for Trim Galore)

- **DOI:** [10.14806/ej.17.1.200](https://doi.org/10.14806/ej.17.1.200) | **Citations:** ~13,000
- **DNase-seq role:** Adapter trimming with Trim Galore. Important for DNase-seq because short sub-nucleosomal fragments frequently read through into adapters, similar to ATAC-seq.