bio-clip-seq-binding-site-annotation

What this skill does

## Version Compatibility

Reference examples tested with: ChIPseeker 1.40+, RCAS 1.30+, GenomicFeatures 1.56+, GenomicRanges 1.56+, rbp-maps (Yeo github), bedtools 2.31+, 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. ChIPseeker 1.40+ changed the priority defaults; verify the priority vector before pipelining.

# Binding Site Annotation

**"Annotate where in transcripts my RBP binds"** -> Map CLIP peaks or single-nucleotide crosslink sites to RNA features and report the per-feature distribution. The interpretation is RBP-class-specific: splicing factors (PTBP1, U2AF2, RBFOX) bind intron-exon junctions; mRNA-stability regulators (HuR, PUM2) bind 3' UTRs; translation factors (EIF3J, RPS19) bind 5' UTRs and CDS; and small-ncRNA-binding RBPs (NSUN2 tRNAs, LARP7 7SK, TROVE2 Y-RNAs) bind specific non-coding transcripts. A correct annotation pipeline (a) resolves overlapping features by priority, (b) preserves transcript-isoform context, (c) generates metagene distributions, and (d) flags repeat-element overlap separately.

- R (peak-level, fast): `ChIPseeker::annotatePeak(peaks, TxDb=txdb, level='gene', tssRegion=c(-100,100))` then `plotAnnoPie(anno)`
- R (transcript-level with RNA-specific regions): RCAS `runReport(queryRegions=peaks, gffData=gencode_gtf, genomeVersion='hg38')`
- CLI (regions only): `bedtools intersect -s -wa -wb -a peaks.bed -b features.bed` (manual hierarchy)
- R (splicing-regulatory map for splice factors): RBP-Maps (`yeolab/rbp-maps`) generates the 1400 nt vectorized cassette-exon map used in Yeo lab ENCODE papers
- CLI (metagene aggregation): `deepTools computeMatrix` or RSeQC `geneBody_coverage.py` for read profiles

The Yeo lab convention for ENCODE eCLIP: ChIPseeker for global region distribution + RBP-Maps for position-resolved splicing context + custom analysis for repeat-element overlap. ChIPseeker alone over-counts intronic binding because TSS-region default extends too far for RNA features.

## Algorithmic Taxonomy

| Tool | Input | Annotation level | Feature priority | Output | Strength | Fails when |
|------|-------|------------------|------------------|--------|----------|------------|
| ChIPseeker (R) | Peak BED + TxDb | Gene-level or transcript-level | Promoter > 5' UTR > 3' UTR > Exon > Intron > Intergenic | Per-peak annotation + pie chart + distance-to-TSS | Mature, widely used, fast | TSS-region default `c(-3000, 3000)` over-extends; promoter category meaningless for RNA |
| RCAS (R) | Peak BED + GFF | Transcript-level | Customizable | HTML report + per-region table | RNA-specific design; ncRNA-aware | Slow; HTML-heavy; less actively maintained |
| RBP-Maps (Yeo) | eCLIP BAM + alternative splicing event tables | Splice junction-level | NA (positional metagene) | Splicing regulatory map (1400nt vector around cassette exon) | The standard for splicing-factor CLIP analysis | Splicing-only; not for 3' UTR or ncRNA binders |
| ASTERIX (Yeo) | Peak BED | Region | NA | Region distribution + region-specific motifs | Companion to RBP-Maps for non-splicing regions | Less granular than ChIPseeker |
| HOMER annotatePeaks.pl | Peak BED | Gene-level | Promoter > UTR > Exon > Intron > Intergenic | Annotation TSV | Familiar from ChIP-seq | Same TSS issue as ChIPseeker; less RNA-aware |
| bedtools intersect | Peak BED + feature BED | User-defined | User-defined | Per-feature overlap | Most flexible | Manual hierarchy logic must be written |
| pyranges (Python) | Peak BED + GTF | Customizable | User-defined | DataFrame | Fast, pythonic, scriptable | Hand-written priority logic |
| pybedtools | Peak BED + GTF | Customizable | User-defined | iterators | Flexible, scripted | Same as pyranges |
| Peakhood (Reuter 2022) | Peak BED + transcripts | Transcript context | NA | Per-peak transcript context | Resolves transcript-isoform ambiguity | Specialized; not a global region tool |

Methodology evolves; verify the GTF source (GENCODE preferred over Ensembl for human; both have feature differences) and the priority hierarchy against published RBP literature.

## Critical Choice: Annotation Hierarchy

RNA features are nested and overlapping. A 3' UTR peak is ALSO in an exon; an intronic peak near a splice site is ALSO in a transcribed-but-spliced region. Without a priority hierarchy the same peak gets counted multiple times.

**Standard CLIP hierarchy (most-specific first):**

1. **3' UTR** (specific to mature mRNA biology)
2. **5' UTR**
3. **CDS** (coding sequence)
4. **Exon** (catch-all for non-UTR exonic)
5. **Splice site** (within 50 nt of donor/acceptor; key for splicing factors)
6. **Intron**
7. **Promoter** (within 1 kb of TSS; usually not relevant for RNA-binding)
8. **5' flanking** / **3' flanking** (within 1 kb of gene)
9. **ncRNA** (lncRNA, snoRNA, snRNA, miRNA host)
10. **Repeat element** (Alu, LINE, LTR, SINE) - separate axis from above
11. **Intergenic**

ChIPseeker's default `level='gene'` reports the FIRST hit by priority. For CLIP, set `level='transcript'` and explicitly customize `tssRegion=c(-100, 100)` (default `c(-3000, 3000)` is for ChIP).

**Why CLIP needs the tight TSS window:** ChIPseeker was designed for ChIP-seq where the "Promoter" category captures transcription-factor binding within 3 kb upstream of a gene. For CLIP-seq, the binding substrate is RNA - the RBP cannot bind upstream of the TSS because there is no RNA upstream of the TSS. A 6 kb TSS window catches 5'-UTR peaks and unrelated upstream sequence, mislabeling 30-50% of CLIP peaks as "Promoter (2-3kb)" - implausible for an RNA-binding study. Tightening to `c(-100, 100)` constrains the "Promoter" category to the mRNA 5' end immediate vicinity and pushes deeper-5'-UTR peaks into the "5UTR" category where they belong.

```r
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

peaks <- readPeakFile('peaks.stringent.bed')

# CLIP-appropriate annotation
anno <- annotatePeak(
    peaks,
    TxDb = txdb,
    level = 'transcript',     # RNA-isoform-aware
    tssRegion = c(-100, 100), # tight; not a meaningful promoter for RNA binding
    genomicAnnotationPriority = c(
        'Promoter', '5UTR', '3UTR', 'Exon',
        'Intron', 'Downstream', 'Intergenic'
    )
)

print(anno)
plotAnnoPie(anno)
plotAnnoBar(anno)
plotDistToTSS(anno)  # CLIP-specific: dist to TSS is dist to mRNA 5' end
```

## RBP Class -> Expected Annotation

| RBP class | Expected dominant region | RBP examples |
|-----------|--------------------------|--------------|
| Splicing factors | Intron (within 500 nt of splice junction) | PTBP1, U2AF2, RBFOX1/2, SRSF1-9, HNRNPC, MBNL1 |
| 3' UTR mRNA stability/decay | 3' UTR | HuR (ELAVL1), PUM2, AUF1, ZFP36, TIA1, NUDT21 |
| 5' UTR translation regulation | 5' UTR, CDS start | EIF4A3, EIF3J, eIF2 subunits, LARP1, PCBP1/2 |
| Ribosome-associated | CDS | RPL/RPS proteins, RACK1 |
| Nuclear export | CDS, 3' UTR | NXF1, ALYREF, SRSF3 |
| m6A reader | 3' UTR, stop codon region | YTHDF1/2/3, IGF2BP1/2/3 |
| miRNA effector | 3' UTR (seed-matched sites) | AGO1/2/3/4, DICER1, TNRC6A/B/C |
| snoRNA-targeting RBP | snoRNA, rRNA | DKC1, FBL, NOP56, NHP2 |
| ncRNA-binding | Specific ncRNA target | TROVE2 (Y RNA), LARP7 (7SK), SSB (vault) |
| Mitochondrial | chrM (mt-mRNA + 12S/16S rRNA) | FASTKD2, LRPPRC, TFAM, MTPAP |
| Histone mRNA | Histone 3' end stem-loop | SLBP |
| Repeat-binding (TE) | Alu, LINE-1, LTR, SINE | MATR3, ZFP36, HNRNPK, HNRNPA2B1, FUS (LINE-1) |
| LCD/IDR / phase separation | Diverse (Alu, intron, 3' UTR) | FUS, TDP-43, EWSR1, TAF15 |

Mismatch between observed and expected sugge