Ellipsoid Gate: Appearance != Performance
Opened this issue · 1 comments
gunthergl commented
A particular ellipsoid gate gives completely off counts.
To reproduce,
- I create a grid of single cells
gate_ellipsoidis visually as intended but results in wrong frequencies/counts (31%)gate_ellipsoid_2andgate_rectangleseem 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_2is 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.7Created on 2024-03-13 with reprex v2.1.0.9000
gunthergl commented
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:
- The gate looks as it should.
- Calculating a polygon from the (wrongly established!)
ellipsoidGateworks.
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.2Created on 2024-03-15 with reprex v2.1.0.9000