/bio-chipseq-qc

name: bio-chipseq-qc
description: ChIP-seq quality control metrics including FRiP (Fraction of Reads in Peaks), cross-correlation analysis (NSC/RSC), library complexity, and IDR (Irreproducibility Discovery Rate) for replicate concordance. Use to assess experiment quality before downstream analysis. Use when assessing ChIP-seq data quality metrics.
tool_type: mixed
primary_tool: deepTools

Version Compatibility

Reference examples tested with: MACS3 3.0+, Subread 2.0+, bedtools 2.31+, deepTools 3.5+, pybedtools 0.9+, pysam 0.22+, samtools 1.19+

Before using code patterns, verify installed versions match. If versions differ:

• Python: `pip show <package>` then `help(module.function)` to check signatures
• R: `packageVersion('<pkg>')` then `?function_name` to verify parameters
• CLI: `<tool> --version` then `<tool> --help` to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

ChIP-seq Quality Control

**"Assess the quality of my ChIP-seq experiment"** → Compute FRiP, cross-correlation (NSC/RSC), library complexity, and IDR replicate concordance to evaluate enrichment success.

• CLI: `deeptools plotFingerprint`, `phantompeakqualtools run_spp.R`
• Python: `pysam` + `pybedtools` for custom QC metrics

Quality metrics for assessing ChIP-seq experiment success and replicate reproducibility.

FRiP (Fraction of Reads in Peaks)

**Goal:** Quantify enrichment strength by measuring the proportion of reads falling within called peaks.

**Approach:** Count reads overlapping peak regions and divide by total mapped reads.

Calculate FRiP with bedtools

# Count reads in peaks
reads_in_peaks=$(bedtools intersect -a chip.bam -b peaks.narrowPeak -u | samtools view -c -)
total_reads=$(samtools view -c -F 260 chip.bam)

# Calculate FRiP
frip=$(echo "scale=4; $reads_in_peaks / $total_reads" | bc)
echo "FRiP: $frip"

Calculate FRiP with featureCounts

# Convert peaks to SAF format
awk 'BEGIN{OFS="\t"} {print $4, $1, $2, $3, "."}' peaks.narrowPeak > peaks.saf

# Count reads in peaks
featureCounts -a peaks.saf -F SAF -o peak_counts.txt chip.bam

# FRiP from summary
grep -v "^#" peak_counts.txt.summary

Calculate FRiP with pysam

import pysam
import pybedtools

def calculate_frip(bam_file, peak_file):
    bam = pysam.AlignmentFile(bam_file, 'rb')
    total_reads = bam.count(read_callback=lambda r: not r.is_unmapped and not r.is_secondary)

    peaks = pybedtools.BedTool(peak_file)
    reads_in_peaks = 0
    for peak in peaks:
        reads_in_peaks += bam.count(peak.chrom, peak.start, peak.end)

    frip = reads_in_peaks / total_reads
    return frip

frip = calculate_frip('chip.bam', 'peaks.narrowPeak')
print(f'FRiP: {frip:.4f}')

FRiP Thresholds

| Target | Minimum FRiP | Good FRiP | |--------|--------------|-----------| | TF (narrow) | 0.01 | > 0.05 | | Histone (broad) | 0.10 | > 0.20 | | H3K4me3 | 0.05 | > 0.15 | | H3K27ac | 0.05 | > 0.10 |

Cross-Correlation Analysis (NSC/RSC)

**Goal:** Assess ChIP enrichment quality by measuring strand cross-correlation signal.

**Approach:** Calculate correlation between forward and reverse strand read coverage at varying shifts to detect fragment-length enrichment.

Run phantompeakqualtools

# Run SPP cross-correlation analysis
Rscript run_spp.R \
    -c=chip.bam \
    -savp=chip_cc.pdf \
    -out=chip_cc.txt \
    -odir=qc/

# Output columns:
# 1: filename
# 2: numReads
# 3: estFragLen (estimated fragment length)
# 4: corr_estFragLen
# 5: phantomPeak
# 6: corr_phantomPeak
# 7: argmin_corr (minimum strand shift)
# 8: min_corr
# 9: NSC (Normalized Strand Coefficient)
# 10: RSC (Relative Strand Coefficient)
# 11: QualityTag

Interpret NSC and RSC

# Parse results
awk -F'\t' '{
    print "Fragment length:", $3
    print "NSC:", $9
    print "RSC:", $10
    print "Quality:", $11
}' chip_cc.txt

NSC/RSC Thresholds

| Metric | Marginal | Acceptable | Ideal | |--------|----------|------------|-------| | NSC | < 1.05 | 1.05 - 1.1 | > 1.1 | | RSC | < 0.8 | 0.8 - 1.0 | > 1.0 | | QualityTag | -2 | 0 | 1 or 2 |

Plot Cross-Correlation in R

library(spp)

chip_data <- read.bam.tags('chip.bam')
binding_characteristics <- get.binding.characteristics(chip_data, srange=c(50, 500), bin=5)

# Cross-correlation plot
pdf('cc_plot.pdf')
plot(binding_characteristics$cross.correlation, type='l',
     xlab='Strand shift', ylab='Cross-correlation')
abline(v=binding_characteristics$peak$x, col='red')
dev.off()

# Extract metrics
print(paste('Fragment length:', binding_characteristics$peak$x))

Library Complexity (NRF, PBC1, PBC2)

**Goal:** Detect PCR amplification artifacts by measuring library complexity metrics.

**Approach:** Calculate the fraction of unique reads and positional redundancy to assess PCR bottlenecking.

Calculate with bedtools

# NRF: Non-Redundant Fraction (unique reads / total reads)
total=$(samtools view -c -F 260 chip.bam)
unique=$(samtools view -F 260 chip.bam | cut -f1-4 | sort -u | wc -l)
nrf=$(echo "scale=4; $unique / $total" | bc)
echo "NRF: $nrf"

# PBC1: PCR Bottleneck Coefficient 1 (M1/Mdistinct)
# M1 = locations with exactly 1 read
# Mdistinct = distinct genomic locations

bedtools bamtobed -i chip.bam | \
    awk '{print $1":"$2"-"$3}' | \
    sort | uniq -c | \
    awk '{
        if($1==1) m1++
        mdist++
    } END {
        print "M1:", m1
        print "Mdistinct:", mdist
        print "PBC1:", m1/mdist
    }'

Library Complexity Thresholds

| Metric | Severe | Mild | None | |--------|--------|------|------| | NRF | < 0.5 | 0.5 - 0.8 | > 0.8 | | PBC1 | < 0.5 | 0.5 - 0.8 | > 0.8 | | PBC2 | < 1 | 1 - 3 | > 3 |

IDR (Irreproducibility Discovery Rate)

**Goal:** Assess replicate concordance by measuring consistency of ranked peak lists.

**Approach:** Compare signal-ranked peaks from two repli