xmlElementsByTagName <- function(...)suppressWarnings(XML::xmlElementsByTagName(...))
compact <- function (x)
{
null <- vapply(x, is.null, logical(1))
x[!null]
}
#' the parameter range from the flow data associated with GatingHierarchy
#' @param ... GatingHierarchy object
#' @param na.rm not used
#' @param type character of "instrument" or "data" indicating whether to retrieve the instrument or the actual data range
#' @param raw.scale logical whether convert the range from transformed scale to raw scale
#' @return matrix
#' @examples
#' \dontrun{
#' range(gh, type = "data")#return data range
#' range(gh) #return instrument range
#' range(gh, raw.scale = TRUE) #inverse transform the range to the raw scale
#' }
range.GatingHierarchy <- function(..., na.rm = FALSE, type = c("instrument", "data"), raw.scale = FALSE){
type <- match.arg(type)
gh <- list(...)[[1]]
fr <- gh_pop_get_data(gh, use.exprs = type == "data")
rng <- range(fr, type = type)
if(raw.scale)
{
chnls <- names(rng)
translist <- gh_get_transformations(gh, only.function = FALSE)
for(chnl in chnls)
{
trans <- translist[[chnl]]
if(!is.null(trans))
{
rng[, chnl] <- trans$inverse(rng[, chnl])
}
}
}
rng
}
is.cytof <- function(gs){
any(grepl("Event_length", colnames(gs)))
}
######################################
#common APIs related to processing comp, trans, gates
#to prepare the GatingML output that can be shared by both cytobank and flowJo modules
######################################
#' @importFrom flowWorkspace gh_pop_get_data
#' @importFrom flowCore compensation identifier identifier<- compensatedParameter asinhtGml2 logicletGml2 logtGml2
export_comp_trans <- function(gs, flowEnv, cytobank.default.scale = FALSE, type = c("cytobank", "flowJo"))
{
type <- match.arg(type)
#parse comp and channel names
comp <- gh_get_compensations(gs[[1]])#assuming the comp is identical across samples
if(is.null(comp)){
#no compensation and channel names in transformation are not prefixed
chnls <- colnames(gs_pop_get_data(gs))
prefix_chnls_orig <- chnls
}else{
chnls <- as.vector(parameters(comp))
#retrieve the prefix for latter trans matching
cmp <- gs_get_compensation_internal(gs@pointer, sampleNames(gs)[[1]])
prefix <- cmp$prefix
suffix <- cmp$suffix
prefix_chnls_orig <- paste0(prefix, chnls, suffix)
}
#add comp
if(!is.null(comp)){
if(type == "cytobank"){
#prefix the channels since cytobank expect that
prefix_chnls <- paste0("Comp_", chnls)
rownames(comp@spillover) <- prefix_chnls #change fluorochromes and leave detector(colnames) unchanged
}else{
#prefix the channels
prefix_chnls <- paste0("Comp-", chnls)
# prefix_chnls <- chnls #flowJo use the raw channel name
# rownames(comp@spillover) <- prefix_chnls #change fluorochromes
# markers <- markernames(gs)
# rownames(comp@spillover) <- markers #change detector to marker names because flowJo seems to expect that
}
comp <- compensation(comp@spillover)
compId <- identifier(comp)
compId <- paste("Spill", compId, sep = "_")
identifier(comp) <- compId
flowEnv[[compId]] <- comp
}else{
compId <- "FCS"
prefix_chnls <- chnls
}
#add trans (assume it is identical across samples)
trans <- gh_get_transformations(gs[[1]], only.function = FALSE)
if(length(trans) == 0)
stop("no transformation is found in GatingSet!")
trans.Gm2objs <- list()
rescale.gate <- FALSE
for(transName in names(trans)){
trans.obj <- trans[[transName]]
type <- trans.obj[["name"]]
# ind <- sapply(chnls, grepl, x = transName, USE.NAMES = FALSE)
ind <- prefix_chnls_orig == transName#do strict match due to the possible mismatch for cases like CD3 vs CD33
nMatched <- sum(ind)
if(nMatched > 1)
stop("More than one channels matched to transformation: ", transName)
else if(nMatched == 1){
trans.func <- trans.obj[["transform"]]
chnl <- chnls[ind]
prefix_chnl <- prefix_chnls[ind]
prefix_chnl_orig <- prefix_chnls_orig[ind]
# if(is.null(comp))
# param.obj <- prefix_chnl
# else
param.obj <- compensatedParameter(prefix_chnl
, spillRefId = compId
, searchEnv = flowEnv
, transformationId = prefix_chnl)
if(type %in% c("asinhtGml2", "flowJo_fasinh")){
#extract parameters
env <- environment(trans.func)
transID <- paste0("Tr_Arcsinh_", prefix_chnl)
flowEnv[[transID]] <- asinhtGml2(parameters = param.obj
, M = env[["m"]]
, T = env[["t"]]
, A = env[["a"]]
, transformationId = transID
)
}else if(type == "logicleGml2"){
#extract parameters
env <- environment(trans.func)
transID <- paste0("Tr_logicleGml2_", prefix_chnl)
flowEnv[[transID]] <- logicletGml2(parameters = param.obj
, M = env[["M"]]
, T = env[["T"]]
, A = env[["A"]]
, W = env[["W"]]
, transformationId = transID
)
}else if(type %in% c("logicle", "flowJo_logicle")){
#extract parameters
env <- environment(trans.func)
transID <- paste0("Tr_logicle_", prefix_chnl)
flowEnv[[transID]] <- logicletGml2(parameters = param.obj
, M = env[["m"]]
, T = env[["t"]]
, A = env[["a"]]
, W = as.vector(env[["w"]])
, transformationId = transID
)
rescale.gate <- TRUE
}else if(type %in% c("flowJo_biexp")){
transID <- paste0("Tr_logicle_", prefix_chnl)
param <- attr(trans.func,"parameters")
flowEnv[[transID]] <- logicletGml2(parameters = param.obj
, M = param[["pos"]]
, T = param[["maxValue"]]
, A = param[["neg"]]
, W = log10(-param[["widthBasis"]]) / 2
, transformationId = transID
)
rescale.gate <- TRUE
}else if(type == "flowJo_flog"){
env <- environment(trans.func)
transID <- paste0("Tr_flog_", prefix_chnl)
flowEnv[[transID]] <- logtGml2(parameters = param.obj
, M = 1
, T = 1
, transformationId = transID
)
rescale.gate <- TRUE
}else if(type == "logtGml2"){
stop("logtGml2 is not supported yet!")#there is some issue associated with gate truncation at negative side during exporting, need to further troubleshoot
env <- environment(trans.func)
transID <- paste0("Tr_flog_", prefix_chnl)
flowEnv[[transID]] <- logtGml2(parameters = param.obj
, M = env[["decade"]]
, T = env[["max_val"]]
, boundMin = env[["min_val"]]
, transformationId = transID
)
rescale.gate <- TRUE
}else{
# browser()
stop("unsupported trans: ", type)
}
if(cytobank.default.scale){
rescale.gate <- TRUE
#overwrite the customed scale with default one
flowEnv[[transID]] <- asinhtGml2(parameters = param.obj
, M = 0.43429448190325176
, T = ifelse(is.cytof(gs), 5.8760059682190064, 176.2801790465702)
, A = 0.0
, transformationId = transID)
}
#save another copy of trans.obj in the list
trans.Gm2objs[[prefix_chnl_orig]] <- flowEnv[[transID]]
}
}
return(list(trans.Gm2objs = trans.Gm2objs, trans = trans, compId = compId))
}
#' @importFrom flowCore eval
processGate <- function(gate, gml2.trans, compId, flowEnv, rescale.gate = FALSE, orig.trans){
params <- as.vector(parameters(gate))
prefixed_chnls_orig <- names(gml2.trans)
for(i in seq_along(params)){
param <- params[i]
# ind <- sapply(chnls, function(chnl)grepl(chnl, param), USE.NAMES = FALSE)
ind <- prefixed_chnls_orig == param
nMatched <- sum(ind)
if(nMatched == 0){
chnl <- gate@parameters[[i]]@parameters
#can't use "uncompensated" because it will cause multiple entries when paring gate back into openCyto through cytobank_to_gatingset
gate@parameters[[i]] <- compensatedParameter(chnl
, spillRefId = compId #"uncompensated"
, searchEnv = flowEnv
, transformationId = chnl
)
}else if(nMatched == 1){
prefixed_chnl_orig <- prefixed_chnls_orig[ind]
orig.trans.obj <- orig.trans[[which(ind)]]
gml2.trans.obj <- gml2.trans[[prefixed_chnl_orig]]
if(rescale.gate){
inv.fun <- orig.trans.obj[["inverse"]]
trans.fun <- eval(gml2.trans.obj)
#rescale
gate <- rescale_gate(gate, inv.fun, param)
gate <- rescale_gate(gate, trans.fun, param)
}
#attach trans reference
gate@parameters[[i]] <- gml2.trans.obj
}else if(nMatched > 1)
stop("multiple trans matched to :", param)
}
# round
gate
}
inverse <- function(gate, boundary){
UseMethod("inverse")
}
inverse.polygonGate <- function(gate, ...){
gate@boundaries <- inverse.polygon(gate@boundaries, ...)
gate
}
inverse.ellipsoidGate <- function(gate, ...){
inverse(as(gate, "polygonGate"), ...)
}
inverse.rectangleGate <- function(gate, ...){
param <- as.vector(parameters(gate))
if(length(param) == 1){
stop("to be implemented!")
}else
inverse(as(gate, "polygonGate"), ...)
}
inverse.polygon <- function(mat, boundary){
dims <- colnames(mat)
stopifnot(all(dims %in% colnames(boundary)))
#enlarge the boundary to ensure the negated polygon includes all events outside of original polygon
boundary["min", ] <- boundary["min", ] - 1e4
boundary["max", ] <- boundary["max", ] + 1e4
x <- dims[1]
y <- dims[2]
#find the top most point as starting point
ind <- which.max(mat[, y])
#reorder mat starting from top most
nRow <- nrow(mat)
orig.seq <- seq_len(nRow)
seq1 <- ind:nRow
seq2 <- setdiff(orig.seq, seq1)
new.seq <- c(seq1, seq2)
new.seq
mat.new <- mat[new.seq, ]
mat.inv <- mat.new[1, ]
#extend to upper bound
pt1 <- c(mat.new[1, x], boundary["max", y])
mat.inv <- rbind(mat.inv, pt1, deparse.level = 0)
#left top
mat.inv <- rbind(mat.inv, c(boundary["min", x], boundary["max",y]))
#left bottom
mat.inv <- rbind(mat.inv, c(boundary["min", x], boundary["min",y]))
#right bottom
mat.inv <- rbind(mat.inv, c(boundary["max", x], boundary["min",y]))
#right top
mat.inv <- rbind(mat.inv, c(boundary["max", x], boundary["max",y]))
#back to first extended point
mat.inv <- rbind(mat.inv, pt1, deparse.level = 0)
#connect to the remain mat.new in reverse order
mat.inv <- rbind(mat.inv, mat.new[rev(seq_len(nRow)[-1]), ])
# plot(mat.inv, type = "n")
# polygon(mat.new, col = "red")
# polygon(mat.inv, col = "blue")
mat.inv
}
