R/mie_transform.R

In rogerswt/flowMie: Modeling a Detector Response to Scattering of (Multi-layer) Spherical Particles

#
# mie_transform.R
#
################################################################################
################################################################################
#                     Copyright Still Pond Cytomics LLC 2019.                 ##
#        All Rights Reserved. No part of this source code may be reproduced   ##
#            without Still Pond Cytomics' express written consent.            ##
################################################################################
################################################################################
#

#' @title Calibrate Detector Response Curve
#' @param detector A description of the detector.  Default = create_detector().
#' @param bead A description of the calibration bead used.  Default is a silica
#' bead with a diameter of 200nm.
#' @param bead_mfi The "MFI" of the signal generated by the calibration bead.  Default = 60000.
#' @param show Display a graph of detector response as a function of EV diameter.  Default = FALSE.
#' @description This function uses the signal from a calibration bead to determine
#' the gain factor for the instrument.  This factor converts the output of
#' \link{calculate_detector_response} such that it reads out in the same units
#' as reported by the detector.
#'
#' This function is a precursor to the use of the \link{mie_transform} function, which
#' calculates vesicle size and appends it to a \link[flowCore]{flowFrame}.
#' @return A list, with elements:
#' \describe{
#'   \item{gain}{The deduced gain factor for the detector}
#'   \item{lut}{A lookup table relating EV size to the corresponding detector signal}
#' }
#' @export
calibrate_ssc_conversion = function(detector = create_detector(),
                                    bead = create_SI(d = 200),
                                    bead_mfi = 60000,
                                    show = FALSE) {
  ssc_raw = calculate_detector_response(particle = bead, detector = detector)
  detector$gain = gain = bead_mfi / ssc_raw

  ssc_ev = vector('numeric')
  diameter = seq(20, 500, by = 2)
  for (i in 1:length(diameter)) {
    P <- create_EV(d = diameter[i])
    ssc_ev[i]  = calculate_detector_response(P, detector)
  }

  detector$lut = lut = data.frame(signal = ssc_ev, diameter = diameter)

  if (show) {
    plot(lut$diameter, lut$signal,
         xlab = "EV Diameter (nm)",
         ylab = "Detector Signal (MFI units)",
         log = 'y', yaxt = 'n',
         cex.lab = 1.5)

    dec.max = ceiling(log10(max(lut$signal)))
    dec.min = floor(log10(min(lut$signal)))
    axis(side = 2, at = 10^(dec.min:dec.max), labels = wadeTools::tick.labels(dec.min, dec.max))

    bead_diameter = 2 * bead@layers[[1]]@r
    ev_sig_at_bead_diameter = lut$signal[which(lut$diameter == bead_diameter)]
    segments(x0 = bead_diameter, y0 = ev_sig_at_bead_diameter, x1 = bead_diameter, y1 = 10 ^ dec.min, col = 'red')
    segments(x0 = 0, y0 = ev_sig_at_bead_diameter, x1 = bead_diameter, y1 = ev_sig_at_bead_diameter, col = 'red')

    x.text = bead_diameter + 50
    y.text = 10 ^ (dec.max - 2)
    text(x = x.text, y = y.text, labels = sprintf("Calibration Bead Diameter = %.1f nm", bead_diameter), pos = 4, cex = 1.25)
    text(x = x.text, y = y.text / 2, labels = sprintf("Calibration Bead Signal = %.1f", bead_mfi), pos = 4, cex = 1.25)
    text(x = x.text, y = y.text / 4, labels = sprintf("EV Signal @ %.1f nm = %.1f", bead_diameter, ev_sig_at_bead_diameter), pos = 4, cex = 1.25)
  }

  return(detector)
}



#' @title Compute the Mie Transform
#' @param ff A \link[flowCore]{flowFrame} containing EV data.
#' @param ssc_param The name of the scattering parameter to use for sizing.
#' **NOTE: should be the same parameter used to detect the calibration bead in \link{calibrate_ssc_conversion()}.**
#' @param transformation What transformation was used on the scattering parameter.  Ones
#' currently supported include c("none", "biexp").
#' @param detector A detector, having been calibrated by \link{calibrate_ssc_conversion()}.
#' @description Calculate the Mie Transform, which converts a scattering signal to
#' a corresponding EV diameter by means of calculated scattering from a calibration particle.
#' @details There are several limitations of this function having to do with whether
#' the scattering data have been transformed or not.  Since the Mie Transform operates in linear
#' space, any non-linear transformation must be inverted prior to calculating the Mie Transform.
#' Currently, the only supported non-linear transformation is link[wadeTools]{biexp}.
#' @return A flowFrame, with a new parameter called "Size (nm)".
#' @export
mie_transform = function(ff,
                         ssc_param = "SSC-A",
                         transformation = c("none", "biexp"),
                         detector) {

  transformation = match.arg(transformation, c("none", "biexp"))

  if (is.null(detector$lut)) {
    stop("Please calibrate this detector using calibrate_ssc_conversion()")
  }
  lut = detector$lut

  # extract the ssc parameter
  ssc = exprs(ff)[, ssc_param]

  # if data were transformed, undo it
  if (transformation == 'biexp') {
    ssc = ibx(ssc)
  }

  # handle negative values
  ssc[ssc < 0] = 0

  # handle large values
  ssc[ssc >= max(lut$signal)] = max(lut$signal)

  # pass it through the lut
  size = approx(lut, xout = ssc)$y

  # tack this onto ff as a new calculated parameter
  tmpmat = exprs(ff)
  tmpmat = cbind(tmpmat, "Size (nm)" = size)
  pdata = pData(parameters(ff))
  pdata = rbind(pdata, list("Size (nm)", "<NA>", 262144, 0, 500))
  last.pname = paste("$P", nrow(pdata), sep = "")
  rownames(pdata)[nrow(pdata)] = last.pname

  ffs = flowFrame(tmpmat, parameters = as(pdata, "AnnotatedDataFrame"))

  ffs
}