diff --git a/.Rbuildignore b/.Rbuildignore
index 0631158..e89f96a 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,4 +1,5 @@
 .github
+^CITATION\.cff$
 ^_pkgdown\.yml$
 ^docs$
 ^pkgdown$
diff --git a/CITATION.cff b/CITATION.cff
index c3fe17c..11fb26d 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -1,5 +1,11 @@
 cff-version: 1.2.0
-message: "If you use this software, please cite it as below."
+message: "If you use this software, please cite the MUON paper below."
+title: "MuDataSeurat"
+abstract: "MuData serialization for Seurat: read and write .h5ad/.h5mu files from R without a Python runtime."
+version: 0.0.1.0000
+url: "https://github.com/zqfang/MuDataSeurat"
+repository-code: "https://github.com/zqfang/MuDataSeurat"
+license: GPL-3.0
 authors:
 - family-names: "Bredikhin"
   given-names: "Danila"
@@ -7,10 +13,9 @@ authors:
 - family-names: "Kats"
   given-names: "Ilia"
   orcid: "https://orcid.org/0000-0001-5220-5671"
-title: "muon"
-version: 1.0.0
-date-released: 2021-06-01
-url: "https://github.com/scverse/muon"
+- family-names: "Fang"
+  given-names: "Zhuoqing"
+  orcid: "https://orcid.org/0000-0002-7418-1313"
 preferred-citation:
   type: article
   authors:
diff --git a/DESCRIPTION b/DESCRIPTION
index 6b4dbc1..f7c523b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -25,7 +25,9 @@ Suggests:
     rmarkdown,
     knitr,
     SeuratData,
-    BiocStyle
+    BiocStyle,
+    testthat (>= 3.0.0),
+    fs
 Config/Needs/website:
     r-lib/downlit
 Remotes: SeuratData=satijalab/seurat-data
@@ -33,4 +35,4 @@ VignetteBuilder: knitr
 License: GPL-3
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.1.2
+RoxygenNote: 7.3.3
diff --git a/R/ReadH5MU.R b/R/ReadH5MU.R
index 11bfcb6..a7792d0 100644
--- a/R/ReadH5MU.R
+++ b/R/ReadH5MU.R
@@ -219,7 +219,8 @@ ReadH5MU <- function(file) {
 
   # Create a Seurat object
   srt <- Seurat::CreateSeuratObject(subset(modalities[[1]], cells = obs_names), assay = names(modalities)[1])
-  for (modality in names(modalities)[2:length(modalities)]) {
+  # NOTE: [-1], not [2:length()], which yields c(NA, ..) for one modality
+  for (modality in names(modalities)[-1]) {
     srt[[modality]] <- subset(modalities[[modality]], cells = obs_names)
   }
 
diff --git a/R/ReadUtils.R b/R/ReadUtils.R
index e22f162..44c9e69 100644
--- a/R/ReadUtils.R
+++ b/R/ReadUtils.R
@@ -51,10 +51,7 @@ read_table_encv1 <- function(dataset, set_index = TRUE) {
         # Make factors out of categorical data
         ref <- values_attr$categories
         values_labels <- ref$dereference(obj = NULL)[[1]]
-        # NOTE: number of labels have to be strictly matching the number of unique integer values.
-        values_notna <- unique(values)
-        values_notna <- values_notna[!is.na(values_notna)]
-        values <- factor(as.integer(values), labels = values_labels$read()[1:length(values_notna)])
+        values <- decode_categorical(as.integer(values), values_labels$read())
       }
     }
     values
@@ -64,20 +61,26 @@ read_table_encv1 <- function(dataset, set_index = TRUE) {
   table
 }
 
+# Reconstruct a factor from AnnData categorical storage.
+# Codes are 0-based indices into `categories`; a negative code denotes NA.
+# The full category list is kept as the factor levels, so categories that are
+# not used by any observation are preserved rather than shifting the labels of
+# the categories that are used.
+decode_categorical <- function(codes, categories) {
+  codes <- as.integer(codes)
+  codes[!is.na(codes) & (codes < 0L | codes >= length(categories))] <- NA
+  factor(categories[codes + 1L], levels = categories)
+}
+
 read_column <- function(column, etype, eversion) {
   values <- NULL
   if (etype == "categorical") {
     if (eversion == "0.2.0") {
       codes <- column[["codes"]]$read()
       categories <- column[["categories"]]$read()
-
-      # NOTE: number of labels have to be strictly matching the number of unique integer values.
-      codes_notna <- unique(codes)
-      codes_notna <- codes_notna[!is.na(codes_notna)]
-
-      values <- factor(as.integer(codes), labels = categories[1:length(codes_notna)])
+      values <- decode_categorical(codes, categories)
     } else {
-      warning(paste0("Cannot recognise encoving-version ", eversion))
+      warning(paste0("Cannot recognise encoding-version ", eversion))
     }
   } else {
     values <- column$read()
@@ -213,9 +216,9 @@ read_layers_to_assay <- function(root, modalityname="") {
   }
 
   obs <- read_table(root[['obs']])
-  if (is("obs", "data.frame"))
-    rownames(obs) <- paste(modalityname, rownames(obs), sep="-")
-
+  # NOTE: obs names must NOT be prefixed with the modality name here.
+  # ReadH5MU takes the intersection of obs names across modalities to build a
+  # single Seurat object, so per-modality prefixes would make it empty.
   colnames(X) <- rownames(obs)
   rownames(X) <- rownames(var)
 
diff --git a/R/WriteH5MU.R b/R/WriteH5MU.R
index 62d8a64..e286582 100644
--- a/R/WriteH5MU.R
+++ b/R/WriteH5MU.R
@@ -110,8 +110,12 @@ WriteH5ADHelper <- function(object, assay, root, scale.data=FALSE, sparse.type="
     write_matrix(layers_group, "counts", x[["counts"]], sparse.type)
     write_matrix(root, "X", x[["data"]], sparse.type)
   } else {
-    which_x <- which(!is.null(x))
-    write_matrix(root, "X", x[[which_x]], sparse.type)
+    # Exactly one of counts/data/scale.data is present: write that one as X.
+    which_x <- which(!vapply(x, is.null, logical(1)))
+    if (length(which_x) == 0) {
+      stop(paste0("Assay ", assay, " has no data in counts, data or scale.data to write."))
+    }
+    write_matrix(root, "X", x[[which_x[1]]], sparse.type)
   }
 
   uns_group <- root$create_group("uns")
@@ -163,8 +167,8 @@ WriteH5ADHelper <- function(object, assay, root, scale.data=FALSE, sparse.type="
         }
 
         # Strip away modality name if the embedding starts with it
-        if (emb_assay == substr(red_name, 1, length(emb_assay))) {
-          red_name <- substr(red_name, length(emb_assay) + 1, length(red_name))
+        if (emb_assay == substr(red_name, 1, nchar(emb_assay))) {
+          red_name <- substr(red_name, nchar(emb_assay) + 1, nchar(red_name))
         }
       }
 
@@ -275,21 +279,25 @@ setMethod("WriteH5AD", "Seurat", function(object, file, assay = NULL, scale.data
   # Do not default to Seurat::DefaultAssay(object)
   # as it is not explicit, is hard to reason about,
   # and does not mean anything for MuData.
-  if (length(object@assays) > 1 && is.null(assay)) {
+  if (is.null(assay)) {
+    if (length(object@assays) > 1) {
+      h5$close()
+      stop(paste0(
+        "An assay to be written has to be provided, one of: ",
+        paste(names(object@assays), collapse = ", "),
+        ".\nUse WriteH5MU() to write all the modalities."
+      ))
+    }
+    assay <- names(object@assays)[1]
+  } else if (!assay %in% names(object@assays)) {
+    # Never fall back to another assay: writing a different assay than the one
+    # that was asked for would silently produce a file with the wrong data.
     h5$close()
     stop(paste0(
-      "An assay to be written has to be provided, one of: ",
+      "Assay ", assay, " not found. Available assays: ",
       paste(names(object@assays), collapse = ", "),
-      ".\nUse WriteH5MU() to write all the modalities."
+      "."
     ))
-  }  
-  else if  (length(object@assays) > 1 && (match(assay, names(object@assays), nomatch = 0) > 0)) {
-    all_assays <- names(object@assays)
-    idx = match(assay, all_assays, nomatch = 0)
-    assay = all_assays[idx]
-  }
-  else {
-    assay <- names(object@assays)[1]
   }
 
   # "Global" attributes such as metadata have to be written
@@ -343,9 +351,27 @@ setMethod("WriteH5MU", "Seurat", function(object, file, scale.data=FALSE, sparse
   names(var_names) <- modalities
   write_data_frame(h5, "var", do.call(c, var_names))
 
+  write_mod_maps(h5, modalities, nrow(obs), var_names)
+
   uns_group <- h5$create_group("uns")
   write_attribute(uns_group, "encoding-type", "dict")
   write_attribute(uns_group, "encoding-version", "0.1.0")
+
+  # obsm/varm/obsp/varp have to exist even when empty: MuData readers expect all
+  # of them, and creating them lazily produced files that could not be opened.
+  obsm_group <- h5$create_group("obsm")
+  write_attribute(obsm_group, "encoding-type", "dict")
+  write_attribute(obsm_group, "encoding-version", "0.1.0")
+  varm_group <- h5$create_group("varm")
+  write_attribute(varm_group, "encoding-type", "dict")
+  write_attribute(varm_group, "encoding-version", "0.1.0")
+  obsp_group <- h5$create_group("obsp")
+  write_attribute(obsp_group, "encoding-type", "dict")
+  write_attribute(obsp_group, "encoding-version", "0.1.0")
+  varp_group <- h5$create_group("varp")
+  write_attribute(varp_group, "encoding-type", "dict")
+  write_attribute(varp_group, "encoding-version", "0.1.0")
+
   # reductions -> .obsm
   # Reductions starting with modality name
   # that corresponds to the assay.used value
@@ -390,8 +416,8 @@ setMethod("WriteH5MU", "Seurat", function(object, file, scale.data=FALSE, sparse
 
 
         # Strip away modality name if the embedding starts with it
-        if (assay_emb == substr(red_name, 1, length(assay_emb))) {
-          red_name <- substr(red_name, length(assay_emb) + 1, length(red_name))
+        if (assay_emb == substr(red_name, 1, nchar(assay_emb))) {
+          red_name <- substr(red_name, nchar(assay_emb) + 1, nchar(red_name))
         }
       }
 
@@ -399,15 +425,7 @@ setMethod("WriteH5MU", "Seurat", function(object, file, scale.data=FALSE, sparse
         next
       }
 
-      if (!"obsm" %in% names(h5)) {
-        obsm <- h5$create_group("obsm")
-        write_attribute(obsm, "encoding-type", "dict")
-        write_attribute(obsm, "encoding-version", "0.1.0")  
-      } else {
-        obsm <- h5[["obsm"]]
-      }
-
-      write_matrix(obsm, paste0("X_", red_name), emb)
+      write_matrix(obsm_group, paste0("X_", red_name), emb)
 
       # loadings -> .varm
       if (!is.null(loadings) && ncol(loadings) == ncol(red)) {
@@ -416,19 +434,6 @@ setMethod("WriteH5MU", "Seurat", function(object, file, scale.data=FALSE, sparse
           varm_key <- OBSM2VARM[[paste0("X_", red_name)]]
         }
 
-        if (modality_specific) {
-          # this should have been written with WriteH5ADHelper
-          next
-        }
-
-        if (!"varm" %in% names(h5)) {
-          varm <- h5$create_group("varm")
-          write_attribute(varm, "encoding-type", "dict")
-          write_attribute(varm, "encoding-version", "0.1.0")   
-        } else {
-          varm <- h5[["varm"]]
-        }
-
         # If only a subset of features was used,
         # this has to be accounted for
         var_names_for_loadings <- do.call(c, var_names)
@@ -446,28 +451,19 @@ setMethod("WriteH5MU", "Seurat", function(object, file, scale.data=FALSE, sparse
           all_loadings <- loadings
         }
 
-        write_matrix(varm, varm_key, t(all_loadings))
+        write_matrix(varm_group, varm_key, t(all_loadings))
       }
 
       # stdev -> .uns[...]['variance']
+      # Modality-specific reductions have already been written by WriteH5ADHelper
+      # and skipped above, so only multimodal reductions reach this point.
       if (length(red@stdev) > 0) {
-        if (modality_specific) {
-          # REMOVE: this should have been written with WriteH5ADHelper
-          if (!red_name %in% names(h5[[paste0("mod/", assay_emb, "/uns")]])) {
-            uns <- h5$create_group(paste0("mod/", assay_emb, "/uns/", red_name))
-            write_attribute(uns, "encoding-type", "dict")
-            write_attribute(uns, "encoding-version", "0.1.0")
-          } else {
-            uns <- uns_group[[paste0("mod/", assay_emb, "/uns/", red_name)]]
-          }
+        if (!red_name %in% names(uns_group)) {
+          uns <- uns_group$create_group(red_name)
+          write_attribute(uns, "encoding-type", "dict")
+          write_attribute(uns, "encoding-version", "0.1.0")
         } else {
-          if (!red_name %in% names(uns_group)) {
-            uns <- uns_group$create_group(red_name)
-            write_attribute(uns, "encoding-type", "dict")
-            write_attribute(uns, "encoding-version", "0.1.0")
-          } else {
-            uns <- uns_group[[red_name]]
-          }
+          uns <- uns_group[[red_name]]
         }
         write_matrix(uns, "variance", red@stdev^2)
       }
@@ -476,9 +472,6 @@ setMethod("WriteH5MU", "Seurat", function(object, file, scale.data=FALSE, sparse
 
   # graphs -> .obsp
   if ('graphs' %in% slotNames(object)) {
-    obsp_group <- h5$create_group("obsp")
-    write_attribute(obsp_group, "encoding-type", "dict")
-    write_attribute(obsp_group, "encoding-version", "0.1.0")
     for (graph_name in names(object@graphs)) {
       graph <- object@graphs[[graph_name]]
 
diff --git a/R/WriteUtils.R b/R/WriteUtils.R
index a061c27..ff2f526 100644
--- a/R/WriteUtils.R
+++ b/R/WriteUtils.R
@@ -40,7 +40,7 @@ finalize_anndata <- function(h5, internal = FALSE) {
   filename <- h5$get_filename()
   h5$close_all()
   h5 <- file(filename, "r+b")
-  writeChar(paste0("anndata (format-version=", .mudataversion, ";creator=", .name, ";creator-version=", .version, ")"), h5)
+  writeChar(paste0("anndata (format-version=", .anndataversion, ";creator=", .name, ";creator-version=", .version, ")"), h5)
   close(h5)
 }
 
@@ -71,6 +71,14 @@ write_attribute <- function(obj, name, value, scalar = TRUE) {
 }
 
 write_matrix <- function(parent, key, mat, storage_sparse_type = "csr_matrix") {
+  # AnnData (>=0.8) has no nullable-string encoding, so a character vector with
+  # NAs is stored as categorical, where missing values are represented by the
+  # code -1. Previously NAs were coerced to the literal string "NaN", silently
+  # corrupting the data.
+  if (is.character(mat) && anyNA(mat)) {
+    return(write_matrix(parent, key, factor(mat), storage_sparse_type))
+  }
+
   if (is.matrix(mat) || is.vector(mat) || is.array(mat)) {
     hasna <- anyNA(mat)
     if (hasna && is.double(mat)) {
@@ -79,12 +87,6 @@ write_matrix <- function(parent, key, mat, storage_sparse_type = "csr_matrix") {
       hasna <- FALSE
     }
 
-    if (hasna && is.character(mat)) {
-      # FIXME: anndata spec to handle character with NAs ?
-      mat[is.na(mat)] <- NaN
-      hasna <- FALSE
-    }
-
     if (!hasna) {
       dset <- write_dataset(parent, key, mat)
       write_attribute(dset, "encoding-type", ifelse(is.character(mat), "string-array", "array"))
@@ -173,62 +175,33 @@ write_data_frame <- function(parent, key, attr_df) {
   }
 }
 
-# Only write _index (obs_names or var_names)
-write_names <- function(attr_group, attr_names) {
-  stype <- H5T_STRING$new(type = "c", size = Inf)
-  stype$set_cset("UTF-8")
-  attr_group$create_dataset("_index", attr_names, dtype = stype)
-
-  # Write attributes
-  attr_group$create_attr("_index", "_index", space = H5S$new("scalar"), dtype = stype)
-  attr_group$create_attr("encoding-type", "dataframe", space = H5S$new("scalar"), dtype = stype)
-  attr_group$create_attr("encoding-version", "0.2.0", space = H5S$new("scalar"), dtype = stype)
-  # When there are no columns, null buffer can't be written to a file.
-  attr_group$create_attr("column-order", dtype = h5types$H5T_NATIVE_DOUBLE, space = H5S$new("simple", 0, 0))
-}
-
-
-
-slot_writer <- function(h5group, mx, name) {
-  xt <- mx
-  if ("i" %in% slotNames(mx)) {
-    sparse_type <- ifelse(class(mx) == "dgCMatrix", "csc_matrix", "csr_matrix")
-    # sparse matrix
-    if (sparse_type == "csc_matrix") {
-      xt <- Matrix::t(mx)
-    } # transpose for anndata
-    mx_group <- h5group$create_group(name)
-    write_sparse_matrix(mx_group, xt, sparse_type)
-  } else {
-    # dense matrix, create a dataset
-    # h5group$create_dataset(name, mx)
-    # h5group in the root.
-    write_dense_matrix(h5group, xt, name)
+# MuData records, for each modality, the 1-based index of every global
+# observation/variable within that modality (0 when it is not present).
+# Readers require /obsmap and /varmap to be present.
+write_mod_maps <- function(h5, modalities, n_obs, var_names) {
+  obsmap_group <- h5$create_group("obsmap")
+  write_attribute(obsmap_group, "encoding-type", "dict")
+  write_attribute(obsmap_group, "encoding-version", "0.1.0")
+  varmap_group <- h5$create_group("varmap")
+  write_attribute(varmap_group, "encoding-type", "dict")
+  write_attribute(varmap_group, "encoding-version", "0.1.0")
+
+  n_var <- length(unlist(var_names, use.names = FALSE))
+  var_offset <- 0L
+  for (mod in modalities) {
+    # A Seurat object shares all of its cells across every assay.
+    write_matrix(obsmap_group, mod, seq_len(n_obs))
+
+    # The global var is the concatenation of the per-modality var_names,
+    # so each modality occupies one contiguous block of it.
+    n_mod_var <- length(var_names[[mod]])
+    varmap <- integer(n_var)
+    varmap[var_offset + seq_len(n_mod_var)] <- seq_len(n_mod_var)
+    write_matrix(varmap_group, mod, varmap)
+    var_offset <- var_offset + n_mod_var
   }
 }
 
-
-write_sparse_matrix <- function(root, x, sparse_type) {
-  stype <- H5T_STRING$new(type = "c", size = Inf)
-  stype$set_cset("UTF-8")
-  root$create_dataset("indices", x@i)
-  root$create_dataset("indptr", x@p)
-  root$create_dataset("data", x@x)
-  h5attr(root, "shape") <- dim(x)
-  root$create_attr("encoding-type", sparse_type, space = H5S$new("scalar"), dtype = stype)
-  root$create_attr("encoding-version", "0.1.0", space = H5S$new("scalar"), dtype = stype)
-}
-
-write_dense_matrix <- function(root, x, name) {
-  stype <- H5T_STRING$new(type = "c", size = Inf)
-  stype$set_cset("UTF-8")
-  dense <- root$create_dataset(name, x)
-  # h5attr(dense, "shape") <- dim(x)
-  dense$create_attr("encoding-type", "array", space = H5S$new("scalar"), dtype = stype)
-  dense$create_attr("encoding-version", "0.2.0", space = H5S$new("scalar"), dtype = stype)
-}
-
-
 reshape_scaled_data <- function(mat, var.meta, mat_name = "scale.data") {
   # If only a subset of features was used,
   # this has to be accounted for
diff --git a/README.md b/README.md
index bcf7806..75f882c 100644
--- a/README.md
+++ b/README.md
@@ -18,9 +18,17 @@ The original repository activity seems quite low, and unfortunately, the bugs ha
 4. Fixed stack overflow issue because of obs column containing NAs
    - skip columns with all NA value
    - fixed string array with NA
-4. Add two new keyword arguments to `WriteH5AD` and `WriteH5MU`: 
+5. Add two new keyword arguments to `WriteH5AD` and `WriteH5MU`: 
    - `scale.data`: whether write `scale.data` to `anndata/mudata` or not.
    - `sparse.type`: store `csc_matrix` or `csr_matrix` in `anndata/mudata`
+6. Correctness fixes
+   - `.h5mu` files are now readable by `mudata`: `obsm`, `varm`, `obsp`, `varp`,
+     `obsmap` and `varmap` are always written.
+   - Categorical (factor) columns keep their labels when a category is unused.
+   - `NA` in string columns is stored as a missing value (as `categorical`,
+     since anndata has no nullable-string encoding) instead of the text `"NaN"`.
+   - `WriteH5AD` now errors on an unknown `assay` instead of silently writing
+     a different one.
 
 
 ## Installation
@@ -56,7 +64,7 @@ WriteH5MU(seu, "export.h5mu", overwrite=TRUE)
 ### Read H5AD to Seurat
 
 ```R
-ReadH5AD()
-ReadH5MU()
+seu <- ReadH5AD("export.h5ad")
+seu <- ReadH5MU("export.h5mu")
 ```
 You may also use the native support of anndata in R: `anndataR::read_h5ad`
diff --git a/man/WriteH5AD.Rd b/man/WriteH5AD.Rd
index d0a4443..8447e66 100644
--- a/man/WriteH5AD.Rd
+++ b/man/WriteH5AD.Rd
@@ -5,9 +5,23 @@
 \alias{WriteH5AD,Seurat-method}
 \title{Write one assay to .h5ad}
 \usage{
-WriteH5AD(object, file, assay = NULL, overwrite = TRUE)
+WriteH5AD(
+  object,
+  file,
+  assay = NULL,
+  scale.data = FALSE,
+  sparse.type = "csr_matrix",
+  overwrite = TRUE
+)
 
-\S4method{WriteH5AD}{Seurat}(object, file, assay = NULL, overwrite = TRUE)
+\S4method{WriteH5AD}{Seurat}(
+  object,
+  file,
+  assay = NULL,
+  scale.data = FALSE,
+  sparse.type = "csr_matrix",
+  overwrite = TRUE
+)
 }
 \arguments{
 \item{object}{\code{Seurat} object.}
@@ -16,6 +30,10 @@ WriteH5AD(object, file, assay = NULL, overwrite = TRUE)
 
 \item{assay}{Assay to write; can be omitted if there is a single assay in the object.}
 
+\item{scale.data}{Boolen, wether to export scale.data.}
+
+\item{sparse.type}{String, save as csr_matrix or csc_matrix.}
+
 \item{overwrite}{Boolean value to indicate if to overwrite the \code{file} if it exists (\code{TRUE} by default).}
 }
 \description{
diff --git a/man/WriteH5ADHelper.Rd b/man/WriteH5ADHelper.Rd
index 42a7f49..70ba905 100644
--- a/man/WriteH5ADHelper.Rd
+++ b/man/WriteH5ADHelper.Rd
@@ -4,7 +4,14 @@
 \alias{WriteH5ADHelper}
 \title{A helper function to write a modality (an assay) to an .h5mu file}
 \usage{
-WriteH5ADHelper(object, assay, root, global = FALSE)
+WriteH5ADHelper(
+  object,
+  assay,
+  root,
+  scale.data = FALSE,
+  sparse.type = "csr_matrix",
+  global = FALSE
+)
 }
 \description{
 A helper function to write a modality (an assay) to an .h5mu file
diff --git a/man/WriteH5MU.Rd b/man/WriteH5MU.Rd
index 71aa099..a793cb0 100644
--- a/man/WriteH5MU.Rd
+++ b/man/WriteH5MU.Rd
@@ -5,15 +5,31 @@
 \alias{WriteH5MU,Seurat-method}
 \title{Create an .h5mu file with data from a \code{\link{Seurat}} object}
 \usage{
-WriteH5MU(object, file, overwrite = TRUE)
+WriteH5MU(
+  object,
+  file,
+  scale.data = FALSE,
+  sparse.type = "csr_matrix",
+  overwrite = TRUE
+)
 
-\S4method{WriteH5MU}{Seurat}(object, file, overwrite = TRUE)
+\S4method{WriteH5MU}{Seurat}(
+  object,
+  file,
+  scale.data = FALSE,
+  sparse.type = "csr_matrix",
+  overwrite = TRUE
+)
 }
 \arguments{
 \item{object}{\code{Seurat} object.}
 
 \item{file}{Path to the .h5mu file.}
 
+\item{scale.data}{Boolen, wether to export scale.data.}
+
+\item{sparse.type}{String, save as csr_matrix or csc_matrix.}
+
 \item{overwrite}{Boolean value to indicate if to overwrite the \code{file} if it exists (\code{TRUE} by default).}
 }
 \description{
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000..8579217
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,4 @@
+library(testthat)
+library(MuDataSeurat)
+
+test_check("MuDataSeurat")
diff --git a/tests/testthat/test_assays.R b/tests/testthat/test_assays.R
new file mode 100644
index 0000000..8fafdff
--- /dev/null
+++ b/tests/testthat/test_assays.R
@@ -0,0 +1,128 @@
+library(Seurat)
+library(MuDataSeurat)
+library(Matrix)
+library(hdf5r)
+library(fs)
+
+# Regression tests for assay selection and layer handling.
+# Seurat v5 (Assay5) is the baseline.
+
+nobs <- 10
+nvar <- 20
+
+obs_names <- paste("obs", 1:nobs, sep = "-")
+var_names <- paste("var", 1:nvar, sep = "-")
+
+make_counts <- function(vars = var_names) {
+  x <- rnbinom(n = nobs * length(vars), prob = .95, size = 10)
+  x <- Matrix(matrix(x, ncol = nobs), sparse = TRUE)
+  colnames(x) <- obs_names
+  rownames(x) <- vars
+  x
+}
+
+make_multimodal <- function() {
+  adt_names <- paste("adt", 1:5, sep = "-")
+  srt <- CreateSeuratObject(counts = make_counts(), assay = "RNA")
+  srt[["ADT"]] <- CreateAssay5Object(counts = make_counts(adt_names))
+  srt
+}
+
+test_that("WriteH5AD errors on an assay that does not exist", {
+  srt <- make_multimodal()
+  file <- paste0(file_temp(), ".h5ad")
+
+  # Must not silently fall back to the first assay.
+  expect_error(WriteH5AD(srt, file, assay = "rna"), "not found")
+})
+
+test_that("WriteH5AD errors when the assay is ambiguous", {
+  srt <- make_multimodal()
+  file <- paste0(file_temp(), ".h5ad")
+
+  expect_error(WriteH5AD(srt, file), "has to be provided")
+})
+
+test_that("WriteH5AD writes the requested assay, not the first one", {
+  srt <- make_multimodal()
+  file <- paste0(file_temp(), ".h5ad")
+
+  expect_true(WriteH5AD(srt, file, assay = "ADT"))
+
+  srt2 <- ReadH5AD(file)
+  expect_equal(nrow(srt2), 5)
+  expect_equal(rownames(srt2), paste("adt", 1:5, sep = "-"))
+})
+
+test_that("an assay with only a data layer can be written", {
+  # The fallback branch used to always pick `counts`, which is NULL here.
+  x <- make_counts()
+  # Seurat warns about the missing counts layer when deriving nCount/nFeature.
+  srt <- suppressWarnings(CreateSeuratObject(CreateAssay5Object(data = x)))
+
+  file <- paste0(file_temp(), ".h5ad")
+  expect_true(WriteH5AD(srt, file))
+
+  srt2 <- ReadH5AD(file)
+  expect_equal(dim(srt2), c(nvar, nobs))
+})
+
+test_that("a multimodal object round-trips through .h5mu", {
+  srt <- make_multimodal()
+  file <- paste0(file_temp(), ".h5mu")
+
+  expect_true(WriteH5MU(srt, file))
+
+  srt2 <- ReadH5MU(file)
+  expect_setequal(Assays(srt2), c("RNA", "ADT"))
+  expect_equal(colnames(srt2), obs_names)
+  expect_equal(nrow(srt2[["RNA"]]), nvar)
+  expect_equal(nrow(srt2[["ADT"]]), 5)
+})
+
+test_that(".h5mu contains every group MuData readers require", {
+  # These used to be created lazily (or not at all), which produced files that
+  # mudata could not open unless multimodal reductions happened to be present.
+  srt <- make_multimodal()
+  file <- paste0(file_temp(), ".h5mu")
+  expect_true(WriteH5MU(srt, file))
+
+  h5 <- H5File$new(file, mode = "r")
+  on.exit(h5$close_all())
+
+  required <- c("mod", "obs", "var", "obsm", "varm",
+                "obsp", "varp", "obsmap", "varmap", "uns")
+  expect_true(all(required %in% names(h5)))
+  for (grp in c("obsm", "varm", "obsp", "varp", "obsmap", "varmap")) {
+    expect_equal(h5attr(h5[[grp]], "encoding-type"), "dict")
+  }
+})
+
+test_that(".h5mu obsmap/varmap index the modalities correctly", {
+  srt <- make_multimodal()
+  file <- paste0(file_temp(), ".h5mu")
+  expect_true(WriteH5MU(srt, file))
+
+  h5 <- H5File$new(file, mode = "r")
+  on.exit(h5$close_all())
+
+  # Every cell is present in every assay
+  expect_equal(as.integer(h5[["obsmap/RNA"]]$read()), 1:nobs)
+  expect_equal(as.integer(h5[["obsmap/ADT"]]$read()), 1:nobs)
+
+  # Global var is RNA (20) followed by ADT (5)
+  expect_equal(as.integer(h5[["varmap/RNA"]]$read()), c(1:nvar, rep(0, 5)))
+  expect_equal(as.integer(h5[["varmap/ADT"]]$read()), c(rep(0, nvar), 1:5))
+})
+
+test_that("a single-modality .h5mu can be read", {
+  # names(modalities)[2:length()] used to yield c(NA, "RNA") here.
+  srt <- CreateSeuratObject(counts = make_counts(), assay = "RNA")
+  file <- paste0(file_temp(), ".h5mu")
+
+  expect_true(WriteH5MU(srt, file))
+
+  srt2 <- ReadH5MU(file)
+  expect_equal(Assays(srt2), "RNA")
+  expect_equal(colnames(srt2), obs_names)
+})
diff --git a/tests/testthat/test_metadata.R b/tests/testthat/test_metadata.R
new file mode 100644
index 0000000..f4e1faf
--- /dev/null
+++ b/tests/testthat/test_metadata.R
@@ -0,0 +1,85 @@
+library(Seurat)
+library(MuDataSeurat)
+library(Matrix)
+library(hdf5r)
+library(fs)
+
+# Regression tests for metadata (obs) round-trips.
+# Seurat v5 (Assay5) is the baseline.
+
+nobs <- 10
+nvar <- 20
+
+obs_names <- paste("obs", 1:nobs, sep = "-")
+var_names <- paste("var", 1:nvar, sep = "-")
+
+make_srt <- function() {
+  x <- rnbinom(n = nobs * nvar, prob = .95, size = 10)
+  x <- Matrix(matrix(x, ncol = nobs), sparse = TRUE)
+  colnames(x) <- obs_names
+  rownames(x) <- var_names
+  CreateSeuratObject(counts = x)
+}
+
+test_that("factor columns round-trip when a category is unused", {
+  # Only "a" and "c" are observed, but "b" is a declared level.
+  # Truncating the category list by the number of observed codes used to
+  # relabel code 2 ("c") as "b".
+  celltype <- factor(rep(c("a", "c"), length.out = nobs),
+                     levels = c("a", "b", "c"))
+
+  srt <- make_srt()
+  srt$celltype <- celltype
+
+  file <- paste0(file_temp(), ".h5ad")
+  expect_true(WriteH5AD(srt, file))
+
+  srt2 <- ReadH5AD(file)
+  expect_equal(as.character(srt2$celltype), as.character(celltype))
+  expect_equal(levels(srt2$celltype), levels(celltype))
+})
+
+test_that("factor columns with NAs round-trip", {
+  celltype <- factor(c("a", NA, "b", NA, rep("a", nobs - 4)),
+                     levels = c("a", "b"))
+
+  srt <- make_srt()
+  srt$celltype <- celltype
+
+  file <- paste0(file_temp(), ".h5ad")
+  expect_true(WriteH5AD(srt, file))
+
+  srt2 <- ReadH5AD(file)
+  expect_equal(as.character(srt2$celltype), as.character(celltype))
+  expect_equal(unname(is.na(srt2$celltype)), is.na(celltype))
+})
+
+test_that("character columns with NAs round-trip instead of becoming 'NaN'", {
+  # AnnData has no nullable-string encoding, so these are stored as
+  # categorical. They previously came back as the literal string "NaN".
+  batch <- c("b1", NA, "b2", rep("b1", nobs - 3))
+
+  srt <- make_srt()
+  srt$batch <- batch
+
+  file <- paste0(file_temp(), ".h5ad")
+  expect_true(WriteH5AD(srt, file))
+
+  srt2 <- ReadH5AD(file)
+  expect_equal(as.character(srt2$batch), batch)
+  expect_true(is.na(as.character(srt2$batch)[2]))
+  expect_false(any(as.character(srt2$batch) %in% "NaN"))
+})
+
+test_that("character columns without NAs round-trip as strings", {
+  batch <- rep(c("b1", "b2"), length.out = nobs)
+
+  srt <- make_srt()
+  srt$batch <- batch
+
+  file <- paste0(file_temp(), ".h5ad")
+  expect_true(WriteH5AD(srt, file))
+
+  srt2 <- ReadH5AD(file)
+  expect_equal(as.character(srt2$batch), batch)
+})
diff --git a/tests/testthat/test_sparse.R b/tests/testthat/test_sparse.R
index 89232e0..2ab6b39 100644
--- a/tests/testthat/test_sparse.R
+++ b/tests/testthat/test_sparse.R
@@ -1,110 +1,80 @@
-context("Creating .h5ad and .h5mu files with sparse matrices")
 library(Seurat)
 library(MuDataSeurat)
 library(Matrix)
 library(hdf5r)
 library(fs)  # for file_temp()
 
-# csc
-fileh5mu_r <- paste0(file_temp(), ".h5mu")
-fileh5ad_r <- paste0(file_temp(), ".h5ad")
-# csc
-fileh5mu_c <- paste0(file_temp(), ".h5mu")
-fileh5ad_c <- paste0(file_temp(), ".h5ad")
-
+# Sparse matrix storage round-trips.
+# Seurat stores features x cells; AnnData stores obs x var, so the shape
+# attribute is reversed on write. The default storage type is csr_matrix.
 
 nobs <- 10
 nvar <- 20
-nvar2 <- 31
 
-obs_names <- paste("obs", 1:nobs, sep="-")
-var_names <- paste("var", 1:nvar, sep="-")
+obs_names <- paste("obs", 1:nobs, sep = "-")
+var_names <- paste("var", 1:nvar, sep = "-")
 
 # Hard-coded value to inject
 true_val <- 0.1234569
 true_val_i <- 3
 true_val_j <- 7
 
-test_that("dgCMatrix can be written to .h5ad", {
-    x <- rnbinom(n = nobs * nvar, prob = .95, size = 10)
-    x <- Matrix(matrix(x, ncol = nobs), sparse = TRUE)  # => dgCMatrix
-    x[true_val_i,true_val_j] <- true_val
-
-    colnames(x) <- obs_names
-    rownames(x) <- var_names
+make_srt <- function() {
+  x <- rnbinom(n = nobs * nvar, prob = .95, size = 10)
+  x <- Matrix(matrix(x, ncol = nobs), sparse = TRUE) # => dgCMatrix
+  x[true_val_i, true_val_j] <- true_val
+  colnames(x) <- obs_names
+  rownames(x) <- var_names
+  CreateSeuratObject(counts = x)
+}
 
-    srt <- CreateSeuratObject(counts = x)
+fileh5ad_r <- paste0(file_temp(), ".h5ad")
+fileh5ad_c <- paste0(file_temp(), ".h5ad")
 
-    expect_true(WriteH5AD(srt, fileh5ad_c))
+test_that("dgCMatrix is written to .h5ad as csr_matrix by default", {
+  expect_true(WriteH5AD(make_srt(), fileh5ad_r))
 
-    h5 <- H5File$new(fileh5ad_c, mode="a")
-    expect_equal(h5attr(h5[["X"]], "encoding-type"), "csc_matrix")
-    h5$close()
+  h5 <- H5File$new(fileh5ad_r, mode = "r")
+  on.exit(h5$close_all())
+  expect_equal(h5attr(h5[["X"]], "encoding-type"), "csr_matrix")
+  expect_equal(as.integer(h5attr(h5[["X"]], "shape")), c(nobs, nvar))
 })
 
-test_that("dgRMatrix can be written to .h5ad", {
-    # Seurat objects would not accept dgRMatrix,
-    # and there is no native dgRMatrix -> dgCMatrix conversion.
-    # We will write a transposed matrix as CSC and then change it to CSR manually
-    # to emulate CSR matrices written by other tools.
-
-    x <- rnbinom(n = nobs * nvar, prob = .95, size = 10)
-    x <- x[x != 0]
-    i <- sample(1:nvar, length(x), replace = T)
-    j <- sample(1:nobs, length(x), replace = T)
-    
-    x_c <- sparseMatrix(i = i, j = j, x = x, dims = c(nvar, nobs), repr = "C")    # => dgCMatrix
-    x_c[true_val_i,true_val_j] <- true_val
-
-    x_r <- sparseMatrix(i = i, j = j, x = x, dims = c(nvar, nobs), repr = "R")  # => dgRMatrix
-    # Value assignment will automatically convert it to dgTMatrix,
-    # and there is no dgTMatrix -> dgRMatrix coercion,
-    # so in this case we won't inject the value.
-    # x_r[true_val_i,true_val_j] <- true_val
+test_that("dgCMatrix can be written to .h5ad as csc_matrix", {
+  expect_true(WriteH5AD(make_srt(), fileh5ad_c, sparse.type = "csc_matrix"))
 
-    colnames(x_c) <- obs_names
-    rownames(x_c) <- var_names
-
-    srt <- CreateSeuratObject(counts = x_c)
-
-    expect_true(WriteH5AD(srt, fileh5ad_r))
-
-    h5 <- H5File$new(fileh5ad_r, mode="a")
-    h5x <- h5[["X"]]
-    expect_equal(h5attr(h5x, "encoding-type"), "csc_matrix")
-
-    h5x$link_delete("indices")
-    h5x[["indices"]] <- x_r@j
-    h5x$link_delete("indptr")
-    h5x[["indptr"]] <- x_r@p
-    h5x$link_delete("data")
-    h5x[["data"]] <- x_r@x
-    h5attr(h5x, "encoding-type") <- "csr_matrix"
-
-    h5$close()
+  h5 <- H5File$new(fileh5ad_c, mode = "r")
+  on.exit(h5$close_all())
+  expect_equal(h5attr(h5[["X"]], "encoding-type"), "csc_matrix")
+  expect_equal(as.integer(h5attr(h5[["X"]], "shape")), c(nobs, nvar))
 })
 
-test_that("dgCMatrix can be read from .h5ad", {
-    srt <- ReadH5AD(fileh5ad_c)
-
-    counts <- srt@assays[[1]]@counts
-    expect_true("dgCMatrix" %in% class(counts))
-    expect_equal(dim(counts), c(nvar, nobs))
-    expect_equal(counts[true_val_i, true_val_j], true_val)
-    expect_equal(rownames(srt), var_names)
-    expect_equal(colnames(srt), obs_names)
+test_that("an unsupported sparse.type is rejected", {
+  expect_error(
+    WriteH5AD(make_srt(), paste0(file_temp(), ".h5ad"), sparse.type = "coo_matrix"),
+    "not supported"
+  )
 })
 
-test_that("dgRMatrix can be read from .h5ad", {
-    srt <- ReadH5AD(fileh5ad_r)
-    counts <- srt@assays[[1]]@counts
-    
-    # Seurat only support dgCMatrix as counts
-    expect_true("dgCMatrix" %in% class(counts))
-    expect_equal(dim(counts), c(nvar, nobs))
-    
-    expect_equal(rownames(srt), var_names)
-    expect_equal(colnames(srt), obs_names)
+test_that("a csr_matrix .h5ad can be read", {
+  srt <- ReadH5AD(fileh5ad_r)
+  counts <- GetAssayData(srt, layer = "counts")
+
+  # Seurat only supports dgCMatrix as counts
+  expect_true("dgCMatrix" %in% class(counts))
+  expect_equal(dim(counts), c(nvar, nobs))
+  expect_equal(counts[true_val_i, true_val_j], true_val)
+  expect_equal(rownames(srt), var_names)
+  expect_equal(colnames(srt), obs_names)
 })
 
+test_that("a csc_matrix .h5ad can be read", {
+  srt <- ReadH5AD(fileh5ad_c)
+  counts <- GetAssayData(srt, layer = "counts")
 
+  expect_true("dgCMatrix" %in% class(counts))
+  expect_equal(dim(counts), c(nvar, nobs))
+  expect_equal(counts[true_val_i, true_val_j], true_val)
+  expect_equal(rownames(srt), var_names)
+  expect_equal(colnames(srt), obs_names)
+})