# Bio Applied Chipseq Pipeline

> FASTQ-to-peaks ChIP-seq pipeline: Bowtie2 align, Picard dedup, MACS2/MACS3 narrow/broad peak calling, FRiP/IDR QC, deepTools bamCoverage/heatmaps. Use for ChIP-seq/CUT&RUN peak calling or FRiP/NRF/IDR QC.

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

---


# ChIP-seq Processing Pipeline

## When to Use

- Processing raw ChIP-seq (or CUT&RUN/CUT&Tag) FASTQ files into aligned, deduplicated BAMs
- Calling narrow peaks (TFs, active marks) or broad peaks (repressive histone marks) with MACS2/MACS3
- Assessing library/experiment quality: FRiP, NRF/PBC, fingerprint plots, replicate IDR
- Generating normalized signal tracks (BigWig) and TSS/peak-centered heatmaps with deepTools
- Deciding read depth, replicate count, or peak-shape expectations before designing an experiment

## Version Compatibility

Bowtie2 ≥2.5, samtools ≥1.19, Picard ≥3.1, MACS3 ≥3.0 (MACS2 ≥2.2 still widely used and interchangeable for standard `callpeak`), deepTools ≥3.5, bedtools ≥2.31. Reference: ENCODE ChIP-seq pipeline v2 conventions (hg38, `hg38-blacklist.v2.bed`).

## Prerequisites

- Tools on `PATH`: `fastqc`, `multiqc`, `trim_galore` (or `fastp`), `bowtie2`, `samtools`, `bedtools`, `picard`, `macs3` (or `macs2`), `deeptools` (`bamCoverage`, `computeMatrix`, `plotHeatmap`, `plotFingerprint`)
- Python: `pandas`, `numpy`, `matplotlib`, optionally `pyBigWig` to read BigWig values directly
- Prior concepts: paired-end FASTQ QC, BAM coordinate filtering (see `bio-read-qc-fastp-workflow`, `bio-read-alignment-bowtie2-alignment`)

## Experiment Design Reference

| Experiment Type | Peak Shape | Recommended Reads | Example Marks |
|---|---|---|---|
| TF | Narrow (< 500 bp) | 20-40 M | CTCF, GATA1, p53 |
| Active histone | Narrow | 30-50 M | H3K4me3, H3K9ac |
| Broad histone | Broad (> 5 kb) | 40-80 M | H3K27me3, H3K9me3 |
| Enhancer mark | Narrow/Mixed | 30-50 M | H3K27ac, H3K4me1 |

**Antibody**: must be ChIP-grade (not just WB-validated). **Replicates**: ENCODE requires ≥2 biological replicates + IDR analysis. **Sequencing**: paired-end strongly preferred (accurate fragment size, better dedup).

## Step 1: QC, Trimming, Alignment, Dedup

**Goal:** turn raw paired-end FASTQ into a clean, deduplicated, blacklist-filtered BAM.
**Approach:** FastQC/MultiQC → Trim Galore (or fastp) → Bowtie2 with strict pairing/MAPQ filters → blacklist removal → Picard `MarkDuplicates`.

```bash
# --- QC ---
fastqc sample_R1.fastq.gz sample_R2.fastq.gz -o qc/raw/ -t 4
multiqc qc/raw/ -o qc/raw/multiqc/

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

# --- Align (Bowtie2), filter unmapped/secondary/low-MAPQ, keep properly paired ---
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/sample_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
samtools flagstat aligned/sample.bam | tee logs/sample_flagstat.txt

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

# --- Mark/remove PCR + optical duplicates (2500 px for NovaSeq, 100 for HiSeq) ---
picard MarkDuplicates \
    I=aligned/sample_filtered.bam O=dedup/sample_dedup.bam \
    M=dedup/sample_dup_metrics.txt REMOVE_DUPLICATES=true \
    OPTICAL_DUPLICATE_PIXEL_DISTANCE=2500 \
    VALIDATION_STRINGENCY=SILENT
samtools index dedup/sample_dedup.bam
```

Targets: alignment rate > 80%; filter flags `-F 4` (unmapped), `-F 256` (secondary), `-q 30` (MAPQ), `-f 2` (properly paired). ENCODE dedup QC: PBC1/NRF > 0.8 acceptable, > 0.9 ideal; total dup rate should be < 30% for TF ChIP.

## Step 2: Peak Calling (MACS3) and FRiP QC

**Goal:** call narrow peaks for TFs/active marks or broad peaks for repressive histone marks against matched input, then quantify signal-to-noise with FRiP.
**Approach:** `macs3 callpeak` with `--broad` only for marks like H3K27me3/H3K9me3; compute FRiP from `bedtools intersect` read counts vs. total mapped reads.

```bash
# --- TF / narrow mark ---
macs3 callpeak -t dedup/sample_dedup.bam -c dedup/input_dedup.bam \
    -f BAMPE -g hs -n tf_sample --outdir peaks/tf/ -q 0.05 --keep-dup all \
    2> logs/macs3_tf.log

# --- Broad histone mark (H3K27me3, H3K9me3) ---
macs3 callpeak -t dedup/h3k27me3_dedup.bam -c dedup/input_dedup.bam \
    -f BAMPE -g hs -n h3k27me3 --outdir peaks/broad/ -q 0.05 --broad \
    --keep-dup all 2> logs/macs3_broad.log

# --- Reads-in-peaks for FRiP (numerator) ---
bedtools intersect -a dedup/sample_dedup.bam -b peaks/tf/tf_sample_peaks.narrowPeak -u \
    | samtools view -c - > reads_in_peaks.txt
```

```python
import pandas as pd


def compute_frip(reads_in_peaks: int, total_mapped: int) -> float:
    """Fraction of Reads in Peaks — core ChIP-seq enrichment QC metric.

    ENCODE thresholds: TF ChIP >= 0.05, histone-mark ChIP >= 0.01.
    """
    if total_mapped <= 0:
        raise ValueError("total_mapped must be > 0")
    return reads_in_peaks / total_mapped


def qc_summary(df: pd.DataFrame) -> pd.DataFrame:
    """Add FRiP and pass/fail flags to a per-sample QC table.

    df columns required: Sample, Type ('TF'|'Histone'|'Input'),
    Mapped_Reads, Reads_in_Peaks.
    """
    df = df.copy()
    df["FRiP"] = df["Reads_in_Peaks"] / df["Mapped_Reads"]
    threshold = df["Type"].map({"TF": 0.05, "Histone": 0.01}).fillna(0)
    df["FRiP_pass"] = df["FRiP"] >= threshold
    return df


if __name__ == "__main__":
    demo = pd.DataFrame({
        "Sample": ["CTCF_rep1", "H3K27ac_rep1"],
        "Type": ["TF", "Histone"],
        "Mapped_Reads": [33_800_000, 45_600_000],
        "Reads_in_Peaks": [8_400_000, 22_300_000],
    })
    out = qc_summary(demo)
    assert out.loc[0, "FRiP_pass"] and out.loc[1, "FRiP_pass"]
    print(out[["Sample", "FRiP", "FRiP_pass"]])
```

## Step 3: Signal Tracks and TSS Heatmaps (deepTools)

**Goal:** convert BAM to normalized BigWig and visualize enrichment around TSS or peak summits.
**Approach:** `bamCoverage` (RPKM/CPM normalization) → `computeMatrix reference-point` → `plotHeatmap`/`plotProfile`; use `plotFingerprint` to sanity-check ChIP vs. input enrichment.

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

plotFingerprint -b dedup/sample_dedup.bam dedup/input_dedup.bam \
    --labels ChIP Input -plot figures/fingerprint.png

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

plotHeatmap -m matrix/tss_matrix.gz -out figures/heatmap_tss.png \
    --colorMap Blues --whatToShow 'heatmap and colorbar' --zMin 0 --zMax 10 \
    --samplesLabel ChIP Input --regionsLabel Genes
plotProfile -m matrix/tss_matrix.gz -out figures/profile_tss.png \
    --samplesLabel ChIP Input --plotTitle 'Signal around TSS (+/-3 kb)' --perGroup
```

```python
import numpy as np


def read_matrix_gz(path: str) -> np.ndarray:
    """Load a deepTools computeMatrix .gz output into a (regions x bins) array.

    Requires deeptools' own gzip/JSON header format; for ad-hoc BigWig
    reads instead use pyBigWig.open(path).values(chrom, start, end, numpy=True).
    """
    import gzip
    import json

    with gzip.open(path, "rt") as fh:
        header = json.loads(fh.readline().lstrip("@"))
        data = np.loadtxt(fh, delimiter="\t", usecols=range(6, 6 + sum(header["sample_boundaries"][1:])))
    return data


def frip_from_narrowpeak(narrowpeak_bed: str, bam_read_count_in_peaks: int, total_mapped: int) -> float:
    """Convenience wrapper: FRiP given a precomputed reads-in-peaks count."""
    return bam_read_count_in_peaks / total_mapped
```

## Pitfalls

- **Coordinate systems**: BED/narrowPeak are 0-based half-open; VCF/GFF are 1-based inclusive — mixing them causes off-by-one errors
- **Input control is mandatory**: never call peaks without a matched IgG/input control; FRiP and fold-enrichment are meaningless without it
- **`--broad` mismatch**: using default (narrow) settings on H3K27me3/H3K9me3 fragments broad domains into many spurious narrow peaks
- **Batch effects**: check for batch confounding (e.g., different sequencing runs per replicate) before interpreting biological signal
- **Multiple testing**: MACS q-values already apply BH-FDR per run, but comparing peaks across many samples/marks still needs correction downstream
- **Blacklist skipped**: leaving ENCODE blacklist regions in inflates peak counts at repetitive/artifact-prone loci

## See Also

- `bio-chip-seq-peak-calling` — deeper MACS2/MACS3 parameter tuning
- `bio-chip-seq-chipseq-qc` — full FRiP/NRF/PBC/cross-correlation QC suite
- `bio-chip-seq-differential-binding` — comparing peaks/signal across conditions
- `bio-chip-seq-peak-annotation` — annotating peaks to genes/genomic features