add first implementation of requantification; see #39

DESCRIPTION

@@ -10,7 +10,7 @@ Maintainer: Laurent Gatto <lg390@cam.ac.uk> and
 Depends: R (>= 3.1.0), methods, MSnbase (>= 1.13.1)
 Imports: hwriter, RColorBrewer, lattice, qvalue, multtest, utils, tools,
         Biobase, knitr, Biostrings, cleaver (>= 1.3.3)
-Suggests: synapterdata, xtable, testthat (>= 0.8), tcltk, tcltk2
+Suggests: synapterdata, xtable, testthat (>= 0.8), tcltk, tcltk2, BRAIN
 Description: The synapter package provides functionality to reanalyse
         label-free proteomics data acquired on a Synapt G2 mass
         spectrometer. One or several runs, possibly processed with

R/requantify-functions.R

@@ -0,0 +1,80 @@
+requantify <- function(msnset, saturationThreshold,
+                       method=c("peptides", "th.mean", "th.median", "th.weighted.mean")) {
+  cn <- fvarLabels(msnset)
+  i <- grep("isotopicDistr", cn)
+
+  if (!length(i)) {
+    stop("Could not find any isotopic distribution information.")
+  }
+
+  f <- fData(msnset)[, i]
+  e <- exprs(msnset)
+
+  method <- match.arg(method)
+
+  if (method == "peptides") {
+    e <- t(apply(f, 1, .requantifyPeptides,
+                 saturationThreshold=saturationThreshold))
+  } else {
+    if (!requireNamespace("BRAIN")) {
+      stop("Please install the BRAIN package via 'biocLite(\"BRAIN\")' to use this method.")
+    }
+    ## remove "th." prefix
+    method <- substring(method, 4, nchar(method))
+    rn <- rownames(e)
+
+    for (i in 1:nrow(e)) {
+      e[i, ] <- .requantifyTheoreticalDistribution(f[i, ], sequence=rn[i],
+                                                   saturationThreshold=saturationThreshold,
+                                                   method=method)
+    }
+  }
+
+  dimnames(e) <- dimnames(exprs(msnset))
+  exprs(msnset) <- e
+  msnset
+}
+
+.requantifyPeptides <- function(x, saturationThreshold=Inf) {
+  x <- .splitIsotopicDistr(unlist(x))
+  commonnm <- .commonIsotopes(x, saturationThreshold)
+
+  if (length(commonnm)) {
+    return(sapply(x, function(x)sum(x[commonnm])))
+  }
+
+  rep.int(NA, length(x))
+}
+
+.requantifyTheoreticalDistribution <- function(x, sequence,
+                                               saturationThreshold=Inf,
+                                               method=c("mean", "median", "weighted.mean")) {
+  x <- .splitIsotopicDistr(unlist(x))
+  commonnm <- .commonIsotopes(x, saturationThreshold)
+
+  if (length(commonnm)) {
+    distr <- lapply(x, function(y).sumIsotopes(y[commonnm]))
+    n <- max(sapply(distr, length))
+    props <- BRAIN::calculateIsotopicProbabilities(BRAIN::getAtomsFromSeq(sequence),
+                                                   nrPeaks=n)
+    distr <- lapply(distr, function(y)y/props)
+
+    return(switch(match.arg(method),
+                  "mean" = sapply(distr, mean, na.rm=TRUE),
+                  "median" = sapply(distr, median, na.rm=TRUE),
+                  "weighted.mean" = sapply(distr, weighted.mean, w=props, na.rm=TRUE)))
+  }
+
+  rep.int(NA, length(x))
+}
+
+.commonIsotopes <- function(x, saturationThreshold=Inf) {
+  distr <- lapply(x, function(y)y[y<saturationThreshold])
+  nm <- lapply(distr, names)
+  Reduce(intersect, nm)
+}
+
+.sumIsotopes <- function(x) {
+  nm <- factor(sapply(strsplit(names(x), "_", fixed=TRUE), "[", 2))
+  tapply(x, nm, sum)
+}