diff --git a/NAMESPACE b/NAMESPACE
index 1458c52..bef684f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,21 +1,24 @@
 # Generated by roxygen2: do not edit by hand
 
-S3method("[",ome_zarr)
 S3method(dimnames,ome_zarr)
 S3method(plot,ome_zarr)
-S3method(print,ome_zarr)
-export(extract_levels)
 export(ome_read)
 export(ome_validate)
 export(ome_write)
+export(print.ome_zarr)
 exportMethods(ome_write)
+exportMethods(print)
+exportMethods(show)
+import(methods)
 importFrom(EBImage,Image)
 importFrom(EBImage,imageData)
 importFrom(EBImage,resize)
 importFrom(Rarr,read_zarr_array)
 importFrom(Rarr,read_zarr_attributes)
+importFrom(S4Vectors,SimpleList)
 importFrom(ZarrArray,ZarrArray)
 importFrom(cli,cli_abort)
+importFrom(methods,show)
 importFrom(stats,setNames)
 importFrom(utils,head)
 importFrom(utils,tail)
diff --git a/R/AllClasses.R b/R/AllClasses.R
new file mode 100644
index 0000000..4dc51c9
--- /dev/null
+++ b/R/AllClasses.R
@@ -0,0 +1,21 @@
+#' @importFrom S4Vectors SimpleList
+setClass(
+  Class = "ome_zarr",
+  contains = "SimpleList",
+  prototype = prototype(elementType = "array")
+)
+
+#' @importFrom utils head
+#' @export
+setMethod("print", "ome_zarr", function(x, level = 1, ...) {
+  cat(
+    "Multiscale OME-Zarr ",
+    x@metadata$type,
+    " (v",
+    x@metadata$version,
+    ") object.\n",
+    sep = ""
+  )
+  cat(sprintf("Scale: %d/%d", level, length(x)), "\n")
+  print(head(x[[level]], rep_len(5, length(dim(x[[level]]))), ...))
+})
diff --git a/R/extract_level.R b/R/extract_level.R
deleted file mode 100644
index ef29a66..0000000
--- a/R/extract_level.R
+++ /dev/null
@@ -1,39 +0,0 @@
-#' Extract specific levels from a multiscale `ome-zarr` object
-#'
-#' @param x An `ome-zarr` object.
-#' @param levels Integer vector specifying the levels to extract.
-#'
-#' @returns
-#' - If `levels` is of length 1, an array
-#' - If `levels` is of length more than 1, an `ome-zarr` object
-#'
-#' @returns An object of `ome_zarr` (OME-Zarr) class representing an
-#'  image or label pyramid.
-#'
-#' @examples
-#' omezarrzip <- system.file("extdata",
-#'                           "test_ngff_image_v04.ome.zarr.zip",
-#'                           package = "romeo")
-#' dir.create(td <- tempfile())
-#' unzip(omezarrzip, exdir = td)
-#' x <- ome_read(td)
-#' extract_levels(x, c(1, 3))
-#' extract_levels(x, 2)
-#'
-#' @export
-extract_levels <- function(x, levels) {
-  stopifnot(
-    inherits(x, "ome_zarr")
-  )
-  if (any(levels < 1) || any(levels > length(x))) {
-    stop("Level must be between 1 and ", length(x))
-  }
-  x <- lapply(levels, function(level) x[[level]])
-
-  if (length(x) == 1) {
-    return(x[[1]])
-  }
-
-  class(x) <- "ome_zarr"
-  x
-}
diff --git a/R/get_dim_names.R b/R/get_dim_names.R
index 26e7f89..e641384 100644
--- a/R/get_dim_names.R
+++ b/R/get_dim_names.R
@@ -14,5 +14,5 @@
 
 #' @export
 dimnames.ome_zarr <- function(x) {
-  attr(x, "dim_names")
+  x@metadata$dim_names
 }
diff --git a/R/print.R b/R/methods.R
similarity index 60%
rename from R/print.R
rename to R/methods.R
index 1897ddc..358bac4 100644
--- a/R/print.R
+++ b/R/methods.R
@@ -3,12 +3,19 @@
 print.ome_zarr <- function(x, level = 1, ...) {
   cat(
     "Multiscale OME-Zarr ",
-    attr(x, "type"),
+    x@metadata$type,
     " (v",
-    attr(x, "version"),
+    x@metadata$version,
     ") object.\n",
     sep = ""
   )
   cat(sprintf("Scale: %d/%d", level, length(x)), "\n")
   print(head(x[[level]], rep_len(5, length(dim(x[[level]]))), ...))
+  invisible(x)
 }
+
+#' @importFrom methods show
+#' @export
+setMethod("show", "ome_zarr", function(object) {
+  print.ome_zarr(object, level = 1)
+})
diff --git a/R/ome_read.R b/R/ome_read.R
index 4fca03d..c21b1fe 100644
--- a/R/ome_read.R
+++ b/R/ome_read.R
@@ -54,7 +54,7 @@ ome_read <- function(path, s3_client = NULL, lazy = TRUE, validate = TRUE) {
     img
   })
 
-  x <- mapply(
+  mapply(
     function(img, scale) {
       attr(img, "scale") <- scale
       img
@@ -64,13 +64,9 @@ ome_read <- function(path, s3_client = NULL, lazy = TRUE, validate = TRUE) {
       unlist(x$coordinateTransformations[[1]]$scale)
     }),
     SIMPLIFY = FALSE
-  )
-  class(x) <- "ome_zarr"
-  attr(x, "type") <- type
-  attr(x, "version") <- ome_version
-  if (!is.null(dim_names)) {
-    attr(x, "dim_names") <- dim_names
-  }
-
-  x
+  ) |>
+    new(
+      "ome_zarr",
+      metadata = list(version = ome_version, type = type, dim_names = dim_names)
+    )
 }
diff --git a/R/romeo-package.R b/R/romeo-package.R
index 3b3c951..a885738 100644
--- a/R/romeo-package.R
+++ b/R/romeo-package.R
@@ -2,6 +2,7 @@
 "_PACKAGE"
 
 ## usethis namespace: start
+#' @import methods
 ## usethis namespace: end
 
 ## mockable bindings: start
diff --git a/R/subset.R b/R/subset.R
deleted file mode 100644
index e28a5d4..0000000
--- a/R/subset.R
+++ /dev/null
@@ -1,50 +0,0 @@
-#' Subset an `ome-zarr` object
-#'
-#' Subset operation is applied on all levels of the multiscale `ome-zarr`
-#' object. The result is an `ome-zarr` object with the same number of levels,
-#' but each level is subsetted according to the provided indices.
-#'
-#' The first image is subsetted using the provided indices, and the
-#' resulting dimensions are used to subset the remaining levels, while
-#' conserving the same scaling factor across levels
-#'
-#' @param x An `ome-zarr` object.
-#' @param ... Indices to subset the `ome-zarr` object.
-#'
-#' @returns A subset of an object of `ome_zarr` (OME-Zarr) class representing
-#'  an image or label pyramid.
-#'
-#' @examples
-#' omezarrzip <- system.file("extdata",
-#'                           "test_ngff_image_v04.ome.zarr.zip",
-#'                           package = "romeo")
-#' dir.create(td <- tempfile())
-#' unzip(omezarrzip, exdir = td)
-#' x <- ome_read(td)
-#' y <- x[1:5,1:5]
-#' plot(y, level = 2)
-#'
-#' @export
-`[.ome_zarr` <- function(x, ...) {
-  y <- lapply(x, function(layer) {
-    scale <- attr(layer, "scale")
-    indices <- list(...)
-    indices <- mapply(
-      function(idx, scaling_factor) {
-        if (is.null(idx)) {
-          return(NULL)
-        }
-        # FIXME: is this the most sensible way to round here?
-        scaled_idx <- unique(ceiling((idx / scaling_factor)))
-        scaled_idx
-      },
-      indices,
-      scale,
-      SIMPLIFY = FALSE
-    )
-    do.call(`[`, c(list(layer), indices))
-  })
-  class(y) <- "ome_zarr"
-  attributes(y) <- attributes(x)
-  y
-}
diff --git a/man/extract_levels.Rd b/man/extract_levels.Rd
deleted file mode 100644
index 3e0387a..0000000
--- a/man/extract_levels.Rd
+++ /dev/null
@@ -1,36 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/extract_level.R
-\name{extract_levels}
-\alias{extract_levels}
-\title{Extract specific levels from a multiscale \code{ome-zarr} object}
-\usage{
-extract_levels(x, levels)
-}
-\arguments{
-\item{x}{An \code{ome-zarr} object.}
-
-\item{levels}{Integer vector specifying the levels to extract.}
-}
-\value{
-\itemize{
-\item If \code{levels} is of length 1, an array
-\item If \code{levels} is of length more than 1, an \code{ome-zarr} object
-}
-
-An object of \code{ome_zarr} (OME-Zarr) class representing an
-image or label pyramid.
-}
-\description{
-Extract specific levels from a multiscale \code{ome-zarr} object
-}
-\examples{
-omezarrzip <- system.file("extdata",
-                          "test_ngff_image_v04.ome.zarr.zip",
-                          package = "romeo")
-dir.create(td <- tempfile())
-unzip(omezarrzip, exdir = td)
-x <- ome_read(td)
-extract_levels(x, c(1, 3))
-extract_levels(x, 2)
-
-}
diff --git a/man/romeo-package.Rd b/man/romeo-package.Rd
index fc55fe0..5ee00f9 100644
--- a/man/romeo-package.Rd
+++ b/man/romeo-package.Rd
@@ -6,7 +6,7 @@
 \alias{romeo-package}
 \title{romeo: Minimal R 'OME-Zarr' Reader}
 \description{
-What the package does (one paragraph).
+A minimal R package to reading, writing and validating multiscale OME-Zarr images.
 }
 \seealso{
 Useful links:
@@ -26,5 +26,11 @@ Authors:
   \item Artür Manukyan \email{artur-man@hotmail.com} (\href{https://orcid.org/0000-0002-0441-9517}{ORCID})
 }
 
+Other contributors:
+\itemize{
+  \item Helmholtz Association (\href{https://ror.org/0281dp749}{ROR}) (Funded by the Helmholtz ScienceServe Initative to improve interoperability of bioimaging, single-cell and spatial omics data, https://www.helmholtz.de) [funder]
+  \item German Network for Bioinformatics Infrastructure - de.NBI (\href{https://ror.org/01vmpm840}{ROR}) [funder]
+}
+
 }
 \keyword{internal}
diff --git a/man/sub-.ome_zarr.Rd b/man/sub-.ome_zarr.Rd
deleted file mode 100644
index 261d60f..0000000
--- a/man/sub-.ome_zarr.Rd
+++ /dev/null
@@ -1,38 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/subset.R
-\name{[.ome_zarr}
-\alias{[.ome_zarr}
-\title{Subset an \code{ome-zarr} object}
-\usage{
-\method{[}{ome_zarr}(x, ...)
-}
-\arguments{
-\item{x}{An \code{ome-zarr} object.}
-
-\item{...}{Indices to subset the \code{ome-zarr} object.}
-}
-\value{
-A subset of an object of \code{ome_zarr} (OME-Zarr) class representing
-an image or label pyramid.
-}
-\description{
-Subset operation is applied on all levels of the multiscale \code{ome-zarr}
-object. The result is an \code{ome-zarr} object with the same number of levels,
-but each level is subsetted according to the provided indices.
-}
-\details{
-The first image is subsetted using the provided indices, and the
-resulting dimensions are used to subset the remaining levels, while
-conserving the same scaling factor across levels
-}
-\examples{
-omezarrzip <- system.file("extdata",
-                          "test_ngff_image_v04.ome.zarr.zip",
-                          package = "romeo")
-dir.create(td <- tempfile())
-unzip(omezarrzip, exdir = td)
-x <- ome_read(td)
-y <- x[1:5,1:5]
-plot(y, level = 2)
-
-}
diff --git a/tests/testthat/test-read.R b/tests/testthat/test-read.R
index d35f6c9..c9e81c7 100644
--- a/tests/testthat/test-read.R
+++ b/tests/testthat/test-read.R
@@ -19,14 +19,12 @@ test_that("parse ome version", {
 
     # image
     x <- ome_read(td)
-    # TODO: why S3 ?
-    expect_s3_class(x, "ome_zarr")
-    expect_equal(attr(x, "type"), "image")
+    expect_s4_class(x, "ome_zarr")
+    expect_identical(x@metadata$type, "image")
 
     # labels
     x <- ome_read(file.path(td, "labels/blobs"))
-    # TODO: why S3 ?
-    expect_s3_class(x, "ome_zarr")
-    expect_equal(attr(x, "type"), "label")
+    expect_s4_class(x, "ome_zarr")
+    expect_identical(x@metadata$type, "label")
   }
 })
diff --git a/tests/testthat/test-write-image.R b/tests/testthat/test-write-image.R
index 751ee4a..7b21ca2 100644
--- a/tests/testthat/test-write-image.R
+++ b/tests/testthat/test-write-image.R
@@ -85,7 +85,7 @@ test_that("writing 0.4 and 0.5 works", {
     )
 
     # check type
-    expect_equal(attr(ome_img, "type"), "image")
+    expect_identical(ome_img@metadata$type, "image")
 
     # zarr exists
     expect_true(zarr_exists(td))
diff --git a/tests/testthat/test-write-label.R b/tests/testthat/test-write-label.R
index 3b95840..a77c2bb 100644
--- a/tests/testthat/test-write-label.R
+++ b/tests/testthat/test-write-label.R
@@ -26,7 +26,7 @@ test_that("writing 0.4 and 0.5 labels works", {
     )
 
     # check type
-    expect_equal(attr(ome_label, "type"), "label")
+    expect_identical(ome_label@metadata$type, "label")
 
     # type is logical in this example
     expect_equal(type(ome_label[[1]]), "logical")