Execute the ENCODE DNase-seq pipeline for chromatin accessibility profiling,\nproducing DNase hypersensitive sites (DHSs) via Hotspot2 and transcription\nfactor footprints.
\nFASTQ -> Trim -> BWA-MEM align -> Filter/dedup -> Hotspot2 -> DHS peaks\n | |\n Signal track Footprinting (HINT)\nENCODE-DCC/dnase-seq-pipelinescripts/Dockerfile in this skill (docker build -t encode-toolkit/pipeline-dnaseseq:1.0.0 scripts/); override with --container| Tool | \nVersion | \nPurpose | \nCitation | \n
|---|---|---|---|
| BWA-MEM | \n0.7.18 | \nAlignment | \nLi & Durbin 2009 | \n
| samtools | \n1.19 | \nBAM operations | \nLi et al. 2009 | \n
| Picard | \n3.1.1 | \nDuplicate marking | \nBroad Institute | \n
| Hotspot2 | \n2.1.2 | \nDHS calling (ENCODE standard) | \nJohn et al. 2011 | \n
| modwt | \n1.0 | \nWavelet smoothing used by Hotspot2 | \nStam Lab | \n
| bedtools | \n2.31.0 | \nGenomic arithmetic | \nQuinlan & Hall 2010 | \n
| BEDOPS | \napt (Ubuntu 22.04) | \nsort-bed and unstarch for the Hotspot2 starch archives | \nNeph et al. 2012 | \n
| HINT (RGT) | \n1.0.2 | \nTF footprinting | \nLi et al. 2019 | \n
| FastQC | \n0.12.1 | \nRead quality | \nAndrews (Babraham) | \n
| Trim Galore | \n0.6.10 | \nAdapter and quality trimming | \nKrueger (Babraham) | \n
| cutadapt | \n4.6 | \nAdapter removal backend for Trim Galore | \nMartin 2011 | \n
| MultiQC | \n1.21 | \nAggregated QC | \nEwels et al. 2016 | \n
John et al. 2011 - \"Chromatin accessibility pre-determines glucocorticoid\nreceptor binding patterns\" (Nature Genetics, ~600 citations)\nDOI: 10.1038/ng.759
\nThurman et al. 2012 - \"The accessible chromatin landscape of the human\ngenome\" (Nature, ~3,000 citations)\nDOI: 10.1038/nature11232
\nVierstra et al. 2020 - \"Global reference mapping of human transcription\nfactor footprints\" (Nature, ~600 citations)\nDOI: 10.1038/s41586-020-2528-x
\nAmemiya et al. 2019 - \"The ENCODE Blacklist\" (Scientific Reports, ~1,372 citations)\nDOI: 10.1038/s41598-019-45839-z
\nLi et al. 2019 - \"Identification of transcription factor binding sites using\nATAC-seq\" (Genome Biology) -- HINT-ATAC footprinting\nDOI: 10.1186/s13059-019-1642-2
\nnextflow run scripts/main.nf \\\n -profile local \\\n --reads '/data/fastq/*_R{1,2}.fastq.gz' \\\n --bwa_index '/ref/bwa_index/genome.fa' \\\n --chrom_sizes '/ref/hg38.chrom.sizes' \\\n --hotspot_center_sites '/ref/hotspot2/hg38.center_sites.n100.starch' \\\n --hotspot_mappable '/ref/hotspot2/hg38.mappable_only.bed' \\\n --rgt_data '/ref/rgtdata' \\\n --blacklist '/ref/hg38-blacklist.v2.bed' \\\n --outdir results/ \\\n -resume\nDrop --rgt_data and add --skip_footprint to stop after hotspot calling.
The slurm profile runs through Singularity, which cannot resolve the default\nDocker image name, so pass the converted .sif with --container:
singularity build pipeline-dnaseseq.sif docker-daemon://encode-toolkit/pipeline-dnaseseq:1.0.0\n\nnextflow run scripts/main.nf \\\n -profile slurm \\\n --container /path/to/pipeline-dnaseseq.sif \\\n --reads '/data/fastq/*_R{1,2}.fastq.gz' \\\n --bwa_index '/ref/bwa_index/genome.fa' \\\n --chrom_sizes '/ref/hg38.chrom.sizes' \\\n --hotspot_center_sites '/ref/hotspot2/hg38.center_sites.n100.starch' \\\n --hotspot_mappable '/ref/hotspot2/hg38.mappable_only.bed' \\\n --rgt_data '/ref/rgtdata' \\\n --blacklist '/ref/hg38-blacklist.v2.bed' \\\n --outdir results/ \\\n -resume\n# Google Cloud Batch\nnextflow run scripts/main.nf -profile gcp \\\n --container us-docker.pkg.dev/<project>/<repo>/pipeline-dnaseseq:1.0.0 \\\n --gcp_project <project> \\\n --gcp_workdir gs://<bucket>/work \\\n --reads 'gs://<bucket>/fastq/*_R{1,2}.fastq.gz' \\\n --bwa_index gs://<bucket>/ref/bwa_index/genome.fa \\\n --chrom_sizes gs://<bucket>/ref/hg38.chrom.sizes \\\n --hotspot_center_sites gs://<bucket>/ref/hotspot2/hg38.center_sites.n100.starch \\\n --hotspot_mappable gs://<bucket>/ref/hotspot2/hg38.mappable_only.bed \\\n --rgt_data gs://<bucket>/ref/rgtdata \\\n --blacklist gs://<bucket>/ref/hg38-blacklist.v2.bed \\\n --outdir gs://<bucket>/results\n\n# AWS Batch\nnextflow run scripts/main.nf -profile aws \\\n --container <account>.dkr.ecr.<region>.amazonaws.com/pipeline-dnaseseq:1.0.0 \\\n --aws_queue <job-queue> \\\n --aws_workdir s3://<bucket>/work \\\n --reads 's3://<bucket>/fastq/*_R{1,2}.fastq.gz' \\\n --bwa_index s3://<bucket>/ref/bwa_index/genome.fa \\\n --chrom_sizes s3://<bucket>/ref/hg38.chrom.sizes \\\n --hotspot_center_sites s3://<bucket>/ref/hotspot2/hg38.center_sites.n100.starch \\\n --hotspot_mappable s3://<bucket>/ref/hotspot2/hg38.mappable_only.bed \\\n --rgt_data s3://<bucket>/ref/rgtdata \\\n --blacklist s3://<bucket>/ref/hg38-blacklist.v2.bed \\\n --outdir s3://<bucket>/results\n--outdir only sets where results are published; Google Batch and AWS Batch\nstage every task through the work directory, and the workflow stops with an\nerror if it or the project/queue is missing.
| Step | \nCPUs | \nRAM | \nTime (per sample) | \n
|---|---|---|---|
| BWA-MEM align | \n8 | \n16 GB | \n1-2 hours | \n
| Filter/dedup | \n4 | \n8 GB | \n30-60 min | \n
| Hotspot2 | \n4 | \n8 GB | \n30-60 min | \n
| Signal generation | \n2 | \n4 GB | \n15-30 min | \n
| Footprinting | \n4 | \n8 GB | \n1-2 hours | \n
| Total | \n8 | \n16 GB | \n3-6 hours | \n
Every process asks for memory { N.GB * task.attempt }, so a task killed for running out\nof memory is retried with more: the second attempt gets twice the figure in the table, the\nthird three times it, bounded by --max_memory (32 GB by default). nextflow.config\nscales the time the same way for BWA_ALIGN, FILTER_DEDUP, HOTSPOT2 and\nFOOTPRINTING, bounded by --max_time; the other processes declare no time limit. A task\nis retried only for exit codes 130-145 and 104 (killed for exceeding a limit); any other\nfailure stops the run.
| Parameter | \nDefault | \nDescription | \n
|---|---|---|
--reads | \nrequired | \nGlob pattern to paired FASTQ files (e.g. '/data/*_R{1,2}.fastq.gz'). Paired-end only | \n
--bwa_index | \nrequired | \nBWA index prefix, i.e. the FASTA path. Every file matching <prefix>* (the .fa plus .amb .ann .bwt .pac .sa) is staged | \n
--chrom_sizes | \nrequired | \nTwo-column chromosome sizes file. Used for the bigWig track and converted to BED for Hotspot2 -c | \n
--hotspot_center_sites | \nrequired | \nHotspot2 center-sites archive (.starch), made once per genome with extractCenterSites.sh (Hotspot2 -C) | \n
--hotspot_mappable | \nnull | \nMappable-regions BED that the center sites were made from (Hotspot2 -M; recommended) | \n
--blacklist | \nrequired | \nENCODE blacklist BED. Applied to the BAM before hotspot calling | \n
--outdir | \n./results | \nOutput directory | \n
--fdr | \n0.05 | \nHotspot2 hotspot FDR (-f). Names every Hotspot2 output file. -F is passed as max(--fdr, 0.05) because it may not be stricter than -f | \n
--skip_footprint | \nfalse | \nSkip footprinting analysis | \n
--organism | \nhg38 | \nGenome name registered in the RGT data directory, used by HINT footprinting | \n
--rgt_data | \nrequired unless --skip_footprint | \nRGT data directory with the genome for --organism set up (see below) | \n
There is no --genome, --single_end, or --fastq_r1/--fastq_r2 parameter.
nextflow.config)| Parameter | \nDefault | \nDescription | \n
|---|---|---|
--container | \nencode-toolkit/pipeline-dnaseseq:1.0.0 | \nImage built from scripts/Dockerfile. Pass a registry image for gcp/aws, or a .sif file for slurm | \n
--max_cpus, --max_memory, --max_time | \n16, 32.GB, 24.h | \nUpper bounds applied to every process | \n
--slurm_queue, --slurm_account | \nnormal, none | \nSLURM partition and account | \n
--gcp_project, --gcp_workdir | \nnone (both required for -profile gcp) | \nGoogle Cloud project and gs:// work directory | \n
--gcp_location, --gcp_disk | \nus-central1, 200.GB | \nGoogle Batch region and per-task disk | \n
--aws_queue, --aws_workdir | \nnone (both required for -profile aws) | \nAWS Batch job queue and s3:// work directory | \n
--aws_region, --aws_cli_path | \nus-east-1, /home/ec2-user/miniconda/bin/aws | \nAWS region, and the AWS CLI path inside the Batch AMI | \n
results/\n fastqc/ # FastQC on the raw reads\n trim_galore/\n {sample}_R1_val_1.fq.gz # Trimmed reads\n {sample}_R2_val_2.fq.gz\n *_trimming_report.txt # Trim Galore reports\n *_fastqc.{html,zip} # FastQC on the trimmed reads\n alignment/\n {sample}.filtered.bam # Filtered, deduplicated, blacklist-free BAM\n {sample}.filtered.bam.bai\n {sample}.flagstat.txt # samtools flagstat on the filtered BAM\n {sample}.dup_metrics.txt # Picard MarkDuplicates metrics\n hotspots/\n {sample}.hotspots.fdr0.05.bed # DHS hotspots (primary output; unstarched)\n {sample}.peaks.narrowPeak # Peaks within hotspots (unstarched)\n {sample}.allcalls.bed # All site calls before FDR filtering (unstarched)\n {sample}.SPOT.txt # SPOT score\n {sample}.density.bw # RPM fragment-coverage signal track (bigWig)\n footprints/\n {sample}.footprints.bed # TF footprints (omitted with --skip_footprint)\n qc/\n {sample}.insert_sizes.txt # samtools stats output (insert sizes in the IS block)\n multiqc/\n multiqc_report.html\n pipeline_info/\n timeline.html\n report.html\n trace.txt\nThe .bed, .narrowPeak and .SPOT.txt files under hotspots/ are the\nunstarch-ed forms of the Hotspot2 .starch archives; the raw archives stay in\nthe Nextflow work directory. {sample}.density.bw is the bedtools RPM fragment\ncoverage track, not the per-base cut-count bigWig Hotspot2 writes internally.
| Metric | \nPass | \nWarning | \nFail | \n
|---|---|---|---|
| SPOT score (Signal Portion of Tags) | \n>0.4 | \n0.2-0.4 | \n<0.2 | \n
| Hotspot count | \n>50,000 | \n20,000-50,000 | \n<20,000 | \n
| Mapping rate | \n>80% | \n60-80% | \n<60% | \n
| Duplication rate | \n<30% | \n30-50% | \n>50% | \n
| NRF (Non-Redundant Fraction) | \n>0.8 | \n0.7-0.8 | \n<0.7 | \n
| PBC1 (PCR Bottleneck Coefficient 1) | \n>0.9 | \n0.7-0.9 | \n<0.7 | \n
| Insert size peak | \n50-150 bp | \nVariable | \nAbnormal | \n
The workflow produces everything the first four rows and the last row need:\nthe SPOT score (hotspots/{sample}.SPOT.txt), the hotspot BED to count, the\nmapping and duplication rates (alignment/{sample}.flagstat.txt and\n{sample}.dup_metrics.txt), and the insert-size distribution\n(qc/{sample}.insert_sizes.txt). NRF, PBC1 and PBC2 are not computed;\nderive them manually from the alignment BAM as shown in\nreferences/03-filtering.md.
The SPOT score (Signal Portion of Tags) is the fraction of reads falling\nwithin hotspots. It is the DNase-seq equivalent of FRiP for ChIP-seq.
\nHigher SPOT = more enrichment in accessible regions = better library quality.
\nIMPORTANT: ENCODE uses Hotspot2 for DNase-seq, NOT MACS2.
\n| Feature | \nHotspot2 | \nMACS2 | \n
|---|---|---|
| Designed for | \nDNase-seq | \nChIP-seq | \n
| Background model | \nLocal tag density + mappability | \nDynamic Poisson | \n
| ENCODE standard | \nYes (DNase-seq) | \nYes (ChIP-seq/ATAC-seq) | \n
| Mappability correction | \nBuilt-in | \nNot available | \n
| Output | \nHotspots + peaks | \nPeaks only | \n
Hotspot2 accounts for mappability variation across the genome, which is\ncritical for DNase-seq because DNase I cuts accessible chromatin regardless\nof whether it is uniquely mappable.
\nThese are different assays measuring the same biology (chromatin accessibility):
\nDNase-seq produces a characteristic fragment size distribution:
\nHotspot2 needs a center-sites file, which is derived from a mappable-regions BED. Both are\nread-length and genome-build specific. Create the center sites once per genome with the script\nthat ships with Hotspot2 (it is on the PATH inside the image):
\n# chrom_sizes.bed is a BED file: chromosome, 0, length\nawk 'BEGIN{OFS=\"\\t\"} {print $1, 0, $2}' hg38.chrom.sizes | sort-bed - > chrom_sizes.bed\nextractCenterSites.sh -c chrom_sizes.bed -M hg38.mappable_only.bed -o hg38.center_sites.n100.starch\nPass the same mappable-regions BED to the workflow as --hotspot_mappable that\nwas used to build the center sites.
Mappable-regions files:
\nThe workflow ends at footprint calling; motif matching against JASPAR is a separate\ndownstream step (see the jaspar-motifs skill).
HINT reads genome sequence and annotation from an RGT data directory, which is several GB and\nis not part of the image. Create it once, then pass it with --rgt_data:
pip install RGT==1.0.2 # creates ~/rgtdata with setupGenomicData.py\ncd ~/rgtdata && python setupGenomicData.py --hg38\nThen run with --rgt_data ~/rgtdata, or use --skip_footprint to stop after hotspot calling.
--blacklist is required and the workflow removes blacklisted reads from the\nBAM (Amemiya et al. 2019) before Hotspot2 runs, so the published peaks are\nalready blacklist-free. Blacklist regions produce artifactual signal in\naccessibility assays.
Filter at the peak level only when the peaks came from a BAM that was not\nfiltered, or when applying an additional list:
\nbedtools intersect -a hotspots.bed -b hg38-blacklist.v2.bed -v > hotspots_filtered.bed\nTranscription factor footprinting detects bound TFs from DNase-seq signal.\nThis is what the workflow runs:
\nrgt-hint footprinting \\\n --dnase-seq \\\n --paired-end \\\n --organism hg38 \\\n --output-location footprints/ \\\n --output-prefix sample \\\n sample.filtered.bam \\\n sample.peaks.narrowPeak\nUse --dnase-seq, not --atac-seq: the two apply different cleavage-bias\nmodels, and the wrong one silently produces wrong footprints.
After pipeline completion, log all outputs:
\nencode_log_derived_file(\n file_path=\"/results/hotspots/sample1.hotspots.fdr0.05.bed\",\n source_accessions=[\"ENCSR...\", \"ENCFF...\"],\n description=\"DNase hypersensitive sites from ENCODE DNase-seq pipeline\",\n file_type=\"DHS_peaks\",\n tool_used=\"BWA 0.7.18 + Hotspot2 2.1.2\",\n parameters=\"FDR 0.05, blacklist filtered, ENCODE hg38 mappability index\"\n)\nDetailed step-by-step documentation is provided in the references/ directory:
01-qc-trimming.md -- Read QC and adapter trimming02-alignment.md -- BWA-MEM alignment for DNase-seq03-filtering.md -- BAM filtering, deduplication, blacklist removal04-hotspot-calling.md -- Hotspot2 DHS detection and signal generation05-footprinting.md -- TF footprint detection with HINTGoal: Process raw DNase-seq FASTQ files through the ENCODE pipeline to generate DNase I hypersensitive site (DHS) peak calls.\nContext: DNase-seq identifies open chromatin via DNase I enzyme digestion. The pipeline uses BWA alignment and Hotspot2 for DHS identification.
\nencode_search_experiments(assay_title=\"DNase-seq\", biosample_term_name=\"K562\", organism=\"Homo sapiens\")\nExpected output:
\n{\n \"results\": [\n {\"accession\": \"ENCSR000DNS\", \"assay_title\": \"DNase-seq\", \"biosample_summary\": \"K562\", \"status\": \"released\"}\n ],\n \"total\": 8,\n \"limit\": 25,\n \"offset\": 0,\n \"has_more\": false,\n \"next_offset\": null\n}\nencode_list_files(experiment_accession=\"ENCSR000DNS\", file_format=\"fastq\")\nExpected output (a JSON array of file records; fields abridged):
\n[\n {\"accession\": \"ENCFF500DN1\", \"file_format\": \"fastq\", \"output_type\": \"reads\", \"file_size_human\": \"2.6 GB\", \"biological_replicates\": [1], \"status\": \"released\"},\n {\"accession\": \"ENCFF501DN2\", \"file_format\": \"fastq\", \"output_type\": \"reads\", \"file_size_human\": \"2.7 GB\", \"biological_replicates\": [1], \"status\": \"released\"}\n]\nencode_download_files(file_accessions=[\"ENCFF500DN1\", \"ENCFF501DN2\"], download_dir=\"/data/dnaseseq/fastq\")\nENCODE names every FASTQ after its accession (ENCFF500DN1.fastq.gz), with no _R1/_R2\nin the name, so the two files of a pair share no prefix and the --reads glob cannot pair\nthem. Link them into the shape the glob expects. Which mate an accession is comes from the\nENCODE file record on encodeproject.org, which carries paired_end (1 or 2) and\npaired_with; the MCP file tools do not return those two fields:
cd /data/dnaseseq/fastq\nln -s ENCFF500DN1.fastq.gz k562_rep1_R1.fastq.gz\nln -s ENCFF501DN2.fastq.gz k562_rep1_R2.fastq.gz\nnextflow run scripts/main.nf \\\n -profile local \\\n --reads '/data/dnaseseq/fastq/k562_*_R{1,2}.fastq.gz' \\\n --bwa_index /ref/bwa_index/genome.fa \\\n --chrom_sizes /ref/hg38.chrom.sizes \\\n --hotspot_center_sites /ref/hotspot2/hg38.center_sites.n100.starch \\\n --hotspot_mappable /ref/hotspot2/hg38.mappable_only.bed \\\n --rgt_data /ref/rgtdata \\\n --blacklist /ref/hg38-blacklist.v2.bed \\\n --outdir results/ \\\n -resume\nKey pipeline steps:
\n| Metric | \nThreshold | \nPurpose | \n
|---|---|---|
| SPOT score | \n> 0.4 | \nSignal portion of tags | \n
| Hotspot count | \n> 50,000 | \nSensitivity | \n
| Duplicate rate | \n< 30% | \nLibrary complexity | \n
encode_search_experiments(assay_title=\"ATAC-seq\", biosample_term_name=\"K562\", organism=\"Homo sapiens\")\nInterpretation: DNase-seq and ATAC-seq both measure accessibility but with different biases. Compare peaks from both assays -- concordant peaks are high confidence.
\nencode_get_facets(assay_title=\"DNase-seq\", organism=\"Homo sapiens\")\nExpected output (facet field names are the top-level keys):
\n{\n \"biosample_ontology.organ_slims\": [\n {\"term\": \"blood\", \"count\": 45},\n {\"term\": \"brain\", \"count\": 30}\n ]\n}\nencode_list_files(experiment_accession=\"ENCSR000DNS\", file_format=\"bed\", output_type=\"peaks\", assembly=\"GRCh38\")\nExpected output (a JSON array of file records; fields abridged):
\n[\n {\"accession\": \"ENCFF800DHS\", \"file_format\": \"bed\", \"file_type\": \"bed narrowPeak\", \"output_type\": \"peaks\", \"assembly\": \"GRCh38\", \"file_size_human\": \"1.5 MB\"}\n]\nencode_track_experiment(accession=\"ENCSR000DNS\", notes=\"K562 DNase-seq for accessibility comparison with ATAC-seq\")\nExpected output (the notes you pass are stored, not echoed back; read them with encode_list_tracked):
{\n \"tracking\": {\"accession\": \"ENCSR000DNS\", \"action\": \"tracked\"},\n \"publications_found\": 0,\n \"publications\": [],\n \"pipelines_found\": 1,\n \"pipelines\": [\n {\"title\": \"DNase-HS pipeline single-end - Version 2\", \"version\": \"2.0\", \"software\": [{\"name\": \"bwa\", \"version\": \"0.7.17\"}], \"status\": \"released\"}\n ]\n}\n| This skill produces... | \nFeed into... | \nPurpose | \n
|---|---|---|
| DHS peaks (narrowPeak) | \naccessibility-aggregation | \nUnion merge with ATAC-seq peaks | \n
| TF footprints | \nmotif-analysis | \nValidate motif predictions with footprint evidence | \n
| Signal tracks (bigWig) | \nvisualization-workflow | \nGenome browser display | \n
| Accessible regions | \nregulatory-elements | \ncCRE classification | \n
| DHS coordinates | \nvariant-annotation | \nAnnotate variants in hypersensitive sites | \n
| QC metrics | \nquality-assessment | \nValidate SPOT score and sensitivity | \n
| Pipeline parameters | \ndata-provenance | \nRecord BWA/Hotspot2 versions | \n
| DHS peak regions | \njaspar-motifs | \nScan accessible sites for known TF motifs | \n
pipeline-guide -- Parent skill with compute resource assessment and cloud setupaccessibility-aggregation -- Aggregate DHS data across samples/tissuesquality-assessment -- Evaluate pipeline output quality metricsdata-provenance -- Track all pipeline inputs, outputs, and parametersdownload-encode -- Download ENCODE DNase-seq FASTQ files for pipeline inputpublication-trust -- Verify literature claims backing analytical decisionsWhen reporting DNase-seq pipeline results:
\nfootprints/{sample}.footprints.bed and note the sequencing depth (>100M reads recommended for reliable footprints)alignment/{sample}.flagstat.txt and alignment/{sample}.dup_metrics.txt, and the insert size peak from the IS block of qc/{sample}.insert_sizes.txt. NRF (>0.8) and PBC1 (>0.9) are not produced by the workflow -- state that they were computed manually (references/03-filtering.md) or that they are unavailablemotif-analysis for TF motif enrichment in DHS peaks, or accessibility-aggregation for merging DHS data across samples