---
title: "Score function test of PSM"
output:
  rmarkdown::html_document:
    toc: true
    toc_float: true
    theme: united
vignette: >
  %\VignetteIndexEntry{Score-function-test-of-PSM}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 10
)
```

```{r, include = FALSE, eval=FALSE}
library(devtools)
build_vignettes()
rmarkdown::render("Score-function-test-of-PSM.Rmd", output_dir = "../doc/")
```

```{r setup}
library(Aerith)
library(dplyr)
library(stringr)
library(ggplot2)
library(tidyr)
```

The Aerith package delivers isotope-aware peptide-spectrum matching (PSM) scoring that keeps pace with modern stable isotope probing (SIP) workflows. This vignette walks through two representative scenarios—an unlabeled peptide and a peptide labeled at roughly 50% $^{13}$C—to illustrate how Aerith quantifies evidence for isotopic incorporation. Each section begins with the experimental context, introduces the relevant Aerith functions, and explains how to interpret the resulting scores for method development or routine quality control.

## Overview of scoring functions

Aerith exposes multiple scoring strategies that capture complementary information:

- **`scorePSM()`** evaluates observed spectra against theoretical expectations across a user-specified isotopic abundance grid by (`WDP`) function in Sipros.
- **`annotatePSM()`** returns cross-correlation (`XcorrScore`) and multivariate hypergeometric (`MVHscore`) metrics after annotating fragment ions.

These scores are designed to support threshold selection and benchmarking against state-of-the-art SIP-aware search engines. The following sections demonstrate how to configure their key parameters and interpret the resulting trends.

### Unlabeled PSM at 1.07% ^13^C 

We begin with an Aerith example file representing an effectively unlabeled spectrum. The goal is to confirm that Aerith reports maximal scores near natural-abundance $^{13}$C levels.

```{r}
demo_file <- system.file("extdata", "107728.FT2", package = "Aerith")
scan2 <- readOneScanMS2(ftFile = demo_file, 107728)
scan1 <- getRealScanFromList(scan2)
spectra <- slot(scan1, "spectra")
charges <- slot(scan1, "charges")
pep <- "HSQVFSTAEDNQSAVTIHVLQGER"
pep2 <- "[HSQVFSTAEDNQSAVTIHVLQGER]"
```

The `readOneScanMS2()` and `getRealScanFromList()` helpers retrieve the centroided fragment spectrum for scan 107728. Throughout this vignette, peptide sequences are provided both as plain strings (`pep`) for annotation and with bracketed notation (`pep2`) for scoring functions that expect capped termini.

```{r}
isoCenter <- scan2$isolationWindowCenterMZ
anno <- annotatePSM(
    spectra$MZ, spectra$Prob,
    spectra$Charge,
    pep, 1:2, "C13",
    0.0107, isoCenter, 5.0, TRUE
)
WDPscores <- sapply((0:100) / 100, simplify = TRUE, function(prob) {
    return(scorePSM(spectra$MZ,
        realIntensity = spectra$Prob,
        realCharge = spectra$Charge, parentCharge = charges[1],
        pepSeq = pep2, Atom = "C13", Prob = prob
    ))
})
scores <- sapply((0:100) / 100, simplify = TRUE, function(prob) {
    anno <- annotatePSM(
        spectra$MZ, spectra$Prob,
        spectra$Charge,
        pep, 1:2, "C13",
        prob, isoCenter, 4.0, TRUE
    )
    return(c(XcorrScore = anno$XcorrScore, MVHscore = anno$MVHscore))
})
scores <- t(scores)
scores <- cbind(WDPscores, scores)
```

Here we scan the full 0–100% $^{13}$C grid to understand score dynamics. Key arguments are:

- `isoCenter`: the isolation window center in $m/z$, sourced from instrument metadata or PSMs file to anchor precursor annotation.
- `prob`: the hypothesized $^{13}$C proportion. The fixed value (0.0107) approximates natural abundance for `annotatePSM()`, while the loop explores alternative hypotheses.
- `parentCharge`: drawn from `scan1@charges` to ensure theoretical spectra match the observed charge state.

`annotatePSM()` produces both cross-correlation and MVH scores, while `scorePSM()` returns the weighted dot product (WDP) score. Evaluating all three allows users to pick thresholds aligned with their downstream filtering strategy.

```{r}
scores_df <- as.data.frame(scores)
scores_df$abundance <- 0:100

# Reshape the data to long format for faceting
scores_long <- scores_df %>%
    pivot_longer(
        cols = c(WDPscores, XcorrScore, MVHscore),
        names_to = "score_type",
        values_to = "score_value"
    )

# Create a publication-quality plot
ggplot(scores_long, aes(x = abundance, y = score_value)) +
    geom_line() +
    facet_wrap(~score_type, scales = "free_y") +
    labs(
        x = expression(paste(~ {}^{
            13
        }, "C %")),
        y = "Score",
        title = expression(paste("Scores Across", ~ {}^{
            13
        }, "C %"))
    ) +
    theme_minimal() +
    theme(
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        strip.text = element_text(size = 12, face = "bold"),
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "white", colour = "black")
    )
```

The resulting faceted plot summarizes how each score reacts to varying $^{13}$C abundance assumptions. For this near-unlabeled spectrum, all metrics peak close to natural abundance, confirming the absence of substantial label incorporation. The free $y$-scales highlight the relative sensitivity of each metric without forcing a shared dynamic range.

### Labeled PSM at 50% ^13^C 

#### Annotate B and Y ion fragments

We now investigate a spectrum from the Aerith labeled dataset, expecting scores to favor higher $^{13}$C proportions consistent with experimental labeling.

```{r}
demo_file <- system.file("extdata", "X13_4068_2596_8182.FT2", package = "Aerith")
ft2 <- readAllScanMS2(demo_file)
scan2 <- ft2[["2596"]]
scan1 <- getRealScanFromList(scan2)
spectra <- slot(scan1, "spectra")
charges <- slot(scan1, "charges")
pep <- "HYAHVDCPGHADYVK"
pep2 <- "[HYAHVDCPGHADYVK]"
pct <- 0.52
```

The full scan list is loaded with `readAllScanMS2()` so we can reference scan 2596 explicitly. 

```{r}
isoCenter <- scan2$isolationWindowCenterMZ
WDPscores <- sapply((0:100) / 100, simplify = TRUE, function(prob) {
    return(scorePSM(spectra$MZ,
        realIntensity = spectra$Prob,
        realCharge = spectra$Charge, parentCharge = charges[1],
        pepSeq = pep2, Atom = "C13", Prob = prob
    ))
})
scores <- sapply((0:100) / 100, simplify = TRUE, function(prob) {
    anno <- annotatePSM(
        spectra$MZ, spectra$Prob,
        spectra$Charge,
        pep, 1:2, "C13",
        prob, isoCenter, 4.0, TRUE
    )
    return(c(XcorrScore = anno$XcorrScore, MVHscore = anno$MVHscore))
})
scores <- t(scores)
scores <- cbind(WDPscores, scores)
```

As before, we iterate across the full probability grid to visualize the likelihood landscape. Practitioners may narrow the grid around anticipated labeling levels (for example, 30–70%) to accelerate large-scale studies, but the full sweep is informative when benchmarking instrument performance or validating new sample types.

```{r}
scores_df <- as.data.frame(scores)
scores_df$abundance <- 0:100

# Reshape the data to long format for faceting
scores_long <- scores_df %>%
    pivot_longer(
        cols = c(WDPscores, XcorrScore, MVHscore),
        names_to = "score_type",
        values_to = "score_value"
    )

# Create a publication-quality plot
ggplot(scores_long, aes(x = abundance, y = score_value)) +
    geom_line() +
    facet_wrap(~score_type, scales = "free_y") +
    labs(
        x = expression(paste(~ {}^{
            13
        }, "C %")),
        y = "Score",
        title = expression(paste("Scores Across", ~ {}^{
            13
        }, "C %"))
    ) +
    theme_minimal() +
    theme(
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12),
        strip.text = element_text(size = 12, face = "bold"),
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "white", colour = "black")
    )
```

In this labeled scenario, each score exhibits a clear maximum near 50% abundance, mirroring the known experimental condition. Divergence among the score shapes highlights how cross-correlation and hypergeometric statistics can provide complementary evidence; Aerith lets analysts combine or prioritize these metrics depending on downstream classification needs.

### Guidance on parameter selection and interpretation

- **Isolation window (`isoCenter`)**: always source from spectrum metadata to avoid systematic shifts in the theoretical isotope envelope.
- **Probability grid**: dense grids (e.g., 0–1% steps) aid exploratory analyses and method validation, while coarser grids reduce compute cost during high-throughput processing. Align step size with the expected precision of labeling efficiency estimates.
- **Score evaluation**: identify the probability that maximizes each metric. Consistent peaks across WDP, Xcorr, and MVH suggest robust assignment. Discrepancies can indicate spectrum quality issues or mixed precursor populations.
- **Thresholding**: use empirical distributions of these scores across known target and decoy sets to set decision thresholds, leveraging Aerith's compatibility with Percolator and other state-of-the-art post-processors.

### Why Aerith stands out

Aerith integrates isotope-aware scoring, fragment annotation, and downstream export utilities within a single R-native workflow. Compared with generic proteomics search engines, Aerith explicitly models SIP isotope envelopes, offers fine-grained control over labeling hypotheses, and interoperates with modern statistical rescoring pipelines. This alignment with current SIP best practices ensures that researchers can bridge method development and routine analysis without leaving the R ecosystem.

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