Bio Chipseq Qc

Version Compatibility

Reference examples tested with: deepTools 3.5+, phantompeakqualtools 1.2.2+, ChIPQC 1.42+, IDR 2.0.4+, samtools 1.19+, bedtools 2.31+, pysam 0.22+, pybedtools 0.9+, MACS2 2.2.9+, MACS3 3.0.4+.

Verify versions before relying on numerical thresholds — phantompeakqualtools has known R-version compatibility issues with R ≥ 4.0 (use kundajelab fork or pin to R 3.6).

ChIP-seq Quality Control

"Should I trust this ChIP-seq experiment?" -> Validate antibody, fragmentation, enrichment, replicate concordance, library complexity, and absence of hyper-ChIPable artifacts before committing to downstream peak calling and differential analysis.

• CLI: Rscript run_spp.R -c=chip.bam -out=cc.txt (NSC/RSC), plotFingerprint -b chip.bam input.bam (enrichment shape), idr --samples rep1.np rep2.np (replicate IDR)
• R: ChIPQC package (Carroll & Stark; computes the full ENCODE metric battery)
• Python: pysam + pybedtools for custom FRiP and library-complexity metrics

ChIP-seq fails for many independent reasons. The QC metrics below probe distinct failure modes — passing one metric does not rescue another. Antibody failure cannot be fixed by sequencing more.

The Antibody Problem is the Real Problem

Every downstream metric is conditional on antibody specificity. "ChIP-grade" on a vendor datasheet is marketing, not validation. Run the cascade:

Histone modification cross-reactivity is universal: H3K9me2 vs H3K9me3, H3K27me2 vs H3K27me3, and H3K4me1 vs H3K4me2 antibodies routinely show 10-30% cross-reactivity. Polyclonals vary lot-to-lot. CRISPR-knockout-validated lots from CST and Epicypher are the modern standard. Always record antibody catalog number + lot in methods.

Fragment-Size Distribution is a Free Diagnostic

The fragment-size distribution from a properly prepared ChIP BAM is itself a quality readout:

# Quick diagnostic — count fragment sizes from properly-paired reads
samtools view -f 0x2 sample.bam | awk '{print $9}' | awk '$1>0' \
    | sort -n | uniq -c | awk '{print $2, $1}' > fragment_sizes.tsv

For CUT&Tag: 25-75 bp characteristic; fragments < 25 bp are Tn5 self-tagmentation noise (see cut-and-run-tag).

QC Metric Battery with ENCODE Thresholds

Practical operational rule: Compute the full battery. Reject any sample failing FRiP OR antibody validation OR fragment-size sanity check, regardless of other metrics. Failing one of NSC/RSC alone with strong FRiP can sometimes be rescued for narrow-peak biology; broad histones are more forgiving on NSC.

Hyper-ChIPable Region Detection

Teytelman 2013 (PNAS): untagged GFP, no antibody, or non-existent targets all produce "binding" signal at highly-transcribed loci (rRNA, tRNA, histone gene clusters, snoRNA hosts, mtDNA, abundant housekeeping genes). ENCODE blacklist v2 (Amemiya 2019) catches repeat-driven artifacts but NOT these hyper-ChIPable transcribed regions.

Detection:

# Top 1% input signal as cell-type-specific custom blacklist
multiBigwigSummary BED-file -b input.bw -o input_signal.npz \
    --BED genes.bed --outRawCounts input_per_gene.tsv
awk 'NR > 1' input_per_gene.tsv | sort -k4,4nr | head -n $(($(wc -l < input_per_gene.tsv) / 100)) \
    > hyper_chipable.bed

# Intersect peaks against this list; flag peaks falling in hyper-ChIPable regions
bedtools intersect -a peaks.narrowPeak -b hyper_chipable.bed -u > suspicious_peaks.bed

Disprove a suspicious peak: Required for any claim at rRNA loci, tRNA clusters, HIST1/2 clusters, mitochondrial DNA:

1. Motif enrichment at peak (artifact has no enrichment)
2. KO/KD signal loss at peak (artifact persists)
3. Untagged-protein control ChIP shows no signal at this locus

Many "novel binding" claims at the rDNA repeat, mtDNA, and HIST1 cluster are spurious artifacts.

Computing the Battery

FRiP

total_reads=$(samtools view -c -F 260 chip.bam)
reads_in_peaks=$(bedtools intersect -a chip.bam -b peaks.narrowPeak -u | samtools view -c -)
frip=$(echo "scale=4; $reads_in_peaks / $total_reads" | bc)

NSC / RSC / fragment length (phantompeakqualtools)

Rscript run_spp.R -c=chip.bam -savp=qc/chip_cc.pdf -out=qc/chip_cc.txt
# Output columns: filename | numReads | estFragLen | corr_estFragLen |
#                 phantomPeak | corr_phantomPeak | argmin_corr | min_corr |
#                 NSC | RSC | QualityTag

Library complexity

# NRF
total=$(samtools view -c -F 260 chip.bam)
unique=$(samtools view -F 260 chip.bam | awk '{print $1, $3, $4}' | sort -u | wc -l)
nrf=$(echo "scale=4; $unique / $total" | bc)

# PBC1, PBC2 (singletons vs distinct positions vs doubletons)
samtools view -F 260 chip.bam | awk '{print $3":"$4}' | sort | uniq -c \
    | awk '{
        if ($1 == 1) m1++;
        if ($1 == 2) m2++;
        mdist++;
      } END {
        print "M1:", m1; print "M2:", m2; print "Mdistinct:", mdist;
        print "PBC1:", m1/mdist; print "PBC2:", m1/m2
      }'

deepTools plotFingerprint

plotFingerprint \
    -b chip.bam input.bam \
    --labels ChIP Input \
    -o qc/fingerprint.pdf \
    --outRawCounts qc/fingerprint_counts.tab \
    --outQualityMetrics qc/fingerprint_qc.txt
# Inspect qc/fingerprint_qc.txt: AUC, JS distance, synthetic JS, X-intercept
# Good ChIP: AUC < 0.6 (TF), JS > 0.3 (TF); Input near diagonal (AUC ~ 0.5)

Replicate Spearman correlation

multiBamSummary bins -b rep1.bam rep2.bam rep3.bam input.bam \
    --binSize 10000 -o results.npz
plotCorrelation -in results.npz --corMethod spearman \
    --whatToPlot heatmap --plotNumbers -o corr.pdf \
    --outFileCorMatrix corr_matrix.tab
# Replicates: > 0.8 (narrow), > 0.6 (broad)
# Cross-condition reps should correlate less than within-condition

ChIPQC R package

library(ChIPQC)
samples <- read.csv('samples.csv')
qc <- ChIPQC(samples, annotation = 'hg38')
ChIPQCreport(qc, reportFolder = 'ChIPQCreport')
# Generates the full ENCODE battery report per sample in one call

ChIPQC remains Bioconductor-maintained but mature; phantompeakqualtools is the canonical NSC/RSC source.

IDR and Replicate Consistency Rules

For TFs: signal-ranked IDR with Nself/Nt consistency check; see chip-seq/peak-calling for full ENCODE workflow. Key thresholds:

• True replicate IDR threshold: 0.05
• Pseudoreplicate IDR threshold: 0.10 (per-rep self-consistency)
• Nself/Nt rule: max(N1self, N2self) / min(N1self, N2self) ≤ 2 AND max(Nt, max(Nself)) / min(Nt, min(Nself)) ≤ 2. Failing both ratios rejects the library.

For histones: naive overlap with ≥ 40% reciprocal overlap (ENCODE default; commonly misquoted as 50%). IDR is too conservative for histone signal dynamic range.

ENCODE 3 vs ENCODE 4 Thresholds (unchanged for most QC)

Most QC thresholds are stable across ENCODE versions; blacklist update is the main practical change.

Per-Tool Failure Modes

phantompeakqualtools / SPP -- R version incompatibility

Trigger: Running with R ≥ 4.0.

Mechanism: spp R package has unmaintained Boost / Rcpp dependencies; some shifts produce NaN cross-correlation values.

Symptom: NSC = NaN, RSC = NaN, or fragment length = 0 in output.

Fix: Pin to R 3.6 + spp 1.16 via conda env; OR use the kundajelab/phantompeakqualtools fork (current); OR substitute deepTools plotFingerprint for enrichment QC and macs3 predictd for fragment length.

deepTools plotFingerprint -- Wrong baseline assumption for broad marks

Trigger: Interpreting JS distance with TF threshold (> 0.3) on broad histone mark.

Mechanism: Broad marks have less concentrated signal; JS distance is naturally lower (0.05-0.15 for H3K27me3) without indicating bad ChIP.

Symptom: Reports "failed JS distance" for high-quality broad-mark ChIP.

Fix: Use mark-specific thresholds: > 0.3 for TFs and sharp histones; > 0.05 for broad histones; check AUC instead (0.6-0.9 for broad; < 0.6 for TF/sharp).