## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 12
)

## ----setup--------------------------------------------------------------------
library(Aerith)
library(dplyr)
library(stringr)
library(ggplot2)

## ----include = FALSE, eval=FALSE----------------------------------------------
# devtools::load_all()
# rmarkdown::render("PSM-annotation-and-visualization.Rmd", output_dir = "../doc/")

## -----------------------------------------------------------------------------
demo_file <- system.file("extdata", "107728.FT2", package = "Aerith")
scan1 <- readOneScanMS2(ftFile = demo_file, 107728)
anno <- annotatePSM(
    scan1$peaks$mz, scan1$peaks$intensity,
    scan1$peaks$charge,
    "HSQVFSTAEDNQSAVTIHVLQGER", 1:2, "C13",
    0.0107, scan1$isolationWindowCenterMZ, 4.0
)
head(anno$ExpectedBYions[anno$ExpectedBYions$matchedIndices != -1, ])
residuePos <- anno$ExpectedBYions$residuePositions[anno$ExpectedBYions$matchedIndices != -1]
table(residuePos)

## -----------------------------------------------------------------------------
set.seed(9527)
p <- plotPSMannotation(
    observedSpect = getRealScanFromList(scan1),
    pep = "HSQVFSTAEDNQSAVTIHVLQGER", Atom = "C13", Prob = 0.01,
    charges = 1:2, isoCenter = 886.65, isoWidth = 4.0,
    ifRemoveNotFoundIon = TRUE
)
p

## -----------------------------------------------------------------------------
a <- getSipBYionSpectra("HSQVFSTAEDNQSAVTIHVLQGER", "C13", 0.01, 1:2)
slot(a, "spectra") <- slot(a, "spectra")[slot(a, "spectra")$MZ < 2000, ]
p <- plot(a)
p <- p + plotSipBYionLabel(a)
demo_file <- system.file("extdata", "107728.FT2", package = "Aerith")
b <- readAllScanMS2(demo_file)
c <- getRealScan(107728, b)
p <- p + plotRealScan(c)
p

## -----------------------------------------------------------------------------
demo_file <- system.file("extdata", "107695.FT1", package = "Aerith")
ft1 <- readOneScanMS1(demo_file, 107695)
precursorScan1 <- getRealScanFromList(ft1)

pep <- "HSQVFSTAEDNQSAVTIHVLQGER"
precursorSP <- getSipPrecursorSpectra(pep, Prob = 0.0107, charges = 3)
slot(precursorSP, "spectra")$Kind <- "Expected"
xlimit <- slot(precursorScan1, "spectra")$MZ > 880 & slot(precursorScan1, "spectra")$MZ < 890
slot(precursorScan1, "spectra") <- slot(precursorScan1, "spectra")[xlimit, ]
slot(precursorScan1, "spectra")$Kind <- "Observed"
maxInt <- max(slot(precursorScan1, "spectra")$Prob)
slot(precursorScan1, "spectra")$Prob <- slot(precursorScan1, "spectra")$Prob / maxInt * 100
p <- plot(precursorSP, linewidth = 0.3) + plotRealScan(precursorScan1, linewidth = 0.3) +
    scale_x_continuous(breaks = seq(880, 890, by = 1)) +
    theme(legend.title = element_blank()) +
    scale_color_manual(values = c("#E7872B", "#F3082F"))
p

## -----------------------------------------------------------------------------
observed_precursor <- getRealScanFromList(ft1)
xlimit <- observed_precursor@spectra$MZ > 880 & observed_precursor@spectra$MZ < 890
observed_precursor@spectra <- observed_precursor@spectra[xlimit, ]

precursor_anno <- annotatePrecursor(
    observed_precursor@spectra$Mass,
    observed_precursor@spectra$Prob,
    observed_precursor@spectra$Charge,
    pepSeq = pep,
    charge = 3,
    Atom = "C13",
    Prob = 0.0107,
    isoCenter = 886.65,
    isoWidth = 4.0
)
matched_precursor <- precursor_anno$ExpectedPrecursorIons[
    precursor_anno$ExpectedPrecursorIons$matchedIndices != -1,
    c("mz", "intensity", "charge", "matchedIndices", "SIPabundances")
]
head(matched_precursor)

plotPrecursorAnnotation(
    observedSpect = observed_precursor,
    pep = pep,
    charge = 3,
    Atom = "C13",
    Prob = 0.0107,
    isoCenter = 886.65,
    isoWidth = 4.0,
    xwidth = 10.0,
    ifRemoveNotFoundIon = TRUE,
    ifAdjustSIPabundance = FALSE
)

## -----------------------------------------------------------------------------
demo_file <- system.file("extdata", "X13_4068_2596_8182.FT2", package = "Aerith")
scan1 <- readAllScanMS2(ftFile = demo_file)[["2596"]]
anno <- annotatePSM(
    scan1$peaks$mz, scan1$peaks$intensity,
    scan1$peaks$charge,
    "HYAHVDCPGHADYVK", 1:2, "C13",
    0.52, scan1$isolationWindowCenterMZ, 5.0
)
head(anno$ExpectedBYions[anno$ExpectedBYions$matchedIndices != -1, ])
residuePos <- anno$ExpectedBYions$residuePositions[anno$ExpectedBYions$matchedIndices != -1]
table(residuePos)

## -----------------------------------------------------------------------------
set.seed(9527)
p <- plotPSMannotation(
    observedSpect = getRealScanFromList(scan1),
    pep = "HYAHVDCPGHADYVK", Atom = "C13", Prob = 0.52,
    charges = 1:2, isoCenter = scan1$isolationWindowCenterMZ, isoWidth = 5.0,
    ifRemoveNotFoundIon = TRUE
)
p

## -----------------------------------------------------------------------------
demo_file <- system.file("extdata", "X13_4068_2596_8182.FT2", package = "Aerith")
ft2 <- readAllScanMS2(demo_file)
a <- getSipBYionSpectra("HYAHVDCPGHADYVK", "C13", 0.52, 1:2)
p <- plot(a)
p <- p + plotSipBYionLabel(a)
c <- getRealScan(2596, ft2)
p <- p + plotRealScan(c)
p

## -----------------------------------------------------------------------------
demo_file <- system.file("extdata", "X13_2559.FT1", package = "Aerith")
ft1 <- readOneScanMS1(demo_file, 2559)
precursorScan1 <- getRealScanFromList(ft1)

pep <- "HYAHVDCPGHADYVK"
precursorSP <- getSipPrecursorSpectra(pep, Prob = 0.5, charges = 3)
slot(precursorSP, "spectra")$Kind <- "Expected"
xlimit <- slot(precursorScan1, "spectra")$MZ > 590 & slot(precursorScan1, "spectra")$MZ < 620
slot(precursorScan1, "spectra") <- slot(precursorScan1, "spectra")[xlimit, ]
slot(precursorScan1, "spectra")$Kind <- "Observed"
maxInt <- max(slot(precursorScan1, "spectra")$Prob)
slot(precursorScan1, "spectra")$Prob <- slot(precursorScan1, "spectra")$Prob / maxInt * 100
p <- plot(precursorSP, linewidth = 0.3) + plotRealScan(precursorScan1, linewidth = 0.3) +
    scale_x_continuous(breaks = seq(590, 620, by = 5)) +
    theme(legend.title = element_blank()) +
    scale_color_manual(values = c("#E7872B", "#F3082F"))
p

## -----------------------------------------------------------------------------
observed_precursor <- getRealScanFromList(ft1)
xlimit <- observed_precursor@spectra$MZ > 590 & observed_precursor@spectra$MZ < 620
observed_precursor@spectra <- observed_precursor@spectra[xlimit, ]

precursor_anno <- annotatePrecursor(
    observed_precursor@spectra$Mass,
    observed_precursor@spectra$Prob,
    observed_precursor@spectra$Charge,
    pepSeq = pep,
    charge = 3,
    Atom = "C13",
    Prob = 0.5,
    isoCenter = 605.0,
    isoWidth = 5.0
)
matched_precursor <- precursor_anno$ExpectedPrecursorIons[
    precursor_anno$ExpectedPrecursorIons$matchedIndices != -1,
    c("mz", "intensity", "charge", "matchedIndices", "SIPabundances")
]
head(matched_precursor)

plotPrecursorAnnotation(
    observedSpect = observed_precursor,
    pep = pep,
    charge = 3,
    Atom = "C13",
    Prob = 0.5,
    isoCenter = 605.0,
    isoWidth = 5.0,
    xwidth = 30.0,
    ifRemoveNotFoundIon = TRUE,
    ifAdjustSIPabundance = FALSE
)

## -----------------------------------------------------------------------------
element <- "C13"
demo_file <- system.file("extdata", "demo.psm.txt", package = "Aerith")
psm <- readPSMtsv(demo_file)
psm <- psm[psm$Filename == "Pan_052322_X13.FT2", ]
psm <- psm[psm$ScanNumber %in% c("4068", "2596", "8182"), ]
demo_file <- system.file("extdata", "X13_4068_2596_8182.FT2", package = "Aerith")
ft2 <- readAllScanMS2(demo_file)
ftFileNames <- psm$Filename
scanNumbers <- psm$ScanNumber
proNames <- psm$ProteinNames
charges <- psm$ParentCharge
pep <- psm$OriginalPeptide
pep <- stringr::str_sub(pep, 2, -2)
pct <- psm$SearchName
pct <- as.numeric(stringr::str_sub(
    stringr::str_split(pct, "_", simplify = TRUE)[, 2], 1, -4
)) / 100 / 1000
realScans <- getRealScans(ft2, scanNumbers)
tmp <- tempdir()
plotPSMs(
    realScans,
    charges,
    element,
    pct,
    BYcharge = 1:2,
    ftFileNames,
    scanNumbers,
    pep,
    proNames,
    path = tmp
)
list.files(tmp, pattern = ".pdf", full.names = TRUE)

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