bio-clip-seq-clip-peak-calling

What this skill does

## Version Compatibility

Reference examples tested with: CLIPper 2.0+, PureCLIP 1.3.1+, Piranha 1.2.1+, omniCLIP 0.2.0+, CTK 1.1.4+, CLAM 1.2+, Paraclu 9+, Skipper (commit 2023.05+), MACS3 3.0+, bedtools 2.31+, samtools 1.19+, idr 2.0.4+.

Before using code patterns, verify installed versions match. If versions differ:
- 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 binary and adapt the example to match the actual CLI rather than retrying. PureCLIP 2.x changed several flag names; Skipper is distributed as a Snakemake workflow with frequently-evolving paths.

# CLIP-seq Peak Calling

**"Call protein-RNA binding sites from my deduplicated CLIP BAM"** -> Identify regions where read pile-up (and, for iCLIP/eCLIP, single-nucleotide truncations) exceed background from a size-matched input (SMInput) control. The choice of peak caller depends on (a) the CLIP variant (HITS-CLIP, iCLIP, eCLIP, PAR-CLIP), (b) whether SMInput is available, (c) the RBP binding mode (narrow motif vs broad zones vs repeat-binding), and (d) the goal (publication-comparable ENCODE peaks vs single-nucleotide crosslink sites vs high-recall site lists).

- CLI (ENCODE eCLIP canonical): `clipper -b dedup.bam -s hg38 -o peaks.bed --save-pickle --FDR-alpha 0.05 --superlocal` then `eclip-norm peaks.bed -i sminput.bam` for log2 fold-change against SMInput (`--FDR-alpha` is the CLIPper flag for the FDR cutoff; older docs sometimes shorten to `--FDR`)
- CLI (single-nucleotide crosslink sites, iCLIP/eCLIP): `pureclip -i dedup.bam -bai dedup.bam.bai -g genome.fa -ibam sminput.bam -ibai sminput.bam.bai -o sites.bed -or regions.bed -nt 8 -dm 8`
- CLI (high-recall, beta-binomial windowed; needs SMInput): `Skipper` Snakemake workflow with config matching cell type and SMInput BAM
- CLI (no SMInput, no truncation): `Piranha -b 50 -p 0.01 -d ZeroTruncatedNegativeBinomial -s -o peaks.bed dedup.bam` (Piranha takes the BAM as a positional argument, not via `-s`)
- CLI (CIMS/CITS single-nt from CTK): `tag2cluster.pl dedup.bed cluster.bed --multi-tag-method coverage`; then `bedExtractCIMS.pl cluster.bed cims.bed`
- CLI (multi-mapper rescue for repeat-binding RBPs): `CLAM peakcaller -i unique.bam multimap.bam -o clam_out_dir/ --gtf gencode.gtf` (CLAM peakcaller writes peaks into an OUTPUT DIRECTORY, not a single file)

The ENCODE eCLIP gold standard for "stringent" peaks is: log2(IP / SMInput) >= 3 AND -log10(p-value) >= 3. "Lenient" peaks use log2 >= 1 AND -log10 >= 2. Both filters operate on CLIPper peak output normalized against SMInput. Without SMInput, neither stringency level can be reproduced - the algorithm has no background term.

## Algorithmic Taxonomy

| Caller | Model | Resolution | SMInput required | Strength | Fails when |
|--------|-------|------------|------------------|----------|------------|
| CLIPper (Yeo) | Poisson with 500 bp local lambda, cubic-spline interpolation | Peak (10-500 nt) | Not for calling; required for ENCODE normalization | ENCODE eCLIP standard; reproducible against published ENCODE peak sets | Assumes most reads not from binding; fails when IP signal dominates a gene; sensitive to highly expressed transcripts |
| PureCLIP (Krakau 2017) | Non-homogeneous HMM jointly modeling enrichment + truncation + sequence biases | Single-nucleotide CL + broad region | Optional via -ibam; recommended | Single-nt resolution; only caller that explicitly models the iCLIP truncation pattern | Misses broad binding zones; very focal (Skipper 2023 benchmark: only ~4 sites per CLIP on test set; F1 ~0.2) |
| Skipper (Boyle 2023) | GC-stratified beta-binomial, 100 bp feature-respecting windows | Window (~100 nt) | Mandatory | 210-320% more sites than CLIPper; 8x faster; properly normalizes against input | Loses single-nt resolution; relatively new (2023); fewer published peak sets to compare |
| Piranha (Smith lab) | Zero-truncated negative binomial regression with optional covariates | Bin (50-200 nt) | Optional; pass as covariate | Mature, widely cited, handles count overdispersion | Biased toward high-expression transcripts; convergence fails with large covariate values (use `-l` log-space) |
| omniCLIP (Drewe-Boss 2018) | Non-homogeneous HMM with Dirichlet-multinomial of variants | Peak | Required | Models replicate variance; integrates background; can call peaks on any CLIP variant | Slow on deep libraries; F1 ~0.4 on standard benchmarks; fails for mitochondrial transcripts (Skipper benchmark: 0 chrM windows for FASTKD2) |
| CTK CIMS/CITS (Shah 2017) | Empirical FDR on crosslink-induced mutations or truncations | Single-nucleotide | No (uses background mutation rate from non-bound transcripts) | Single-nt resolution; works on HITS-CLIP deletions, PAR-CLIP T->C, iCLIP truncations | Empirical FDR less principled than HMM/beta-binomial; perl-based pipeline harder to integrate |
| CLAM peakcaller (Zhang & Xing 2017) | EM-assigned multi-mapper count + Piranha-like negative binomial | Peak | Optional | Only solution for repeat-binding RBPs; 10-30% additional sites in repeats | Inherits Piranha limitations on coverage-based stats; slower than Piranha |
| Paraclu (Kawaji) | Parametric clustering with min-density and max-density thresholds | Variable | No | Simple, parameter-tunable, works on bedgraph | Heuristic; no statistical significance; cluster boundaries sensitive to thresholds |
| PIPE-CLIP | Online CLIP pipeline with custom peak caller | Peak | Optional | Web interface; integrated end-to-end | Slow on cloud; less customizable; community has moved on |
| CLIPick (Park 2018) | Expression-deconvolved peak caller | Peak | No (RNA-seq used instead) | Models RNA-seq abundance as background | RNA-seq cannot capture nonspecific binding; less popular post-Skipper |
| Flipper (Tu 2024) | Negative-binomial differential test downstream of Skipper | Window | Yes | Companion differential tool to Skipper | Only meaningful in differential context; see clip-seq/differential-clip |
| MACS3 callpeak | Local Poisson | Peak | Optional (treats SMInput as ChIP-seq input) | Familiar, fast | Not designed for CLIP; misses truncation signal; produces wider-than-typical peaks |

Methodology evolves; verify the current ENCODE eCLIP standard operating procedure (encodeproject.org/eclip) and the nf-core/clipseq pipeline configuration before locking on a single caller. The 2023 Skipper benchmark (Boyle Cell Genomics) is the most recent comprehensive comparison and is the rationale for many recent eCLIP reanalyses.

## Critical Choice: Coverage-Based vs Crosslink-Site-Based vs Window-Based

Three fundamentally different statistical frameworks:

**Coverage-based (CLIPper, Piranha, MACS3, Paraclu):** Tally reads in a window; significance from local Poisson or negative binomial. Produces multi-nt peaks. Best when binding zones are broader than the read footprint (PUM2 3' UTR, HuR ARE elements).

**Crosslink-site-based (PureCLIP, CTK CIMS/CITS, PARalyzer):** Identify single nucleotides enriched in truncations (iCLIP/eCLIP) or specific mutations (PAR-CLIP T->C, HITS-CLIP deletions). Produces single-nt sites that can be aggregated into footprints. Best when single-nt resolution is essential (motif registration with mCross; allele-specific binding).

**Window-based (Skipper):** Tile transcriptome a priori into fixed-size feature-respecting windows; test each window's IP/IN ratio with beta-binomial. Produces window-level calls. Best when SMInput is available and the goal is high-recall, comparable-across-RBPs site sets.

| Goal | Recommended caller | Why |
|------|-------------------|-----|
| Reproduce ENCODE published peaks | CLIPper + IP/SMInput normalization | The Yeo lab pipeline is the canonical reference |
| Maximum sensitivity, modern statistics, SMInput available | Skipper | 210-320% more sites; explicit input normalization |
| Single-nucleotide crosslink-site ma