feat: add flowCore conversions for gates

CHANGELOG.md

@@ -16,6 +16,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - The `authenticate` function now supports two-factor authentication.
 - Support `byName()` for the `gateId` argument in `updateGate()`.
 - Introduce `lintr` as a linter with an Actions workflow
+- Add `fromFlowCore` to convert flowCore objects to CellEngine
+- Add `toFlowCore` to convert CellEngine objects to flowCore
 
 ### Changed
 - Fix for [confusing bulk entity retrieval](https://github.com/primitybio/cellengine-r-toolkit/issues/48)

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: cellengine
 Type: Package
 Title: 'CellEngine' API Toolkit
-Version: 0.2.0
+Version: 0.3.0
 Authors@R: c(
     person("Zach", "Bjornson", email = "[email protected]", role = c("aut", "cre")),
     person("Gerrit", "Egnew", email = "[email protected]", role = c("aut")))
@@ -12,5 +12,6 @@ Encoding: UTF-8
 LazyData: true
 Depends: R (>= 3.0.0)
 Imports: httr(>= 1.3.0), jsonlite, utils
-RoxygenNote: 7.1.1
-Suggests: httptest(>= 3.2.2), testthat, getPass, rstudioapi
+RoxygenNote: 7.1.2
+Suggests: httptest(>= 3.2.2), testthat, mockthat, getPass, rstudioapi, mockery
+Enhances: flowCore, flowDensity

NAMESPACE

@@ -24,6 +24,7 @@ export(deleteGates)
 export(deletePopulationAuto)
 export(downloadAttachment)
 export(downloadFcsFiles)
+export(fromFlowCore)
 export(generateId)
 export(getAttachments)
 export(getCompensations)
@@ -41,6 +42,7 @@ export(getScaleSets)
 export(getStatistics)
 export(setFcsFilePanel)
 export(setServer)
+export(toFlowCore)
 export(updateExperiment)
 export(updateGate)
 export(updateGateFamily)

R/annotateFcsFile.R

@@ -16,8 +16,10 @@
 #' @export
 #' @examples
 #' \dontrun{
-#' annotations <- list(list(name = "annotations 1", value = 1),
-#'   list(name = "annotation 2", value = "myValue"))
+#' annotations <- list(
+#'   list(name = "annotations 1", value = 1),
+#'   list(name = "annotation 2", value = "myValue")
+#' )
 #' annotateFcsFile(experimentId, fcsFileId, annotations)
 #' }
 annotateFcsFile <- function(experimentId, fcsFileId, annotations) {

R/convertFromFlowCore.R

@@ -0,0 +1,81 @@
+library("stats")
+
+#' Convert a flowCore object to a CellEngine object
+#'
+#' Creates a CellEngine object from a flowCore object, saving it to CellEngine.
+#'
+#' @param flowObject A flowCore object.
+#' @param experimentId The ID of the experiment in which to create the object.
+#' @param name The name to use for the CellEngine object. If not set, it will
+#' default to the identifier of the flowCore object.
+#' @param ... Other parameters passed to \code{createRectangleGate}, \code{createPolygonGate}, etc.
+#' @export
+#' @examples
+#' \dontrun{
+#' rectGate <- rectangleGate(
+#'   filterId = "Gate 1",
+#'   "FL1-H" = c(0, 12), "FL2-H" = c(0, 12)
+#' )
+#' experimentId <- "5d2f8b4b21fd0676fb3a6a8c"
+#' fromFlowCore(rectGate, experimentId, name = "my gate")
+#' }
+fromFlowCore <- function(flowObject, experimentId, name = NULL, ...) {
+  if (is.null(name)) {
+    name <- flowCore::identifier(flowObject)
+  }
+
+  switch(class(flowObject)[1],
+    "ellipsoidGate" = {
+      convertEllipsoidGate(flowObject, experimentId, name, ...)
+    },
+    "polygonGate" = {
+      convertPolygonGate(flowObject, experimentId, name, ...)
+    },
+    "rectangleGate" = {
+      convertRectangleGate(flowObject, experimentId, name, ...)
+    }
+  )
+}
+
+convertPolygonGate <- function(flowGate, experimentId, name = NULL, ...) {
+  boundaries <- flowGate@boundaries
+  createPolygonGate(
+    experimentId,
+    colnames(boundaries)[1],
+    colnames(boundaries)[2],
+    name,
+    vertices = list(boundaries),
+    ...
+  )
+}
+
+convertEllipsoidGate <- function(flowGate, experimentId, name = NULL, ...) {
+  params <- covarToParameters(flowGate@cov)
+
+  createEllipseGate(
+    experimentId,
+    colnames(flowGate@cov)[1],
+    colnames(flowGate@cov)[2],
+    name,
+    flowGate@mean[[1]],
+    flowGate@mean[[2]],
+    params$angle,
+    2 * params$major,
+    2 * params$minor,
+    ...
+  )
+}
+
+convertRectangleGate <- function(flowGate, experimentId, name = NULL, ...) {
+  createRectangleGate(
+    experimentId,
+    names(flowGate@min)[1],
+    names(flowGate@min)[2],
+    name,
+    flowGate@min[[1]],
+    flowGate@max[[1]],
+    flowGate@min[[2]],
+    flowGate@max[[2]],
+    ...
+  )
+}

R/convertToFlowCore.R

@@ -0,0 +1,94 @@
+library("stats")
+
+#' Convert a gate to flowCore
+#'
+#' Converts a CellEngine object to its flowCore analogue.
+#'
+#' @param cellengineObject The CellEngine object to be converted.
+#' @export
+toFlowCore <- function(cellengineObject) {
+  if (!requireNamespace("flowCore")) {
+    stop("This function requires the 'flowCore' package.")
+  }
+
+  if ("scales" %in% names(cellengineObject)) {
+    class_ <- "ScaleSet"
+  } else if ("type" %in% names(cellengineObject)) {
+    class_ <- cellengineObject$type
+  }
+
+  switch(class_,
+    "RectangleGate" = toFlowCoreRectangleGate(cellengineObject),
+    "EllipseGate" = toFlowCoreEllipsoidGate(cellengineObject),
+    "PolygonGate" = toFlowCorePolygonGate(cellengineObject),
+    "QuadrantGate" = stop("This gate representation is not yet implemented"),
+    "SplitGate" = stop("This gate representation is not yet implemented"),
+    "RangeGate" = stop("This gate representation is not yet implemented"),
+    "ScaleSet" = scaleSetToTransformList(cellengineObject)
+  )
+}
+
+toFlowCoreRectangleGate <- function(gate) {
+  flowCore::rectangleGate(
+    filterId = gate$name,
+    stats::setNames(
+      list(
+        c(gate$model$rectangle$x1, gate$model$rectangle$x2),
+        c(gate$model$rectangle$y1, gate$model$rectangle$y2)
+      ),
+      c(gate$xChannel, gate$yChannel)
+    )
+  )
+}
+
+toFlowCoreEllipsoidGate <- function(gate) {
+  ellipse <- gate$model$ellipse
+  x <- unlist(ellipse$center)[1]
+  y <- unlist(ellipse$center)[2]
+  points <- getEllipsePoints(ellipse$angle, 0.5 * ellipse$major, 0.5 * ellipse$minor, x, y)
+  points <- t(data.frame(points))
+
+  result <- fitEllipsePoints(points)
+
+  cov <- result$covar
+  colnames(cov) <- c(gate$xChannel, gate$yChannel)
+  rownames(cov) <- c(gate$xChannel, gate$yChannel)
+
+  return(flowCore::ellipsoidGate(filterId = gate$name, .gate = cov, mean = c(result$x, result$y)))
+}
+
+toFlowCorePolygonGate <- function(gate) {
+  m <- gate$model$polygon$vertices
+  if (is.null(dim(m))) {
+    m <- m[[1]]
+  }
+  colnames(m) <- c(gate$xChannel, gate$yChannel)
+  flowCore::polygonGate(filterId = gate$name, m)
+}
+
+#' Convert ScaleSet
+#'
+#' Convert a CellEngine ScaleSet to a flowCore transformList
+#'
+#' @param scaleSet The CellEngine scaleSet to be converted
+scaleSetToTransformList <- function(scaleSet) {
+  if (!requireNamespace("flowCore")) {
+    stop("This function requires the 'flowCore' package to be installed.")
+  }
+  scales <- scaleSet$scales[[1]]
+
+  funs <- sapply(seq_len(nrow(scales)), function(i) {
+    x <- scales$scale[i, ]
+    switch(x$type,
+      "LinearScale" = function(a) a,
+      "LogScale" = function(a) log10(pmax(1, a)),
+      "ArcSinhScale" = function(a) asinh(a / x$cofactor)
+    )
+  })
+
+  flowCore::transformList(
+    from = scales$channelName,
+    tfun = funs,
+    transformationId = scaleSet$name
+  )
+}

R/covarToParameters.R

@@ -0,0 +1,26 @@
+#' Covert covariance matrix to ellipse parameters
+#'
+#' Returns list("major" = major, "minor" = minor, "angle" = angle)
+#'
+#' @param covar Covariance matrix
+covarToParameters <- function(covar) {
+  eigens <- eigen(covar, symmetric = TRUE)
+  eigenvalues <- eigens$values
+
+  a <- covar[1]
+  b <- covar[2]
+  c <- covar[4]
+
+  major <- sqrt(eigenvalues[1])
+  minor <- sqrt(eigenvalues[2])
+
+  if ((b == 0) & a >= c) {
+    alpha <- 0
+  } else if ((b == 0) & a < c) {
+    alpha <- pi / 2
+  } else {
+    alpha <- atan2(eigenvalues[1] - a, b)
+  }
+
+  return(list("major" = major, "minor" = minor, "angle" = alpha))
+}

R/fitEllipsePoints.R

@@ -0,0 +1,86 @@
+#' Fit ellipse points
+#'
+#' Calculate the minimum volume enclosing ellipsoid for a set of points
+#'
+#' Uses the Khachiyan Algorithm to find the the minimum volume enclosing
+#' ellipsoid (MVEE) of a set of points. In two dimensions, this is just
+#' the bounding ellipse (this function is limited to two dimensions).
+#' Adapted from https://stat.ethz.ch/pipermail/r-help/2011-March/272997.html
+#'
+#' @param xy a two-column data frame containing x and y coordinates.
+#' @param tolerance a tolerance value (default = 0.005).
+#' @param max.iter Maximum algorithm iterations.
+#'
+#' @return A named list containing:
+#'  - "covar", the 2x2 covariance matrix of the ellipse
+#'  - "x", the x-coordinate of the ellipse center
+#'  - "y", the y-coordinate of the ellipse center
+#'  - "major", the length of the major axis
+#'  - "minor", the length of the minor axis
+#'  - "angle", the angle of rotation
+fitEllipsePoints <- function(xy, tolerance = 1e-20, max.iter = 1000) {
+  if (ncol(xy) != 2) {
+    stop("xy must be a two-column data frame")
+  }
+  n <- nrow(xy)
+  d <- ncol(xy)
+  if (n <= d) stop("There cannot be more columns than rows")
+
+  ## Add a column of 1s to the (n x 2) matrix xy - so it is now (n x 3)
+  Q <- t(cbind(xy, rep(1, n))) # nolint
+
+  ## Initialize
+  count <- 1
+  err <- 1
+  u <- rep(1 / n, n)
+  ## Khachiyan Algorithm
+  while (err > tolerance) {
+    X <- Q %*% diag(u) %*% t(Q) # nolint
+    M <- diag(t(Q) %*% solve(X) %*% Q) # nolint
+    maximum <- max(M)
+    j <- which(M == maximum)
+    step_size <- (maximum - d - 1) / ((d + 1) * (maximum - 1))
+    new_u <- (1 - step_size) * u
+    new_u[j] <- new_u[j] + step_size
+    err <- sqrt(sum((new_u - u)^2))
+    count <- count + 1
+    if (count > max.iter) {
+      message <- paste(
+        "Iterated ", max.iter, " times and
+       still can't find the bounding
+       ellipse. \n Either increase the
+       tolerance or the maximum number of
+       iterations. \n",
+        sep = ""
+      )
+      stop(message)
+    }
+    u <- new_u
+  }
+  # Put the elements of the vector u into the diagonal of a matrix
+  U <- diag(u) # nolint
+  # Take the transpose of xy
+  P <- t(xy) # nolint
+  # Compute the center, adding back the shifted values
+  c <- as.vector((P %*% u))
+  x <- c[1]
+  y <- c[2]
+  # Compute the A-matrix
+  A <- (1 / d) * solve(P %*% U %*% t(P) - (P %*% u) %*% t(P %*% u)) # nolint
+  # covar is the covariance matrix of the ellipse
+  covar <- solve(A)
+  # Find the Eigenvalues of matrix A which will be used to get the major and minor axes
+  A.eigen <- eigen(A) # nolint
+  # Calculate the length of the semi-major and semi-minor axes
+  # from the Eigenvalues of A.
+  semi.axes <- sort(1 / sqrt(A.eigen$values), decreasing = TRUE) # nolint
+  major <- semi.axes[1]
+  minor <- semi.axes[2]
+  # Calculate the rotation angle from the first Eigenvector
+  alpha <- atan2(A.eigen$vectors[2, 1], A.eigen$vectors[1, 1]) - pi / 2
+
+
+  params <- covarToParameters(covar)
+
+  list("covar" = covar, "x" = x, "y" = y, "major" = params$major, "minor" = params$minor, angle = params$angle)
+}

R/getEllipsePoints.R

@@ -0,0 +1,24 @@
+#' Get ellipse points
+#'
+#' Returns 8 perimeter points at angular increments of pi / 4.
+#'
+#' @param x X-coordinate of the ellipse center.
+#' @param y Y-coordinate of the ellipse center.
+#' @param angle Angle of the ellipse.
+#' @param major Length of the major axis.
+#' @param minor Length of the minor axis.
+getEllipsePoints <- function(angle, major, minor, x, y) {
+  points <- list()
+  phi <- 0
+
+  while (phi < (pi * 2)) {
+    p <- c(
+      x + major * cos(phi) * cos(angle) - minor * sin(phi) * sin(angle),
+      y + major * cos(phi) * sin(angle) + minor * sin(phi) * cos(angle)
+    )
+
+    points[length(points) + 1] <- list(p)
+    phi <- phi + (pi / 4)
+  }
+  points
+}

R/getEvents.R

@@ -87,7 +87,8 @@ getEvents <- function(experimentId,
 
   fullURL <- paste(
     paste(pkg.env$baseURL, "experiments", experimentId, "fcsfiles", fcsFileId, sep = "/"),
-    format, sep = "."
+    format,
+    sep = "."
   )
 
   params <- list(