bio-chipseq-spike-in-normalization

What this skill does

## Version Compatibility

Reference examples tested with: DiffBind 3.20+, DESeq2 1.42+, edgeR 4.0+, csaw 1.36+, ChIPseqSpikeInFree 1.6+, SpikChIP 1.0+, SpikeFlow (NAR Genom Bioinform 2024), samtools 1.19+, bowtie2 2.5+.

# ChIP-seq Spike-In Normalization

**"Account for global signal changes that defeat standard normalization"** -> Add exogenous reference chromatin (Drosophila for human/mouse ChIP-Rx; E. coli carryover for CUT&RUN/CUT&Tag) at fixed concentration BEFORE IP, derive scaling factors from spike-in read counts, and apply at the read or size-factor level (never to peak counts) to enable quantitative cross-condition comparison.

- CLI: align reads to combined target + spike genome; count spike reads via `samtools view -c`
- R (DiffBind integration): `dba.normalize(obj, spikein = TRUE)`
- R (DESeq2 / edgeR): `sizeFactors(dds) <- 1 / scale_factors` (note inverse)
- CLI (deepTools tracks): `bamCoverage --scaleFactor <derived>` (mutually exclusive with `--normalizeUsing`)
- Wrapper: SpikeFlow (Snakemake; 2024) automates end-to-end
- Post-hoc detection: ChIPseqSpikeInFree (when no spike-in was added)

The fundamental rule: spike-in scaling is applied at the READ level (via size factors or `--scaleFactor`), never multiplied into peak counts. ~25% of published spike-in ChIP papers violate this (Patel L et al 2024 *Nat Biotechnol* 42:1343).

## When Spike-In Is Required

| Experimental design | Spike-in needed? |
|---------------------|------------------|
| HDAC inhibitor -> global H3K27ac increase | Yes |
| BET inhibitor (JQ1, OTX015) -> global BRD4 / H3K27ac decrease | Yes |
| EZH2 inhibitor -> global H3K27me3 loss | Yes |
| DNMT inhibitor -> global 5mC loss; downstream histone mark shifts | Yes |
| Target factor knockdown / degron | Yes (or matched-input subtraction) |
| Cell-cycle synchronization / arrest | Yes |
| Dosage titration | Yes |
| Standard TF perturbation, local rebinding expected | No (reads-in-peaks RLE works) |
| Histone mark cross-cell-type comparison | Recommended |
| CUT&RUN/CUT&Tag standard | E. coli carryover (automatic); deliberate Drosophila for high-stakes |
| Replicate-only experiment, no condition comparison | No |

**Why this is necessary:** Standard normalization (RLE on reads-in-peaks, TMM on bins) assumes most regions don't change. When the perturbation IS the change-everything-globally biology, these methods force the median log2FC to zero, hiding the real effect.

## Spike-In Protocol Taxonomy

| Protocol | Spike organism | Added when | Notes |
|----------|----------------|------------|-------|
| **ChIP-Rx** (Orlando 2014) | Drosophila S2 nuclei | After lysis, before IP | 50k Drosophila nuclei per 5M target cells (Egan 2016) |
| **ChIP-Rx variant** (Bonhoure 2014) | Drosophila chromatin | After fragmentation, before IP | Different normalization layer |
| **CUT&RUN/Tag E. coli** | E. coli (carryover) | Automatic from bacterial pA-MNase/Tn5 | Free; variable across enzyme batches |
| **Heterologous spike-in** (Hu 2015) | Defined yeast / E. coli chromatin | Added at lysis | Less common; defined concentration |
| **xenoChIP** | Species swap (mouse cells + human chromatin spike) | Before IP | Niche; specific cancer xenograft contexts |

**The dominant standard for human/mouse ChIP is Drosophila (ChIP-Rx).** Drosophila is genetically distinct enough that mapping is unambiguous, and the genome size (~140 Mb) gives adequate read depth at small chromatin input.

## Scaling Factor Calculation

**RRPM (Orlando 2014):** reference-adjusted reads per million.

```
scale_factor_i = min(N_spike) / N_spike_i
```

Apply at the read level. The sample with the fewest spike reads gets scale_factor = 1; others get > 1 (compensating for less observed spike chromatin).

**Rx-Input (Fursova 2019):** additionally scales by input spike-in to correct IP efficiency variation.

```
RxInput_i = (N_spike_chip_i / N_total_chip_i) / (N_spike_input_i / N_total_input_i)
```

This is more rigorous when input controls are available; required for some inhibitor experiments where IP efficiency itself changes.

## Workflow: Drosophila ChIP-Rx Spike-In

**Goal:** Compute per-sample scaling factors from Drosophila spike-in reads and apply at the read level (not peak counts) to enable quantitative cross-condition ChIP-seq comparison.

**Approach:** Align reads to combined target + Drosophila genome, count spike reads at high mapq after deduplication, derive RRPM scaling factors (min/each), then apply via DESeq2 sizeFactors, DiffBind spike-in flag, or bamCoverage scaleFactor. Validate against internal-control regions (blacklist).

### Step 1: Alignment to combined genome

```bash
# Build combined index (target + Drosophila)
cat hg38.fa dm6.fa > hg38_dm6.fa
bowtie2-build hg38_dm6.fa hg38_dm6

# Align reads
bowtie2 -x hg38_dm6 -1 R1.fq -2 R2.fq -S aln.sam --very-sensitive --no-mixed
samtools view -bS aln.sam | samtools sort -o aln.bam
samtools index aln.bam
```

### Step 2: Filter, deduplicate, count spike reads

```bash
# Apply ENCODE filter (-F 1804 -q 30) BEFORE counting spike reads
samtools view -F 1804 -q 30 -b aln.bam > aln.filt.bam
samtools index aln.filt.bam

# Count Drosophila reads (NOT total reads)
DROSO_READS=$(samtools view -c aln.filt.bam chr2L chr2R chr3L chr3R chr4 chrX chrY)
echo "$SAMPLE: Drosophila reads = $DROSO_READS"

# Separate into target-only BAM for peak calling
samtools view -b aln.filt.bam chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 \
    chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22 chrX chrY \
    > aln.filt.hg38.bam
samtools index aln.filt.hg38.bam
```

### Step 3: Compute scaling factors

```bash
# Per-sample Drosophila counts (assume saved in droso_counts.tsv)
# sample_id, droso_reads
# ctrl_1, 145000
# ctrl_2, 132000
# treat_1, 98000
# treat_2, 85000

awk 'BEGIN{min=1e10} NR>1{if($2<min) min=$2} END{print "min:", min}' droso_counts.tsv
# Use min as numerator: scale_factor_i = min / droso_reads_i
```

### Step 4: Apply scaling — three layers

**Layer 1: bigWig tracks**

```bash
SCALE=$(echo "scale=6; $MIN_DROSO / $SAMPLE_DROSO" | bc)
bamCoverage -b sample.bam -o sample.scaled.bw \
    --scaleFactor $SCALE --binSize 10 --extendReads 200
# DO NOT also pass --normalizeUsing; mutually exclusive
```

**Layer 2: DiffBind**

```r
library(DiffBind)
dba_obj <- dba(sampleSheet = 'samples.csv')   # spike-in BAM in sample sheet
dba_obj <- dba.count(dba_obj, summits = 250, bParallel = TRUE)

# Spike-in normalization
dba_obj <- dba.normalize(dba_obj, spikein = TRUE,
                          library = DBA_LIBSIZE_FULL,
                          normalize = DBA_NORM_LIB)

# Verify what was applied
dba.normalize(dba_obj, bRetrieve = TRUE)
```

**Layer 3: DESeq2 / edgeR direct**

```r
library(DESeq2)
# Read spike-in counts into a vector aligned with sample order
spike_reads <- c(ctrl_1 = 145000, ctrl_2 = 132000, treat_1 = 98000, treat_2 = 85000)
scale_factors <- min(spike_reads) / spike_reads

dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ condition)
# DESeq2 expects sizeFactors in INVERSE convention (sample with smallest factor gets largest sizeFactor)
sizeFactors(dds) <- 1 / scale_factors
dds <- DESeq(dds, fitType = 'parametric')
```

## Workflow: E. coli Spike-In (CUT&RUN/CUT&Tag Automatic)

E. coli DNA from bacterial pA-MNase/pA-Tn5 production is automatic spike-in carryover.

```bash
# Combined index
cat hg38.fa ecoli_k12.fa > hg38_ecoli.fa
bowtie2-build hg38_ecoli.fa hg38_ecoli

# Align as in ChIP-Rx; count E. coli reads
ECOLI_READS=$(samtools view -c aln.filt.bam ecoli_chr1)
TOTAL_READS=$(samtools view -c aln.filt.bam)
echo "E. coli fraction: $(echo "scale=4; $ECOLI_READS / $TOTAL_READS" | bc)"
# Target: 0.005-0.02 (0.5-2%); IgG: 0.02-0.05 (2-5%)

# Scale factor same as ChIP-Rx: min(ecoli) / per_sample_ecoli
# Apply at read or sizeFactors level
```

**E. coli carryover is variable** between enzyme production batches. For publication-grade cross-condition claims, supplement with deliberate Droso