---
title: "Scoring simulated SAINT input with the pure-R engine"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Scoring simulated SAINT input with the pure-R engine}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
set.seed(42)
```

This vignette walks through scoring a simulated affinity-purification experiment
with `saintexpress::run_saint()`. No native binary is involved — everything runs
in pure R.

## Simulate a small AP-MS experiment

We define two real baits (`BaitA`, `BaitB`) with two replicates each, plus two
controls. Six preys exist; `BaitA` enriches for `Prey1`/`Prey2` and `BaitB` for
`Prey3`/`Prey4`. The remaining preys are background.

```{r simulate}
simulate_si <- function(seed = 42, mode = c("spc", "int")) {
  mode <- match.arg(mode)
  set.seed(seed)
  preys  <- paste0("Prey", 1:6)
  baits  <- c("BaitA", "BaitA", "BaitB", "BaitB", "Ctrl1", "Ctrl2")
  ips    <- paste0("IP", seq_along(baits))
  cort   <- c("T", "T", "T", "T", "C", "C")

  # One row per IP x prey; recover the bait for each IP.
  grid <- expand.grid(ipId = ips, preyId = preys, stringsAsFactors = FALSE)
  grid$baitId <- baits[match(grid$ipId, ips)]

  # A prey is enriched ("signal") only for its true bait; everything else is background.
  signal <- (grid$baitId == "BaitA" & grid$preyId %in% c("Prey1", "Prey2")) |
            (grid$baitId == "BaitB" & grid$preyId %in% c("Prey3", "Prey4"))

  # Draw every quantity in one vectorized call.
  if (mode == "spc") {
    lambda <- ifelse(signal, ifelse(grid$baitId == "BaitA", 20, 15), 0.5)
    grid$quant <- stats::rpois(nrow(grid), lambda)
  } else {
    mu <- ifelse(signal, ifelse(grid$baitId == "BaitA", 1e7, 5e6), 1e4)
    grid$quant <- stats::rexp(nrow(grid), 1 / mu)
  }

  inter <- grid[grid$quant > 0, c("ipId", "baitId", "preyId", "quant")]
  rownames(inter) <- NULL
  list(
    inter = inter,
    prey  = data.frame(preyId = preys, preyLength = 500L, preyGeneId = preys,
                       stringsAsFactors = FALSE),
    bait  = data.frame(ipId = ips, baitId = baits, CorT = cort,
                       stringsAsFactors = FALSE)
  )
}

si_spc <- simulate_si(mode = "spc")
str(si_spc, max.level = 2)
```

## Validate the input tables

```{r validate}
saintexpress::validate_saint_input(si_spc)
```

## Spectral-count scoring

```{r score-spc}
scores_spc <- saintexpress::run_saint(si_spc, mode = "spc", optimizer = "base")
scores_spc[, c("Bait", "Prey", "AvgP", "BFDR", "SaintScore")]
```

The true interactors (`Prey1`/`Prey2` for `BaitA`, `Prey3`/`Prey4` for `BaitB`)
should land at the top of each bait by `AvgP`:

```{r top-hits-spc}
top_per_bait <- function(df) {
  df <- df[order(df$Bait, -df$AvgP), ]
  do.call(rbind, by(df, df$Bait, head, 2))
}
top_per_bait(scores_spc)[, c("Bait", "Prey", "AvgP", "BFDR")]
```

## Volcano plot

The `FoldChange` column summarizes enrichment over the estimated background.
We plot it as an effect-size coordinate (`EFC`, here `log2(FoldChange)`) against
the `BFDR` evidence coordinate.

```{r volcano-spc, fig.width = 6, fig.height = 4}
plot_volcano <- function(scores, title) {
  volcano_data <- scores
  volcano_data$EFC <- log2(volcano_data$FoldChange)
  volcano_data$neg_log10_bfdr <- -log10(pmax(volcano_data$BFDR, 1e-16))
  volcano_data$is_hit <- volcano_data$BFDR <= 0.05 & volcano_data$EFC >= 1

  plot(
    volcano_data$EFC,
    volcano_data$neg_log10_bfdr,
    pch = 19,
    col = ifelse(volcano_data$is_hit, "#B23A48", "#3D6F99"),
    xlab = "EFC (log2 fold change)",
    ylab = "-log10(BFDR)",
    main = title
  )
  abline(h = -log10(0.05), lty = 2, col = "grey40")
  abline(v = 1, lty = 2, col = "grey40")
  text(
    volcano_data$EFC,
    volcano_data$neg_log10_bfdr,
    labels = paste(volcano_data$Bait, volcano_data$Prey, sep = ":"),
    pos = 3,
    cex = 0.7
  )
}

plot_volcano(scores_spc, "SAINTexpress spectral-count volcano plot")
```

## Intensity scoring

The same simulator, with `mode = "int"`, produces continuous abundances.

```{r score-int}
si_int <- simulate_si(mode = "int")
scores_int <- saintexpress::run_saint(si_int, mode = "int", optimizer = "base")
top_per_bait(scores_int)[, c("Bait", "Prey", "AvgP", "BFDR")]
```

```{r volcano-int, fig.width = 6, fig.height = 4}
plot_volcano(scores_int, "SAINTexpress intensity volcano plot")
```

## See also

For running the original native SAINTexpress C++ binary on the same input,
see the companion package
[`saintexpressbin`](https://github.com/prolfqua/saintexpressbin).

## Session information

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