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).

# 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.

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.

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.

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.

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.

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

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