fix: map counts to layers['counts'] across all Seurat->h5ad paths

Counts were written to the legacy /raw/X group, which modern scverse tools
(scvi-tools, current scanpy workflows) do not look at. The in-memory direct
writers were migrated to /layers/counts, but the h5seurat-file and CLI paths
still emitted /raw/X, producing structurally different h5ad depending on the
entry point and silently dropping counts on a CLI round-trip.

Migrate the remaining writers and teach the CLI reader the modern convention:

- H5SeuratToH5AD (R hub): counts -> /layers/counts (dict 0.1.0 + csr_matrix
  0.1.0 + shape); drop the /raw and /raw/var machinery. Layers share /var.
- CLI writer (sc_h5ad.c): counts -> /layers/counts instead of /raw/X + raw/var.
- CLI reader (sc_h5ad.c): read /layers/counts into the assay counts layer and
  route X to data when a counts slot exists; /raw/X still read for back-compat.
- sc_stream_layers: skip "counts"/"data" so the primary matrices are not
  duplicated into a stray top-level /layers or collided on round-trip.
- Rename internal raw_list -> counts_list in .writeH5AD_c for accuracy.

Encoding verified byte-for-byte against anndata 0.12.7 (layers dict 0.1.0,
csr_matrix 0.1.0, shape; no dtype attr - anndata does not emit one). Full
round-trips validated for the C writer, R writer, R hub writer, and CLI both
directions; legacy /raw/X reading preserved. Test suites: cli-integration
117/0, roundtrip-property 193/0, saveload 19/0, misc-roundtrip 36/0, varp 18/0.

Author: mianaz

R/Convert.R

@@ -1291,6 +1291,18 @@ H5ADToH5Seurat <- function(
   # would otherwise be decoded per assay-copy pass.
   .categorical_cache <- new.env(parent = emptyenv())
 
+  # Counts source priority: raw/X (legacy h5ad) > layers/counts (modern
+  # anndata). When neither is present, X itself is treated as counts —
+  # matching scanpy's behaviour when only X exists.
+  counts_source <- if (source$exists(name = 'raw') &&
+                       source[['raw']]$exists(name = 'X')) {
+    'raw/X'
+  } else if (source$exists(name = 'layers') &&
+             source[['layers']]$exists(name = 'counts')) {
+    'layers/counts'
+  } else {
+    'X'
+  }
   ds.map <- c(
     scale.data = if (inherits(x = source[['X']], what = 'H5D')) {
       'X'
@@ -1299,11 +1311,7 @@ H5ADToH5Seurat <- function(
     },
     # Always use X for data slot - X contains normalized data when raw/X has counts
     data = 'X',
-    counts = if (source$exists(name = 'raw')) {
-      'raw/X'
-    } else {
-      'X'
-    }
+    counts = counts_source
   )
   # Add assay data
   assay.group <- dfile[['assays']]$create_group(name = assay)
@@ -3254,164 +3262,42 @@ H5SeuratToH5AD <- function(
     )
   }
 
-  # Add raw
+  # Add counts as layers/counts (modern anndata convention). Layers share the
+  # main /var, so no separate var group is written (unlike the legacy /raw).
   if (!is.null(x = raw.data)) {
     if (verbose) {
-      message("Adding ", raw.data, " from ", assay, " as raw")
-    }
-
-    # Get shape from source BEFORE copying (most reliable)
-    raw_source_dims <- NULL
-    if (grepl("/", raw.data)) {
-      # Handle nested paths like "layers/counts"
-      parts <- strsplit(raw.data, "/")[[1]]
-      src_obj <- assay.group
-      for (part in parts) {
-        src_obj <- src_obj[[part]]
-      }
-      if (src_obj$attr_exists(attr_name = 'dims')) {
-        raw_source_dims <- h5attr(x = src_obj, which = 'dims')
-      }
-    } else {
-      # Simple path
-      if (assay.group$exists(name = raw.data) &&
-          assay.group[[raw.data]]$attr_exists(attr_name = 'dims')) {
-        raw_source_dims <- h5attr(x = assay.group[[raw.data]], which = 'dims')
-      }
+      message("Adding ", raw.data, " from ", assay, " as layers/counts")
     }
 
-    dfile$create_group(name = 'raw')
-    CopyH5MatrixData(
-      src_group = assay.group,
-      src_path = raw.data,
-      dst_loc = dfile[['raw']],
-      dst_name = 'X',
-      verbose = verbose
-    )
-
-    # Ensure shape attribute exists for raw/X (required by anndata)
-    if (!dfile[['raw/X']]$attr_exists(attr_name = 'shape')) {
-      # Try to get from copied dims attribute first
-      if (dfile[['raw/X']]$attr_exists(attr_name = 'dims')) {
-        dims <- h5attr(x = dfile[['raw/X']], which = 'dims')
-        dfile[['raw/X']]$create_attr(
-          attr_name = 'shape',
-          robj = rev(x = dims),
-          dtype = GuessDType(x = dims)
-        )
-        dfile[['raw/X']]$attr_delete(attr_name = 'dims')
-      } else if (!is.null(raw_source_dims)) {
-        # Use dims from source
-        dfile[['raw/X']]$create_attr(
-          attr_name = 'shape',
-          robj = rev(x = raw_source_dims),
-          dtype = GuessDType(x = raw_source_dims)
-        )
-      }
-    }
-
-    AddEncoding(dname = 'raw/X')
-    # Get number of features from the ACTUAL raw/X matrix dimensions
-    raw_x_shape <- NULL
-    if (dfile[['raw/X']]$attr_exists(attr_name = 'shape')) {
-      raw_x_shape <- h5attr(x = dfile[['raw/X']], which = 'shape')
-    } else if (dfile[['raw/X']]$attr_exists(attr_name = 'dims')) {
-      raw_x_shape <- rev(h5attr(x = dfile[['raw/X']], which = 'dims'))
-    }
-    if (!is.null(raw_x_shape) && length(raw_x_shape) >= 2) {
-      n_raw_features <- raw_x_shape[2]
-    } else {
-      n_raw_features <- prod(assay.group[['features']]$dims)
-    }
-    # Add meta.features with validation
-    if (assay.group$exists(name = 'meta.features')) {
-      meta.features.src <- assay.group[['meta.features']]
-      if (inherits(x = meta.features.src, what = c('H5D', 'H5Group'))) {
-        TransferDF(
-          src = meta.features.src,
-          dname = 'raw/var',
-          index = seq.default(from = 1, to = n_raw_features)
-        )
-      } else {
-        message("meta.features is not a valid H5D or H5Group in raw, creating empty var")
-        dfile[['raw']]$create_group(name = 'var')
-      }
-    } else {
-      dfile[['raw']]$create_group(name = 'var')
-    }
-    # Add feature names
-    if (Exists(x = dfile[['raw/var']], name = rownames)) {
-      dfile[['raw/var']]$link_delete(name = rownames)
-    }
-    # Get feature names from correct location based on structure
-    if (has_layers && assay.group$exists(name = 'meta.data/_index')) {
-      # V5 structure: feature names are in meta.data/_index
-      features_data_raw <- SafeH5DRead(assay.group[['meta.data/_index']])
-    } else {
-      # Legacy structure: feature names are in features dataset
-      features_data_raw <- SafeH5DRead(assay.group[['features']])
-    }
-    # Subset to match actual raw/X dimensions
-    features_data_raw <- as.character(features_data_raw[seq_len(n_raw_features)])
-    dfile[['raw/var']]$create_dataset(
-      name = rownames,
-      robj = features_data_raw,
-      dtype = GuessDType(x = features_data_raw[1])
-    )
-    # Add encoding attributes to raw/var/_index (required for AnnData to read gene names correctly)
-    dfile[['raw/var']][[rownames]]$create_attr(
+    layers.grp <- dfile$create_group(name = 'layers')
+    layers.grp$create_attr(
       attr_name = 'encoding-type',
-      robj = 'string-array',
-      dtype = CachedGuessDType(x = 'string-array'),
+      robj = 'dict',
+      dtype = CachedGuessDType(x = 'dict'),
       space = ScalarSpace()
     )
-    dfile[['raw/var']][[rownames]]$create_attr(
+    layers.grp$create_attr(
       attr_name = 'encoding-version',
-      robj = '0.2.0',
-      dtype = CachedGuessDType(x = '0.2.0'),
+      robj = '0.1.0',
+      dtype = CachedGuessDType(x = '0.1.0'),
       space = ScalarSpace()
     )
-    dfile[['raw/var']]$create_attr(
-      attr_name = rownames,
-      robj = rownames,
-      dtype = GuessDType(x = rownames),
-      space = ScalarSpace()
+
+    CopyH5MatrixData(
+      src_group = assay.group,
+      src_path = raw.data,
+      dst_loc = layers.grp,
+      dst_name = 'counts',
+      verbose = verbose
     )
-    # Add dataframe encoding attributes to raw/var group (required for AnnData)
-    # Only add column-order if it doesn't already exist (TransferDF may have added it)
-    if (!dfile[['raw/var']]$attr_exists(attr_name = 'column-order')) {
-      # Get list of columns (everything except _index)
-      raw_var_cols <- setdiff(dfile[['raw/var']]$names, rownames)
-      if (length(raw_var_cols) > 0) {
-        dfile[['raw/var']]$create_attr(
-          attr_name = 'column-order',
-          robj = raw_var_cols,
-          dtype = GuessDType(x = raw_var_cols)
-        )
-      } else {
-        # Empty column-order for dataframes with only _index
-        dfile[['raw/var']]$create_attr(
-          attr_name = 'column-order',
-          robj = character(0),
-          dtype = StringType('utf8')
-        )
-      }
-    }
-    # Add encoding-type and encoding-version
-    encoding.info <- c('type' = 'dataframe', 'version' = '0.2.0')
-    names(x = encoding.info) <- paste0('encoding-', names(x = encoding.info))
-    for (i in seq_along(along.with = encoding.info)) {
-      attr.name <- names(x = encoding.info)[i]
-      attr.value <- encoding.info[i]
-      if (dfile[['raw/var']]$attr_exists(attr_name = attr.name)) {
-        dfile[['raw/var']]$attr_delete(attr_name = attr.name)
-      }
-      dfile[['raw/var']]$create_attr(
-        attr_name = attr.name,
-        robj = attr.value,
-        dtype = GuessDType(x = attr.value),
-        space = ScalarSpace()
-      )
+
+    # AddEncoding sets encoding-type/version and recomputes the anndata shape
+    # (n_obs, n_vars) from indptr plus the copied `dims` attribute. Drop the
+    # leftover h5seurat `dims` attribute afterwards for byte-parity with anndata.
+    AddEncoding(dname = 'layers/counts')
+    counts.grp <- layers.grp[['counts']]
+    if (counts.grp$attr_exists(attr_name = 'dims')) {
+      counts.grp$attr_delete(attr_name = 'dims')
     }
   }
   # Add cell-level metadata with validation

R/WriteH5AD.R

@@ -225,20 +225,21 @@ writeH5AD <- function(
   has_counts <- layer_has("counts")
   has_data   <- layer_has("data")
 
-  # anndata convention: X holds normalized data, raw/X holds counts. When both
-  # layers exist, route counts to the raw group; when only one exists, it
-  # becomes X directly.
+  # anndata convention: X holds normalized data, counts live in
+  # layers['counts']. When both layers exist, counts is passed to the C writer
+  # (which emits /layers/counts); when only one exists, it becomes X directly.
   mat <- if (has_data) GetAssayData(assay_obj, layer = "data")
          else if (has_counts) GetAssayData(assay_obj, layer = "counts")
          else NULL
   if (is.null(mat)) return(FALSE)
 
-  raw_list <- NULL
+  # `counts_list` carries the counts CSR arrays to C_write_h5ad's 8th argument.
+  counts_list <- NULL
   if (has_data && has_counts) {
-    raw_mat <- GetAssayData(assay_obj, layer = "counts")
-    raw_list <- list(i = raw_mat@i, p = raw_mat@p, x = raw_mat@x,
-                     dim = dim(raw_mat))
-    rm(raw_mat)
+    counts_mat <- GetAssayData(assay_obj, layer = "counts")
+    counts_list <- list(i = counts_mat@i, p = counts_mat@p, x = counts_mat@x,
+                        dim = dim(counts_mat))
+    rm(counts_mat)
   }
 
   mat_list <- list(
@@ -282,7 +283,7 @@ writeH5AD <- function(
 
   result <- .Call(C_write_h5ad,
     filename, mat_list, meta, reductions, graphs,
-    assay_name, as.integer(gzip_level), raw_list
+    assay_name, as.integer(gzip_level), counts_list
   )
 
   if (isTRUE(result) && verbose) {
@@ -374,22 +375,33 @@ DirectSeuratToH5AD <- function(
     !is.null(d) && prod(dim(d)) > 0
   }, error = function(e) FALSE)
 
+  # X is the primary expression matrix: normalized data when present,
+  # otherwise raw counts (keeps counts-only objects valid for readH5AD's
+  # required /X skeleton).
   if (has_data) {
     data_mat <- GetAssayData(object, assay = assay, layer = 'data')
     write_csr_group(dfile, "X", data_mat)
     rm(data_mat); gc(verbose = FALSE)
-    if (has_counts) {
-      counts_mat <- GetAssayData(object, assay = assay, layer = 'counts')
-      raw_grp <- dfile$create_group("raw")
-      write_csr_group(raw_grp, "X", counts_mat)
-      rm(counts_mat); gc(verbose = FALSE)
-    }
   } else if (has_counts) {
     counts_mat <- GetAssayData(object, assay = assay, layer = 'counts')
     write_csr_group(dfile, "X", counts_mat)
     rm(counts_mat); gc(verbose = FALSE)
   }
 
+  # counts slot -> data.layers["counts"] (modern anndata convention).
+  # Only when a separate data slot exists; counts-only objects already
+  # placed counts in /X above to preserve round-trip with the reader.
+  if (has_data && has_counts) {
+    counts_mat <- GetAssayData(object, assay = assay, layer = 'counts')
+    layers_grp <- dfile$create_group("layers")
+    layers_grp$create_attr(attr_name = 'encoding-type', robj = 'dict',
+                           dtype = CachedGuessDType('dict'), space = ScalarSpace())
+    layers_grp$create_attr(attr_name = 'encoding-version', robj = '0.1.0',
+                           dtype = CachedGuessDType('0.1.0'), space = ScalarSpace())
+    write_csr_group(layers_grp, "counts", counts_mat)
+    rm(counts_mat); gc(verbose = FALSE)
+  }
+
   # ========== obs (cell metadata) ==========
   if (verbose) message("  Writing obs...")
   meta <- object[[]]
@@ -420,12 +432,6 @@ DirectSeuratToH5AD <- function(
   }
   write_df_group(dfile, "var", meta_features, gene_names)
 
-  # Write raw/var if raw/X exists
-  if (has_data && has_counts) {
-    write_df_group(dfile[["raw"]], "var",
-                   data.frame(row.names = gene_names), gene_names)
-  }
-
   # ========== obsm (dimensional reductions) ==========
   reducs <- Reductions(object)
   if (length(reducs) > 0) {

src/sc_groups.c

@@ -636,6 +636,12 @@ int sc_stream_layers(hid_t src, hid_t dst, sc_direction_t dir, int gzip_level) {
         H5Gget_objname_by_idx(src_grp, i, name, sizeof(name));
         int obj_type = H5Gget_objtype_by_idx(src_grp, i);
 
+        /* "counts" and "data" are the primary matrices, handled explicitly as
+         * X / layers-counts. Skip them here so they are neither duplicated into
+         * a stray top-level /layers entry nor collide on round-trip. */
+        if (strcmp(name, "counts") == 0 || strcmp(name, "data") == 0)
+            continue;
+
         if (obj_type == H5G_GROUP) {
             /* Fast path: H5Ocopy for sparse layers (CSR structure preserved) */
             int copy_rc = sc_copy_group_h5ocopy(src_grp, name, dst_grp, name);