R/apply.sonar.calibration.R

In cbroeckl/calibrate.sonar: Functions for calibrating the quadrupole isolation range for waters sonar acquisition.

#' generate.sonar.calibration
#'
#' Create calibration parameters from a low collision energy SONAR scan - traditional or 'hybrid' Sonar
#' @details This function uses Apex3D64.exe, distributed by Waters Corp, to perform 4d peak detection in SONAR data, reads that output data, filters, if desired, and calibrates the linear relationships between sonar drift bin and m/z value.  An R-object is returned, and a .csv file exported.  These data contain details on the instrumentation and acquisition method used to collect the calibration data, to ensure that the calibration is only applied to compatible data files.  Masslynx Version, Instrument Serial number, and sonar settings must all match.
#'
#' @param raw.file full path (or path within working directory) to the .RAW file to which the calibration is to be applied  An entire directory can be input instead, and if a directory is chosen all appropriate .raw files in that directory will be calibrated.
#' @param sonar.cal either an R object generated by 'generate.sonar.calibration, or a path to a csv calibration file written by generate.sonar.calibration.
#' @param recursive logical.  if 'raw.file' is a path to a directory of raw files, should the function also look for raw files in all subdirectories? see ?list.files
#' @return returns nothing to the R environment
#' @return selected .RAW files will have been modified with a new _sonar.INF and optionally some file copying and renaming if using hybrid sonar.
#' @concept Waters SONAR
#' @concept mass spectrometry
#' @concept data-independent MS/MS
#' @author Corey Broeckling
#' @export

apply.sonar.calibration <- function(
  raw.file = "C:/Users/cbroeckl/Documents/temp/hybridSonar.raw",
  sonar.cal = my.cal,
  recursive = FALSE
) {

  if(!file.exists(raw.file)) {
    stop("raw file: ", raw.file, "does not exist")
  }
  
  if(tolower(substring(raw.file, first = nchar(raw.file)-2, last = nchar(raw.file))) != 'raw') {
    raw.file <- list.files(raw.file, pattern = ".raw", ignore.case = TRUE, recursive = recursive)
  }
  
  if(!is.data.frame(sonar.cal)) {
    if(file.exists(sonar.cal)) {
      sonar.cal <- read.csv(sonar.cal)
    }
  }
  
  needed.names <- c("ms.model", "ms.serial.number", "mass.lynx.version", "mass.lynx.scn", 
                    "interscan.delay", "pusher.period", "scan.time", "q.start", "q.stop", "q.width", "intercept.start", 
                    "intercept.mid", "intercept.stop", "slope.start", "slope.mid", 
                    "slope.stop")
  
  if(!is.data.frame(sonar.cal)) stop(cat("sonar.cal must be a dataframe with column names", paste(needed.names, collapse = " "), '\n'))
  if(nrow(sonar.cal) == 0) stop("the provided sonar.cal has zero rows (and zero calibrations)", '\n')
  
  ## SONAR INF file to be exported.
  ## file named _sonar.INF
  sonar.inf <-
    c("[VERSION#1]",
      "MassLynx v4.2 SCN1018",
      "", "[DATE]",
      "acquire.date",   ## must update for each output file, original "Acquired Date=dd-mmmm-yyyy"
      "<---- comment",
      "acquire.time",       ## must update for each output file, original "Acquired Time=hh:mm:ss"
      "[INSTRUMENT]",
      "Quad Resolution=15.0",
      "Collision Cell Pressure=3.0",
      "Collision Cell Gradient=5.0",
      "Ion Energy=0.2",
      "Post Filter=2",
      "Resolution Mode=Resolution",
      "Pusher Period=76.25",         ## must update for each output file
      "[SONAR#1]",
      "Scan Time=scan.time",         ## must update for each output file
      "Start M/Z=q.start",           ## must update for each output file
      "End M/Z=q.stop",             ## must update for each output file
      "Quadrupole Window=q.width",  ## must update for each output file
      "Calibration=q.start,y.low.1,y.mid.1,y.high.1",  ## must update for each output file
      "Calibration=q.mid,y.low.2,y.mid.2,y.high.2",    ## must update for each output file
      "Calibration=q.stop,y.low.3,y.mid.3,y.high.3"    ## must update for each output file
    )

  for(i in 1:length(raw.file)) {
    h <- readLines(paste0(raw.file[i], "/_HEADER.TXT"))
    ms.model <- h[grep("Instrument: ", h)]
    ms.serial.number <- strsplit(ms.model, "#", fixed = TRUE)[[1]][2]
    ms.model <- strsplit(strsplit(ms.model, "Instrument:")[[1]][2], "#", fixed = TRUE)[[1]][1]
    acquire.date <- h[grep("Acquired Date", h)]
    acquire.date <- gsub("$$ ", "", acquire.date, fixed = TRUE)
    acquire.date <- gsub(": ", "=", acquire.date, fixed = TRUE)
    acquire.time <- h[grep("Acquired Time", h)]
    acquire.time <- gsub("$$ ", "", acquire.time, fixed = TRUE)
    acquire.time <- gsub(": ", "=", acquire.time, fixed = TRUE)
    
    ## get SONAR parameters:
    method <- readLines(paste0(raw.file[i], "/_extern.inf"))
    
    f1 <- method[grep("Function 1", method):(grep("Function 2", method)-1)]
    f1.sonar <- unlist(strsplit(f1[grep("UseSONARMode", f1)], "\t", fixed = TRUE))
    f1.sonar <- as.logical(f1.sonar[length(f1.sonar)])
    f2 <- method[grep("Function 2", method):(grep("Function 3", method)-1)]
    f2.sonar <- unlist(strsplit(f2[grep("UseSONARMode", f2)], "\t", fixed = TRUE))
    f2.sonar <- as.logical(f2.sonar[length(f2.sonar)])
    
    if(!(f1.sonar | f2.sonar)) {
      next(raw.file[i], "did not use SONAR", '\n')
    }
    
    hybrid <-  !(f1.sonar && f2.sonar)
    
    mass.lynx.version <- as.numeric(strsplit(strsplit(method[2], " v")[[1]][2], " ")[[1]][1])
    mass.lynx.scn <- as.numeric(strsplit(method[2], " SCN")[[1]][2])
    # "Function Parameters - Function 2 - TOF MS FUNCTION"
    # "SONARQuadrupoleStartMass\t\t\t60.0"
    # "SONARQuadrupoleStopMass\t\t\t\t440.0"
    # "Scan Time (sec)\t\t\t\t\t0.300"
    # "SONARQuadrupolePeakWidth12"
    method <- method[grep("Function 2", method):(grep("Function 3", method)-1)]
    method <- gsub('\t', "", method)
    scan.time <- method[grep("Scan Time", method)]
    scan.time <- as.numeric(unlist(strsplit(scan.time, ")"))[2])
    q.start <- method[grep("SONARQuadrupoleStartMass", method)]
    q.start <- as.numeric(gsub("SONARQuadrupoleStartMass", "", q.start))
    q.stop <- method[grep("SONARQuadrupoleStopMass", method)]
    q.stop <- as.numeric(gsub("SONARQuadrupoleStopMass", "", q.stop))
    q.mid <- round(mean(c(q.stop, q.start)))
    q.width <- method[grep("SONARQuadrupolePeakWidth", method)]
    q.width <- as.numeric(gsub("SONARQuadrupolePeakWidth", "", q.width))
    interscan.delay <- method[grep("Interscan Time", method)]
    interscan.delay <- as.numeric(gsub("Interscan Time (sec)", "", interscan.delay, fixed = TRUE))
    pusher.period <- method[grep("ADC Pusher Frequency", method)]
    pusher.period <- as.numeric(strsplit(pusher.period, "s)")[[1]][2])
    
    current.file <- as.data.frame(t(as.matrix(c(
      "ms.model" = ms.model,
      "ms.serial.number" = ms.serial.number,
      "mass.lynx.version" = mass.lynx.version,
      "mass.lynx.scn" = mass.lynx.scn,
      "interscan.delay" = interscan.delay,
      "pusher.period" = pusher.period,
      "scan.time" = scan.time,
      "q.start" = round(q.start, digits = 1),
      "q.stop" = round(q.stop, digits = 1),
      "q.width" = round(q.width, digits = 1)
    ))))
    
    ## if hybrid sonar, we need to calibrate of function 2, which requires some hacking to ensure 
    ## that the data can be converted and read properly
    if(hybrid) {
      # you could move to _func003.cdt and _func003.ind to a separate folder
      suppressWarnings(dir.create(paste0(dirname(raw.file[i]), "/SonarCalTmp")))
      file.rename(from = paste0(raw.file[i], c("/_func003.cdt", "/_func003.ind")),
                to = paste0(raw.file[i], c("/_func003_.cdt", "/_func003_.ind")))
    }
    
    
    match.cal <- which(sapply(1:nrow(sonar.cal), FUN = function(x) identical(as.character(current.file[1,]), as.character(sonar.cal[x,names(current.file)]))))
    if(length(match.cal) == 0) {cat(basename(raw.file[i]), ": no calibration available with identical acquisition parameters to this file.", '\n')}
    if(length(match.cal) >= 1) {
      sonar.cal <- sonar.cal[order(sonar.cal[,"date.time"]),]
      use.cal <- sonar.cal[min(match.cal),]
      cat(basename(raw.file[i]), ": calibration successful", '\n')
      }
    
    ### update sonar.inf for export
    sonar.inf.out <- sonar.inf
    sonar.inf.out <- gsub("acquire.date", acquire.date, sonar.inf.out)
    sonar.inf.out <- gsub("acquire.time", acquire.time, sonar.inf.out)
    sonar.inf.out <- gsub("q.start", q.start, sonar.inf.out)
    sonar.inf.out <- gsub("q.stop", q.stop, sonar.inf.out)
    sonar.inf.out <- gsub("q.width", q.width, sonar.inf.out)
    sonar.inf.out <- gsub("q.width", q.width, sonar.inf.out)
    sonar.inf.out <- gsub("q.mid", q.mid, sonar.inf.out)
    sonar.inf.out <- gsub("scan.time", scan.time, sonar.inf.out)
    
    ### calculate calibration values
    for(j in c(
      "scan.time", "q.start", "q.stop", "intercept.start", "intercept.mid",
      "intercept.stop", "slope.start", "slope.mid", "slope.stop")) {
      use.cal[,j] <- as.numeric(use.cal[,j])
    }
    
    ## quad low end cal
    y.low.1 <- round(use.cal$slope.start*q.start + use.cal$intercept.start, digits = 4)
    y.mid.1 <- round(use.cal$slope.mid*q.start + use.cal$intercept.mid, digits = 4)
    y.high.1 <- round(use.cal$slope.stop*q.start + use.cal$intercept.stop, digits = 4)
    
    ## quad mid end cal
    y.low.2 <- round(use.cal$slope.start*q.mid + use.cal$intercept.start, digits = 4)
    y.mid.2 <- round(use.cal$slope.mid*q.mid + use.cal$intercept.mid, digits = 4)
    y.high.2 <- round(use.cal$slope.stop*q.mid + use.cal$intercept.stop, digits = 4)
    
    ## quad high end cal
    y.low.3 <- round(use.cal$slope.start*q.stop + use.cal$intercept.start, digits = 4)
    y.mid.3 <- round(use.cal$slope.mid*q.stop + use.cal$intercept.mid, digits = 4)
    y.high.3 <- round(use.cal$slope.stop*q.stop + use.cal$intercept.stop, digits = 4)
    
    ## insert calculated y values into sonar.inf.out
    sonar.inf.out <- gsub("y.low.1", y.low.1, sonar.inf.out)
    sonar.inf.out <- gsub("y.mid.1", y.mid.1, sonar.inf.out)
    sonar.inf.out <- gsub("y.high.1", y.high.1, sonar.inf.out)
    sonar.inf.out <- gsub("y.low.2", y.low.2, sonar.inf.out)
    sonar.inf.out <- gsub("y.mid.2", y.mid.2, sonar.inf.out)
    sonar.inf.out <- gsub("y.high.2", y.high.2, sonar.inf.out)
    sonar.inf.out <- gsub("y.low.3", y.low.3, sonar.inf.out)
    sonar.inf.out <- gsub("y.mid.3", y.mid.3, sonar.inf.out)
    sonar.inf.out <- gsub("y.high.3", y.high.3, sonar.inf.out)
    
    sink(paste0(raw.file, "/_sonar.INF"))
    cat(paste(sonar.inf.out, collapse = '\n'))
    sink()
    
  }
}