End-to-end CLIP-seq pipeline from FASTQ to ENCODE-compliant binding sites, single-nucleotide crosslink maps, annotation, motifs, and (optionally) differential binding. Use when running the full Yeo lab eCLIP / iCLIP / iCLIP2 / iCLIP3 / irCLIP / PAR-CLIP analysis with SMInput control, protocol-specific UMI extraction, ENCODE STAR parameters, CLIPper or Skipper peak calling with stringent log2 FC and -log10 p thresholds, IDR rescue and self-consistency QC, and downstream motif registration with mCross or PEKA.
End-to-end CLIP-seq pipeline from FASTQ to ENCODE-compliant binding sites, single-nucleotide crosslink maps, annotation, motifs, and (optionally) differential binding. Use when running the full Yeo lab eCLIP / iCLIP / iCLIP2 / iCLIP3 / irCLIP / PAR-CLIP analysis with SMInput control, protocol-specific UMI extraction, ENCODE STAR parameters, CLIPper or Skipper peak calling with stringent log2 FC and -log10 p thresholds, IDR rescue and self-consistency QC, and downstream motif registration with mCross or PEKA.
Before using code patterns, verify installed versions match. If versions differ:
CLI: <tool> --version then <tool> --help to confirm flags
Python: pip show <package> then help(module.function) to check signatures
R: packageVersion('<pkg>') then ?function_name to verify parameters
If code throws unexpected errors, introspect the installed tool and adapt the example rather than retrying.
CLIP-seq End-to-End Pipeline
"Analyze my CLIP-seq data from raw FASTQ to ENCODE-compliant binding sites" -> Orchestrate protocol-specific UMI extraction, 3'-only adapter trimming (preserving the R2 5' truncation = crosslink site -1), ENCODE STAR alignment, UMI-based deduplication, library complexity QC, peak calling against SMInput with stringent thresholds (log2 FC >= 3 AND -log10 p >= 3), single-nucleotide crosslink-site detection, ChIPseeker annotation with CLIP-appropriate tssRegion, motif discovery with GC-matched background and CL-position registration, and optional differential binding between conditions.
This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step.
The governing principle
A CLIP callset is decided at four seams, not inside the peak caller.
The protocol variant is the master commitment made once at the top. It fixes the UMI pattern, the STAR mismatch ceiling, and the crosslink signal (truncation vs PAR-CLIP T->C vs STAMP C->U edit). The CLIP Variant Selection table below is that decision; everything downstream inherits it.
Every IP is normalized against its SMInput with ENCODE-stringent thresholds (log2 FC >= 3 AND -log10 p >= 3). Peaks called without SMInput normalization are enrichment-uncontrolled and dominated by abundance.
The R2 5' end IS the crosslink site (-1), so preprocessing and alignment must PRESERVE it. Trim 3'-only and permissively (-q 6, never -g on R1), and align with STAR --alignEndsType EndToEnd - soft-clipping or aggressive 5' trimming destroys the truncation base and with it single-nucleotide resolution.
40-70% PCR duplication is BY DESIGN; low duplication signals a FAILED IP, not a clean library. The IP enriches a small molecule pool, so the real quality metric is the UNIQUE-fragment count after UMI dedup - never the raw duplication rate.
Pipeline Overview
FASTQ + SMInput
-> [clip-preprocessing] UMI extract + 3' adapter trim (-q 6 -m 18) + two-pass for eCLIP
-> [clip-alignment] STAR ENCODE block (alignEndsType EndToEnd, mismatch 0.04 or 0.07 for PAR-CLIP) + UMI dedup
-> [clip-qc] preseq, FRiP, IDR rescue + self-consistency, read distribution
-> [clip-peak-calling] CLIPper + SMInput log2 norm (stringent: log2 FC >= 3, -log10 p >= 3) OR Skipper (substantially more sites)
-> [crosslink-site-detection] PureCLIP or CTK CITS for single-nt CL positions
-> [binding-site-annotation] ChIPseeker (tssRegion=c(-100,100), level=transcript) + RBP-Maps for splicing factors
-> [clip-motif-analysis] HOMER + mCross (registered) + RBNS Kd cross-check
-> [differential-clip] DEWSeq window-level NB with type:condition interaction (optional)
CLIP Variant Selection
Variant
When to use
UMI pattern
STAR mismatch ceiling
Detection signal
eCLIP (Van Nostrand 2016)
ENCODE comparability; SMInput available
10 nt R1
0.04
R2 5' truncation
iCLIP / iCLIP2 / iCLIP3
Single-end; high motif specificity
NNNXXXXNN (3+4+2; demux first)
0.04
R1 5' truncation
irCLIP / FLASH
Non-radioactive; fast
Protocol-specific
0.04
Truncation
PAR-CLIP
Photoactivatable nucleoside (4SU); HEK293/K562
4 nt typical
0.07 (raised for T->C)
T->C transitions
miCLIP / miCLIP2
m6A modification
iCLIP-style
0.04
Truncation + C->T at m6A
STAMP / scSTAMP
Antibody-free; in vivo or single-cell
NA (no UV)
0.04 (RNA-seq mode)
C->U editing (RBP-APOBEC1 fusion)
chimeric eCLIP / miR-eCLIP
Direct miRNA-target pairs
10 nt R1
0.04
Chimeric reads
Step 1: Quality Control of Raw FASTQ
# Initial QC
fastqc raw_R1.fq.gz raw_R2.fq.gz -o qc/raw/
# Inspect first 12 bases of 100 reads to verify UMI pattern matches the prep
zcat raw_R1.fq.gz | awk 'NR%4==2' | head -100 | cut -c1-12 | sort | uniq -c | sort -rn | head# Random barcode positions show ~25% per base; library barcodes are fixed
Step 2: Preprocessing (Protocol-Specific)
Goal: Convert raw CLIP FASTQ into UMI-deduplicated, alignment-ready FASTQ while preserving the R2 5' end (= crosslink site -1) that drives single-nucleotide resolution downstream.
Approach: Use the protocol-matched UMI pattern (10 nt eCLIP, NNNXXXXNN iCLIP, 4 nt PAR-CLIP), run umi_tools extract to move random barcodes to read names, then apply cutadapt with 3'-only adapter trimming at -q 6 -m 18 (permissive 5' to protect the truncation base). eCLIP uses two-pass trimming to remove read-through inline adapters from R2 5' only; iCLIP and PAR-CLIP use single-pass.
# eCLIP: 10 nt UMI on R1; two-pass adapter trim for read-through# See clip-seq/clip-preprocessing for protocol-specific patterns
umi_tools extract \
--bc-pattern=NNNNNNNNNN \
--stdin=raw_R1.fq.gz --read2-in=raw_R2.fq.gz \
--stdout=R1.umi.fq.gz --read2-out=R2.umi.fq.gz \
--log=qc/umi_extract.log
# Pass 1: 3' adapter on both reads# -q 6 is intentionally permissive; aggressive trimming destroys R2 5' = CL site -1
cutadapt \
-a AGATCGGAAGAGCACACGTCT \
-A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \
--quality-base 33 -q 6 -m 18 \
-j 8 \
-o R1.p1.fq.gz -p R2.p1.fq.gz \
R1.umi.fq.gz R2.umi.fq.gz \
> qc/cutadapt_pass1.log 2>&1
# Pass 2: strip read-through 5' adapter from R2 only (NEVER -g on R1)
cutadapt \
-G GATCGTCGGACTGTAGAACTCTGAAC \
--quality-base 33 -q 6 -m 18 \
-j 8 \
-o R1.trim.fq.gz -p R2.trim.fq.gz \
R1.p1.fq.gz R2.p1.fq.gz \
>> qc/cutadapt_pass2.log 2>&1
For PAR-CLIP: same UMI extraction but downstream alignment raises --outFilterMismatchNoverReadLmax from 0.04 to 0.07 (the T->C signature would otherwise be filtered as sequencing error). See clip-seq/clip-preprocessing for full per-protocol guidance.
For PAR-CLIP: change --outFilterMismatchNoverReadLmax 0.04 to 0.07. For repeat-binding RBPs (MATR3, ZFP36, FUS at LINE-1, HNRNPK at SINEs): change --outFilterMultimapNmax 1 to 100 and add --outSAMmultNmax -1, then run CLAM downstream for EM-based multi-mapper assignment. See clip-seq/clip-alignment for full guidance.
CLIP libraries have 40-70% PCR duplication BY DESIGN (the IP enriches a small molecule pool). Low duplication usually means failed IP, not a good library. The unique-fragment count after UMI dedup is the actual quality metric. See clip-seq/clip-qc for full five-gate diagnostic.
Step 5: Peak Calling
# CLIPper (ENCODE canonical) + SMInput log2 normalization
clipper \
-b sample_dedup.bam \
-s GRCh38 \
-o peaks/sample.clipper.bed \
--FDR 0.05 \
--save-pickle \
--processors 8 # super-local p-values are hard-coded ON in current CLIPper; the --superlocal flag was removed# ENCODE stringent: log2(IP/SMInput) >= 3 AND -log10 p >= 3# (Yeo lab eclip-pipeline scripts implement the normalization; see clip-seq/clip-peak-calling)
python overlap_peakfi_with_bam_PE.py \
peaks/sample.clipper.bed \
sample_dedup.bam sminput_dedup.bam \
sample_dedup.bam.readnum.txt sminput_dedup.bam.readnum.txt \
peaks/sample.normed.bed
python compress_l2foldenrpeakfi_for_replicate_overlapping_bedformat.py \
peaks/sample.normed.bed \
peaks/sample.compressed.bed
# Stringent filter
awk 'BEGIN{FS=OFS="\t"} $5 >= 3 && $6 >= 3' peaks/sample.compressed.bed > peaks/sample.stringent.bed
For maximum sensitivity (substantially more sites than CLIPper for mRNA-binding RBPs), use the Skipper Snakemake workflow with the same SMInput control. Mandatory for FASTKD2 / mt-RBPs which CLIPper misses on chrM. See clip-seq/clip-peak-calling for the full caller taxonomy.
# PureCLIP: HMM jointly modeling enrichment + truncation + CL motif.# -iv learns HMM parameters on a CHROMOSOME SUBSET (semicolon-delimited) to cut memory/runtime# (per PureCLIP docs); it is NOT a BED. To limit the callset to expressed regions, pre-filter the input BAM.
pureclip \
-i sample_dedup.bam -bai sample_dedup.bam.bai \
-g genome.fa \
-ibam sminput_dedup.bam -ibai sminput_dedup.bam.bai \
-o crosslinks/sample.sites.bed \
-or crosslinks/sample.regions.bed \
-nt 8 -dm 8 \
-iv 'chr1;chr2;chr3;'
Single-nt CL sites feed mCross motif registration and allele-specific binding analyses. They are NOT a replacement for the broad peak list; complementary outputs. See clip-seq/crosslink-site-detection.
Step 7: IDR Across Replicates
# Sort each replicate's compressed BED by signal (log2 FC, column 5)sort -k5,5gr peaks/rep1.compressed.bed > peaks/rep1.sorted.bed
sort -k5,5gr peaks/rep2.compressed.bed > peaks/rep2.sorted.bed
# True replicates threshold 0.05
idr --samples peaks/rep1.sorted.bed peaks/rep2.sorted.bed \
--input-file-type bed --rank 5 \
--output-file qc/idr.true.out \
--idr-threshold 0.05 \
--plot --log-output-file qc/idr.log
# ENCODE rule: rescue + self-consistency ratios both < 2 to pass# Pseudo-replicate IDR (split BAM in half) at threshold 0.10
Default ChIPseeker tssRegion=c(-3000, 3000) over-extends for CLIP (would label 30-50% peaks as "Promoter"). Splicing factors additionally need RBP-Maps (Yeo lab) for the 1400 nt cassette-exon regulatory metagene. See clip-seq/binding-site-annotation.
Step 9: Motif Analysis (De Novo + CL-Registered)
# Extract peak sequences (strand-preserving)
bedtools getfasta -fi genome.fa -bed peaks/sample.stringent.bed -s -fo motifs/peaks.fa
# GC-matched 3' UTR background (NOT auto-shuffled, which biases to AU)
bedtools shuffle -i peaks/sample.stringent.bed -g chrom.sizes \
-incl expressed_3utr.bed -seed 42 > motifs/background.bed
bedtools getfasta -fi genome.fa -bed motifs/background.bed -s -fo motifs/background.fa
# HOMER de novo
findMotifs.pl motifs/peaks.fa fasta motifs/homer \
-rna -len 5,6,7,8 -p 8 -fasta motifs/background.fa
# mCross for CL-position-registered motif. mCross.pl takes a POSITIONAL FASTA of sequences# pre-extracted/registered around the CL sites and an output stem (not a BED + genome + -i/-g/-k/-o):# bedtools slop -i crosslinks/sample.sites.bed -g genome.sizes -b 10 | bedtools getfasta -fi genome.fa -bed - -s > motifs/peakseqs.fa
mCross.pl motifs/peakseqs.fa motifs/mcross # see clip-seq/clip-motif-analysis for options
UV254 crosslinking has a strong U bias at CL sites; naive logos centered on CL positions are U-enriched even for non-U-binding RBPs. mCross corrects this by registering motif relative to the CL offset. See clip-seq/clip-motif-analysis.
Step 10: Differential Binding (Optional, Across Conditions)
# DEWSeq window-level NB with the interaction-term design# The interaction `~ type + condition + type:condition` tests whether IP/SMInput ratio shifts;# naive `~ condition` confounds binding with expression changes.
library(DEWSeq)
counts <- read.table('counts/merged.tsv', sep='\t', header=TRUE, row.names=1)
colData <- data.frame(
type = relevel(factor(c('ip','ip','ip','ip','sminput','sminput','sminput','sminput')), ref='sminput'),
condition = relevel(factor(c('treat','treat','ctrl','ctrl','treat','treat','ctrl','ctrl')), ref='ctrl'))
dds <- DESeqDataSetFromSlidingWindows(
countData=counts, colData=colData,
annotObj='annotation.txt',# htseq-clip TAB annotation table (named columns), NOT a plain BED
design =~ type + condition + type:condition
)
dds <- DESeq(dds)# with sminput/ctrl as the references, the interaction coefficient is typeip.conditiontreat
res <- results(dds, name='typeip.conditiontreat')
See clip-seq/differential-clip for full DEWSeq workflow and the htseq-clip preprocessing required upstream.
Default tssRegion=c(-3000,3000) over-extends for CLIP
Tight tssRegion=c(-100,100), level='transcript'
Differential binding confounded with expression
~ condition design
~ type + condition + type:condition interaction (DEWSeq)
References
Van Nostrand EL, Pratt GA, Shishkin AA, et al (2016) Robust transcriptome-wide discovery of RNA-binding protein binding sites with enhanced CLIP (eCLIP). Nature Methods 13:508-514. DOI 10.1038/nmeth.3810.
Van Nostrand EL, Freese P, Pratt GA, et al (2020) A large-scale binding and functional map of human RNA-binding proteins. Nature 583:711-719. DOI 10.1038/s41586-020-2077-3. (ENCODE RBP; SMInput + IDR practice.)
Krakau S, Richard H, Marsico A (2017) PureCLIP: capturing target-specific protein-RNA interaction footprints from single-nucleotide CLIP-seq data. Genome Biology 18:240. DOI 10.1186/s13059-017-1364-2.
Li Q, Brown JB, Huang H, Bickel PJ (2011) Measuring reproducibility of high-throughput experiments. Annals of Applied Statistics 5:1752-1779. DOI 10.1214/11-AOAS466. (IDR.)
Related Skills
clip-seq/clip-preprocessing - UMI extraction and adapter trimming details
clip-seq/clip-alignment - STAR ENCODE block + multi-mapper rescue
