Nothing
#' @title CorMID.
#' @description \code{CorMID} will compute a MID (Mass Isotopologues
#' 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 \ifelse{html}{\out{<sup>13</sup>C}}{\eqn{^13^C}}
#' vs total carbon which will amount to about 1.1\% (the natural
#' \ifelse{html}{\out{<sup>13</sup>C}}{\eqn{^13^C}} 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. Set to 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.
#' @param algo The algorithm used to estimate the isotopic distribution of a
#' chemical formula.
#' @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)
#' out
#' 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, algo = c("CorMID", "Rdisop")) {
# 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"))
algo <- match.arg(algo)
if (algo=="Rdisop") verify_suggested("Rdisop")
# 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))
# compute theoretical distribution matrix from formula (assuming 1% 13C abundance)
td <- CalcTheoreticalMDV(fml = fml, nbio = attr(fml, "nbio"), nmz = attr(fml, "nmz"), algo = algo)
# 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
# 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(ifelse(is.na(attr(x, "err")), NA, formatC(attr(x, "err"), format = "g")), "\n\n")
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.