Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Ellipsoid Gate: Appearance != Performance #268

Open
gunthergl opened this issue Mar 13, 2024 · 1 comment
Open

Ellipsoid Gate: Appearance != Performance #268

gunthergl opened this issue Mar 13, 2024 · 1 comment

Comments

@gunthergl
Copy link

A particular ellipsoid gate gives completely off counts.

To reproduce,

  • I create a grid of single cells
  • gate_ellipsoid is visually as intended but results in wrong frequencies/counts (31%)
  • gate_ellipsoid_2 and gate_rectangle seem to work fine (1.9%)

Expected behavior
gate_ellipsoid should have a similar percentage of cells as the other two gates.

Additional info

  • I checked that the covariance matrices are positive definite with eigen().
  • gate_ellipsoid_2 is similar, but not identical, so that is no workaround

Any idea what is happening?

example_expression <- expand.grid(
    x = seq(1, 1e4, length.out = 1e4 / 150),
    y = seq(1, 2e3, length.out = 2e3 / 20)
)
ff <- flowCore::flowFrame(exprs = as.matrix(example_expression))

my_covmat <- matrix(c(774761, -309675, -12387, 16687), nrow = 2)
colnames(my_covmat) <- c("x", "y")
rownames(my_covmat) <- c("x", "y")
eigen(my_covmat)
#> eigen() decomposition
#> $values
#> [1] 779787.79  11660.21
#> 
#> $vectors
#>            [,1]       [,2]
#> [1,]  0.9266082 0.01623032
#> [2,] -0.3760282 0.99986828


my_covmat_2 <- matrix(c(774761, -100000, -10000, 16687), nrow = 2)
colnames(my_covmat_2) <- c("x", "y")
rownames(my_covmat_2) <- c("x", "y")
eigen(my_covmat_2)
#> eigen() decomposition
#> $values
#> [1] 776077.84  15370.16
#> 
#> $vectors
#>            [,1]       [,2]
#> [1,]  0.9914408 0.01316731
#> [2,] -0.1305574 0.99991331

gate_ellipsoid <- flowCore::ellipsoidGate(mean = c(5000, 1000), cov = my_covmat)
gate_ellipsoid_2 <- flowCore::ellipsoidGate(mean = c(5000, 1000), cov = my_covmat_2)
gate_rectangle <- flowCore::rectangleGate("x" = c(4000, 6000), "y" = c(900, 1100))
plotlist <- lapply(list(gate_ellipsoid, gate_ellipsoid_2, gate_rectangle), function(gate_x) {
    flowViz::xyplot(
        y ~ x,
        data = ff,
        nrpoints = 10000,
        smooth = FALSE,
        stats = TRUE,
        filter = gate_x
    )
})
print(plotlist)
#> [[1]]

#> 
#> [[2]]

#> 
#> [[3]]

sessionInfo()
#> R version 4.3.3 (2024-02-29)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04.6 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3;  LAPACK version 3.9.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Europe/Berlin
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> loaded via a namespace (and not attached):
#>  [1] compiler_4.3.3      reprex_2.1.0.9000   Rcpp_1.0.12        
#>  [4] Biobase_2.60.0      cytolib_2.12.0      png_0.1-8          
#>  [7] yaml_2.3.8          fastmap_1.1.1       lattice_0.22-5     
#> [10] deldir_2.0-2        hexbin_1.28.3       RProtoBufLib_2.12.0
#> [13] latticeExtra_0.6-30 knitr_1.45          MASS_7.3-60        
#> [16] BiocGenerics_0.46.0 tibble_3.2.1        interp_1.1-6       
#> [19] R.cache_0.15.0      RColorBrewer_1.1-3  pillar_1.9.0       
#> [22] R.utils_2.11.0      rlang_1.1.3         utf8_1.2.4         
#> [25] flowCore_2.12.0     xfun_0.42           fs_1.6.3           
#> [28] IDPmisc_1.1.21      cli_3.6.2           withr_3.0.0        
#> [31] magrittr_2.0.3      digest_0.6.34       grid_4.3.3         
#> [34] lifecycle_1.0.4     R.methodsS3_1.8.1   R.oo_1.24.0        
#> [37] S4Vectors_0.38.1    vctrs_0.6.5         KernSmooth_2.23-22 
#> [40] evaluate_0.23       glue_1.7.0          flowViz_1.64.0     
#> [43] styler_1.7.0        stats4_4.3.3        fansi_1.0.6        
#> [46] rmarkdown_2.26      purrr_1.0.2         jpeg_0.1-10        
#> [49] matrixStats_1.2.0   tools_4.3.3         pkgconfig_2.0.3    
#> [52] htmltools_0.5.7

Created on 2024-03-13 with reprex v2.1.0.9000

@gunthergl
Copy link
Author

gunthergl commented Mar 15, 2024

Summary:
Non-symmetrical matrices defining the ellipse do not make sense (obviously).
Everything is fine within flowCore. Sorry for the confusion.

Potential enhancement:
Raise an error if the covariance matrix given to ellipsoidGate is not symmetric + positive definite. I will leave this issue open in case you want to have this error added. I could also propose a pull request - but I am unsure if that's something you want to fit in the package.

Long:
From previous (wrong) calculations, I came to the described matrices. I didn't think about them being non-symmetric after I "checked" that they were positive definite according to the Eigenvalues (which is wrong because it's non-Hermitian matrices).

Given a matrix which is not positive definite, ellipsoidGate does something funny, shown in the next plots.

Interestingly:

  1. The gate looks as it should.
  2. Calculating a polygon from the (wrongly established!) ellipsoidGate works.
example_expression <- expand.grid(
    x = seq(1, 1e4, length.out = 1e4 / 150),
    y = seq(1, 2e3, length.out = 2e3 / 20)
)
ff <- flowCore::flowFrame(exprs = as.matrix(example_expression))

my_covmat <- matrix(c(774761, -309675, -12387, 16687), nrow = 2)
colnames(my_covmat) <- c("x", "y")
rownames(my_covmat) <- c("x", "y")
gate_ellipsoid <- flowCore::ellipsoidGate(mean = c(5000, 1000), cov = my_covmat, filterId = "gate_ellipsoid")
gate_ellipsoid_pg <- as(gate_ellipsoid, "polygonGate")
gate_ellipsoid_pg@filterId <- "gate_ellipsoid_pg"

all_gates <- list(
    gate_ellipsoid,
    gate_ellipsoid_pg
)

flowCore::write.FCS(ff, "removeme.fcs")
#> [1] "removeme.fcs"
cs <- flowWorkspace::load_cytoset_from_fcs("removeme.fcs")
gs <- flowWorkspace::GatingSet(cs)
for (gate_x in all_gates) {
    flowWorkspace::gs_pop_add(gs, gate_x)
}
gated <- flowWorkspace::gh_apply_to_cs(gs, cs)
#> generating new GatingSet from the gating template...
flowWorkspace::gs_get_pop_paths(gated)
#> [1] "root"               "/gate_ellipsoid"    "/gate_ellipsoid_pg"
flowWorkspace::recompute(gated)
#> done!
print(flowWorkspace::gs_pop_get_count_fast(gated))
#>            name         Population Parent Count ParentCount
#>          <char>             <char> <char> <int>       <int>
#> 1: removeme.fcs    /gate_ellipsoid   root  2052        6700
#> 2: removeme.fcs /gate_ellipsoid_pg   root   114        6700
pacman::p_load(ggplot2)
# pdf("removeme.pdf")
for (gate_x in all_gates) {
    print(
        ggcyto::ggcyto(gated, ggplot2::aes(x = x, y = y)) +
            # ggplot2::geom_point(size = .25, aes(col = flowWorkspace::gh_pop_get_indices(gated, gate_x@filterId))) +
            ggplot2::geom_point(aes(col = flowWorkspace::gh_pop_get_indices(gated, gate_x@filterId))) +
            ggcyto::ggcyto_par_set(limits = list(
                x = c(4000, 6000),
                y = c(800, 1200)
            )) +
            ggcyto::geom_gate(gate_x@filterId) +
            ggplot2::ggtitle(gate_x@filterId) +
            ggpubr::theme_pubr()
    )
    # print(ggcyto::autoplot(gated, gate_x@filterId))
    # ggcyto::autoplot(gated, "/gate_rectangle")
}
#> Coordinate system already present. Adding new coordinate system, which will
#> replace the existing one.
#> Coordinate system already present. Adding new coordinate system, which will
#> replace the existing one.

sessionInfo()
#> R version 4.3.3 (2024-02-29)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04.6 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3;  LAPACK version 3.9.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Europe/Berlin
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_3.5.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.4         xfun_0.42            rstatix_0.7.2       
#>  [4] Biobase_2.60.0       lattice_0.22-5       vctrs_0.6.5         
#>  [7] tools_4.3.3          generics_0.1.3       stats4_4.3.3        
#> [10] tibble_3.2.1         flowWorkspace_4.12.0 fansi_1.0.6         
#> [13] pacman_0.5.1         pkgconfig_2.0.3      R.oo_1.24.0         
#> [16] data.table_1.15.2    RColorBrewer_1.1-3   S4Vectors_0.38.1    
#> [19] graph_1.78.0         lifecycle_1.0.4      R.cache_0.15.0      
#> [22] compiler_4.3.3       farver_2.1.1         munsell_0.5.0       
#> [25] carData_3.0-5        htmltools_0.5.7      yaml_2.3.8          
#> [28] flowCore_2.12.0      pillar_1.9.0         hexbin_1.28.3       
#> [31] car_3.1-2            ggpubr_0.6.0         tidyr_1.3.1         
#> [34] R.utils_2.11.0       abind_1.4-5          RProtoBufLib_2.12.0 
#> [37] styler_1.7.0         tidyselect_1.2.0     digest_0.6.35       
#> [40] dplyr_1.1.4          purrr_1.0.2          labeling_0.4.3      
#> [43] fastmap_1.1.1        grid_4.3.3           colorspace_2.1-0    
#> [46] cli_3.6.2            magrittr_2.0.3       ncdfFlow_2.46.0     
#> [49] XML_3.99-0.16.1      utf8_1.2.4           broom_1.0.5         
#> [52] withr_3.0.0          scales_1.3.0         backports_1.4.1     
#> [55] rmarkdown_2.26       matrixStats_1.2.0    gridExtra_2.3       
#> [58] ggsignif_0.6.4       cytolib_2.12.0       R.methodsS3_1.8.1   
#> [61] evaluate_0.23        knitr_1.45           rlang_1.1.3         
#> [64] Rcpp_1.0.12          glue_1.7.0           Rgraphviz_2.44.0    
#> [67] BiocGenerics_0.46.0  reprex_2.1.0.9000    R6_2.5.1            
#> [70] plyr_1.8.9           fs_1.6.3             zlibbioc_1.46.0     
#> [73] ggcyto_1.30.2

Created on 2024-03-15 with reprex v2.1.0.9000

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant