bio-applied-chipseq-pipeline

name	bio-applied-chipseq-pipeline
description	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.
tool_type	bash
primary_tool	MACS2/MACS3

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.

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

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

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)
     out.loc[, ]  out.loc[, ]
    (out[[, , ]])

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.

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

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