R/CorMID.R

In CorMID: Correct Mass Isotopologue Distribution Vectors

#'@title CorMID.
#'@description \code{CorMID} will compute a MID (Mass Isotopomer Distribution)
#'  based on measured ion intensities in GC-APCI-MS.
#'@details Let's assume we measured the ion intensities of all 3 isotopes of an
#'  individual compound containing 2 carbons and observe a vector of {978,22,0}.
#'  We may calculate the enrichment \strong{E} out of this data, i.e. the relative proportion
#'  of 13C vs total carbon which will amount to about 1.1\% (the natural 13C abundance)
#'  under standard conditions.
#'  The equivalent corMID vector would be {1,0,0}, indicating that the non-labeled
#'  isotopologue (where non-labeled means non-labeled above the natural 1.1\%)
#'  is the only component observed.
#'  During a labeling experiment we may change the measurement values in different
#'  ways (either labeling only one carbon or both), which potentially can translate
#'  into similar values for \strong{E} being larger 1.1\%.
#'  The MIDs will provide additional information about the isotopologue fraction
#'  which gave rise to the observed \strong{E}'s (cf. examples). The \emph{r} parameter
#'  indicates an overlay of chemical rearrangements which may occur.
#'@param int Named numeric vector of measured ion intensities of a fragment. Names will give position of values relative to M+H (see details).
#'@param fml Chemical formula of the fragment as string.
#'@param r Either a character vector giving fragments to be considered
#'  OR a named numeric giving relative amounts of fragments
#'  OR NULL (all known fragments will  be estimated)
#'  OR a 2-row matrix giving the lower and upper allowed ratio (see examples).
#'@param penalize Numeric exponent penalizing solutions with low M+H occurrence. Formula is 1+3*(1-x)^penalty. NA to omit penalizing.
#'@param mid_fix May provide a numeric vector used as a given MID. Allows to estimate \emph{r} individually.
#'@param trace_steps For testing purposes. Print the results of intermediate steps to console.
#'@param prec Precision of the estimation of MID, set to 1\% as default.
#'@return Estimated percent representation of each isotopologue measured (corMID).
#'@references <doi:10.3390/metabo12050408>
#'@examples
#'# make up some fake measurement data for Pyruvic acid 2TMS with 3 biological carbon
#'# assuming 10% labeling at M3 and 2 fragments
#'fml <- "C9H20O3Si2"
#'mid <- c(0.9,0,0,0.1)
#'r <- unlist(list("M+H"=0.8, "M+H2O-CH4"=0.2))
#'int <- CorMID::recMID(mid=mid, r=r, fml=fml)
#'plot(int)
#'
#'# full estimation of M and r
#'out <- CorMID::CorMID(int=int, fml=fml)
#'plot(out)
#'
#'# get an improved result setting r to the correct values
#'CorMID::CorMID(int=int, fml=fml, r=r, prec=0.0001)
#'
#'# provoke a wrong estimation using a fixed r
#'CorMID::CorMID(int=int, fml=fml, r=unlist(list("M+H"=1)))
#'
#'# calculate r if you know the true corMID for a compound
#'r <- attr(CorMID::CorMID(int=int, fml=fml, mid_fix=c(0.9,0,0,0.1)), "ratio")
#'round(CorMID::CorMID(int=int, fml=fml, r=r, prec=0.0001),3)
#'
#'# deal with missing intensity values
#'CorMID::CorMID(int=int[-3], fml=fml)
#'
#'\donttest{
#'# perform estimation with banded r and observation of optimization steps
#'r <- matrix(c(0.5,1,0,0.5,0,0.5), nrow=2, dimnames=list(NULL,c("M+H","M+","M+H2O-CH4")))
#'CorMID::CorMID(int=int, fml=fml, r=r, trace=TRUE)
#'
#'#process Gln data from publication
#'utils::data("prep", package = "CorMID")
#'int <- prep[[24]][["int"]][,6]
#'fml <- prep[[24]]$fml
#'CorMID::CorMID(int=int, fml=fml, trace=TRUE)
#'
#'# check the effect of the penalize parameter on selection of adducts
#'int <- c(1560, 119203, 41927, 16932, 4438)
#'names(int) <- c(-2, 0, 1, 2, 3)
#'fml <- "C19H37NO4Si3"
#'CorMID::CorMID(int=int, fml=fml, r=NULL, trace=TRUE)
#'CorMID::CorMID(int=int, fml=fml, r=NULL, trace=TRUE, penalize=7)
#'}
#'@export
CorMID <- function(int=NULL, fml="", r=NULL, penalize=7, mid_fix=NULL, trace_steps=FALSE, prec = 0.01) {

  # potential parameters
  # specify known fragments; I use this named vector to check in the fitting function where a certain fragment starts,
  # to make the function really flexible this would need to be user definable potentially in a data.frame together with parameter 'r'
  known_frags <- unlist(list("M+H"=0,"M+"=-1,"M-H"=-2,"M+H2O-CH4"=+2))
  # this is the maximum number of carbons computable (otherwise matrices get to big)
  lim_nbio <- min(12, attr(fml, "nmz"))

  # QC for fml
  # we need to know the biological carbon within fml and specify nmz based on our known frags
  if (is.null(attr(fml, "nbio"))) {
    cce <- CountChemicalElements(x=fml, ele=c("C","Si"))
    attr(fml, "nbio") <- unname(cce["C"]-3*cce["Si"])
  }
  attr(fml, "nbio") <- min(lim_nbio, attr(fml, "nbio"))
  attr(fml, "nmz") <- attr(fml, "nbio")+diff(range(known_frags))

  # QC for r
  # if r is unspecified
  if (is.null(r)) r <- matrix(rep(c(0,1), length(known_frags)), nrow=2, dimnames=list(c("low","high"), names(known_frags)))
  # if r specifies only fragments
  if (is.character(r) && all(r %in% names(known_frags))) r <- matrix(rep(c(0,1),length(r)), nrow=2, dimnames=list(c("low","high"), r))
  # convert r into matrix if numeric vector is provided
  if (is.null(nrow(r)) || nrow(r)==1) {
    r <- r/sum(r) # ensure that sum(r)==1
    r <- matrix(c(r,r), byrow=T, nrow=2, dimnames=list(c("low","high"),names(r)))
  }

  # QC for int and frag
  frag <- min(known_frags):(max(known_frags)+attr(fml, "nbio"))
  rawMID <- rep(0, length(frag))
  names(rawMID) <-  paste0("M", formatC(x = frag, format = "d", flag = "+"))

  # set non-finite intensity values to zero
  if (any(!is.finite(int))) int[!is.finite(int)] <- 0
  # adjust names of int to match rawMID
  if (is.null(names(int))) {
    names(int) <- paste0("M", formatC(x = 0:(length(int)-1), format = "d", flag = "+"))
  } else {
    tmp <- as.numeric(gsub("M", "", names(int)))
    if (all(is.finite(tmp))) {
      names(int) <- paste0("M", formatC(x = tmp, format = "d", flag = "+"))
    }
  }
  if (!all(names(int) %in% names(rawMID))) stop("rawMID specified without names indicating position relative to [M+H].")
  rawMID[names(int)] <- int
  # ensure that intensity vector is normalized to sum
  rawMID <- rawMID/sum(rawMID)

  # compute theoretical distribution matrix from formula (assuming 1% 13C abundance)
  td <- CalcTheoreticalMDV(fml=fml, nbio=attr(fml, "nbio"), nmz=attr(fml, "nmz"))

  # QC for mid_fix
  if (!is.null(mid_fix)) {
    if (length(mid_fix)>nrow(td)) mid_fix <- mid_fix[1:nrow(td)]
    if (length(mid_fix)<nrow(td)) warning("length(mid_fix)<nrow(td)")
    names(mid_fix) <- rownames(td)
    mid_fix <- mid_fix/sum(mid_fix)
  }

  # fitted return value
  out <- FitMID(md=rawMID, td=td, r=r, mid_fix=mid_fix, prec=prec, trace_steps=trace_steps, penalize=penalize)

  # set the class to allow dedicated plotting and printing
  class(out) <- "CorMID"

  # return relative representation of each isotopologue measured (= corrected MIDs in %)
  return(out)
}

#'@rdname CorMID
#'@param x Object of class CorMID.
#'@param ... Further plotting parameters.
#'@importFrom graphics axis par
#'@importFrom grDevices grey
#'@export
#'@method plot CorMID
plot.CorMID <- function(x, ...) {
  stopifnot(is.numeric(x))
  stopifnot(all(names(x) %in% paste0("M", 0:12)))
  argg <- c(as.list(environment()), list(...))
  if (!"xlab" %in% names(argg)) xlab <- "" else xlab <- argg$xlab
  if (!"ylab" %in% names(argg)) ylab <- "Relative intensity [%]" else ylab <- argg$ylab
  if (!"lwd" %in% names(argg)) lwd <- 7 else lwd <- argg$lwd
  if (!"lend" %in% names(argg)) lend <- 1 else lend <- argg$lend
  if (!"xlim" %in% names(argg)) xlim <- c(0.5,length(x)+0.5) else xlim <- argg$xlim
  if (!"ylim" %in% names(argg)) ylim <- c(0,100) else ylim <- argg$ylim
  if (!"las" %in% names(argg)) las <- 2 else las <- argg$las
  default_cols <- c(1:7, rep(grDevices::grey(0.9), 6))
  names(default_cols) <- paste0("M", 0:12)
  col <- default_cols[names(x)]
  if ("col" %in% names(argg)) { col <- argg$col }
  opar <- graphics::par(no.readonly = TRUE)
  on.exit(graphics::par("mar"=opar$mar))
  par(mar=c(ifelse(xlab=="", 3, 5), 4, 0, 0) + 0.3)
  plot(as.numeric(x), type="h", axes=FALSE, lwd=lwd, lend=lend, col=col, xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim)
  graphics::axis(1, at=1:length(x), labels=names(x), las=las)
  graphics::axis(2, las=las)
}

#'@rdname CorMID
#'@param x Object of class CorMID.
#'@export
#'@method print CorMID
print.CorMID <- function(x, ...) {
  stopifnot(is.numeric(x))
  # $JL$ I had implemented a colored version of the print to console
  # however, I decided to keep it in standard message style for time being
  cat_or_message <- function(txt) {
    if (interactive() & FALSE) { cat(paste0("\033[0;34m", txt, "\033[0m", "\n")) } else { message(txt) }
  }
  cat_or_message("[class] 'CorMID'")
  cat_or_message(paste0("MID [%] (", attr(x, "mid_status"), ")"))
  cat("    ", paste(names(x), collapse="    "), "\n", sep="")
  cat(" ", paste(formatC(x, digits=2, format="f", width=5, flag="0"), collapse=" "), "\n", sep="")
  cat_or_message(paste0("[attr] 'r' (", attr(x, "ratio_status"), ")"))
  r <- attr(x, "ratio")
  cat("  ", paste(names(r), collapse="   "), "\n", sep="")
  cat(sapply(1:length(r), function(i) { formatC(r[i], digits=2, format="f", width=2+nchar(names(r)[i])) }), "\n")
  cat_or_message("[attr] 'err'")
  cat(formatC(attr(x, "err"), format="g"),"\n\n")
}