Bio Workflows Chipseq Pipeline

Version Compatibility

Reference examples tested with: Bowtie2 2.5.3+, MACS3 3.0+, HOMER 4.11+, bedtools 2.31+, deepTools 3.5+, fastp 0.23+, samtools 1.19+, ChIPseeker 1.38+

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Note: macs3 callpeak -f BAMPE uses real fragment lengths and IGNORES --shift/--extsize/--nomodel (those apply to single-end -f BAM); the -g shortcut (hs/mm) sets the effective genome size and must match the build/read length. Confirm in-tool before quoting.

ChIP-seq Pipeline

"Process my ChIP-seq data from FASTQ to annotated peaks" -> Chain QC/trim, alignment, pre-dedup complexity QC, dedup + blacklist filtering, control-matched peak calling, reproducibility, signal tracks, and annotation.

• CLI + R: fastp -> bowtie2 -> (NRF/PBC pre-dedup) -> samtools markdup -> chrM/blacklist filter -> macs3 callpeak (IP vs input) -> IDR -> bamCoverage -> ChIPseeker

This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step. Every step below cross-references the component skill that teaches its mechanism.

The governing principle

A ChIP-seq peakset is decided at four seams, not inside the caller.

1. The reference build + blacklist version + effective genome size is one coordinate commitment made once and inherited by everything downstream — peak coordinates, signal-track scaling, and every overlap. Blacklist filtering is a committed pipeline step, not optional cleanup: ENCODE-blacklisted regions (satellite/rDNA/high-signal artifacts) produce reproducible false peaks in every dataset regardless of biology, so they are removed before calling (Amemiya 2019).
2. An IP is only interpretable against its matched control — the control IS the enrichment background. Calling peaks without the right input/IgG fabricates peaks at open/accessible and copy-number-amplified regions. Pair each IP with its control at the calling step.
3. Library-complexity QC (NRF/PBC1/PBC2) is computed on the PRE-dedup BAM. After markdup -r the duplicates are gone, so computing complexity afterward reads ~1.0 and is meaningless. Compute it on the filtered, position-sorted BAM before removing duplicates.
4. Normalization must not silently undo the experiment. deepTools RPKM/CPM rescales every library to the same depth, which ERASES a spike-in global-shift signal (the whole point of ChIP-Rx). For spike-in experiments use --scaleFactor + --normalizeUsing None; for standard experiments RPKM/CPM is fine (chip-seq/spike-in-normalization).

Reproducibility corollary: pool replicates for a consensus peakset, but keep PER-REPLICATE peaks — IDR needs individual replicates plus pooled pseudo-replicates; running IDR on an already-pooled peakset is not IDR.

Pipeline map

FASTQ (IP + matched Input, replicates)
  | [1] QC & trim -----------------> fastp                (read-qc/fastp-workflow)
  v
  | [2] Align ---------------------> bowtie2 (-q30 unique) (read-alignment/bowtie2-alignment)
  v     ^-- commitment: build + blacklist version + effective genome size
  | [3] Complexity QC (PRE-dedup) -> NRF/PBC1/PBC2         (chip-seq/chipseq-qc)
  v
  | [4] Dedup + filter ------------> markdup -r; drop chrM; SUBTRACT ENCODE blacklist  (alignment-files/duplicate-handling)
  v
  | [5] Peak calling (IP vs input)-> macs3 callpeak (narrow | --broad)  (chip-seq/peak-calling)
  v     ^-- keep PER-REPLICATE peaks for IDR
  | [6] Reproducibility -----------> IDR (per-rep + pooled pseudo-reps)  (chip-seq/peak-calling)
  v
  | [7] Signal tracks -------------> bamCoverage (RPKM | spike-in scaleFactor)  (chip-seq/chipseq-visualization)
  v
  | [8] QC + Annotate -------------> FRiP/NSC/RSC/fingerprint; ChIPseeker  (chip-seq/chipseq-qc, peak-annotation)
  v
Blacklist-filtered, annotated, reproducible peaks

Made-once commitments

The canonical order and why

1. QC/trim (fastp) both IP and input.
2. Align (bowtie2), keep uniquely-mapped (samtools view -q 30), coordinate-sort.
3. Compute NRF/PBC1/PBC2 on the PRE-dedup BAM — order-trap: after dedup they are meaningless.
4. Mark/remove duplicates (collate -> fixmate -m -> sort -> markdup -r), then drop chrM and subtract the ENCODE blacklist — order-trap: skipping the blacklist leaves reproducible artifact peaks.
5. Call peaks against the matched control (narrow or --broad).
6. IDR on per-replicate peaks (+ pooled pseudo-replicates) — order-trap: IDR on a pooled peakset is not IDR.
7. Signal tracks — RPKM/CPM for standard; --scaleFactor + --normalizeUsing None for spike-in (order-trap: RPKM erases the spike-in global shift).
8. QC (FRiP/NSC/RSC/fingerprint) and annotate (ChIPseeker).

Choosing the caller and peak shape

Pipeline-level selection only; mechanism lives in the component skills.

Primary path: Bowtie2 + MACS3 + ChIPseeker

Goal: turn IP+input FASTQ into a blacklist-filtered, control-matched, annotated peakset.

Approach: align and keep unique reads, measure complexity before dedup, dedup + drop chrM + subtract the blacklist, call against the control, then annotate. Full runnable script: examples/narrow_peak_workflow.sh; annotation: examples/peak_annotation.R.

bowtie2 -p 8 -x bt2_index/genome -1 trimmed/${s}_R1.fq.gz -2 trimmed/${s}_R2.fq.gz \
    --no-mixed --no-discordant --maxins 1000 2> aligned/${s}.log \
  | samtools view -@4 -bS -q 30 - | samtools sort -@4 -o aligned/${s}.sorted.bam
samtools index aligned/${s}.sorted.bam

# Complexity QC on the PRE-dedup BAM (NRF = distinct positions / total; PBC1 = singletons / distinct).
# Counted per-mate here (close to ENCODE fragment-level values); use `bamtobed -bedpe` for exact parity.
bedtools bamtobed -i aligned/${s}.sorted.bam | awk 'BEGIN{OFS="\t"}{print $1,$2,$3,$6}' | sort | uniq -c \
  | awk '{tot+=$1; dist++; if($1==1) one++} END{printf "NRF=%.3f PBC1=%.3f\n", dist/tot, one/dist}'

# Dedup, drop chrM, then SUBTRACT the ENCODE blacklist (committed step, not optional)
samtools collate -@8 -O -u aligned/${s}.sorted.bam | samtools fixmate -m -u - - \
  | samtools sort -@8 -u - | samtools markdup -r -@8 - aligned/${s}.dedup.bam
samtools index aligned/${s}.dedup.bam
samtools idxstats aligned/${s}.dedup.bam | cut -f1 | grep -v -e '^chrM$' -e '^MT$' \
  | xargs samtools view -b aligned/${s}.dedup.bam > aligned/${s}.nochrM.bam
bedtools intersect -v -a aligned/${s}.nochrM.bam -b ENCODE_blacklist.bed > aligned/${s}.final.bam
samtools index aligned/${s}.final.bam

# Narrow (TFs, sharp marks) vs broad (spreading marks). -f BAMPE uses real fragment sizes.
macs3 callpeak -t aligned/IP_rep1.final.bam aligned/IP_rep2.final.bam \
    -c aligned/Input_rep1.final.bam aligned/Input_rep2.final.bam \
    -f BAMPE -g hs -n experiment --outdir peaks -q 0.01 --keep-dup all   # dedup done upstream (markdup -r); tell MACS3 to keep all
# Broad marks: add  --broad --broad-cutoff 0.1  (do NOT call H3K27me3 with narrow settings)

For IDR, call peaks PER REPLICATE (and on pooled pseudo-replicates) with a relaxed -q, then run idr across them (chip-seq/peak-calling). For higher confidence, intersect a second caller (HOMER -style histone for all histone marks).