---
tags: Source
---
# 1.1 Peakpicking with IPO

#install IPO package
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("IPO")
#browseVignettes("IPO")

#load libraries-----------------------------------------------
library(IPO)
library(xcms) # main algorithm for peak picking
# Carl brunius has written a script of how to sort out the data, eg function "getFiles" CMSITools: https://gitlab.com/CarlBrunius/CMSITools
#library(CMSITools,lib.loc = "/home/elise/R/x86_64-pc-linux-gnu-library/3.6/") #please check where you put your library
library(tidyverse) # for function %>%
library(stringi) 
library(StatTools) # For the imputation function, 'mvImp', script written by Carl Brunius, https://gitlab.com/CarlBrunius/StatTools
library(BiocParallel) # parallell computing
register(bpstart(SnowParam(12))) #for parallel computing in Windows
#register(bpstart(MulticoreParam(12))) #for parallel computing in BIANCA
#source(/tmp/src/getFilesMAX.R) #to source getFilesMAX.R function. This function reads the samples type embedded in the file names. otherwise, prepare a table with the sample groups, file names, and injection orders.

#Input MS data------------------------------------------------------------------------------
LCMS_data <- getFiles("/castor/project/proj/Elise/Negative mzML/") %>% cleanFiles(c("WASH","sQCcond","blank", "cond", "iterative"))
sQCs <- LCMS_data %>% getRP() %>% getNEG() %>% getQC("sQC")
Batches<- stri_extract_first_regex(sQCs$batch, "[0-9]+")
sQCs_batch_specific <- sQCs[Batches==1,] #run IPO per batch
sQCs_batch_names <- sQCs_batch_specific %>% extractNames() %>% sample(5)`

#optimize_peak_picking-----------------------------------------------------------------
peakpickingParameters <- getDefaultXcmsSetStartingParams('centWave')
peakpickingParameters$min_peakwidth <- c(5, 9)
peakpickingParameters$max_peakwidth <- c(20, 40)
peakpickingParameters$noise <- 500 
peakpickingParameters$ppm <- c(10, 30) 
peakpickingParameters$prefilter <- 3
peakpickingParameters$value_of_prefilter <-3000   
peakpickingParameters$snthresh <- 10
time.xcmsSet <- system.time({ # measuring time, time to do optimization
  resultPeakpicking <- 
    optimizeXcmsSet(files = sQCs_batch_names, 
                    params = peakpickingParameters,
                    isotopeIdentification= "IPO",
                    BPPARAM = MulticoreParam(12),
                    subdir = "location", #creating folder called "location" to save the IPO plots
                    plot = TRUE)
})

#print the best settings to a .txt file--------------------------------------
sink(paste(ppath, "IPO/RPpos/QCs_batch1_bestsettings.txt", sep = ""), append=TRUE)
print(resultPeakpicking$best_settings$result)
sink()
    
#save the results to an rds object if necessary------------------------------
optimizedXcmsSetObject <- resultPeakpicking$best_settings$xset
saveRDS(optimizedXcmsSetObject, file = paste0(ppath,"IPO/RPpos/date_IPO_QCs.rds"))
dev.off()

resultPeakpicking$best_settings$result

# time for optimizing peak picking parameters
time.xcmsSet 

#version 1.1 checked by Stef - Dec 2, 2021