---
title: "Introduction to RBPSpecificity"
author: "Soon Yi"
date: "`r Sys.Date()`"
output: 
  BiocStyle::html_document:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Introduction to RBPSpecificity}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
    echo = TRUE,
    collapse = TRUE,
    comment = "#>",
    fig.width = 7,
    fig.height = 5,
    warning = FALSE,
    message = FALSE
)
```

# Overview

**RBPSpecificity** provides tools to analyze RNA-Binding Protein (RBP) binding 
specificities from high-throughput sequencing data. This package provides tools for:

- Calculating **Specificity** from k-mer enrichment data
- Calculating **Sensitivity** for top motifs
- Visualizing specificity distributions and sensitivity profiles
- Computing enrichment scores from CLIP peak data

# Introduction

## Background

The package implements two key metrics for characterizing RBP binding behavior:

1. **Specificity**: Quantifies how specifically an RBP binds its 
   top motif relative to all other motifs. Higher specificity indicates more specific binding.
2. **Sensitivity**: Measures how much the binding affinity 
   changes upon single nucleotide variations in the motif. Higher sensitivity indicates
   greater sensitivity to sequence changes.

These metrics enable quantitative comparison of RBP binding properties across
different experimental conditions or between different RBPs.

## Reference

For a detailed description of these metrics and their biological applications, 
please refer to:

> Yi S, Singh SS, Ye X, Krishna R, Kothwela V, Jankowsky E, Luna JM. (2025). 
> *Inherent Specificity and Mutational Sensitivity as Quantitative Metrics for RBP Binding.* 
> bioRxiv. https://www.biorxiv.org/content/10.1101/2025.03.28.646018v4

## Workflow Integration

**RBPSpecificity** is designed to seamlessly integrate into standard Bioconductor 
sequence and genomic analysis workflows:

- **Input Integration**: Raw peak files from CLIP-seq or similar experiments can be 
  imported as `GRanges` objects using `rtracklayer::import()` and standard annotation 
  packages (e.g. `GenomicFeatures`).
- **Sequence Retrieval**: The `motifEnrichment()` function interfaces directly with 
  `BSgenome` genome packages (e.g., `BSgenome.Hsapiens.UCSC.hg38`) to extract genomic 
  sequences for enrichment calculations.
- **Downstream Analysis**: The Specificity and Sensitivity scores and matrices returned 
  by the package can be integrated with other downstream analyses such as differential 
  binding studies (e.g., `DiffBind`, `DESeq2`), regulatory network modeling, or motif 
  annotation pipelines.

## Sample Data

The sample data included in this package were obtained from the ENCODE portal.

### RBNS Data

RNA Bind-n-Seq (RBNS) datasets provide in vitro RNA-RBP affinity measurements. 
For each K-mer enrichment dataset, the RBP concentration generating the highest 
R-value was selected as the representative sample. The collected RBNS R scores 
were feature-scaled to [1, e] followed by natural log transformation to normalize 
to a 0-to-1 scale.

### eCLIP Data

eCLIP datasets provide in vivo RBP binding site information. Peak coordinates 
were extended 25 nucleotides upstream from the 5'-end. Motif enrichment was 
calculated by counting K-mer appearances across extended peaks and subtracting 
background counts from randomly shifted genomic regions (up to 500 nt away). 
Background generation was repeated 100 times for robust averaging.

## Setup

```{r install, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("RBPSpecificity")
```

```{r load-package}
library(RBPSpecificity)
library(ggplot2)
```

# In Vitro Specificity from Pre-computed Enrichment Data (RBNS)

This workflow demonstrates how to calculate specificity and sensitivity from pre-computed 
enrichment scores, such as those from RBNS (RNA Bind-n-Seq) experiments.

## Loading RBNS Data

The sample data contains normalized 5-mer enrichment scores for four RBPs:
HNRNPC, PCBP2 (structural context independent), RBFOX2, and EIF4G2 
(structural context dependent).

```{r load-rbns}
# Load sample RBNS data from package
rbns_file <- system.file("extdata", "RBNS_normalized_5mer.csv",
    package = "RBPSpecificity"
)
rbns_data <- read.csv(rbns_file, stringsAsFactors = FALSE)

# Preview the data
head(rbns_data)
```

## Calculating Specificity and Sensitivity for All Four RBPs

The package expects a data.frame with `MOTIF` and `Score` columns:

```{r calc-rbns-metrics}
rbps <- c("HNRNPC", "PCBP2", "RBFOX2", "EIF4G2")

rbns_metrics <- lapply(rbps, function(rbp) {
    rbns_enrichment <- data.frame(
        MOTIF = rbns_data$Motif,
        Score = rbns_data[[rbp]]
    )
    is_val <- returnSpecificity(rbns_enrichment)
    vs_val <- returnSensitivity(rbns_enrichment, output_type = "number")
    list(
        enrichment = rbns_enrichment,
        Specificity = is_val,
        Sensitivity = vs_val
    )
})
names(rbns_metrics) <- rbps

# Display in vitro metrics
invitro_summary <- data.frame(
    RBP = rbps,
    Specificity = sapply(rbns_metrics, function(x) round(x$Specificity, 2)),
    Sensitivity = sapply(rbns_metrics, function(x) round(x$Sensitivity, 4))
)
print(invitro_summary)
```

## Visualization

### Specificity Distribution Plots

The `plotSpecificity()` function displays the distribution of k-mer scores, highlighting 
the top motif and its specificity value:

```{r plot-is-rbns, fig.cap="Distribution of k-mer scores from RBNS data for four RBPs."}
plotSpecificity(rbns_metrics[["HNRNPC"]]$enrichment)
plotSpecificity(rbns_metrics[["PCBP2"]]$enrichment)
plotSpecificity(rbns_metrics[["RBFOX2"]]$enrichment)
plotSpecificity(rbns_metrics[["EIF4G2"]]$enrichment)
```

### Sensitivity Profile Plots

The `plotSensitivity()` function visualizes how sensitive binding is to mutations at each 
position:

```{r plot-vs-rbns, fig.cap="Sensitivity profiles from RBNS data for four RBPs."}
plotSensitivity(rbns_metrics[["HNRNPC"]]$enrichment)
plotSensitivity(rbns_metrics[["PCBP2"]]$enrichment)
plotSensitivity(rbns_metrics[["RBFOX2"]]$enrichment)
plotSensitivity(rbns_metrics[["EIF4G2"]]$enrichment)
```

# Cellular Specificity from De Novo Enrichment (eCLIP)

This workflow demonstrates how to calculate enrichment scores directly from 
peak data using the `motifEnrichment()` function.

> **Note**: While this example uses eCLIP-derived peaks, `motifEnrichment()` 
> works with any peak data in BED format (e.g., iCLIP, PAR-CLIP, or other 
> ChIP-seq derived peaks).

> **Requirement**: This workflow requires `BSgenome.Hsapiens.UCSC.hg38` to be installed.

When calculating enrichment from cellular peak data, the resulting metrics are 
referred to as **Cellular Specificity** and **Cellular Sensitivity**,
as they reflect in vivo binding behavior rather than intrinsic biochemical affinity.

## Loading Peak Data

```{r load-eclip, eval = requireNamespace("BSgenome.Hsapiens.UCSC.hg38", quietly = TRUE)}
# Load eCLIP peaks (10-column narrowPeak format) from package
bed_cols <- c(
    "chr", "start", "end", "name", "score", "strand",
    "signalValue", "pValue", "qValue", "peak"
)

rbp_names <- c("HNRNPC", "PCBP2", "RBFOX2", "EIF4G2")
cell_lines <- c("K562", "HepG2", "K562", "K562")

peak_data <- lapply(seq_along(rbp_names), function(i) {
    bed_file <- system.file("extdata",
        sprintf("ENCODE_eCLIP_%s_%s_narrowPeak.bed", rbp_names[i], cell_lines[i]),
        package = "RBPSpecificity"
    )
    peaks <- read.table(bed_file, header = FALSE, sep = "\t")
    colnames(peaks) <- bed_cols
    message("Loaded ", nrow(peaks), " ", rbp_names[i], " peaks.")
    peaks
})
names(peak_data) <- rbp_names
```

## Running motifEnrichment Pipeline

In addition to standard input variables, `motifEnrichment()` provides several 
parameters for customization:

| Parameter | Description |
|-----------|-------------|
| `extension` | Numeric vector of length 2: `c(five_prime, three_prime)` indicating shifting of 5' and 3' ends. Positive extends, negative trims. |
| `method` | Enrichment model: `"anr"` (default, any number of repetitions, conditional Binomial test) or `"zoops"` (zero-or-one occurrence per sequence, Hypergeometric test) |
| `bkg_iter` | Number of background iterations for averaging (higher = more robust) |
| `bkg_min_dist` | Minimum distance (bp) to shift peaks when generating background regions |
| `bkg_max_dist` | Maximum distance (bp) for shifting peaks in background generation |
| `scramble_bkg` | Logical, scramble background sequences to control for nucleotide composition (default: FALSE) |
| `nucleic_acid_type`| Character string, `"DNA"` or `"RNA"` (default: `"DNA"`) |

ZOOPS method can also be used instead. Note that ANR and ZOOPS calculation can 
result in different motif enrichment values and thus different specificity and 
sensitivity calculations. For more information on ANR vs. ZOOPS, please check 
the GitHub README and Further Reading section. In addition, the scramble 
background option can be turned on (`scramble_bkg = TRUE`) to control for 
nucleotide composition.

```{r run-enrichment, eval = requireNamespace("BSgenome.Hsapiens.UCSC.hg38", quietly = TRUE)}
# NOTE: This requires BSgenome.Hsapiens.UCSC.hg38

cellular_results <- lapply(rbp_names, function(rbp) {
    message("Running motifEnrichment for ", rbp, "...")
    enrichment <- motifEnrichment(
        coordinates = peak_data[[rbp]],
        species_or_build = "hg38",
        K = 5,
        extension = c(25, 0),
        method = "anr",
        bkg_iter = 100,
        scramble_bkg = FALSE
    )

    cs_val <- returnSpecificity(enrichment)
    cvs_val <- returnSensitivity(enrichment, output_type = "number")

    list(
        enrichment = enrichment,
        Specificity = cs_val,
        Sensitivity = cvs_val
    )
})
names(cellular_results) <- rbp_names

# Display cellular metrics
cellular_summary <- data.frame(
    RBP = rbp_names,
    Specificity = sapply(cellular_results, function(x) round(x$Specificity, 2)),
    Sensitivity = sapply(cellular_results, function(x) round(x$Sensitivity, 4))
)
print(cellular_summary)
```

## Visualization

The same `plotSpecificity()` and `plotSensitivity()` functions used for RBNS data can be applied
to visualize cellular specificity and cellular sensitivity from eCLIP enrichment results:

```{r eclip-plots, eval = requireNamespace("BSgenome.Hsapiens.UCSC.hg38", quietly = TRUE)}
# Specificity Distribution plots
plotSpecificity(cellular_results[["HNRNPC"]]$enrichment)
plotSpecificity(cellular_results[["PCBP2"]]$enrichment)
plotSpecificity(cellular_results[["RBFOX2"]]$enrichment)
plotSpecificity(cellular_results[["EIF4G2"]]$enrichment)

# Sensitivity Profile plots
plotSensitivity(cellular_results[["HNRNPC"]]$enrichment)
plotSensitivity(cellular_results[["PCBP2"]]$enrichment)
plotSensitivity(cellular_results[["RBFOX2"]]$enrichment)
plotSensitivity(cellular_results[["EIF4G2"]]$enrichment)
```

# Comparing In Vitro (RBNS) vs. Cellular (eCLIP) Specificity

A key application is comparing specificity and sensitivity across in vitro and cellular contexts. RBPs that bind directly to single-stranded sequence motifs (e.g., HNRNPC, PCBP2) are expected to show correlated specificity and sensitivity values. RBPs whose binding depends on structural context (e.g., RBFOX2, EIF4G2) may show divergence between in vitro and cellular metrics.

```{r comparison, eval = requireNamespace("BSgenome.Hsapiens.UCSC.hg38", quietly = TRUE)}
comparison <- data.frame(
    RBP = rbp_names,
    Structural_Context = c("Independent", "Independent", "Dependent", "Dependent"),
    In_Vitro_Specificity = sapply(rbns_metrics[rbp_names], function(x) round(x$Specificity, 2)),
    Cellular_Specificity = sapply(cellular_results, function(x) round(x$Specificity, 2)),
    In_Vitro_Sensitivity = sapply(rbns_metrics[rbp_names], function(x) round(x$Sensitivity, 4)),
    Cellular_Sensitivity = sapply(cellular_results, function(x) round(x$Sensitivity, 4))
)
print(comparison)
```

# Session Info

```{r session}
sessionInfo()
```