Quicker subsetting of points using categorized 'im' than polygons and…

… 'over'

In gating() function, creates a categorized 'im' based on critical p-value, assigns that value to every point in a 'ppp' object, and subsets points by category.

Removed 'maptools' and 'sp' packages from Imports

Updated links in package.R for change in packages

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: gateR
 Type: Package
 Title: Flow/Mass Cytometry Gating via Spatial Kernel Density Estimation
-Version: 0.1.7-1
-Date: 2021-02-02
+Version: 0.1.7
+Date: 2021-02-04
 Authors@R: 
     c(person(given = "Ian D.",
              family = "Buller",
@@ -48,9 +48,7 @@ Imports:
     graphics,
     grDevices,
     lifecycle,
-    maptools,
     raster,
-    sp,
     sparr,
     SpatialPack,
     spatstat.geom,

NAMESPACE

@@ -20,17 +20,15 @@ importFrom(lifecycle,badge)
 importFrom(lifecycle,deprecate_warn)
 importFrom(lifecycle,deprecated)
 importFrom(lifecycle,is_present)
-importFrom(maptools,unionSpatialPolygons)
 importFrom(raster,cut)
 importFrom(raster,extent)
 importFrom(raster,raster)
-importFrom(raster,rasterToPolygons)
 importFrom(raster,values)
-importFrom(sp,coordinates)
-importFrom(sp,over)
 importFrom(sparr,OS)
 importFrom(sparr,risk)
 importFrom(spatstat.geom,as.im)
+importFrom(spatstat.geom,cut.im)
+importFrom(spatstat.geom,marks)
 importFrom(spatstat.geom,owin)
 importFrom(spatstat.geom,ppp)
 importFrom(stats,na.omit)

R/gating.R

@@ -54,20 +54,17 @@
 #'
 #' @importFrom grDevices chull
 #' @importFrom lifecycle badge deprecate_warn deprecated is_present
-#' @importFrom maptools unionSpatialPolygons
-#' @importFrom raster rasterToPolygons values
-#' @importFrom sp coordinates over
-#' @importFrom spatstat.geom owin
+#' @importFrom spatstat.geom cut.im marks owin ppp
 #' @importFrom tibble add_column
 #' @export
 #' 
 #' @examples
 #' if (interactive()) {
 #' ## Single condition, no multiple testing correction
-#'   test_gate <- gateR::gating(dat = randCyto,
-#'                              vars = c("arcsinh_CD4", "arcsinh_CD38",
-#'                                       "arcsinh_CD8", "arcsinh_CD3"),
-#'                              n_condition = 1)
+#'   test_gate <- gating(dat = randCyto,
+#'                       vars = c("arcsinh_CD4", "arcsinh_CD38",
+#'                                "arcsinh_CD8", "arcsinh_CD3"),
+#'                       n_condition = 1)
 #' }
 #' 
 gating <- function(dat,
@@ -226,15 +223,15 @@ gating <- function(dat,
                c1n = c1n,
                ...)
     }
-
-    # Convert p-value surface into a categorized raster
+    
+    # Convert p-value surface into a categorized 'im'
     ## v == 2 : significant T1 numerator;  v == 1: not
-    Ps <- pval_plot(input = out$P, alpha = out$alpha)
+    Ps <- spatstat.geom::cut.im(out$P, breaks = c(-Inf, out$alpha / 2, 1 - out$alpha / 2, Inf), labels = c("1", "2", "3")) 
     list_gate[[k]] <- out # save for output
-    rm(out, df, win_gate) # conserve memory
+    rm(out, df) # conserve memory
 
     # Go back one gate if current gate has no significant area and produce output of previous gate
-    if (all(raster::values(Ps)[!is.na(raster::values(Ps))] == 2) | all(is.na(raster::values(Ps)))) {
+    if (all(Ps$v == "2", na.rm = TRUE) | all(is.na(Ps$v))) {
       if (k > 1) {
       message(paste("Gate", k, "yielded no significant", type_cluster, "cluster(s)...",
                 "Returning results from Gate", k-1,
@@ -243,19 +240,13 @@ gating <- function(dat,
       out_list <- list("obs" = output, "n" = n_out, "gate" = list_gate)
       return(out_list)
       } else {
-        stop(paste("Gate 1 yielded no significant", type_cluster, "cluster(s)... Returning no results", sep = " "))
+        stop(paste("Gate 1 yielded no significant clustering... Returning no results", sep = " "))
       }
     }
 
-    # convert categorized raster to gridded polygons
-    out_pol <- raster::rasterToPolygons(Ps)
-    rm(Ps) # conserve memory
+    if (numerator == TRUE) { v <- "1" } else { v <- "3" }
 
-    if (numerator == TRUE) { v <- 1 } else { v <- 3 }
-    # combine gridded polygons of similar categories
-    pols <- try(maptools::unionSpatialPolygons(out_pol[out_pol$layer == v, ],
-                                               IDs = rep(1, length(out_pol[out_pol$layer == v, ]))), silent = TRUE)
-    if("try-error" %in% class(pols)) {
+    if (!any(Ps$v == v, na.rm = TRUE)) {
       if (k > 1) {
       message(paste("Gate", k, "yielded no significant", type_cluster, "cluster(s)...",
                 "Returning results from Gate", k-1,
@@ -267,15 +258,16 @@ gating <- function(dat,
         stop(paste("Gate 1 yielded no significant", type_cluster, "cluster(s)... Returning no results", sep = " "))
       }
     }
-    rm(out_pol) # conserve memory
 
     # Overlay points
-    sp::coordinates(dat_gate) <- ~ cbind(dat_gate[,which(colnames(dat_gate) %in% vars[j])],
-                                         dat_gate[,which(colnames(dat_gate) %in% vars[j + 1])])
+    suppressMessages(suppressWarnings(full_ppp <- spatstat.geom::ppp(x = dat_gate[,which(colnames(dat_gate) %in% vars[j])],
+                                                                     y = dat_gate[,which(colnames(dat_gate) %in% vars[j + 1])],
+                                                                     marks = dat_gate,
+                                                                     window = win_gate)))
 
-    # extract points within significant cluster
-    dat_gate <- as.data.frame(dat_gate[!is.na(sp::over(dat_gate, pols)), ])
-    rm(pols) # conserve memory
+    # extract points within significant cluster(s)
+    spatstat.geom::marks(full_ppp)$gate <- Ps[full_ppp, drop = FALSE]
+    dat_gate <- spatstat.geom::marks(full_ppp)[spatstat.geom::marks(full_ppp)$gate == v, ]
 
     # Output for the final gate
     if (k == n_gate) {

R/package.R

@@ -24,7 +24,7 @@
 #' @aliases gateR-package gateR
 #' @docType package
 #' 
-#' @section Dependencies: The 'gateR' package relies heavily upon \code{\link{sparr}}, \code{\link{sp}}, and \code{\link{raster}}. For a two-level comparison, the spatial relative risk function uses the \code{\link[sparr]{risk}} function. The \code{gating} function extracts cells within significant clusters using the \code{\link[raster]{rasterToPolygons}}, \code{\link[maptools]{unionSpatialPolygons}}, and \code{\link[sp]{over}} functions. The calculation of a Bonferroni correction for multiple testing accounting for the spatial correlation of the estimated surface uses the \code{\link[SpatialPack]{modified.ttest}} function. Basic visualizations rely on the \code{\link[fields]{image.plot}} function.
+#' @section Dependencies: The 'gateR' package relies heavily upon \code{\link{sparr}}, \code{\link{spatstat.geom}}, and \code{\link{raster}}. For a two-level comparison, the spatial relative risk function uses the \code{\link[sparr]{risk}} function. The calculation of a Bonferroni correction for multiple testing accounting for the spatial correlation of the estimated surface uses the \code{\link[SpatialPack]{modified.ttest}} function. Basic visualizations rely on the \code{\link[fields]{image.plot}} function.
 #' 
 #' @author Ian D. Buller\cr \emph{Occupational and Environmental Epidemiology Branch, Division of Cancer Epidemiology and Genetics, National Cancer Institute, National Institutes of Health, Rockville, Maryland, USA.}
 #' 
@@ -37,12 +37,10 @@ NULL
 #' @importFrom graphics close.screen par screen split.screen 
 #' @importFrom grDevices chull colorRampPalette dev.off png
 #' @importFrom lifecycle badge deprecate_warn deprecated is_present
-#' @importFrom maptools unionSpatialPolygons
-#' @importFrom raster cut extent raster rasterToPolygons values
-#' @importFrom sp coordinates over
+#' @importFrom raster cut extent raster values
 #' @importFrom sparr OS risk
 #' @importFrom SpatialPack modified.ttest
-#' @importFrom spatstat.geom as.im owin ppp 
+#' @importFrom spatstat.geom as.im cut.im marks owin ppp
 #' @importFrom stats na.omit pnorm relevel
 #' @importFrom tibble add_column
 #' @importFrom fields image.plot

README.md

@@ -196,7 +196,6 @@ test_gating <- gateR::gating(dat = obs_dat,
                              vars = c("arcsinh_CD4", "arcsinh_CD38",
                                       "arcsinh_CD8", "arcsinh_CD3"),
                              n_condition = 1,
-                             p_correct = "FDR",
                              plot_gate = TRUE,
                              upper_lrr = 1,
                              lower_lrr = -1)