R/gating-methods.R

In openCyto: Hierarchical Gating Pipeline for flow cytometry data

#' @templateVar old gating
#' @templateVar new gt_gating
#' @template template-depr_pkg
NULL

#' @importFrom flowWorkspace gh_pop_is_negated
#' @export
gt_gating <- function(x, y, ...){
  UseMethod("gt_gating")
}
#' @export
gating <- function(x, y, ...){
  .Deprecated("gt_gating")
  gt_gating(x,y,...)
}

#' Applies a gatingTemplate to a GatingSet.
#'
#' It loads the gating methods by topological order and applies them to \code{GatingSet}.
#' 
#' @name gt_gating
#' @usage gt_gating(x, y, ...)
#' @aliases 
#' gating
#' gating,gatingTemplate,GatingSet-method
#' gt_gating,gatingTemplate,GatingSet-method
#' @param x a \code{gatingTemplate} object
#' @param y a \code{GatingSet} object
#' @param ... 
#'  \itemize{
#'      \item{start}{ a \code{character} that specifies the population (correspoding to 'alias' column in csv template) where the gating process will start from. It is useful to quickly skip some gates and go directly to the target population in the testing run. Default is "root".}
#'      \item{stop.at}{ a \code{character} that specifies the population (correspoding to 'alias' column in csv template) where the gating prcoess will stop at. Default is NULL, indicating the end of gating tree.}
#'      \item{keep.helperGates}{a \code{logical} flag indicating whether to keep the intermediate helper gates that are automatically generated by openCyto. Default is TRUE.}
#'      \item{mc.cores}{ passed to \code{multicore} package for parallel computing}
#'      \item{parallel_type}{  \code{character} specifying the parallel type. The valid options are "none", "multicore", "cluster".}
#'      \item{cl}{ \code{cluster} object passed to \code{parallel} package (when \code{parallel_type} is "cluster")}
#'  }  
#' 
#' 
#' @return 
#' Nothing. As the side effect, gates generated by gating methods are saved in \code{GatingSet}. 
#' @examples
#' \dontrun{
#'  gt <- gatingTemplate(file.path(path, "data/ICStemplate.csv"), "ICS")
#'  gs <- GatingSet(fs) #fs is a flowSet/ncdfFlowSet
#'  gt_gating(gt, gs)
#'  gt_gating(gt, gs, stop.at = "v") #proceed the gating until population 'v'
#'  gt_gating(gt, gs, start = "v") # start from 'v'
#'  gt_gating(gt, gs, parallel_type = "multicore", mc.cores = 8) #parallel gating using multicore
#'  #parallel gating by using cluster
#'  cl1 <- makeCluster (8, type = "MPI")
#'  gt_gating(gt, gs, parallel_type = "cluster", cl = cl1)
#'  stopCluster ( cl1 )
#' }
#' @export
gt_gating.gatingTemplate <- function(x, y, ...) {
  .gating_gatingTemplate(x, y, ...)
}
# gt_gating.gatingTemplate <- function(x, y, env_fct = NULL, ...) {
#       .gating_gatingTemplate(x,y,env_fct,...)
#     }

#' internal function (gating_gatingTemplate)
#'  
#' @param stop.at a \code{character} that specifies the population (correspoding to 'alias' column in csv template) where the gating prcoess will stop at.
#' @param start a \code{character} that specifies the population (correspoding to 'alias' column in csv template) where the gating prcoess will start from. It is useful to quickly skip some gates and go directly to the target population in the testing run. 
#' @param env_fct a \code{environment} that contains \code{fcTree} object named as 'fct'. If NULL (by default), no \code{fcTree} will be constructed. It is currently reserved for the internal debugging.
#' @param ... other arguments passed to the gatingMethod-specific \code{gating} methods.
#' @importFrom RBGL tsort
#' @importFrom plyr ldply
#' @keywords internal
#' @noRd 
.gating_gatingTemplate <- function(x, y, env_fct = NULL, start = "root", stop.at = NULL, keep.helperGates = TRUE, ...) {
  gt <- x
  if (!is.null(env_fct)) {
    # use the fcTree if already exists
    if (exists("fct", env_fct)) {
      fct <- get("fct", env_fct)
    } else {
      # create one from gt if not
      fct <- fcTree(gt)
      assign("fct", fct, env_fct)
    }
    
  }
  
  nodePaths <- names(sapply(gt_get_nodes(gt),alias))
    #validity check for stop.at argument
  if(!is.null(stop.at)){
    #escape meta character
    stop.pat <- gsub("\\+", "\\\\+", stop.at)
    stop.pat <- paste0(stop.pat, "$")
    if(substr(stop.at, 1,1) != "/")
      stop.pat <- paste0("/", stop.pat) #prepend path delimiter to avoid partial node name matching
    
    matchInd <- grep(stop.pat, nodePaths)
    if(length(matchInd) == 0)
      stop("Can't find stop point: ", stop.at)
    else if(length(matchInd) > 1)
      stop("ambiguous stop point: ", stop.at)
  }
  # gate each node 
  gt_nodes <- tsort(gt)[-1]#by the topological order

  #try to skip some nodes to save time
  if(start != "root"){
    #escape meta character
    start.pat <- gsub("\\+", "\\\\+", start)
    start.pat <- paste0(start.pat, "$") #treat it as terminal node
    if(substr(start, 1,1) != "/")
      start.pat <- paste0("/", start.pat) #prepend path delimiter to avoid partial node name matching
    matchInd <- grep(start.pat, nodePaths)
    
    if(length(matchInd) == 0)
      stop("Can't find start point: ", start)
    else if(length(matchInd) > 1)
      stop("ambiguous start point: ", start)
    
    gt_nodes <- gt_nodes[matchInd:length(gt_nodes)]
  }  
    
                   
  for (node in gt_nodes) {
    
    # get parent node to gate
    gt_node_pop <- gt_get_nodes(gt, node)
    parent <- gt_get_parent(gt, node)
    
    if(!is.null(stop.at)){
      nodeName <- alias(gt_node_pop)
      nodePath <- file.path(parent, nodeName)
      matchInd <- grep(stop.pat, nodePath)
      if(length(matchInd) == 1)
      {
        message("stop at: ",stop.at)
        break
      }
      
    }    
    # extract gate method from one edge(since multiple edge to the same node is
    # redudant)
    this_gate <- gt_get_gate(gt, parent, node)
    
    #get preprocessing method
    this_ppm <- ppMethod(gt, parent, node)
    
    parentInd <- match(parent, gs_get_pop_paths(y[[1]], showHidden = TRUE))
    if (is.na(parentInd)) 
      stop("parent node '", parent, "' not gated yet!")
    
    #preprocessing
    pp_res <- NULL
#    browser()
    if(class(this_ppm) == "ppMethod")
      pp_res <- preprocessing(x = this_ppm, y, parent = parent, gtPop = gt_node_pop, gm = this_gate, ...)
#    browser()
    # pass the pops and gate to gating routine
    filterObj <- gt_gating(x = this_gate, y, parent = parent, gtPop = gt_node_pop, pp_res = pp_res, ...)
    if(!keep.helperGates)
    {
      gt_delete_helpergates(x, y)
    }
    # update fct
    if (!is.null(env_fct) && !is.null(filterObj)) {
      nodeData(env_fct$fct, node, "fList")[[1]] <- filterObj
    }  
    
    
  }
  message("finished.")
}


#' apply a \link[openCyto:gtMethod-class]{gtMethod} to the \code{GatingSet}
#' 
#' The actual workhorse of most gating methods
#' 
#' @param x \code{gtMethod}
#' @param y \code{GatingSet}
#' @param ... other arguments 
#' @aliases
#' gating,gtMethod,GatingSet-method 
#' gating,gtMethod,GatingSetList-method
#' @keywords internal
#' @noRd
gt_gating.gtMethod <- function(x, y, ...) {
  .gating_gtMethod(x,y,...)
}

#' evalRd tag is currently not working because the output of eval string isn't positioned properly in Rd file
#' @noRd
roxygen_parameter <- function() {
  paste(" 
          gtPop a \\code{gtPopulation} object that contains the information of the cell population that is going to be generated by this gating method.
          parent\\code{character} specifying the parent node within\\code{GatingSet} to which the new popoulation is to be attached.
          pp_res\\code{list} the preprocessing results to be used by the gating method. The names of the list has to be matched to the sample names of the data.
          mc.cores passed to\\code{multicore} package for parallel computing
          parallel_type\\code{character} specifying the parallel type. The valid options are 'none', 'multicore', 'cluster'.
          cl \\code{cluster} object passed to \\code{parallel} package")
}
#' internal function (gating_gtMethod)
#' 
#' It is a generic gating function that does:
#' 1. parse the gating parameters
#' 2. group the data when applicable 
#' 3. apply parallelism when applicable  
#' 4. pass the flow data(maybe grouped) and preprocessing results to the adaptor function ".gating_adaptor"
#' 5. collect the gates and add to GatingSet object
#' 
#' 
#' @param x \code{gtMethod}
#' @param y \code{GatingSet}
#' @evalRd roxygen_parameter()
#' @import flowWorkspace
#' @importFrom flowCore getChannelMarker
#' @noRd
.gating_gtMethod <- function(x, y, gtPop, parent, pp_res 
            , mc.cores = getOption("mc.cores", 2L), parallel_type = c("none", "multicore", "cluster"), cl = NULL
            ,  ...) {
  
  requireNamespace("parallel")
  gFunc_args <- parameters(x)

  # HOTFIX: This resolve an error when args is a named list with name NA and object NA.
  # The resulting error occurs down below and is:
  # Error in thisCall[[arg]] <- args[[arg]] : subscript out of bounds
  if (!is.null(names(gFunc_args))) {
    gFunc_args <- gFunc_args[!is.na(names(gFunc_args))]
  }
  
  ## Optionally omit arguments that should not be used in this context
  if(x@name == "gate_tail"){
    to_omit <- na.omit(match("positive", names(gFunc_args)))
    if(length(to_omit) > 0)
      gFunc_args <- gFunc_args[-to_omit]
  }
  ##  
    
  gm <- paste0(".", names(x))
  
  dims <- dims(x)
  
  is_1d_gate <- length(dims) == 1
  
  popAlias <- alias(gtPop)
  popName <- names(gtPop)
  popId <- gtPop@id
  gs_nodes <- basename(gs_pop_get_children(y[[1]], parent))
  
  if (length(gs_nodes) == 0)
    isGated <- FALSE
  else
    isGated <- any(popAlias %in% gs_nodes)
    
  if(!isGated)
  {
    message("Gating for '", popAlias, "'")
    
    parent_data <- gs_pop_get_data(y, parent)
    parallel_type <- match.arg(parallel_type)
    ## get the accurate channel name by matching to the fr
    frm <- parent_data[[1, use.exprs = FALSE]]
    channels <-  unname(sapply(dims, function(channel)as.character(getChannelMarker(frm, channel)$name)))
    
    parent_data <- parent_data[, channels] #it is more efficient to only pass the channels of interest
    # Splits the flow set into a list.
    # By default, each element in the list is a flowSet containg one flow frame,
    # corresponding to the invidual sample names.
    # If 'split' is given, we split using the unique combinations within pData.
    # In this case 'split' is specified with column names of the pData.
    # For example, "PTID:VISITNO"
    # when split is numeric, do the grouping by every N samples
    groupBy <- groupBy(x)
    isCollapse <- isCollapse(x)
    if (groupBy != "" && isCollapse) {
      #when x@collapse == FALSE,then ignore groupBy argument since grouping is only used for collapsed gating
      split_by <- as.character(groupBy)
      split_by_num <- as.numeric(split_by)
      #split by every N samples
      if(!is.na(split_by_num)){
            nSamples <- length(parent_data)
            if(nSamples==1){
              split_by <- 1
            }else{
              split_by <-  sample(rep(1:nSamples, each = split_by_num, length.out= nSamples))  
            }
            
      }else{
        #split by study variables
        split_by <- strsplit(split_by, ":")[[1]]
        split_by <- apply(pData(parent_data)[, split_by, drop = FALSE], 1, paste, collapse = ":")
        split_by <- as.character(split_by)
      }
    } else {
      split_by <- sampleNames(parent_data)
    }

    fslist <- split(parent_data, split_by)
    
    if(is.null(pp_res))
      pp_res <- sapply(names(fslist),function(i)pp_res)
    else
      pp_res <- pp_res[names(fslist)] #reorder pp_res to make sure it is consistent with fslist
    # construct method call
    thisCall <- substitute(f1(fslist,pp_res))
    thisCall[["FUN"]] <- as.symbol(".gating_adaptor")
    # args to be passed to gating_adaptor
    args <- list()
    args[["gFunc"]] <- gm  #set gating method
    args[["popAlias"]] <- popAlias  
    args[["channels"]] <- channels 
    
    negated <- popName=="-"
    args[["gFunc_args"]] <- gFunc_args
    
    thisCall[["MoreArgs"]] <- args
    
    ## choose serial or parallel mode
    
      
    if (parallel_type == "multicore") {
      message("Running in parallel mode with ", mc.cores, " cores.")
      thisCall[[1]] <- quote(mcmapply)
      thisCall[["mc.cores"]] <- mc.cores
    }else if(parallel_type == "cluster"){
      if(is.null(cl))
          stop("cluster object 'cl' is empty!")
        thisCall[[1]] <- quote(clusterMap)
        thisCall[["cl"]] <- cl
        thisCall[["fun"]] <- thisCall[["FUN"]] 
        thisCall[["FUN"]] <- NULL
        thisCall[["SIMPLIFY"]] <- TRUE
     }else {
      thisCall[[1]] <- quote(mapply)  #select loop mode
  
     }
    
    flist <- eval(thisCall)

    # Handles the case that 'flist' is a list of lists.
    #   The outer lists correspond to the split by pData factors.
    #   The inner lists contain the actual gates.
    #     The order do not necessarily match up with sampleNames()
    # Unforunately, we cannot simply use 'unlist' because the list element names
    # are mangled.
    
    if (all(sapply(flist, is.list))) {
      flist <- do.call(c, unname(flist))
    }
    #check failed sample 
    #eventually we want to handle this properly (like inserting dummy gates)
    #in order for the other samples proceed the gating
    failed <- sapply(flist, function(i)extends(class(i), "character"))
    if(any(failed)){
      print(flist[failed])
      stop("some samples failed!")
      
    }
    
    #this is flowClust-specific operation, which
    # be abstracted out of this framework
    if (extends(class(flist[[1]]), "fcFilter")) {
      fcflist <- try(fcFilterList(flist), silent = TRUE)
      if(!is(fcflist, "try-error"))
        flist <- fcflist
      }
    
    
    
    
    if(length(popAlias) == 1){
      #when Alias is meta character, then pass NULL
      # to add method which uses filterId slot to name the populations
      if(popAlias == "*")
        popAlias <- NULL  
    }

    # For gate_quad methods, need to filter down to just the gates that were asked for
    if(names(x) %in% c("quadGate.seq", "gate_quad_sequential", "quadGate.tmix", "gate_quad_tmix")){
      pops <- gtPop@name
      pops <- gsub("([\\+-])([^/$])", "\\1&\\2", pops)
      pops <- strsplit(pops, "&")[[1]]
      pops <- strsplit(pops, "/")
      pops <- paste0(rep(pops[[1]], each=length(pops[[2]])), pops[[2]])
      pops <- match(pops, c("-+", "++", "+-", "--"))
      flist <- lapply(flist, function(sublist) filters(sublist[pops]))
    }
    
    gs_node_id <- gs_pop_add(y, flist, parent = parent, name = popAlias, validityCheck = FALSE, negated = negated)
    if(!is.null(popAlias))
      recompute(y, file.path(parent, popAlias))
    else
      recompute(y, parent)
	message("done.")
    
  }else{
    
    message("Skip gating! Population '", paste(popAlias, collapse = ","), "' already exists.")
    flist <- NULL
  }
  
  flist
}

# apply a \code{boolMethod} to the \code{GatingSet}
#' @rdname gt_gating
#' @aliases
#' gating,boolMethod,GatingSet-method
#' gating,boolMethod,GatingSetList-method
#' @noRd
gt_gating.boolMethod <- function(x, y, ...) {
  .gating_boolMethod(x,y,...)
}

#' internal function (gating_boolMethod)
#' 
#'  @param y \code{GatingSet}
#'  @param gtPop a \code{gtPopulation} object that contains the information of the cell population that is going to be generated by this gating method.
#'  @param  parent \code{character} specifying the parent node within \code{GatingSet} to which the new popoulation is to be attached.
#' @noRd 
.gating_boolMethod <- function(x, y, gtPop, parent, ...) {
  
  args <- parameters(x)[[1]]
  gm <- paste0(".", names(x))
  popAlias <- alias(gtPop)
  popName <- names(gtPop)
  popId <- gtPop@id
  
  gs_nodes <- basename(gs_pop_get_children(y[[1]], parent))
  tNodes <- deparse(args)
  if (!(popAlias %in% gs_nodes)) {
    message(popAlias, " gating...")
    bf <- eval(substitute(booleanFilter(x), list(x = args)))
    bf@filterId <- tNodes
    #set recompute to FALSE because we want recompute method take over the
    #computing job since it is smart on determining whehter flow data needs to be loaded
    #for boolean gates
    invisible(gs_node_id <- gs_pop_add(y, bf, parent = parent, name = popAlias))
    newNode <- file.path(parent, popAlias)
    invisible(recompute(y, newNode))
    message("done.")
  } else {
    message("Skip gating! Population '", popAlias, "' already exists.")
  }
  
  # gs_node_id
  NULL
}

# apply a \code{polyFunctions} gating method to the \code{GatingSet}
# 
# It generates a batch of \code{boolMethod}s based on the expression defined in \code{polyFunctions} objects.
# It is a convenience way to generate different boolean combinations of cytokine gates. 
#' @rdname gt_gating
#' @aliases
#' gating,polyFunctions,GatingSet-method
#' gating,polyFunctions,GatingSetList-method
#' @noRd
gt_gating.polyFunctions <- function(x, y, ...) {
  .gating_polyFunctions(x,y,...)
}

#' internal function (gating_polyFunctions)
#' 
#'  @param y \code{GatingSet}
#'  @param gtPop a \code{gtPopulation} object that contains the information of the cell population that is going to be generated by this gating method.
#'  @param  parent \code{character} specifying the parent node within \code{GatingSet} to which the new popoulation is to be attached.
#' @importFrom gtools permutations 
#' @noRd 
.gating_polyFunctions <- function(x, y, gtPop, parent, ...) {
  
  refNodes <- x@refNodes
  popAlias <- alias(gtPop)
  message("Population '", paste(popAlias, collapse = ","), "'")
  
  nMarkers <- length(refNodes)

  ## all the comibnations of A & B & C
  opList <- permutations(n = 1, r = nMarkers - 1, c("&"), repeats.allowed = TRUE)
  isNotList <- permutations(n = 2, r = nMarkers, c("!", ""), repeats.allowed = TRUE)
  polyExprsList <- apply(opList, 1, function(curOps) {
    apply(isNotList, 1, function(curIsNot) {
      polyExprs <- curIsNot
      polyExprs[-1] <- paste0(curOps, curIsNot[-1])
      
      paste(paste0(polyExprs, refNodes), collapse = "")
    })
  })
  polyExprsList <- as.vector(polyExprsList)
  
  gs_nodes <- basename(gs_pop_get_children(y[[1]], parent))
  # actual gating
  lapply(polyExprsList, function(polyExpr) {
    
    #replace the slash with colon 
    #since forward slash is reserved for gating path
    if(grepl("/",polyExpr)){
      old_name <- polyExpr
      new_name <- gsub("/",":",polyExpr)
      warning(old_name, " is replaced with ", new_name)
    }else
      new_name <- polyExpr
#    browser()
    
    isExist <- new_name %in% gs_nodes
    if (!isExist) {
      message("adding ", new_name, " ...")
      
      bf <- eval(substitute(booleanFilter(v), list(v = as.symbol(polyExpr))))

      invisible(gs_node_id <- .addGate_fast(y, bf, parent = parent, name = polyExpr))
      
  
    } else {
      message("Skip!Population '", new_name, "' already exists.")
    }        
  })
  #to reduce overhead,compute from parent node once instead of do it multiple times for each individual new bool gate
    
  invisible(recompute(y, parent))
  message("done.")
  
  list()
}

# apply a \code{refGate} to the \code{GatingSet}
#' @rdname gt_gating
#' 
#' @aliases
#' gating,refGate,GatingSet-method
#' gating,refGate,GatingSetList-method
#' gating,dummyMethod,GatingSet-method
#' gating,dummyMethod,GatingSetList-method
#' @noRd
gt_gating.refGate <- function(x, y, ...) {
  .gating_refGate(x, y, ...)
}

gt_gating.dummyMethod <- function(x, y, ...) {
  #do nothing
}

#' internal function (gating_refGate)
#' 
#'  @param y \code{GatingSet}
#'  @param gtPop a \code{gtPopulation} object that contains the information of the cell population that is going to be generated by this gating method.
#'  @param  parent \code{character} specifying the parent node within \code{GatingSet} to which the new popoulation is to be attached.
#' @noRd 
.gating_refGate <- function(x, y, gtPop, parent, ...) {
#  negated <- FALSE
  refNodes <- x@refNodes
  popAlias <- alias(gtPop)
  popName <- names(gtPop)
  dims <- dims(x)
  
  my_gh <- y[[1]] 
  gs_nodes <- basename(gs_pop_get_children(my_gh, parent))
  if (length(gs_nodes) == 0 || !popAlias %in% gs_nodes) {
    
    message("Population '", paste(popAlias, collapse = ","), "'")
    if (length(refNodes) > 2) {
      stop("Not sure how to construct gate from more than 2 reference nodes!")
    }
    
    fr <- gh_pop_get_data(my_gh, use.exprs = FALSE)
    
    #check if parent node is shared 
    #to determine whether simply grab the indices from ref nodes without recompute
    ref_parents <- sapply(refNodes, function(refNode)gs_pop_get_parent(my_gh, refNode))
    isSameParent <- all(ref_parents == parent)
    flist <- flowWorkspace::lapply(y, function(gh) {
          
       glist <- lapply(refNodes, function(refNode) {
          
          gh_pop_get_gate(gh, refNode)
        })
       
      # standardize the names for the gate parameters and dims
      gate_params <- unlist(lapply(glist, function(g) {
        cur_param <- parameters(g)
        getChannelMarker(fr, cur_param)["name"]
      }))

      if(length(glist)==1){
        #1d ref gate
        dims <- dims[!is.na(dims)]
        nDims <- length(dims)
        
        
        if(nDims==2){
          # #before flowCore and flowViz support the negated filter
          # #we use refGate+boolGate in csv template as the workaround
          # if (popName == "-") {
          #   stop("negated 2d gate is not supported yet!")
          # }
          #pass the gate as it is 
          glist[[1]]
        }else{
          dim_params <-  getChannelMarker(fr, dims)["name"]
          y_g <- glist[[1]]               
          y_coord <- c(y_g@min, y_g@max)
          cut.y <- y_coord[!is.infinite(y_coord)]
          
          
          if (popName == "+") {
            #handle all infinite coordinates
            if(length(cut.y) == 0)
              cut.y <- -Inf
            gate_coordinates <- list(c(cut.y, Inf))
          } else if(popName=="-"){
            if(length(cut.y) == 0)
              cut.y <- Inf
            gate_coordinates <- list(c(-Inf, cut.y))
          }else{
            stop("unknown population pattern, ",popName)
          }
          names(gate_coordinates) <- as.character(dim_params)
          
          rectangleGate(gate_coordinates)  
        }
        
      }else{
          #2d quad gate
          dim_params <- unlist(lapply(dims, function(dim) {
                    getChannelMarker(fr, dim)["name"]
                  }))
          
          # match the gate param to dims to find gates for x, y dimensions
          x_ref_id <- match(dim_params[1], gate_params)
          y_ref_id <- match(dim_params[2], gate_params)
          x_g <- glist[[x_ref_id]]
          y_g <- glist[[y_ref_id]]
                        
          
          # pick the non-infinite coordinate as the cut points
          x_coord <- c(x_g@min, x_g@max)
          y_coord <- c(y_g@min, y_g@max)
          x_inf_vec <- is.infinite(x_coord)
          y_inf_vec <- is.infinite(y_coord)
          cut.x <- x_coord[!x_inf_vec]
          cut.y <- y_coord[!y_inf_vec]
          
          # In order, the following vector has the patterns:
          # 1. top left (-+)
          # 2. top right (++)
          # 3. bottom right (+-)
          # 4. bottom left (--)
          quadPatterns <- c("-+", "++", "+-", "--")
          
          quadInd <- match(popName, quadPatterns)
          
          is_negated <- lapply(refNodes, function(refNode) {
            
            gh_pop_is_negated(gh, refNode)
          })          
          x_ref <- refNodes[x_ref_id] 
          y_ref <- refNodes[y_ref_id]
          x_negated <- is_negated[[x_ref_id]]
          y_negated <- is_negated[[y_ref_id]]
          #take into account of (-Inf, x) for bool operations in ocRectRefGate
          if(x_coord[x_inf_vec] < 0)
            x_negated <- !x_negated
          if(y_coord[y_inf_vec] < 0)
            y_negated <- !y_negated
          #standardize event ind to x+ and y+
          #by checking the positive sign of both reference gates
          #to prepare for the bool operation later
          if(length(cut.x) == 0){
          #handle dummy ref gate that has both boudnary as Inf
#            x_event_ind <- TRUE           
            x_ref <- NULL
          }
          if(length(cut.y) == 0){
#            y_event_ind <- TRUE
              y_ref <- NULL
          }
          
#          browser()
          # construct rectangleGate from reference cuts
          if (quadInd == 1) {#-+
            #handle all infinite coordinates
            if(length(cut.x) == 0)
              cut.x <- Inf
            else
            {
              if(!x_negated)
                x_ref <- paste0("!", x_ref) 
            }
              
            if(length(cut.y) == 0)
              cut.y <- -Inf
            else
            {
              if(y_negated)
                y_ref <- paste0("!", y_ref)
            }
            coord <- list(c(-Inf, cut.x), c(cut.y, Inf))
            
          } else if (quadInd == 2) {#++
            #handle all infinite coordinates
            if(length(cut.x) == 0)
              cut.x <- -Inf
            else
            {
              if(x_negated)
                x_ref <- paste0("!", x_ref) 
            }
            if(length(cut.y) == 0)
              cut.y <- -Inf
            else
            {
              if(y_negated)
                y_ref <- paste0("!", y_ref) 
            }
            coord <- list(c(cut.x, Inf), c(cut.y, Inf))
            
          } else if (quadInd == 3) {#+-
            #handle all infinite coordinates
            if(length(cut.x) == 0)
              cut.x <- -Inf
            else
            {
              if(x_negated)
                x_ref <- paste0("!", x_ref) 
            }
            if(length(cut.y) == 0)
              cut.y <- Inf
            else
            {
              if(!y_negated)
                y_ref <- paste0("!", y_ref) 
            }
            coord <- list(c(cut.x, Inf), c(-Inf, cut.y))
             
          } else if (quadInd == 4) {#--
            #handle all infinite coordinates
            if(length(cut.x) == 0)
              cut.x <- Inf
            else
            {
              if(!x_negated)
                x_ref <- paste0("!", x_ref) 
            }
            if(length(cut.y) == 0)
              cut.y <- Inf
            else
            {
              if(!y_negated)
                y_ref <- paste0("!", y_ref) 
            }
            coord <- list(c(-Inf, cut.x), c(-Inf, cut.y))
            
          } else stop("Pop names does not match to any quadrant pattern!")
          
          names(coord) <- as.character(dim_params)
          fres <- rectangleGate(coord)
          if(isSameParent&&!is.null(x_ref)&&!is.null(y_ref)){
            fres <- ocRectRefGate(fres, paste0(x_ref, "&", y_ref))
            fres
          }
          
        fres  
      }
      
    })
    
    flist <- filterList(flist)
    gs_node_id <- gs_pop_add(y, flist, parent = parent, name = popAlias, validityCheck = FALSE)
	if(!is(flist[[1]], "ocRectRefGate"))
		recompute(y, file.path(parent, popAlias))
	
   
  } else {
    message("Skip gating! Population '", popAlias, "' already exists.")
    flist <- NULL
  }

  message("done.")
  flist
}