GPTomics/bio-clip-seq-differential-clip

Version Compatibility

Reference examples tested with: DEWSeq 1.18+, htseq-clip 2.0+, DESeq2 1.44+, edgeR 4.2+, limma 3.60+, Flipper (commit 2024.04+), Skipper (commit 2023.05+), pybedtools 0.10+, pyranges 0.0.129+.

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

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

If code throws unexpected errors, introspect the installed package and adapt the example to match the actual API rather than retrying.

Differential CLIP-seq Analysis

"Identify regions with changed RBP binding across conditions" -> Test for condition-level differences in IP enrichment relative to SMInput, accounting for replicate variance and (where available) sequencing depth normalization. Three statistical scales are possible: peak-level (test each peak as a unit), window-level (test fixed transcriptome windows; DEWSeq, Flipper), or crosslink-site level (test single-nt positions). The choice depends on the biology (narrow regulatory shift vs broad binding-mode change) and on which upstream peak caller was used (CLIPper -> peak-level; Skipper -> window-level; PureCLIP -> CL-level).

• R (window-level, DEWSeq + htseq-clip): library(DEWSeq); dds <- DESeqDataSetFromSlidingWindows(counts, colData, design=~condition); dds <- DESeq(dds); res <- results(dds)
• R (Skipper downstream, Flipper): flipper differential -i skipper_out/ --design design.tsv --contrast treatment vs control -o flipper_out/
• R (peak-level, edgeR): dge <- DGEList(counts=peak_counts, group=condition); dge <- calcNormFactors(dge); design <- model.matrix(~condition); dge <- estimateDisp(dge, design); fit <- glmQLFit(dge, design); res <- glmQLFTest(fit)
• R (peak-level, limma-voom): v <- voom(dge, design); fit <- lmFit(v, design); fit <- eBayes(fit); res <- topTable(fit, coef=2, number=Inf)
• CLI (htseq-clip preprocessing): htseq-clip extract -i annotation.gff -o annotation_windows.bed -w 50 -s 20 && htseq-clip count -i sample.bam -w annotation_windows.bed -o sample.counts.txt

DEWSeq is the EMBL/Hentze-group windowed-NB approach designed specifically for differential CLIP. Flipper is the Skipper-companion tool (Tu 2024) for the modern Skipper workflow. Peak-level edgeR/limma-voom work when the upstream peak caller produced a comparable peak BED across conditions (e.g., consensus peaks from CLIPper).

Algorithmic Taxonomy

| Tool | Scale | Statistical model | Replicate requirement | Strength | Fails when | |------|-------|-------------------|----------------------|----------|------------| | DEWSeq (Schwarzl 2020) | 50-100 nt sliding window | Negative binomial GLM (DESeq2 internals) | >= 2 reps per condition | Designed specifically for CLIP; integrates SMInput | Slow on dense libraries; output window-resolution | | Flipper (Tu 2024) | Skipper window (100 nt) | Negative binomial; designed for Skipper output | >= 2 reps per condition | Modern; pairs with Skipper peak caller | Only useful if upstream is Skipper | | ASpeak | Peak | Negative binomial | >= 2 reps | Peak-level for CLIPper output | Less popular; legacy | | edgeR (general) | Peak or window | Quasi-likelihood F-test on NB | >= 2 reps | Mature, widely cited | Generic; not CLIP-aware; needs careful normalization | | limma-voom | Peak or window | Linear model with mean-variance trend | >= 2 reps | Fast; well-validated; handles small samples | Generic; treats counts as continuous after voom | | DESeq2 (direct) | Peak or window | Negative binomial GLM | >= 2 reps | Mature; same engine as DEWSeq | Same as edgeR caveats | | MASTR-seq (Banks 2023) | Cell-resolved (scCLIP) | Cell-mixture model | Cells | Single-cell CLIP differential | Specialized; few published applications | | diffbind (CLIP adaptation) | Peak | DESeq2 or edgeR backend | >= 2 reps | Familiar from ChIP/ATAC | Designed for ChIP; needs CLIP-specific normalization | | MAnorm2 | Peak | Hierarchical empirical Bayes | >= 2 reps | Tested on ChIP; less on CLIP | Less CLIP-specific |

Methodology evolves; verify DEWSeq vignette and Flipper paper for current best practice. The DEWSeq + htseq-clip pipeline (Schwarzl 2020) is the most-published CLIP-specific differential framework; Flipper (Tu 2024) is the modern Skipper-coupled alternative.

Critical Decision: The Interaction-Term Design

Use ~ type + condition + type:condition and test the interaction-term coefficient. This is the single most consequential statistical choice in differential CLIP.

• type = ip vs sminput (whether the library is IP or size-matched input)
• condition = treated vs untreated (or KD vs control)
• type:condition interaction = "Does the IP-vs-input ratio shift across conditions?"

A naive ~ condition design tests whether read counts differ regardless of whether they come from IP or SMInput - this confounds binding changes with expression changes. The interaction term explicitly tests for differential binding (the biologically meaningful signal) rather than differential expression at peak loci.

When testing, extract the interaction-term coefficient: results(dds, name = 'typeip.conditiontreat'). The log2FoldChange returned is the change in IP/input ratio in treated vs untreated - this is what "differential binding" means.

Critical Choice: Peak-Level vs Window-Level vs Crosslink-Level

Three scales differ in resolution and statistical power:

Peak-level (edgeR, limma-voom, DESeq2 on CLIPper peaks): Test each consensus peak's IP/SMInput log2 FC across conditions. Pro: peak boundaries are biologically meaningful; output interpretable. Con: peak set changes across conditions (a new peak in treatment but missing in control complicates testing); SMInput normalization must be applied consistently.

Window-level (DEWSeq, Flipper): Tile transcriptome into fixed 50-100 nt windows; test each window. Pro: comparable across conditions (windows are pre-defined); handles binding-mode shifts within a peak; high statistical power. Con: window-resolution; multiple-testing burden (millions of windows); biological meaning of a window needs translation.

Crosslink-level (custom): Test each single-nt CL position. Pro: nucleotide resolution; captures motif-level shifts. Con: very low coverage per position; massive multiple-testing burden; rarely used in published differential CLIP.

| Goal | Scale | Tool | |------|-------|------| | ENCODE-style peak-level differential (CLIPper upstream) | Peak | DESeq2 / edgeR / limma-voom on CLIPper consensus peaks | | Maximum sensitivity windowed differential | Window | DEWSeq (with htseq-clip) or Flipper (with Skipper) | | Modern Skipper-coupled workflow | Window | Flipper | | Single-cell scCLIP differential | Cell | MASTR-seq or scCLAM (sparse) | | Compare binding-mode shifts within a peak | Window | DEWSeq | | Allele-specific differential | CL site | BEAPR + custom logistic | | RBP-KD effect on binding profile | Peak/window | DEWSeq (handles KD-effect on RBP itself) |

DEWSeq Workflow (Window-Level Differential)

DEWSeq is the EMBL/Hentze-group differential framework specifically designed for CLIP. The pipeline is:

Goal: Identify transcriptome windows where the IP-vs-SMInput ratio shifts across conditions, accounting for replicate variance with the negative-binomial GLM.

Approach: Use htseq-clip to extract sliding 50 nt windows across annotated features, count reads per window per sample, build a DESeqDataSetFromSlidingWindows object with the ~ type + condition + type:condition interaction design, extract the typeip.conditiontreat interaction coefficient as the differential-binding effect size, and aggregate adjacent significant windows with bedtools merge -d 100.

# Step 1: htseq-clip generates sliding-window count matrices
htseq-clip extract \
    -i gencode.v38.annotation.gff \
    -o annotation_windows.bed \
    --window-size 50 \
    --window-step 20 \
    --feature-type CDS,UTR

# Step 2: count IP and SMInput reads per window per sample
for sample in ip_rep1 ip_rep2 sminput_rep1 sminput_rep2; do
    htseq-clip count \
        -i ${sample}.dedup.bam \
        -a annotation_windows.bed \
        -o ${sample}.counts.txt \
        --mate 2
done
# (For eCLIP, --mate 2 because R2 5' is the truncation site; for iCLIP single-end use --mate 1)

# Step 3: DEWSeq differential testing
htseq-clip mergeCounts \
    -i ip_rep1.counts.txt ip_rep2.counts.txt sminput_rep1.counts.txt sminput_rep2.counts.txt \
    -o merged_counts.tsv

library(DEWSeq)

counts <- read.table('merged_counts.tsv', sep='\t', header=TRUE, row.names=1)
colData <- data.frame(
    sample = c('ip_rep1','ip_rep2','sminput_rep1','sminput_rep2'),
    type = c('ip','ip','sminput','sminput'),
    condition = c('treated','treated','untreated','untreated')
)

dds <- DESeqDataSetFromSlidingWindows(
    countData = counts,
    colData = colData,
    annotObj = 'annotation_windows.bed',
    design = ~ type + condition
)

dds <- DESeq(dds)
# In a model with `~ type + condition`, the IP-vs-input contrast is the simple condition
# main effect; to test the differential CLIP signal between conditions use the interaction
# term name from the design matrix (matches the skill's interaction-term guidance below):
res <- results(dds, name='typeip.conditiontreated')

# Window-level FDR adjustment
res_filtered <- res[!is.na(res$padj) & res$padj < 0.05 & abs(res$log2FoldChange) > 1, ]

# Aggregate adjacent significant windows into differential regions
sig_windows <- as.data.frame(res_filtered)
sig_windows$chr <- gsub('_.*', '', rownames(sig_windows))
# Custom reduce: combine adjacent windows within 100 nt

Flipper Workflow (Skipper-Coupled)

Flipper (Tu 2024) is the differential companion to Skipper, using the same 100 nt feature-respecting windows and the same beta-binomial framework.

# Assume Skipper has been run on all samples; Skipper output is at skipper_out/
flipper differential \
    -i skipper_out/ \
    --design design.tsv \
    --contrast treatment vs control \
    -o flipper_out/

# design.tsv format:
# sample_id   condition   replicate   ip_or_input
# ip_treat_r1 treatment   1           ip
# ip_treat_r2 treatment   2           ip
# in_treat_r1 treatment   1           input
# ...

Output: differential window BED with log2 FC, p, padj per window.

Peak-Level Differential (CLIPper Upstream)

library(DESeq2)
library(GenomicRanges)
library(Rsubread)

# Step 1: union of CLIPper peaks across conditions
# (See bedtools merge upstream)
peaks <- read.table('consensus_peaks.bed', sep='\t', col.names=c('chr','start','end','name','score','strand'))

# Step 2: count reads per peak per sample with featureCounts
saf <- data.frame(
    GeneID = peaks$name,
    Chr = peaks$chr,
    Start = peaks$start + 1,  # 1-based for featureCounts
    End = peaks$end,
    Strand = peaks$strand
)
counts_ip <- featureCounts(
    files = c('ip_treat_r1.bam','ip_treat_r2.bam','ip_ctrl_r1.bam','ip_ctrl_r2.bam'),
    annot.ext = saf,
    isGTFAnnotationFile = FALSE,
    strandSpecific = 1,
    isPairedEnd = TRUE
)$counts

counts_in <- featureCounts(
    files = c('in_treat_r1.bam','in_treat_r2.bam','in_ctrl_r1.bam','in_ctrl_r2.bam'),
    annot.ext = saf,
    isGTFAnnotationFile = FALSE,
    strandSpecific = 1,
    isPairedEnd = TRUE
)$counts

# Step 3: DESeq2 with IP vs SMInput interaction
all_counts <- cbind(counts_ip, counts_in)
colData <- data.frame(
    type = rep(c('ip','input'), each=4),
    condition = rep(c('treat','treat','ctrl','ctrl'), 2),
    replicate = rep(c('r1','r2','r1','r2'), 2)
)

dds <- DESeqDataSetFromMatrix(countData = all_counts, colData = colData,
                               design = ~ type + condition + type:condition)
dds <- DESeq(dds)

# The interaction term `typeip.conditiontreat` tests:
# does the IP/input ratio differ in treatment vs control?
res <- results(dds, name = 'typeip.conditiontreat')

# Filter
res_sig <- res[!is.na(res$padj) & res$padj < 0.05 & abs(res$log2FoldChange) > 1, ]

The interaction-term design (type:condition) is the correct statistical model for differential CLIP: it tests whether the IP-vs-input ratio differs across conditions, which is what "differential binding" means. Naive testing of just condition (ignoring SMInput) confounds binding changes with expression changes.

RBP Knockdown Experiment Design

The canonical differential CLIP design is to knock down the RBP and observe what binding sites are lost. Caveats:

| Issue | Implication | Mitigation | |-------|-------------|------------| | RBP KD also depletes the RBP protein in cells | siRNA/shRNA reduces RBP -> reduces IP yield -> reduces unique fragments | Normalize against SMInput WITHIN each condition; the relative IP/SMInput captures binding, not protein level | | RBP KD changes transcript abundance | mRNA stability regulators (HuR, PUM2) when knocked down change target abundance | Both IP and SMInput see the change; ratio still works | | Off-target effects of siRNA | Multiple binding profiles change | Use multiple independent siRNAs; require concordance | | KD efficiency varies | Lower KD -> smaller binding-loss signal | Validate KD by WB on the same IP lysate; > 70% protein loss target | | Rescue requires re-introducing RBP | siRNA-resistant RBP cDNA for rescue | The standard differential validation experiment |

Per-Tool Failure Modes

DEWSeq -- Slow on dense libraries

Trigger: Whole-genome window tiling at 20 nt step; dense library (50M unique fragments); 4+ samples.

Mechanism: DEWSeq runs DESeq2 internals on millions of windows; the dispersion fit on this many features is slow.

Symptom: Runtime > 6 h; out-of-memory; "size of object exceeds vector limit".

Fix: Increase window step size to 50 nt; pre-filter windows with low counts; or restrict to expressed transcripts only. DEWSeq vignette suggests keep <- rowSums(counts(dds)) >= 30; dds <- dds[keep,] before testing.

DEWSeq -- Custom adjacency aggregation needed

Trigger: Windows are 50 nt; biological binding sites are 50-500 nt; user expects DEWSeq to output continuous "differential regions" but gets individual windows.

Mechanism: DEWSeq outputs per-window results; aggregating adjacent significant windows into regions is a separate step.

Symptom: Output has 10,000 individual windows; user expects 1,000 biological regions.

Fix: Use the DEWSeq utility resultsDEWSeq() then bedtools merge -d 100 on the significant-window BED. Or use the top_hits_to_bed.R script from DEWSeq examples.

Peak-level differential -- Peak set differs between conditions

Trigger: CLIPper called peaks separately per condition; treatment has peaks at sites missing in control (and vice versa).

Mechanism: Peak unification requires a consensus peakset; testing on a "treatment-only" peak underestimates evidence in control (zero reads) and produces spurious DE.

Symptom: "Treatment-specific" peaks dominate DE results; biologically implausible.

Fix: Generate consensus peakset across all conditions (bedtools merge of all per-condition peak BEDs); count reads per consensus peak across all samples; THEN run differential. The Yeo lab convention is consensus peakset across all samples.

Interaction term forgotten

Trigger: DESeq2 design ~ condition instead of ~ type + condition + type:condition.

Mechanism: Simple ~ condition tests whether read counts differ between treatment and control regardless of whether reads are from IP or SMInput. A condition-driven expression change in SMInput is detected as DE binding.

Symptom: DE results dominated by transcripts with global expression changes (housekeeping shifts).

Fix: Always use the interaction-term design. The biologically meaningful test is the interaction type:condition p-value.

Normalization assumptions

Trigger: edgeR calcNormFactors(method='TMM') on CLIP-seq data.

Mechanism: TMM assumes most features (genes) are not differentially expressed. CLIP-seq peak counts can be globally shifted if the RBP itself is knocked down; TMM normalization would force the shift to be invisible.

Symptom: Knockdown experiment shows ~0 DE peaks; expected hundreds.

Fix: Use SMInput as the control library; spike-in normalization if available; or skip TMM and use library-size normalization only. Some CLIP-specific tools (DEWSeq, Flipper) handle this internally.

Flipper requires Skipper upstream

Trigger: Flipper called on CLIPper output.

Mechanism: Flipper expects Skipper's window-level output format with beta-binomial estimates.

Symptom: Flipper crashes or produces nonsense.

Fix: Use DEWSeq with htseq-clip for CLIPper-upstream workflows; use Flipper only with Skipper.

Decision Tree by Scenario

| Scenario | Tool + design | Why | |----------|---------------|-----| | KD vs control eCLIP, CLIPper upstream | DEWSeq + htseq-clip | CLIP-specific NB GLM with interaction term | | KD vs control eCLIP, Skipper upstream | Flipper | Skipper companion; matches windowing | | Treatment vs vehicle (small effect) | DEWSeq (window-level higher power) | Sliding windows capture small shifts | | Multiple time points | DEWSeq with time as covariate | Continuous design with time vector | | Allele-specific differential | BEAPR per-allele + custom logistic | See clip-seq/clip-alignment for WASP | | Single-cell CLIP differential | MASTR-seq or scCLAM | Specialized; few options | | Differential motif occupancy | Window-level + DEWSeq + motif overlap | Combine differential windows with motif BED | | RBP overexpression vs control | Same as KD reversed | Same statistical framework | | Compare two RBPs | NOT differential CLIP; use SPIDR or separate CLIPs | Different RBPs need separate IPs | | Spike-in normalization needed | DEWSeq + spike-in size factors | For global occupancy shifts | | chimeric eCLIP differential miRNA targets | Custom; treat each miRNA-target chimera as feature | Specialized; see clip-seq/ago-clip-mirna-targets |

Reconciliation: When Differential Tools Disagree

| Pattern | Likely cause | Action | |---------|--------------|--------| | DEWSeq finds many DE windows; edgeR peak-level finds few | Window-level higher power for narrow shifts | Aggregate DEWSeq windows; cross-check | | edgeR many DE; DEWSeq few | edgeR not accounting for SMInput | Re-run edgeR with interaction term | | DE peaks dominated by expression changes | No interaction term | Use ~ type + condition + type:condition | | KD experiment shows ~0 DE | TMM over-corrects global shift | Switch to library-size norm only; use SMInput | | siRNA replicates discordant | Off-target effects vary | Use multiple independent siRNAs; require concordance | | Treatment-only peaks dominate DE | No consensus peakset | Generate consensus first; then test on unified set | | Significant windows scattered | Window aggregation step skipped | bedtools merge -d 100 on significant-window BED | | Same gene appears in many DE windows | Multiple binding sites per gene differential | Report at gene-level too; not just window-level | | Flipper fails with non-Skipper input | Upstream mismatch | Use DEWSeq for non-Skipper workflows | | DESeq2 dispersion fit fails | Too few replicates (n=2 per condition); too few features after filtering | Increase replicates; or relax filtering |

Operational rule for high-confidence differential reporting: (a) Use SMInput-aware design (~ type + condition + type:condition); (b) generate consensus peakset across conditions; (c) require padj < 0.05 AND |log2FC| > 1; (d) require concordance with at least one orthogonal method (e.g., DEWSeq + edgeR peak-level on same data); (e) for KD experiments, validate KD efficiency by WB and require multiple independent siRNAs.

Common Errors

| Error / symptom | Cause | Solution | |-----------------|-------|----------| | DESeq2 ~ condition finds many DE | No interaction with type | Use ~ type + condition + type:condition and test interaction | | DEWSeq output is gene-level not window-level | keep <- ... filter too aggressive | Loosen pre-filter | | Adjacent significant windows not aggregated | Forgot bedtools merge | bedtools merge -d 100 on sig windows | | edgeR TMM fits all libraries to one value | Most-features-not-DE assumption violated | Use library-size norm; or DEWSeq for CLIP-specific | | Flipper crashes on CLIPper input | Tool mismatch | Switch to DEWSeq | | Few replicates -> unstable estimates | n=2 not enough for dispersion | Increase n; or use limma-voom (more tolerant) | | Peaks differ across conditions | Per-condition peak calls | Unify with consensus peakset | | KD experiment yields 0 DE | Normalization over-corrected; OR KD efficiency too low | Validate KD WB; check normalization | | Global shift in IP relative to SMInput | RBP itself knocked down so IP yield lower | Normalize WITHIN each condition | | Lots of "treatment-only" peaks | Caller stringency higher in one condition | Use consensus peakset for fairness |

References

• Schwarzl T et al 2020 Bioconductor pkg DEWSeq (windowed CLIP differential)
• Sahadevan S et al 2022 Methods Mol Biol 2404:189 (DEWSeq + htseq-clip pipeline)
• Tu Y et al 2024 PMC PMC13060198 (Flipper, Skipper-companion differential)
• Boyle EA et al 2023 Cell Genomics 3:100317 (Skipper, parent of Flipper)
• Love MI et al 2014 Genome Biol 15:550 (DESeq2)
• McCarthy DJ et al 2012 Nucleic Acids Res 40:4288 (edgeR)
• Ritchie ME et al 2015 Nucleic Acids Res 43:e47 (limma)
• Wu B et al 2018 Nat Commun 9:5117 (BEAPR allele-specific differential)
• Van Nostrand EL et al 2020 Nature 583:711 (ENCODE 150 RBP shRNA + eCLIP comparison)

Related Skills

• clip-seq/clip-peak-calling - CLIPper / Skipper outputs feed differential
• clip-seq/clip-qc - Replicate QC required for valid differential
• clip-seq/binding-site-annotation - Annotate differential regions
• clip-seq/clip-motif-analysis - Motif analysis on differential windows
• clip-seq/ago-clip-mirna-targets - Differential miRNA targeting from chimeric eCLIP
• differential-expression/deseq2-basics - Underlying NB GLM model
• differential-expression/de-results - DE results interpretation
• differential-expression/edger-basics - edgeR for peak counts
• chip-seq/differential-binding - DNA-protein analogue