Bio Clip Seq Clip Peak Calling

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

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.

ENCODE eCLIP Stringency Thresholds

The Van Nostrand et al. ENCODE eCLIP standards define two stringency levels for CLIPper peaks normalized against SMInput:

Both thresholds are applied to the IDR-passing reproducible peak set. ENCODE requires:

• = 2 biological replicates

• = 1M unique fragments per replicate (or saturated peak detection)

• IDR rescue ratio and self-consistency ratio both < 2
• Narrow-binding RBPs: FRiP >= 0.005 (atypical-binding RBPs exempt)

# CLIPper -> SMInput normalization -> ENCODE thresholds
clipper -b dedup.bam -s hg38 -o peaks.bed --save-pickle --FDR-alpha 0.05 --superlocal

# Normalize against SMInput; eclip-norm is the Yeo lab tool
python overlap_peakfi_with_bam_PE.py peaks.bed dedup.bam sminput.bam dedup.bam.readnum.txt sminput.bam.readnum.txt peaks.normed.bed
python compress_l2foldenrpeakfi_for_replicate_overlapping_bedformat.py peaks.normed.bed peaks.compressed.bed

# Stringent filter: log2 FC >= 3, -log10 p >= 3
awk 'BEGIN{FS=OFS="\t"} $4 >= 3 && $5 >= 3' peaks.compressed.bed > peaks.stringent.bed

# Lenient
awk 'BEGIN{FS=OFS="\t"} $4 >= 1 && $5 >= 2' peaks.compressed.bed > peaks.lenient.bed

The Yeo lab scripts (overlap_peakfi_with_bam_PE.py, compress_l2foldenrpeakfi_for_replicate_overlapping_bedformat.py) live in the eclip-pipeline repository; they implement the exact log2 fold-change computation used for ENCODE data downloads.

Per-Caller Failure Modes

CLIPper -- High-expression transcript bias

Trigger: RBP that binds rare transcripts (snoRNAs, tRNAs, mitochondrial mRNAs) such as TROVE2 (Y RNA), NSUN2 (tRNA), or FASTKD2 (chrM).

Mechanism: CLIPper's Poisson local-lambda model assumes only a minority of reads come from binding. For rare-transcript binders, the IP reads ARE the majority signal on that transcript; the model treats them as background.

Symptom: Per-transcript peak count plummets compared to qualitatively-visible IP signal; chrM peaks called by Skipper but missed by CLIPper for FASTKD2.

Fix: Use Skipper for rare-transcript binders (beta-binomial against SMInput respects the actual IP/IN ratio independently of local lambda). Or pre-restrict CLIPper to the target transcript with --gene custom.bed.

CLIPper -- Requires SMInput downstream

Trigger: Submitting CLIPper output as "peaks" without IP/SMInput log2 normalization.

Mechanism: CLIPper's own p-value is uncorrected for SMInput background. ENCODE peak BED files distributed on the portal are CLIPper + SMInput log2 normalization combined.

Symptom: Peak count from the analysis is 5-10x higher than ENCODE values for the same RBP; the peak set includes obvious housekeeping-gene noise.

Fix: Always normalize CLIPper output against SMInput. The Yeo lab scripts (overlap_peakfi_with_bam_PE.py + compress_l2foldenrpeakfi_for_replicate_overlapping_bedformat.py) produce the canonical ENCODE-style peak.

PureCLIP -- Too focal; misses broad binding zones

Trigger: RBP with broad binding mode (PUM2 in 3' UTR clusters; SR proteins across exonic enhancer regions); using PureCLIP and disappointed.

Mechanism: PureCLIP's HMM emits single-nt crosslink-site state at high stringency. In the Boyle 2023 Skipper benchmark, PureCLIP was highly focal with low recall on broad-binding RBPs (few true-positive windows).

Symptom: Site count 100x lower than expected; sites cluster within known binding regions but vast majority of region is "background" in PureCLIP output.

Fix: Use PureCLIP for single-nt CL maps (motif registration, ASB) but pair with CLIPper or Skipper for the broader binding-site list. PureCLIP's -or regions file gives broader output but is still focal compared to CLIPper.

Piranha -- Top-expression-decile bias

Trigger: Highly expressed transcript (GAPDH, ACTB, rRNA-flanking) appears in most-significant Piranha peaks even when RBP biology is elsewhere.

Mechanism: ZTNB regression weights peak significance by absolute count, not relative enrichment. Genes with 10x baseline expression dominate the top of the peak list.

Symptom: Top-100 Piranha peaks include housekeeping genes; motif enrichment in top peaks is weak; functional GO term analysis of peak-bearing genes returns generic "cellular metabolism."

Fix: Use Piranha with SMInput as a covariate (-c sminput_counts.bed -l) so the regression accounts for input pile-up. Or switch to Skipper for principled input normalization.

Piranha -- Convergence failure with covariates

Trigger: -c rna_seq.bed for expression normalization; large per-bin RNA-seq counts (> 10000).

Mechanism: ZTNB optimization fails on extreme covariate values; the iterative algorithm diverges.

Symptom: "Maximum iterations reached" or "Singularity in regression" error; or output BED file empty.

Fix: Pass -l to convert covariates to log-space; or pre-filter covariate BED to remove top-1% outlier bins; or switch to Skipper which handles input normalization internally.

omniCLIP -- Mitochondrial transcript blind spot

Trigger: Mitochondrial RBP (FASTKD2, LRPPRC) submitted to omniCLIP.

Mechanism: omniCLIP's HMM is trained on autosomal/X transcripts; chrM has different background statistics that the model treats as outlier (or excludes if blacklisted).

Symptom: chrM peaks return 0; FASTKD2 peak file has no peaks at all.

Fix: Use Skipper for chrM-binding RBPs (GC-stratified beta-binomial is per-bin and indifferent to chromosome). Or call peaks on chrM separately with PureCLIP.

CTK CIMS -- HITS-CLIP only; iCLIP/eCLIP misuse

Trigger: CTK CIMS (-mutation flag) applied to iCLIP/eCLIP data.

Mechanism: iCLIP/eCLIP use truncation (RT stops at adduct), not mutations. CIMS expects deletions/substitutions inside the read; CITS is the truncation-based companion.

Symptom: CIMS returns very few sites (mutation rate of CL is ~3-7% for HITS-CLIP, but only ~1-2% for iCLIP since most reads truncate before reaching the adduct).

Fix: For iCLIP/eCLIP, use CTK CITS (tag2cluster.pl ... -cs5 5 truncation-mode). For HITS-CLIP and PAR-CLIP, CIMS is correct (HITS deletions, PAR T->C).

CLAM peakcaller -- Multi-mapper BAM absent

Trigger: Ran STAR with --outFilterMultimapNmax 1 then tried CLAM; CLAM fails or returns identical output to Piranha.

Mechanism: CLAM's EM rescue requires the multi-mapper BAM from STAR (--outFilterMultimapNmax 100 --outSAMmultNmax -1). Without it, CLAM falls back to unique-only Piranha-style peak calling.

Symptom: CLAM output looks identical to Piranha; repeat peaks not rescued.

Fix: Re-align with multi-mapper-aware STAR parameters (see clip-seq/clip-alignment). Verify the multi-mapper BAM has alignments at non-unique positions before running CLAM.