search

RGLab / cytolib

c++ library for representing and interacting the gated cytometry data structure
GNU Affero General Public License v3.0
12 stars 11 forks source link

Error Gating Empty Populations #22

Closed DillonHammill closed 4 years ago

DillonHammill commented 4 years ago

Hi @mikejiang & @jacobpwagner ,

Not sure if this is the best place to post this issue. It looks like there are issues gating empty populations with the latest cytolib, flowWorkspace and openCyto.

> gt_gating(gt, gs)
Preprocessing for 'Cells'
Gating for 'Cells'
done!
done.
Preprocessing for 'Single Cells'
Gating for 'Single Cells'
done!
done.
Preprocessing for 'Dead Cells'
Gating for 'Dead Cells'
done!
done.
Live Cells gating...
done!
done.
Preprocessing for 'T Cells'
Gating for 'T Cells'
done!
done.
Preprocessing for 'CD8 T Cells'

error: Mat::max(): object has no elements
Error in get_cytoset_from_node(obj@pointer, y) : 
  Mat::max(): object has no elements

Thanks for your help.

Dillon

DillonHammill commented 4 years ago

Here is the traceback:

> traceback()
14: stop(list(message = "Mat::max(): object has no elements", call = get_cytoset_from_node(obj@pointer, 
        y), cppstack = NULL))
13: get_cytoset_from_node(obj@pointer, y)
12: initialize(value, ...)
11: initialize(value, ...)
10: new("cytoset", pointer = get_cytoset_from_node(obj@pointer, y))
9: gs_pop_get_data(y, parent)
8: .preprocessing(x, y, ...)
7: preprocessing(x = this_ppm, y, parent = parent, gtPop = gt_node_pop, 
       gm = this_gate, ...)
6: preprocessing(x = this_ppm, y, parent = parent, gtPop = gt_node_pop, 
       gm = this_gate, ...)
5: .gating_gatingTemplate(x, y, ...)
4: gt_gating.gatingTemplate(gt, gs)
3: gt_gating(gt, gs)
2: withCallingHandlers(expr, warning = function(w) invokeRestart("muffleWarning"))
1: gt_gating(gt, gs)
jacobpwagner commented 4 years ago

Thanks for the traceback @DillonHammill. I'm looking in to it now.

DillonHammill commented 4 years ago

Thanks @jacobpwagner. I will give it a try later today.

DillonHammill commented 4 years ago

Looks good! Thanks @jacobpwagner.

DillonHammill commented 4 years ago

Sorry @jacobpwagner, looks like there are still some problems...

> gt_gating(gt,gs)
Preprocessing for 'Cells'
Gating for 'Cells'
done!
done.
Preprocessing for 'Single Cells'
Gating for 'Single Cells'
done!
done.
Preprocessing for 'Dead Cells'
Gating for 'Dead Cells'
done!
done.
Live Cells gating...
done!
done.
Preprocessing for 'T Cells'
Gating for 'T Cells'
done!
done.
Preprocessing for 'CD8 T Cells'
Gating for 'CD8 T Cells'
 Error in cf_getData(object@pointer) : 
  length of 'dimnames' [2] not equal to array extent 
> traceback()
22.
cf_getData(object@pointer) 
21.
exprs(fr) 
20.
exprs(fr) 
19.
cytoframe_to_flowFrame(fr) 
18.
.local(x, i, j, ...) 
17.
from[[1]] 
16.
from[[1]] 
15.
asMethod(object) 
14.
as(fs, "flowFrame") 
13.
(function (fs, pp_res, gFunc, popAlias, channels, gFunc_args) 
{
    require(openCyto)
    sn <- sampleNames(fs) ... 
12.
mapply(list(Activation_1.fcs = new("cytoset", pointer = <pointer: 0x00000249a1b6f610>, 
    frames = <environment>, phenoData = new("AnnotatedDataFrame", 
        varMetadata = structure(list(labelDescription = character(0)), row.names = character(0), class = "data.frame"), 
        data = structure(list(), .Names = character(0), row.names = integer(0), class = "data.frame"),  ... 
11.
eval(thisCall) 
10.
eval(thisCall) 
9.
.gating_gtMethod(x, y, ...) 
8.
gt_gating.gtMethod(x = this_gate, y, parent = parent, gtPop = gt_node_pop, 
    pp_res = pp_res, ...) 
7.
gt_gating(x = this_gate, y, parent = parent, gtPop = gt_node_pop, 
    pp_res = pp_res, ...) 
6.
.gating_gatingTemplate(x, y, ...) 
5.
gt_gating.gatingTemplate(gt, x, ...) 
4.
gt_gating(gt, x, ...) at cyto_gatingTemplate-helpers.R#315
3.
withCallingHandlers(expr, warning = function(w) invokeRestart("muffleWarning")) 
2.
suppressWarnings(gt_gating(gt, x, ...)) at cyto_gatingTemplate-helpers.R#315
1.
cyto_gatingTemplate_apply(gs, "Activation-gatingTemplate.csv")
jacobpwagner commented 4 years ago

I can look in to it, but because I added data members to the CytoFrameView class, you'll have to rebuild every package that builds against cytolib. Otherwise you will get some very weird behavior. So (if you haven't already done it), please try rebuilding flowCore, flowWorkspace, openCyto, and CytoML (after rebuilding cytolib of course) and seeing if that fixes the problem.

DillonHammill commented 4 years ago

I re-built cytolib and then flowCore and flowWorkspace but I did not re-build openCyto, I will do this now...

DillonHammill commented 4 years ago

Still seeing the same error after:

devtools::install_github("RGLab/cytolib", ref = "trunk", force = TRUE)
devtools::install_github("RGLab/flowCore", ref = "trunk", force = TRUE)
devtools::install_github("RGLab/flowWorkspace", ref = "trunk", force = TRUE)
devtools::install_github("RGLab/openCyto", ref = "trunk", force = TRUE)
devtools::install_github("RGLab/CytoML", ref = "trunk", force = TRUE)
jacobpwagner commented 4 years ago

I'm trying to reproduce it. Is the starting GatingSet from CytoRSuiteData or CytoExploreRData? Can you point me to a reprex from one of your docs or tests?

jacobpwagner commented 4 years ago

I'm checking this snippet, based on CytoRSuite/tests/testthat/helper-lib.R and just adding in the required library calls and it is successful in my hands:

library(flowCore)
library(flowWorkspace)
library(CytoRSuite)

library(CytoRSuiteData)

# Load in vdiffr for image comparison -
library(vdiffr)

# Load in grDevices for recordPlot() -
library(grDevices)

# Load in robustbase (colMedians) -
library(robustbase)

# Data directory CytoRSuiteData ------------------------------------------------
datadir <- system.file("extdata", package = "CytoRSuiteData")

# Assign Activation dataset to fs -
fs <- read.flowSet(path = paste0(datadir,"/Activation"), pattern = ".fcs")

# Sample names
nms <- sampleNames(fs)

# Sample for speed -
fs <- flowSet(lapply(seq_len(length(fs)), function(x) {
    fs[[x]][1:2000, ]
}))
sampleNames(fs) <- nms
pData(fs)$name <- nms

# Samplenames
nms <- sampleNames(fs)

# Extract channels
chans <- colnames(fs)

# pData information -
pData(fs)$OVAConc <- c(0, 0.005, 0.05, 0.5)
pData(fs)$Treatment <- c("A", "A", "B", "B")
pData(fs)$Treatment <- factor(pData(fs)$Treatment, levels = c("B", "A"))
chnls <- c("Alexa Fluor 405-A", 
           "Alexa Fluor 430-A", 
           "APC-Cy7-A", "PE-A", 
           "Alexa Fluor 488-A", 
           "Alexa Fluor 700-A", 
           "Alexa Fluor 647-A", 
           "7-AAD-A")
markers <- c("Hoechst-405", 
             "Hoechst-430", 
             "CD11c", 
             "Va2", 
             "CD8", 
             "CD4", 
             "CD44", 
             "CD69")
names(markers) <- chnls
markernames(fs) <- markers

# GatingSet -
gs <- GatingSet(fs)

# Compensation -
gs <- compensate(gs, fs[[1]]@description$SPILL)

# Transformation -
trans <- estimateLogicle(gs[[4]], cyto_fluor_channels(gs))
gs <- transform(gs, trans)

# gatingTemplate -
gt <- gatingTemplate(paste0(datadir,"/Activation-gatingTemplate.csv"))

# gatingTemplate file -
gtf <- read.csv(system.file("extdata", 
                            "Activation-gatingTemplate.csv", 
                            package = "CytoRSuiteData"))

# Gating -
gating(gt, gs)

Checking the stats also shows it getting past the empty parent issue that was the crux of the issue before (the child of the empty CD8 T Cells subpop):

> gs_pop_get_stats(gs)
             sample                                                                  pop count
 1: Activation1.fcs                                                                 root  2000
 2: Activation1.fcs                                                               /Cells  1458
 3: Activation1.fcs                                                  /Cells/Single Cells  1351
 4: Activation1.fcs                                       /Cells/Single Cells/Live Cells   619
 5: Activation1.fcs                       /Cells/Single Cells/Live Cells/Dendritic Cells    20
 6: Activation1.fcs                               /Cells/Single Cells/Live Cells/T Cells   218
 7: Activation1.fcs                   /Cells/Single Cells/Live Cells/T Cells/CD8 T Cells   117
 8: Activation1.fcs /Cells/Single Cells/Live Cells/T Cells/CD8 T Cells/CD69+ CD8 T Cells     0
 9: Activation1.fcs                   /Cells/Single Cells/Live Cells/T Cells/CD4 T Cells    85
10: Activation1.fcs /Cells/Single Cells/Live Cells/T Cells/CD4 T Cells/CD69+ CD4 T Cells     0
11: Activation2.fcs                                                                 root  2000
12: Activation2.fcs                                                               /Cells  1456
13: Activation2.fcs                                                  /Cells/Single Cells  1350
14: Activation2.fcs                                       /Cells/Single Cells/Live Cells   658
15: Activation2.fcs                       /Cells/Single Cells/Live Cells/Dendritic Cells    31
16: Activation2.fcs                               /Cells/Single Cells/Live Cells/T Cells   232
17: Activation2.fcs                   /Cells/Single Cells/Live Cells/T Cells/CD8 T Cells   116
18: Activation2.fcs /Cells/Single Cells/Live Cells/T Cells/CD8 T Cells/CD69+ CD8 T Cells     1
19: Activation2.fcs                   /Cells/Single Cells/Live Cells/T Cells/CD4 T Cells   105
20: Activation2.fcs /Cells/Single Cells/Live Cells/T Cells/CD4 T Cells/CD69+ CD4 T Cells     6
21: Activation3.fcs                                                                 root  2000
22: Activation3.fcs                                                               /Cells  1460
23: Activation3.fcs                                                  /Cells/Single Cells  1353
24: Activation3.fcs                                       /Cells/Single Cells/Live Cells   714
25: Activation3.fcs                       /Cells/Single Cells/Live Cells/Dendritic Cells    32
26: Activation3.fcs                               /Cells/Single Cells/Live Cells/T Cells   260
27: Activation3.fcs                   /Cells/Single Cells/Live Cells/T Cells/CD8 T Cells   157
28: Activation3.fcs /Cells/Single Cells/Live Cells/T Cells/CD8 T Cells/CD69+ CD8 T Cells     9
29: Activation3.fcs                   /Cells/Single Cells/Live Cells/T Cells/CD4 T Cells    89
30: Activation3.fcs /Cells/Single Cells/Live Cells/T Cells/CD4 T Cells/CD69+ CD4 T Cells     9
31: Activation4.fcs                                                                 root  2000
32: Activation4.fcs                                                               /Cells  1435
33: Activation4.fcs                                                  /Cells/Single Cells  1321
34: Activation4.fcs                                       /Cells/Single Cells/Live Cells   673
35: Activation4.fcs                       /Cells/Single Cells/Live Cells/Dendritic Cells    36
36: Activation4.fcs                               /Cells/Single Cells/Live Cells/T Cells   249
37: Activation4.fcs                   /Cells/Single Cells/Live Cells/T Cells/CD8 T Cells   140
38: Activation4.fcs /Cells/Single Cells/Live Cells/T Cells/CD8 T Cells/CD69+ CD8 T Cells    18
39: Activation4.fcs                   /Cells/Single Cells/Live Cells/T Cells/CD4 T Cells    71
40: Activation4.fcs /Cells/Single Cells/Live Cells/T Cells/CD4 T Cells/CD69+ CD4 T Cells    24
             sample                                                                  pop count

It's entirely possible that there's another separate issue, however. It's just hard to track it down without a reproducible example.

DillonHammill commented 4 years ago

I can send through the compensated and transformed GatingSet. I will see if I can generate a reprex without needing the gatingTemplate.

jacobpwagner commented 4 years ago

Well, if you can make a reprex for making the GatingSet (like in that snippet above), that's fine too.

DillonHammill commented 4 years ago

@jacobpwagner the raw data files come from CytoExploreRData which is a public repo. If you have trouble installing it, the individual files can be downloaded from https://github.com/DillonHammill/CytoExploreRData/tree/master/inst/extdata/Activation. Activation_33.fcs is an unstained control and I suspect this is where the problems are coming from. The initial gates are fine, it happens downstream...

# Load Packages ----------------------------------------------------------------
devtools::install_github("DillonHammill/CytoExploreRData")
library(CytoExploreRData)

# Data directory CytoExploreRData ----------------------------------------------
datadir <- system.file("extdata", package = "CytoExploreRData")

# Activation GatingSet ---------------------------------------------------------

# Assign Activation dataset to fs -
files <- list.files(paste0(datadir, "/Activation"), full.names = TRUE)
fs <- read.ncdfFlowSet(files)

# pData information -
pData(fs)$OVAConc <- c(rep(c(0, 0, 0.005, 0.005, 0.05, 0.05, 0.5, 0.5), 4), 0
                       )
pData(fs)$Treatment <- c(
  rep("Stim-A", 8),
  rep("Stim-B", 8),
  rep("Stim-C", 8),
  rep("Stim-D", 8),"NA"
)
pData(fs)$Treatment <- factor(pData(fs)$Treatment, levels = c(
  "Stim-A",
  "Stim-B",
  "Stim-C",
  "Stim-D", "NA"
))
chans <- c(
  "Alexa Fluor 405-A",
  "Alexa Fluor 430-A",
  "APC-Cy7-A", 
  "PE-A",
  "Alexa Fluor 488-A",
  "Alexa Fluor 700-A",
  "Alexa Fluor 647-A",
  "7-AAD-A"
)
markers <- c(
  "Hoechst-405",
  "Hoechst-430",
  "CD11c",
  "Va2",
  "CD8",
  "CD4",
  "CD44",
  "CD69"
)
names(markers) <- chans
markernames(fs) <- markers

# GatingSet -
gs <- GatingSet(fs)

# Compensation -
spill <- fs[[1]]@description$SPILL
gs <- compensate(gs, spill)

# Transformations -
chans <- c("Alexa Fluor 488-A", "PE-A", "7-AAD-A", "Alexa Fluor 405-A", "Alexa Fluor 430-A", "Alexa Fluor 647-A", "Alexa Fluor 700-A", "APC-Cy7-A")
trans <- estimateLogicle(gs[[1]], channels = chans)
gs <- transform(gs, trans)

# gatingTemplate -
gt <- gatingTemplate(paste0(datadir,"/Activation-gatingTemplate.csv"))

# gatingTemplate file -
gtf <- read.csv(system.file("extdata",  "Activation-gatingTemplate.csv", package = "CytoRSuiteData"))

# Gating -
gating(gt, gs)
DillonHammill commented 4 years ago

Hopefully the error pops up when you try to apply the CytoRSuite gatingTemplate as in the snippet you posted.

jacobpwagner commented 4 years ago

Other than adding the required library calls at the top:

library(flowCore)
library(flowWorkspace)
library(ncdfFlow)
library(openCyto)

I ran your snippet as it's written, but I get this:

> gating(gt, gs)
Loading required package: parallel
Error in .preprocessing(x, y, ...) : 
  Can't gate using unregistered method .pp_cyto_gate_draw

I thought that might be in CytoRSuite, so I added that library call, but still got the same message. So I'm guessing this is something you added to CytoExploreR, which I don't have access to. Can you give me read permissions for that repo?

DillonHammill commented 4 years ago

I have granted you access. You will need to apply the new gatingTemplate from CytoExploreRData. You will need to install the devel branch.

# Load Packages ----------------------------------------------------------------
devtools::install_github("DillonHammill/CytoExploreRData")
devtools::install_github("DillonHammill/CytoExploreR", ref = "devel")
library(CytoExploreRData)

# Data directory CytoExploreRData ----------------------------------------------
datadir <- system.file("extdata", package = "CytoExploreRData")

# Activation GatingSet ---------------------------------------------------------

# Assign Activation dataset to fs -
files <- list.files(paste0(datadir, "/Activation"), full.names = TRUE)
fs <- read.ncdfFlowSet(files)

# pData information -
pData(fs)$OVAConc <- c(rep(c(0, 0, 0.005, 0.005, 0.05, 0.05, 0.5, 0.5), 4), 0
                       )
pData(fs)$Treatment <- c(
  rep("Stim-A", 8),
  rep("Stim-B", 8),
  rep("Stim-C", 8),
  rep("Stim-D", 8),"NA"
)
pData(fs)$Treatment <- factor(pData(fs)$Treatment, levels = c(
  "Stim-A",
  "Stim-B",
  "Stim-C",
  "Stim-D", "NA"
))
chans <- c(
  "Alexa Fluor 405-A",
  "Alexa Fluor 430-A",
  "APC-Cy7-A", 
  "PE-A",
  "Alexa Fluor 488-A",
  "Alexa Fluor 700-A",
  "Alexa Fluor 647-A",
  "7-AAD-A"
)
markers <- c(
  "Hoechst-405",
  "Hoechst-430",
  "CD11c",
  "Va2",
  "CD8",
  "CD4",
  "CD44",
  "CD69"
)
names(markers) <- chans
markernames(fs) <- markers

# GatingSet -
gs <- GatingSet(fs)

# Compensation -
spill <- fs[[1]]@description$SPILL
gs <- compensate(gs, spill)

# Transformations -
chans <- c("Alexa Fluor 488-A", "PE-A", "7-AAD-A", "Alexa Fluor 405-A", "Alexa Fluor 430-A", "Alexa Fluor 647-A", "Alexa Fluor 700-A", "APC-Cy7-A")
trans <- estimateLogicle(gs[[1]], channels = chans)
gs <- transform(gs, trans)

# gatingTemplate -
gt <- gatingTemplate(paste0(datadir,"/Activation-gatingTemplate.csv"))

# Gating -
gt_gating(gt, gs)
jacobpwagner commented 4 years ago

Okay. I see it now. I'll get to work on it.

DillonHammill commented 4 years ago

Thanks @jacobpwagner.

jacobpwagner commented 4 years ago

After 1b8d9a36ce68f18cfc61f8d3a2f2de284f11a356, your snippet works fine for me and gt_gating succeeds. You may have to do a full rebuild of the relevant packages like before unfortunately. It looks like a few of the data driven methods complain about not having enough events from a few subpopulations, but that's a different problem.

2: In (function (fs, pp_res, gFunc, popAlias, channels, gFunc_args)  :
  Activation_33.fcs: Not enough events to proceed the data-driven gating!Returning a dummy gate instead.
3: In (function (fs, pp_res, gFunc, popAlias, channels, gFunc_args)  :
  Activation_33.fcs: Not enough events to proceed the data-driven gating!Returning a dummy gate instead.
4: In (function (fs, pp_res, gFunc, popAlias, channels, gFunc_args)  :
  Activation_33.fcs: Not enough events to proceed the data-driven gating!Returning a dummy gate instead.
5: In (function (fs, pp_res, gFunc, popAlias, channels, gFunc_args)  :
  Activation_33.fcs: Not enough events to proceed the data-driven gating!Returning a dummy gate instead.
DillonHammill commented 4 years ago

Thanks, these are the same warnings that I used to see. So as far I can tell it is now performing as well as it used to. I will try running on my computer now...

DillonHammill commented 4 years ago

Thanks @jacobpwagner. Really appreciate your help. This is working now as expected now. I will keep this open in case we want to address the warnings. From what I can tell they don't seem to impact on anything.

jacobpwagner commented 4 years ago

Yeah. I can look in to it, but they are probably appropriate warnings for automated gating methods applied to some of those minimal populations.

DillonHammill commented 4 years ago

I think it is because I am using the gatingTemplate in an unconventional way, by passing gates manually through the template. There aren't any data driven gates in these templates.

mikejiang commented 4 years ago

Not enough events warning is from openCyto, which is irrelevant to cytolib. So I am closing it for now.