Clip qc

Version Compatibility

Reference examples tested with: preseq 3.2+, picard 3.1+, samtools 1.19+, bedtools 2.31+, deeptools 3.5+, idr 2.0.4+, MultiQC 1.21+, RSeQC 5.0+, pysam 0.22+, fastp 0.23+.

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.

CLIP-seq Quality Control

"Did my CLIP library pass?" -> Assess preprocessing retention, alignment rate, library complexity, replicate reproducibility (IDR), fraction reads in peaks (FRiP), read-distribution metagene, rRNA/snoRNA contamination, fragment-length distribution, and SMInput vs IP enrichment. ENCODE eCLIP compliance is the canonical bar: >= 1M unique fragments per replicate, IDR rescue and self-consistency ratios both < 2, FRiP >= 0.005 (narrow-binding), library complexity rising linearly with depth on preseq lc_extrap. A library can fail at any of these stages, and the failure mode determines whether the data is salvageable.

• CLI (library complexity, primary QC): preseq lc_extrap -B -P aligned.bam -o complexity.txt
• CLI (FRiP after peak calling): bedtools intersect -c -s -a peaks.bed -b dedup.bam | awk '{s+=$NF} END{print s}' then divide by total reads
• CLI (IDR, ENCODE convention): idr --samples rep1.sorted.bed rep2.sorted.bed --input-file-type bed --rank 5 --output-file idr.out --idr-threshold 0.05 --plot
• CLI (read distribution / metagene): RSeQC read_distribution.py -i dedup.bam -r gencode.v38.bed + geneBody_coverage.py -i dedup.bam -r housekeeping.bed -o gb
• CLI (rRNA contamination check): samtools idxstats dedup.bam | awk '$1 ~ /rRNA|45S|18S|28S/ { sum+=$3 } END {print sum}'
• CLI (consolidated report): multiqc <run_dir> aggregates FastQC + cutadapt + STAR + umi_tools + preseq + samtools stats

The ENCODE eCLIP standards (encodeproject.org/eclip) define: >= 2 biological replicates with >= 1M unique fragments each (or saturated peak detection); IDR rescue and self-consistency ratios both < 2; narrow-binding RBPs FRiP >= 0.005. CLIP libraries should have 40-70% PCR duplication BY DESIGN - the IP enriches a small molecule pool, so high pre-dedup duplication is normal; low duplication suggests failed IP.

QC Stage Hierarchy

CLIP QC progresses through five gates; failure at an earlier gate makes later gates meaningless.

Gate	Metric	Tool	ENCODE threshold	Failure interpretation
1. Preprocessing retention	% reads retained after UMI + adapter trim	cutadapt log	>= 70%	Adapter pattern wrong; degraded RNA
2. Alignment rate	% reads aligned to genome (unique)	STAR Log.final.out	>= 60% for eCLIP, 70% for iCLIP/PAR-CLIP	Wrong genome; rRNA pre-map missing
3. Library complexity	Predicted unique fragments at sequenced depth	preseq lc_extrap	>= 1M unique	Over-amplified or under-input library
4. IP enrichment	log2(IP/SMInput) at expected sites; FRiP	bedtools + idr	FRiP >= 0.005; log2 >= 3 at top peaks	Failed antibody / antibody not IP-grade
5. Reproducibility	IDR rescue + self-consistency ratios	idr	both < 2	Biological variation too high; or low complexity

A library failing at Gate 3 (complexity) cannot be rescued analytically; gates 4-5 fail downstream of complexity by construction.

Library Complexity with preseq

Goal: Determine whether the CLIP library captured enough independent molecules to support genome-wide peak calling (ENCODE: >= 1M unique fragments per replicate).

Approach: Run preseq lc_extrap on the PRE-dedup BAM (preseq counts PCR duplicates to extrapolate); also compute picard ESTIMATED_LIBRARY_SIZE at sequenced depth. Flag any library predicted to plateau below 3M unique fragments at infinite depth.

# After alignment, BEFORE UMI dedup (preseq counts PCR duplicates)
preseq lc_extrap \
    -B -P \
    -o sample_complexity.txt \
    sample_aligned.bam

# Output columns:
# TOTAL_READS  EXPECTED_DISTINCT  LOWER_0.95CI  UPPER_0.95CI
# At 100M reads, EXPECTED_DISTINCT:
#   >= 10M  = excellent complexity
#   3-10M   = acceptable; restrict to high-expression transcripts
#   < 3M    = library failed; cannot rescue analytically

# picard direct estimate at current depth
picard EstimateLibraryComplexity \
    I=sample_aligned.bam \
    O=picard_complexity.txt
# ESTIMATED_LIBRARY_SIZE > 5M = healthy CLIP library

A linear plateau on preseq's curve at low depth indicates over-amplification; a curve still climbing at sequenced depth means more sequencing would yield more unique fragments.

FRiP (Fraction Reads in Peaks)

FRiP measures how much of the IP signal falls into the called peak set. ENCODE eCLIP narrow-binding RBP minimum: FRiP >= 0.005. Atypical-binding RBPs (rare-transcript binders like TROVE2 on Y RNAs) are exempt.

# Reads in peaks (use stringent peaks: log2 FC >= 3, -log10 p >= 3)
reads_in_peaks=$(bedtools intersect -c -s -a peaks.stringent.bed -b dedup.bam | awk '{s+=$NF} END {print s}')
total_reads=$(samtools view -c -F 4 dedup.bam)
frip=$(echo "scale=4; $reads_in_peaks / $total_reads" | bc)
echo "FRiP: $frip"

# Per-region FRiP breakdown
for region in three_utr exon intron; do
    rip=$(bedtools intersect -c -s -a peaks_${region}.bed -b dedup.bam | awk '{s+=$NF} END {print s}')
    echo "${region}: $(echo "scale=4; $rip / $total_reads" | bc)"
done

RBP class	Expected FRiP (ENCODE eCLIP)
Splicing factors (PTBP1, U2AF2)	0.01 - 0.10
3' UTR mRNA stability (HuR, PUM2)	0.02 - 0.20
Translation factors (EIF3J)	0.01 - 0.05
Repeat binders (MATR3)	0.05 - 0.30 (high; concentrated in repeats)
Mitochondrial (FASTKD2)	0.05 - 0.40 (very high; chrM is small)
snoRNA binders (DKC1)	0.10 - 0.50 (high; snoRNA is rare)
Failed IP (any RBP)	< 0.005

IDR for CLIP Reproducibility

IDR (Li et al 2011) measures peak-rank reproducibility across replicates. ENCODE eCLIP convention applies IDR identically to ChIP-seq, using CLIPper + SMInput log2 FC + -log10 p as the ranking signal. The CLIP-specific consideration: rank by signalValue (log2 FC) or p-value, NOT by score column (CLIPper score is sparse and tied).

# Sort each replicate's peaks by signal
sort -k5,5gr rep1.compressed.bed > rep1.sorted
sort -k5,5gr rep2.compressed.bed > rep2.sorted

# True-replicates IDR (threshold 0.05)
idr --samples rep1.sorted rep2.sorted \
    --input-file-type bed --rank 5 \
    --output-file idr_true.out \
    --idr-threshold 0.05 \
    --plot --log-output-file idr.log

# Pseudo-replicates from each individual replicate (split BAM in half)
samtools view -b -h -s 1.5 rep1.dedup.bam > rep1.psr1.bam   # seed 1, fraction 0.5
samtools view -b -h -s 2.5 rep1.dedup.bam > rep1.psr2.bam   # seed 2 (different)
# Re-run peak calling on each pseudoreplicate, then IDR at threshold 0.10

ENCODE consistency rules for eCLIP:

• Nt = peaks passing IDR on true replicates
• Nself = peaks passing IDR on pseudo-replicates of each rep
• Library passes if: max(Nt, Nself) / min(Nt, Nself) <= 2
• If both ratios > 2: library rejected

SMInput vs IgG Control: Which?

The eCLIP design uses a size-matched input (SMInput) from the SAME lysate, treated identically (UV, IP buffer, RNase, ligation, IP, but with NO antibody addition - just bead-only control or a non-specific control IP). This is fundamentally different from IgG controls and from RNA-seq.

Control	What it measures	Pros	Cons
SMInput	Background from non-specific binding + ligation/RT/gel biases at same size	Captures all CLIP-specific biases; ENCODE standard	Requires same-day prep; cannot use a previous IgG library
IgG-IP	Non-specific antibody binding to the same RBP-naive lysate	Direct nonspecificity measure	Yields very low (~3-10x less reads); high PCR dup; hard to normalize
Empty beads (mock)	Bead surface non-specificity only	Cleanest baseline	Misses real CLIP background (IP buffer + ligation step bias)
RNA-seq (matched cell type)	Transcript abundance	Easy to obtain	Misses CLIP-specific biases entirely; not a real CLIP control
Total RNA / nuclear RNA	Cellular RNA distribution	Easy	Same as RNA-seq above

Consensus (ENCODE / Hentze / Yeo): SMInput. The bead-only IgG/mock alternatives are ill-suited for quantification because their library yields are 5-10x lower, dominated by PCR duplicates, and produce sparse read-density tracks. RNA-seq cannot replace SMInput because it does not capture the non-specific binding that occurs during IP.

# SMInput preparation - same as IP except no antibody added during incubation
# Same UV dose, same lysate aliquot, same RNase, same library prep
# Critical: same SDS-PAGE size cut from membrane as the IP

# Check SMInput vs IP enrichment globally
# Whole-genome log2(IP_RPKM / SMInput_RPKM):
# > 0.5 globally indicates IP enrichment of expressed transcripts (good)
# ~ 0 indicates failed IP (SMInput == IP)

Read Distribution Metagene

# RSeQC read_distribution.py shows fractional read placement
read_distribution.py -i dedup.bam -r gencode.v38.bed > read_dist.txt

# Output reports:
#   CDS_Exons, 5'UTR_Exons, 3'UTR_Exons, Introns, TES_down_10kb, TES_down_1kb,
#   TSS_up_10kb, TSS_up_1kb, Intergenic_region
# CLIP-seq biology-specific patterns:
#   Splicing factor: > 60% Introns + 5' UTR exons (containing 5' splice sites)
#   3' UTR factor: > 50% 3'UTR_Exons
#   m6A reader: 3'UTR_Exons + Stop_codon region
#   Failed IP: matches RNA-seq distribution (no enrichment)

# GeneBody coverage for 5' vs 3' bias
geneBody_coverage.py \
    -i dedup.bam \
    -r housekeeping.bed \
    -o sample_gb
# Flat curve = no positional bias (normal for most RBPs)
# 3' end bias = polyA-dependent enrichment (suspect)
# 5' end bias = nascent / TSS-proximal (only sensible for some RBPs)

Fragment-Length Distribution (Paired-End)

# Paired-end fragment-length distribution
samtools view -f 2 dedup.bam | awk '$9 > 0 && $9 < 500 {print $9}' | sort -n | uniq -c > fragment_lengths.txt

# CLIP normal range: 20-75 nt insert (rises from short trimmed reads to ~75 nt)
# Wider range (20-200) seen in high-RNase / long-fragment protocols
# Narrow peak at one length (e.g., all reads at 30 nt) = over-trimmed
# Bimodal at 30 nt and 150 nt = library preparation artifact

Pre-Map rRNA / snoRNA Contamination Check

# rRNA reads as fraction of total
total=$(samtools view -c -F 4 dedup.bam)
rRNA=$(samtools idxstats dedup.bam | awk '$1 ~ /rRNA|45S|18S|28S|5_8S/ { sum+=$3 } END {print sum}')
rrna_frac=$(echo "scale=4; $rRNA / $total" | bc)
echo "rRNA fraction: $rrna_frac"

# eCLIP without pre-map: rRNA 5-30% normal
# After pre-map: < 2% expected
# > 30% indicates rRNA dominance - IP captured ribosomes preferentially
# Some RBPs (RPL/RPS) expected to bind ribosomes; verify against RBP biology

Antibody Validation Sanity Check

The antibody is the single most common point of failure in CLIP. Even commercial "IP-grade" antibodies fail at 30-50% rates in practice. Cross-check IP enrichment against expected biology:

import pandas as pd

# Load peak file
peaks = pd.read_csv('peaks.stringent.bed', sep='\t', header=None,
                    names=['chr','start','end','name','log2fc','strand'])

# Filter top 100 peaks
top_peaks = peaks.nlargest(100, 'log2fc')

# Check enrichment in expected biology
# 1. If RBP is splicing factor: top peaks should be intronic or splice-site flanking
# 2. If RBP is HuR: > 70% top peaks should be 3' UTR
# 3. If RBP is FASTKD2: top peaks should be chrM
# 4. If RBP is FUS: GUGGU motif should appear in top motif enrichment

# Quick check
chrom_dist = top_peaks['chr'].value_counts(normalize=True)
print(chrom_dist.head(10))
# Healthy RBP: chromosomes represented in proportion to expression
# Failed IP: top chromosome chrM (mt artifact) or chr21/22 (housekeeping bias) > 20%

GO term sanity: top-peak genes should enrich for the expected biology (e.g., HuR -> immune / inflammation / mRNA stability terms; PTBP1 -> RNA splicing terms). If the top GO term is "cellular metabolism" or "translation" for a splicing factor, the IP failed.

Per-Stage Failure Modes

Gate 1: Preprocessing retention < 70%

Trigger: cutadapt log shows > 30% reads filtered (too short or no adapter).

Mechanism: Adapter pattern wrong; or RNA was degraded before fragmentation (all reads adapter-only).

Symptom: Catastrophic loss at the adapter trim step.

Fix: Verify adapter sequence in library prep documentation; check library quality control (Bioanalyzer / TapeStation) for RIN value; re-prep if RIN < 7.

Gate 2: Alignment rate < 60%

Trigger: STAR Log.final.out reports Uniquely mapped reads % < 60% for eCLIP, < 70% for iCLIP.

Mechanism: rRNA contamination dominating (no pre-map); wrong species genome; or chimeric library (contamination).

Symptom: Low alignment rate; samtools idxstats shows most reads as rRNA-aligned.

Fix: Run bowtie2 pre-map to rRNA + repeats index; verify genome species; check for sample swap with samtools view -h sample.bam | grep '^@SQ'.

Gate 3: Library complexity < 1M unique

Trigger: preseq lc_extrap predicts plateau < 1M unique at sequenced depth; or picard reports ESTIMATED_LIBRARY_SIZE < 1M.

Mechanism: Over-amplification (> 25 PCR cycles); low input (< 5M cells); failed IP capturing only a few molecules.

Symptom: preseq curve flattens early; UMI clusters have median size > 8.

Fix: No analytic rescue. Re-prep with more input cells (10-20M for eCLIP), fewer PCR cycles (14-18 for eCLIP, 16-20 for iCLIP2). If forced to use this library, restrict downstream analysis to high-expression transcripts and acknowledge the limitation.

Gate 4: FRiP < 0.005 for narrow-binding RBP

Trigger: FRiP value below ENCODE threshold for the RBP class.

Mechanism: IP failed to enrich specific binding sites; antibody cross-reactive or not IP-grade.

Symptom: Top peaks dominated by abundant transcripts (GAPDH, ACTB); GO enrichment generic; motif analysis returns AU-rich background.

Fix: Re-test antibody on a knockdown lysate (siRNA against the RBP) - if WB signal does not decrease, antibody is non-specific; switch antibody. ENCODE-validated antibodies are at encodeproject.org/biosamples.

Gate 4 variant: IP/SMInput global log2 ~ 0

Trigger: Whole-genome log2(IP/SMInput) is near 0 instead of positive.

Mechanism: IP and SMInput are equivalent - the antibody is not enriching any RBP-bound RNA.

Symptom: Per-peak log2 FC scaled around 0; no obvious enrichment at known motif sites.

Fix: Same as above - antibody failure. Consider tagged-RBP system (Halo-CLIP, GoldCLIP, or knock-in endogenous tag).

Gate 5: IDR fails - rescue/self-consistency > 2

Trigger: ENCODE IDR test reports rescue ratio > 2 OR self-consistency > 2.

Mechanism: Replicates inconsistent. Biological variation high; OR one replicate has lower library complexity; OR one replicate failed IP.

Symptom: Per-replicate peak counts differ > 2x; IDR plot shows the two replicates have very different peak rank distributions.

Fix: Run preseq and FRiP on each replicate independently; identify the failing one; re-sequence to higher depth or re-prep. ENCODE-style: down-sample both replicates to common depth, re-test IDR.

High PCR duplication confused with low complexity

Trigger: Saw 60% PCR duplication and panicked.

Mechanism: CLIP libraries have 40-70% PCR duplication BY DESIGN - the IP enriches a small pool. UMI dedup recovers the unique molecules underneath.

Symptom: "60% duplication" headline. Actual unique count is what matters.

Fix: Confirm UMI dedup has been applied and unique fragment count >= 1M. The duplication rate metric is meaningless without UMI context for CLIP.

Decision Tree by Failure

Symptom	Likely failure	Action
> 50% reads "too short" in cutadapt	Adapter wrong or RNA degraded	Verify adapter; check RIN
Most reads align to rRNA	No pre-map; or IP captured ribosomes	bowtie2 pre-map; or accept if RBP is ribosomal
Unique frag count < 1M	Over-amplified or under-input	No rescue; re-prep
FRiP < 0.005 (narrow RBP)	Failed IP	Switch antibody
IP/SMInput global log2 ~ 0	Antibody non-specific	Switch antibody or use tagged-RBP
IDR rescue > 2	Replicate inconsistency	Identify failing replicate; re-sequence
Top GO terms generic	Failed IP or contamination	Check antibody on KD lysate WB
chrM peaks abundant for non-mt-RBP	Mitochondrial contamination	Filter chrM or accept if mt-RBP
Fragment length all 30 nt	Over-trimmed	Loosen quality trim from -q 20 to -q 6
Strand-specific peaks lost	bedtools without -s	Add -s strand flag

Standard QC Report (MultiQC)

# Aggregate all QC tools into a single report
multiqc \
    fastqc_output/ \
    cutadapt_logs/ \
    star_logs/ \
    umi_tools_logs/ \
    preseq_output/ \
    samtools_stats/ \
    -o multiqc_report/

MultiQC reads logs from FastQC, cutadapt, STAR, umi_tools, preseq, samtools stats, picard, and others, and produces a single HTML report. For CLIP-seq, additionally run RSeQC read_distribution.py and document FRiP / IDR results manually.

ENCODE-Style Comprehensive QC Table

Metric	Threshold	Source
Biological replicates	>= 2	ENCODE eCLIP
Read length	>= 50 nt (PE)	ENCODE
Unique fragments per replicate	>= 1M (or saturated)	ENCODE
Preprocessing retention	>= 70%	Practitioner consensus
Genome alignment rate (unique)	>= 60% eCLIP, 70% iCLIP	Practitioner consensus
Library complexity (preseq)	>= 10M expected unique at 100M	ENCODE
FRiP (narrow-binding RBP)	>= 0.005	ENCODE
FRiP (atypical-binding RBP)	Exempt	ENCODE
log2(IP/SMInput) at stringent peaks	>= 3	ENCODE
-log10 p at stringent peaks	>= 3	ENCODE
IDR rescue ratio	< 2	ENCODE
IDR self-consistency ratio	< 2	ENCODE
rRNA fraction post pre-map	< 2%	Practitioner consensus
Read distribution match to RBP class	Y	Sanity check
Top motif consistent with literature	Y	Sanity check

Common Errors

Error / symptom	Cause	Solution
preseq "all reads have duplicates"	Pre-deduplicated BAM passed to preseq	Re-run on PRE-dedup BAM
FRiP value > 1.0	Peak set overlaps reads counted twice	Use unique-reads BAM; verify peak BED is non-redundant
IDR plot empty	Ranking column wrong; all peaks tied at same value	Rank by --rank 5 (log2 FC); sort BED first
MultiQC missing modules	Logs not in expected directory structure	Verify log paths; rerun multiqc with explicit -d
Read distribution: > 50% intergenic for mRNA RBP	TxDb missing transcripts; rRNA contamination	Update GENCODE; pre-map rRNA
GeneBody coverage 3' biased	polyA-selected library; not CLIP	Verify library prep was random hexamer
Fragment-length distribution narrow at 30 nt	Adapter aggressively trimmed at 5'	Loosen trim parameters
Antibody KD WB shows no decrease	Antibody non-specific	Switch antibody; use ENCODE-validated

References

• Van Nostrand EL et al 2016 Nat Methods 13:508 (eCLIP, ENCODE QC standards)
• Li Q et al 2011 Ann Appl Stat 5:1752 (IDR framework)
• Landt SG et al 2012 Genome Res 22:1813 (ChIP/CLIP QC guidelines, IDR Nself rule)
• Daley T & Smith AD 2013 Nat Methods 10:325 (preseq library complexity)
• Wang L et al 2012 Bioinformatics 28:2184 (RSeQC)
• Ewels P et al 2016 Bioinformatics 32:3047 (MultiQC)
• Smith T et al 2017 Genome Res 27:491 (UMI-tools)
• ENCODE eCLIP Data Standards (encodeproject.org/eclip) - canonical thresholds
• Van Nostrand EL et al 2020 Nature 583:711 (ENCODE 150 RBP QC patterns)

Related Skills

• clip-seq/clip-preprocessing - Gate 1 preprocessing logs
• clip-seq/clip-alignment - Gate 2 alignment metrics
• clip-seq/clip-peak-calling - Gates 4-5 peak QC + FRiP
• clip-seq/differential-clip - QC for differential analyses
• read-qc/quality-reports - General FastQC / MultiQC
• read-qc/contamination-screening - Cross-species contamination
• chip-seq/chipseq-qc - DNA-protein QC analogue