--- a
+++ b/tests/testthat/test-ondisk.R
@@ -0,0 +1,262 @@
+# packages 
+skip_if_not_installed("rhdf5")
+skip_if_not_installed("pizzarr")
+skip_if_not_installed("HDF5Array")
+skip_if_not_installed("HDF5DataFrame")
+skip_if_not_installed("ZarrArray")
+skip_if_not_installed("ZarrDataFrame")
+skip_if_not_installed("ImageArray")
+skip_if_not_installed("BPCells")
+
+# create dir
+dir.create(td <- tempfile())
+output_zarr <- paste0(td, "/xenium_data_zarr_test")
+output_h5ad <- paste0(td, "/xenium_data_h5_test")
+output_h5ad_2 <- paste0(td, "/xenium_data_h5_test")
+output_h5ad_merged <- paste0(td, "/xenium_data_merged_h5_test")
+
+test_that("save and load on disk", {
+  
+  # get data
+  data("xenium_data")
+  
+  # HDF5
+  xenium_data2 <- saveVoltRon(xenium_data, 
+                              output = output_h5ad, 
+                              format = "HDF5VoltRon", 
+                              replace = TRUE, 
+                              verbose = FALSE)
+  xenium_data2 <- loadVoltRon(dir = output_h5ad)
+
+  # TODO: output zarr problem with path
+  # zarr
+  xenium_data2 <- saveVoltRon(xenium_data,
+                              output = output_zarr,
+                              format = "ZarrVoltRon",
+                              replace = TRUE,
+                              verbose = FALSE)
+  xenium_data2 <- loadVoltRon(dir = output_zarr)
+
+  # remove files
+  unlink(output_h5ad, recursive = TRUE)
+  unlink(output_zarr, recursive = TRUE)
+
+  expect_equal(1,1L)
+})
+
+test_that("double write and merging", {
+  
+  # get data
+  data("visium_data")
+  data("xenium_data")
+  
+  # TODO: output zarr problem with path
+  # write merged data
+  xenium_data2 <- xenium_data
+  xenium_data2$Sample <- "XeniumR2"
+  xenium_data2 <- merge(xenium_data, xenium_data2, verbose = FALSE)
+  xenium_data3 <- saveVoltRon(xenium_data2,
+                              output = output_zarr,
+                              format = "ZarrVoltRon",
+                              replace = TRUE,
+                              verbose = FALSE)
+  
+  # merged ondisk data
+  xenium_data_disk <- saveVoltRon(xenium_data, 
+                              output = output_h5ad, 
+                              format = "HDF5VoltRon", 
+                              replace = TRUE, 
+                              verbose = FALSE)
+  xenium_data2 <- xenium_data 
+  xenium_data2$Sample <- "XeniumR2"
+  xenium_data2_disk <- saveVoltRon(xenium_data2, 
+                                   output = output_h5ad_2, 
+                                   format = "HDF5VoltRon", 
+                                   replace = TRUE, 
+                                   verbose = FALSE)
+  xenium_data_merged <- merge(xenium_data_disk, xenium_data2_disk, verbose = FALSE)
+  xenium_data_merged_disk <- saveVoltRon(xenium_data_merged, 
+                                         output = output_h5ad_merged, 
+                                         format = "HDF5VoltRon", 
+                                         replace = TRUE, 
+                                         verbose = FALSE)
+  
+  # remove files
+  unlink(output_h5ad, recursive = TRUE)
+  unlink(output_h5ad_2, recursive = TRUE)
+  unlink(output_h5ad_merged, recursive = TRUE)
+  unlink(output_zarr, recursive = TRUE)
+  
+  expect_equal(1,1L)
+})
+
+test_that("subsetting", {
+  
+  # get data
+  data("xenium_data")
+  
+  # HDF5
+  xenium_data2 <- saveVoltRon(xenium_data, 
+                              output = output_h5ad, 
+                              format = "HDF5VoltRon", 
+                              replace = TRUE, 
+                              verbose = FALSE)
+
+  # by spatialpoints
+  spatpoints <- vrSpatialPoints(xenium_data2)[1:30]
+  xenium_data2_subset <- subset(xenium_data2, spatialpoints = spatpoints)
+  expect_equal(vrSpatialPoints(xenium_data2_subset), spatpoints)
+  
+  # visualize
+  vrSpatialPlot(xenium_data2_subset)
+  vrSpatialFeaturePlot(xenium_data2_subset, features = "KRT15")
+  
+  # by image
+  xenium_data2_subset <- subset(xenium_data2, image = "290x202+98+17")
+  expect_equal(length(vrSpatialPoints(xenium_data2_subset)), 392)
+  expect_equal(is(vrImages(xenium_data2_subset)), "magick-image")
+  expect_equal(is(vrImages(xenium_data2_subset, as.raster = TRUE)), "Image_Array")
+  
+  # visualize
+  vrSpatialPlot(xenium_data2_subset)
+  vrSpatialFeaturePlot(xenium_data2_subset, features = "KRT15")
+  
+  # remove files
+  unlink(output_h5ad, recursive = TRUE)
+  
+  expect_equal(1,1L)
+})
+
+test_that("metadata", {
+  
+  # get data
+  data("xenium_data")
+  
+  # HDF5
+  xenium_data2 <- saveVoltRon(xenium_data, 
+                              output = output_h5ad, 
+                              format = "HDF5VoltRon", 
+                              replace = TRUE,
+                              verbose = FALSE)
+  
+  # add column to metadata
+  xenium_data2$new_column <- vrSpatialPoints(xenium_data2)
+  
+  # save updated metadata
+  xenium_data2 <- saveVoltRon(xenium_data2, verbose = FALSE)
+  meta.data <- Metadata(xenium_data2)
+  expect_true(is(meta.data@listData[["new_column"]], "HDF5ColumnVector"))
+  
+  # load after update
+  xenium_data2 <- loadVoltRon(output_h5ad)
+  meta.data <- Metadata(xenium_data2)
+  expect_true(is(meta.data@listData[["new_column"]], "HDF5ColumnVector"))
+  
+  # remove files
+  unlink(output_h5ad, recursive = TRUE)
+  
+  expect_equal(1,1L)
+})
+
+
+test_that("embeddings", {
+  
+  # get data
+  data("xenium_data")
+  
+  # HDF5
+  xenium_data2 <- saveVoltRon(xenium_data, 
+                              output = output_h5ad, 
+                              format = "HDF5VoltRon", 
+                              replace = TRUE, 
+                              verbose = FALSE)
+  xenium_data2 <- getPCA(xenium_data2, features = vrFeatures(xenium_data2), type = "pca_key")
+  
+  # check embeddings
+  emb <- vrEmbeddings(xenium_data2, type = "pca_key")
+  expect_true(all(dim(emb) > 1))
+  expect_true(all(vrEmbeddingNames(xenium_data2) == c("pca", "umap", "pca_key")))
+  
+  # remove files
+  unlink(output_h5ad, recursive = TRUE)
+  
+  expect_equal(1,1L)
+})
+
+test_that("visualization", {
+  
+  # get data
+  data("visium_data")
+  
+  # HDF5
+  visium_data2 <- saveVoltRon(visium_data, 
+                              output = output_h5ad, 
+                              format = "HDF5VoltRon", 
+                              replace = TRUE, 
+                              verbose = FALSE)
+
+  # check spatial plots
+  vrSpatialPlot(visium_data2, group.by = "Sample")
+  vrSpatialFeaturePlot(visium_data2, features = "Count")
+  
+  # check multilayer spatial plots
+  data("merged_object")
+  merged_object2 <- saveVoltRon(merged_object, 
+                                output = output_h5ad, 
+                                format = "HDF5VoltRon", 
+                                replace = TRUE, 
+                                verbose = FALSE)
+  vrSpatialPlot(merged_object2, plot.segments = FALSE, n.tile = 100) |>
+    addSpatialLayer(merged_object2, assay = "Assay3", group.by = "Sample", alpha = 0.4, colors = list(Block = "blue"))
+  vrSpatialPlot(merged_object2, plot.segments = TRUE, n.tile = 100) |>
+    addSpatialLayer(merged_object2, assay = "Assay3", group.by = "Sample", alpha = 0.4, colors = list(Block = "blue"))
+  vrSpatialPlot(merged_object2, assay = "Assay2", group.by = "gene", alpha = 1, colors = list(KRT15 = "blue", KRT14 = "green"), n.tile = 100) |>
+    addSpatialLayer(merged_object2, assay = "Assay3", group.by = "Sample", alpha = 0.4, colors = list(Block = "blue"))
+  
+  # remove files
+  unlink(output_h5ad, recursive = TRUE)
+  
+  expect_equal(1,1L)
+})
+
+test_that("neighbors", {
+  
+  # get data
+  data("xenium_data")
+  
+  # HDF5
+  xenium_data2 <- saveVoltRon(xenium_data, 
+                              output = output_h5ad, 
+                              format = "HDF5VoltRon", 
+                              replace = TRUE, 
+                              verbose = FALSE)
+  
+  # spatial neighbors, delaunay
+  xenium_data2 <- getSpatialNeighbors(xenium_data2, method = "delaunay")
+  graphs <- vrGraph(xenium_data2, graph.type = "delaunay")
+  expect_true(inherits(graphs,"igraph"))
+  expect_true(length(igraph::E(graphs)) > 0)
+  
+  # profile neighbors, kNN
+  xenium_data2 <- getProfileNeighbors(xenium_data2, method = "kNN", k = 10, data.type = "pca")
+  graphs <- vrGraph(xenium_data2, graph.type = "kNN")
+  expect_true(inherits(graphs,"igraph"))
+  expect_true(length(igraph::E(graphs)) > 0)
+  xenium_data2 <- getClusters(xenium_data2, graph = "kNN", label = "cluster_knn")
+  expect_true(is.numeric(unique(xenium_data2$cluster_knn)))
+  
+  # update metadata
+  xenium_data2 <- saveVoltRon(xenium_data2, verbose = FALSE)
+  meta.data <- Metadata(xenium_data2)
+  expect_true(is(meta.data@listData[["cluster_knn"]], "HDF5ColumnVector"))
+  
+  # load after update
+  xenium_data2 <- loadVoltRon(output_h5ad)
+  meta.data <- Metadata(xenium_data2)
+  expect_true(is(meta.data@listData[["cluster_knn"]], "HDF5ColumnVector"))
+  
+  # remove files
+  unlink(output_h5ad, recursive = TRUE)
+  
+  expect_equal(1,1L)
+})
\ No newline at end of file