---
tags: Source
---

# 1.3. PeakPicking_preBW.R

#load libraries
library(xcms)
library(CMSITools)
library(tidyverse)
library(stringi)
library(StatTools)
register(bpstart(MulticoreParam(12))) #in BIANCA, in Windows SnowParam
source('/tmp/src/getFilesMAX.R')

# Summary:
# This scripts performs Peak Picking using xcms, based on IPO fitted parameters

# Inputs:
# mzml files location
LCMS_filepath <- "" 

#based on parameters optimized using PeakPicking_wIPO.R
peakwidth <- c(8.05,61.11)
ppm <- 20.68
mzdiff <- 0.00174353
prefilter <- c(3,1000)
noise <- 500

# Outputs:
# directory path to the peak picking file
PP_filepath <- "" 
# name of the peak picking file
filename <- "" 

# Before in pipeline:
# - PeakPicking_wIPO.R

# Next in pipeline:
# - IPO_RTcorrection.R #to optimize several parameters for alignment using Obiwarp method. 
#Otherwise, use PeakGroup method in
# - BWalignment.R #to align peaks based on multiple bandwidth values

filepath <- file.path(filepath, filename)
LCMS_data <- getFilesMAX(files_location)


#----------------------------------define groups-------------------------------------------------------------------------
LCMS_data$group <- rep(NA,nrow(LCMS_data))
LCMS_data$group[LCMS_data$sample %>% grep("sQC",.)] <- "sQC"
LCMS_data$group[LCMS_data$sample %>% grep("ltQC",.)] <- "ltQC"
LCMS_data$group[-c(LCMS_data$sample %>% grep("sQC",.), LCMS_data$sample %>% grep("ltQC",.))] <- "sample"
LCMS_data_names <- LCMS_data %>% extractNames()



#----------------------------------peakpicking---------------------------------------------------------------------------
pd <- data.frame(
  sample_name = sub(basename(LCMS_data_names), pattern = ".mzML",
                    replacement = "", fixed = TRUE),
  sample_group = LCMS_data$group, stringsAsFactors = FALSE)

# TEAM: not sure if the rt parameter here should be changed
raw_data <- readMSData(files = LCMS_data_names, pdata = new("NAnnotatedDataFrame", pd),
                              mode = "onDisk") %>% filterRt(.,rt = c(0, 630)) 
                              #for HILIC and c(0,930) for RP, or depending on the length of the column
#if there is any missing spectra
raw_data <- filterEmptySpectra(raw_data)
#if necessary to filter based on the Mz, too
raw_data <- filterMz(raw_data, c(0, 1000))

# runs the peak picking
cwp <- CentWaveParam(peakwidth=peakwidth,ppm=ppm, mzdiff=mzdiff, prefilter=prefilter, noise=noise)
xdata <- findChromPeaks(raw_data, param = cwp) 
saveRDS(xdata, file=filepath)

version 1.1 - Stef - Dec 2, 2021