Quality control for the target decoy approach for peptide identification

Ghent University

1 load libraries

suppressPackageStartupMessages({
library(tidyverse)
library(TargetDecoy)
library(RCurl)
})

2 Download data from Zenodo

# Download data from zenodo
options(timeout=300)
url <- "https://zenodo.org/record/7308022/files/search-results.zip?download=1"
destFile <- "searchResults.zip"
if (!file.exists(destFile)) download.file(url, destFile)
unzip(destFile, exdir = "./data", overwrite = TRUE)

3 Pyrococcus

3.1 Import Pyrococcus Data in R

allTsvFiles <- list.files(
  path = "data", 
  pattern = ".tsv$",
  full.names = TRUE)
msgfFiles <- allTsvFiles[grepl("msgf",allTsvFiles)&grepl("PXD001077",allTsvFiles)]
dfsmsgf <- lapply(msgfFiles, read_tsv)

## Rows: 136982 Columns: 34
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr  (8): peptidoform, spectrum_id, run, protein_list, source, provenance:mz...
## dbl (22): score, precursor_mz, retention_time, rank, meta:calculatedMassToCh...
## lgl  (4): collection, is_decoy, qvalue, pep
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 126277 Columns: 34
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr  (8): peptidoform, spectrum_id, run, protein_list, source, provenance:mz...
## dbl (22): score, precursor_mz, retention_time, rank, meta:calculatedMassToCh...
## lgl  (4): collection, is_decoy, qvalue, pep
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

crapEntries <- scan("db-gpm-crap-entries.txt", what = "character")
dfsmsgf <- lapply(dfsmsgf, function(db) 
  db[rowSums(sapply(crapEntries, grepl,fixed=TRUE,x = db$protein_list)) == 0,])

tandemFiles <- allTsvFiles[grepl("xtandem",allTsvFiles)&grepl("PXD001077",allTsvFiles)]
dfsTandem <- lapply(tandemFiles, read_tsv)

## Rows: 12586 Columns: 18
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (5): peptidoform, spectrum_id, protein_list, source, provenance:xtandem_...
## dbl (7): score, precursor_mz, retention_time, provenance:xtandem_id, meta:xt...
## lgl (6): run, collection, is_decoy, qvalue, pep, rank
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 13949 Columns: 18
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (5): peptidoform, spectrum_id, protein_list, source, provenance:xtandem_...
## dbl (7): score, precursor_mz, retention_time, provenance:xtandem_id, meta:xt...
## lgl (6): run, collection, is_decoy, qvalue, pep, rank
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

dfsTandem <- lapply(dfsTandem, function(db) 
  db[rowSums(sapply(crapEntries, grepl,fixed=TRUE,x = db$protein_list)) == 0,])

3.2 FDR function

Calculate FDR on tibble with search results. It is assumed that higher scores are better.

# FDR function
fdrDb <- function(db)
{
  db <- db %>% 
    arrange(desc(score)) %>% 
    mutate(FP = cumsum(is_decoy),
           FDR = cumsum(is_decoy)/cumsum(!is_decoy))
  FDR <- db$FDR
  FDRmin <- FDR[length(FDR)]
  for (j in (length(FDR)-1):1)
  {
    if (FDR[j] < FDRmin)
      FDRmin <- FDR[j] else
      FDR[j] <- FDRmin   
  }
  db$FDR <- FDR
  db <- db %>% 
    arrange(spectrum_id) 
  return(db)
}

3.3 Generate plots for MSGF+

Histograms and P-P plots for rank1, rank 2 and search against human DB.

# Generate MSGF+ plots pyrococcus
histsMsgf <- lapply(
  dfsmsgf,
  function(db, score, decoy, log10) 
    evalTargetDecoysHist(
      db %>% filter(rank==1), 
      decoy, 
      score, 
      log10) + 
    xlab("Score") + 
    ggtitle(NULL) +   
    geom_histogram(
      bins = 50, 
      position = "identity",
      alpha = .9), 
  score = "score", 
  decoy = "is_decoy", 
  log10 = TRUE)
histsMsgfR2 <- lapply(
  dfsmsgf,
  function(db, score, decoy, log10) 
    evalTargetDecoysHist(db %>% filter(rank==2), 
                         decoy, 
                         score, 
                         log10) + 
    xlab("Score") + 
    ggtitle(NULL) +  
    geom_histogram(bins=50, position="identity",alpha=.9), 
  score = "score", 
  decoy = "is_decoy", 
  log10 = TRUE)
ppPlotsMsgf <- lapply(dfsmsgf,
                      function(db, score, decoy, log10) 
                        evalTargetDecoysPPPlot(db %>% filter(rank==1), 
                                               decoy, 
                                               score, 
                                               log10) + 
                        ggtitle(NULL) +
                        xlab("Fd") +
                        ylab("Ft"), 
                      score = "score", 
                      decoy = "is_decoy", 
                      log10=TRUE)

pyroId <- which(
  grepl(
    pattern = "swissprot",
    msgfFiles,
    fixed = TRUE) & 
    grepl(
      pattern = "pfuriosus",
      msgfFiles,
      fixed = TRUE)
  ) 

dfHlp <- dfsmsgf[[pyroId]]  %>% 
  filter(rank==1) %>% 
  mutate(score = -log10(score)) %>%
  fdrDb()
thresh <- dfHlp %>%  
  filter(FDR < 0.01) %>% 
  pull(score) %>% 
  min
nSig <- dfHlp %>% filter(FDR < 0.01) %>% pull(is_decoy) %>% `!` %>% sum

humanId <- which(
  grepl(pattern="swissprot",msgfFiles,fixed = TRUE) & grepl(pattern="hsapiens",msgfFiles,fixed = TRUE)
  ) 

figMsgfSwissHistsPyroR1R2_Human <- gridExtra::grid.arrange(
      histsMsgf[[pyroId]] +
        geom_histogram(
          bins = 50, 
          position = "identity",
          alpha = .9) +
         ggtitle(paste(nSig,"target PSMs at 1% FDR")) +
        annotate("rect", 
                 xmin = thresh, 
                 xmax = 35, 
                 ymin = -10, 
                 ymax = 1600, 
                 alpha = .2) + 
      annotate(geom = "text", 
               x = thresh+1, 
               y = 1500, 
               label = "x > t['1% FDR']",
              color = "black",
              hjust = 0 ,
              parse = TRUE) +
        annotate(geom = "rect",
                 xmin = thresh,
                 xmax = thresh,
                 ymin = 0,
                 ymax = 1600,
                 color = "red") +
         theme(legend.position = c(0.75,0.75)),
    histsMsgfR2[[pyroId]] + 
         geom_histogram(
           bins = 50, 
           position = "identity",
           alpha = .9) +
         ggtitle("Rank 2 PSMs") +
         theme(legend.position = c(0.75,0.75)),
    histsMsgf[[humanId]] + 
         geom_histogram(
           bins = 50, 
           position = "identity",
           alpha = .9) +
         ggtitle("Spectra matched to H. sapiens") +
         theme(legend.position = c(0.75,0.75)),
    ncol=3)

ggsave( 
  "./figs/figMsgfSwissHistsPyroR1R2_Human.png", 
  plot = figMsgfSwissHistsPyroR1R2_Human,
  device = "png", 
  width = 11.7,
  height = 3.9)

ggsave( 
  "./tiffs/figMsgfSwissHistsPyroR1R2_Human.tiff", 
  plot = figMsgfSwissHistsPyroR1R2_Human,
  device = "tiff", 
  width = 11.7,
  height = 3.9)

## Warning in grDevices::dev.off(): unable to open TIFF file './tiffs/
## figMsgfSwissHistsPyroR1R2_Human.tiff'

figMsgfSwissPPplotsPyro_Human <- 
    gridExtra::grid.arrange(
      ppPlotsMsgf[[humanId]] + 
         ggtitle(NULL),
       ppPlotsMsgf[[pyroId]] + 
         ggtitle(NULL),
       ncol=2)

ggsave( 
  "./figs/figMsgfSwissPPplotsPyro_Human.png", 
  plot = figMsgfSwissPPplotsPyro_Human,
  device = "png", 
  width = 7.8, 
  height = 3.9)

ggsave( 
  "./tiffs/figMsgfSwissPPplotsPyro_Human.tiff", 
  plot = figMsgfSwissPPplotsPyro_Human,
  device = "tiff", 
  width = 7.8, 
  height = 3.9)

## Warning in grDevices::dev.off(): unable to open TIFF file './tiffs/
## figMsgfSwissPPplotsPyro_Human.tiff'

3.4 X-tandem searches

3.4.1 Preprocess tandem data

Convert spectrum id into double for sorting + add FDR

# Preprocess Tandem pyrococcus results
for (i in 1:length(dfsTandem))
{
  
  dfsTandem[[i]] <- dfsTandem[[i]] %>%
    mutate(spectrum_id_orig = spectrum_id,
           spectrum_id = sapply(
             spectrum_id %>% strsplit(split=" "), 
             function(x) substr(x[3],6,1000)) %>% as.double
    )
}
dfsTandem <- lapply(dfsTandem, fdrDb)

3.4.2 Plots for search with and without refinement.

# plots xtandem Results
histsTandem <- lapply(dfsTandem,
                      evalTargetDecoysHist, 
                      score = "score", 
                      decoy = "is_decoy", 
                      log10 = FALSE)
ppPlotsTandem <- lapply(dfsTandem,
                        evalTargetDecoysPPPlot, 
                        score = "score", 
                        decoy = "is_decoy", 
                        log10 = FALSE)

noRefineId <- which(grepl("no-refine",tandemFiles))
refineId <- which(!grepl("no-refine",tandemFiles))


figTandemRefineSwissHistPP <- 
    gridExtra::grid.arrange(
       histsTandem[[refineId]] + 
         geom_histogram(
           bins = 50, 
           position = "identity",
           alpha=.9) +
         ylim(0, 1500) +
         xlim(-6.2, 35) +
         ggtitle(NULL) +
         theme(legend.position = c(0.75,0.75)),
       ppPlotsTandem[[refineId]] + 
         ggtitle(NULL) +
         ylim(0,.6), 
       ncol = 2)

## Warning: Removed 4 rows containing missing values (`geom_bar()`).
## Removed 4 rows containing missing values (`geom_bar()`).

ggsave( 
  "./figs/figTandemRefineSwissHistPP.png", 
  plot =  figTandemRefineSwissHistPP,
  device = "png", 
  width = 7.8, 
  height = 3.9)

ggsave( 
  "./tiffs/figTandemRefineSwissHistPP.tiff", 
  plot =  figTandemRefineSwissHistPP,
  device = "tiff", 
  width = 7.8, 
  height = 3.9)

## Warning in grDevices::dev.off(): unable to open TIFF file './tiffs/
## figTandemRefineSwissHistPP.tiff'

figTandemNoRefineSwissHistPP <- 
    gridExtra::grid.arrange(
       histsTandem[[noRefineId]] + 
        geom_histogram(
          bins = 50, 
          position = "identity",
          alpha = .9) +
         ylim(0, 1500) + 
         xlim(-6.2, 35) +
         ggtitle(NULL) + 
         theme(legend.position = c(0.75,0.75)),
       ppPlotsTandem[[noRefineId]] + 
         ggtitle(NULL) +
         ylim(0,.6), 
       ncol=2)

## Warning: Removed 4 rows containing missing values (`geom_bar()`).
## Removed 4 rows containing missing values (`geom_bar()`).

ggsave(
  "./figs/figTandemNoRefineSwissHistPP.png", 
  plot = figTandemNoRefineSwissHistPP,
  device = "png", 
  width = 7.8, 
  height = 3.9)

ggsave(
  "./tiffs/figTandemNoRefineSwissHistPP.tiff", 
  plot = figTandemNoRefineSwissHistPP,
  device = "tiff", 
  width = 7.8, 
  height = 3.9)

## Warning in grDevices::dev.off(): unable to open TIFF file './tiffs/
## figTandemNoRefineSwissHistPP.tiff'

3.4.3 Plot to stratify refinement searched peptides

Calculate 1% FDR threshold for plot

thresh <- dfsTandem[[refineId]] %>%  
  filter(FDR < 0.01) %>% 
  pull(score) %>% 
  min

Identify type of PSM: same match, swiched match or new spectrum upon refinement.

dfsTandem[[refineId]]$refine <- "switched"
dfsTandem[[refineId]]$inNoRefine <- dfsTandem[[refineId]]$spectrum_id %in% 
  dfsTandem[[noRefineId]]$spectrum_id

dfsTandem[[refineId]]$refine[!dfsTandem[[refineId]]$inNoRefine] <- "new"

dfRefineCommon <- dfsTandem[[refineId]] %>% 
  filter(inNoRefine)

dfsTandem[[noRefineId]]$inRefine <- dfsTandem[[noRefineId]]$spectrum_id %in% 
  dfsTandem[[refineId]]$spectrum_id
dfNoRefineCommon <- dfsTandem[[noRefineId]] %>% 
  filter(inRefine)
dfRefineCommon$refine[dfRefineCommon$peptidoform == dfNoRefineCommon$peptidoform] <-"same"

dfsTandem[[refineId]]$refine[dfsTandem[[refineId]]$inNoRefine] <- dfRefineCommon$refine
dfsTandem[[refineId]]$refineTD <- paste(
  dfsTandem[[refineId]]$refine,
  ifelse(
    dfsTandem[[refineId]]$is_decoy,
    "decoy",
    "target")
  ) %>% 
  factor(levels = c("same target", 
                    "same decoy",
                    "new target",
                    "switched target",
                    "new decoy",
                    "switched decoy")
         )

Construct figure with histograms stratefied according to PSM type.

figsTandemSwissRefinementHistPsmType <- dfsTandem[[refineId]] %>% 
  ggplot(aes(x = score,
             fill = refineTD, 
             col = I("black"))) +
  geom_histogram(bins = 50, 
                 position = "identity",
                 alpha = .9) +   
  scale_fill_manual(
    "values" = c(`same target` = "#009900",
                 `same decoy` = "#FF9900", 
                 `new target` = "#5df542", 
                 `new decoy`="#c03600",
                 `switched target`="#42f5c8",
                 `switched decoy` = "#f54281")
    ) +
  theme_bw() + 
  theme(plot.title = element_text(size = rel(1.5)), 
        axis.title = element_text(size = rel(1.2)), 
        axis.text = element_text(size = rel(1.2)), 
        axis.title.y = element_text(angle = 0))

psmType <- dfsTandem[[refineId]] %>% 
  pull(refineTD) %>% 
  table
pi0 <- sum(dfsTandem[[refineId]]$is_decoy)/ sum(!dfsTandem[[refineId]]$is_decoy)
piSame <- psmType["same decoy"]/sum(!dfsTandem[[refineId]]$is_decoy)
piSameNew <- (psmType["same decoy"]+psmType["new decoy"])/sum(!dfsTandem[[refineId]]$is_decoy)
piSameNewSwitch <- (psmType["same decoy"]+psmType["new decoy"] + psmType["switched decoy"])/sum(!dfsTandem[[refineId]]$is_decoy)
piSameNewSwitchNewTarget <- (psmType["same decoy"]+psmType["new decoy"] + psmType["switched decoy"] + psmType["new target"])/sum(!dfsTandem[[refineId]]$is_decoy)


figTandemRefineSwissPSMtypeHistPP <- 
    gridExtra::grid.arrange(
         figsTandemSwissRefinementHistPsmType +
           ggtitle(NULL) +
           ylab(NULL) +
           xlim(-6.2, 35) +
           theme(
             legend.position = c(0.73,0.68),
             legend.title = element_blank()) +
           annotate("rect", 
                    xmin = thresh,
                    xmax = 35, 
                    ymin = -10, 
                    ymax = 1500, 
                    alpha = .2) + 
      annotate(geom = "text", 
               x = thresh+1, 
               y = 1400, 
               label = "x > t['1% FDR']",
               color = "black",
               hjust = 0 ,
               parse = TRUE) +
        annotate(geom = "rect",
                 xmin = thresh,
                 xmax = thresh,
                 ymin = 0,
                 ymax = 1500,
                 color = "red"), 
         ppPlotsTandem[[refineId]] + 
           geom_polygon(
             aes(x = x, y = y),
             data = data.frame(
               x = c(0, 1, 1),
               y = c(0, 0, piSame)),
             fill = "#FF9900",
             inherit.aes = FALSE) +
           geom_polygon(
             aes(x = x,y = y),
             data = data.frame(
               x = c(0, 1, 1),
               y = c(0, piSame, piSameNew)),
             fill = "#c03600",
             inherit.aes = FALSE) +
           geom_polygon(
             aes(x = x,y = y),
             data = data.frame(
               x = c(0, 1, 1),
               y = c(0, piSameNew, piSameNewSwitch)),
             fill = "#f54281",
             inherit.aes = FALSE) +
           geom_polygon(
             aes(x = x,y = y),
             data = data.frame(
               x = c(0, 1, 1),
               y = c(0, piSameNewSwitch, piSameNewSwitchNewTarget)),
             fill = "#5df542",
             inherit.aes = FALSE) +
           geom_abline(slope = pi0) + 
           geom_abline(slope = piSameNewSwitchNewTarget) + 
           ylim(0, 0.6) + 
           ggtitle(NULL), 
      ncol = 2)

## Warning: Removed 12 rows containing missing values (`geom_bar()`).

ggsave("./figs/figTandemRefineSwissPSMtypeHistPP.png", 
       plot = figTandemRefineSwissPSMtypeHistPP,
       device = "png", 
       width = 7.8, 
       height = 3.9)

ggsave("./tiffs/figTandemRefineSwissPSMtypeHistPP.tiff", 
       plot = figTandemRefineSwissPSMtypeHistPP,
       device = "tiff", 
       width = 7.8, 
       height = 3.9)

## Warning in grDevices::dev.off(): unable to open TIFF file './tiffs/
## figTandemRefineSwissPSMtypeHistPP.tiff'

4 Human sample

Make diagnostic plots for rank 1 and rank 2.

# P-P plot for human sample 
path2File <- "./data/PXD028735-LFQ_Orbitrap_DDA_Human_01_uniprot-hsapiens-canonical-isoforms-crap_msgfplus-phospho.msgf.mzid.tsv"
df <- read_tsv(path2File)

## Rows: 1058347 Columns: 21
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (8): peptidoform, spectrum_id, run, protein_list, source, provenance:mzi...
## dbl (9): score, precursor_mz, retention_time, rank, meta:calculatedMassToCha...
## lgl (4): collection, is_decoy, qvalue, pep
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

df <- df[
  rowSums(
    sapply(crapEntries, 
           grepl,
           fixed = TRUE,
           x = df$protein_list)
    ) == 0,]



figHumanMsgfPlus <- gridExtra::grid.arrange(
    evalTargetDecoysHist(
      df %>% filter(rank==1),
      "is_decoy",
      "score",
      TRUE) + 
      ggtitle(NULL) +
      geom_histogram(bins = 50, 
                     alpha = .9, 
                     position = "identity") +
      theme(
             legend.position = c(0.75,0.75),
             ),
    evalTargetDecoysPPPlot(
      df %>% filter(rank==1),
      "is_decoy", 
      "score",
      TRUE) + 
      ggtitle(NULL),
    ncol = 2
  )

ggsave(
  filename = "./figs/figHumanMsgfPlus.png",
  plot = figHumanMsgfPlus,
  device = "png", 
  width = 7.8, 
  height = 3.9)

ggsave(
  filename = "./tiffs/figHumanMsgfPlus.tiff",
  plot = figHumanMsgfPlus,
  device = "tiff", 
  width = 7.8, 
  height = 3.9)

## Warning in grDevices::dev.off(): unable to open TIFF file './tiffs/
## figHumanMsgfPlus.tiff'

figHumanMsgfPlusR2 <- 
    gridExtra::grid.arrange(
     evalTargetDecoysHist(
      df %>% filter(rank==2),
      "is_decoy",
      "score",
      TRUE) + 
      ggtitle(NULL) +
      geom_histogram(bins = 50, 
                     alpha = .9,
        position = "identity") +
      theme(legend.position = c(0.75,0.75)),
    evalTargetDecoysPPPlot(df %>% filter(rank==2),
                           "is_decoy", 
                           "score",
                           TRUE) + 
      ggtitle(NULL),
    ncol = 2
  )

ggsave(
  filename = "./figs/figHumanMsgfPlusR2.png",
  plot = figHumanMsgfPlusR2,
  device = "png", 
  width = 7.8, 
  height = 3.9)

ggsave(
  filename = "./tiffs/figHumanMsgfPlusR2.tiff",
  plot = figHumanMsgfPlusR2,
  device = "tiff", 
  width = 7.8, 
  height = 3.9)

## Warning in grDevices::dev.off(): unable to open TIFF file './tiffs/
## figHumanMsgfPlusR2.tiff'

5 Immunopeptidomics

Make diagnostic plots

# plot immunopeptidomics plot
path2File <- "./data/immunopeptidomics_msms_IAA.txt"
df <- read_tsv(path2File)

## Rows: 7622632 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr  (8): peptide, run, protein_list, source, provenance:msms_filename, meta...
## dbl (10): spectrum_id, score, precursor_mz, retention_time, meta:Delta Score...
## lgl  (4): collection, is_decoy, qvalue, pep
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

figPeptidomics <- 
    gridExtra::grid.arrange(
    evalTargetDecoysHist(df,
                         "is_decoy", 
                         "score",
                         FALSE) + 
      geom_histogram(bins = 50,
                     alpha = .9, 
                     position = "identity") +
      ggtitle(NULL) +
      theme(legend.position = c(0.3,0.7)),
    evalTargetDecoysPPPlot(df,
                           "is_decoy", 
                           "score",
                           FALSE) + 
      xlab("Fd") + 
      ylab("Ft") +
      ggtitle(NULL),
    ncol = 2
  )

ggsave(
  filename = "./figs/figPeptidomics.png",
  plot = figPeptidomics,
  device = "png", 
  width = 7.8, 
  height = 3.9)

ggsave(
  filename = "./tiffs/figPeptidomics.tiff",
  plot = figPeptidomics,
  device = "tiff", 
  width = 7.8, 
  height = 3.9)

## Warning in grDevices::dev.off(): unable to open TIFF file './tiffs/
## figPeptidomics.tiff'

6 Joint P-Plot

Put searches of MSGF+ and XTandem with refinement on one plot.

# Joint PP-plot
msgfId <- which(
  grepl(
    pattern = "swissprot",
    msgfFiles,
    fixed = TRUE) & 
    grepl(pattern = "pfuriosus",
          msgfFiles,
          fixed = TRUE)
  ) 

h <- TargetDecoy:::processObjects(
  list(dfsTandem[[refineId]],
       dfsmsgf[[msgfId]] %>% filter(rank==1)),
  decoy = c("is_decoy", "is_decoy"),
  score = c("score", "score"),
  log10=c(FALSE, TRUE)
  ) %>% 
  TargetDecoy:::ppScoresData()

jointPlot <-  h$df %>%
    ggplot(aes(Fdp, z, color = id %>% as.factor)) + 
    geom_point() +
    geom_abline(slope = 0) + 
    theme_bw() + 
    theme(legend.title = element_blank(),
          legend.position = c(0.2,.85)) +
    ylab("Ft-pi0") +
    scale_color_discrete(labels=c("X!Tandem","MSGF+")) +
    xlab("Fd")

jointPlot

ggsave(
  file = "./figs/figTandemMsGfSwissCombindedPP.png",
  plot = jointPlot,
  device = "png", 
  width = 3.9, 
  height = 3.9)

ggsave(
  file = "./tiffs/figTandemMsGfSwissCombindedPP.tiff",
  plot = jointPlot,
  device = "tiff", 
  width = 3.9, 
  height = 3.9)

## Warning in grDevices::dev.off(): unable to open TIFF file './tiffs/
## figTandemMsGfSwissCombindedPP.tiff'

7 Graphical abstract

histsTandem25 <- lapply(dfsTandem,
                        evalTargetDecoysHist, 
                        score="score", 
                        decoy="is_decoy", 
                        log10 = FALSE,
                        nBins = 25)

histsMsgf25 <- lapply(dfsmsgf,
                      function(db, score, decoy, log10) 
                        evalTargetDecoysHist(
                          db %>% filter(rank==1), 
                          decoy, 
                          score, 
                          log10) + 
                        xlab("Score") + 
                        ggtitle(NULL) +   
                        geom_histogram(bins = 25, 
                                       position = "identity",
                                       alpha = .9), 
                      score = "score", 
                      decoy = "is_decoy", 
                      log10 = TRUE)

  set.seed(15123)
  df <- data.frame(
    xtheo = c(
      sample(
        seq(0.5, 2.1, length = 30), 10),
      rep(NA,5)
      ),
    y0 = rep(0.8, 15), 
    ytheo =rep(2.2, 15)
    ) %>% 
  mutate(
      xreal = c(sample(xtheo,8), 
                sample(seq(0.7, 2.1, length = 20),
                       7)
                ),
      yreal = c(runif(15, .1, 1.5))
  )
theoPlot <- df %>%
  ggplot() + 
  annotate("rect",
           xmin = 0,
           xmax = 2,
           ymin = 1,
           ymax = 3, 
           col = "black",
           fill = "white") + 
  annotate("rect",
           xmin = 0.2,
           xmax = 2.2,
           ymin = .8,
           ymax = 2.8, 
           col = "black", 
           fill = "white") +
    annotate("rect",
             xmin = 0.4,
             xmax = 2.4,
             ymin = .6,
             ymax = 2.6, 
             col = "black",
             fill = "white") +
     annotate("rect",
              xmin = 0,
              xmax = 2.0,
              ymin = -1.6,
              ymax = 0.4, 
              col = "black",
              fill = "white") +
   annotate("rect",
            xmin = 0.2,
            xmax = 2.2,
            ymin = -1.8,
            ymax = 0.2, 
            col = "black",
            fill = "white") +
   annotate("rect",
            xmin = 0.4,
            xmax = 2.4,
            ymin = -2,
            ymax = 0, 
            col = "black",
            fill = "white") +
  geom_segment(aes(x = xtheo,
                   xend = xtheo,
                   y = y0,
                   yend = ytheo)) +
  annotate("text",
           x = 1, 
           y=3.2, 
           label="Database") + 
  annotate("text",
           x = 1.4, 
           y=-2.2, 
           label = "MS2 spectra") +
  geom_segment(aes(
    x = xreal,
    xend = xreal,
    y = y0-2.6,
    yend = yreal-1.8)) +
  theme(axis.line = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position = "none",
        panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.background = element_blank()) +
    annotate("segment", 
             x = 2.8,
             xend = 3.8, 
             y = 0.3,
             yend = 0.3, 
             arrow = arrow()) +
    annotate("text", 
             x = 3.3, 
             y = 1.3, 
             label = "Search\n Engine") +
    annotate("text", 
             x = 3.3, 
             y = -.5, 
             label = "TDA") 
    
tocPlot <- gridExtra::grid.arrange(
    theoPlot,
    histsMsgf25[[pyroId]] +
        geom_histogram(
          bins = 25, 
          position = "identity",
          alpha = .9) +
        ggtitle("Valid TDA") +
        theme(
            legend.position = "none", 
            plot.title = element_text(hjust = 0.5), 
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank()
            ),
    ppPlotsMsgf[[pyroId]] + 
        ggtitle(NULL) + 
        xlab("Fd") + 
        ylab("Ft") +
        ylim(0, .6) + 
        scale_x_continuous(breaks = c(0, 0.5, 1.0)) +
        theme(
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank()
            ),
    histsTandem25[[refineId]] + 
         geom_histogram(bins = 25, 
                        position = "identity",
                        alpha = .9) +
        scale_fill_manual(
            values = c(`FALSE` = "#009900", 
                       `TRUE` = "#FF9900"),
            labels=c("Target", "Decoy")) +
        xlim(-6.2,35) +
        ggtitle("Invalid TDA") +
        theme(
            legend.position = c(0.75,0.6),
            legend.title = element_blank(), 
            plot.title = element_text(hjust = 0.5), 
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank()),
    ppPlotsTandem[[refineId]] + 
        ggtitle(NULL) +
        ylim(0, .6) + 
        xlab("Fd") + 
        ylab("Ft") + 
        scale_x_continuous(breaks = c(0, 0.5, 1.0)) +
        theme(panel.grid.major = element_blank(), 
              panel.grid.minor = element_blank()),
    nrow = 2,
    ncol = 3, 
    layout_matrix = rbind(c(1, 2, 4),
                          c(1, 3, 5))
    )

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

## Warning: Removed 5 rows containing missing values (`geom_segment()`).

## Warning: Removed 4 rows containing missing values (`geom_segment()`).

## Warning: Removed 4 rows containing missing values (`geom_bar()`).
## Removed 4 rows containing missing values (`geom_bar()`).

tocPlot

## TableGrob (2 x 3) "arrange": 5 grobs
##   z     cells    name           grob
## 1 1 (1-2,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,2-2) arrange gtable[layout]
## 4 4 (1-1,3-3) arrange gtable[layout]
## 5 5 (2-2,3-3) arrange gtable[layout]

ggsave("./figs/figGraphicalAbstract.png",
       device = "png",
       plot = tocPlot, 
       width = 7,
       height = 3.9)

ggsave("./tiffs/figGraphicalAbstract.tiff",
       device ="tiff",
       plot = tocPlot, 
       width = 7,
       height = 3.9)

## Warning in grDevices::dev.off(): unable to open TIFF file './tiffs/
## figGraphicalAbstract.tiff'

8 Clean up files

f <- list.files("./data", recursive = TRUE, full.names = TRUE)
file.remove(f)

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

file.remove("data")

## [1] TRUE

file.remove(destFile)

## [1] TRUE

9 Session Info

sessionInfo()

## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## attached base packages:
## [1] stats     graphics  grDevices datasets  utils     methods   base     
## 
## other attached packages:
##  [1] RCurl_1.98-1.9    TargetDecoy_1.4.0 forcats_0.5.2     stringr_1.4.1    
##  [5] dplyr_1.0.10      purrr_0.3.5       readr_2.1.3       tidyr_1.2.1      
##  [9] tibble_3.1.8      ggplot2_3.4.0     tidyverse_1.3.2  
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-7        ProtGenerics_1.30.0 fs_1.5.2           
##  [4] bit64_4.0.5         lubridate_1.9.0     doParallel_1.0.17  
##  [7] httr_1.4.4          tools_4.2.2         backports_1.4.1    
## [10] bslib_0.4.1         utf8_1.2.2          R6_2.5.1           
## [13] DBI_1.1.3           BiocGenerics_0.44.0 colorspace_2.0-3   
## [16] withr_2.5.0         gridExtra_2.3       tidyselect_1.2.0   
## [19] bit_4.0.4           compiler_4.2.2      cli_3.4.1          
## [22] rvest_1.0.3         Biobase_2.58.0      xml2_1.3.3         
## [25] labeling_0.4.2      sass_0.4.2          scales_1.2.1       
## [28] digest_0.6.30       rmarkdown_2.18      pkgconfig_2.0.3    
## [31] htmltools_0.5.3     highr_0.9           dbplyr_2.2.1       
## [34] fastmap_1.1.0       rlang_1.0.6         readxl_1.4.1       
## [37] shiny_1.7.3         farver_2.1.1        jquerylib_0.1.4    
## [40] generics_0.1.3      jsonlite_1.8.3      vroom_1.6.0        
## [43] mzID_1.36.0         googlesheets4_1.0.1 magrittr_2.0.3     
## [46] Rcpp_1.0.9          munsell_0.5.0       fansi_1.0.3        
## [49] lifecycle_1.0.3     stringi_1.7.8       yaml_2.3.6         
## [52] plyr_1.8.7          grid_4.2.2          parallel_4.2.2     
## [55] promises_1.2.0.1    crayon_1.5.2        miniUI_0.1.1.1     
## [58] haven_2.5.1         hms_1.1.2           mzR_2.32.0         
## [61] knitr_1.40          pillar_1.8.1        codetools_0.2-18   
## [64] reprex_2.0.2        XML_3.99-0.12       glue_1.6.2         
## [67] evaluate_0.18       renv_0.16.0         modelr_0.1.10      
## [70] vctrs_0.5.0         tzdb_0.3.0          httpuv_1.6.6       
## [73] foreach_1.5.2       cellranger_1.1.0    gtable_0.3.1       
## [76] assertthat_0.2.1    cachem_1.0.6        xfun_0.34          
## [79] mime_0.12           xtable_1.8-4        broom_1.0.1        
## [82] later_1.3.0         ncdf4_1.19          googledrive_2.0.0  
## [85] gargle_1.2.1        iterators_1.0.14    timechange_0.1.1   
## [88] ellipsis_0.3.2