R/idFunctionsNeg.R

Defines functions idSMneg idBAneg idCLneg idAcylCerneg idCerPneg idCerneg idSphPneg idSphneg idPSneg idPIneg idPGneg idPEpneg idPEoneg idPEneg idPCpneg idPConeg idPCneg idLPSneg idLPIneg idLPGneg idLPEneg idLPCneg idFAHFAneg idFAneg idNEG

Documented in idAcylCerneg idBAneg idCerneg idCerPneg idCLneg idFAHFAneg idFAneg idLPCneg idLPEneg idLPGneg idLPIneg idLPSneg idNEG idPCneg idPConeg idPCpneg idPEneg idPEoneg idPEpneg idPGneg idPIneg idPSneg idSMneg idSphneg idSphPneg

# idNEG
#' Lipids annotation for ESI-
#'
#' Lipids annotation based on fragmentation patterns for LC-MS/MS DIA or DDA data
#' acquired in negative mode. This function compiles all functions writen for
#' ESI- annotations.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 5 seconds.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param lipidClasses classes of interest to run the identification functions.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#' 
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification); and the annotatedPeaklist element shows the original 
#' MS1 peaklist with the annotations on it.
#'
#' @examples
#' \dontrun{
#' msobject <- idNEG(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maribel_alcoriza@iislafe.es>
idNEG <- function(msobject,
                  ppm_precursor = 5,
                  ppm_products = 10,
                  rttol = 5,
                  coelCutoff = 0.8,
                  lipidClasses = c("FA", "FAHFA", "LPC", "LPE", "LPG", "LPI",
                                   "LPS", "PC", "PCo", "PCp", "PE", "PEo", "PEp", 
                                   "PG", "PI", "PS", "Sph", "SphP", "Cer", "CerP", 
                                   "AcylCer", "SM", "CL", "BA"),
                  dbs,
                  verbose = TRUE){

  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if ("results" %in% names(msobject$annotation)){
    if(verbose){cat("\n Removing previous results...")}
    msobject$annotation$results <- NULL
    msobject$annotation$detailsAnnotation <- NULL
    msobject$annotatedPeaklist <- NULL
    if(verbose){cat("OK")}
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(lipidClasses %in% c("FA", "FAHFA", "LPC", "LPE", "LPG", "LPI",
                               "LPS", "PC", "PCo", "PCp", "PE", "PEo", "PEp", 
                               "PG", "PI", "PS", "Sph", "SphP", "Cer", "CerP", 
                               "AcylCer", "SM", "CL", "BA"))){
    stop("Lipid classes allowed for negative annotation are: FA, FAHFA, LPC, LPE,
         LPG, LPI, LPS, PC, PCo, PCp, PE, PEo, PEp, PG, PI, PS, Sph, SphP, Cer, 
         CerP, AcylCer, SM, CL and BA")
  }

  if(verbose){cat("\n Starting annotation...")}
  if ("FA" %in% lipidClasses){
    if(verbose){cat("\n  Searching for FA...")}
    msobject <-  idFAneg(msobject = msobject, ppm_precursor = ppm_precursor,
                         ppm_products = ppm_products, rttol = rttol,
                         coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("FAHFA" %in% lipidClasses){
    if(verbose){cat("\n  Searching for FAHFA...")}
    msobject <-  idFAHFAneg(msobject = msobject, ppm_precursor = ppm_precursor,
                         ppm_products = ppm_products, rttol = rttol,
                         coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("LPC" %in% lipidClasses){
    if(verbose){cat("\n  Searching for LPC...")}
    msobject <-  idLPCneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("LPE" %in% lipidClasses){
    if(verbose){cat("\n  Searching for LPE...")}
    msobject <-  idLPEneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("LPG" %in% lipidClasses){
    if(verbose){cat("\n  Searching for LPG...")}
    msobject <-  idLPGneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("LPI" %in% lipidClasses){
    if(verbose){cat("\n  Searching for LPI...")}
    msobject <-  idLPIneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("LPS" %in% lipidClasses){
    if(verbose){cat("\n  Searching for LPS...")}
    msobject <-  idLPSneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("PC" %in% lipidClasses){
    if(verbose){cat("\n  Searching for PC...")}
    msobject <-  idPCneg(msobject = msobject, ppm_precursor = ppm_precursor,
                         ppm_products = ppm_products, rttol = rttol,
                         coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("PCo" %in% lipidClasses){
    if(verbose){cat("\n  Searching for PCo...")}
    msobject <-  idPConeg(msobject = msobject, ppm_precursor = ppm_precursor,
                         ppm_products = ppm_products, rttol = rttol,
                         coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("PCp" %in% lipidClasses){
    if(verbose){cat("\n  Searching for PCp...")}
    msobject <-  idPCpneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("PE" %in% lipidClasses){
    if(verbose){cat("\n  Searching for PE...")}
    msobject <-  idPEneg(msobject = msobject, ppm_precursor = ppm_precursor,
                         ppm_products = ppm_products, rttol = rttol,
                         coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("PEo" %in% lipidClasses){
    if(verbose){cat("\n  Searching for PEo...")}
    msobject <-  idPEoneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("PEp" %in% lipidClasses){
    if(verbose){cat("\n  Searching for PEp...")}
    msobject <-  idPEpneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("PG" %in% lipidClasses){
    if(verbose){cat("\n  Searching for PG...")}
    msobject <-  idPGneg(msobject = msobject, ppm_precursor = ppm_precursor,
                         ppm_products = ppm_products, rttol = rttol,
                         coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("PI" %in% lipidClasses){
    if(verbose){cat("\n  Searching for PI...")}
    msobject <-  idPIneg(msobject = msobject, ppm_precursor = ppm_precursor,
                         ppm_products = ppm_products, rttol = rttol,
                         coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("PS" %in% lipidClasses){
    if(verbose){cat("\n  Searching for PS...")}
    msobject <-  idPSneg(msobject = msobject, ppm_precursor = ppm_precursor,
                         ppm_products = ppm_products, rttol = rttol,
                         coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("Sph" %in% lipidClasses){
    if(verbose){cat("\n  Searching for Sph...")}
    msobject <-  idSphneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("SphP" %in% lipidClasses){
    if(verbose){cat("\n  Searching for SphP...")}
    msobject <-  idSphPneg(msobject = msobject, ppm_precursor = ppm_precursor,
                           ppm_products = ppm_products, rttol = rttol,
                           coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("Cer" %in% lipidClasses){
    if(verbose){cat("\n  Searching for Cer...")}
    msobject <-  idCerneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("CerP" %in% lipidClasses){
    if(verbose){cat("\n  Searching for CerP...")}
    msobject <-  idCerPneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("AcylCer" %in% lipidClasses){
    if(verbose){cat("\n  Searching for AcylCer...")}
    msobject <-  idAcylCerneg(msobject = msobject, ppm_precursor = ppm_precursor,
                           ppm_products = ppm_products, rttol = rttol,
                           coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("SM" %in% lipidClasses){
    if(verbose){cat("\n  Searching for SM...")}
    msobject <-  idSMneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("CL" %in% lipidClasses){
    if(verbose){cat("\n  Searching for CL...")}
    msobject <-  idCLneg(msobject = msobject, ppm_precursor = ppm_precursor,
                         ppm_products = ppm_products, rttol = rttol,
                         coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }
  if ("BA" %in% lipidClasses){
    if(verbose){cat("\n  Searching for Bile acids...")}
    msobject <-  idBAneg(msobject = msobject, ppm_precursor = ppm_precursor,
                          ppm_products = ppm_products, rttol = rttol,
                          coelCutoff = coelCutoff, dbs = dbs, verbose = verbose)
    if(verbose){cat("OK")}
  }

  if(verbose){cat("\n Preparing output...")}
  msobject <- crossTables(msobject,
                          ppm = ppm_precursor, 
                          rttol = rttol,
                          dbs = dbs)
  if(verbose){cat("OK\n")}
  return(msobject)
}

# idFAneg
#' Fatty Acids (FA) annotation for ESI-
#'
#' FA identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for FA in ESI-. Adducts allowed can
#' be modified in addutcsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idFAneg} function involves 2 steps. 1) FullMS-based
#' identification of candidate FA as M-H or 2M-H. 2) Search of FA class
#' fragments: neutral loss of H2O coeluting with the precursor ion or the
#' molecular ion.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (in this
#' case, just MS-only or Subclass level (if any class fragment is defined) are
#' possible) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idFAneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maribel_alcoriza@iislafe.es>
idFAneg <- function(msobject,
                    ppm_precursor = 5,
                    ppm_products = 10,
                    rttol = 3,
                    rt,
                    adducts = c("M-H", "2M-H"),
                    clfrags = c("fa_M-H", "fa_M-H-H2O"),
                    clrequired = c(FALSE, FALSE),
                    ftype = c("BB", "BB"),
                    coelCutoff = 0.8,
                    dbs,
                    verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "FA",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("FA" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious FA annotations removed")}
      msobject$annotation$detailsAnnotation$FA <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$fadb, ppm = ppm_precursor,
                               rt = rt, adducts = adducts, rttol = rttol,
                               dbs = dbs, rawData = rawData,
                               coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb = list(),
                           intrules  = c(), intConf = list(), nchains = 0,
                           class="FA",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$FA <- list()
    msobject$annotation$detailsAnnotation$FA$candidates <- candidates
    msobject$annotation$detailsAnnotation$FA$classfragments <- classConf$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$FA$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$FA <- list()
  }
  return(msobject)
}

# idFAHFAneg
#' FAHFA annotation for ESI-
#'
#' FAHFA identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for FAHFA in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idFAHFAneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate FAHFA as M-H. 2) Search of FAHFA class fragments:
#' there is't any class fragment by default. 3) Search of specific fragments
#' that inform about chain composition in sn1 (HFA as M-H resulting from the
#' loss of the FA chain) and sn2 (FA chain as M-H). 4) Look for possible
#' chains structure based on the combination of chain fragments. 5) Check
#' intensity rules to confirm chains position. In this case, HFA intensity has
#' to be higher than FA.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idFAHFAneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idFAHFAneg <- function(msobject,
                       ppm_precursor = 5,
                       ppm_products = 10,
                       rttol = 3,
                       rt,
                       adducts = c("M-H"),
                       clfrags = c(),
                       clrequired = c(),
                       ftype = c(),
                       chainfrags_sn1 = c("hfa_M-H"),
                       chainfrags_sn2 = c("fa_M-H"),
                       intrules = c("hfa_sn1/fa_sn2"),
                       rates = c("3/1"),
                       intrequired = c(T),
                       coelCutoff = 0.8,
                       dbs,
                       verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous FAHFA annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "FAHFA",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("FAHFA" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious FAHFA annotations removed")}
      msobject$annotation$detailsAnnotation$FAHFA <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$fahfadb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="FAHFA",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$FAHFA <- list()
    msobject$annotation$detailsAnnotation$FAHFA$candidates <- candidates
    msobject$annotation$detailsAnnotation$FAHFA$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$FAHFA$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$FAHFA$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$FAHFA <- list()
  }
  return(msobject)
}

# idLPCneg
#' Lysophosphocholines (LPC) annotation for ESI-
#'
#' LPC identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for LPC in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments. See \link{chainFrags} for details.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idLPCneg} function involves 3 steps. 1) FullMS-based
#' identification of candidate LPC as M+CH3COO, M-CH3 and M+CH3COO-CH3. To avoid
#' incorrect annotations of PE as PC, candidates which are present just as M-CH3
#' will be ignored. 2) Search of LPC class fragments: 168.0426, 224.0688, lysoPA
#' as M-H or lysoPC as M-CH3 coeluting with the precursor ion. 3) Search of
#' specific fragments that confirm chain composition (FA as M-H).
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (in this
#' case, as LPC only have one chain, only Subclass and FA level are possible)
#' and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idLPCneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idLPCneg <- function(msobject,
                     ppm_precursor = 5,
                     ppm_products = 10,
                     rttol = 3,
                     rt,
                     adducts = c("M+CH3COO", "M-CH3", "M+CH3COO-CH3"),
                     clfrags = c(168.0426, 224.0688, "lysopa_M-H", "lysopc_M-CH3"),
                     clrequired = c(F, F, F, F),
                     ftype = c("F", "F", "BB", "BB"),
                     chainfrags_sn1 = c("fa_M-H"),
                     coelCutoff = 0.8,
                     dbs,
                     verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous LPC annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "LPC",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("LPC" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious LPC annotations removed")}
      msobject$annotation$detailsAnnotation$LPC <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$lysopcdb, ppm = ppm_precursor,
                               rt = rt, adducts = adducts, rttol = rttol,
                               dbs = dbs, rawData = rawData,
                               coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=1, sn1)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules = c(), rates = c(),
                                   intrequired = c(), nchains=1,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb,
                           intrules  = c(), intConf, nchains = 1, class="LPC",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$LPC <- list()
    msobject$annotation$detailsAnnotation$LPC$candidates <- candidates
    msobject$annotation$detailsAnnotation$LPC$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$LPC$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$LPC$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$LPC <- list()
  }
  return(msobject)
}

# idLPEneg
#' Lysophosphoethanolamines (LPE) annotation for ESI-
#'
#' LPE identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for LPE in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments. See \link{chainFrags} for details.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idLPEneg} function involves 3 steps. 1) FullMS-based
#' identification of candidate LPE as M-H. 2) Search of
#' LPE class fragments: 140.0115, 196.038 and 214.048 coeluting with the
#' precursor ion. If a loss of CH3 group is found coeluting with any candidate,
#' this will be excluded as it is a characteristic fragment of LPC.3) Search of
#' specific fragments that confirm chain composition (FA as M-H).
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (in this
#' case, as LPE only have one chain, only Subclass and FA level are possible)
#' and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idLPEneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idLPEneg <- function(msobject, ppm_precursor = 5,
                     ppm_products = 10,
                     rttol = 3,
                     rt,
                     adducts = c("M-H"),
                     clfrags = c(140.0115, 196.038, 214.048, "lysope_M-CH3"),
                     clrequired = c(F, F, F, "excluding"),
                     ftype = c("F", "F", "F", "BB"),
                     chainfrags_sn1 = c("fa_M-H"),
                     coelCutoff = 0.8,
                     dbs,
                     verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "LPE",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("LPE" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious LPE annotations removed")}
      msobject$annotation$detailsAnnotation$LPE <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$lysopedb, ppm = ppm_precursor,
                               rt = rt, adducts = adducts, rttol = rttol,
                               dbs = dbs, rawData = rawData,
                               coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=1, sn1)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules = c(), rates = c(),
                                   intrequired = c(), nchains=1,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb,
                           intrules  = c(), intConf, nchains = 1, class="LPE",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$LPE <- list()
    msobject$annotation$detailsAnnotation$LPE$candidates <- candidates
    msobject$annotation$detailsAnnotation$LPE$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$LPE$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$LPE$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$LPE <- list()
  }
  return(msobject)
}

# idLPGneg
#' Lysophosphoglycerols (LPG) annotation for ESI-
#'
#' LPG identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for LPG in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments. See \link{chainFrags} for details.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idLPGneg} function involves 3 steps. 1) FullMS-based
#' identification of candidate LPG as M-H. 2) Search of LPG class fragments:
#' 152.9958, 227.0326, 209.022 and neutral loss of 74.0359 coeluting with the
#' precursor ion. 3) Search of specific fragments that confirm chain composition
#' (FA as M-H).
#'
#' Results data frame shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (in this
#' case, as LPG only have one chain, only Subclass and FA level are possible)
#' and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idLPGneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idLPGneg <- function(msobject,
                     ppm_precursor = 5,
                     ppm_products = 10,
                     rttol = 3,
                     rt,
                     adducts = c("M-H"),
                     clfrags = c(152.9958, 227.0326, 209.022, 74.0359),
                     clrequired = c(F, F, F, F),
                     ftype = c("F", "F", "F", "NL"),
                     chainfrags_sn1 = c("fa_M-H"),
                     coelCutoff = 0.8,
                     dbs,
                     verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "LPG",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("LPG" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious LPG annotations removed")}
      msobject$annotation$detailsAnnotation$LPG <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$lysopgdb, ppm = ppm_precursor,
                               rt = rt, adducts = adducts, rttol = rttol,
                               dbs = dbs, rawData = rawData,
                               coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=1, sn1)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules = c(), rates = c(),
                                   intrequired = c(), nchains=1,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb,
                           intrules  = c(), intConf, nchains = 1, class="LPG",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$LPG <- list()
    msobject$annotation$detailsAnnotation$LPG$candidates <- candidates
    msobject$annotation$detailsAnnotation$LPG$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$LPG$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$LPG$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$LPG <- list()
  }
  return(msobject)
}

# idLPIneg
#' Lysophosphoinositols (LPI) annotation for ESI-
#'
#' LPI identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for LPI in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments. See \link{chainFrags} for details.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idLPIneg} function involves 3 steps. 1) FullMS-based
#' identification of candidate LPI as M-H. 2) Search of
#' LPI class fragments: 241.0115, 223.0008, 259.0219 and 297.0375 coeluting
#' with the precursor ion. 3) Search of specific fragments that confirm chain
#' composition (FA as M-H).
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (in this
#' case, as LPI only have one chain, only Subclass and FA level are possible)
#' and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idLPIneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idLPIneg <- function(msobject,
                     ppm_precursor = 5,
                     ppm_products = 10,
                     rttol = 3,
                     rt,
                     adducts = c("M-H"),
                     clfrags = c(241.0115, 223.0008, 259.0219, 297.0375),
                     clrequired = c(F, F, F, F),
                     ftype = c("F", "F", "F", "F"),
                     chainfrags_sn1 = c("fa_M-H"),
                     coelCutoff = 0.8,
                     dbs, 
                     verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "LPI",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("LPI" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious LPI annotations removed")}
      msobject$annotation$detailsAnnotation$LPI <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$lysopidb, ppm = ppm_precursor,
                               rt = rt, adducts = adducts, rttol = rttol,
                               dbs = dbs, rawData = rawData,
                               coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=1, sn1)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules = c(), rates = c(),
                                   intrequired = c(), nchains=1,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb,
                           intrules  = c(), intConf, nchains = 1, class="LPI",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$LPI <- list()
    msobject$annotation$detailsAnnotation$LPI$candidates <- candidates
    msobject$annotation$detailsAnnotation$LPI$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$LPI$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$LPI$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$LPI <- list()
  }
  return(msobject)
}

# idLPSneg
#' Lysophosphoserines (LPS) annotation for ESI-
#'
#' LPS identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for LPS in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments. See \link{chainFrags} for details.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idLPSneg} function involves 3 steps. 1) FullMS-based
#' identification of candidate LPS as M-H and M+Na-2H. 2) Search of
#' LPS class fragments: neutral loss of 87.032 coeluting with the precursor ion.
#' 3) Search of specific fragments that confirm chain composition (FA as M-H).
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (in this
#' case, as LPS only have one chain, only Subclass and FA level are possible)
#' and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idLPSneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idLPSneg <- function(msobject,
                     ppm_precursor = 5,
                     ppm_products = 10,
                     rttol = 3,
                     rt,
                     adducts = c("M-H", "M+Na-2H"),
                     clfrags = c(87.032),
                     clrequired = c(F),
                     ftype = c("NL"),
                     chainfrags_sn1 = c("fa_M-H"),
                     coelCutoff = 0.8,
                     dbs,
                     verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "LPS",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("LPS" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious LPS annotations removed")}
      msobject$annotation$detailsAnnotation$LPS <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$lysopsdb, ppm = ppm_precursor,
                               rt = rt, adducts = adducts, rttol = rttol,
                               dbs = dbs, rawData = rawData,
                               coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=1, sn1)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules = c(), rates = c(),
                                   intrequired = c(), nchains=1,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb,
                           intrules  = c(), intConf, nchains = 1, class="LPS",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$LPS <- list()
    msobject$annotation$detailsAnnotation$LPS$candidates <- candidates
    msobject$annotation$detailsAnnotation$LPS$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$LPS$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$LPS$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$LPS <- list()
  }
  return(msobject)
}

# idPCneg
#' Phosphocholines (PC) annotation for ESI-
#'
#' PC identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for PC in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idPCneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate PC as M+CH3COO, M-CH3 or M+CH3COO-CH3. To avoid
#' incorrect annotations of PE as PC, candidates which are present just as M-CH3
#' will be ignored. 2) Search of PC class fragments: 168.0426, 224.0688 or loss
#' of CH3 coeluting with the precursor ion. 3) Search of specific fragments that
#' inform about chain composition in sn1 (lysoPC as M-CH3 resulting from the
#' loss of the FA chain at sn2) and sn2 (lysoPC as M-CH3 resulting from the loss
#' of sn1 or FA as M-H). 4) Look for possible chains structure based on the
#' combination of chain fragments. 5) Check intensity rules to confirm chains
#' position. In this case, lysoPC from sn1 is at least 3 times more intense than
#' lysoPC from sn2.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idPCneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idPCneg <- function(msobject,
                    ppm_precursor = 5,
                    ppm_products = 10,
                    rttol = 3,
                    rt,
                    adducts = c("M+CH3COO", "M-CH3", "M+CH3COO-CH3"),
                    clfrags = c(168.0426, 224.0688, "pc_M-CH3"),
                    clrequired = c(F, F, F),
                    ftype = c("F", "F", "BB"),
                    chainfrags_sn1 = c("lysopc_M-CH3"),
                    chainfrags_sn2 = c("fa_M-H", "lysopc_M-CH3"),
                    intrules = c("lysopc_sn1/lysopc_sn2"),
                    rates = c("3/1"),
                    intrequired = c(T),
                    coelCutoff = 0.8,
                    dbs,
                    verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "PC",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("PC" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious PC annotations removed")}
      msobject$annotation$detailsAnnotation$PC <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$pcdb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)
  # remove PC which ony appear as M-CH3
  if(length(adducts) > 1 & "M-CH3" %in% adducts){
    candidates <- candidates[candidates$adducts != "M-CH3",]
  }

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="PC",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PC <- list()
    msobject$annotation$detailsAnnotation$PC$candidates <- candidates
    msobject$annotation$detailsAnnotation$PC$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$PC$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$PC$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PC <- list()
  }
  return(msobject)
}

# idPConeg
#' Plasmanyl Phosphocholines (PCo) annotation for ESI-
#'
#' PCo identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for PCo in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idPConeg} function involves 5 steps. 1) FullMS-based
#' identification of candidate PCo as M+CH3COO, M-CH3 or M+CH3COO-CH3. To avoid
#' incorrect annotations of PEo as PCo, candidates which are present just as M-CH3
#' will be ignored. 2) Search of PCo class fragments: 168.0426, 224.0688 or loss
#' of CH3 coeluting with the precursor ion. 3) Search of specific fragments that
#' inform about chain composition in sn1 (LPCo as M-CH3 and M-CH3-H2O resulting 
#' from the loss of the FA chain at sn2) and sn2 (FA as M-H and M-CO2-H). 
#' 4) Look for possible chains structure based on the  combination of chain 
#' fragments. 5) Check intensity rules to confirm chains position. In this case, 
#' FA fragments from sn2 are at least 3 times more intense than LPCo from sn1.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idPCneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idPConeg <- function(msobject,
                    ppm_precursor = 5,
                    ppm_products = 10,
                    rttol = 3,
                    rt,
                    adducts = c("M+CH3COO", "M-CH3", "M+CH3COO-CH3"),
                    clfrags = c(168.0426, 224.0688, "pco_M-CH3"),
                    clrequired = c(F, F, F),
                    ftype = c("F", "F", "BB"),
                    chainfrags_sn1 = c("lysopco_M-CH3", "lysopco_M-CH3-H2O"),
                    chainfrags_sn2 = c("fa_M-H", "fa_M-CO2-H"),
                    intrules = c("lysopco_sn1/fa_sn2"),
                    rates = c(1/3),
                    intrequired = c(T),
                    coelCutoff = 0.8,
                    dbs,
                    verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "PCo",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("PCo" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious PCo annotations removed")}
      msobject$annotation$detailsAnnotation$PCo <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps
  
  # candidates search
  candidates <- findCandidates(MS1, dbs$pcodb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)
  # remove PCo which ony appear as M-CH3 to avoid wrong annotations of PEo
  if(length(adducts) > 1 & "M-CH3" %in% adducts){
    candidates <- candidates[candidates$adducts != "M-CH3",]
  }
  
  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }
    
    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)
    
    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)
    
    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)
    
    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)
    
    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="PCo",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)
    
    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PCo <- list()
    msobject$annotation$detailsAnnotation$PCo$candidates <- candidates
    msobject$annotation$detailsAnnotation$PCo$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$PCo$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$PCo$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PCo <- list()
  }
  return(msobject)
}

# idPCpneg
#' Plasmenyl Phosphocholines (PCp) annotation for ESI-
#'
#' PCp identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for PCp in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idPCpneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate PCp as M+CH3COO, M-CH3 or M+CH3COO-CH3. To avoid
#' incorrect annotations of PEp as PCp, candidates which are present just as M-CH3
#' will be ignored. 2) Search of PCp class fragments: 168.0426, 224.0688 or loss
#' of CH3 coeluting with the precursor ion. 3) Search of specific fragments that
#' inform about chain composition in sn1 (LPCp as M-CH3 and M-CH3-H2O resulting 
#' from the loss of the FA chain at sn2) and sn2 (FA as M-H and M-CO2-H). 
#' 4) Look for possible chains structure based on the  combination of chain 
#' fragments. 5) Check intensity rules to confirm chains position. In this case, 
#' FA fragments from sn2 are at least 3 times more intense than LPCp from sn1.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idPCpneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idPCpneg <- function(msobject,
                     ppm_precursor = 5,
                     ppm_products = 10,
                     rttol = 3,
                     rt,
                     adducts = c("M+CH3COO", "M-CH3", "M+CH3COO-CH3"),
                     clfrags = c(168.0426, 224.0688, "pcp_M-CH3"),
                     clrequired = c(F, F, F),
                     ftype = c("F", "F", "BB"),
                     chainfrags_sn1 = c("lysopcp_M-CH3", "lysopcp_M-CH3-H2O"),
                     chainfrags_sn2 = c("fa_M-H", "fa_M-CO2-H"),
                     intrules = c("lysopcp_sn1/fa_sn2"),
                     rates = c(1/3),
                     intrequired = c(T),
                     coelCutoff = 0.8,
                     dbs,
                     verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "PCp",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("PCp" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious PC annotations removed")}
      msobject$annotation$detailsAnnotation$PCp <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps
  
  # candidates search
  candidates <- findCandidates(MS1, dbs$pcpdb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)
  # remove PCp which ony appear as M-CH3 to avoid wrong annotations of PEp
  if(length(adducts) > 1 & "M-CH3" %in% adducts){
    candidates <- candidates[candidates$adducts != "M-CH3",]
  }
  
  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }
    
    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)
    
    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)
    
    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)
    
    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)
    
    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="PCp",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)
    
    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PCp <- list()
    msobject$annotation$detailsAnnotation$PCp$candidates <- candidates
    msobject$annotation$detailsAnnotation$PCp$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$PCp$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$PCp$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PCp <- list()
  }
  return(msobject)
}

# idPEneg
#' Phosphoethanolamines (PE) annotation for ESI-
#'
#' PE identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for PE in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idPEneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate PE as M-H. 2) Search of PE class fragments:
#' 140.0115, 196.038, 214.048 ion coeluting with the precursor ion. If a loss of
#' CH3 group is found coeluting with any candidate, this will be excluded as it
#' is a characteristic fragment of PC. 3) Search of specific fragments that
#' inform about chain composition in sn1 (lysoPE as M-H resulting from the loss
#' of the FA chain at sn2) and sn2 (lysoPE as M-H resulting from the loss
#' of the FA chain at sn1 or FA chain as M-H). 4) Look for possible
#' chains structure based on the combination of chain fragments. 5) Check
#' intensity rules to confirm chains position. In this case, lysoPE from sn1 is
#' at least 3 times more intense than lysoPE from sn2.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idPEneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idPEneg <- function(msobject,
                    ppm_precursor = 5,
                    ppm_products = 10,
                    rttol = 5,
                    rt,
                    adducts = c("M-H"),
                    clfrags = c(140.0118, 196.038, 214.048, "pe_M-CH3"),
                    clrequired = c(F, F, F, "excluding"),
                    ftype = c("F", "F", "F", "BB"),
                    chainfrags_sn1 = c("lysope_M-H"),
                    chainfrags_sn2 = c("lysope_M-H", "fa_M-H"),
                    intrules = c("lysope_sn1/lysope_sn2"),
                    rates = c("3/1"),
                    intrequired = c(T),
                    coelCutoff = 0.8,
                    dbs,
                    verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "PE",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("PE" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious PE annotations removed")}
      msobject$annotation$detailsAnnotation$PE <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$pedb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="PE",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PE <- list()
    msobject$annotation$detailsAnnotation$PE$candidates <- candidates
    msobject$annotation$detailsAnnotation$PE$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$PE$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$PE$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PE <- list()
  }
  return(msobject)
}

# idPEoneg
#' Plasmanyl Phosphoethanolamines (PEo) annotation for ESI-
#'
#' PEo identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for PEo in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idPEoneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate PEo as M-H and M+NaCH3COO. 2) Search 
#' of PEo class fragments: 140.0115, 196.038, 214.048 ion coeluting with the 
#' precursor ion. If a loss of CH3 group is found coeluting with any candidate, 
#' this will be excluded as it is a characteristic fragment of PCo. 3) Search of 
#' specific fragments that inform about chain composition in sn1 (lysoPEo as M-H 
#' and M-H-H2O resulting from the loss of the FA chain at sn2) and sn2 (FA chain 
#' as M-H). 4) Look for possible chains structure based on the combination of 
#' chain fragments. 5) Check intensity rules to confirm chains position. In this 
#' case, FA fragments from sn2 are at least 3 times more intense than LPEo from 
#' sn1.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idPEoneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idPEoneg <- function(msobject,
                    ppm_precursor = 5,
                    ppm_products = 10,
                    rttol = 5,
                    rt,
                    adducts = c("M-H", "M+NaCH3COO"),
                    clfrags = c(140.0118, 196.038, 214.048, "peo_M-CH3"),
                    clrequired = c(F, F, F, "excluding"),
                    ftype = c("F", "F", "F", "BB"),
                    chainfrags_sn1 = c("lysopeo_M-H", "lysopeo_M-H-H2O"),
                    chainfrags_sn2 = c("fa_M-H"),
                    intrules = c("lysopeo_sn1/fa_sn2"),
                    rates = c(1/3),
                    intrequired = c(T),
                    coelCutoff = 0.8,
                    dbs,
                    verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "PEo",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("PEo" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious PEo annotations removed")}
      msobject$annotation$detailsAnnotation$PEo <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps
  
  # candidates search
  candidates <- findCandidates(MS1, dbs$peodb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)
  
  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }
    
    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)
    
    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)
    
    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)
    
    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)
    
    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="PEo",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)
    
    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PEo <- list()
    msobject$annotation$detailsAnnotation$PEo$candidates <- candidates
    msobject$annotation$detailsAnnotation$PEo$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$PEo$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$PEo$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PEo <- list()
  }
  return(msobject)
}

# idPEpneg
#' Plasmenyl Phosphoethanolamines (PEp) annotation for ESI-
#'
#' PEp identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for PEp in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idPEpneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate PEp as M-H and M+NaCH3COO. 2) Search 
#' of PEp class fragments: 140.0115, 196.038, 214.048 ion coeluting with the 
#' precursor ion. If a loss of CH3 group is found coeluting with any candidate, 
#' this will be excluded as it is a characteristic fragment of PCp. 3) Search of 
#' specific fragments that inform about chain composition in sn1 (lysoPEp as M-H 
#' and M-H-H2O resulting from the loss of the FA chain at sn2) and sn2 (FA chain 
#' as M-H). 4) Look for possible chains structure based on the combination of 
#' chain fragments. 5) Check intensity rules to confirm chains position. In this 
#' case, FA fragments from sn2 are at least 3 times more intense than LPEp from 
#' sn1.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idPEoneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idPEpneg <- function(msobject,
                     ppm_precursor = 5,
                     ppm_products = 10,
                     rttol = 5,
                     rt,
                     adducts = c("M-H", "M+NaCH3COO"),
                     clfrags = c(140.0118, 196.038, 214.048, "pep_M-CH3"),
                     clrequired = c(F, F, F, "excluding"),
                     ftype = c("F", "F", "F", "BB"),
                     chainfrags_sn1 = c("lysopep_M-H", "lysopep_M-H-H2O"),
                     chainfrags_sn2 = c("fa_M-H"),
                     intrules = c("lysopep_sn1/fa_sn2"),
                     rates = c(1/3),
                     intrequired = c(T),
                     coelCutoff = 0.8,
                     dbs,
                     verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "PEp",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("PEp" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious PEp annotations removed")}
      msobject$annotation$detailsAnnotation$PEp <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps
  
  # candidates search
  candidates <- findCandidates(MS1, dbs$pepdb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)
  
  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }
    
    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)
    
    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)
    
    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)
    
    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)
    
    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="PEp",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)
    
    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PEp <- list()
    msobject$annotation$detailsAnnotation$PEp$candidates <- candidates
    msobject$annotation$detailsAnnotation$PEp$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$PEp$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$PEp$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PEp <- list()
  }
  return(msobject)
}

# idPGneg
#' Phosphoglycerols (PG) annotation for ESI-
#'
#' PG identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for PG in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idPGneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate PG as M-H. 2) Search of PG class fragments:
#' 152.9958, 227.0326, 209.022 and neutral loss of 74.0359 coeluting with the
#' precursor ion. 3) Search of specific fragments that inform about chain
#' composition at sn1 (lysoPG as M-H resulting from the loss of the FA chain
#' at sn2) and sn2 (lysoPG as M-H resulting from the loss of the FA chain
#' at sn1 or FA chain as M-H). 4) Look for possible chains structure
#' based on the combination of chain fragments. 5) Check intensity rules to
#' confirm chains position. In this case, lysoPG from sn1 is at least 3 times
#' more intense than lysoPG from sn2.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idPGneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idPGneg <- function(msobject, ppm_precursor = 5,
                    ppm_products = 10,
                    rttol = 3,
                    rt,
                    adducts = c("M-H"),
                    clfrags = c(152.9958, 227.0326, 209.022, 74.0359),
                    clrequired = c(F, F, F, F),
                    ftype = c("F", "F", "F", "NL"),
                    chainfrags_sn1 = c("lysopg_M-H"),
                    chainfrags_sn2 = c("lysopg_M-H", "fa_M-H"),
                    intrules = c("lysopg_sn1/lysopg_sn2"),
                    rates = c("2/1"),
                    intrequired = c(T),
                    coelCutoff = 0.8,
                    dbs,
                    verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "PG",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("PG" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious PG annotations removed")}
      msobject$annotation$detailsAnnotation$PG <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$pgdb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="PG",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PG <- list()
    msobject$annotation$detailsAnnotation$PG$candidates <- candidates
    msobject$annotation$detailsAnnotation$PG$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$PG$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$PG$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PG <- list()
  }
  return(msobject)
}

# idPIneg
#' Phosphoinositols (PI) annotation for ESI-
#'
#' PI identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for PI in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idPIneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate PI as M-H. 2) Search of PI class fragments:
#' 241.0115, 223.0008, 259.0219 and 297.0375 coeluting with the precursor
#' ion. 3) Search of specific fragments that inform about chain composition at
#' sn1 (lysoPI as M-H resulting from the loss of the FA chain at sn2 or lysoPA
#' as M-H if it also losses the head group) and sn2 (lysoPI or lysoPA as M-H
#' resulting from the loss of the FA chain at sn1 or FA chain as M-H). 4) Look
#' for possible chains structure based on the combination of chain fragments.
#' 5) Check intensity rules to confirm chains position. In this case, lysoPI or
#' lysoPA from sn1 is at least 3 times more intense than lysoPI or lysoPA from
#'  sn2.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idPIneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idPIneg <- function(msobject,
                    ppm_precursor = 5,
                    ppm_products = 10,
                    rttol = 3,
                    rt,
                    adducts = c("M-H"),
                    clfrags = c(241.0115, 223.0008, 259.0219, 297.0375),
                    clrequired = c(F, F, F, F),
                    ftype = c("F", "F", "F", "F"),
                    chainfrags_sn1 = c("lysopi_M-H", "lysopa_M-H"),
                    chainfrags_sn2 = c("lysopi_M-H", "lysopa_M-H", "fa_M-H"),
                    intrules = c("lysopi_sn1/lysopi_sn2", "lysopa_sn1/lysopa_sn2"),
                    rates = c("3/1", "3/1"),
                    intrequired = c(F, F),
                    coelCutoff = 0.8,
                    dbs,
                    verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "PI",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("PI" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious PI annotations removed")}
      msobject$annotation$detailsAnnotation$Cer <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$pidb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="PI",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PI <- list()
    msobject$annotation$detailsAnnotation$PI$candidates <- candidates
    msobject$annotation$detailsAnnotation$PI$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$PI$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$PI$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PI <- list()
  }
  return(msobject)
}

# idPSneg
#' Phosphoserines (PS) annotation for ESI-
#'
#' PS identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for PS in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idPSneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate PS as M-H or M+Na-2H. 2) Search of PS class
#' fragments: neutral loss of 87.032 (serine) coeluting with the precursor ion.
#' 3) Search of specific fragments that inform about chain composition at sn1
#' (lysoPA as M-H or M-H-H2O resulting from the loss of the FA chain at sn2 and
#' the head group) and sn2 (lysoPA as M-H or M-H-H2O resulting from the loss of
#' the FA chain at sn1 and the head group or FA chain as M-H). 4) Look for
#' possible chains structure based on the combination of chain fragments.
#' 5) Check intensity rules to confirm chains position. In this case, lysoPA
#' from sn1 is at least 3 times more intense than lysoPA from sn2.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idPSneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idPSneg <- function(msobject,
                    ppm_precursor = 5,
                    ppm_products = 10,
                    rttol = 3,
                    rt,
                    adducts = c("M-H", "M+Na-2H"),
                    clfrags = c(87.032, 152.9958),
                    clrequired = c(F, F),
                    ftype = c("NL", "F"),
                    chainfrags_sn1 = c("lysopa_M-H", "lysopa_M-H-H2O"),
                    chainfrags_sn2 = c("lysopa_M-H", "lysopa_M-H-H2O", "fa_M-H"),
                    intrules = c("lysopa_sn1/lysopa_sn2"),
                    rates = c("3/1"),
                    intrequired = c(T),
                    coelCutoff = 0.8,
                    dbs,
                    verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "PS",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("PS" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious PS annotations removed")}
      msobject$annotation$detailsAnnotation$PS <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$psdb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="PS",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PS <- list()
    msobject$annotation$detailsAnnotation$PS$candidates <- candidates
    msobject$annotation$detailsAnnotation$PS$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$PS$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$PS$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$PS <- list()
  }
  return(msobject)
}

# idSphneg
#' Sphingoid bases (Sph) annotation for ESI-
#'
#' Sph identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for Sph in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idSphneg} function involves 2 steps. 1) FullMS-based
#' identification of candidate Sph as M-H. 2) Search of Sph class fragments:
#' neutral loss of 1 or 2 H2O molecules.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (in this
#' case, as Sph only have one chain, only Subclass and FA level are possible)
#' and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idSphneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idSphneg <- function(msobject,
                     ppm_precursor = 5,
                     ppm_products = 10,
                     rttol = 3, rt,
                     adducts = c("M-H"),
                     clfrags = c("sph_M-H-H2O", "sph_M-H-2H2O"),
                     clrequired = c(F, F),
                     ftype = c("BB", "BB"),
                     coelCutoff = 0.8,
                     dbs,
                     verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "Sph",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("Sph" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious Sph annotations removed")}
      msobject$annotation$detailsAnnotation$Sph <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$sphdb, ppm = ppm_precursor,
                               rt = rt, adducts = adducts, rttol = rttol,
                               dbs = dbs, rawData = rawData,
                               coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb = list(),
                           intrules  = c(), intConf = list(), nchains = 0,
                           class="Sph",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$Sph <- list()
    msobject$annotation$detailsAnnotation$Sph$candidates <- candidates
    msobject$annotation$detailsAnnotation$Sph$classfragments <- classConf$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$Sph$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$Sph <- list()
  }
  return(msobject)
}

# idSphPpneg
#' Sphingoid bases phosphate (SphP) annotation for ESI-
#'
#' SphP identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for SphP in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idSphpos} function involves 2 steps. 1) FullMS-based
#' identification of candidate SphP as M-H. 2) Search of SphP class fragments:
#' 78.9585, 96.969 or neutral loss of 1 H2O molecule.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (in this
#' case, as SphP only have one chain, only Subclass and FA level are possible)
#' and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idSphPneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idSphPneg <- function(msobject,
                      ppm_precursor = 5,
                      ppm_products = 10,
                      rttol = 3,
                      rt,
                      adducts = c("M-H"),
                      clfrags = c(78.9585, 96.9691, "sphP_M-H-H2O"),
                      clrequired = c(F, F, F),
                      ftype = c("F", "F", "BB"),
                      coelCutoff = 0.8,
                      dbs,
                      verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "SphP",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("SphP" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious SphP annotations removed")}
      msobject$annotation$detailsAnnotation$SphP <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$sphPdb, ppm = ppm_precursor,
                               rt = rt, adducts = adducts, rttol = rttol,
                               dbs = dbs, rawData = rawData,
                               coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb = list(),
                           intrules  = c(), intConf = list(), nchains = 0,
                           class="SphP",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$SphP <- list()
    msobject$annotation$detailsAnnotation$SphP$candidates <- candidates
    msobject$annotation$detailsAnnotation$SphP$classfragments <- classConf$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$SphP$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$SphP <- list()
  }
  return(msobject)
}

# idCerneg
#' Ceramides (Cer) annotation for ESI-
#'
#' Cer identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for Cer in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected ratesbetween fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idCerneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate Cer as M-H and M+CH3COO. 2) Search of Cer class
#' fragments: there are no class fragment by default. 3) Search of specific
#' fragments that inform about the sphingoid base (Sph as M-H-2H2O resulting
#' from the loss of the FA chain or loss of part of the sphingoid base) and the
#' FA chain (FA as M-H but with a N instead of an O, what means a mass difference
#' of 1.9918 from the exact mass of the FA or FA as M-H-H2O). 4) Look for 
#' possible chains structure based on the combination of chain fragments. 
#' 5) Check intensity rules to confirm chains position. In this case, there are 
#' no intensity rules by default.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idCerneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idCerneg <- function(msobject,
                     ppm_precursor = 5,
                     ppm_products = 10,
                     rttol = 3,
                     rt,
                     adducts = c("M-H","M+CH3COO"),
                     clfrags = c(),
                     clrequired = c(),
                     ftype = c(),
                     chainfrags_sn1 = c("NL-nlsph_M-H", "sph_M-H-2H2O", "sph_M-H-H2O"),
                     chainfrags_sn2 = c("fa_Mn-1.9918", "fa_M-H-H2O"),
                     intrules = c(),
                     rates = c(),
                     intrequired = c(),
                     coelCutoff = 0.8,
                     dbs,
                     verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "Cer",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("Cer" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious Ceramide annotations removed")}
      msobject$annotation$detailsAnnotation$Cer <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$cerdb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="Cer",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$Cer <- list()
    msobject$annotation$detailsAnnotation$Cer$candidates <- candidates
    msobject$annotation$detailsAnnotation$Cer$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$Cer$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$Cer$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$Cer <- list()
  }
  return(msobject)
}

# idCerPneg
#' Ceramides phosphate (CerP) annotation for ESI-
#'
#' CerP identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for CerP in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected ratesbetween fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idCerPneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate CerP as M-H. 2) Search of CerP class
#' fragments: 78.9585 and 96.9691. 3) Search of specific fragments that inform 
#' about the sphingoid base (SphP as M-H resulting from the loss of the FA chain) 
#' and the FA chain (FA as M-H but with a N instead of an O, what results in a 
#' mass difference of 1.9918 from the exact mass of the FA, or the difference 
#' between precursor and sn1 chain fragments). 4) Look for possible chains
#' structure based on the combination of chain fragments. 5) Check intensity
#' rules to confirm chains position. In this case, there are no intensity rules
#' by default.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idCerPneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idCerPneg <- function(msobject,
                     ppm_precursor = 5,
                     ppm_products = 10,
                     rttol = 3,
                     rt,
                     adducts = c("M-H"),
                     clfrags = c(78.9585, 96.9691),
                     clrequired = c(F, F),
                     ftype = c("F", "F"),
                     chainfrags_sn1 = c("sphP_M-H"),
                     chainfrags_sn2 = c("fa_Mn-1.9918", ""),
                     intrules = c(),
                     rates = c(),
                     intrequired = c(),
                     coelCutoff = 0.8,
                     dbs,
                     verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "CerP",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("CerP" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious CeramideP annotations removed")}
      msobject$annotation$detailsAnnotation$CerP <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps
  
  # candidates search
  candidates <- findCandidates(MS1, dbs$cerPdb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)
  
  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }
    
    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)
    
    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)
    
    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)
    
    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)
    
    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="CerP",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)
    
    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$CerP <- list()
    msobject$annotation$detailsAnnotation$CerP$candidates <- candidates
    msobject$annotation$detailsAnnotation$CerP$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$CerP$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$CerP$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$CerP <- list()
  }
  return(msobject)
}

# idAcylCerneg
#' Acylceramides (AcylCer) annotation for ESI-
#'
#' AcylCer identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for AcylCer in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the sphingoid base. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param chainfrags_sn3 character vector containing the fragmentation rules for
#' the acyl chain. See \link{chainFrags} for details.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected ratesbetween fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idAcylCerneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate AcylCer as M-H and M+CH3COO. 2) Search of AcylCer 
#' class fragments: no class fragments by default. 3) Search of specific fragments 
#' that inform about the acyl chain (Cer as M-H), the sphingoid base (neutral 
#' loss of 62.0600 of the Sph) and the FA chain (FA as M-H and M-H2O but with a N 
#' instead of an O, what results in a mass differences of 1.9918 and 19.0179 
#' respectively). 4) Look for possible chains structure based on the combination 
#' of chain fragments. 5) Check intensity rules to confirm chains position. In 
#' this case, the fragment coming from the loss of the acyl chain must be at least 
#' 5 times more intense the fragment from the sphingoid base and this one, two 
#' times more intense than the FA chain from sn3.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idAcylCerneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idAcylCerneg <- function(msobject,
                      ppm_precursor = 5,
                      ppm_products = 10,
                      rttol = 3,
                      rt,
                      adducts = c("M-H", "M+CH3COO"),
                      clfrags = c(),
                      clrequired = c(),
                      ftype = c(),
                      chainfrags_sn1 = c("cbdiff-cer_M-H"),
                      chainfrags_sn2 = c("sph_Mn-62.06001", "sph_M-H-H2O"),
                      chainfrags_sn3 = c("fa_Mn-1.9918", "fa_Mn-19.0179"),
                      intrules = c("cbdiff-cer_sn1/sph_sn2", "sph_sn2/fa_sn3"),
                      rates = c("5/1", "2/1"),
                      intrequired = c(T, T),
                      coelCutoff = 0.8,
                      dbs,
                      verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "AcylCer",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("AcylCer" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious AcylCer annotations removed")}
      msobject$annotation$detailsAnnotation$AcylCer <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps
  
  # candidates search
  candidates <- findCandidates(MS1, dbs$acylcerdb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)
  
  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }
    
    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)
    
    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)
    sn3 <- chainFrags(coelfrags, chainfrags_sn3, ppm_products, candidates, sn1,
                      dbs)
    
    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=3, sn1, sn2, sn3)
    
    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)
    
    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 3, class="AcylCer",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)
    
    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$AcylCer <- list()
    msobject$annotation$detailsAnnotation$AcylCer$candidates <- candidates
    msobject$annotation$detailsAnnotation$AcylCer$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$AcylCer$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$AcylCer$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$AcylCer <- list()
  }
  return(msobject)
}

# idCLneg
#' Cardiolipines (CL) annotation for ESI-
#'
#' CL identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for CL in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details.
#' @param chainfrags_sn3 character vector containing the fragmentation rules for
#' the chain fragments in sn3 position. See \link{chainFrags} for details.
#' @param chainfrags_sn4 character vector containing the fragmentation rules for
#' the chain fragments in sn4 position. See \link{chainFrags} for details.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}. If some intensity rules should be employed to
#' identify the chains position but they are't known yet, use "Unknown". If it
#' isn't required, leave an empty vector.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idCLneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate CL as M-H or M-2H. 2) Search of CL class fragments:
#' no class fragments are searched by defaults as they use to have bad coelution
#' scores. 3) Search of specific fragments that inform about chain composition
#' at sn1 (lysoPA as M-H-H2O), sn2 (lysoPA as M-H-H2O), sn3 (lysoPA as M-H-H2O)
#' and sn4 (lysoPA as M-H-H2O). 4) Look for possible chains structure based on
#' the combination of chain fragments. 5) Check intensity rules to confirm
#' chains position. For CL there are no intensity rules by default.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idCLneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idCLneg <- function(msobject,
                    ppm_precursor = 5,
                    ppm_products = 10,
                    rttol = 5,
                    rt,
                    adducts = c("M-H", "M+Na-2H"),
                    clfrags = c(),
                    clrequired = c(),
                    ftype = c(),
                    chainfrags_sn1 = c("lysopa_M-H-H2O"),
                    chainfrags_sn2 = c("lysopa_M-H-H2O"),
                    chainfrags_sn3 = c("lysopa_M-H-H2O"),
                    chainfrags_sn4 = c("lysopa_M-H-H2O"),
                    intrules = c("Unknown"),
                    rates = c(),
                    intrequired = c(),
                    coelCutoff = 0.8,
                    dbs,
                    verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "CL",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("CL" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious CL annotations removed")}
      msobject$annotation$detailsAnnotation$CL <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$cldb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)

  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)

    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)
    sn3 <- chainFrags(coelfrags, chainfrags_sn3, ppm_products, candidates, sn1,
                      dbs)
    sn4 <- chainFrags(coelfrags, chainfrags_sn4, ppm_products, candidates, sn1,
                      dbs)

    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=4, sn1, sn2, sn3, sn4)

    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=4,
                                   chainsComb)

    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 4, class="CL",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)


    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$CL <- list()
    msobject$annotation$detailsAnnotation$CL$candidates <- candidates
    msobject$annotation$detailsAnnotation$CL$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$CL$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$CL$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$CL <- list()
  }
  return(msobject)
}

# idBAneg
#' Bile Acids (BA) annotation for ESI-
#'
#' BA identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for BA in ESI-. Adducts allowed can
#' be modified in the adducsTable (dbs argument).
#' @param conjfrag character vector containing the fragmentation rules for
#' the BA-conjugates. By default just taurine and glycine are considered,
#' but baconjdb can be modified to add more possible conjugates.
#' See \link{chainFrags} for details. It can also be an empty vector.
#' @param bafrag character vector containing the fragmentation rules for
#' other BA fragments. See \link{chainFrags} for details. It can be an empty
#' vector.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idBAneg} function involves 3 steps. 1) FullMS-based
#' identification of candidate BA as M-H. 2) Search of BA-conjugate fragments if
#' required. 3) Search of fragments coming from the loss of H2O.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (MS-only
#' if no rules are defined, or Subclass level if they are supported by fragments)
#' and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idBAneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idBAneg <- function(msobject,
                    ppm_precursor = 5,
                    ppm_products = 10,
                    rttol = 3,
                    rt,
                    adducts = c("M-H"),
                    conjfrag = c("baconj_M-H"),
                    bafrag = c("ba_M-H-H2O", "ba_M-H-2H2O"),
                    coelCutoff = 0.8,
                    dbs,
                    verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "BA",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("BA" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious BA annotations removed")}
      msobject$annotation$detailsAnnotation$BA <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps

  # candidates search
  candidates <- findCandidates(MS1, dbs$badb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)


  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }

    # search fragments
    check <- rep(list(vector()), nrow(candidates))
    conjugates <- rep(list(vector()), nrow(candidates))
    bas <- rep(list(vector()), nrow(candidates))
    if (length(conjfrag) > 0){
      conjugates <- chainFrags(coelfrags, conjfrag, ppm_products, dbs = dbs,
                               candidates = candidates)
      for (c in 1:nrow(candidates)){
        if (nrow(conjugates[[c]]) > 0){
          if (dbs$badb[dbs$badb$total ==
                       candidates$cb[c], "conjugate"] %in%
              conjugates[[c]]$cb | dbs$badb[dbs$badb$total ==
                                            candidates$cb[c],
                                            "conjugate"] == ""){
            check[[c]] <- append(check[[c]], TRUE)
          } else {
            check[[c]] <- append(check[[c]], FALSE)
          }
        } else if (dbs$badb[dbs$badb$total ==
                            candidates$cb[c],
                            "conjugate"] == ""){
          check[[c]] <- append(check[[c]], TRUE)
          conjugates[[c]] <- data.frame(0,0,0,0,0,0,0,0)
          colnames(conjugates[[c]]) <- c("cb", "mz", "RT", "int",
                                         "peakID", "coelScore", "db", "adduct")
        } else {
          check[[c]] <- append(check[[c]], FALSE)
          conjugates[[c]] <- data.frame(0,0,0,0,0,0,0,0)
          colnames(conjugates[[c]]) <- c("cb", "mz", "RT", "int",
                                         "peakID", "coelScore", "db", "adduct")
        }
      }
      classConf <- list()
    }
    if (length(bafrag) > 0){
      bas <- chainFrags(coelfrags, bafrag, ppm_products, dbs = dbs,
                        candidates = candidates)
      for (c in 1:nrow(candidates)){
        if(nrow(bas[[c]]) > 0){
          if (candidates$cb[c] %in% bas[[c]]$cb ||
              dbs$badb[dbs$badb$total == candidates$cb[c], "base"]
              %in% bas[[c]]$cb){
            check[[c]] <- append(check[[c]], TRUE)
          } else {
            check[[c]] <- append(check[[c]], FALSE)
          }
        } else {
          check[[c]] <- append(check[[c]], FALSE)
          bas[[c]] <- data.frame(0,0,0,0,0,NA,0,0)
          colnames(bas[[c]]) <- c("cb", "mz", "RT", "int",
                                  "peakID", "coelScore", "db", "adduct")
        }
      }
    }
    check <- matrix(unlist(check), nrow = nrow(candidates), byrow = T)
    classConf$presence <- check
    classConf$passed <- unlist(apply(check, 1, sum)) > 0
    classConf$fragments <- Map(rbind, conjugates, bas)
    classConf$fragments <- lapply(classConf$fragments, function(x){
      x[x$peakID != "0",,drop = FALSE]
    })

    
    # prepare output
    if (length(bafrag) > 0 | length(conjfrag) > 0){
      clfrags <- c("fragment")
    } else {
      clfrags <- c()
    }

    res <- organizeResults(candidates, clfrags = clfrags, classConf = classConf,
                           chainsComb = c(), intrules = c(),
                           intConf = c(), nchains = 0, class="BA",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)

    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$BA <- list()
    msobject$annotation$detailsAnnotation$BA$candidates <- candidates
    msobject$annotation$detailsAnnotation$BA$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$BA$chainfragments <- bas$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$BA$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$BA <- list()
  }
  return(msobject)
}

# idSMneg
#' Sphingomyelines (SM) annotation for ESI-
#'
#' SM identification based on fragmentation patterns for LC-MS/MS DIA or DDA
#' data acquired in negative mode.
#'
#' @param msobject an msobject returned by \link{dataProcessing}.
#' @param ppm_precursor mass tolerance for precursor ions. By default, 5 ppm.
#' @param ppm_products mass tolerance for product ions. By default, 10 ppm.
#' @param rttol total rt window for coelution between precursor and product
#' ions. By default, 3 seconds.
#' @param rt rt range where the function will look for candidates. By default,
#' it will search within all RT range in MS1.
#' @param adducts expected adducts for PC in ESI-. Adducts allowed can
#' be modified in adductsTable (dbs argument).
#' @param clfrags vector containing the expected fragments for a given lipid
#' class. See \link{checkClass} for details.
#' @param ftype character vector indicating the type of fragments in clfrags.
#' It can be: "F" (fragment), "NL" (neutral loss) or "BB" (building block).
#' See \link{checkClass} for details.
#' @param clrequired logical vector indicating if each class fragment is
#' required or not. If any of them is required, at least one of them must be
#' present within the coeluting fragments. See \link{checkClass} for details.
#' @param chainfrags_sn1 character vector containing the fragmentation rules for
#' the chain fragments in sn1 position. See \link{chainFrags} for details.
#' @param chainfrags_sn2 character vector containing the fragmentation rules for
#' the chain fragments in sn2 position. See \link{chainFrags} for details. If
#' empty, it will be estimated based on the difference between precursors and
#' sn1 chains.
#' @param intrules character vector specifying the fragments to compare. See
#' \link{checkIntensityRules}.
#' @param rates character vector with the expected rates between fragments given
#' as a string (e.g. "3/1"). See \link{checkIntensityRules}.
#' @param intrequired logical vector indicating if any of the rules is required.
#' If not, at least one must be verified to confirm the structure.
#' @param coelCutoff coelution score threshold between parent and fragment ions.
#' Only applied if rawData info is supplied. By default, 0.8.
#' @param dbs list of data bases required for annotation. By default, dbs
#' contains the required data frames based on the default fragmentation rules.
#' If these rules are modified, dbs may need to be supplied. See \link{createLipidDB}
#' and \link{assignDB}.
#' @param verbose print information messages.
#'
#' @return annotated msobject (list with several elements). The results element
#' is a data frame that shows: ID, lipid class, CDB (total number of carbons
#' and double bounds), FA composition (specific chains composition if it has
#' been confirmed), mz, RT (in seconds), I (intensity), Adducts, ppm (mz error),
#' confidenceLevel (Subclass, FA level, where chains are known but not their
#' positions, or FA position level), peakID, and Score (parent-fragment coelution 
#' score mean in DIA data or relative sum intensity in DDA of all fragments used 
#' for the identification).
#'
#' @details \code{idSMneg} function involves 5 steps. 1) FullMS-based
#' identification of candidate SM as M+CH3COO, M-CH3 or M+CH3COO-CH3. 2) Search 
#' of SM class fragments: 168.0426, 224.0688 or loss of CH3 coeluting with the 
#' precursor ion. 3) Search of specific fragments that inform about chain 
#' composition in sn1 (Sph+phosphocholine as M-CH3-H2O which results in a mass 
#' difference of Sph+150.032) and sn2 (difference between precursor and sn1 chain
#' fragments). 4) Look for possible chains structure based on the
#' combination of chain fragments. 5) Check intensity rules to confirm chains 
#' position. In this case, there are no intensity rules by default.
#'
#' Results data frame shows: ID, lipid class, CDB (total number
#' of carbons and double bounds), FA composition (specific chains composition if
#' it has been confirmed), mz, RT (in seconds), I (intensity, which comes
#' directly from de input), Adducts, ppm (mz error), confidenceLevel (Subclass,
#' FA level, where chains are known but not their positions, or FA position
#' level) and Score (parent-fragment coelution score mean in DIA data or relative 
#' sum intensity in DDA of all fragments used for the identification).
#'
#' @note This function has been writen based on fragmentation patterns
#' observed for three different platforms (QTOF 6550 from Agilent, Synapt G2-Si
#' from Waters and Q-exactive from Thermo), but it may need to be customized for
#' other platforms or acquisition settings.
#'
#' @examples
#' \dontrun{
#' msobject <- idSMneg(msobject)
#' }
#'
#' @author M Isabel Alcoriza-Balaguer <maialba@alumni.uv.es>
idSMneg <- function(msobject,
                    ppm_precursor = 5,
                    ppm_products = 10,
                    rttol = 3,
                    rt,
                    adducts = c("M+CH3COO", "M-CH3", "M+CH3COO-CH3"),
                    clfrags = c(168.0426, 224.0688, "sm_M-CH3"),
                    clrequired = c(F, F, F),
                    ftype = c("F", "F", "BB"),
                    chainfrags_sn1 = c("sph_Mn+150.032"),
                    chainfrags_sn2 = c("fa_Mn-1.9918", ""),
                    intrules = c(),
                    rates = c(),
                    intrequired = c(),
                    coelCutoff = 0.8,
                    dbs,
                    verbose = TRUE){
  ##############################################################################
  # check arguments
  if (msobject$metaData$generalMetadata$polarity != "negative"){
    stop("Data wasn't acquired in negative mode")
  }
  if (missing(dbs)){
    dbs <- assignDB()
  }
  if (!all(c("metaData", "processing", "rawData", "peaklist") %in% names(msobject))){
    stop("Wrong msobject format")
  }
  if (!all(c("MS1", "MS2") %in% names(msobject$rawData))){
    stop("Wrong msobject format")
  }
  if (!msobject$metaData$generalMetadata$acquisitionmode %in% c("DIA", "DDA")){
    stop("Acquisition mode must be DIA or DDA")
  }
  if (!all(adducts %in% dbs[["adductsTable"]]$adduct)){
    stop("Some adducts can't be found at the adductsTable. Add them.")
  }
  if (length(clfrags) > 0){
    if (length(clfrags) != length(clrequired) | length(clfrags) !=
        length(ftype)){
      stop("clfrags, clrequired and ftype should have the same length")
    }
    if (!all(ftype %in% c("F", "NL", "BB"))){
      stop("ftype values allowed are: \"F\", \"NL\" or\"BB\"")
    }
    strfrag <- which(grepl("_", clfrags))
    if (length(strfrag) > 0){
      d <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 1))
      a <- unlist(lapply(strsplit(clfrags[strfrag], "_"), "[[", 2))
      if (!all(a %in% dbs[["adductsTable"]]$adduct)){
        stop("Adducts employed in clfrags also need to be at adductsTable.")
      }
      if (!all(paste(d, "db", sep="") %in% names(dbs))){
        stop("All required dbs must be supplied through dbs argument.")
      }
    }
  }
  ##############################################################################
  # extract data from msobject
  # Peaklist MS1: remove isotopes
  MS1 <- msobject$peaklist$MS1
  MS1 <- MS1[MS1$isotope %in% c("[M+0]"),
             !colnames(MS1) %in% c("isotope", "isoGroup")]
  # Peaklist MS2:
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    MS2 <- msobject$rawData$MS2[,c("mz", "RT", "int", "peakID")]
  } else {
    MS2 <- msobject$peaklist$MS2[,c("mz", "RT", "int", "peakID")]
  }
  rawData <- rbind(msobject$rawData$MS1, msobject$rawData$MS2)
  # if acquisition mode is DDA, extract precursors
  if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
    precursors <- msobject$metaData$scansMetadata[msobject$metaData$scansMetadata$collisionEnergy > 0 &
                                                    msobject$metaData$scansMetadata$msLevel == 2,
                                                  c("RT", "precursor", "Scan")]
  }
  ##############################################################################
  # Remove previous ceramide annotations
  if ("results" %in% names(msobject$annotation)){
    if (nrow(msobject$annotation$results) > 0){
      msobject$annotation$results <- msobject$annotation$results[msobject$annotation$results$Class != "SM",]
    }
  }
  if ("detailsAnnotation" %in% names(msobject$annotation)){
    if("SM" %in% names(msobject$annotation$detailsAnnotation)){
      if(verbose){cat("\nPrevious SM annotations removed")}
      msobject$annotation$detailsAnnotation$SM <- list()
    }
  }
  ##############################################################################
  # set rt limits
  if (missing(rt)){
    rt <- c(min(MS1$RT), max(MS1$RT))
  }
  ##############################################################################
  # Start identification steps
  
  # candidates search
  candidates <- findCandidates(MS1, dbs$smdb, ppm = ppm_precursor, rt = rt,
                               adducts = adducts, rttol = rttol, dbs = dbs,
                               rawData = rawData, coelCutoff = coelCutoff)
  
  if (nrow(candidates) > 0){
    if (msobject$metaData$generalMetadata$acquisitionmode == "DIA"){
      if (nrow(rawData) == 0){
        coelCutoff <- 0 # if no rawData is supplied, coelution score between precursors and fragments will be ignored
      }
      # isolation of coeluting fragments
      coelfrags <- coelutingFrags(candidates, MS2, rttol, rawData,
                                  coelCutoff = coelCutoff)
    } else if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      coelCutoff <- 0
      coelfrags <- ddaFrags(candidates, precursors, rawData, ppm = ppm_products)
    }
    
    # check class fragments
    classConf <- checkClass(candidates, coelfrags, clfrags, ftype, clrequired,
                            ppm_products, dbs)
    
    # search chains fragments
    sn1 <- chainFrags(coelfrags, chainfrags_sn1, ppm_products, dbs = dbs,
                      candidates = candidates)
    sn2 <- chainFrags(coelfrags, chainfrags_sn2, ppm_products, candidates, sn1,
                      dbs)
    
    # combine chain fragments
    chainsComb <- combineChains(candidates, nchains=2, sn1, sn2)
    
    # check chains position based on intensity ratios
    intConf <- checkIntensityRules(intrules, rates, intrequired, nchains=2,
                                   chainsComb)
    
    # prepare output
    res <- organizeResults(candidates, clfrags, classConf, chainsComb, intrules,
                           intConf, nchains = 2, class="SM",
                           acquisitionmode = msobject$metaData$generalMetadata$acquisitionmode)
    
    # update msobject
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$SM <- list()
    msobject$annotation$detailsAnnotation$SM$candidates <- candidates
    msobject$annotation$detailsAnnotation$SM$classfragments <- classConf$fragments
    msobject$annotation$detailsAnnotation$SM$chainfragments <- chainsComb$fragments
    if (msobject$metaData$generalMetadata$acquisitionmode == "DDA"){
      msobject$annotation$detailsAnnotation$SM$coelfrags <- coelfrags
    }
  } else {
    res <- data.frame()
    if ("results" %in% names(msobject$annotation)){
      msobject$annotation$results <- rbind(msobject$annotation$results, res)
    } else {
      msobject$annotation$results <- res
    }
    msobject$annotation$detailsAnnotation$SM <- list()
  }
  return(msobject)
}

Try the LipidMS package in your browser

Any scripts or data that you put into this service are public.

LipidMS documentation built on March 18, 2022, 7:14 p.m.