diff --git a/NAMESPACE b/NAMESPACE
index 9a5a9987..ad4a13cb 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -160,6 +160,7 @@ export(importCosMx)
 export(importGenePS)
 export(importGeoMx)
 export(importImageData)
+export(importPhenoCycler)
 export(importSTOmics)
 export(importVisium)
 export(importXenium)
diff --git a/R/import.R b/R/import.R
index 6ed2b913..8ece3348 100644
--- a/R/import.R
+++ b/R/import.R
@@ -1358,6 +1358,318 @@ importSTOmics <- function(h5ad.path, assay_name = "STOmics", sample_name = NULL,
   formVoltRon(rawdata, metadata = metadata, coords = coords, main.assay = assay_name, params = params, assay.type = "spot", image_name = image_name, main_channel = channel_name, sample_name = sample_name, ...)
 }
 
+####
+# Akoya ####
+####
+
+####
+## PhenoCycler ####
+####
+
+#' readPhenoCyclerMat
+#' 
+#' Read and Load Akoya CODEX data, adapted from \code{ReadAkoya} function \code{Seurat} package
+#'
+#' @param filename Path to matrix generated by upstream processing.
+#' @param type Specify which type matrix is being provided.
+#' \itemize{
+#'  \item \dQuote{\code{processor}}: matrix generated by CODEX Processor
+#'  \item \dQuote{\code{inform}}: matrix generated by inForm
+#'  \item \dQuote{\code{qupath}}: matrix generated by QuPath
+#' }
+#' @param filter A pattern to filter features by; pass \code{NA} to skip feature filtering
+#' @param inform.quant When \code{type} is \dQuote{\code{inform}}, the quantification level to read in
+#' 
+#' @noRd
+readPhenoCyclerMat <- function(
+    filename,
+    type = c('inform', 'processor', 'qupath'),
+    filter = 'DAPI|Blank|Empty',
+    inform.quant = c('mean', 'total', 'min', 'max', 'std')
+) {
+  if (!requireNamespace("data.table", quietly = TRUE)) {
+    stop("Please install 'data.table' for this function")
+  }
+  # Check arguments
+  if (!file.exists(filename)) {
+    stop(paste("Can't file file:", filename))
+  }
+  type <- tolower(x = type[1L])
+  type <- match.arg(arg = type)
+  ratio <- getOption(x = 'Seurat.input.sparse_ratio', default = 0.4)
+  # Preload matrix
+  sep <- switch(EXPR = type, 'inform' = '\t', ',')
+  mtx <- data.table::fread(
+    file = filename,
+    sep = sep,
+    data.table = FALSE,
+    verbose = FALSE
+  )
+  # Assemble outputs
+  outs <- switch(
+    EXPR = type,
+    'processor' = {
+      # Create centroids data frame
+      centroids <- data.frame(
+        x = mtx[['x:x']],
+        y = mtx[['y:y']],
+        cell = as.character(x = mtx[['cell_id:cell_id']]),
+        stringsAsFactors = FALSE
+      )
+      rownames(x = mtx) <- as.character(x = mtx[['cell_id:cell_id']])
+      # Create metadata data frame
+      md <- mtx[, !grepl(pattern = '^cyc', x = colnames(x = mtx)), drop = FALSE]
+      colnames(x = md) <- vapply(
+        X = strsplit(x = colnames(x = md), split = ':'),
+        FUN = '[[',
+        FUN.VALUE = character(length = 1L),
+        2L
+      )
+      # Create expression matrix
+      mtx <- mtx[, grepl(pattern = '^cyc', x = colnames(x = mtx)), drop = FALSE]
+      colnames(x = mtx) <- vapply(
+        X = strsplit(x = colnames(x = mtx), split = ':'),
+        FUN = '[[',
+        FUN.VALUE = character(length = 1L),
+        2L
+      )
+      if (!is.na(x = filter)) {
+        mtx <- mtx[, !grepl(pattern = filter, x = colnames(x = mtx)), drop = FALSE]
+      }
+      mtx <- t(x = mtx)
+      if ((sum(mtx == 0) / length(x = mtx)) > ratio) {
+        mtx <- as.sparse(x = mtx)
+      }
+      list(matrix = mtx, centroids = centroids, metadata = md)
+    },
+    'inform' = {
+      inform.quant <- tolower(x = inform.quant[1L])
+      inform.quant <- match.arg(arg = inform.quant)
+      expr.key <- c(
+        mean = 'Mean',
+        total = 'Total',
+        min = 'Min',
+        max = 'Max',
+        std = 'Std Dev'
+      )[inform.quant]
+      expr.pattern <- '\\(Normalized Counts, Total Weighting\\)'
+      rownames(x = mtx) <- mtx[['Cell ID']]
+      mtx <- mtx[, setdiff(x = colnames(x = mtx), y = 'Cell ID'), drop = FALSE]
+      # Create centroids
+      centroids <- data.frame(
+        x = mtx[['Cell X Position']],
+        y = mtx[['Cell Y Position']],
+        cell  = rownames(x = mtx),
+        stringsAsFactors = FALSE
+      )
+      # Create metadata
+      cols <- setdiff(
+        x = grep(
+          pattern = expr.pattern,
+          x = colnames(x = mtx),
+          value = TRUE,
+          invert = TRUE
+        ),
+        y = paste('Cell', c('X', 'Y'), 'Position')
+      )
+      md <- mtx[, cols, drop = FALSE]
+      # Create expression matrices
+      exprs <- data.frame(
+        cols = grep(
+          pattern = paste(expr.key, expr.pattern),
+          x = colnames(x = mtx),
+          value = TRUE
+        )
+      )
+      exprs$feature <- vapply(
+        X = trimws(x = gsub(
+          pattern = paste(expr.key, expr.pattern),
+          replacement = '',
+          x = exprs$cols
+        )),
+        FUN = function(x) {
+          x <- unlist(x = strsplit(x = x, split = ' '))
+          x <- x[length(x = x)]
+          return(gsub(pattern = '\\(|\\)', replacement = '', x = x))
+        },
+        FUN.VALUE = character(length = 1L)
+      )
+      exprs$class <- tolower(x = vapply(
+        X = strsplit(x = exprs$cols, split = ' '),
+        FUN = '[[',
+        FUN.VALUE = character(length = 1L),
+        1L
+      ))
+      classes <- unique(x = exprs$class)
+      outs <- vector(
+        mode = 'list',
+        length = length(x = classes) + 2L
+      )
+      names(x = outs) <- c(
+        'matrix',
+        'centroids',
+        'metadata',
+        setdiff(x = classes, y = 'entire')
+      )
+      outs$centroids <- centroids
+      outs$metadata <- md
+      # browser()
+      for (i in classes) {
+        df <- exprs[exprs$class == i, , drop = FALSE]
+        expr <- mtx[, df$cols]
+        colnames(x = expr) <- df$feature
+        if (!is.na(x = filter)) {
+          expr <- expr[, !grepl(pattern = filter, x = colnames(x = expr)), drop = FALSE]
+        }
+        expr <- t(x = expr)
+        if ((sum(expr == 0, na.rm = TRUE) / length(x = expr)) > ratio) {
+          expr <- as.sparse(x = expr)
+        }
+        outs[[switch(EXPR = i, 'entire' = 'matrix', i)]] <- expr
+      }
+      outs
+    },
+    'qupath' = {
+      rownames(x = mtx) <- as.character(x = seq_len(length.out = nrow(x = mtx)))
+      # Create centroids
+      xpos <- sort(
+        x = grep(pattern = 'Centroid X', x = colnames(x = mtx), value = TRUE),
+        decreasing = TRUE
+      )[1L]
+      ypos <- sort(
+        x = grep(pattern = 'Centroid Y', x = colnames(x = mtx), value = TRUE),
+        decreasing = TRUE
+      )[1L]
+      centroids <- data.frame(
+        x = mtx[[xpos]],
+        y = mtx[[ypos]],
+        cell = rownames(x = mtx),
+        stringsAsFactors = FALSE
+      )
+      # Create metadata
+      cols <- setdiff(
+        x = grep(
+          pattern = 'Cell: Mean',
+          x = colnames(x = mtx),
+          ignore.case = TRUE,
+          value = TRUE,
+          invert = TRUE
+        ),
+        y = c(xpos, ypos)
+      )
+      md <- mtx[, cols, drop = FALSE]
+      # Create expression matrix
+      idx <- which(x = grepl(
+        pattern = 'Cell: Mean',
+        x = colnames(x = mtx),
+        ignore.case = TRUE
+      ))
+      mtx <- mtx[, idx, drop = FALSE]
+      colnames(x = mtx) <- vapply(
+        X = strsplit(x = colnames(x = mtx), split = ':'),
+        FUN = '[[',
+        FUN.VALUE = character(length = 1L),
+        1L
+      )
+      if (!is.na(x = filter)) {
+        mtx <- mtx[, !grepl(pattern = filter, x = colnames(x = mtx)), drop = FALSE]
+      }
+      mtx <- t(x = mtx)
+      if ((sum(mtx == 0) / length(x = mtx)) > ratio) {
+        mtx <- as.sparse(x = mtx)
+      }
+      list(matrix = mtx, centroids = centroids, metadata = md)
+    },
+    stop("Unknown matrix type: ", type)
+  )
+  return(outs)
+}
+
+
+#' Title
+#'
+#' @param dir.path path to PhenoCycler output folder
+#' @param assay_name the assay name of the SR object
+#' @param sample_name the name of the sample
+#' @param image_name the image name of the Xenium assay, Default: main
+#' @param filename Path to matrix generated by upstream processing.
+#' @param type Specify which type matrix is being provided.
+#' \itemize{
+#'  \item \dQuote{\code{processor}}: matrix generated by CODEX Processor
+#'  \item \dQuote{\code{inform}}: matrix generated by inForm
+#'  \item \dQuote{\code{qupath}}: matrix generated by QuPath
+#' }
+#' @param filter A pattern to filter features by; pass \code{NA} to skip feature filtering
+#' @param inform.quant When \code{type} is \dQuote{\code{inform}}, the quantification level to read in 
+#' @param ... additional parameters passed to \link{formVoltRon} 
+#' 
+#' @importFrom magick image_info image_read
+#'
+#' @export
+importPhenoCycler <- function(dir.path, assay_name = "PhenoCycler", sample_name = NULL, image_name = "main", 
+                              type = c('inform', 'processor', 'qupath'), filter = 'DAPI|Blank|Empty', inform.quant = c('mean', 'total', 'min', 'max', 'std'), ...){
+  
+  # raw counts, metadata and coordinates 
+  listoffiles <- list.files(paste0(dir.path, "/processed/segm/segm-1/fcs/compensated/"), full.names = TRUE)
+  datafile <- listoffiles[grepl("_compensated.csv$", listoffiles)][1]
+  if(file.exists(datafile)){
+    rawdata <- readPhenoCyclerMat(filename = datafile, type = type, filter = filter, inform.quant = inform.quant)
+  } else {
+    stop("There are no files named ending with '_compensated.csv' in the processed/segm/segm-1/fcs/compensated/ subfolder")
+  }
+  
+  # cell id 
+  cellid <- paste0("cell", colnames(rawdata$matrix))
+  
+  # coordinates
+  coords <- rawdata$centroids[,c("x", "y")]
+  rownames(coords) <- cellid
+  coords <- as.matrix(coords)
+  
+  # metadata
+  metadata <- rawdata$metadata
+  rownames(metadata) <- cellid
+  
+  # data
+  rawdata <- rawdata$matrix
+  colnames(rawdata) <- cellid
+  rownames(rawdata) <- gsub("(\t|\r|\n)", "", rownames(rawdata))
+  
+  # images 
+  image_dir <- paste0(dir.path, "/processed/stitched/reg001/")
+  list_files <- list.files(image_dir)
+  if(!dir.exists(image_dir)){
+    message("There are no images of channels!")
+    image_list <- NULL
+  } else {
+    if(!any(grepl(".tif$", list_files))){
+      stop("The folder doesnt have any images associated with channels!")
+    } else{
+      image_channel_names <- sapply(list_files, function(x) {
+        name <- strsplit(x, split = "_")[[1]]
+        name <- gsub(".tif$", "", name[length(name)])
+        name <- make.unique(name)
+        return(name)
+      })
+      image_channel_names <- c(image_channel_names[grepl("DAPI", image_channel_names)][1], 
+                               image_channel_names[image_channel_names %in% rownames(rawdata)])
+      list_files <- names(image_channel_names)
+      image_list <- lapply(list_files, function(x){
+        magick::image_read(paste0(image_dir, "/", x))
+      })
+      names(image_list) <- image_channel_names   
+      coords[,2] <- magick::image_info(image_list[[1]])$height - coords[,2]
+    }
+  }
+  
+  # voltron object
+  object <- formVoltRon(data = rawdata, metadata = metadata, image = image_list, coords = coords, assay.type = "cell", 
+                        sample_name = sample_name, main.assay = assay_name, image_name = image_name, ...)
+  
+  # return
+  object
+}
+
 ####
 # Image Data ####
 ####
diff --git a/R/visualization.R b/R/visualization.R
index 241403c2..c7666b22 100644
--- a/R/visualization.R
+++ b/R/visualization.R
@@ -1253,9 +1253,11 @@ vrEmbeddingFeaturePlot <- function(object, embedding = "pca", features = NULL, n
 
     # get data
     if(feat %in% vrFeatures(object)){
-      datax$score <- normdata[feat,]
+      # datax$score <- normdata[feat,]
+      datax$score <- normdata[feat, rownames(datax)]
     } else {
-      datax$score <- metadata[,feat]
+      # datax$score <- metadata[,feat]
+      datax$score <- metadata[rownames(datax),feat]
     }
 
     # get image information and plotting features
diff --git a/man/importPhenoCycler.Rd b/man/importPhenoCycler.Rd
new file mode 100644
index 00000000..bdf1bcbb
--- /dev/null
+++ b/man/importPhenoCycler.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/import.R
+\name{importPhenoCycler}
+\alias{importPhenoCycler}
+\title{Title}
+\usage{
+importPhenoCycler(
+  dir.path,
+  assay_name = "PhenoCycler",
+  sample_name = NULL,
+  image_name = "main",
+  type = c("inform", "processor", "qupath"),
+  filter = "DAPI|Blank|Empty",
+  inform.quant = c("mean", "total", "min", "max", "std"),
+  ...
+)
+}
+\arguments{
+\item{dir.path}{path to PhenoCycler output folder}
+
+\item{assay_name}{the assay name of the SR object}
+
+\item{sample_name}{the name of the sample}
+
+\item{image_name}{the image name of the Xenium assay, Default: main}
+
+\item{type}{Specify which type matrix is being provided.
+\itemize{
+ \item \dQuote{\code{processor}}: matrix generated by CODEX Processor
+ \item \dQuote{\code{inform}}: matrix generated by inForm
+ \item \dQuote{\code{qupath}}: matrix generated by QuPath
+}}
+
+\item{filter}{A pattern to filter features by; pass \code{NA} to skip feature filtering}
+
+\item{inform.quant}{When \code{type} is \dQuote{\code{inform}}, the quantification level to read in}
+
+\item{...}{additional parameters passed to \link{formVoltRon}}
+
+\item{filename}{Path to matrix generated by upstream processing.}
+}
+\description{
+Title
+}
diff --git a/tests/testthat/test-coordinates.R b/tests/testthat/test-coordinates.R
new file mode 100644
index 00000000..6ff0169e
--- /dev/null
+++ b/tests/testthat/test-coordinates.R
@@ -0,0 +1,27 @@
+# Testing functions of manipulating coordinates ####
+test_that("coordinates", {
+  
+  # get data
+  data("visium_data")
+  
+  # coordinates
+  coords <- vrCoordinates(visium_data)
+  coords <- vrCoordinates(visium_data, image_name = "main")
+  coords <- vrCoordinates(visium_data, spatial_name = "main")
+  expect_warning(coords <- vrCoordinates(visium_data, reg = TRUE))
+  expect_warning(coords <- vrCoordinates(visium_data, assay = "Assay1", reg = TRUE))
+  
+  # update coordinates
+  vrCoordinates(visium_data) <- coords*2
+  expect_error(vrCoordinates(visium_data, reg = TRUE) <- coords*3)
+  
+  # flip coordinates
+  visium_data <- flipCoordinates(visium_data)
+  
+  # segments
+  segments <- vrSegments(visium_data)
+  expect_warning(segments <- vrSegments(visium_data, reg = TRUE))
+  expect_warning(segments <- vrSegments(visium_data, assay = "Assay1", reg = TRUE))
+  
+  expect_equal(1,1L)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-embeddings.R b/tests/testthat/test-embeddings.R
new file mode 100644
index 00000000..6947230e
--- /dev/null
+++ b/tests/testthat/test-embeddings.R
@@ -0,0 +1,16 @@
+# Testing functions of manipulating embeddings ####
+test_that("embeddings", {
+  
+  # get data
+  data("visium_data")
+  
+  # write embedding
+  vrEmbeddings(visium_data, type = "pca") <- vrCoordinates(visium_data)
+  
+  # check overwrite embeddings
+  expect_error(vrEmbeddings(visium_data, type = "pca") <- vrCoordinates(visium_data))
+  vrEmbeddings(visium_data, type = "pca", overwrite = TRUE) <- vrCoordinates(visium_data)
+  
+  # return
+  expect_equal(1,1L)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-image.R b/tests/testthat/test-image.R
new file mode 100644
index 00000000..108e446c
--- /dev/null
+++ b/tests/testthat/test-image.R
@@ -0,0 +1,36 @@
+# Testing functions of manipulating images ####
+test_that("image", {
+  
+  # get data
+  data("visium_data")
+  
+  # get image
+  images <- vrImages(visium_data)
+  images <- vrImages(visium_data, name = "main")
+  expect_error(images <- vrImages(visium_data, name = "main2"))
+  images <- vrImages(visium_data, name = "main", channel = "H&E")
+  expect_warning(images <- vrImages(visium_data, name = "main", channel = "H&E2"))
+  
+  # manipulate image
+  visium_data_resize <- resizeImage(visium_data, size = 400)
+  visium_data_modulate <- modulateImage(visium_data, brightness = 400)
+  
+  # add new image
+  visium_data[["Assay1"]]@image[["new_image"]] <- vrImages(visium_data_resize)
+  
+  # get image names
+  expect_equal(vrImageNames(visium_data), c("main", "new_image"))
+  
+  # get main image
+  expect_equal(vrMainImage(visium_data[["Assay1"]]), "main")
+  expect_equal(vrMainSpatial(visium_data[["Assay1"]]), "main")
+  
+  # change main image
+  vrMainImage(visium_data[["Assay1"]]) <- "new_image"
+  vrMainSpatial(visium_data[["Assay1"]]) <- "new_image"
+  expect_equal(vrMainImage(visium_data[["Assay1"]]), "new_image")
+  expect_equal(vrMainSpatial(visium_data[["Assay1"]]), "new_image")
+  
+  # return
+  expect_equal(1,1L)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-imagedata.R b/tests/testthat/test-imagedata.R
new file mode 100644
index 00000000..058ecd5f
--- /dev/null
+++ b/tests/testthat/test-imagedata.R
@@ -0,0 +1,18 @@
+# Testing functions of image datasets ####
+test_that("imagedata", {
+  
+  # import image data
+  img <- system.file(package = "VoltRon", "extdata/DAPI.tif")
+  vrimgdata <- importImageData(img, tile.size = 10, stack.id = 1)
+  vrimgdata <- importImageData(img, tile.size = 10, stack.id = 2)
+  
+  # import image data when image is a stack
+  img_magick <- magick::image_read(img)
+  img_stack <- c(img_magick, img_magick)
+  vrimgdata <- importImageData(img_stack, tile.size = 10, stack.id = 1)
+  vrimgdata <- importImageData(img_stack, tile.size = 10, stack.id = 2)
+  expect_error(vrimgdata <- importImageData(img_stack, tile.size = 10, stack.id = 3))
+  
+  # return
+  expect_equal(1,1L)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
new file mode 100644
index 00000000..f3f7b484
--- /dev/null
+++ b/tests/testthat/test-plots.R
@@ -0,0 +1,68 @@
+# Testing plotting functions
+test_that("plots", {
+
+  # get data
+  data("visium_data")
+  data("xenium_data")
+
+  # get custom colors
+  colors <- scales::hue_pal()(length(unique(xenium_data$clusters)))
+  names(colors) <- unique(xenium_data$clusters)
+
+  # embedding plot
+  vrEmbeddingPlot(xenium_data, group.by = "clusters", embedding = "umap", label = T)
+  vrEmbeddingPlot(xenium_data, group.by = "clusters", embedding = "umap", group.ids = c(1,3,4), label = T)
+  vrEmbeddingPlot(xenium_data, group.by = "clusters", embedding = "umap", colors = colors, label = T)
+  vrEmbeddingPlot(xenium_data, group.by = "clusters", embedding = "umap", group.ids = c(1,3,4), colors = colors[c(1,3,4)], label = T)
+  vrEmbeddingPlot(xenium_data, group.by = "clusters", ncol = 3, split.by = "clusters")
+  vrEmbeddingPlot(xenium_data, group.by = "clusters", ncol = 3, split.by = "Sample")
+  expect_error(vrEmbeddingPlot(xenium_data, group.by = "clusters", ncol = 3, split.by = "art"))
+
+  # spatial plot
+  vrSpatialPlot(xenium_data, group.by = "clusters", plot.segments = TRUE)
+  vrSpatialPlot(xenium_data, group.by = "clusters", group.ids = c(1,3,4), plot.segments = TRUE)
+  vrSpatialPlot(xenium_data, group.by = "clusters", colors = colors, plot.segments = TRUE)
+  vrSpatialPlot(xenium_data, group.by = "clusters", group.ids = c(1,3,4), colors = colors[c(1,3,4)], plot.segments = TRUE)
+  vrSpatialPlot(xenium_data, group.by = "clusters", background = "black")
+  vrSpatialPlot(xenium_data, group.by = "clusters", background = "white")
+  vrSpatialPlot(xenium_data, group.by = "clusters", background = "main")
+  expect_warning(vrSpatialPlot(xenium_data, group.by = "clusters", background = c("main", "DAPI2")))
+  
+  # spatial plot without segmentation
+  vrSpatialPlot(xenium_data, group.by = "clusters", plot.segments = FALSE)
+
+  # spatial plot of visium
+  vrSpatialPlot(visium_data)
+
+  # spatial plot of melc data
+  vrSpatialPlot(melc_data, group.by = "Clusters")
+  expect_error(vrSpatialPlot(melc_data, group.by = "Clusters_new"))
+
+  # feature plots
+  vrSpatialFeaturePlot(visium_data, features = "Count")
+  vrSpatialFeaturePlot(visium_data, features = "Stat1", norm = TRUE, log = TRUE)
+  expect_error(vrSpatialFeaturePlot(visium_data, features = "Count_new"))
+  
+  # return
+  expect_equal(1,1L)
+})
+
+# Testing plotting functions
+test_that("missing_embedding_values", {
+  
+  # get data
+  data("visium_data")
+  data("xenium_data")
+  
+  # change embeddings 
+  vrEmbeddings(xenium_data, type = "new_umap") <- vrEmbeddings(xenium_data, type = "umap")[sample(1:length(vrSpatialPoints(xenium_data)), 500),]
+  vrEmbeddingPlot(xenium_data, embedding = "new_umap")
+  expect_error(vrEmbeddingPlot(xenium_data, embedding = "new_umap", group.by = "cluster"))
+  vrEmbeddingFeaturePlot(xenium_data, embedding = "new_umap", features = "Count")
+  expect_error(vrEmbeddingFeaturePlot(xenium_data, embedding = "new_umap", features = "Counts"))
+  vrEmbeddingFeaturePlot(xenium_data, embedding = "new_umap", features = "REXO4")
+  expect_error(vrEmbeddingFeaturePlot(xenium_data, embedding = "new_umap", features = "REXO4s"))
+  
+  # return
+  expect_equal(1,1L)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-spatialpoints.R b/tests/testthat/test-spatialpoints.R
new file mode 100644
index 00000000..971a8b29
--- /dev/null
+++ b/tests/testthat/test-spatialpoints.R
@@ -0,0 +1,16 @@
+# Testing functions of manipulating spatialpoints ####
+test_that("spatialpoints", {
+  
+  # get data
+  data("visium_data")
+  
+  # get spatial points
+  print(head(vrSpatialPoints(visium_data)))
+  print(head(vrSpatialPoints(visium_data, assay = "Assay1")))
+  
+  # subset on spatial points
+  spatialpoints <- vrSpatialPoints(visium_data)
+  visium_data_sub <- subset(visium_data, spatialpoints = spatialpoints[1:5])
+  
+  expect_equal(1,1L)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-subset.R b/tests/testthat/test-subset.R
new file mode 100644
index 00000000..c56212e5
--- /dev/null
+++ b/tests/testthat/test-subset.R
@@ -0,0 +1,26 @@
+# Testing functions of subsetting ####
+test_that("subset", {
+  
+  # get data
+  data("visium_data")
+  
+  # subset based on assay
+  subset(visium_data, assays = "Assay1")
+  subset(visium_data, assays = "Visium")
+  expect_error(subset(visium_data, assays = "Visium2"))
+  
+  # subset based on samples
+  subset(visium_data, samples = "Anterior1")
+  expect_error(subset(visium_data, samples = "Anterior2"))
+  
+  # subset based on assay
+  subset(visium_data, spatialpoints = c("GTTATATTATCTCCCT-1_Assay1", "GTTTGGGTTTCGCCCG-1_Assay1"))
+  expect_error(subset(visium_data, spatialpoints = c("point")))
+  
+  # subset based on features
+  subset(visium_data, features = c("Map3k19", "Rab3gap1"))
+  expect_error(subset(visium_data, features = c("feature")))
+  
+  # return
+  expect_equal(1,1L)
+})
\ No newline at end of file
diff --git a/tests/testthat/test-voltronobjects.R b/tests/testthat/test-voltronobjects.R
index 1d870f3f..0f5e0a3e 100644
--- a/tests/testthat/test-voltronobjects.R
+++ b/tests/testthat/test-voltronobjects.R
@@ -43,198 +43,3 @@ test_that("sample", {
 
   expect_equal(1,1L)
 })
-
-# Testing functions of manipulating spatialpoints ####
-test_that("spatialpoints", {
-
-  # get data
-  data("visium_data")
-
-  # get spatial points
-  print(head(vrSpatialPoints(visium_data)))
-  print(head(vrSpatialPoints(visium_data, assay = "Assay1")))
-
-  # subset on spatial points
-  spatialpoints <- vrSpatialPoints(visium_data)
-  visium_data_sub <- subset(visium_data, spatialpoints = spatialpoints[1:5])
-
-  expect_equal(1,1L)
-})
-
-# Testing functions of manipulating coordinates ####
-test_that("coordinates", {
-
-  # get data
-  data("visium_data")
-
-  # coordinates
-  coords <- vrCoordinates(visium_data)
-  coords <- vrCoordinates(visium_data, image_name = "main")
-  coords <- vrCoordinates(visium_data, spatial_name = "main")
-  expect_warning(coords <- vrCoordinates(visium_data, reg = TRUE))
-  expect_warning(coords <- vrCoordinates(visium_data, assay = "Assay1", reg = TRUE))
-
-  # update coordinates
-  vrCoordinates(visium_data) <- coords*2
-  expect_error(vrCoordinates(visium_data, reg = TRUE) <- coords*3)
-
-  # flip coordinates
-  visium_data <- flipCoordinates(visium_data)
-
-  # segments
-  segments <- vrSegments(visium_data)
-  expect_warning(segments <- vrSegments(visium_data, reg = TRUE))
-  expect_warning(segments <- vrSegments(visium_data, assay = "Assay1", reg = TRUE))
-
-  expect_equal(1,1L)
-})
-
-# Testing functions of manipulating images ####
-test_that("image", {
-
-  # get data
-  data("visium_data")
-
-  # get image
-  images <- vrImages(visium_data)
-  images <- vrImages(visium_data, name = "main")
-  expect_error(images <- vrImages(visium_data, name = "main2"))
-  images <- vrImages(visium_data, name = "main", channel = "H&E")
-  expect_warning(images <- vrImages(visium_data, name = "main", channel = "H&E2"))
-
-  # manipulate image
-  visium_data_resize <- resizeImage(visium_data, size = 400)
-  visium_data_modulate <- modulateImage(visium_data, brightness = 400)
-
-  # add new image
-  visium_data[["Assay1"]]@image[["new_image"]] <- vrImages(visium_data_resize)
-
-  # get image names
-  expect_equal(vrImageNames(visium_data), c("main", "new_image"))
-
-  # get main image
-  expect_equal(vrMainImage(visium_data[["Assay1"]]), "main")
-  expect_equal(vrMainSpatial(visium_data[["Assay1"]]), "main")
-  
-  # change main image
-  vrMainImage(visium_data[["Assay1"]]) <- "new_image"
-  vrMainSpatial(visium_data[["Assay1"]]) <- "new_image"
-  expect_equal(vrMainImage(visium_data[["Assay1"]]), "new_image")
-  expect_equal(vrMainSpatial(visium_data[["Assay1"]]), "new_image")
-
-  # return
-  expect_equal(1,1L)
-})
-
-# Testing functions of manipulating embeddings ####
-test_that("embeddings", {
-
-  # get data
-  data("visium_data")
-
-  # write embedding
-  vrEmbeddings(visium_data, type = "pca") <- vrCoordinates(visium_data)
-
-  # check overwrite embeddings
-  expect_error(vrEmbeddings(visium_data, type = "pca") <- vrCoordinates(visium_data))
-  vrEmbeddings(visium_data, type = "pca", overwrite = TRUE) <- vrCoordinates(visium_data)
-
-  # return
-  expect_equal(1,1L)
-})
-
-# Testing functions of image datasets ####
-test_that("importimagedata", {
-
-  # import image data
-  img <- system.file(package = "VoltRon", "extdata/DAPI.tif")
-  vrimgdata <- importImageData(img, tile.size = 10, stack.id = 1)
-  vrimgdata <- importImageData(img, tile.size = 10, stack.id = 2)
-
-  # import image data when image is a stack
-  img_magick <- magick::image_read(img)
-  img_stack <- c(img_magick, img_magick)
-  vrimgdata <- importImageData(img_stack, tile.size = 10, stack.id = 1)
-  vrimgdata <- importImageData(img_stack, tile.size = 10, stack.id = 2)
-  expect_error(vrimgdata <- importImageData(img_stack, tile.size = 10, stack.id = 3))
-
-  # return
-  expect_equal(1,1L)
-})
-
-# Testing functions of plots ####
-test_that("plots", {
-
-  # get data
-  data("visium_data")
-  data("xenium_data")
-
-  # get custom colors
-  colors <- scales::hue_pal()(length(unique(xenium_data$clusters)))
-  names(colors) <- unique(xenium_data$clusters)
-
-  # embedding plot
-  vrEmbeddingPlot(xenium_data, group.by = "clusters", embedding = "umap", label = T)
-  vrEmbeddingPlot(xenium_data, group.by = "clusters", embedding = "umap", group.ids = c(1,3,4), label = T)
-  vrEmbeddingPlot(xenium_data, group.by = "clusters", embedding = "umap", colors = colors, label = T)
-  vrEmbeddingPlot(xenium_data, group.by = "clusters", embedding = "umap", group.ids = c(1,3,4), colors = colors[c(1,3,4)], label = T)
-  vrEmbeddingPlot(xenium_data, group.by = "clusters", ncol = 3, split.by = "clusters")
-  vrEmbeddingPlot(xenium_data, group.by = "clusters", ncol = 3, split.by = "Sample")
-  expect_error(vrEmbeddingPlot(xenium_data, group.by = "clusters", ncol = 3, split.by = "art"))
-
-  # spatial plot
-  vrSpatialPlot(xenium_data, group.by = "clusters", plot.segments = TRUE)
-  vrSpatialPlot(xenium_data, group.by = "clusters", group.ids = c(1,3,4), plot.segments = TRUE)
-  vrSpatialPlot(xenium_data, group.by = "clusters", colors = colors, plot.segments = TRUE)
-  vrSpatialPlot(xenium_data, group.by = "clusters", group.ids = c(1,3,4), colors = colors[c(1,3,4)], plot.segments = TRUE)
-  vrSpatialPlot(xenium_data, group.by = "clusters", background = "black")
-  vrSpatialPlot(xenium_data, group.by = "clusters", background = "white")
-  vrSpatialPlot(xenium_data, group.by = "clusters", background = "main")
-  expect_warning(vrSpatialPlot(xenium_data, group.by = "clusters", background = c("main", "DAPI2")))
-  
-  # spatial plot without segmentation
-  vrSpatialPlot(xenium_data, group.by = "clusters", plot.segments = FALSE)
-
-  # spatial plot of visium
-  vrSpatialPlot(visium_data)
-
-  # spatial plot of melc data
-  vrSpatialPlot(melc_data, group.by = "Clusters")
-  expect_error(vrSpatialPlot(melc_data, group.by = "Clusters_new"))
-
-  # feature plots
-  vrSpatialFeaturePlot(visium_data, features = "Count")
-  vrSpatialFeaturePlot(visium_data, features = "Stat1", norm = TRUE, log = TRUE)
-  expect_error(vrSpatialFeaturePlot(visium_data, features = "Count_new"))
-  
-  # return
-  expect_equal(1,1L)
-})
-
-# Testing functions of subsetting ####
-test_that("subset", {
-  
-  # get data
-  data("visium_data")
-
-  # subset based on assay
-  subset(visium_data, assays = "Assay1")
-  subset(visium_data, assays = "Visium")
-  expect_error(subset(visium_data, assays = "Visium2"))
-  
-  # subset based on samples
-  subset(visium_data, samples = "Anterior1")
-  expect_error(subset(visium_data, samples = "Anterior2"))
-  
-  # subset based on assay
-  subset(visium_data, spatialpoints = c("GTTATATTATCTCCCT-1_Assay1", "GTTTGGGTTTCGCCCG-1_Assay1"))
-  expect_error(subset(visium_data, spatialpoints = c("point")))
-  
-  # subset based on features
-  subset(visium_data, features = c("Map3k19", "Rab3gap1"))
-  expect_error(subset(visium_data, features = c("feature")))
-  
-  # return
-  expect_equal(1,1L)
-})
-