R/graphGML_methods.R

In RGLab/CytoML: A GatingML Interface for Cross Platform Cytometry Data Sharing

#' @include read.gatingML.cytobank.R
NULL

#' get nodes from {graphGML} object
#'
#' @param x \code{graphGML}
#' @param y \code{character} node index. When \code{missing}, return all the nodes
#' @param order \code{character} specifying the order of nodes. options are "default", "bfs", "dfs", "tsort"
#' @param only.names \code{logical} specifiying whether user wants to get the entire \code{nodeData} or just the name of the population node
#' @return It returns the node names and population names by default. Or return the entire nodeData associated with each node.
#' @importFrom flowWorkspace getNodes
#' @importFrom graph nodeData
#' @examples
#' \dontrun{
#  acsfile <- system.file("extdata/cytobank_experiment.acs", package = "CytoML")
#  ce <- open_cytobank_experiment(acsfile)
#  xmlfile <- ce$gatingML
#' g <- read.gatingML.cytobank(xmlfile)
#' getNodes(g)
#' getNodes(g, only.names = FALSE)
#' }
setMethod("getNodes", signature = c("graphGML"),
          definition = function(x, y
                                  , order = c("default", "bfs", "dfs", "tsort")
                                  , only.names = TRUE) {

  if (missing(y)){
    res <- nodeData(x)
    order <- match.arg(order)
    if(order != "default"){
      nodeIds <- eval(substitute(f1(x),list(f1=as.symbol(order))))
      if(order == "dfs")
        nodeIds <- nodeIds$discovered
      res <- res[nodeIds]
    }
  }else
  {
    res <- nodeData(x, y)

  }
  if(only.names){
    res <- sapply(res,`[[`,"popName")
  }
  if(length(res) == 1 && class(res) == "list")
      res <- res[[1]]
   res
})

#' get full path of the parent
#' @param x \code{graphGML}
#' @param y \code{character} node index. When \code{missing}, return all the nodes
#' @noRd
.getPath <- function(x, y){
  #get full path
  nodeIds <- y
  thisNodeID <- y
  while(length(thisNodeID) > 0){
    thisNodeID <- getParent(x, thisNodeID)
    nodeIds <- c(thisNodeID,nodeIds)
  }

  pops <- lapply(nodeIds, function(i)nodeData(x,i)[[1]][["popName"]])
  path <- paste(pops, collapse = "/")
  paste0("/", path)
}

#' get children nodes
#'
#' @param obj \code{graphGML}
#' @param y \code{character} parent node path
#' @return a graphNEL node
#' @examples
#' \dontrun{
#  acsfile <- system.file("extdata/cytobank_experiment.acs", package = "CytoML")
#  ce <- open_cytobank_experiment(acsfile)
#  xmlfile <- ce$gatingML
#' g <- read.gatingML.cytobank(xmlfile)
#' getChildren(g, "GateSet_722326")
#' getParent(g, "GateSet_722326")
#' }
#' @importClassesFrom methods character ANY data.frame environment list logical matrix missing numeric oldClass
#' @importFrom flowWorkspace getChildren
setMethod("getChildren", signature = c("graphGML", "character"),
          definition = function(obj, y) {
  edges(obj, y)[[1]]
})

#' get parent nodes
#'
#' @param obj \code{graphGML}
#' @param y \code{character} child node path
#' @return a graphNEL node
#' @importFrom flowWorkspace gs_pop_get_parent getParent
setMethod("getParent", signature = c("graphGML", "character"),
          definition = function(obj, y) {

   inEdges(y, obj)[[1]]

})

#' get gate from the node
#'
#' @param obj \code{graphGML}
#' @param y \code{character} node path
#' @return the gate information associated with the node
#' @importFrom flowWorkspace getGate
setMethod("getGate", signature = c("graphGML", "character"),
          definition = function(obj, y) {

          nodeData(obj, y)[["gateInfo"]]
})


#' show method for graphGML
#'
#' show method for graphGML
#'
#' @param object \code{graphGML}
#' @return nothing
#' @importFrom methods show
setMethod("show", signature = c("graphGML"),
          definition = function(object) {
  cat("--- Gating hieararchy parsed from GatingML: ")

  cat("\n")
  cat("\twith ", length(object@nodes), " populations defined\n")
})


#' plot the population tree stored in graphGML.
#'
#' The node with dotted order represents the population that has tailored gates (sample-specific gates) defined.
#'
#' @param x a graphNEL generated by constructTree function
#' @param y not used
#' @param label specifies what to be dispaled as node label. Can be either 'popName' (population name parsed from GateSets) or 'gateName'(the name of the actual gate associated with each node)
#' @return nothing
#' @importFrom graph nodeData nodes<- nodeRenderInfo<-
#' @importFrom Rgraphviz renderGraph layoutGraph
#' @examples
#' \dontrun{
#  acsfile <- system.file("extdata/cytobank_experiment.acs", package = "CytoML")
#  ce <- open_cytobank_experiment(acsfile)
#  xmlfile <- ce$gatingML
#' g <- read.gatingML.cytobank(xmlfile)
#' plot(g)
#'}
setMethod("plot", signature = c(x = "graphGML", y = "missing"), definition = function(x, y = "missing", label = c("popName", "gateName")){
  label <- match.arg(label, c("popName", "gateName"))
  if(label == "popName")
    nodeLabel  <- sapply(nodeData(x), `[[`, "popName")
  else
    nodeLabel  <- sapply(nodeData(x), function(i)i[["gateInfo"]][["gateName"]])


  #annotate the node with tailor gate info
  nTailoredGate <- sapply(nodeData(x), function(i)length(i[["gateInfo"]][["tailored_gate"]]))

  nAttrs <- list()

  nAttrs$label <- nodeLabel

  nAttrs$lty <- sapply(nTailoredGate
                       ,function(i)
                       {
                         ifelse(i>0,"dotted","solid")
                       })

  nodeRenderInfo(x) <- nAttrs
  lay <- layoutGraph(x
                     ,attrs=list(graph=list(rankdir="LR",page=c(8.5,11))
                                 ,node=list(fixedsize=FALSE
                                            ,fontsize = 12
                                            ,shape="ellipse"
                                 )
                     )
  )
  renderGraph(lay)

})

#' Apply the gatingML graph to a GatingSet
#'
#' It applies the gates to the GatingSet based on the population tree described in graphGML.
#'
#' @param x graphGML
#' @param y GatingSet
#' @param ... other arguments
#' @return
#' Nothing. As the side effect, gates generated by gating methods are saved in \code{GatingSet}.
#' @noRd 
#' @importFrom flowWorkspace gs_pop_set_name gs_pop_get_children recompute sampleNames gs_pop_add
#' @importFrom RBGL tsort
gating_graphGML <- function(x, y, trans = NULL, ...) {
  if(is.null(trans))
  trans <- getTransformations(x)

  gt_nodes <- tsort(x)
  for (nodeID in gt_nodes) {

    # get parent node to gate
    gt_node <- getNodes(x, nodeID, only.names = FALSE)
    popName <- gt_node[["popName"]]


    parentID <- getParent(x, nodeID)

    if(length(parentID) == 0)
      parent <- "root"
    else{
      parent <- .getPath(x, parentID)
    }


    gs_nodes <- basename(gs_pop_get_children(y[[1]], parent))

    if (length(gs_nodes) == 0)
      isGated <- FALSE
    else
      isGated <- any(popName %in% gs_nodes)

    #TODO: rename the node name with path in order to match against gs
#     parentInd <- match(parent, getNodes(y[[1]], showHidden = TRUE))
#     if (is.na(parentInd))
#       stop("parent node '", parent, "' not gated yet!")
    if(isGated){
      message("Skip gating! Population '", paste(popName, collapse = ","), "' already exists.")
      next
    }
    message(popName)
    gateInfo <- gt_node[["gateInfo"]]
    this_gate <- gateInfo[["gate"]]



    # transform bounds if applicable
    bound <- gateInfo[["bound"]]
    if(!is.null(trans))
    {
      for(rn in rownames(bound)){
        thisTrans <- trans[[rn]]
        if(!is.null(thisTrans))
          bound[rn, ] <- thisTrans[["transform"]](unlist(bound[rn, ]))
      }
    }

        # if(popName == "MDSC(gran-cd15+)")
        #   browser()
    this_gate <- extend(this_gate,bound = bound)

    sn <- sampleNames(y)
    this_gate <- sapply(sn, function(i)this_gate)

    #update gates that are tailored for specific samples
    tailor_gate <- gateInfo[["tailored_gate"]]
    #lookup by fcs name|fileid
    tg_idx <- tailor_gate[["file_vs_gateid"]][sn]
    tg_idx <- tg_idx[!is.na(tg_idx)]
    dup <- duplicated(tg_idx)
    if(any(dup))
      stop("Unexpected behavior!The same tailor gate is matched by both file id and file name!", paste(names(tg_idx[dup]), collapse = " "))
    if(length(tg_idx) > 0){
      this_tgs <- lapply(tailor_gate[["gateid_vs_gate"]][tg_idx], extend,bound = bound)
      tg_sn <- names(tg_idx)
      this_gate[tg_sn] <- this_tgs
    }




    gs_pop_add(y, this_gate, parent = parent, name = popName)

  }
  recompute(y)
}

#' Extract compensation from graphGML object.
#' @param x graphGML
#' @return compensation object or "FCS" when compensation comes from FCS keywords
#' @importFrom flowWorkspace getCompensationMatrices
#' @method getCompensationMatrices graphGML
getCompensationMatrices.graphGML <- function(x){
  x@graphData[["compensation"]]

}

#' Extract transformations from graphGML object.
#' @param x graphGML
#' @param ... not used
#' @return transformerList object
#' @importFrom flowCore eval parameters colnames
#' @importFrom flowWorkspace transformerList asinh_Gml2 flow_trans asinhtGml2_trans logicleGml2_trans logtGml2_trans getTransformations
#' @importFrom methods extends
#' @method getTransformations graphGML
getTransformations.graphGML <- function(x, ...){
  trans <- x@graphData[["transformations"]]
  if(!is.null(trans)){
    chnls <- names(trans)
    trans <- sapply(trans, function(thisTrans){
      trans
      #convert from transform object to function since transform has empty function in .Data slot
      #which is not suitable for transformList constructor
      # trans.fun <- eval(thisTrans)
      trans.type <- class(thisTrans)
      if(methods::extends(trans.type, "asinhtGml2")){
        # inv.func <- asinh_Gml2(thisTrans@T, thisTrans@M, thisTrans@A, inverse = TRUE)
        trans.obj <- asinhtGml2_trans(thisTrans@T, thisTrans@M, thisTrans@A)
      }else if(methods::extends(trans.type, "logicletGml2")){
        trans.obj <- logicleGml2_trans(thisTrans@T, thisTrans@M, thisTrans@W, thisTrans@A)
      }else if(methods::extends(trans.type, "logtGml2")){
        trans.obj <- logtGml2_trans(thisTrans@T, thisTrans@M)
      }else
        stop("Don't know how to inverse transformation: ", trans.type)

      # trans.obj <- flow_trans(trans.type, trans.fun, inv.func)
      trans.obj
    }
    , USE.NAMES = FALSE, simplify = FALSE)

    trans <- transformerList(chnls, trans)
  }
  trans
}

#' compensate a GatingSet based on the compensation information stored in graphGML object
#'
#'
#' @param x GatingSet
#' @param spillover graphGML
#' @param ... unused.
#' @return compensated GatingSet
#' @importFrom flowCore compensate keyword
#' @importFrom flowWorkspace gs_pop_get_data
setMethod("compensate", signature = c("GatingSet", "graphGML"), function(x, spillover, ...){

  comp <- getCompensationMatrices(spillover)
  if(is(comp, "compensation")){
    # prefix <- TRUE
    skip <- FALSE
  }else if(comp == "FCS"){
    # prefix <- FALSE
    fs <- gs_pop_get_data(x)
    fr <- fs[[1, use.exprs = FALSE]]
    #can't use spillover method directly because it will error out when none is found
    mat <- keyword(fr, c("spillover", "SPILL"))
    mat <- compact(mat)
    if(length(mat) == 0){
      skip <- TRUE
      warning("Compensation is skipped!Because gates refer to 'FCS' for compensation but no spillover is found in FCS.")
    }else{
      skip <- FALSE
      mat <- mat[[1]]
      comp <- compensation(mat)
    }
  }else if(comp == "NONE")
    skip <- TRUE
  if(skip)
    return(x)
  else{
    x <- compensate(x, comp)

#     if(prefix){
#
#       comp_param <- colnames(comp@spillover)
#       #strip prefix
#       comp_param <- sapply(comp_param, function(i)sub("(^Comp_)(.*)", "\\2", i), USE.NAMES = FALSE)
#       #match to chnls
#       chnls <- colnames(x)
#       ind <- match(comp_param, chnls)
#       chnls[ind] <- paste0("Comp_", chnls[ind])
#       colnames(x) <- chnls
#     }
    return(x)
  }

})