R/get_xic.R

In EMSL-Computing/pspecterlib: PSpecteR Library - Visualization tools for top-down and bottom-up proteomics data

#' Get the Extracted Ion Chromatogram object "xic"
#'
#' @description Pulls the intensity and retention times for each isotope and adjacent charge state.
#'
#' @param ScanMetadata Object of the scan_metadata class from get_scan_metadata. Required.
#' @param MZ The mz value to extract the ion chromatogram. Required.
#' @param RTRange The retention time range to pull isotope and charge values from. It
#'    should be a vector of length two with the min and max RT in minutes. Required.
#' @param IsotopeNumber A vector of isotopes to plot. 1 indicates M+1, etc.
#'    The max isotopes is 5, ranging from 0 to 5 where 0 is the non-isotope peak. Default is 0:5.
#' @param Charges A vector indicating which charges to plot. There is a max of 3,
#'   ranging from 1-3. Default is 1:3.
#' @param XICTol The tolerance for XIC, used only with raw data.
#' @param Messages Should progress messages be printed? Default is TRUE.
#'
#' @importMethodsFrom MSnbase filterMz
#'
#' @details
#' The data.table outputted by this function contains 16 main columns of data from both the MS and ID file.
#' \tabular{ll}{
#' Charge \tab The charge (Z) state for each trace of the XIC. \cr
#' \tab \cr
#' Isotope \tab Written as "M+n", the n is indicative of each successive isotope after the fragment, where the fragment is M+0 or 0. \cr
#' \tab \cr
#' Intensity \tab The intensity at the given m/z and retention time range. \cr
#' \tab \cr
#' RT \tab The retention time \cr
#' \tab \cr
#' }
#'
#' @examples
#' \dontrun{
#'
#' # Test with bottom up data
#' xictest1 <- get_xic(ScanMetadata = BU_ScanMetadata, MZ = 659.8596, RTRange = c(58, 60))
#'
#' # Test with top down raw data
#' xictest2 <- get_xic(ScanMetadata = RAW_ScanMetadata, MZ = 669.8381, RTRange = c(0, 1))
#'
#' }
#'
#' @export
get_xic <- function(ScanMetadata,
                    MZ,
                    RTRange,
                    IsotopeNumber = c(0:5),
                    Charges = c(1:3),
                    XICTol = 1000,
                    Messages = TRUE) {

  ##################
  ## CHECK INPUTS ##
  ##################

  # Assert that ScanMetadata is a ScanMetadata object.
  if ("scan_metadata" %in% class(ScanMetadata) == FALSE) {
    stop("ScanMetadata must be a scan_metadata object generated by get_scan_metadata.")
  }

  # Assert that MZ is numeric
  if (!is.numeric(MZ)) {
    stop("MZ must be a number.")
  }

  # Assert that RT is a range
  if (!is.numeric(RTRange) || length(RTRange) != 2 || RTRange[1] > RTRange[2]) {
    stop("RT range must be written as a vector with two numerics, a min and max.")
  }

  # Assert that the isotope number is in range 0:5
  if (!is.numeric(IsotopeNumber) || min(IsotopeNumber) < 0 || max(IsotopeNumber) > 5) {
    stop("The range of isotope values must be a vector of integers with values between 0 and 5.")
  } else {
    IsotopeNumber <- floor(IsotopeNumber)
  }

  # Assert that the charge is in range 1:3
  if (!is.numeric(Charges) || min(Charges) < 1 || max(Charges) > 3) {
    stop("The range of charge values must be a vector of integers between 1 and 3.")
  } else {
    Charges <- floor(Charges)
  }

  # Assert that messages is TRUE/FALSE
  if (!is.logical(Messages)) {
    stop("Messages must be TRUE or FALSE.")
  }
  if (is.na(Messages)) {Messages <- FALSE}

  ###############
  ## PULL DATA ##
  ###############

  # Determine all MZ values to check
  XIC <- do.call(rbind, lapply(Charges, function(Charge) {
    MZValue <- (MZ + IsotopeNumber) / Charge
    return(data.frame("Charge" = Charge, "Isotope" = IsotopeNumber,
                      "Min" = MZValue - ((1/Charge) / 2), "Max" = MZValue + ((1/Charge) / 2)))
  }))

  # Different methods have to be used whether the data is XML or RAW
  if (attr(ScanMetadata, "pspecter")$MSFileType == "XML") {

    # Latest MSnbase update requires a library load
    library(MSnbase)

    # Filter down data to MS1 scans, the MS range of all XIC data, and the retention time range
    MetadataForXIC <- attr(ScanMetadata, "pspecter")$MSnbase %>%
      MSnbase::filterMsLevel(msLevel = 1) %>%
      MSnbase::filterRt(rt = RTRange*60) %>%
      MSnbase::filterMz(mz = c(min(XIC$Min), max(XIC$Max)))

    # Generate the XIC data.table
    XIC_Traces <- do.call(rbind, lapply(1:nrow(XIC), function(row) {

      # Add message if enable
      if (Messages) {message(paste("...extracting data for charge", XIC$Charge[row], "and isotope", XIC$Isotope[row]))}

      # Filter Metadata down to exact data range
      Chromatogram <- MetadataForXIC %>%
        MSnbase::filterMz(mz = c(XIC$Min[row], XIC$Max[row])) %>%
        MSnbase::chromatogram(aggregationFun = "max")

      # Create a chromatogram data table
      Chromatogram <- data.table::data.table(
        "Charge" = XIC$Charge[row],
        "Isotope" = XIC$Isotope[row],
        "Intensity" = unlist(Chromatogram@.Data)[[1]]@intensity,
        "RT" = unlist(Chromatogram@.Data)[[1]]@rtime / 60
      )

      # Filter out any NAs
      Chromatogram <- Chromatogram[is.na(Chromatogram$Intensity) == FALSE,]

      # Return Chromatograms
      return(Chromatogram)

    }))

  } else {

    # Use rawrr's readChromatogram
    XICread <- rawrr::readChromatogram(
      rawfile = attr(ScanMetadata, "pspecter")$MSPath,
      mass = lapply(1:nrow(XIC), function(row) {mean(c(XIC$Min[row], XIC$Max[row]))}) %>% unlist(),
      tol = XICTol
    )

    # Make XIC data table
    XIC_Traces <- do.call(rbind, lapply(1:length(XICread), function(el) {

      # If no intensity, return NULL
      #if (is.null(XICread[[el]]$intensities)) {return(NULL)}

      # Pull chromatogram
      Chromatogram <- data.table::data.table(
        "Charge" = XIC$Charge[el],
        "Isotope" = XIC$Isotope[el],
        "Intensity" = XICread[[el]]$intensities,
        "RT" = XICread[[el]]$times
      )

      # Subset by time
      Chromatogram <- subset(Chromatogram, RT >= min(RTRange) & RT <= max(RTRange))

      return(Chromatogram)

    }))

  }

  # Remove row names
  row.names(XIC_Traces) <- 1:nrow(XIC_Traces)

  ##################
  ## BUILD OBJECT ##
  ##################

  # Add attributes of input data
  attr(XIC_Traces, "pspecter")$MSPath = attr(ScanMetadata, "pspecter")$MSPath
  attr(XIC_Traces, "pspecter")$MZ = MZ
  attr(XIC_Traces, "pspecter")$RTRange = RTRange
  attr(XIC_Traces, "pspecter")$IsotopeNumber = IsotopeNumber
  attr(XIC_Traces, "pspecter")$Charges = Charges

  # Add class information to XIC_Traces
  class(XIC_Traces) <- c(class(XIC_Traces), "xic_pspecter")

  # Return all traces
  return(XIC_Traces)

}