🧠 ChIP-Seq Fundamentals

Q: What is ChIP-Seq?

ChIP-Seq (Chromatin Immunoprecipitation Sequencing) combines chromatin immunoprecipitation with next-generation sequencing to identify genome-wide binding sites of proteins (transcription factors, histones) and their modifications. It reveals regulatory elements, epigenetic marks, and protein-DNA interactions.

Q: What is the difference between transcription factor and histone modification ChIP-Seq?

Transcription factor (TF) ChIP-Seq produces sharp, narrow peaks (100-500 bp) at specific binding sites. Histone modification ChIP-Seq produces broader, often bimodal peaks: active marks like H3K4me3 are sharp, while repressive marks like H3K27me3 and H3K9me3 are broad (several kb). Different peak callers are needed.

Q: Which peak caller should I use?

Use MACS2 or MACS3 for transcription factors and sharp histone marks (H3K4me3, H3K27ac). Use SICER or epic for broad histone marks (H3K27me3, H3K9me3, H3K36me3). For differential binding, use DiffBind or csaw. For ATAC-Seq, MACS2 works well.

Q: What is a control (input) sample in ChIP-Seq?

An input control is chromatin that has not been immunoprecipitated. It's used to correct for background, copy number, and open chromatin biases. MACS2 uses the input to model local background and calculate fold enrichment. Without input, peak calling is less accurate.

Q: What is the difference between ChIP-Seq and ATAC-Seq?

ChIP-Seq uses antibodies to pull down specific proteins or histone modifications, giving targeted binding maps. ATAC-Seq (Assay for Transposase-Accessible Chromatin) uses Tn5 transposase to map open chromatin regions genome-wide, revealing regulatory elements without antibodies. Both use similar peak calling workflows.

🎯 What is ChIP-Seq?

ChIP-Seq (Chromatin Immunoprecipitation Sequencing) combines chromatin immunoprecipitation with next-generation sequencing to identify genome-wide binding sites of proteins (transcription factors, histones) and their modifications. It reveals regulatory elements, epigenetic marks, and protein-DNA interactions — the foundation of epigenomics.

Why ChIP-Seq Matters

TF: Binding sites. Histone: Marks. Enhancers: H3K27ac. Promoters: H3K4me3. Repression: H3K27me3. Disease: Regulatory. Clinical: Epigenetics.

📌 TF Binding📌 Histone Marks📌 Enhancers📌 Peaks

# ChIP-Seq Analysis — the big picture
# Mapping protein-DNA interactions genome-wide.

# CONCEPTUAL OVERVIEW
print("""
=== ChIP-Seq Analysis ===

What it is:
  - Chromatin immunoprecipitation
  - Sequencing of bound DNA
  - Identify binding sites
  - Genome-wide map

Why it's powerful:
  - Whole genome
  - Quantitative
  - Unbiased
  - Regulatory insight

Types of ChIP-Seq:
  - TF ChIP-Seq: Sharp peaks
  - Histone ChIP-Seq: Sharp or broad
  - ATAC-Seq: Open chromatin
  - CUT&RUN: Low input
  - CUT&Tag: Low input, low background

Histone marks:
  - H3K4me3: Active promoters
  - H3K4me1: Enhancers
  - H3K27ac: Active enhancers
  - H3K27me3: Repressed
  - H3K9me3: Heterochromatin
  - H3K36me3: Gene bodies

Workflow:
  1. Experimental design
  2. ChIP + sequencing
  3. QC (FastQC)
  4. Trimming (fastp)
  5. Alignment (Bowtie2, BWA)
  6. Peak calling (MACS2)
  7. Annotation (ChIPseeker)
  8. Motif (MEME, HOMER)
  9. Differential binding (DiffBind)
  10. Visualization (IGV, deepTools)

Applications:
  - TF binding maps
  - Enhancer discovery
  - Epigenetic marks
  - Disease mechanisms
  - Drug response
""")

# WHY IT MATTERS
print("=== Why It Matters ===")
print("""
Biological questions:
  - Where does this TF bind?
  - What are the active enhancers?
  - How do histone marks change?
  - What motifs are enriched?
  - How does binding differ between conditions?

ChIP-Seq answers:
  - Peak locations
  - Peak strengths
  - Target genes
  - Motifs
  - Differential binding
""")

# HISTORICAL CONTEXT
print("=== Historical Context ===")
print("""
1980s: ChIP (chromatin immunoprecipitation)
2000s: ChIP-chip (microarrays)
2007: First ChIP-Seq (Johnson et al.)
2008: MACS (peak caller)
2010: ENCODE project
2012: DiffBind
2013: MACS2
2014: CUT&RUN
2019: CUT&Tag
2020s: Single-cell epigenomics
""")

# CHALLENGES
print("=== Challenges ===")
print("""
1. ANTIBODY QUALITY
   - Specificity
   - Efficiency
   - Validation

2. INPUT CONTROL
   - Background
   - Copy number
   - Open chromatin bias

3. PEAK CALLING
   - Narrow vs broad
   - Parameter tuning
   - False positives

4. REPRODUCIBILITY
   - Replicates
   - Batch effects
   - IDR

5. INTERPRETATION
   - Target genes
   - Motifs
   - Biological meaning
""")

🎯 Why ChIP-Seq Matters

🎯 The Value of ChIP-Seq

ChIP-Seq is the primary tool for mapping protein-DNA interactions and epigenetic marks, essential for understanding gene regulation.

Biological Insights

TF: Binding. Enhancers: Discovery. Promoters: Activity. Repression: Marks. Disease: Regulatory. Drug: Response.

🎯 TF Binding🎯 Enhancers🎯 Histone Marks🎯 Disease

# Why ChIP-Seq Matters
print("=== 1. Transcription Factor Binding ===")
print("""
TF binding:
  - Where TFs bind
  - Target genes
  - Regulatory networks

Examples:
  - p53 (tumor suppressor)
  - MYC (oncogene)
  - OCT4 (stem cells)
  - NF-κB (immune)

Output:
  - Peak locations
  - Binding motifs
  - Target genes
""")

print("\n=== 2. Enhancer Discovery ===")
print("""
Enhancers:
  - Distal regulatory elements
  - H3K4me1 (poised)
  - H3K27ac (active)
  - Cell-type specific

Workflow:
  1. H3K27ac ChIP-Seq
  2. Peak calling
  3. Super-enhancers
  4. Target genes

Tools:
  - ROSE (super-enhancers)
  - ChIPseeker
""")

print("\n=== 3. Promoter Activity ===")
print("""
Promoters:
  - Active: H3K4me3
  - Poised: H3K4me3 + H3K27me3
  - Inactive: H3K27me3

Examples:
  - Developmental genes
  - Housekeeping genes
  - Cell cycle genes
""")

print("\n=== 4. Heterochromatin ===")
print("""
Heterochromatin:
  - Repressed regions
  - H3K9me3
  - H3K27me3
  - Silencing

Examples:
  - Centromeres
  - Transposons
  - Inactive X
""")

print("\n=== 5. Differential Binding ===")
print("""
Differential binding:
  - TF binding changes
  - Histone mark changes
  - Between conditions

Examples:
  - Treatment vs control
  - Disease vs normal
  - Time course

Tools:
  - DiffBind
  - csaw
  - DBChIP
""")

print("\n=== 6. Motif Discovery ===")
print("""
Motif:
  - Sequence pattern
  - TF binding site
  - Co-factor

Tools:
  - MEME
  - HOMER
  - DREME
  - JASPAR

Output:
  - Motif logo
  - E-value
  - Enrichment
""")

print("\n=== 7. Disease Mechanisms ===")
print("""
Disease:
  - Regulatory variants
  - Enhancer hijacking
  - Epigenetic changes

Examples:
  - Cancer
  - Autoimmune
  - Neurodegenerative

Workflow:
  1. ChIP-Seq
  2. Peaks
  3. GWAS overlap
  4. Mechanism
""")

print("\n=== 8. Drug Response ===")
print("""
Drug response:
  - Epigenetic changes
  - TF binding changes
  - Resistance

Examples:
  - HDAC inhibitors
  - BET inhibitors
  - DNMT inhibitors

Workflow:
  1. Treat
  2. ChIP-Seq
  3. Differential binding
  4. Mechanism
""")

⚙️ Workflow Overview

🎯 From Reads to Peaks

The ChIP-Seq workflow has well-defined steps: QC, alignment, peak calling, annotation, and motif analysis.

Workflow

QC: FastQC. Trim: fastp. Align: Bowtie2. Peak: MACS2. Annotate: ChIPseeker. Motif: MEME. Visualize: IGV.

⚙️ QC⚙️ Align⚙️ Peak⚙️ Motif

# ChIP-Seq Workflow
print("=== Complete Pipeline ===")
print("""
1. EXPERIMENTAL DESIGN
   - Biological replicates (n ≥ 2, ideally 3)
   - Input control
   - Antibody validation
   - Spike-in (optional)

2. LIBRARY PREP
   - ChIP
   - Reverse crosslinks
   - Library prep
   - Sequencing (single or paired)

3. SEQUENCING
   - Illumina
   - 20-40M reads per sample
   - 50-75 bp single-end
   - Or 2×50 bp paired-end

4. QUALITY CONTROL
   - FastQC
   - MultiQC
   - Check: quality, GC, adapters

5. TRIMMING
   - fastp / Trimmomatic
   - Remove adapters
   - Quality filtering

6. ALIGNMENT
   - Bowtie2 (standard)
   - BWA-MEM (alternative)
   - Filter: MAPQ ≥ 30
   - Remove duplicates

7. PEAK CALLING
   - MACS2 (narrow)
   - SICER / epic (broad)
   - Input control
   - q-value < 0.05

8. PEAK ANNOTATION
   - ChIPseeker (R)
   - HOMER
   - Genomic features
   - Nearest genes

9. MOTIF DISCOVERY
   - MEME
   - HOMER
   - DREME
   - JASPAR

10. DIFFERENTIAL BINDING
    - DiffBind (R)
    - csaw (R)
    - DESeq2 / edgeR

11. VISUALIZATION
    - IGV
    - deepTools
    - ChIPseeker
""")

print("\n=== Two Main Approaches ===")
print("""
A. Narrow peaks (TF, H3K4me3, H3K27ac):
   Reads → Bowtie2 → BAM → MACS2 → narrowPeak
   - Sharp, well-defined
   - Usually < 1 kb

B. Broad peaks (H3K27me3, H3K9me3, H3K36me3):
   Reads → Bowtie2 → BAM → SICER/epic → broadPeak
   - Broad, diffuse
   - Several kb to Mb
""")

print("\n=== Key Decisions ===")
print("""
1. Input control?
   → Yes, always (background correction)

2. Paired vs single?
   → Paired-end (better for repeats)

3. Narrow vs broad?
   → TF: narrow
   → Histone: depends on mark

4. Peak caller?
   → MACS2 (narrow)
   → SICER (broad)

5. Replicates?
   → n ≥ 2, ideally 3
""")

🧹 QC & Preprocessing

🎯 Clean Before You Call Peaks

Quality control and preprocessing are critical for accurate peak calling. ChIP-Seq has specific QC metrics.

QC

FastQC: Raw reads. MultiQC: Aggregate. fastp: Trim. FRiP: Fraction in peaks. Cross-correlation: Strand shift. IDR: Reproducibility.

🧹 FastQC🧹 MultiQC🧹 FRiP🧹 IDR

# QC & Preprocessing
print("=== FastQC ===")
print("""
FastQC:
  - Per-base quality
  - GC content
  - Adapter contamination
  - Duplication
  - Overrepresented sequences

Command:
  fastqc reads.fastq.gz -o qc/

Output:
  reads_fastqc.html
  reads_fastqc.zip
""")

print("\n=== MultiQC ===")
print("""
MultiQC:
  - Aggregate all QC
  - One report
  - Compare samples

Command:
  multiqc qc/ -o multiqc_report/
""")

print("\n=== fastp ===")
print("""
fastp:
  - Fast
  - All-in-one
  - Adapter detection
  - Quality filtering

Command:
  fastp \
    -i reads.fastq.gz \
    -o clean.fastq.gz \
    --detect_adapter_for_pe \
    --qualified_quality_phred 20 \
    --length_required 36 \
    --json report.json \
    --html report.html

Key options:
  -q 20 (quality)
  -l 36 (min length)
  -u 40 (max uncalled)
""")

print("\n=== ChIP-Seq Specific QC ===")
print("""
1. FRiP (Fraction of Reads in Peaks):
   - FRiP = reads in peaks / total reads
   - > 1%: Good
   - > 5%: Excellent
   - < 1%: Problem

2. Cross-correlation (phantompeakqualtools):
   - NSC (Normalized Strand Coefficient)
   - RSC (Relative Strand Coefficient)
   - NSC > 1.05: Good
   - RSC > 0.8: Good
   - Also: Quality tag

3. IDR (Irreproducible Discovery Rate):
   - Reproducibility between replicates
   - IDR < 0.05: Good
   - Number of peaks at IDR < 0.05

4. Library complexity:
   - Number of unique fragments
   - NRF, PBC1, PBC2
   - Preseq
""")

print("\n=== phantompeakqualtools ===")
print("""
phantompeakqualtools:
  - Cross-correlation
  - NSC, RSC
  - Quality tag

Command:
  Rscript run_spp.R \
    -c=sample.bam \
    -savp=sample.pdf \
    -out=sample.spp.out

Output:
  sample.spp.out
  sample.pdf
""")

print("\n=== IDR ===")
print("""
IDR:
  - Reproducibility
  - Between replicates

Command:
  idr --samples rep1.narrowPeak rep2.narrowPeak \
    --input-file-type narrowPeak \
    --output-file idr.txt \
    --plot \
    --log-output-file idr.log

Output:
  idr.txt
  idr.png
""")

print("\n=== Best Practices ===")
print("""
1. Always run FastQC
2. Aggregate with MultiQC
3. Trim adapters
4. Check FRiP
5. Check cross-correlation
6. Check IDR
7. Check library complexity
8. Remove duplicates (optional)
9. Filter MAPQ ≥ 30
10. Document everything
""")

🗺️ Alignment

🎯 Mapping Reads to the Genome

ChIP-Seq alignment is similar to other NGS workflows, but with specific considerations for duplicates and multi-mapping reads.

Alignment

Bowtie2: Standard. BWA: Alternative. MAPQ: Filter. Duplicates: Remove. Blacklist: Remove. BAM: Output.

🗺️ Bowtie2🗺️ BWA🗺️ MAPQ🗺️ Blacklist

# Alignment
print("=== Bowtie2 ===")
print("""
Bowtie2:
  - Fast
  - Sensitive
  - Standard for ChIP-Seq

Build index:
  bowtie2-build genome.fa bowtie2_index

Align (single-end):
  bowtie2 -p 8 \
    -x bowtie2_index \
    -U clean.fastq.gz \
    -S sample.sam

Align (paired-end):
  bowtie2 -p 8 \
    -x bowtie2_index \
    -1 clean_1.fq.gz \
    -2 clean_2.fq.gz \
    -S sample.sam

Options:
  --very-sensitive (more sensitive)
  -k 1 (report 1 alignment)
  --no-mixed
  --no-discordant
""")

print("\n=== BWA-MEM ===")
print("""
BWA-MEM:
  - Fast
  - Accurate
  - Alternative

Build index:
  bwa index genome.fa

Align:
  bwa mem -t 8 genome.fa clean_1.fq.gz clean_2.fq.gz > sample.sam

Options:
  -t 8 (threads)
  -M (mark shorter split)
""")

print("\n=== Post-alignment ===")
print("""
1. Convert to BAM:
   samtools view -bS sample.sam > sample.bam

2. Sort:
   samtools sort -@ 8 -o sample.sorted.bam sample.bam

3. Index:
   samtools index sample.sorted.bam

4. Filter MAPQ:
   samtools view -q 30 -b sample.sorted.bam > sample.filtered.bam

5. Remove duplicates:
   samtools markdup -r sample.sorted.bam sample.dedup.bam
   # or
   picard MarkDuplicates \
     I=sample.sorted.bam \
     O=sample.dedup.bam \
     M=metrics.txt \
     REMOVE_DUPLICATES=true

6. Remove blacklist:
   bedtools intersect -v \
     -a sample.dedup.bam \
     -b blacklist.bed > sample.clean.bam
""")

print("\n=== Alignment Stats ===")
print("""
Check:
  - Alignment rate > 80%
  - Unique alignment > 70%
  - Duplication rate < 20%
  - Mitochondrial reads < 5%

samtools flagstat sample.bam

Output:
  total: 20,000,000
  mapped: 18,000,000 (90%)
  properly paired: 17,500,000 (87.5%)
  duplicates: 2,000,000 (10%)
""")

print("\n=== ENCODE Blacklist ===")
print("""
Blacklist:
  - Artifact regions
  - High signal
  - Remove from analysis

Download:
  wget https://github.com/Boyle-Lab/Blacklist/raw/master/lists/hg38-blacklist.v2.bed.gz

Use:
  bedtools intersect -v -a sample.bam -b blacklist.bed > clean.bam
""")

print("\n=== Best Practices ===")
print("""
1. Use Bowtie2 or BWA
2. Filter MAPQ ≥ 30
3. Remove duplicates
4. Remove blacklist
5. Check alignment stats
6. Check FRiP
7. Check cross-correlation
""")

📈 Peak Calling (MACS2/3)

🎯 Finding Binding Sites

Peak calling identifies regions of enrichment (peaks) from ChIP-Seq data. MACS2 is the gold standard for narrow peaks.

Peak Calling

MACS2: Model-based. Input: Control. q-value: FDR. Fold: Enrichment. Narrow: TF. Broad: Histone.

📈 MACS2📈 q-value📈 Fold enrichment📈 Input control

# Peak Calling with MACS2
print("=== MACS2 ===")
print("""
MACS2 (Model-based Analysis of ChIP-Seq):
  - Most popular
  - Model-based
  - Input control
  - Narrow + broad

Command (narrow peaks):
  macs2 callpeak \
    -t sample.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n sample \
    --outdir macs2_output \
    -q 0.05 \
    --nomodel \
    --extsize 147

Command (broad peaks):
  macs2 callpeak \
    -t sample.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n sample \
    --outdir macs2_output \
    --broad \
    --broad-cutoff 0.1

Options:
  -t: treatment
  -c: control (input)
  -f: format (BAM, BED)
  -g: genome size (hs, mm, ce, dm)
  -q: q-value cutoff
  --broad: broad peaks
  --nomodel: skip model building
  --extsize: extension size
""")

print("\n=== MACS2 Output ===")
print("""
Output files:
  sample_peaks.narrowPeak: Peak calls
  sample_peaks.xls: Tabular peaks
  sample_summits.bed: Summit positions
  sample_model.r: Model plot
  sample_peaks.bed: BED format

narrowPeak columns:
  chrom, start, end, name, score, strand,
  fold_enrichment, -log10(pvalue),
  -log10(qvalue), summit_position
""")

print("\n=== MACS3 ===")
print("""
MACS3:
  - Updated MACS2
  - Python 3
  - Faster
  - Same interface

Command:
  macs3 callpeak \
    -t sample.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n sample \
    -q 0.05

Install:
  pip install MACS3
""")

print("\n=== Peak Calling Parameters ===")
print("""
Genome size (-g):
  - hs: human (2.7e9)
  - mm: mouse (1.87e9)
  - ce: C. elegans (9e7)
  - dm: Drosophila (1.2e8)
  - Custom: 2.7e9

q-value (-q):
  - 0.05: Standard
  - 0.01: Stringent
  - 0.1: Lenient

Fold enrichment:
  - > 5: Strong
  - 2-5: Moderate
  - < 2: Weak

Summit:
  - Peak center
  - TF binding site
""")

print("\n=== Peak Stats ===")
print("""
Number of peaks:
  - TF: 1,000-50,000
  - H3K4me3: 10,000-30,000
  - H3K27ac: 20,000-100,000
  - H3K27me3: 1,000-10,000 (broad)

Peak width:
  - TF: 100-500 bp
  - H3K4me3: 500-2000 bp
  - H3K27ac: 500-5000 bp
  - H3K27me3: 5-100 kb
""")

print("\n=== Best Practices ===")
print("""
1. Always use input control
2. Use appropriate genome size
3. Use q-value < 0.05
4. Check peak number
5. Check FRiP
6. Check peak width
7. Use IDR for replicates
8. Remove blacklist
""")

🌊 Broad Peaks (SICER, epic)

🎯 Histone Marks Need Special Handling

Broad histone marks (H3K27me3, H3K9me3, H3K36me3) require different peak callers that account for diffuse enrichment.

Broad Peaks

SICER: Window-based. epic: HMM. Broad: Diffuse. Histone: Marks. H3K27me3: Repressed. H3K9me3: Heterochromatin.

🌊 SICER🌊 epic🌊 H3K27me3🌊 H3K9me3

# Broad Peak Calling
print("=== Why Broad Peaks? ===")
print("""
Broad marks:
  - H3K27me3: Repressed (5-100 kb)
  - H3K9me3: Heterochromatin (10-1000 kb)
  - H3K36me3: Gene bodies (1-100 kb)
  - H3K4me1: Enhancers (1-5 kb)

MACS2 limitations:
  - Designed for sharp peaks
  - May fragment broad domains
  - May miss diffuse enrichment

Solution:
  - SICER: Window-based
  - epic: HMM-based
  - Broad mode in MACS2
""")

print("\n=== SICER ===")
print("""
SICER (Spatial Clustering for Identification of ChIP-Enriched Regions):
  - Window-based
  - Designed for broad marks
  - Popular

Command:
  sicer \
    -t sample.bed \
    -c input.bed \
    -g hg38 \
    -w 1000 \
    -f 200 \
    -g 600 \
    -e 5 \
    -o sicer_output

Options:
  -w: window size (1000)
  -f: fragment size (200)
  -g: gap size (600)
  -e: E-value (5)
""")

print("\n=== epic ===")
print("""
epic (Enrichment of Peaks with Internal Control):
  - HMM-based
  - Accurate
  - Broad + narrow

Command:
  epic \
    -t sample.bed \
    -c input.bed \
    -g hg38 \
    --output epic_output

Install:
  pip install epic
""")

print("\n=== MACS2 Broad Mode ===")
print("""
MACS2 broad:
  macs2 callpeak \
    -t sample.bam \
    -c input.bam \
    -f BAM \
    -g hs \
    -n sample \
    --broad \
    --broad-cutoff 0.1

Options:
  --broad: broad peaks
  --broad-cutoff: 0.1 (q-value)
""")

print("\n=== Histone Marks Guide ===")
print(f"{'Mark':<15} {'Type':<20} {'Caller':<15} {'Width':<15}")
print("-" * 70)
print(f"{'H3K4me3':<15} {'Active promoter':<20} {'MACS2':<15} {'500-2000':<15}")
print(f"{'H3K4me1':<15} {'Enhancer':<20} {'MACS2':<15} {'1-5 kb':<15}")
print(f"{'H3K27ac':<15} {'Active enhancer':<20} {'MACS2':<15} {'500-5000':<15}")
print(f"{'H3K27me3':<15} {'Repressed':<20} {'SICER':<15} {'5-100 kb':<15}")
print(f"{'H3K9me3':<15} {'Heterochromatin':<20} {'SICER':<15} {'10-1000 kb':<15}")
print(f"{'H3K36me3':<15} {'Gene body':<20} {'SICER':<15} {'1-100 kb':<15}")

print("\n=== Best Practices ===")
print("""
1. Choose caller by mark
2. Use input control
3. Check peak width
4. Check FRiP
5. Validate with known loci
6. Use IDR for replicates
7. Visualize in IGV
""")

⚖️ Differential Binding (DiffBind)

🎯 Comparing Conditions

Differential binding identifies peaks that change between conditions — the ChIP-Seq equivalent of differential expression.

DiffBind

DiffBind: R package. csaw: Window-based. Counts: Peak matrix. DESeq2: Statistical. edgeR: Alternative. Fold: Change.

⚖️ DiffBind⚖️ csaw⚖️ Fold change⚖️ FDR

# Differential Binding
print("=== DiffBind ===")
print("""
DiffBind:
  - R package
  - Peak matrix
  - DESeq2/edgeR
  - Standard

R code:
  library(DiffBind)
  
  # Sample sheet
  samples <- data.frame(
    SampleID = c("ctrl1","ctrl2","trt1","trt2"),
    Tissue = c("control","control","treated","treated"),
    Factor = c("H3K27ac","H3K27ac","H3K27ac","H3K27ac"),
    Replicate = c(1,2,1,2),
    bamReads = c("ctrl1.bam","ctrl2.bam","trt1.bam","trt2.bam"),
    Peaks = c("ctrl1_peaks.narrowPeak","ctrl2_peaks.narrowPeak",
              "trt1_peaks.narrowPeak","trt2_peaks.narrowPeak"),
    PeakCaller = c("macs","macs","macs","macs")
  )
  
  # Load
  dba <- dba(sampleSheet=samples)
  
  # Count
  dba <- dba.count(dba)
  
  # Contrast
  dba <- dba.contrast(dba, categories=DBA_TREATMENT)
  
  # Analyze
  dba <- dba.analyze(dba)
  
  # Report
  report <- dba.report(dba)
  
  # Plot
  dba.plotMA(dba)
  dba.plotVolcano(dba)
  dba.plotHeatmap(dba)
""")

print("\n=== csaw ===")
print("""
csaw:
  - Window-based
  - No peak calling needed
  - Flexible

R code:
  library(csaw)
  
  # Load BAMs
  bam.files <- c("ctrl1.bam","ctrl2.bam","trt1.bam","trt2.bam")
  
  # Count windows
  windows <- windowCounts(bam.files, width=150, spacing=50)
  
  # Filter
  keep <- filterWindowsGlobal(windows, type="background")
  windows <- windows[keep,]
  
  # Normalize
  windows <- normFactors(windows)
  
  # Design
  design <- model.matrix(~ condition)
  
  # Test
  results <- diffBind(windows, design=design)
""")

print("\n=== Results Interpretation ===")
print("""
Output columns:
  - seqnames: Chromosome
  - start: Start
  - end: End
  - width: Width
  - strand: Strand
  - Conc: Mean concentration
  - Conc_control: Control mean
  - Conc_treated: Treated mean
  - Fold: Fold change
  - p.value: Raw p-value
  - FDR: Adjusted p-value

Thresholds:
  - FDR < 0.05
  - |Fold| > 1 (2-fold)
  - Or |Fold| > 0.58 (1.5-fold)
""")

print("\n=== Visualization ===")
print("""
Plots:
  - MA plot: Fold vs mean
  - Volcano: Fold vs significance
  - Heatmap: Peak matrix
  - PCA: Sample clustering
  - Boxplot: Peak distribution
""")

print("\n=== Best Practices ===")
print("""
1. Use biological replicates (n ≥ 2)
2. Use input control
3. Normalize properly
4. Use consistent peak calling
5. Multiple testing correction
6. Effect size threshold
7. Validate key changes
8. Visualize results
""")

📝 Peak Annotation

🎯 What Genes Are Near My Peaks?

Peak annotation assigns genomic features and nearest genes to peaks, revealing regulatory targets.

Annotation

ChIPseeker: R package. HOMER: CLI. Feature: Promoter, exon. Nearest: Gene. Distance: TSS. Enrichment: GO, pathway.

📝 ChIPseeker📝 HOMER📝 Nearest gene📝 GO

# Peak Annotation
print("=== ChIPseeker ===")
print("""
ChIPseeker:
  - R package
  - Comprehensive
  - Visualization
  - Standard

R code:
  library(ChIPseeker)
  library(TxDb.Hsapiens.UCSC.hg38.knownGene)
  library(org.Hs.eg.db)
  
  # Load peaks
  peaks <- readPeakFile("sample_peaks.narrowPeak")
  
  # Annotate
  txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
  anno <- annotatePeak(peaks, tssRegion=c(-3000, 3000),
                       TxDb=txdb, annoDb="org.Hs.eg.db")
  
  # Results
  as.data.frame(anno)
  
  # Plot
  plotAnnoPie(anno)
  plotAnnoBar(anno)
  plotDistToTSS(anno)
  
  # GO enrichment
  library(clusterProfiler)
  genes <- as.data.frame(anno)$geneId
  ego <- enrichGO(genes, OrgDb=org.Hs.eg.db, ont="BP")
  dotplot(ego)
""")

print("\n=== HOMER ===")
print("""
HOMER:
  - CLI
  - Fast
  - Motif + annotation

Command:
  annotatePeaks.pl \
    sample_peaks.narrowPeak \
    hg38 \
    > sample_annotated.txt

Options:
  -gtf annotation.gtf
  -go GO_analysis
  -genome hg38
""")

print("\n=== Genomic Features ===")
print("""
Features:
  - Promoter (< 3 kb from TSS)
  - 5' UTR
  - Exon
  - Intron
  - 3' UTR
  - Intergenic

Distribution:
  - TF: Often promoter + enhancer
  - H3K4me3: Promoter
  - H3K27ac: Enhancer
  - H3K27me3: Intergenic
""")

print("\n=== Nearest Gene ===")
print("""
Nearest gene:
  - Distance to TSS
  - Gene symbol
  - Strand

Output:
  - Gene ID
  - Gene name
  - Distance
  - Feature
""")

print("\n=== Functional Enrichment ===")
print("""
GO:
  - BP: Biological process
  - MF: Molecular function
  - CC: Cellular component

Pathways:
  - KEGG
  - Reactome

Tools:
  - clusterProfiler
  - DAVID
  - Enrichr
  - GREAT
""")

print("\n=== Best Practices ===")
print("""
1. Use appropriate TSS region
2. Check feature distribution
3. Check nearest genes
4. Use multiple databases
5. Visualize distribution
6. Validate key targets
7. Report parameters
""")

🔍 Motif Discovery (MEME, HOMER)

🎯 What Sequence Patterns Are Enriched?

Motif discovery identifies enriched sequence patterns in peaks, revealing TF binding sites and co-factors.

Motif

MEME: De novo. HOMER: Known + de novo. DREME: Short motifs. JASPAR: Database. Logo: Visualization. E-value: Significance.

🔍 MEME🔍 HOMER🔍 JASPAR🔍 Motif logo

# Motif Discovery
print("=== MEME ===")
print("""
MEME:
  - De novo motif discovery
  - Comprehensive suite
  - Standard

Command:
  meme peaks.fa \
    -dna \
    -oc meme_output \
    -nmotifs 10 \
    -minw 6 \
    -maxw 20 \
    -mod zoops

Options:
  -dna: DNA sequences
  -nmotifs: Number (10)
  -minw: Min width (6)
  -maxw: Max width (20)
  -mod: Distribution (zoops, anr, oops)
""")

print("\n=== HOMER ===")
print("""
HOMER:
  - Known + de novo
  - Fast
  - Popular

Command:
  findMotifsGenome.pl \
    sample_peaks.narrowPeak \
    hg38 \
    homer_output \
    -size 200 \
    -mask

Options:
  -size: Peak size (200)
  -mask: Repeat masking
  -bg: Background
""")

print("\n=== DREME ===")
print("""
DREME:
  - Short motifs (4-8 bp)
  - Fast
  - Part of MEME

Command:
  dreme -p peaks.fa -o dreme_output
""")

print("\n=== JASPAR ===")
print("""
JASPAR:
  - Database of TF motifs
  - Known motifs
  - Compare de novo

Tools:
  - Tomtom (compare)
  - FIMO (scan)

Command (Tomtom):
  tomtom meme_output/meme.txt \
    JASPAR2020_CORE_vertebrates.meme \
    -oc tomtom_output
""")

print("\n=== Motif Output ===")
print("""
Output:
  - Motif logo (PWM)
  - E-value
  - Number of sites
  - Enrichment
  - Best match

Example:
  Motif 1: E-value 1e-50
  Match: CTCF (JASPAR)
  Sites: 5000
""")

print("\n=== Motif Logo ===")
print("""
Logo:
  - Height = information content
  - Letter size = frequency
  - Shows conservation

Tools:
  - WebLogo
  - Skylign
  - ggseqlogo
""")

print("\n=== Best Practices ===")
print("""
1. Use appropriate peak size
2. Use repeat masking
3. Use good background
4. Multiple tools
5. Compare to JASPAR
6. Check E-value
7. Validate key motifs
8. Visualize logos
""")

🐍 Python Integration

🎯 ChIP-Seq from Python

You can run ChIP-Seq tools from Python using subprocess and parse outputs with pandas and pybedtools.

Python

subprocess: Run CLI. pysam: BAM. pybedtools: BED. pandas: Peaks. matplotlib: Plot. Pipeline: Automate.

🐍 subprocess🐍 pysam🐍 pybedtools🐍 pandas

# Python Integration
print("=== Run Bowtie2 ===")
print("""
import subprocess

def run_bowtie2(index, reads, output_bam, threads=8):
    \"\"\"Run Bowtie2 alignment.\"\"\"
    cmd = (
        f"bowtie2 -p {threads} -x {index} -U {reads} "
        f"| samtools view -bS - "
        f"| samtools sort -@ {threads} -o {output_bam} -"
    )
    
    result = subprocess.run(cmd, shell=True, capture_output=True, text=True)
    
    if result.returncode != 0:
        raise RuntimeError(f"Bowtie2 failed: {result.stderr}")
    
    # Index BAM
    subprocess.run(['samtools', 'index', output_bam], check=True)
    
    return output_bam

# Use
# bam = run_bowtie2("bowtie2_index", "reads.fq.gz", "sample.bam")
""")

print("\n=== Run MACS2 ===")
print("""
def run_macs2(treatment_bam, control_bam, output_dir,
              genome="hs", qvalue=0.05, broad=False):
    \"\"\"Run MACS2 peak calling.\"\"\"
    cmd = [
        'macs2', 'callpeak',
        '-t', treatment_bam,
        '-c', control_bam,
        '-f', 'BAM',
        '-g', genome,
        '-n', 'sample',
        '--outdir', output_dir,
        '-q', str(qvalue)
    ]
    
    if broad:
        cmd.extend(['--broad', '--broad-cutoff', '0.1'])
    
    result = subprocess.run(cmd, capture_output=True, text=True)
    
    if result.returncode != 0:
        raise RuntimeError(f"MACS2 failed: {result.stderr}")
    
    return f"{output_dir}/sample_peaks.narrowPeak"

# Use
# peaks = run_macs2("sample.bam", "input.bam", "macs2_output")
""")

print("\n=== Parse Peaks ===")
print("""
import pandas as pd

def parse_peaks(peak_file):
    \"\"\"Parse narrowPeak file.\"\"\"
    columns = [
        'chrom', 'start', 'end', 'name', 'score', 'strand',
        'fold_enrichment', 'neg_log10_pvalue',
        'neg_log10_qvalue', 'summit'
    ]
    
    peaks = pd.read_csv(peak_file, sep='\\t', header=None,
                        names=columns)
    
    print(f"Peaks: {len(peaks)}")
    print(f"Mean width: {peaks['end'].sub(peaks['start']).mean():.0f} bp")
    print(f"Mean fold: {peaks['fold_enrichment'].mean():.2f}")
    
    return peaks

# Use
# peaks = parse_peaks("sample_peaks.narrowPeak")
""")

print("\n=== pybedtools ===")
print("""
import pybedtools

def analyze_peaks(peak_file, blacklist_file):
    \"\"\"Analyze peaks with pybedtools.\"\"\"
    peaks = pybedtools.BedTool(peak_file)
    blacklist = pybedtools.BedTool(blacklist_file)
    
    # Remove blacklist
    clean = peaks.intersect(blacklist, v=True)
    
    # Statistics
    print(f"Original peaks: {peaks.count()}")
    print(f"Clean peaks: {clean.count()}")
    print(f"Total bp: {clean.total_coverage()}")
    
    return clean

# Use
# clean = analyze_peaks("sample_peaks.narrowPeak", "blacklist.bed")
""")

print("\n=== Peak Annotation ===")
print("""
import pybedtools

def annotate_peaks(peak_file, gtf_file):
    \"\"\"Annotate peaks with genes.\"\"\"
    peaks = pybedtools.BedTool(peak_file)
    genes = pybedtools.BedTool(gtf_file)
    
    # Intersect
    annotated = peaks.intersect(genes, wa=True, wb=True)
    
    # Parse
    results = []
    for feature in annotated:
        results.append({
            'peak_chrom': feature[0],
            'peak_start': int(feature[1]),
            'peak_end': int(feature[2]),
            'gene_chrom': feature[9],
            'gene_start': int(feature[10]),
            'gene_end': int(feature[11]),
            'gene_name': feature[15] if len(feature) > 15 else 'N/A'
        })
    
    return pd.DataFrame(results)

# Use
# anno = annotate_peaks("peaks.narrowPeak", "genes.gtf")
""")

print("\n=== Motif Analysis ===")
print("""
import subprocess

def run_meme(peak_file, output_dir, genome="hg38"):
    \"\"\"Run MEME motif discovery.\"\"\"
    # Extract sequences
    fasta = f"{output_dir}/peaks.fa"
    subprocess.run([
        'bedtools', 'getfasta',
        '-fi', f"{genome}.fa",
        '-bed', peak_file,
        '-fo', fasta
    ], check=True)
    
    # Run MEME
    subprocess.run([
        'meme', fasta,
        '-dna',
        '-oc', f"{output_dir}/meme",
        '-nmotifs', '10',
        '-minw', '6',
        '-maxw', '20',
        '-mod', 'zoops'
    ], check=True)
    
    return f"{output_dir}/meme/meme.txt"

# Use
# motifs = run_meme("peaks.narrowPeak", "motif_output")
""")

print("\n=== Complete Pipeline ===")
print("""
import subprocess
import os

def chipseq_pipeline(reads, input_reads, genome_index,
                     genome, output_dir="chipseq_output",
                     threads=8):
    \"\"\"Complete ChIP-Seq pipeline.\"\"\"
    os.makedirs(output_dir, exist_ok=True)
    
    # 1. Trim
    print("Trimming...")
    subprocess.run([
        'fastp', '-i', reads,
        '-o', f'{output_dir}/clean.fq.gz',
        '--detect_adapter_for_pe'
    ], check=True)
    
    # 2. Align
    print("Aligning...")
    bam = run_bowtie2(genome_index,
                      f'{output_dir}/clean.fq.gz',
                      f'{output_dir}/sample.bam', threads)
    
    # 3. Align input
    print("Aligning input...")
    input_bam = run_bowtie2(genome_index, input_reads,
                            f'{output_dir}/input.bam', threads)
    
    # 4. Peak calling
    print("Calling peaks...")
    peaks = run_macs2(bam, input_bam, f'{output_dir}/macs2',
                      genome=genome)
    
    # 5. Parse
    peak_df = parse_peaks(peaks)
    
    print("Pipeline complete!")
    return peak_df

# Use
# peaks = chipseq_pipeline("sample.fq.gz", "input.fq.gz",
#                          "bowtie2_index", "hs")
""")

🎨 Visualization

🎯 Seeing the Peaks

Visualization is critical for ChIP-Seq interpretation: genome browser tracks, heatmaps, and profile plots.

Visualization

IGV: Genome browser. deepTools: Heatmaps. ChIPseeker: Annotation. UCSC: Browser. Profile: TSS. Heatmap: Peaks.

🎨 IGV🎨 deepTools🎨 ChIPseeker🎨 Heatmap

# Visualization
print("=== IGV ===")
print("""
IGV (Integrative Genomics Viewer):
  - Genome browser
  - Interactive
  - Standard

Load:
  - BAM (coverage)
  - BigWig (coverage)
  - BED (peaks)
  - GTF (genes)

Command line:
  igv.sh -g hg38 sample.bam sample_peaks.narrowPeak
""")

print("\n=== deepTools ===")
print("""
deepTools:
  - Heatmaps
  - Profile plots
  - Correlation
  - Standard

bamCoverage:
  bamCoverage -b sample.bam \
    -o sample.bw \
    --normalizeUsing RPKM \
    --binSize 10

computeMatrix:
  computeMatrix reference-point \
    -S sample.bw \
    -R genes.bed \
    -a 3000 -b 3000 \
    -o matrix.gz

plotHeatmap:
  plotHeatmap -m matrix.gz \
    -o heatmap.png \
    --colorMap RdBu_r

plotProfile:
  plotProfile -m matrix.gz \
    -o profile.png
""")

print("\n=== ChIPseeker Plots ===")
print("""
R code:
  library(ChIPseeker)
  
  # Pie chart
  plotAnnoPie(anno)
  
  # Bar chart
  plotAnnoBar(anno)
  
  # Distance to TSS
  plotDistToTSS(anno)
  
  # Profile
  plotProfile(anno)
  
  # Heatmap
  plotHeatmap(anno)
""")

print("\n=== UCSC Genome Browser ===")
print("""
UCSC:
  - Web browser
  - Custom tracks
  - Public data

Upload:
  - BED
  - BigWig
  - BAM (via URL)

URL:
  https://genome.ucsc.edu/
""")

print("\n=== Profile Plot ===")
print("""
Profile:
  - Average signal
  - Around TSS
  - Around peaks

Tools:
  - deepTools
  - ChIPseeker
  - ngsplot
""")

print("\n=== Heatmap ===")
print("""
Heatmap:
  - Signal around regions
  - Clustering
  - Patterns

Tools:
  - deepTools
  - pheatmap
  - ComplexHeatmap
""")

print("\n=== Best Practices ===")
print("""
1. Use BigWig for coverage
2. Use IGV for inspection
3. Use deepTools for heatmaps
4. Use ChIPseeker for annotation
5. Colorblind-friendly
6. High resolution
7. Vector format
8. Document parameters
""")

⚖️ Tool Comparison

🎯 Which Tool?

Choosing the right ChIP-Seq tools depends on your target, peak type, and analysis goals.

Comparison

Bowtie2: Aligner. MACS2: Narrow. SICER: Broad. DiffBind: Differential. MEME: Motif. ChIPseeker: Annotation.

⚖️ Bowtie2⚖️ MACS2⚖️ SICER⚖️ DiffBind

Tool	Type	Speed	Accuracy	Best For
Bowtie2	Aligner	Fast	Excellent	Standard alignment
BWA-MEM	Aligner	Fast	Excellent	Alternative
MACS2	Peak caller	Fast	Excellent	Narrow peaks (TF)
MACS3	Peak caller	Fast	Excellent	Updated MACS2
SICER	Peak caller	Medium	Very good	Broad peaks (histone)
epic	Peak caller	Medium	Excellent	Broad + narrow
DiffBind	Differential	Medium	Excellent	Differential binding
csaw	Differential	Fast	Very good	Window-based
MEME	Motif	Slow	Excellent	De novo motifs
HOMER	Motif	Fast	Very good	Known + de novo
ChIPseeker	Annotation	Fast	Excellent	Peak annotation

# Decision Guide
print("=== Choose Peak Caller ===")
print("""
Transcription factors (sharp):
  → MACS2 / MACS3
  → narrowPeak

H3K4me3, H3K27ac (sharp):
  → MACS2 / MACS3
  → narrowPeak

H3K27me3, H3K9me3 (broad):
  → SICER
  → epic
  → MACS2 --broad

H3K36me3 (gene body):
  → SICER
  → epic

ATAC-Seq:
  → MACS2
  → Genrich
  → HMMRATAC
""")

print("\n=== Choose Aligner ===")
print("""
Bowtie2:
  ✓ Standard
  ✓ Fast
  ✓ Sensitive
  ✓ Good for most

BWA-MEM:
  ✓ Alternative
  ✓ Good for longer reads
  ✓ SV detection

Recommendation:
  → Bowtie2 for most ChIP-Seq
  → BWA for ATAC-Seq
""")

print("\n=== Speed Comparison ===")
print(f"{'Tool':<20} {'10M reads':<15} {'50M reads':<15}")
print("-" * 50)
print(f"{'Bowtie2':<20} {'2 min':<15} {'10 min':<15}")
print(f"{'BWA-MEM':<20} {'3 min':<15} {'15 min':<15}")
print(f"{'MACS2':<20} {'1 min':<15} {'5 min':<15}")
print(f"{'MACS3':<20} {'1 min':<15} {'5 min':<15}")
print(f"{'SICER':<20} {'5 min':<15} {'25 min':<15}")
print(f"{'DiffBind':<20} {'2 min':<15} {'10 min':<15}")

print("\n=== When to Use Which ===")
print("""
USE MACS2:
  ✓ TF ChIP-Seq
  ✓ H3K4me3
  ✓ H3K27ac
  ✓ Narrow peaks

USE SICER:
  ✓ H3K27me3
  ✓ H3K9me3
  ✓ H3K36me3
  ✓ Broad peaks

USE DiffBind:
  ✓ Differential binding
  ✓ Multiple conditions

USE MEME:
  ✓ De novo motifs
  ✓ Novel TFs

USE HOMER:
  ✓ Known motifs
  ✓ Fast annotation
""")

💼 Use Cases

🎯 Where ChIP-Seq Is Used

ChIP-Seq is used in virtually every area of biology: gene regulation, epigenetics, cancer, development, and more.

Use Cases

TF: Binding. Enhancers: Discovery. Histone: Marks. Cancer: Epigenetics. Development: Differentiation. Drug: Response.

💼 TF Binding💼 Enhancers💼 Cancer💼 Development

# Use Cases
print("=== 1. TF Binding Maps ===")
print("""
Use case:
  - Where TFs bind
  - Target genes
  - Regulatory networks

Workflow:
  1. ChIP-Seq (TF)
  2. Peak calling (MACS2)
  3. Annotation (ChIPseeker)
  4. Motif (MEME)
  5. Network

Examples:
  - p53
  - MYC
  - OCT4
  - NF-κB

Tools:
  - MACS2
  - ChIPseeker
  - MEME
""")

print("\n=== 2. Enhancer Discovery ===")
print("""
Use case:
  - Active enhancers
  - Super-enhancers
  - Cell-type specific

Workflow:
  1. H3K27ac ChIP-Seq
  2. Peak calling (MACS2)
  3. Super-enhancers (ROSE)
  4. Target genes

Examples:
  - Cancer
  - Development
  - Immune

Tools:
  - MACS2
  - ROSE
  - ChIPseeker
""")

print("\n=== 3. Histone Modifications ===")
print("""
Use case:
  - Active/repressed
  - Chromatin states
  - Epigenetic changes

Workflow:
  1. Histone ChIP-Seq
  2. Peak calling (MACS2/SICER)
  3. ChromHMM
  4. Differential

Examples:
  - H3K4me3
  - H3K27ac
  - H3K27me3
  - H3K9me3

Tools:
  - MACS2
  - SICER
  - ChromHMM
""")

print("\n=== 4. Cancer Epigenetics ===")
print("""
Use case:
  - Epigenetic changes
  - Enhancer hijacking
  - Drug targets

Workflow:
  1. Tumor ChIP-Seq
  2. Peak calling
  3. Differential binding
  4. Targets

Examples:
  - Enhancer hijacking
  - Super-enhancers
  - Epigenetic drugs

Tools:
  - MACS2
  - DiffBind
  - ROSE
""")

print("\n=== 5. Development ===")
print("""
Use case:
  - Differentiation
  - Cell fate
  - Time course

Workflow:
  1. Time points
  2. ChIP-Seq
  3. Differential binding
  4. Trajectory

Examples:
  - Stem cells
  - Embryonic
  - Organoids

Tools:
  - MACS2
  - DiffBind
  - ChromHMM
""")

print("\n=== 6. Drug Response ===")
print("""
Use case:
  - Mechanism
  - Resistance
  - Biomarkers

Workflow:
  1. Treat
  2. ChIP-Seq
  3. Differential binding
  4. Mechanism

Examples:
  - HDAC inhibitors
  - BET inhibitors
  - DNMT inhibitors

Tools:
  - MACS2
  - DiffBind
""")

print("\n=== 7. Immune Cells ===")
print("""
Use case:
  - TF networks
  - Activation
  - Differentiation

Workflow:
  1. ChIP-Seq
  2. Peaks
  3. Motifs
  4. Networks

Examples:
  - T cells
  - B cells
  - Macrophages

Tools:
  - MACS2
  - HOMER
  - MEME
""")

print("\n=== 8. Stem Cells ===")
print("""
Use case:
  - Pluripotency
  - Differentiation
  - Reprogramming

Workflow:
  1. ChIP-Seq
  2. Peaks
  3. Super-enhancers
  4. Targets

Examples:
  - OCT4
  - SOX2
  - NANOG

Tools:
  - MACS2
  - ROSE
  - ChIPseeker
""")

print("\n=== Comparison ===")
print(f"{'Use Case':<25} {'Tool':<20} {'Output':<20}")
print("-" * 65)
print(f"{'TF binding':<25} {'MACS2':<20} {'Peaks':<20}")
print(f"{'Enhancers':<25} {'MACS2 + ROSE':<20} {'Super-enhancers':<20}")
print(f"{'Histone marks':<25} {'SICER':<20} {'Broad peaks':<20}")
print(f"{'Cancer':<25} {'DiffBind':<20} {'Differential':<20}")
print(f"{'Development':<25} {'ChromHMM':<20} {'States':<20}")
print(f"{'Drug':<25} {'DiffBind':<20} {'Changes':<20}")
print(f"{'Immune':<25} {'HOMER':<20} {'Motifs':<20}")
print(f"{'Stem cells':<25} {'ROSE':<20} {'Super-enhancers':<20}")

❓ FAQ

🎯 Frequently Asked Questions

Answers to common questions about ChIP-Seq analysis.

What is ChIP-Seq?

What is the difference between transcription factor and histone modification ChIP-Seq?

Transcription factor (TF) ChIP-Seq produces sharp, narrow peaks (100-500 bp) at specific binding sites. Histone modification ChIP-Seq produces broader, often bimodal peaks: active marks like H3K4me3 are sharp, while repressive marks like H3K27me3 and H3K9me3 are broad (several kb). Different peak callers are needed.

Which peak caller should I use?

Use MACS2 or MACS3 for transcription factors and sharp histone marks (H3K4me3, H3K27ac). Use SICER or epic for broad histone marks (H3K27me3, H3K9me3, H3K36me3). For differential binding, use DiffBind or csaw. For ATAC-Seq, MACS2 works well.

What is a control (input) sample in ChIP-Seq?

An input control is chromatin that has not been immunoprecipitated. It's used to correct for background, copy number, and open chromatin biases. MACS2 uses the input to model local background and calculate fold enrichment. Without input, peak calling is less accurate.

What is the difference between ChIP-Seq and ATAC-Seq?

ChIP-Seq uses antibodies to pull down specific proteins or histone modifications, giving targeted binding maps. ATAC-Seq (Assay for Transposase-Accessible Chromatin) uses Tn5 transposase to map open chromatin regions genome-wide, revealing regulatory elements without antibodies. Both use similar peak calling workflows.

How many reads do I need for ChIP-Seq?

For TF ChIP-Seq: 20-40 million reads. For histone marks: 20-40 million. For broad marks: 40-60 million. For input control: 20-40 million. More reads improve peak detection, especially for weak binding sites.

What is FRiP and why does it matter?

FRiP (Fraction of Reads in Peaks) is the proportion of reads that fall within called peaks. It's a key quality metric: > 1% is good, > 5% is excellent, < 1% indicates a problem. Low FRiP suggests poor antibody, insufficient signal, or wrong peak calling parameters.

What is IDR and why is it important?

IDR (Irreproducible Discovery Rate) measures reproducibility between replicates. It's the ENCODE standard for ChIP-Seq. Use IDR < 0.05 to define a reproducible peak set. It's more robust than simple peak overlap and is required for publication-quality data.

How do I find motifs in my peaks?

Use MEME for de novo motif discovery or HOMER for known + de novo motifs. Extract peak sequences with bedtools getfasta, then run MEME/HOMER. Compare discovered motifs to JASPAR using Tomtom. Check E-values and enrichment.

What is differential binding analysis?

Differential binding identifies peaks that change between conditions (e.g., treatment vs control). Use DiffBind or csaw in R. They count reads in a consensus peak set and use DESeq2/edgeR to test for significant changes. Requires biological replicates (n ≥ 2, ideally 3).

# Quick Reference
print("=== Quick Answers ===")
print("""
Q: What is ChIP-Seq?
A: Map protein-DNA interactions genome-wide.

Q: Best aligner?
A: Bowtie2.

Q: Best peak caller (TF)?
A: MACS2 / MACS3.

Q: Best peak caller (broad)?
A: SICER / epic.

Q: Input control?
A: Yes, always.

Q: How many reads?
A: 20-40M per sample.

Q: What is FRiP?
A: Fraction of reads in peaks.

Q: What is IDR?
A: Reproducibility metric.

Q: Motif discovery?
A: MEME, HOMER.

Q: Differential binding?
A: DiffBind.

Q: Free tools?
A: Yes — all major tools are free.
""")

🚀 Future

🎯 Where ChIP-Seq Is Going

The future is low-input, single-cell, multi-omics, AI-augmented, and clinical-grade.

Future Directions

CUT&RUN: Low input. CUT&Tag: Low background. Single-cell: Per-cell. Multi-omics: Integration. AI: Prediction. Clinical: Epigenetics.

🚀 CUT&RUN🚀 CUT&Tag🚀 Single-cell🚀 AI

# Future of ChIP-Seq
print("=== 1. CUT&RUN & CUT&Tag ===")
print("""
Trend: Low-input epigenomics

CUT&RUN:
  - Cleavage Under Targets and Release Using Nuclease
  - Low input (100-1000 cells)
  - Low background
  - High resolution

CUT&Tag:
  - Cleavage Under Targets and Tagmentation
  - Tn5 transposase
  - Very low input
  - Low background
  - Single-cell compatible

Advantages:
  - Fewer cells
  - Less sequencing
  - Better signal
  - Single-cell

Tools:
  - MACS2
  - SEACR
  - GoPeaks
""")

print("\n=== 2. Single-Cell Epigenomics ===")
print("""
Trend: Per-cell epigenomics

Current:
  - Bulk ChIP-Seq
  - Average signal

Single-cell:
  - scATAC-Seq
  - scChIP-Seq
  - scCUT&Tag
  - Per-cell profiles

Tools:
  - ArchR
  - Signac
  - SnapATAC
  - Scanpy

Future:
  - Cell-type specific
  - Rare populations
  - Developmental
""")

print("\n=== 3. Multi-Omics Integration ===")
print("""
Trend: Integrate everything

Current:
  - Single modality
  - Separate analysis

Multi-omics:
  - ChIP-Seq + RNA-Seq
  - ChIP-Seq + ATAC-Seq
  - ChIP-Seq + Hi-C
  - Single-cell multi-omics

Tools:
  - Seurat
  - MOFA
  - mixOmics
  - ArchR

Future:
  - Unified models
  - Regulatory networks
  - Systems biology
""")

print("\n=== 4. AI-Augmented Analysis ===")
print("""
Trend: Deep learning

Ideas:
  - Peak calling with CNNs
  - Motif discovery with transformers
  - Enhancer prediction
  - Chromatin state

Examples:
  - Enformer (DeepMind)
  - Basenji
  - Sei
  - ChromBPNet

Future:
  - Learned representations
  - Zero-shot prediction
  - Variant effect
  - Clinical
""")

print("\n=== 5. Clinical Epigenomics ===")
print("""
Trend: Routine clinical

Current:
  - Research
  - Specialized

Future:
  - Clinical routine
  - Cancer diagnostics
  - Epigenetic drugs
  - Biomarkers

Examples:
  - Enhancer hijacking
  - Super-enhancers
  - Methylation
  - Histone marks

Requirements:
  - Accuracy
  - Speed
  - Cost
  - Regulation
""")

print("\n=== 6. Spatial Epigenomics ===")
print("""
Trend: Location + epigenetics

Current:
  - Dissociated cells
  - Lose location

Spatial:
  - Tissue context
  - Spatial epigenomics
  - Spatial ATAC
  - Spatial CUT&Tag

Platforms:
  - Spatial-ATAC
  - Spatial-CUT&Tag
  - Slide-tags

Future:
  - 3D epigenomics
  - Tissue atlases
  - Clinical
""")

# TIMELINE
print("\n=== Timeline ===")
print("""
Past:
  2007: First ChIP-Seq
  2008: MACS
  2010: ENCODE
  2012: DiffBind
  2013: MACS2
  2014: CUT&RUN
  2019: CUT&Tag
  2020: Single-cell ATAC

Present:
  2025: Multi-omics
  2026: AI-augmented

Future:
  2030: Clinical routine?
  2035: Spatial epigenomics?
  2040: Predictive epigenomics?
""")

print("\n=== Final Thoughts ===")
print("""
ChIP-Seq is:
  - Fundamental to epigenomics
  - Needed for gene regulation
  - Needed for disease
  - Rapidly improving

The future is:
  - Low-input
  - Single-cell
  - Multi-omics
  - AI-augmented
  - Clinical-grade
  - Spatial

"ChIP-Seq maps the regulatory genome.
 Analysis is how we read it."

Keep exploring. Keep building.
""")