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

## ----install, eval=FALSE------------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# 
# BiocManager::install("RBPSpecificity")

## ----load-package-------------------------------------------------------------
library(RBPSpecificity)
library(ggplot2)

## ----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)

## ----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)

## ----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)

## ----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)

## ----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

## ----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)

## ----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)

## ----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------------------------------------------------------------------
sessionInfo()