script ipo e xcms.R

#IPO e XCMS script - Elisa R. (elisarm.antunes@gmail.com)

##-------------------- Package References ----------------------##

# Websites:

# https://bioconductor.org/packages/release/bioc/html/xcms.html
# https://www.bioconductor.org/packages/release/bioc/html/IPO.html


# Citation

citation("xcms")
citation("IPO")

##-------------------- installation ----------------------##

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

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

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

BiocManager::install("xcms") # install xcms

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

BiocManager::install("CAMERA") # install CAMERA.

##-------------------- files ----------------------##

datafiles <- list.files("D:/...", recursive = TRUE, full.names = TRUE)
datafiles # lists files names. Check if the correct files were found

# OBS: files need to be .mzxML format -  Check MSConvert from Proteowizard.

##-------------------- loading packages ----------------------##

library(xcms)
library(Rmpi)
library(CAMERA)
library(IPO)
library(beepr) 

##-------------------- IPO ----------------------##

# peak picking parameters optimization
peakpickingParameters <- getDefaultXcmsSetStartingParams('matchedFilter') #only matchedFilter (lowres)

peakpickingParameters$step <- c(0.1, 0.2)
peakpickingParameters$fwhm <- c(1, 2)
peakpickingParameters$steps <- 2

time.xcmsSet <- system.time({
  resultPeakpicking <-
    optimizeXcmsSet(files = datafiles[x:x], # check datafiles object. Input only the files to be tested
                    params = peakpickingParameters,
                    nSlaves = 1,# only change if your computer can handle it.
                    subdir = 'C:/...', #can be null if no plot is to be saved
                    plot = TRUE)
})

beep()
resultPeakpicking$best_settings$result
optimizedXcmsSetObject <- resultPeakpicking$best_settings$xset

# retention time correction optimization
# OBS: both chunks take a long time, depending on the amount of files.
# Save the R image or object (save.image() + load () | saveRDS() + readRDS()) to keep the objects stored

retcorGroupParameters <- getDefaultRetGroupStartingParams()
retcorGroupParameters$profStep <- 1
retcorGroupParameters$gapExtend <- 2.7
retcorGroupParameters$mzwid <-
  time.RetGroup <- system.time({
    resultRetcorGroup <-
      optimizeRetGroup(xset = optimizedXcmsSetObject,
                       params = retcorGroupParameters,
                       nSlaves = 1, #only increase if the PC can handle it
                       subdir = 'C:/...', #can be set no NULL
                       plot = TRUE)
  })

beep()

writeRScript(resultPeakpicking$best_settings$parameters,
             resultRetcorGroup$best_settings)

time.xcmsSet
time.RetGroup

sessionInfo()

##-------------------- XCMS ----------------------##

# Separate the QC and other files on subdirectories - based on groups.

setwd ("C:/...")

# the script generated by IPO is similar to the one below. Substitute this chunk
# with the one generated by IPO and rename the objects so they don't overwrite each other.


xset <- xcmsSet(
  method   = "matchedFilter",
  fwhm     = x,
  snthresh = x,
  step     = x,
  steps    = x,
  sigma    = x,
  max      = x,
  mzdiff   = x,
  index    = x,)

xset2 <- retcor(
  xset,
  method         = "obiwarp",
  plottype       = "none",
  distFunc       = "cor_opt",
  profStep       = x, #default =1 possible bug, if so, change it to .99
  center         = x,
  response       = x,
  gapInit        = x,
  gapExtend      = x,
  factorDiag     = x,
  factorGap      = x,
  localAlignment = x,)

xset3 <- group(
  xset2,
  method  = "density",
  bw      = 10,
  mzwid   = 0.1,
  minfrac = 1,
  minsamp = 1,
  max     = 100)

xset4 <- fillPeaks(xset3)

# The IPO script ends here

# Substitute the object names inside the ( ) accordingly.

an <- xsAnnotate(xset4)

#Creation of an xsAnnotate object

anF <- groupFWHM(an, perfwhm = 0.6)

#Perfwhm = parameter defines the window width, which is used for matching

anI <- findIsotopes(anF, mzabs=0.01)

#Mzabs = the allowed m/z error

anIC <- groupCorr(anI, cor_eic_th=0.75)
anFA <- findAdducts(anIC, polarity="positive") #change polarity accordingly

write.csv(getPeaklist(anIC), file="Name.csv") # generates a table of features
beep()