BCH709 Introduction to Bioinformatics: ChIP-Seq Tutorial

Paper Reading

Please read these papers before class:

• Kharchenko et al. (2008) Design and analysis of ChIP-seq experiments for DNA-binding proteins. Nature Biotechnology 26:1351–1359
• Landt et al. (2012) ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Research 22:1813–1831

Overview: ChIP-Sequencing

ChIP-Seq (Chromatin Immunoprecipitation followed by Sequencing) identifies genomic regions bound by a protein of interest (transcription factor, histone modification) across the entire genome.

How ChIP-Seq Works

1. Crosslink — Fix protein-DNA interactions with formaldehyde
2. Shear — Fragment chromatin by sonication or MNase digestion
3. Immunoprecipitate — Pull down target protein with a specific antibody
4. Reverse crosslinks — Release DNA from protein
5. Sequence — Sequence the enriched DNA fragments

Types of ChIP-Seq Targets

ChIP-Seq Controls

The ChIP-Seq Workflow

Raw FASTQ (ChIP + Input/Control)
    │
    ├── FastQC → MultiQC (QC)
    ├── fastp (trimming)
    │
    ├── Minimap2 (alignment to reference genome)
    │       │
    │       └── SAMtools sort + index
    │
    ├── Picard MarkDuplicates (remove PCR duplicates)
    │
    ├── deepTools bamCoverage (generate bigWig tracks)
    │
    ├── MACS3 (peak calling)
    │       │
    │       ├── Narrow peaks (TF, H3K4me3)
    │       └── Broad peaks (H3K27me3, H3K36me3)
    │
    ├── deepTools plotFingerprint (IP quality)
    ├── deepTools plotCorrelation (replicate concordance)
    ├── deepTools computeMatrix + plotHeatmap (signal profiles)
    │
    ├── HOMER / MEME-ChIP (motif analysis)
    │
    └── ChIPseeker / DiffBind (peak annotation, differential binding)

1. Environment Setup

Create Conda Environment

conda create -n chipseq -c bioconda -c conda-forge python=3.11 -y
conda activate chipseq

# Core alignment + QC tools (via conda)
conda install -c bioconda -c conda-forge fastqc fastp minimap2 samtools -y
conda install -c bioconda -c conda-forge openjdk=17 picard bedtools -y

# MACS3 and deepTools install cleanest via pip (conda has dependency conflicts)
pip install macs3 deeptools multiqc

# HOMER (optional, for motif analysis — large download)
conda install -c bioconda homer -y

# IDR (optional, for replicate reproducibility analysis)
# NOTE: the pip package named "idr" is unrelated; install from GitHub:
pip install git+https://github.com/nboley/idr.git

Note on conda errors: if you see TypeError: 'NoneType' object is not iterable from conda install, your conda solver is too old. Either upgrade conda (conda update -n base conda) or add --solver=libmamba to each install command.

Note: R packages (ChIPseeker, DiffBind) are installed separately within R (see Sections 13–14).

libcrypto fix (if samtools reports missing library):

# Check which libcrypto is present first
ls ${CONDA_PREFIX}/lib/libcrypto.so.*
# Symlink it (use the actual version you see, often .so.3 on modern installs)
ln -sf ${CONDA_PREFIX}/lib/libcrypto.so.3 ${CONDA_PREFIX}/lib/libcrypto.so.1.0.0

Verify Installations

fastp --version
minimap2 --version
samtools --version
macs3 --version
deeptools --version

Create Working Directory

mkdir -p ~/bch709/chipseq
cd ~/bch709/chipseq

2. Experimental Design Considerations

Replicates

Replicate concordance targets (ENCODE):

• IDR (Irreproducibility Discovery Rate) < 0.05
• Pearson correlation of read counts: > 0.9 (same antibody)

Sequencing Depth

Antibody Validation

Use only antibodies validated for ChIP:

• Check manufacturer’s ChIP validation data
• Run western blot to confirm specificity
• Use ENCODE validated antibodies when possible

3. Quality Control

Download Example Data

We use CTCF ChIP-Seq on human K562 cells (ENCODE experiment ENCSR000DWE) together with its matching input control (ENCSR000DWA). All three files are single-end 36 bp Illumina reads aligned to hg19.

cd ~/bch709/chipseq

# CTCF ChIP-Seq on human K562 (ENCODE ENCSR000DWE)
wget https://www.encodeproject.org/files/ENCFF001HTP/@@download/ENCFF001HTP.fastq.gz -O chip.fastq.gz

# Matching input control for K562 (ENCODE ENCSR000DWA)
wget https://www.encodeproject.org/files/ENCFF001HTT/@@download/ENCFF001HTT.fastq.gz -O input.fastq.gz

ls -lh *.fastq.gz   # expect ~910MB and ~690MB respectively

Verify the download: ENCODE’s S3 URLs are large (~1 GB each). If wget gets interrupted, the file will be truncated and later steps (fastp, minimap2) fail with unexpected end of file. Always validate:

gunzip -t chip.fastq.gz && echo "chip OK"   # should print "chip OK"
gunzip -t input.fastq.gz && echo "input OK"

If the test fails, resume with wget -c <url> (the -c flag continues from where it stopped).

For ENCODE data, browse experiments at encodeproject.org and download FASTQ files directly. Each experiment page lists its “Controls” — always pair ChIP files with their matching input, not a random control.

Run FastQC

fastqc -t 4 chip_R1.fastq.gz input.fastq.gz
multiqc .

ChIP-Seq Specific QC Metrics

4. Read Trimming

fastp performs adapter detection, quality trimming, and QC reporting in a single pass.

mkdir -p trim

# Trim ChIP sample
fastp \
  --in1 chip.fastq.gz \
  --out1 trim/chip_trimmed.fq.gz \
  --qualified_quality_phred 20 \
  --length_required 25 \
  --thread 4 \
  --html trim/chip_fastp.html \
  --json trim/chip_fastp.json

# Trim Input control
fastp \
  --in1 input.fastq.gz \
  --out1 trim/input_trimmed.fq.gz \
  --qualified_quality_phred 20 \
  --length_required 25 \
  --thread 4 \
  --html trim/input_fastp.html \
  --json trim/input_fastp.json

multiqc trim/ -n trim_report

For paired-end data, add --in2, --out2, and --detect_adapter_for_pe options.

5. Reference Genome Preparation

Download Reference

Our example FASTQs (ENCFF001HTP / ENCFF001HTT) are from early ENCODE and were aligned to hg19. For full-genome mapping we download the full assembly; for a quick tutorial run, chr19 alone (~58 MB) is enough.

# Option A — full hg19 genome (~3 GB after decompression, recommended)
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz
gunzip hg19.fa.gz
mv hg19.fa reference.fasta

# Option B — chr19 only, for a fast tutorial run (small disk / short time)
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/chr19.fa.gz
gunzip chr19.fa.gz
mv chr19.fa reference.fasta

Using Option B (single chromosome)? Only reads that map to chr19 will align (~2–6% of total reads). This is fine for testing pipeline commands but won’t give biologically meaningful peak counts. For real analysis use Option A.

Create FASTA Index

samtools faidx reference.fasta

Index the Reference for Minimap2

Minimap2 can build the index on the fly, but pre-computing it speeds up repeated alignments:

# -d : save index to file
minimap2 -d reference.mmi reference.fasta

Why Minimap2 for ChIP-Seq? Minimap2 (-ax sr mode) handles short-read DNA alignment efficiently and is suitable for both short (< 100 bp) and longer Illumina reads. Unlike HISAT2/STAR, it is not splice-aware, making it appropriate for ChIP-Seq where reads come from genomic DNA.

6. Alignment with Minimap2

Minimap2 uses the -ax sr preset for short Illumina reads.

Align ChIP Sample

set -o pipefail   # so any failure in the pipe is caught

minimap2 -ax sr -t 4 reference.mmi trim/chip_trimmed.fq.gz \
  2> chip.minimap2.log \
  | samtools sort -@ 4 -o chip.bam

samtools index chip.bam
samtools quickcheck chip.bam && echo "chip BAM OK"
samtools flagstat chip.bam

Align Input Control

minimap2 -ax sr -t 4 reference.mmi trim/input_trimmed.fq.gz \
  2> input.minimap2.log \
  | samtools sort -@ 4 -o input.bam

samtools index input.bam

Paired-end ChIP-Seq reads? Supply both FASTQ files:

minimap2 -ax sr -t 4 reference.mmi \
  trim/chip_R1_trimmed.fq.gz trim/chip_R2_trimmed.fq.gz \
  2> chip.minimap2.log | samtools sort -@ 4 -o chip.bam

Filter Low-Quality and Unmapped Reads

samtools view -b -q 30 -F 4 chip.bam > chip.filt.bam
samtools index chip.filt.bam

samtools view -b -q 30 -F 4 input.bam > input.filt.bam
samtools index input.filt.bam

7. Mark and Remove Duplicates

For ChIP-Seq, PCR duplicates must be removed (unlike variant calling where they are only marked).

picard MarkDuplicates \
  I=chip.filt.bam \
  O=chip.dedup.bam \
  M=chip.dedup.metrics \
  REMOVE_DUPLICATES=true \
  VALIDATION_STRINGENCY=SILENT

samtools index chip.dedup.bam

picard MarkDuplicates \
  I=input.filt.bam \
  O=input.dedup.bam \
  M=input.dedup.metrics \
  REMOVE_DUPLICATES=true \
  VALIDATION_STRINGENCY=SILENT

samtools index input.dedup.bam

8. Generate BigWig Signal Tracks

BigWig files are used for visualizing ChIP-Seq signal in genome browsers.

Normalize by Sequencing Depth (RPKM)

bamCoverage \
  --bam chip.dedup.bam \
  --outFileName chip.bw \
  --normalizeUsing RPKM \
  --binSize 10 \
  --numberOfProcessors 4

bamCoverage \
  --bam input.dedup.bam \
  --outFileName input.bw \
  --normalizeUsing RPKM \
  --binSize 10 \
  --numberOfProcessors 4

ChIP vs. Input Ratio (log2 fold change)

bamCompare \
  --bamfile1 chip.dedup.bam \
  --bamfile2 input.dedup.bam \
  --outFileName chip_vs_input.bw \
  --scaleFactorsMethod None \
  --normalizeUsing RPKM \
  --operation log2 \
  --binSize 10 \
  --numberOfProcessors 4

Common Error — --normalizeUsing RPKM is only valid if you also use --scaleFactorsMethod None!

Cause: bamCompare applies scaling in two places by default; combining both with RPKM is a conflict. Fix: Always pass --scaleFactorsMethod None when using --normalizeUsing RPKM (shown above).

Visualize in IGV

Load chip.bw, input.bw, and chip_vs_input.bw into IGV for visual inspection.

9. IP Quality Assessment with deepTools

Fingerprint Plot

The fingerprint plot shows how well the ChIP enrichment worked. A perfect ChIP shows a steep curve; input should be diagonal.

plotFingerprint \
  --bamfiles chip.dedup.bam input.dedup.bam \
  --labels ChIP Input \
  --plotFile fingerprint.png \
  --outRawCounts fingerprint.tab \
  --numberOfProcessors 4

Interpreting the fingerprint plot: A successful ChIP enrichment shows a steep curve (reads concentrated at a few genomic loci). The input control should be near-diagonal (reads evenly distributed). Poor enrichment = ChIP curve resembles the input.

Replicate Correlation

# Bin the genome and count reads in each bin
multiBamSummary bins \
  --bamfiles chip_rep1.dedup.bam chip_rep2.dedup.bam input.dedup.bam \
  --labels Rep1 Rep2 Input \
  --outFileName multibam.npz \
  --numberOfProcessors 4

# Plot correlation heatmap
plotCorrelation \
  --corData multibam.npz \
  --corMethod pearson \
  --skipZeros \
  --whatToPlot heatmap \
  --colorMap RdYlBu \
  --plotFile correlation.png \
  --outFileCorMatrix correlation.txt

Common Error — error: the following arguments are required: --whatToPlot/-p

Older documentation often says --plotType but current plotCorrelation uses --whatToPlot (valid values: heatmap, scatterplot).

Target: Pearson correlation between biological replicates > 0.9

10. Peak Calling with MACS3

MACS3 (Model-based Analysis of ChIP-Seq) is the standard tool for identifying enriched genomic regions.

MACS3 Key Options

Narrow Peak Calling (TF, H3K4me3, H3K9ac)

macs3 callpeak \
  -t chip.dedup.bam \
  -c input.dedup.bam \
  --format BAM \
  --gsize hs \
  --name chip_narrow \
  --outdir macs3_narrow \
  --qvalue 0.05 \
  2>&1 | tee macs3_narrow.log

Common Error — MACS3 needs at least 100 paired peaks at + and - strand to build the model, but can only find 0

Cause: MACS3 tries to build a fragment-size model from strand-paired peaks. With low-depth data (or a single chromosome / small reference), it can’t find enough peaks to model. Fix: Bypass the model and set a fixed fragment size (147 bp is a reasonable default for nucleosome-associated reads):

macs3 callpeak -t chip.dedup.bam -c input.dedup.bam --format BAM --gsize hs \
    --name chip_narrow --outdir macs3_narrow --qvalue 0.05 \
    --nomodel --extsize 147

Broad Peak Calling (H3K27me3, H3K36me3, H3K9me3)

macs3 callpeak \
  -t chip.dedup.bam \
  -c input.dedup.bam \
  --format BAM \
  --gsize hs \
  --name chip_broad \
  --outdir macs3_broad \
  --broad \
  --broad-cutoff 0.1 \
  2>&1 | tee macs3_broad.log

MACS3 Output Files

NarrowPeak Format

Count and Filter Peaks

# Count total peaks
wc -l macs3_narrow/chip_narrow_peaks.narrowPeak

# View top peaks by fold enrichment (column 7)
sort -k7,7rn macs3_narrow/chip_narrow_peaks.narrowPeak | head -20

# Filter by FDR q-value < 0.01 (column 9 = -log10 qvalue, so > 2)
awk '$9 > 2' macs3_narrow/chip_narrow_peaks.narrowPeak > chip_narrow_filtered.bed

11. Signal Profiles and Heatmaps

Visualize ChIP signal around features (e.g., TSS, peak centers).

Download Gene Annotation and Build TSS BED

# Download refGene table (UCSC format: tab-separated text)
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz
gunzip refGene.txt.gz

# Build a BED file of transcription start sites (TSS = txStart on +, txEnd on -)
awk 'BEGIN{OFS="\t"} {
    if ($4=="+") print $3, $5, $5+1, $2, "0", $4;
    else         print $3, $6-1, $6, $2, "0", $4
}' refGene.txt | sort -k1,1 -k2,2n | uniq > TSS.bed

head TSS.bed
wc -l TSS.bed

The resulting TSS.bed has one row per transcript, with coordinates pointing to the exact TSS base.

Compute Signal Matrix Around TSS

computeMatrix reference-point \
  --scoreFileName chip.bw \
  --regionsFileName TSS.bed \
  --referencePoint TSS \
  --upstream 3000 \
  --downstream 3000 \
  --binSize 50 \
  --numberOfProcessors 4 \
  -o tss_matrix.gz

plotHeatmap \
  --matrixFile tss_matrix.gz \
  --outFileName tss_heatmap.png \
  --colorMap Blues \
  --whatToShow 'heatmap and colorbar' \
  --zMin 0 --zMax 5

Profile Plot

plotProfile \
  --matrixFile tss_matrix.gz \
  --outFileName tss_profile.png \
  --plotTitle "ChIP Signal at TSS"

12. Motif Analysis

Motif analysis identifies DNA sequence motifs enriched in ChIP-Seq peaks — typically the binding motif of the immunoprecipitated protein.

HOMER Motif Analysis

# Install HOMER genome (run once)
configureHomer.pl -install hg38

findMotifsGenome.pl usage: findMotifsGenome.pl <peaks.bed> <genome> <output_dir/> [options]

# Prepare peak file (convert narrowPeak to BED6)
awk '{print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6}' \
  chip_narrow_filtered.bed > peaks_homer.bed

# Find motifs (de novo + known)
findMotifsGenome.pl peaks_homer.bed hg38 motif_output/ -size 200 -mask -p 4

HOMER Output

13. Peak Annotation with ChIPseeker (R)

ChIPseeker annotates peaks relative to genomic features (promoter, exon, intron, etc.).

Install ChIPseeker

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("ChIPseeker", "TxDb.Hsapiens.UCSC.hg38.knownGene", "clusterProfiler"))

Annotate Peaks

library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(clusterProfiler)

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

# Read peaks
peaks <- readPeakFile("chip_narrow_filtered.bed")

# Annotate
peakAnno <- annotatePeak(peaks, tssRegion=c(-3000, 3000),
                          TxDb=txdb, annoDb="org.Hs.eg.db")

# Plot annotation pie chart
plotAnnoPie(peakAnno)

# Plot distance to TSS
plotDistToTSS(peakAnno, title="Distribution of ChIP peaks\nrelative to TSS")

# Gene ontology of nearby genes
genes <- seq2gene(peaks, tssRegion=c(-1000, 1000),
                   flankDistance=3000, TxDb=txdb)
pathway <- enrichPathway(genes)
dotplot(pathway)

14. Differential Binding Analysis with DiffBind (R)

Compare ChIP-Seq signal between two conditions (e.g., treated vs. untreated).

Install DiffBind

BiocManager::install("DiffBind")

Sample Sheet Format

The sample sheet CSV must contain these columns:

Run DiffBind

library(DiffBind)

# Load sample sheet
samples <- read.csv("samplesheet.csv")
dba_obj <- dba(sampleSheet=samples)

# Count reads in peaks
dba_obj <- dba.count(dba_obj, bUseSummarizeOverlaps=TRUE)

# Set contrast
dba_obj <- dba.contrast(dba_obj, categories=DBA_CONDITION)

# Run differential analysis (DESeq2 by default)
dba_obj <- dba.analyze(dba_obj)

# Get significant peaks
db_peaks <- dba.report(dba_obj, th=0.05)

15. IDR Analysis (Irreproducibility Discovery Rate)

IDR measures reproducibility of peak calls across replicates. ENCODE requires IDR < 0.05.

# Sort peaks by score (column 5)
sort -k5,5rn chip_rep1_peaks.narrowPeak > rep1_sorted.bed
sort -k5,5rn chip_rep2_peaks.narrowPeak > rep2_sorted.bed

# Run IDR
idr \
  --samples rep1_sorted.bed rep2_sorted.bed \
  --input-file-type narrowPeak \
  --rank p.value \
  --output-file idr_peaks.txt \
  --plot \
  --log-output-file idr.log

Common Error — AttributeError: module 'numpy' has no attribute 'int'

IDR 2.0.3 is not compatible with NumPy ≥ 1.20 (where np.int was removed). Fix: Downgrade NumPy in the chipseq env, or use the bioconda build (which patches this):

pip install "numpy<1.20"
# OR use the bioconda-packaged idr:
conda install -c bioconda idr

Common Error — idr: Image Data Repository access library (wrong package)

pip install idr pulls a completely unrelated package (for microscopy images). Always install from the Kundaje/Boley repo:

pip uninstall idr -y
pip install git+https://github.com/nboley/idr.git

16. Full Workflow Summary

17. Troubleshooting — Quick Reference

18. Cleanup

# Remove large intermediate files when done
rm -f chip.bam chip.bam.bai chip.filt.bam chip.filt.bam.bai
rm -f input.bam input.bam.bai input.filt.bam input.filt.bam.bai
rm -rf trim/

conda deactivate
# Optional: remove environment entirely
conda env remove --name chipseq

References