/deeptools-ngs-analysis

name: deeptools-ngs-analysis
description: "NGS CLI for ChIP/RNA/ATAC-seq. BAM→bigWig with RPGC/CPM/RPKM, sample correlation/PCA, heatmaps/profiles around features, fingerprints. For alignment use STAR/BWA; for peak calling use MACS2."
license: BSD-3-Clause

deepTools — NGS Data Analysis Toolkit

Overview

deepTools is a command-line toolkit for processing and visualizing high-throughput sequencing data. It converts BAM alignments to normalized coverage tracks (bigWig), performs quality control (correlation, PCA, fingerprint), and generates publication-quality heatmaps and profile plots around genomic features. Supports ChIP-seq, RNA-seq, ATAC-seq, and MNase-seq.

When to Use

• Converting BAM files to normalized bigWig coverage tracks
• Comparing ChIP-seq treatment vs input control (log2 ratio tracks)
• Assessing sample quality: replicate correlation, PCA, coverage depth
• Evaluating ChIP enrichment strength (fingerprint plots)
• Creating heatmaps and profile plots around TSS, peaks, or other genomic regions
• Analyzing ATAC-seq data with Tn5 offset correction
• Generating strand-specific RNA-seq coverage tracks
• Use **omics-plotting** SKILL for custom figures from exported score matrices; standard signal heatmaps/profiles use deeptools' own `plotHeatmap`/`plotProfile`
• For **read alignment**, use STAR, BWA, or bowtie2 instead
• For **peak calling**, use MACS2 or HOMER instead
• For **BAM/VCF file manipulation**, use pysam instead

Prerequisites

pip install deeptools
# Verify installation
bamCoverage --version

**Input requirements**: BAM files must be sorted and indexed (`.bai` file present). Generate index with `samtools index input.bam`. BED files for genomic regions (genes, peaks) in standard 3+ column format.

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: normalizationMethod
    kind: required
    source: user
    ask: "How should coverage be scaled so tracks from different samples can be compared - by read depth, by genome coverage, or not at all?"
    default: "none - raw coverage, comparable only within one sample"

  - id: D2
    param: effectiveGenomeSize
    kind: derived
    source: upstream
    depends_on: [D1]
    ask: "Which assembly's mappable size should the coverage normalization use?"
    default: "looked up for the alignment reference"
    skip_if: "normalization does not require a genome size"

  - id: D3
    param: readExtension
    kind: required
    source: user
    ask: "Should reads be extended to the sequenced fragment length? Correct for ChIP and ATAC, wrong for spliced RNA."
    default: "not extended"

  - id: D4
    param: duplicateHandling
    kind: required
    source: user
    ask: "Should reads flagged as duplicates be ignored when building coverage?"
    default: "included"

  - id: D5
    param: mappingQualityFloor
    kind: required
    source: user
    ask: "Below what mapping confidence should reads be excluded from the signal?"
    default: "no filter"

  - id: D6
    param: binSize
    kind: optional
    source: user
    ask: "At what resolution should signal be summarized?"
    default: "50 bp"

  - id: D7
    param: comparisonOperation
    kind: optional
    source: user
    ask: "When comparing two samples, should the track hold a log ratio, a plain ratio, or a difference?"
    default: "log2 ratio"

  - id: D8
    param: profileAnchor
    kind: optional
    source: user
    ask: "Should profiles be anchored on feature starts, ends, centres, or scaled across whole features?"
    default: "feature start, 500 bp either side"

  - id: D9
    param: processors
    kind: never_ask
    source: data
    reason: "Affects runtime only, not the signal"
    default: "min(8, available_cores)"

D1 and D3 are the pair that makes tracks comparable or quietly misleading. Unnormalized tracks put the deeper-sequenced sample higher everywhere, which reads as biology on a genome browser; extending RNA-seq reads to a fragment length fills introns with signal that was never sequenced.

Quick Start

# Convert BAM to normalized bigWig
bamCoverage --bam sample.bam --outFileName sample.bw \
    --normalizeUsing RPGC --effectiveGenomeSize 2913022398 \
    --binSize 10 --numberOfProcessors 8

# Create heatmap around TSS
computeMatrix reference-point -S sample.bw -R genes.bed \
    -b 3000 -a 3000 --referencePoint TSS -o matrix.gz
plotHeatmap -m matrix.gz -o heatmap.png --colorMap RdBu

Core API

1. BAM to Coverage Conversion

Convert BAM alignments to normalized coverage tracks (bigWig or bedGraph).

# Basic conversion with RPGC normalization
bamCoverage --bam input.bam --outFileName output.bw \
    --normalizeUsing RPGC --effectiveGenomeSize 2913022398 \
    --binSize 10 --numberOfProcessors 8 \
    --extendReads 200 --ignoreDuplicates

# CPM normalization (simpler, no genome size needed)
bamCoverage --bam input.bam --outFileName output.bw \
    --normalizeUsing CPM --binSize 10 -p 8

# RNA-seq: strand-specific coverage
bamCoverage --bam rnaseq.bam --outFileName forward.bw \
    --filterRNAstrand forward --normalizeUsing CPM -p 8
# IMPORTANT: Never use --extendReads for RNA-seq (spans splice junctions)

2. Sample Comparison

Compare treatment vs control or generate ratio tracks.

# Log2 ratio: treatment / control
bamCompare -b1 treatment.bam -b2 control.bam -o log2ratio.bw \
    --operation log2 --scaleFactorsMethod readCount \
    --extendReads 200 -p 8

# Subtract control from treatment
bamCompare -b1 treatment.bam -b2 control.bam -o subtract.bw \
    --operation subtract --scaleFactorsMethod readCount

3. Quality Control

Assess sample quality, replicate concordance, and enrichment strength.

# Sample correlation heatmap
multiBamSummary bins --bamfiles rep1.bam rep2.bam rep3.bam \
    -o counts.npz --binSize 10000 -p 8
plotCorrelation -in counts.npz --corMethod pearson \
    --whatToShow