End-to-end ChIP-seq workflow from FASTQ files to annotated peaks. Covers QC, alignment, peak calling with MACS3, and peak annotation with ChIPseeker. Use when processing ChIP-seq data from alignment through peak annotation.
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:
packageVersion('<pkg>') then ?function_name to verify parameters<tool> --version then <tool> --help to confirm flagsIf 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.
"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.
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.
A ChIP-seq peakset is decided at four seams, not inside the caller.
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.--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.
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
| Commitment | Choice | Consequence inherited downstream |
|------------|--------|----------------------------------|
| Build + blacklist + effective genome size | One genome build; the matching ENCODE blacklist BED; -g hs/mm/numeric | Mixed builds mis-place peaks; skipping the blacklist plants reproducible false peaks; wrong -g mis-scales p-values |
| Control pairing | Each IP has its input/IgG | No control => peaks at open chromatin / CN-amplified loci |
| Peak shape | Narrow (TF, H3K4me3, H3K27ac) vs broad (H3K27me3, H3K36me3, H3K9me3) | Broad marks called with narrow settings fragment into many small peaks |
| Fragment model | PE: -f BAMPE (real fragments); SE: -f BAM + --nomodel --extsize from predictd/xcorr | BAMPE silently ignores --shift/--extsize |
samtools view -q 30), coordinate-sort.--broad).--scaleFactor + --normalizeUsing None for spike-in (order-trap: RPKM erases the spike-in global shift).Pipeline-level selection only; mechanism lives in the component skills.
| Fork | Lean toward | Hand off to |
|------|-------------|-------------|
| Caller | MACS3 (standard IP+input); SEACR (CUT&RUN/CUT&Tag, low background); Genrich (some ChIP/ATAC, built-in blacklist/replicate handling) | chip-seq/peak-calling, chip-seq/cut-and-run-tag |
| Narrow vs broad | Narrow: TFs, H3K4me3, H3K27ac. Broad (--broad --broad-cutoff 0.1): H3K27me3, H3K36me3, H3K9me3 | chip-seq/peak-calling |
| Reproducibility | ENCODE IDR (per-rep + pooled pseudo-reps) for TFs; naive overlap acceptable for exploratory histone | chip-seq/peak-calling |
| Consensus set | Pool for a union/consensus set AFTER IDR selects the reproducible threshold | chip-seq/differential-binding |
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).
# Standard experiment: RPKM/CPM is fine. SPIKE-IN experiment: this would ERASE the global shift.
bamCoverage -b aligned/IP_rep1.final.bam -o bigwig/IP_rep1.bw --normalizeUsing RPKM -p 8
# Spike-in (ChIP-Rx): bamCoverage --scaleFactor <spike-in factor> --normalizeUsing None (chip-seq/spike-in-normalization)
Annotation uses a project GTF via makeTxDbFromGFF() when provided, else a pre-built TxDb. overlap='all' couples gene assignment with feature overlap (host-gene convention); default overlap='TSS' assigns the nearest-TSS gene independently. Full code: examples/peak_annotation.R.
| After | Gate | Interpretation | |-------|------|----------------| | QC/trim | Q30 >85%, adapter <5% | DNA higher quality than RNA | | Alignment | Mapping >80%, unique >70% | Low unique = repeats/contamination/wrong build | | PRE-dedup | NRF >0.8, PBC1 >0.8 | Low complexity = over-amplification/low input; MUST be computed before dedup | | Peaks | FRiP >1% (TF) / >5% (sharp histone; broad marks run lower); NSC >1.05; RSC >0.8; fingerprint separates IP/input | Low FRiP/flat fingerprint = weak antibody or failed enrichment (chip-seq/chipseq-qc) | | IDR | rescue ratio and self-consistency ratio both <=2 | Poor replicate consistency; run IDR on PER-replicate peaks |
| Symptom | Cause | Fix |
|---------|-------|-----|
| Reproducible peaks over satellite/rDNA/high-signal regions | ENCODE blacklist never subtracted | bedtools intersect -v the blacklist BED before calling (committed step) |
| NRF/PBC ~1.0 and uninformative | Computed after markdup -r | Compute complexity on the PRE-dedup, filtered BAM |
| Peaks at open chromatin / CN-amplified loci | Called without a matched control | Pair each IP with its input/IgG in callpeak -c |
| H3K27me3/H3K9me3 fragmented into many tiny peaks | Broad mark called with narrow settings | Add --broad --broad-cutoff 0.1 |
| --shift/--extsize had no effect | Used with -f BAMPE (ignored for PE) | Use -f BAM + --nomodel for SE; BAMPE derives fragments |
| Spike-in global shift disappears in tracks | bamCoverage RPKM/CPM re-equalized depth | --scaleFactor + --normalizeUsing None (chip-seq/spike-in-normalization) |
| "IDR" numbers look too good | IDR run on a pooled peakset | Run IDR on per-replicate peaks + pooled pseudo-replicates |
The complete runnable scripts are in this skill's examples/ (narrow_peak_workflow.sh, peak_annotation.R).
npx skills add GPTomics/bio-workflows-chipseq-pipeline下载完整 Skill 目录，包含 SKILL.md 及所有相关文件
Search for places (restaurants, cafes, etc.) via Google Places API proxy on localhost.
Interact with GitHub using the `gh` CLI. Use `gh issue`, `gh pr`, `gh run`, and `gh api` for issues, PRs, CI runs, and advanced queries.
Create or update AgentSkills. Use when designing, structuring, or packaging skills with scripts, references, and assets.
Start voice calls via the OpenClaw voice-call plugin.
Notion API for creating and managing pages, databases, and blocks.
Gemini CLI for one-shot Q&A, summaries, and generation.
Category:developer