# Chipseq Epigenomics

> Align/call ChIP-seq peaks with Bowtie2+MACS3, run DiffBind binding and ChIPseeker annotation in R, visualize with deepTools. Use for TF/histone peak calling, narrow vs broad peaks, FRiP/NSC/RSC QC, or peak annotation.

- Skill: `pavel-kravchenko/chipseq-epigenomics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/chipseq-epigenomics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/chipseq-epigenomics/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/chipseq-epigenomics

---


# ChIP-seq Epigenomics

## When to Use

- Processing raw ChIP-seq FASTQ (TF binding or histone modification) into aligned, deduplicated BAM files
- Calling narrow peaks (TF, H3K4me3, H3K27ac) or broad peaks (H3K27me3, H3K9me3) with MACS3
- Running differential binding analysis between conditions (DiffBind) or annotating peaks to genomic features (ChIPseeker)
- Computing ENCODE-style QC metrics: FRiP, NSC/RSC, duplication rate, IDR reproducibility
- Visualizing enrichment as BigWig tracks, TSS heatmaps, or average profile plots (deepTools)

## Version Compatibility

- MACS3 ≥ 3.0, Bowtie2 ≥ 2.5, samtools ≥ 1.19, Picard ≥ 3.1, deepTools ≥ 3.5
- R ≥ 4.3, Bioconductor ≥ 3.18: DiffBind ≥ 3.12, ChIPseeker ≥ 1.38
- Python ≥ 3.10: pandas ≥ 2.0, pybedtools ≥ 0.10, pyBigWig ≥ 0.3

## Prerequisites

- conda/mamba environment with: `trim-galore`, `bowtie2`, `samtools`, `picard`, `macs3`, `deeptools`, `bedtools`, `pybedtools`, `pyBigWig`
- R/Bioconductor packages: `DiffBind`, `ChIPseeker`, `TxDb.Hsapiens.UCSC.hg38.knownGene`, `org.Hs.eg.db`
- ENCODE blacklist BED for the genome build (e.g. `hg38-blacklist.v2.bed`)
- Familiarity with BAM/SAM handling (`bio-applied-variant-calling-and-snp-analysis`) and BED intervals

**Goal:** Turn paired-end ChIP-seq FASTQ into a filtered, deduplicated BAM and call peaks with MACS3.
**Approach:** trim adapters, align with Bowtie2, filter to properly-paired high-MAPQ reads, remove ENCODE blacklist regions, mark/remove duplicates with Picard, then call peaks with MACS3 (`--broad` for repressive/heterochromatin marks, default narrow mode otherwise).

```bash
# 1. Trim adapters (Trim Galore auto-detects adapter type)
trim_galore --paired --fastqc --cores 4 --quality 20 --length 20 \
    --output_dir trimmed/ sample_R1.fastq.gz sample_R2.fastq.gz

# 2. Align, keep only properly-paired, high-confidence, primary alignments
bowtie2 -x /data/indices/hg38 \
    -1 trimmed/sample_R1_val_1.fq.gz -2 trimmed/sample_R2_val_2.fq.gz \
    -p 8 --no-mixed --no-discordant 2> logs/bowtie2.log | \
    samtools view -bS -F 4 -F 256 -q 30 -f 2 | \
    samtools sort -@ 4 -o aligned/sample.bam
samtools index aligned/sample.bam

# 3. Remove ENCODE blacklist regions
bedtools intersect -a aligned/sample.bam -b hg38-blacklist.v2.bed -v \
    > aligned/sample_filtered.bam
samtools index aligned/sample_filtered.bam

# 4. Remove PCR/optical duplicates
picard MarkDuplicates I=aligned/sample_filtered.bam O=dedup/sample_dedup.bam \
    M=dedup/sample_dup_metrics.txt REMOVE_DUPLICATES=true \
    VALIDATION_STRINGENCY=SILENT 2> logs/picard_dedup.log
samtools index dedup/sample_dedup.bam

# 5a. Narrow peaks (TF, H3K4me3, H3K27ac) — needs matched input/IgG control
macs3 callpeak -t dedup/sample_dedup.bam -c dedup/input_dedup.bam \
    -f BAMPE -g hs -n sample --outdir peaks/ -q 0.05 --keep-dup all

# 5b. Broad peaks (H3K27me3, H3K9me3)
macs3 callpeak -t dedup/h3k27me3_dedup.bam -c dedup/input_dedup.bam \
    -f BAMPE -g hs -n h3k27me3 --outdir peaks/ \
    --broad --broad-cutoff 0.1 --keep-dup all
```

**Goal:** Confirm peak calls pass ENCODE QC thresholds before downstream analysis.
**Approach:** compute FRiP (fraction of reads in peaks) from the deduplicated BAM and the called peak set; flag against the standard TF (≥5%) vs. histone (≥1%) thresholds.

```python
import subprocess

def frip_score(bam: str, peaks_bed: str) -> float:
    """Fraction of Reads in Peaks (FRiP) — quick estimate via bedtools.

    For production use, prefer featureCounts -a peaks.narrowPeak -F SAF,
    which correctly counts fragments rather than raw BAM lines.
    """
    total = int(subprocess.check_output(
        ["samtools", "view", "-c", "-F", "4", bam]).decode().strip())
    in_peaks = int(subprocess.check_output(
        ["bedtools", "intersect", "-a", bam, "-b", peaks_bed, "-u"]
    ).decode().count("\n"))
    return in_peaks / total

frip = frip_score("sample_dedup.bam", "sample_peaks.narrowPeak")
threshold = 0.05  # use 0.01 for histone marks
print(f"FRiP = {frip:.3f} ({'PASS' if frip >= threshold else 'FAIL'}, threshold={threshold})")
```

**Goal:** Parse and filter MACS3 narrowPeak output, and compare peak sets across conditions.
**Approach:** read narrowPeak as a DataFrame (BED6+4), filter by q-value, then use pybedtools for reciprocal-overlap comparisons between replicates or conditions.

```python
import pandas as pd
import pybedtools

def read_narrowpeak(path: str) -> pd.DataFrame:
    """Load a MACS3 narrowPeak file into a DataFrame.

    Column 9 (qvalue) is -log10(q); qvalue >= 1.3 corresponds to q <= 0.05.
    """
    cols = ["chrom", "start", "end", "name", "score", "strand",
            "signalValue", "pvalue", "qvalue", "peak"]
    return pd.read_csv(path, sep="\t", header=None, names=cols)

peaks = read_narrowpeak("sample_peaks.narrowPeak")
sig_peaks = peaks[peaks["qvalue"] >= 1.3]  # q <= 0.05
print(f"{len(sig_peaks)} of {len(peaks)} peaks pass q <= 0.05")

# Compare two conditions with reciprocal 50% overlap
a = pybedtools.BedTool("condition_A.narrowPeak")
b = pybedtools.BedTool("condition_B.narrowPeak")
shared = a.intersect(b, f=0.5, r=True)
a_only = a.intersect(b, f=0.5, r=True, v=True)
b_only = b.intersect(a, f=0.5, r=True, v=True)
print(f"Shared: {len(shared)}, A-only: {len(a_only)}, B-only: {len(b_only)}")
```

**Goal:** Generate normalized BigWig tracks and a TSS enrichment heatmap.
**Approach:** normalize BAM to RPKM-scaled BigWig with deepTools, then compute and plot a signal matrix around TSS.

```bash
bamCoverage -b dedup/sample_dedup.bam -o bigwig/sample.bw \
    --normalizeUsing RPKM --binSize 10 --numberOfProcessors 8 --extendReads

computeMatrix reference-point -S bigwig/sample.bw bigwig/input.bw \
    -R genes_hg38.bed --referencePoint TSS -b 3000 -a 3000 --binSize 10 \
    -o matrix/tss_matrix.gz --outFileSortedRegions matrix/tss_regions_sorted.bed

plotHeatmap -m matrix/tss_matrix.gz -out figures/heatmap_tss.png \
    --colorMap Blues --samplesLabel 'ChIP' 'Input' --regionsLabel 'Genes'
plotProfile -m matrix/tss_matrix.gz -out figures/profile_tss.png \
    --samplesLabel 'ChIP' 'Input' --perGroup
```

**Goal:** Test for statistically significant differential binding between two conditions.
**Approach:** build a `DiffBind` sample sheet (SampleID, Condition, bamReads, bamControl, Peaks), count reads over the consensus peak set, normalize (RLE, same as DESeq2), and run the built-in DESeq2 contrast.

```r
library(DiffBind)

# samplesheet.csv columns: SampleID, Condition, Replicate,
#   bamReads, bamControl, Peaks, PeakCaller
dba_obj <- dba(sampleSheet = "samplesheet.csv")
dba_obj <- dba.count(dba_obj, bUseSummarizeOverlaps = TRUE)
dba_obj <- dba.normalize(dba_obj, normalize = DBA_NORM_RLE)
dba_obj <- dba.contrast(dba_obj, categories = DBA_CONDITION, minMembers = 2)
dba_obj <- dba.analyze(dba_obj, method = DBA_DESEQ2)

# FDR < 0.05, |log2FC| > log2(1.5)
db_peaks <- dba.report(dba_obj, th = 0.05, fold = log2(1.5))
dba.plotVolcano(dba_obj)
dba.plotPCA(dba_obj, DBA_CONDITION, label = DBA_ID)
rtracklayer::export(db_peaks, "diffbind_results.bed")
```

**Goal:** Annotate peaks to genomic features (promoter, exon, intron, intergenic) and distance-to-TSS.
**Approach:** load a `TxDb` for the genome build and run `ChIPseeker::annotatePeak`.

```r
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
peaks <- readPeakFile("peaks/sample_peaks.narrowPeak")
anno <- annotatePeak(peaks, tssRegion = c(-2000, 200),
                      TxDb = txdb, annoDb = "org.Hs.eg.db")
plotAnnoPie(anno)
plotDistToTSS(anno)

# Promoter-proximal peaks (within 2 kb of TSS)
anno_df <- as.data.frame(anno)
promo_peaks <- anno_df[abs(anno_df$distanceToTSS) < 2000, ]
```

## Quality Metrics

| Metric | Recommendation | Source |
|--------|---------------|--------|
| FRiP (TF) | ≥ 5% | reads in peaks / mapped reads |
| FRiP (histone) | ≥ 1% | reads in peaks / mapped reads |
| NSC (normalized strand coefficient) | > 1.05 | deepTools/SPP |
| RSC (relative strand correlation) | > 0.8 | SPP |
| Duplication rate | < 30% (TF) | Picard `MarkDuplicates` metrics |
| IDR (replicate reproducibility) | < 0.05 | ENCODE IDR pipeline |

## Pitfalls

- **Input control is mandatory** — always subtract input (or IgG) for peak calling; without it MACS3 cannot distinguish signal from open-chromatin background
- **Paired-end vs single-end** — use `-f BAMPE` for PE, `-f BAM` for SE in MACS3; mismatching this silently halves effective fragment length
- **Broad vs narrow peaks** — H3K4me3/H3K27ac are narrow; H3K27me3/H3K9me3 need `--broad --broad-cutoff`, or peaks fragment into many small calls
- **Blacklist regions** — always filter ENCODE blacklisted regions before dedup/peak calling; they generate spurious high-signal artifacts
- **Cross-condition normalization** — RPKM/CPM normalization is invalid for comparing histone mark abundance across conditions with global changes; use spike-in (Orlando method) normalization instead
- **Minimum depth** — TF ChIP-seq needs ≥20M uniquely mapped reads; histone marks need 40-80M due to broader genomic coverage

## See Also

- `bio-applied-tf-footprinting` — TF footprinting from ATAC-seq/DNase-seq
- `atac-seq-analysis` — open chromatin peak calling, ATAC-specific QC
- `rnaseq` — differential expression, DESeq2 normalization
- `python-bio-data-visualization` — heatmaps, genome browser tracks, volcano plots