Question

custom gatingMethod not functioning / registering as expected (openCyto) (minDensity)

0

Entering edit mode

k.teijgeler • 0

@kteijgeler-15773

Last seen 5.2 years ago

Dear reader,

I've been happily using openCyto for a while now and i really like the mindensity gatingMethod. It's very robust in panels/samples that have 2 peaks for a marker (+- 80-90% of the samples for this particular marker, figure 1) .

figure 1

However i've experienced some cases in which i would like an alternative logic when the peakfinder finds just 1 peak or finds no valleys (+- 10-20% of the samples for this particular marker, figure 2).

figure 2

Im experimenting with different strategies how to tackle this problem. Currently instead of using GateRange1 or GateRange2 when there is 1 peak or no valleys im attempting to get the median and add 1*SD and put a gate there (figure 2).

I've got this working in a small script. The gating method will register, but when i attempt to gate from a .csv gatingfile i get the following error:

"[232] "Gating for 'CD4posCD45RApos'"                                                                                                                                             
[233] "Warning in .gating_gtMethod(x, y, ...) : NAs introduced by coercion"                                                                                                      
[234] "Running in parallel mode with 1 cores."                                                                                                                                   
[235] "Error in (function (fs, pp_res, gFunc, popAlias, channels, gFunc_args)  : "                                                                                               
[236] "  failed at 650583660906 mock_a8_a08_075.fcs650583660906 mock_b8_b08_076.fcs650583660906 rsv_d8_d08_095.fcs650583660906 rsv_e8_e08_096.fcs650583660906 seb_g8_g08_108.fcs"
[237] "Error in abs(i) : non-numeric argument to mathematical function"                                                                                                          
[238] ""

I've attached the code below, in case the uploaded file is no longer available.: If possible i would like some help figuring out why this function is not working as expected (atleast how i expected ;) ). Please tell me if you need more information.

Links:

csv Gatingfile, compensation file, fcs files and gatingfunction

SessionInfo:

Kind regards, Kasper

mindensityNoValley <-function (fr, channel, filterId = "", positive = TRUE, pivot = FALSE, 
                                gate_range = NULL, min = NULL, max = NULL, peaks = NULL, 
                                ...) 
{
  if (missing(channel) || length(channel) != 1) {
    stop("A single channel must be specified.")
  }
  if (!(is.null(min) && is.null(max))) {
    fr <- openCyto:::.truncate_flowframe(fr, channels = channel, min = min, 
                              max = max)
  }
 x <- exprs(fr)[, channel]
  if (is.null(peaks)) 
    peaks <- openCyto:::.find_peaks(x, ...)[, "x"]
  if (is.null(gate_range)) {
    gate_range <- c(min(x), max(x))
  }
  else {
    gate_range <- sort(gate_range)
  }
  if (length(peaks) == 1) {
    ### Original content of mindensity
    #cutpoint <- ifelse(positive, gate_range[1], gate_range[2])

    ### new content:
    ### calculate KDE
    ### get median and SD
    ### gate = median + 1 * SD
    KDEdata <- density(x)                       # new
    KDEdatamedian <- median(KDEdata)            # new                                          
    KDEdataStdev <- Stdev <- sd( x = x)         # new         
    cutpoint <- KDEdatamedian + KDEdataStdev    # new  
    }

  else {
    valleys <- try(.find_valleys(x, ...), silent = TRUE)
    valleys <- openCyto:::.between_interval(x = valleys, interval = gate_range)
    if anyis.na(valleys))) {
      ### Original content of mindensity
      #cutpoint <- ifelse(positive, gate_range[1], gate_range[2])

      ### new content:
      ### calculate KDE
      ### get median and SD
      ### gate = median + 1 * SD
      KDEdata <- density(x)                       # new
      KDEdatamedian <- median(KDEdata)            # new                                          
      KDEdataStdev <- Stdev <- sd( x = x)         # new         
      cutpoint <- KDEdatamedian + KDEdataStdev    # new
    }
    else if (length(valleys) == 1) {
      cutpoint <- as.vector(valleys)
    }
    else if (length(valleys) > 1) {
      peaks <- sort(peaks[1:2])
      cutpoint <- openCyto:::.between_interval(valleys, peaks)[1]
      if is.na(cutpoint)) {
        cutpoint <- valleys[1]
      }
    }
  }
  gate_coordinates <- list(c(cutpoint, Inf))
  names(gate_coordinates) <- channel
  rectangleGate(gate_coordinates, filterId = filterId)
}

.mindensityNoValley <- function(fr, pp_res, channels=NA, ...){

  my_gate <- mindensityNoValley(fr[pp_res,],channel=channels, filter_id=filterId, ...)
  return(my_gate)
}
registerPlugins( fun = .mindensityNoValley, methodName = "mindensityNoValley" , type = "gating")

openCyto registerPlugin gatingMethod minDensity • 1.8k views

ADD COMMENT • link 5.2 years ago k.teijgeler • 0

score 1 · Answer 1 · 2019-08-08

1

Entering edit mode

Jake Wagner ▴ 310

@jake-wagner-19995

Last seen 4.0 years ago

I'll gladly look in to your function quickly to see what the issue might be, but have you also checked out tailgate (it will be renamed gate_tail in the upcoming API changeover). It uses the derivative of the KDE to handle this situation (one peak, so no valley). One of the vignettes briefly touches on it (https://www.bioconductor.org/packages/release/bioc/vignettes/openCyto/inst/doc/HowToAutoGating.html) but I'm happy to discuss parameter tuning a bit more. I also discussed it a bit on a Github issue recently: https://github.com/RGLab/openCyto/issues/191.

ADD COMMENT • link 5.2 years ago Jake Wagner ▴ 310

0

Entering edit mode

Thanks for the quick answer!

Good idea! I use the tailgate for other nodes but it might be better to just incoorporate a tailgate when only 1 peak and no valleys are found. Great :).

To test whether i could register a new plugin i did the following: - Made a copy function of mindensity called mindensity3 - Changed nothing

This does register and work from the gating.csv file as expected.

Conclusion: There is definitely something in the code i added myself that is causing the error.

Used Code:

mindensity3 <- function (fr, channel, filterId = "", positive = TRUE, pivot = FALSE, 
          gate_range = NULL, min = NULL, max = NULL, peaks = NULL, 
          ...) 
{
  if (missing(channel) || length(channel) != 1) {
    stop("A single channel must be specified.")
  }
  if (!(is.null(min) && is.null(max))) {
    fr <- openCyto:::.truncate_flowframe(fr, channels = channel, min = min, 
                              max = max)
  }
  x <- exprs(fr)[, channel]
  if (is.null(peaks)) 
    peaks <- openCyto:::.find_peaks(x, ...)[, "x"]
  if (is.null(gate_range)) {
    gate_range <- c(min(x), max(x))
  }
  else {
    gate_range <- sort(gate_range)
  }
  if (length(peaks) == 1) {
    cutpoint <- ifelse(positive, gate_range[1], gate_range[2])
  }
  else {
    valleys <- try(openCyto:::.find_valleys(x, ...), silent = TRUE)
    valleys <- openCyto:::.between_interval(x = valleys, interval = gate_range)
    if anyis.na(valleys))) {
      cutpoint <- ifelse(positive, gate_range[1], gate_range[2])
    }
    else if (length(valleys) == 1) {
      cutpoint <- as.vector(valleys)
    }
    else if (length(valleys) > 1) {
      peaks <- sort(peaks[1:2])
      cutpoint <- openCyto:::.between_interval(valleys, peaks)[1]
      if is.na(cutpoint)) {
        cutpoint <- valleys[1]
      }
    }
  }
  gate_coordinates <- list(c(cutpoint, Inf))
  names(gate_coordinates) <- channel
  rectangleGate(gate_coordinates, filterId = filterId)
}

.mindensity3 <- function (fr, pp_res = NULL, channels, ...) 
{
  if (length(channels) != 1) 
    stop("invalid number of channels for mindensity!")
  gate <- mindensity3(fr, channel = channels, ...)
  gate
}

registerPlugins( fun = .mindensity3, methodName = "mindensity3" )

ADD REPLY • link 5.2 years ago k.teijgeler • 0

0

Entering edit mode

To further diagnose my peak patterns i thought of adding the following (see code below). Is there a way to get the name of the alias or other information from the gatingFile into context in the function?

KDEdata <- density(x)                        # new
KDEdatamedian <- median(x)
KDEdatamean <- mean(x)
KDEdataStdev <- sd( x = x)          # new         

line <- paste0( "Samplename =: " , fr@description$sampleNames[1] )
write(line,file="C:/temp/output.txt",append=TRUE)
line <- paste0( "Channel =: " , channel )
write(line,file="C:/temp/output.txt",append=TRUE)
line <- paste0( "Peaks =: " , peaks )
write(line,file="C:/temp/output.txt",append=TRUE)
line <- paste0( "Mean =: " , KDEdatamean )
write(line,file="C:/temp/output.txt",append=TRUE)
line <- paste0( "Median =: " , KDEdatamedian )
write(line,file="C:/temp/output.txt",append=TRUE)

ADD REPLY • link 5.2 years ago k.teijgeler • 0

score 0 · Answer 2 · 2019-08-09

A collegue told me to put some print() calls in the function and guess what i was referring to an incorrect object. Oeps;)

KDEdata <- density(x)                       
KDEdatamedian <- median(KDEdata)                     # median of KDEdata instead of median of (x)                                           
KDEdataStdev <- Stdev <- sd( x = x)                
cutpoint <- KDEdatamedian + KDEdataStdev

New code would then be:

KDEdata <- density(x)                       
KDEdatamedian <- median(x)                           #new                 
KDEdataStdev <- sd( x = x)                   
cutpoint <- KDEdatamedian + KDEdataStdev

Now its working :). But im still more a fan of your idea @Jake Wagner. I will see if i can implement a tailGate like approach instead. Will come back on that later!