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

## ----load-packages------------------------------------------------------------
library(RBPEqBind)
library(data.table)
library(ggplot2)

## ----load-models--------------------------------------------------------------
# Load raw models from package sample data
model_file <- system.file("extdata", "model_RBP.csv", package = "RBPEqBind")
raw_models <- loadModel(model_file)

# View available RBPs
names(raw_models)

## ----view-raw-models, fig.cap="Raw score distributions for model RBPs"--------
viewModel(raw_models, rbp = c("HH", "HL", "LH", "LL"), bins = 100, alpha = 0.4)

## ----set-models---------------------------------------------------------------
# High affinity (H) = high Ka -> low Kd
# Low affinity (L) = low Ka -> high Kd
max_affinities <- c(
  "HH" = 100, "HL" = 100,  # Ka_max = 100 -> Kd_min = 0.01 nM
  "LH" = 10,  "LL" = 10    # Ka_max = 10  -> Kd_min = 0.1 nM
)

min_affinities <- c(
  "HH" = 0.001, "HL" = 0.001,  # Ka_min = 0.001 -> Kd_max = 1000 nM
  "LH" = 0.001, "LL" = 0.001
)

rbp_models <- setModel(raw_models, 
                        max_affinity = max_affinities, 
                        min_affinity = min_affinities)

## ----view-models-all, fig.cap="Ka distributions for all RBPs"-----------------
viewModel(rbp_models, rbp = c("HH", "HL", "LH", "LL"), metric = "Ka", 
           bins = 100, alpha = 0.4)

## ----view-models-subset, fig.cap="Comparing HH vs LL affinity distributions"----
viewModel(rbp_models, rbp = c("HH", "LL"), metric = "Ka", bins = 100, alpha = 0.5)

## ----generate-rna-------------------------------------------------------------
# Generate a random 100-nt RNA sequence
random_seq <- generateRNA(100)
cat("Generated sequence:", substr(random_seq, 1, 50), "...\n")

# Or use a specific reference sequence
seq <- "UAGCGGUGCGAUUGGCCCGUGGACCGCGUUUUUGCACUCAUCGUUUCGCACUAAGUACAUAUAGUUGCGACAAAGCCGCUUAUGAGUUGGGGGUAUAUUC"

## ----scenario1-simulate-------------------------------------------------------
# Define concentrations (nM) - competitive binding of all 4 RBPs
prot_concs <- c("HH" = 100, "HL" = 100, "LH" = 100, "LL" = 100)
rna_conc <- 10.0  # nM

# Detect k-mer size from models
k_size <- nchar(rbp_models$HH$motif[1])

# Run simulation
res <- simulateBinding(
  sequence = seq,
  rbp_models = rbp_models,
  protein_concs = prot_concs,
  rna_conc = rna_conc,
  k = k_size
)

# Add transcript column for visualization
res$transcript <- "Custom RNA"

head(res)

## ----scenario1-binding-raw, fig.cap="Raw density_fc (positions 20-60)"--------
plotBinding(res, rbp = c("HH", "HL", "LH", "LL"), metric = "density_fc", 
             transcript = "Custom RNA", xlim = c(20, 60))

## ----scenario1-binding-window, fig.cap="Density FC with 5-mer smoothing window (positions 20-60)"----
plotBinding(res, rbp = c("HH", "HL", "LH", "LL"), metric = "density_fc", 
             transcript = "Custom RNA", window = 5, xlim = c(20, 60))

## ----scenario1-heatmap, fig.cap="Heatmap of occupancy_fc (zoomed region)"-----
plotHeatmap(res, transcript = "Custom RNA", xlim = c(20, 60), 
             xaxis_type = "both", metric = "occupancy_fc")

## ----scenario1-export---------------------------------------------------------
# Export to JSON
tmp_json <- tempfile(fileext = ".json")
exportResults(res, output_file = tmp_json, format = "json")

# Export to CSV
tmp_csv <- tempfile(fileext = ".csv")
exportResults(res, output_file = tmp_csv, format = "csv")

# Clean up
unlink(tmp_json)
unlink(tmp_csv)

## ----scenario1-se-------------------------------------------------------------
se <- makeSE(res, rbp_models = rbp_models)
se

## ----scenario2-future-demo, eval = FALSE--------------------------------------
# # Optional: enable parallel processing across 2 workers
# library(future)
# plan(multisession, workers = 2)

## ----scenario2-grid-----------------------------------------------------------
# Define concentration grids
prot_grid <- list(
  HH = c(10, 100), 
  HL = c(10, 100),
  LH = c(10, 100),
  LL = c(10, 100)
)
rna_grid <- c(10)

res_grid <- simulateGrid(
  sequence = seq,
  rbp_models = rbp_models,
  protein_conc_grid = prot_grid,
  rna_conc_grid = rna_grid,
  k = k_size
)

## ----scenario2-binding--------------------------------------------------------
# Filter to specific concentration combination
plotBinding(
  results = res_grid,
  rbp = c("HH", "HL", "LH", "LL"),
  rna_conc = 10,
  protein_conc = c(HH=10, HL=10, LH=100, LL=100),
  metric = "density"
)

## ----scenario2-plot-grid------------------------------------------------------
# First run a 2-RBP simulation for cleaner visualization
rbp_models_2rbp <- rbp_models[c("HH", "LL")]

res_grid_2rbp <- simulateGrid(
  sequence = seq,
  rbp_models = rbp_models_2rbp,
  protein_conc_grid = list(HH = c(10, 50), LL = c(10, 50)),
  rna_conc_grid = c(10),
  k = k_size
)

plotGrid(
  results = res_grid_2rbp,
  rbp1 = "HH",
  rbp2 = "LL",
  rbp1_concs = c(10, 50),
  roi_range = c(29, 33),
  metric = "density"
)

## ----scenario3-simulate-------------------------------------------------------
# Load sample FASTA from package
fasta_file <- system.file("extdata", "test_transcripts.fa", package = "RBPEqBind")

res_fasta <- simulateBindingFasta(
  fasta_file = fasta_file,
  rbp_models = rbp_models,
  protein_concs = c(HH = 100, HL = 100, LH = 100, LL = 100),
  rna_conc = 10,
  k = 5
)
head(res_fasta)

## ----scenario4-grid-----------------------------------------------------------
fasta_file <- system.file("extdata", "test_transcripts.fa", package = "RBPEqBind")

res_grid_fasta <- simulateGridFasta(
  fasta_file = fasta_file,
  rbp_models = rbp_models,
  protein_conc_grid = list(
    HH = c(10, 100),
    HL = c(10, 100),
    LH = c(10, 100),
    LL = c(10, 100)
  ),
  rna_conc_grid = c(10),
  k = 5
)
head(res_grid_fasta)

## ----session------------------------------------------------------------------
sessionInfo()