diff --git a/DESCRIPTION b/DESCRIPTION index f2a4415..9107392 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -80,6 +80,7 @@ Collate: 'Transpose.R' 'V5Compatibility.R' 'V5LayerSupport.R' + 'VerifyH5AD.R' 'WriteH5AD.R' 'ZarrRemote.R' 'ZarrStore.R' diff --git a/NEWS.md b/NEWS.md index 6b360ae..ac92fb8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,72 @@ # scConvert 0.3.0 (development) +## Reverse-conversion (h5ad -> Seurat) integrity overhaul + +`readH5AD()` now applies the same read-back discipline to the reverse +direction that the write direction already learned the hard way. New +helpers live in `R/VerifyH5AD.R`; regression coverage in +`tests/testthat/test-reverse-conversion.R`. + +- **Version/layout-aware counts resolution.** Where the raw counts live is + now resolved explicitly, in priority order: the `/uns/scConvert/ + counts_location` stamp (see below), `/layers/counts` (modern + anndata/scanpy convention), `/raw/X` (legacy scanpy), then `/X`. The + resolved slot is loaded into the Seurat `counts` layer and X into `data` + in every path -- in-memory, `components` subsets, the C-reader fallback, + and BPCells on-disk mode (which previously only knew about `raw/X` and + silently served log-normalized X as on-disk "counts" for + `layers['counts']` files). Relocating counts out of `layers/counts` is + now a hard `scConvert_data_error` on failure instead of a swallowed + message that left log-normalized values sitting in the counts slot. + +- **Counts integrality guard.** If the values that end up in the counts + layer are non-integer, `readH5AD()` raises a classed + `scConvert_counts_warning` naming the slot they came from, instead of + silently handing normalized data to downstream steps that assume raw + counts. + +- **Writer provenance stamp.** All h5ad writers (`.writeH5AD_c`, + `DirectSeuratToH5AD`, `H5SeuratToH5AD`) now stamp + `/uns/scConvert/{version, counts_location}` by inspecting the file just + written -- not the writer's intent -- so future reads branch on recorded + fact rather than layout heuristics. + +- **Post-read verification.** Before returning, the loaded object is + asserted against the file: dims must match the file's obs/var counts + (transpose/orientation check), cell names must equal the obs index in + order, feature names must equal the var index modulo Seurat's documented + underscore-to-dash replacement, and no duplicate barcodes or genes may + survive. Violations raise `scConvert_data_error`. What the reader did + (counts source, where X went, scConvert version, dedup flags) is + recorded in `misc$scConvert_read`. + +- **Duplicate names are loud.** Duplicate cell barcodes or feature names + in the file are made unique with a classed `scConvert_names_warning` + (previously features were deduplicated silently and duplicate barcodes + fell through to Seurat with undefined downstream alignment). + +- **Gene identity is never silently lost.** When final feature names + differ from the file's var index (deduplication or Seurat's underscore + mangling), the original identifiers are preserved in a new + `orig_var_index` feature-metadata column. var metadata columns are now + keyed by the object's final rownames, fixing silent misalignment (and + layer drops) for files with underscores in gene names. + +- **Reduction key collisions and remapping.** obsm keys that clean to the + same reduction name (e.g. `X_pca` and `pca`) no longer silently + overwrite each other; later claimants keep their raw obsm key (made + unique if needed) with a `scConvert_reduction_warning`. A new + `reductions` argument to `readH5AD()` selects a subset of obsm keys + and/or renames them (`reductions = c(scvi = "X_scVI")`), replacing blind + full-object import with explicit key remapping. + +- **Categorical order and orderedness round-trip.** `DecodeCategorical()` + keeps the stored category order as the factor level order (now pinned by + tests against non-alphabetical orders) and gains an `ordered` argument; + both hdf5r reader paths read the AnnData `ordered` flag, and the + C-reader path restores it post-hoc, so ordered pandas categoricals come + back as ordered factors. + ## New features - **`readZarr()` index slicing: `obs_idx` and `var_idx` push down to diff --git a/R/AnnDataEncoding.R b/R/AnnDataEncoding.R index d75b0c4..33c54c0 100644 --- a/R/AnnDataEncoding.R +++ b/R/AnnDataEncoding.R @@ -8,14 +8,22 @@ NULL #' Decode AnnData categorical encoding to R factor #' +#' The factor levels follow the stored category order verbatim (pandas +#' preserves an explicit category order; re-sorting it alphabetically would +#' silently reorder positional palettes and any order-dependent downstream +#' logic). When the AnnData categorical carries \code{ordered = TRUE}, pass +#' it here to get an ordered factor back. +#' #' @param codes Integer vector of 0-based category codes (-1 = NA) #' @param categories Character vector of category labels +#' @param ordered Logical; produce an ordered factor (AnnData's +#' \code{ordered} categorical flag). Default \code{FALSE}. #' #' @return A factor vector #' #' @keywords internal #' -DecodeCategorical <- function(codes, categories) { +DecodeCategorical <- function(codes, categories, ordered = FALSE) { codes[codes == -1L] <- NA_integer_ valid <- !is.na(codes) & codes >= 0L & codes < length(categories) if (!all(valid[!is.na(codes)])) { @@ -26,7 +34,7 @@ DecodeCategorical <- function(codes, categories) { )) codes[!is.na(codes) & !valid] <- NA_integer_ } - factor(categories[codes + 1L], levels = categories) + factor(categories[codes + 1L], levels = categories, ordered = isTRUE(ordered)) } #' Encode R factor as AnnData categorical diff --git a/R/Convert.R b/R/Convert.R index 0e6fa43..c2cd848 100644 --- a/R/Convert.R +++ b/R/Convert.R @@ -4023,6 +4023,11 @@ H5SeuratToH5AD <- function( } } + # Provenance stamp: /uns/scConvert/{version, counts_location}, derived by + # inspecting the file just written, so readers resolve the counts layer + # from a version-aware fact instead of layout heuristics. + .h5ad_stamp_provenance(dfile) + dfile$flush() return(dfile) } diff --git a/R/LoadH5AD.R b/R/LoadH5AD.R index f59ade9..ae36d25 100644 --- a/R/LoadH5AD.R +++ b/R/LoadH5AD.R @@ -109,10 +109,31 @@ #' @param use.c Use compiled C reader when available (default: TRUE). Set to #' FALSE to force the pure-R hdf5r path. #' @param verbose Show progress messages +#' @param reductions Which \code{obsm} entries to load as dimensional +#' reductions. \code{NULL} (default) loads all of them under cleaned names +#' (leading \code{X_} stripped). Pass a character vector of obsm keys to +#' load a subset, optionally named to control the Seurat reduction names: +#' \code{reductions = c(scvi = "X_scVI", umap = "X_umap")} loads only those +#' two keys as reductions \code{scvi} and \code{umap}. Name collisions +#' (e.g. \code{X_pca} and \code{pca} both present) are resolved by keeping +#' the raw obsm key for later claimants, with a warning, instead of the +#' previous silent overwrite. #' #' @return A \code{Seurat} object. If \code{use.bpcells} is set, the count matrix #' is stored on disk in BPCells format and the object uses minimal memory. #' +#' @section Post-read verification: +#' Before returning, the loaded object is verified against the file: dims +#' must match the file's obs/var counts (orientation check), cell names must +#' equal the obs index in order, feature names must equal the var index +#' modulo Seurat's documented underscore-to-dash replacement, and no +#' duplicate barcodes/features may survive. Violations raise +#' \code{scConvert_data_error}. Duplicate names in the file are made unique +#' with a \code{scConvert_names_warning}; a counts layer left holding +#' non-integer values raises a \code{scConvert_counts_warning}. What the +#' reader did (which slot became the counts layer, where X went, the +#' scConvert version) is recorded in \code{misc$scConvert_read}. +#' #' @importFrom hdf5r H5File h5attr h5attr_names #' @importFrom Matrix sparseMatrix #' @importFrom Seurat CreateSeuratObject SetAssayData CreateDimReducObject @@ -121,7 +142,8 @@ #' @export #' readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, - components = NULL, use.c = TRUE, verbose = TRUE) { + components = NULL, use.c = TRUE, verbose = TRUE, + reductions = NULL) { if (!file.exists(file)) { stop("File not found: ", file, call. = FALSE) } @@ -148,7 +170,7 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, if (c_available) { c_result <- tryCatch( .readH5AD_c(file, assay.name = assay.name, components = components, - verbose = verbose), + reductions = reductions, verbose = verbose), error = function(e) { if (verbose) message("C reader failed (", conditionMessage(e), "), using R reader") NULL @@ -167,150 +189,14 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, message("Loading H5AD file: ", file) } - # Helper function to read H5AD sparse or dense matrix - ReadH5ADMatrix <- function(h5_obj, transpose = TRUE) { - if (inherits(h5_obj, "H5Group")) { - # Sparse matrix (CSR or CSC format in h5ad) - if (h5_obj$exists("data") && h5_obj$exists("indices") && h5_obj$exists("indptr")) { - # Pre-read structural sanity check: indptr[length(indptr)] is the - # number of stored non-zeros, which must equal both length(data) and - # length(indices). A truncated /X/data would otherwise read in - # quietly and the downstream sparseMatrix call may segfault inside - # the C-level constructor. Compare HDF5-reported dataset sizes - # before pulling data into memory. - data_d <- h5_obj[["data"]] - indices_d <- h5_obj[["indices"]] - indptr_d <- h5_obj[["indptr"]] - n_data_h5 <- if (!is.null(data_d$dims)) data_d$dims[1] else NA_integer_ - n_indices_h5 <- if (!is.null(indices_d$dims)) indices_d$dims[1] else NA_integer_ - n_indptr_h5 <- if (!is.null(indptr_d$dims)) indptr_d$dims[1] else NA_integer_ - if (!is.na(n_data_h5) && !is.na(n_indices_h5) && - n_data_h5 != n_indices_h5) { - cond <- structure( - class = c("scConvert_data_error", "error", "condition"), - list(message = sprintf( - "Malformed h5ad sparse group: data has %d entries but indices has %d", - n_data_h5, n_indices_h5), - call = NULL)) - stop(cond) - } - # We can also cross-check against indptr[-1] once indptr is read. - data_vals <- h5_obj[["data"]][] - indices <- h5_obj[["indices"]][] # 0-based - indptr <- h5_obj[["indptr"]][] - if (length(indptr) >= 1L) { - declared_nnz <- as.integer(indptr[length(indptr)]) - if (declared_nnz != length(data_vals) || - declared_nnz != length(indices)) { - cond <- structure( - class = c("scConvert_data_error", "error", "condition"), - list(message = sprintf( - "Malformed h5ad sparse group: indptr declares %d nonzeros but data has %d and indices has %d", - declared_nnz, length(data_vals), length(indices)), - call = NULL)) - stop(cond) - } - } - - # Detect encoding type: CSR vs CSC - encoding <- tryCatch(h5attr(h5_obj, "encoding-type"), error = function(e) "csr_matrix") - is_csc <- identical(encoding, "csc_matrix") - - # Get dimensions from shape attribute - if (h5_obj$attr_exists("shape")) { - shape <- h5attr(h5_obj, "shape") - n_rows <- shape[1] - n_cols <- shape[2] - } else if (is_csc) { - n_cols <- length(indptr) - 1L - n_rows <- if (length(indices) > 0) max(indices) + 1L else 0L - } else { - n_rows <- length(indptr) - 1L - n_cols <- if (length(indices) > 0) max(indices) + 1L else 0L - } + # Read h5ad sparse/dense matrices via the shared package-level reader + # (.h5ad_read_matrix in R/VerifyH5AD.R), which this closure used to inline. + ReadH5ADMatrix <- .h5ad_read_matrix - # Helper: sort row indices within each column for valid dgCMatrix - # dgCMatrix requires @i to be increasing within each column (defined by @p). - # scipy CSR/CSC matrices may have unsorted indices (e.g. scanpy pbmc3k raw/X). - .sort_dgc_indices <- function(i, p, x) { - needs_sort <- FALSE - n_col <- length(p) - 1L - for (ci in seq_len(n_col)) { - start <- p[ci] + 1L - end <- p[ci + 1L] - if (end > start && is.unsorted(i[start:end])) { - needs_sort <- TRUE - break - } - } - if (!needs_sort) return(list(i = i, x = x)) - # Sort indices within each column - for (ci in seq_len(n_col)) { - start <- p[ci] + 1L - end <- p[ci + 1L] - if (end > start) { - seg <- start:end - ord <- order(i[seg]) - i[seg] <- i[seg][ord] - x[seg] <- x[seg][ord] - } - } - list(i = i, x = x) - } - - if (is_csc) { - # CSC format: indptr = column pointers, indices = row indices - # dgCMatrix is natively CSC, so construct directly - sorted <- .sort_dgc_indices(as.integer(indices), as.integer(indptr), - as.numeric(data_vals)) - mat <- new("dgCMatrix", - i = sorted$i, - p = as.integer(indptr), - x = sorted$x, - Dim = c(as.integer(n_rows), as.integer(n_cols)) - ) - if (transpose) { - mat <- Matrix::t(mat) - } - } else if (transpose) { - # CSR->CSC reinterpretation: CSR of (n_rows x n_cols) == CSC of (n_cols x n_rows) - # h5ad CSR indices become dgCMatrix @i (row indices in transposed view). - # scipy CSR may have unsorted column indices within rows, so we must - # sort after reinterpretation to produce a valid dgCMatrix. - i_int <- as.integer(indices) - p_int <- as.integer(indptr) - x_num <- as.numeric(data_vals) - sorted <- .sort_dgc_indices(i_int, p_int, x_num) - mat <- new("dgCMatrix", - i = sorted$i, - p = p_int, - x = sorted$x, - Dim = c(as.integer(n_cols), as.integer(n_rows)) - ) - } else { - # Keep original CSR orientation (e.g. obsp graphs) - indices_1based <- indices + 1L - row_indices <- rep(seq_len(n_rows), diff(indptr)) - mat <- sparseMatrix( - i = row_indices, - j = indices_1based, - x = data_vals, - dims = c(n_rows, n_cols), - index1 = TRUE - ) - } - return(mat) - } - } else if (inherits(h5_obj, "H5D")) { - # Dense matrix - mat <- h5_obj[,] - if (transpose) { - mat <- t(mat) - } - return(mat) - } - stop("Unknown matrix format", call. = FALSE) - } + # Resolve where the raw counts live before anything is loaded: the + # /uns/scConvert stamp when present (version-aware branch), else + # /layers/counts (modern convention), /raw/X (legacy scanpy), or /X. + counts_source <- .h5ad_resolve_counts_source(h5ad) # 1. Read cell names if (verbose) message("Reading cell names...") @@ -400,10 +286,22 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, feature.names <- paste0("Gene", seq_len(n_features)) } - # Deduplicate feature names (some datasets have duplicates, e.g. squidpy four_i) - if (anyDuplicated(feature.names)) { + # Deduplicate names (some datasets have duplicates, e.g. squidpy four_i). + # Duplicates are a data-integrity hazard (silent misalignment on any + # name-based join), so the rename is loud (scConvert_names_warning) and + # recorded in misc$scConvert_read. The original var index survives in the + # 'orig_var_index' feature metadata column via .h5ad_preserve_var_identity. + var_index_original <- feature.names + dedup_features <- anyDuplicated(feature.names) > 0L + if (dedup_features) { + .h5ad_warn_duplicate_names("feature", sum(duplicated(feature.names))) feature.names <- make.unique(feature.names) } + dedup_cells <- anyDuplicated(cell.names) > 0L + if (dedup_cells) { + .h5ad_warn_duplicate_names("cell", sum(duplicated(cell.names))) + cell.names <- make.unique(cell.names) + } # 3. Read main expression matrix if (verbose) message("Reading expression matrix...") @@ -418,14 +316,16 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, "Install with: remotes::install_github('bnprks/BPCells')", call. = FALSE) } - # Prefer raw/X (raw counts) over X (often normalized) -- matches non-BPCells path - has_raw <- h5ad$exists("raw") && h5ad[["raw"]]$exists("X") - bp_group <- if (has_raw) "raw/X" else "X" + # Load the resolved counts source on disk -- matches the non-BPCells + # path's convention. Previously only raw/X was considered, so files + # following the modern layers['counts'] convention silently served + # log-normalized X values as on-disk "counts". + bp_group <- counts_source if (verbose) { - if (has_raw) { - message("Loading raw counts (raw/X) via BPCells (on-disk)...") - } else { + if (identical(bp_group, "X")) { message("Loading expression matrix (X) via BPCells (on-disk)...") + } else { + message("Loading raw counts (", bp_group, ") via BPCells (on-disk)...") } } @@ -489,90 +389,24 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, min.features = 0 ) - # 5. Add raw counts if present - # In scanpy convention: X = normalized/processed, raw/X = raw counts - # When raw/X exists, X should become the "data" layer and raw/X the "counts" layer - if (h5ad$exists("raw") && h5ad[["raw"]]$exists("X")) { - if (use_bpcells) { - # In BPCells mode, skip raw/X to preserve on-disk matrix. - # Loading raw/X in-memory would defeat the purpose of on-disk mode. - if (verbose) message("Skipping raw counts (BPCells on-disk mode preserves X as counts)") - } else { - if (verbose) message("Adding raw counts...") - - raw_features <- NULL - if (h5ad[["raw"]]$exists("var")) { - raw_var <- h5ad[["raw/var"]] - if (raw_var$exists("_index")) { - raw_features <- as.character(raw_var[["_index"]][]) - } else if (raw_var$exists("index")) { - raw_features <- as.character(raw_var[["index"]][]) - } - } - - if (!is.null(raw_features)) { - raw_matrix <- ReadH5ADMatrix(h5ad[["raw/X"]], transpose = TRUE) - - # Handle dimension mismatches (dense matrices may need additional transpose) - n_raw_features <- length(raw_features) - if (nrow(raw_matrix) == n_cells && ncol(raw_matrix) == n_raw_features && - nrow(raw_matrix) != n_raw_features) { - raw_matrix <- t(raw_matrix) - } - - # Match dimensions - raw_features <- raw_features[seq_len(min(length(raw_features), nrow(raw_matrix)))] - rownames(raw_matrix) <- raw_features - colnames(raw_matrix) <- cell.names - - # Find common features - common_features <- intersect(feature.names, raw_features) - if (length(common_features) > 0) { - raw_subset <- raw_matrix[common_features, , drop = FALSE] - # raw/X -> counts layer (actual raw counts) - seurat_obj[[assay.name]] <- SetAssayData( - object = seurat_obj[[assay.name]], - layer = "counts", - new.data = raw_subset - ) - # X (already loaded as counts in step 4) -> data layer (normalized) - # expr_matrix contains the X values which are normalized when raw exists - x_subset <- expr_matrix[common_features, , drop = FALSE] - seurat_obj[[assay.name]] <- SetAssayData( - object = seurat_obj[[assay.name]], - layer = "data", - new.data = x_subset - ) - if (verbose) message(" Set raw/X as counts, X as data (normalized)") - } - } - } - } - - # 5b. Ensure the default assay has a "data" layer even in the simple - # single-X case. Seurat 5's CreateSeuratObject(counts=) populates only - # the counts layer; FeaturePlot / FetchData require "data" and will - # fail with "layer 'data' is not found in the object". Populate data - # from counts as a copy (the user is expected to call NormalizeData() - # later when they actually need normalised values). - if (!use_bpcells) { - tryCatch({ - data_layer <- tryCatch( - Seurat::GetAssayData(seurat_obj, assay = assay.name, layer = "data"), - error = function(e) NULL - ) - if (is.null(data_layer) || length(data_layer) == 0L) { - counts_layer <- Seurat::GetAssayData(seurat_obj, assay = assay.name, - layer = "counts") - seurat_obj[[assay.name]] <- SetAssayData( - object = seurat_obj[[assay.name]], - layer = "data", - new.data = counts_layer - ) - } - }, error = function(e) { - if (verbose) message(" Could not copy counts -> data layer: ", e$message) - }) + # 5. Apply the counts/data convention for the resolved counts source + # (.h5ad_apply_counts_convention in R/VerifyH5AD.R): + # raw/X -> counts <- raw/X (common features), data <- X + # layers/counts -> counts <- layers/counts, data <- X + # X -> counts stays X; data <- copy (call NormalizeData() + # when real normalized values are needed) + # In BPCells mode the chosen source is already the on-disk matrix, so the + # in-memory relocation is skipped; when that source is not X itself, X was + # never loaded at all. + if (use_bpcells) { + x_mapped_to <- if (identical(counts_source, "X")) "counts" else "not_loaded" + } else { + counts_conv <- .h5ad_apply_counts_convention( + seurat_obj, h5ad, assay.name, expr_matrix, feature.names, cell.names, + counts_source, verbose = verbose + ) + seurat_obj <- counts_conv$object + x_mapped_to <- counts_conv$x_mapped_to } # 6. Add layers if present @@ -584,6 +418,14 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, layer_names <- names(h5ad[["layers"]]) for (layer_name in layer_names) { + # layers/counts was already placed in the counts layer by the + # counts-convention step (5); re-reading it here would be wasted IO + # and, worse, its failure mode used to be a swallowed message that + # left log-normalized X sitting in the counts slot. + if (identical(layer_name, "counts") && + identical(counts_source, "layers/counts")) { + next + } if (verbose) message(" Adding layer: ", layer_name) # hdf5r reads a dense h5py (cells x features) dataset as an R matrix @@ -618,8 +460,11 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, layer = layer_name)) stop(cond) } - rownames(layer_matrix) <- feature.names - colnames(layer_matrix) <- cell.names + # Label with the object's final dimnames (Seurat may have replaced + # underscores with dashes in feature names; a file-index label would + # make SetAssayData fail and silently drop the layer) + rownames(layer_matrix) <- rownames(seurat_obj) + colnames(layer_matrix) <- colnames(seurat_obj) # Map layer names to Seurat slots seurat_slot <- switch(layer_name, @@ -701,7 +546,9 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, }) } - # Read modern categoricals (groups with codes/categories) + # Read modern categoricals (groups with codes/categories). The stored + # category order becomes the factor level order verbatim, and the + # AnnData `ordered` flag round-trips into an ordered factor. for (col in cat_cols) { tryCatch({ col_obj <- obs_group[[col]] @@ -709,7 +556,10 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, if (encoding_type == "categorical" && col_obj$exists("categories") && col_obj$exists("codes")) { codes <- col_obj[["codes"]]$read() categories <- as.character(col_obj[["categories"]]$read()) - obs_batch[[col]] <- DecodeCategorical(codes, categories) + is_ordered <- tryCatch(isTRUE(as.logical(h5attr(col_obj, "ordered"))[1]), + error = function(e) FALSE) + obs_batch[[col]] <- DecodeCategorical(codes, categories, + ordered = is_ordered) } }, error = function(e) { if (verbose) message("Could not add metadata column '", col, "': ", e$message) @@ -734,11 +584,17 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, obsm_obj <- h5ad[["obsm"]] if (inherits(obsm_obj, "H5D")) { - # Legacy compound dataset: obsm is a structured array with named fields + # Legacy compound dataset: obsm is a structured array with named fields. + # exclude_spatial = FALSE: the modern spatial pipeline (step 12) only + # sees H5Group obsm, so a legacy 'spatial' field must stay a reduction. obsm_df <- obsm_obj$read() + obsm_plan <- .h5ad_plan_reductions(names(obsm_df), + reductions = reductions, + exclude_spatial = FALSE) n_cells <- length(cell.names) - for (reduc_name in names(obsm_df)) { - clean_name <- gsub("^X_", "", reduc_name) + for (plan_i in seq_along(obsm_plan)) { + reduc_name <- obsm_plan[[plan_i]] + clean_name <- names(obsm_plan)[plan_i] if (verbose) message(" Adding reduction: ", clean_name) tryCatch({ vals <- obsm_df[[reduc_name]] @@ -755,11 +611,14 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, }) } } else { - # Modern format: obsm is an H5Group with named datasets - for (reduc_name in names(obsm_obj)) { - clean_name <- gsub("^X_", "", reduc_name) - # Skip 'spatial' -- handled separately in step 12 - if (clean_name == "spatial") next + # Modern format: obsm is an H5Group with named datasets. + # 'spatial' is excluded from the plan -- handled separately in step 12. + obsm_plan <- .h5ad_plan_reductions(names(obsm_obj), + reductions = reductions, + exclude_spatial = TRUE) + for (plan_i in seq_along(obsm_plan)) { + reduc_name <- obsm_plan[[plan_i]] + clean_name <- names(obsm_plan)[plan_i] if (verbose) message(" Adding reduction: ", clean_name) tryCatch({ @@ -825,10 +684,13 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, } if (length(meta_values) == nrow(seurat_obj)) { - names(meta_values) <- feature.names + # Name by the object's final rownames, not the file var index: + # Seurat's underscore-to-dash replacement would otherwise leave + # the names unmatched and the assignment misaligned. + names(meta_values) <- rownames(seurat_obj) seurat_obj[[assay.name]][[col]] <- meta_values if (col == "highly_variable" && is.logical(meta_values)) { - VariableFeatures(seurat_obj) <- feature.names[meta_values] + VariableFeatures(seurat_obj) <- rownames(seurat_obj)[meta_values] } } }, error = function(e) { @@ -857,12 +719,16 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, col_obj <- var_group[[col]] if (inherits(col_obj, "H5Group")) { - # Modern h5ad categorical format + # Modern h5ad categorical format; category order and the `ordered` + # flag both round-trip into the factor. encoding_type <- tryCatch(h5attr(col_obj, "encoding-type"), error = function(e) "") if (encoding_type == "categorical" && col_obj$exists("categories") && col_obj$exists("codes")) { codes <- col_obj[["codes"]]$read() categories <- as.character(col_obj[["categories"]]$read()) - meta_values <- DecodeCategorical(codes, categories) + is_ordered <- tryCatch(isTRUE(as.logical(h5attr(col_obj, "ordered"))[1]), + error = function(e) FALSE) + meta_values <- DecodeCategorical(codes, categories, + ordered = is_ordered) } } else if (inherits(col_obj, "H5D")) { # Check legacy categorical format @@ -887,14 +753,16 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, } } - # Ensure length matches and name with feature names for Seurat v5 compatibility + # Ensure length matches and name with the object's final rownames for + # Seurat v5 compatibility (the file var index may differ after + # Seurat's underscore-to-dash replacement). if (length(meta_values) == nrow(seurat_obj)) { - names(meta_values) <- feature.names + names(meta_values) <- rownames(seurat_obj) seurat_obj[[assay.name]][[col]] <- meta_values # Set variable features if highly_variable column exists if (col == "highly_variable" && is.logical(meta_values)) { - VariableFeatures(seurat_obj) <- feature.names[meta_values] + VariableFeatures(seurat_obj) <- rownames(seurat_obj)[meta_values] } } }, error = function(e) { @@ -904,6 +772,12 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, } # end else (non-compound var) } + # 9b. Gene identity: if the final rownames differ from the file's var index + # (make.unique dedup or Seurat's underscore mangling), keep the original + # identifiers as feature metadata so identity is never silently lost. + seurat_obj <- .h5ad_preserve_var_identity(seurat_obj, assay.name, + var_index_original) + # 10. Add neighbor graphs from obsp. Auxiliary slot: a single corrupt # graph should warn and skip rather than abort the whole load, matching # the varp pattern below. @@ -1018,6 +892,19 @@ readH5AD <- function(file, assay.name = "RNA", use.bpcells = NULL, # object keyed by library name. Closes the CosMx/Xenium silent-loss gap. seurat_obj <- .rebuild_fovs_from_h5ad(h5ad, seurat_obj, verbose = verbose) + # 13. Post-read verification (read-back discipline): assert shape / + # orientation, cell/feature identity and order, and name uniqueness + # against the file rather than trusting the load; surface a counts layer + # that ended up holding non-integer values; record what the reader did. + .h5ad_verify_read(seurat_obj, cell.names, feature.names, file) + if (!use_bpcells) { + .h5ad_warn_noninteger_counts(seurat_obj, assay.name, counts_source) + } + seurat_obj <- .h5ad_record_provenance(seurat_obj, file, counts_source, + x_mapped_to, + dedup_cells = dedup_cells, + dedup_features = dedup_features) + # Store source path for deferred loading (Optimization 4) if (!setequal(components, all_components)) { seurat_obj@misc[[".__h5ad_path__"]] <- normalizePath(file) @@ -1137,7 +1024,8 @@ scLoadMeta <- function(object, components = NULL, verbose = TRUE) { obj } -.readH5AD_c <- function(file, assay.name = "RNA", components = NULL, verbose = TRUE) { +.readH5AD_c <- function(file, assay.name = "RNA", components = NULL, + reductions = NULL, verbose = TRUE) { if (verbose) message("Loading H5AD file (C reader): ", file) # Call C reader for requested components @@ -1145,7 +1033,7 @@ scLoadMeta <- function(object, components = NULL, verbose = TRUE) { if (is.null(result)) { if (verbose) message("C reader failed, falling back to R reader") return(readH5AD(file, assay.name = assay.name, components = components, - use.c = FALSE, verbose = verbose)) + use.c = FALSE, verbose = verbose, reductions = reductions)) } # 1. Construct expression matrix from C result @@ -1168,12 +1056,23 @@ scLoadMeta <- function(object, components = NULL, verbose = TRUE) { ) } - # Set dimnames + # Set dimnames; duplicates are made unique loudly (scConvert_names_warning) + # and recorded in misc$scConvert_read, mirroring the R path. cell.names <- mat_data$colnames feature.names <- mat_data$rownames if (is.null(cell.names)) cell.names <- paste0("Cell", seq_len(ncol(expr_matrix))) if (is.null(feature.names)) feature.names <- paste0("Gene", seq_len(nrow(expr_matrix))) - if (anyDuplicated(feature.names)) feature.names <- make.unique(feature.names) + var_index_original <- feature.names + dedup_features <- anyDuplicated(feature.names) > 0L + if (dedup_features) { + .h5ad_warn_duplicate_names("feature", sum(duplicated(feature.names))) + feature.names <- make.unique(feature.names) + } + dedup_cells <- anyDuplicated(cell.names) > 0L + if (dedup_cells) { + .h5ad_warn_duplicate_names("cell", sum(duplicated(cell.names))) + cell.names <- make.unique(cell.names) + } rownames(expr_matrix) <- feature.names colnames(expr_matrix) <- cell.names @@ -1181,6 +1080,21 @@ scLoadMeta <- function(object, components = NULL, verbose = TRUE) { if (verbose) message("Creating Seurat object...") seurat_obj <- .fast_create_seurat(expr_matrix, assay.name = assay.name) + # 2b. Counts/data convention. The C reader hands back /X only; where the + # raw counts actually live (layers/counts, raw/X, or X itself) is resolved + # against the file and the matrices are relocated accordingly. Runs + # unconditionally -- counts identity is part of the primary-matrix + # contract, not an optional component. + h5ad <- H5File$new(file, mode = "r") + on.exit(tryCatch(h5ad$close_all(), error = function(e) NULL), add = TRUE) + counts_source <- .h5ad_resolve_counts_source(h5ad) + counts_conv <- .h5ad_apply_counts_convention( + seurat_obj, h5ad, assay.name, expr_matrix, feature.names, cell.names, + counts_source, verbose = verbose + ) + seurat_obj <- counts_conv$object + x_mapped_to <- counts_conv$x_mapped_to + # 3. Add obs metadata if ("obs" %in% components && !is.null(result[["obs"]])) { if (verbose) message("Adding cell metadata...") @@ -1209,10 +1123,10 @@ scLoadMeta <- function(object, components = NULL, verbose = TRUE) { else if (is.numeric(meta_values)) meta_values <- as.logical(meta_values) } if (length(meta_values) == nrow(seurat_obj)) { - names(meta_values) <- feature.names + names(meta_values) <- rownames(seurat_obj) seurat_obj[[assay.name]][[col]] <- meta_values if (col == "highly_variable" && is.logical(meta_values)) { - VariableFeatures(seurat_obj) <- feature.names[meta_values] + VariableFeatures(seurat_obj) <- rownames(seurat_obj)[meta_values] } } }, error = function(e) { @@ -1221,12 +1135,23 @@ scLoadMeta <- function(object, components = NULL, verbose = TRUE) { } } - # 5. Add obsm reductions + # 4b. The compiled reader preserves category order but drops the AnnData + # `ordered` flag; restore it from the file's attributes. Also keep the + # original var index when feature names were deduplicated or mangled. + seurat_obj <- .h5ad_restore_ordered_factors(seurat_obj, h5ad, assay.name) + seurat_obj <- .h5ad_preserve_var_identity(seurat_obj, assay.name, + var_index_original) + + # 5. Add obsm reductions ('spatial' excluded from the plan: the spatial + # pipeline in the hdf5r fallback below handles it) if ("obsm" %in% components && !is.null(result[["obsm"]])) { if (verbose) message("Adding dimensional reductions...") - for (reduc_name in names(result[["obsm"]])) { - clean_name <- gsub("^X_", "", reduc_name) - if (clean_name == "spatial") next + obsm_plan <- .h5ad_plan_reductions(names(result[["obsm"]]), + reductions = reductions, + exclude_spatial = TRUE) + for (plan_i in seq_along(obsm_plan)) { + reduc_name <- obsm_plan[[plan_i]] + clean_name <- names(obsm_plan)[plan_i] tryCatch({ embeddings <- result[["obsm"]][[reduc_name]] if (!is.matrix(embeddings)) embeddings <- as.matrix(embeddings) @@ -1278,69 +1203,24 @@ scLoadMeta <- function(object, components = NULL, verbose = TRUE) { } } - # 7. Handle remaining components via hdf5r fallback (varp, layers, uns, spatial, raw) - # These are less performance-critical, so use the existing R reader + # 7. Handle remaining components via the already-open hdf5r handle. + # raw/X and layers/counts were consumed by the counts-convention step + # (2b), which replaced the hand-rolled raw/X reconstruction that used to + # live here (and, unlike it, validates encoding, sorts indices, and covers + # the layers/counts convention). needs_hdf5r <- any(c("varp", "layers", "uns") %in% components) if (needs_hdf5r) { - h5ad <- H5File$new(file, mode = "r") - on.exit(tryCatch(h5ad$close_all(), error = function(e) NULL)) - - # raw/X handling: set counts and data layers - if ("layers" %in% components || "X" %in% components) { - if (h5ad$exists("raw") && h5ad[["raw"]]$exists("X")) { - if (verbose) message("Adding raw counts...") - raw_features <- NULL - if (h5ad[["raw"]]$exists("var")) { - raw_var <- h5ad[["raw/var"]] - if (raw_var$exists("_index")) raw_features <- as.character(raw_var[["_index"]][]) - else if (raw_var$exists("index")) raw_features <- as.character(raw_var[["index"]][]) - } - if (!is.null(raw_features)) { - # Use the same ReadH5ADMatrix helper pattern - raw_obj <- h5ad[["raw/X"]] - if (inherits(raw_obj, "H5Group") && raw_obj$exists("data")) { - data_vals <- raw_obj[["data"]][] - indices <- raw_obj[["indices"]][] - indptr <- raw_obj[["indptr"]][] - encoding <- tryCatch(h5attr(raw_obj, "encoding-type"), error = function(e) "csr_matrix") - if (identical(encoding, "csc_matrix")) { - raw_matrix <- new("dgCMatrix", - i = as.integer(indices), p = as.integer(indptr), x = as.numeric(data_vals), - Dim = c(as.integer(length(raw_features)), as.integer(length(cell.names))) - ) - } else { - # CSR of (n_cells x n_features) == CSC of (n_features x n_cells) - raw_matrix <- new("dgCMatrix", - i = as.integer(indices), p = as.integer(indptr), x = as.numeric(data_vals), - Dim = c(as.integer(length(raw_features)), as.integer(length(cell.names))) - ) - } - } else if (inherits(raw_obj, "H5D")) { - raw_matrix <- t(raw_obj[,]) - } - if (exists("raw_matrix")) { - common <- intersect(feature.names, raw_features) - if (length(common) > 0) { - rownames(raw_matrix) <- raw_features - colnames(raw_matrix) <- cell.names - seurat_obj[[assay.name]] <- SetAssayData( - seurat_obj[[assay.name]], layer = "counts", new.data = raw_matrix[common, , drop = FALSE] - ) - seurat_obj[[assay.name]] <- SetAssayData( - seurat_obj[[assay.name]], layer = "data", new.data = expr_matrix[common, , drop = FALSE] - ) - } - } - } - } - } - # Layers. Mirror the R-reader behaviour: a shape mismatch between a # layer and X is structurally malformed; raise scConvert_data_error # instead of silently dropping the layer. if ("layers" %in% components && h5ad$exists("layers")) { if (verbose) message("Adding layers...") for (layer_name in names(h5ad[["layers"]])) { + # Already placed into the counts layer by the counts-convention step + if (identical(layer_name, "counts") && + identical(counts_source, "layers/counts")) { + next + } layer_obj <- h5ad[["layers"]][[layer_name]] if (inherits(layer_obj, "H5Group") && layer_obj$exists("data")) { ld <- layer_obj[["data"]][]; li <- layer_obj[["indices"]][]; lp <- layer_obj[["indptr"]][] @@ -1363,8 +1243,8 @@ scLoadMeta <- function(object, components = NULL, verbose = TRUE) { stop(cond) } tryCatch({ - rownames(layer_matrix) <- feature.names - colnames(layer_matrix) <- cell.names + rownames(layer_matrix) <- rownames(seurat_obj) + colnames(layer_matrix) <- colnames(seurat_obj) seurat_slot <- switch(layer_name, "counts" = "counts", "data" = "data", "log_normalized" = "data", "scale.data" = "scale.data", @@ -1419,16 +1299,27 @@ scLoadMeta <- function(object, components = NULL, verbose = TRUE) { } } - # Spatial - if ("obsm" %in% components && h5ad$exists("obsm") && "spatial" %in% names(h5ad[["obsm"]])) { - seurat_obj <- H5ADSpatialToSeurat(h5ad_file = h5ad, seurat_obj = seurat_obj, - assay_name = assay.name, verbose = verbose) - } - # FOV rebuild (mirrors the non-components path above) seurat_obj <- .rebuild_fovs_from_h5ad(h5ad, seurat_obj, verbose = verbose) } + # Spatial: gated on obsm alone. The handle is open regardless of + # needs_hdf5r now, so a components subset like c("X", "obs", "obsm") no + # longer silently drops the spatial image data. + if ("obsm" %in% components && h5ad$exists("obsm") && + "spatial" %in% names(h5ad[["obsm"]])) { + seurat_obj <- H5ADSpatialToSeurat(h5ad_file = h5ad, seurat_obj = seurat_obj, + assay_name = assay.name, verbose = verbose) + } + + # Post-read verification + provenance (mirrors the R reader) + .h5ad_verify_read(seurat_obj, cell.names, feature.names, file) + .h5ad_warn_noninteger_counts(seurat_obj, assay.name, counts_source) + seurat_obj <- .h5ad_record_provenance(seurat_obj, file, counts_source, + x_mapped_to, + dedup_cells = dedup_cells, + dedup_features = dedup_features) + # Store source path for deferred loading all_components <- c("X", "obs", "var", "obsm", "obsp", "varp", "layers", "uns") if (!setequal(components, all_components)) { diff --git a/R/VerifyH5AD.R b/R/VerifyH5AD.R new file mode 100644 index 0000000..e89fea8 --- /dev/null +++ b/R/VerifyH5AD.R @@ -0,0 +1,754 @@ +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +# Reverse-conversion (h5ad -> Seurat) integrity helpers +# +# readH5AD() historically trusted its own bookkeeping: whatever /X held was +# loaded into the Seurat counts layer, obsm keys were cleaned into reduction +# names with silent last-writer-wins collisions, and callers had no way to +# confirm that the returned object matched the file. The helpers here close +# those gaps: +# +# * .h5ad_resolve_counts_source() version/layout-aware detection of where +# the raw counts live (uns stamp, +# /layers/counts, /raw/X, /X) +# * .h5ad_apply_counts_convention() single place that maps the resolved +# source onto Seurat counts/data layers +# * .h5ad_warn_noninteger_counts() loud, classed warning when the counts +# layer ends up holding normalized values +# * .h5ad_plan_reductions() obsm -> reduction-name plan with +# collision handling and user remapping +# * .h5ad_verify_read() post-read shape/orientation and +# name-identity assertions +# * .h5ad_record_provenance() what-the-reader-did record in misc +# * .h5ad_stamp_provenance() writer-side /uns/scConvert stamp that +# future reads can branch on +#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +#' Read an h5ad-encoded matrix (sparse CSR/CSC group or dense dataset) +#' +#' Package-level version of the reader previously defined as a closure inside +#' \code{readH5AD}, so the counts-convention helper and the C-path fallback can +#' share one implementation. \code{transpose = TRUE} returns the Seurat +#' orientation (features x cells) for an h5ad (cells x genes) matrix. +#' +#' @param h5_obj hdf5r H5Group (sparse) or H5D (dense) +#' @param transpose Transpose into features x cells orientation +#' @return dgCMatrix or dense matrix +#' @keywords internal +#' @noRd +.h5ad_read_matrix <- function(h5_obj, transpose = TRUE) { + if (inherits(h5_obj, "H5Group")) { + # Sparse matrix (CSR or CSC format in h5ad) + if (h5_obj$exists("data") && h5_obj$exists("indices") && h5_obj$exists("indptr")) { + # Pre-read structural sanity check: indptr[length(indptr)] is the + # number of stored non-zeros, which must equal both length(data) and + # length(indices). A truncated /X/data would otherwise read in + # quietly and the downstream sparseMatrix call may segfault inside + # the C-level constructor. Compare HDF5-reported dataset sizes + # before pulling data into memory. + data_d <- h5_obj[["data"]] + indices_d <- h5_obj[["indices"]] + indptr_d <- h5_obj[["indptr"]] + n_data_h5 <- if (!is.null(data_d$dims)) data_d$dims[1] else NA_integer_ + n_indices_h5 <- if (!is.null(indices_d$dims)) indices_d$dims[1] else NA_integer_ + n_indptr_h5 <- if (!is.null(indptr_d$dims)) indptr_d$dims[1] else NA_integer_ + if (!is.na(n_data_h5) && !is.na(n_indices_h5) && + n_data_h5 != n_indices_h5) { + cond <- structure( + class = c("scConvert_data_error", "error", "condition"), + list(message = sprintf( + "Malformed h5ad sparse group: data has %d entries but indices has %d", + n_data_h5, n_indices_h5), + call = NULL)) + stop(cond) + } + # We can also cross-check against indptr[-1] once indptr is read. + data_vals <- h5_obj[["data"]][] + indices <- h5_obj[["indices"]][] # 0-based + indptr <- h5_obj[["indptr"]][] + if (length(indptr) >= 1L) { + declared_nnz <- as.integer(indptr[length(indptr)]) + if (declared_nnz != length(data_vals) || + declared_nnz != length(indices)) { + cond <- structure( + class = c("scConvert_data_error", "error", "condition"), + list(message = sprintf( + "Malformed h5ad sparse group: indptr declares %d nonzeros but data has %d and indices has %d", + declared_nnz, length(data_vals), length(indices)), + call = NULL)) + stop(cond) + } + } + + # Detect encoding type: CSR vs CSC + encoding <- tryCatch(h5attr(h5_obj, "encoding-type"), error = function(e) "csr_matrix") + is_csc <- identical(encoding, "csc_matrix") + + # Get dimensions from shape attribute + if (h5_obj$attr_exists("shape")) { + shape <- h5attr(h5_obj, "shape") + n_rows <- shape[1] + n_cols <- shape[2] + } else if (is_csc) { + n_cols <- length(indptr) - 1L + n_rows <- if (length(indices) > 0) max(indices) + 1L else 0L + } else { + n_rows <- length(indptr) - 1L + n_cols <- if (length(indices) > 0) max(indices) + 1L else 0L + } + + # Helper: sort row indices within each column for valid dgCMatrix + # dgCMatrix requires @i to be increasing within each column (defined by @p). + # scipy CSR/CSC matrices may have unsorted indices (e.g. scanpy pbmc3k raw/X). + .sort_dgc_indices <- function(i, p, x) { + needs_sort <- FALSE + n_col <- length(p) - 1L + for (ci in seq_len(n_col)) { + start <- p[ci] + 1L + end <- p[ci + 1L] + if (end > start && is.unsorted(i[start:end])) { + needs_sort <- TRUE + break + } + } + if (!needs_sort) return(list(i = i, x = x)) + # Sort indices within each column + for (ci in seq_len(n_col)) { + start <- p[ci] + 1L + end <- p[ci + 1L] + if (end > start) { + seg <- start:end + ord <- order(i[seg]) + i[seg] <- i[seg][ord] + x[seg] <- x[seg][ord] + } + } + list(i = i, x = x) + } + + if (is_csc) { + # CSC format: indptr = column pointers, indices = row indices + # dgCMatrix is natively CSC, so construct directly + sorted <- .sort_dgc_indices(as.integer(indices), as.integer(indptr), + as.numeric(data_vals)) + mat <- new("dgCMatrix", + i = sorted$i, + p = as.integer(indptr), + x = sorted$x, + Dim = c(as.integer(n_rows), as.integer(n_cols)) + ) + if (transpose) { + mat <- Matrix::t(mat) + } + } else if (transpose) { + # CSR->CSC reinterpretation: CSR of (n_rows x n_cols) == CSC of (n_cols x n_rows) + # h5ad CSR indices become dgCMatrix @i (row indices in transposed view). + # scipy CSR may have unsorted column indices within rows, so we must + # sort after reinterpretation to produce a valid dgCMatrix. + i_int <- as.integer(indices) + p_int <- as.integer(indptr) + x_num <- as.numeric(data_vals) + sorted <- .sort_dgc_indices(i_int, p_int, x_num) + mat <- new("dgCMatrix", + i = sorted$i, + p = p_int, + x = sorted$x, + Dim = c(as.integer(n_cols), as.integer(n_rows)) + ) + } else { + # Keep original CSR orientation (e.g. obsp graphs) + indices_1based <- indices + 1L + row_indices <- rep(seq_len(n_rows), diff(indptr)) + mat <- sparseMatrix( + i = row_indices, + j = indices_1based, + x = data_vals, + dims = c(n_rows, n_cols), + index1 = TRUE + ) + } + return(mat) + } + } else if (inherits(h5_obj, "H5D")) { + # Dense matrix + mat <- h5_obj[,] + if (transpose) { + mat <- t(mat) + } + return(mat) + } + stop("Unknown matrix format", call. = FALSE) +} + +# Slots an h5ad counts stamp may legally point at. +.h5ad_counts_locations <- c("layers/counts", "raw/X", "X") + +#' Does a nested h5 path ("layers/counts") exist in an open handle? +#' @keywords internal +#' @noRd +.h5ad_slot_exists <- function(h5ad, path) { + parts <- strsplit(path, "/", fixed = TRUE)[[1]] + node <- h5ad + for (p in parts) { + ok <- tryCatch(node$exists(p), error = function(e) FALSE) + if (!isTRUE(ok)) return(FALSE) + node <- tryCatch(node[[p]], error = function(e) NULL) + if (is.null(node)) return(FALSE) + } + TRUE +} + +#' Decide which slot of an h5ad file holds the raw counts +#' +#' Version/layout-aware resolution, in priority order: +#' \enumerate{ +#' \item \code{/uns/scConvert/counts_location} -- stamp written by +#' scConvert's own h5ad writers recording where counts were placed. +#' Files written by older scConvert (or other tools) carry no stamp and +#' fall through to the structural heuristics. +#' \item \code{/layers/counts} -- modern anndata/scanpy convention +#' (X = normalized data). +#' \item \code{/raw/X} -- legacy scanpy convention. +#' \item \code{/X} -- no dedicated counts slot; X is all there is. +#' } +#' +#' @param h5ad Open hdf5r H5File handle +#' @return character(1): "layers/counts", "raw/X", or "X" +#' @keywords internal +#' @noRd +.h5ad_resolve_counts_source <- function(h5ad) { + stamped <- tryCatch({ + if (.h5ad_slot_exists(h5ad, "uns/scConvert/counts_location")) { + as.character(h5ad[["uns"]][["scConvert"]][["counts_location"]]$read())[1] + } else { + NULL + } + }, error = function(e) NULL) + if (!is.null(stamped) && stamped %in% .h5ad_counts_locations && + .h5ad_slot_exists(h5ad, stamped)) { + return(stamped) + } + # A stamp of "none" (data-only export) still resolves to X below: the + # non-integer counts warning is what surfaces the missing-counts state. + if (.h5ad_slot_exists(h5ad, "layers/counts")) return("layers/counts") + if (.h5ad_slot_exists(h5ad, "raw/X")) return("raw/X") + "X" +} + +#' Move X / raw/X / layers/counts into the Seurat layers the convention implies +#' +#' \code{CreateSeuratObject()} is always fed /X, whatever /X holds. This step +#' relocates matrices so the Seurat layers match the resolved counts source: +#' \itemize{ +#' \item \code{"raw/X"}: counts <- raw/X (common features), data <- X +#' \item \code{"layers/counts"}: counts <- layers/counts, data <- X +#' \item \code{"X"}: counts stays X; data <- copy of counts +#' } +#' Runs on every load, not gated on \code{components}: counts identity is part +#' of the primary-matrix contract, exactly like the raw/X handling always was. +#' Failure to relocate X out of the counts layer is the silent +#' lognorm-in-counts failure this step exists to prevent, so the +#' layers/counts branch raises \code{scConvert_data_error} instead of +#' downgrading to a message. +#' +#' @param seurat_obj Seurat object freshly built from /X +#' @param h5ad Open hdf5r H5File handle +#' @param assay.name Assay to populate +#' @param expr_matrix The matrix read from /X (features x cells, dimnames set) +#' @param feature.names Final (deduplicated) feature names +#' @param cell.names Final (deduplicated) cell names +#' @param counts_source Result of \code{.h5ad_resolve_counts_source()} +#' @param verbose Emit progress messages +#' @return list(object = , x_mapped_to = "counts"|"data") +#' @keywords internal +#' @noRd +.h5ad_apply_counts_convention <- function(seurat_obj, h5ad, assay.name, + expr_matrix, feature.names, cell.names, + counts_source, verbose = TRUE) { + x_mapped_to <- "counts" + + # Label relocation matrices with the object's final dimnames, not the file + # index: Seurat's underscore-to-dash replacement would otherwise make + # SetAssayData fail with "no feature overlap" (or, in the pre-fix layers + # loop, silently skip -- leaving log-normalized X in the counts slot for + # any file with underscore gene names). Positional relabeling is safe: + # the row/col order is the file order, which .h5ad_verify_read asserts. + object_features <- rownames(seurat_obj) + object_cells <- colnames(seurat_obj) + x_relabeled <- expr_matrix + if (nrow(x_relabeled) == length(object_features) && + ncol(x_relabeled) == length(object_cells)) { + rownames(x_relabeled) <- object_features + colnames(x_relabeled) <- object_cells + } + + if (identical(counts_source, "raw/X")) { + if (verbose) message("Adding raw counts...") + + raw_features <- NULL + if (h5ad[["raw"]]$exists("var")) { + raw_var <- h5ad[["raw/var"]] + if (raw_var$exists("_index")) { + raw_features <- as.character(raw_var[["_index"]][]) + } else if (raw_var$exists("index")) { + raw_features <- as.character(raw_var[["index"]][]) + } + } + + if (!is.null(raw_features)) { + raw_matrix <- .h5ad_read_matrix(h5ad[["raw/X"]], transpose = TRUE) + + # Handle dimension mismatches (dense matrices may need additional transpose) + n_raw_features <- length(raw_features) + n_cells <- length(cell.names) + if (nrow(raw_matrix) == n_cells && ncol(raw_matrix) == n_raw_features && + nrow(raw_matrix) != n_raw_features) { + raw_matrix <- t(raw_matrix) + } + + # Match dimensions + raw_features <- raw_features[seq_len(min(length(raw_features), nrow(raw_matrix)))] + # Same Seurat name transformation the object went through, so the + # intersection happens in the object's namespace. + rownames(raw_matrix) <- gsub("_", "-", raw_features, fixed = TRUE) + colnames(raw_matrix) <- object_cells + + # Find common features + common_features <- intersect(object_features, rownames(raw_matrix)) + if (length(common_features) > 0) { + # raw/X -> counts layer (actual raw counts) + seurat_obj[[assay.name]] <- SetAssayData( + object = seurat_obj[[assay.name]], + layer = "counts", + new.data = raw_matrix[common_features, , drop = FALSE] + ) + # X (loaded as counts at creation) -> data layer (normalized) + seurat_obj[[assay.name]] <- SetAssayData( + object = seurat_obj[[assay.name]], + layer = "data", + new.data = x_relabeled[common_features, , drop = FALSE] + ) + x_mapped_to <- "data" + if (verbose) message(" Set raw/X as counts, X as data (normalized)") + } + } + } else if (identical(counts_source, "layers/counts")) { + if (verbose) message("Adding counts from layers/counts...") + + lyr_obj <- h5ad[["layers"]][["counts"]] + # Dense layers arrive from hdf5r already shaped (features x cells); + # sparse groups need the transpose into Seurat orientation. Mirrors the + # generic layers loop in readH5AD(). + counts_matrix <- if (inherits(lyr_obj, "H5D")) { + lyr_obj[, ] + } else { + .h5ad_read_matrix(lyr_obj, transpose = TRUE) + } + + if (nrow(counts_matrix) != nrow(expr_matrix) || + ncol(counts_matrix) != ncol(expr_matrix)) { + stop(.scconvert_data_error(sprintf( + "Malformed h5ad: layer '/layers/counts' has shape (%d x %d) but X has shape (%d x %d)", + nrow(counts_matrix), ncol(counts_matrix), + nrow(expr_matrix), ncol(expr_matrix)), + layer = "counts")) + } + rownames(counts_matrix) <- object_features + colnames(counts_matrix) <- object_cells + + # X first moves into data (it currently sits in the counts layer), then + # the real counts replace it. + seurat_obj[[assay.name]] <- SetAssayData( + object = seurat_obj[[assay.name]], + layer = "data", + new.data = x_relabeled + ) + seurat_obj[[assay.name]] <- SetAssayData( + object = seurat_obj[[assay.name]], + layer = "counts", + new.data = counts_matrix + ) + x_mapped_to <- "data" + if (verbose) message(" Set layers/counts as counts, X as data (normalized)") + } + + # Safety net for every source: guarantee a "data" layer exists. Seurat 5's + # CreateSeuratObject(counts=) populates only the counts layer; FeaturePlot / + # FetchData require "data" and fail with "layer 'data' is not found" without + # it. In the counts_source == "X" case this is the documented copy (the user + # is expected to call NormalizeData() when they need normalized values). + tryCatch({ + # suppressWarnings: probing a missing data layer makes Seurat 5 warn + # "Layer 'data' is empty" -- emptiness is exactly what is being tested. + data_layer <- tryCatch( + suppressWarnings( + Seurat::GetAssayData(seurat_obj, assay = assay.name, layer = "data") + ), + error = function(e) NULL + ) + if (is.null(data_layer) || length(data_layer) == 0L) { + counts_layer <- Seurat::GetAssayData(seurat_obj, assay = assay.name, + layer = "counts") + seurat_obj[[assay.name]] <- SetAssayData( + object = seurat_obj[[assay.name]], + layer = "data", + new.data = counts_layer + ) + } + }, error = function(e) { + if (verbose) message(" Could not copy counts -> data layer: ", e$message) + }) + + list(object = seurat_obj, x_mapped_to = x_mapped_to) +} + +#' Warn (classed) when the counts layer holds non-integer values +#' +#' The silent failure mode behind counts-layer ambiguity: a normalized /X +#' read into the Seurat counts slot breaks every downstream step that assumes +#' integer counts, without any signal. This check makes it loud. The warning +#' carries class \code{scConvert_counts_warning} so pipelines can +#' \code{withCallingHandlers()} on it specifically. +#' +#' @param object Seurat object after counts-convention application +#' @param assay.name Assay to check +#' @param counts_source Where the counts layer came from +#' @param max_check Cap on the number of values inspected +#' @return invisible(TRUE) if a warning was raised, invisible(FALSE) otherwise +#' @keywords internal +#' @noRd +.h5ad_warn_noninteger_counts <- function(object, assay.name, counts_source, + max_check = 1e6L) { + vals <- tryCatch({ + m <- Seurat::GetAssayData(object, assay = assay.name, layer = "counts") + if (inherits(m, "dgCMatrix")) m@x else as.numeric(m) + }, error = function(e) NULL) + if (is.null(vals) || length(vals) == 0L) return(invisible(FALSE)) + if (length(vals) > max_check) vals <- vals[seq_len(max_check)] + fractional <- vals != floor(vals) + if (!isTRUE(any(fractional, na.rm = TRUE))) return(invisible(FALSE)) + + msg <- if (identical(counts_source, "X")) { + paste0( + "readH5AD: the counts layer holds non-integer values taken from /X. ", + "The file has no dedicated counts slot (/layers/counts or /raw/X), so ", + "/X was used as counts even though it appears to contain normalized ", + "data. Downstream steps that assume integer counts should not trust ", + "this layer.") + } else { + sprintf(paste0( + "readH5AD: the counts layer (read from /%s) contains non-integer ", + "values; that slot may hold normalized data rather than raw counts."), + counts_source) + } + warning(warningCondition(msg, class = "scConvert_counts_warning")) + invisible(TRUE) +} + +#' Emit a classed warning for name deduplication +#' @keywords internal +#' @noRd +.h5ad_warn_duplicate_names <- function(what, n_dup) { + extra <- if (identical(what, "feature")) { + " The original index is preserved in the 'orig_var_index' feature metadata column." + } else { + "" + } + warning(warningCondition(sprintf( + "readH5AD: %d duplicate %s name(s) in the file were made unique with make.unique().%s", + n_dup, what, extra), class = "scConvert_names_warning")) +} + +#' Plan which obsm keys become which Seurat reductions +#' +#' Replaces the old inline \code{gsub("^X_", "", key)} mapping, which silently +#' overwrote reductions when two obsm keys cleaned to the same name (e.g. +#' \code{X_pca} and \code{pca}). Also implements the user-facing remapping +#' interface: \code{reductions = c(scvi = "X_scVI")} loads only \code{X_scVI} +#' and names it \code{scvi} in the Seurat object. +#' +#' @param obsm_keys Character vector of keys present under /obsm +#' @param reductions NULL (load all), or a character vector of obsm keys to +#' load; names, when present, set the Seurat reduction names +#' @param exclude_spatial Drop spatial keys from the default plan (they are +#' handled by the spatial/image pipeline instead) +#' @return Named character vector: values = obsm keys, names = reduction names +#' @keywords internal +#' @noRd +.h5ad_plan_reductions <- function(obsm_keys, reductions = NULL, + exclude_spatial = TRUE) { + if (is.null(reductions)) { + keys <- obsm_keys + if (exclude_spatial) { + keys <- keys[gsub("^X_", "", keys) != "spatial"] + } + targets <- gsub("^X_", "", keys) + } else { + if (!is.character(reductions) || length(reductions) == 0L) { + stop("`reductions` must be NULL or a character vector of obsm keys ", + "(optionally named with the Seurat reduction names to use)", + call. = FALSE) + } + missing_keys <- setdiff(reductions, obsm_keys) + if (length(missing_keys) > 0L) { + warning(warningCondition(sprintf( + "readH5AD: requested obsm key(s) not present in file: %s. Available: %s", + paste(missing_keys, collapse = ", "), + paste(obsm_keys, collapse = ", ")), + class = "scConvert_reduction_warning")) + } + keep <- reductions %in% obsm_keys + keys <- unname(reductions[keep]) + targets <- if (is.null(names(reductions))) { + rep("", length(keys)) + } else { + names(reductions)[keep] + } + targets[!nzchar(targets)] <- gsub("^X_", "", keys[!nzchar(targets)]) + } + if (length(keys) == 0L) { + return(stats::setNames(character(0), character(0))) + } + + dup <- duplicated(targets) + if (any(dup)) { + collided <- unique(targets[dup]) + # First claimant keeps the clean name; later ones fall back to their raw + # obsm key, then make.unique as a last resort. + targets[dup] <- keys[dup] + if (anyDuplicated(targets)) targets <- make.unique(targets) + warning(warningCondition(sprintf( + "readH5AD: multiple obsm keys map to the same reduction name (%s); kept all as: %s", + paste(collided, collapse = ", "), + paste(sprintf("%s -> %s", keys, targets), collapse = ", ")), + class = "scConvert_reduction_warning")) + } + stats::setNames(keys, targets) +} + +#' Restore the `ordered` flag on categorical metadata (C-reader path) +#' +#' The compiled reader builds plain factors (category order preserved, ordered +#' flag dropped). Re-read just the cheap `ordered` attributes via hdf5r and +#' upgrade the matching columns. +#' +#' @param seurat_obj Seurat object with metadata already attached +#' @param h5ad Open hdf5r H5File handle +#' @param assay.name Assay carrying feature metadata +#' @return The Seurat object with ordered factors restored +#' @keywords internal +#' @noRd +.h5ad_restore_ordered_factors <- function(seurat_obj, h5ad, assay.name) { + ordered_cols <- function(grp_name) { + out <- character(0) + if (!isTRUE(tryCatch(h5ad$exists(grp_name), error = function(e) FALSE))) { + return(out) + } + grp <- h5ad[[grp_name]] + if (!inherits(grp, "H5Group")) return(out) + for (nm in names(grp)) { + child <- tryCatch(grp[[nm]], error = function(e) NULL) + if (is.null(child) || !inherits(child, "H5Group")) next + enc <- tryCatch(h5attr(child, "encoding-type"), error = function(e) "") + if (!identical(enc, "categorical")) next + ord <- tryCatch(isTRUE(as.logical(h5attr(child, "ordered"))[1]), + error = function(e) FALSE) + if (ord) out <- c(out, nm) + } + out + } + + for (col in ordered_cols("obs")) { + v <- seurat_obj@meta.data[[col]] + if (is.factor(v) && !is.ordered(v)) { + seurat_obj@meta.data[[col]] <- factor(v, levels = levels(v), ordered = TRUE) + } + } + for (col in ordered_cols("var")) { + tryCatch({ + fmeta <- seurat_obj[[assay.name]][[]] + v <- fmeta[[col]] + if (is.factor(v) && !is.ordered(v)) { + v <- factor(v, levels = levels(v), ordered = TRUE) + names(v) <- rownames(seurat_obj[[assay.name]]) + seurat_obj[[assay.name]][[col]] <- v + } + }, error = function(e) NULL) + } + seurat_obj +} + +#' Preserve the file's original var index when feature names had to change +#' +#' Feature names can drift from the file's var index through +#' \code{make.unique()} deduplication or Seurat's documented underscore-to-dash +#' replacement. When that happens the original identifiers (often the only +#' stable gene identity, e.g. Ensembl IDs) are kept as feature-level metadata +#' column \code{orig_var_index}, so gene identity is never silently lost and +#' round-trips can re-join on it. +#' +#' @param seurat_obj Seurat object with final feature names +#' @param assay.name Assay to annotate +#' @param var_index_original var index exactly as stored in the file +#' @return The Seurat object, possibly with an orig_var_index column added +#' @keywords internal +#' @noRd +.h5ad_preserve_var_identity <- function(seurat_obj, assay.name, + var_index_original) { + final_names <- rownames(seurat_obj) + if (length(var_index_original) != length(final_names)) return(seurat_obj) + if (identical(final_names, var_index_original)) return(seurat_obj) + existing <- tryCatch(colnames(seurat_obj[[assay.name]][[]]), + error = function(e) character(0)) + if ("orig_var_index" %in% existing) return(seurat_obj) + tryCatch({ + vals <- var_index_original + names(vals) <- final_names + seurat_obj[[assay.name]][["orig_var_index"]] <- vals + }, error = function(e) NULL) + seurat_obj +} + +#' Post-read structural verification +#' +#' Read-back discipline for the reverse conversion: instead of trusting that +#' the load succeeded, assert what can be asserted. Failures raise +#' \code{scConvert_data_error}. +#' \itemize{ +#' \item shape/orientation: object dims must equal the file's (n_var, n_obs) +#' \item cell identity/order: colnames must be exactly the obs index +#' \item feature identity/order: rownames must be the var index modulo +#' Seurat's documented underscore-to-dash replacement (and the +#' make.unique that can follow it) +#' \item uniqueness: no duplicate barcodes or genes may survive +#' } +#' +#' @param object Loaded Seurat object +#' @param cell.names Expected cell names (file obs index, post-dedup) +#' @param feature.names Expected feature names (file var index, post-dedup) +#' @param file Source path, for error messages +#' @return invisible(TRUE); raises scConvert_data_error otherwise +#' @keywords internal +#' @noRd +.h5ad_verify_read <- function(object, cell.names, feature.names, file) { + n_obs <- length(cell.names) + n_var <- length(feature.names) + if (ncol(object) != n_obs || nrow(object) != n_var) { + stop(.scconvert_data_error(sprintf( + paste0("Post-read verification failed for '%s': loaded object is ", + "%d features x %d cells but the file declares %d var x %d obs. ", + "This indicates a transpose/orientation fault or a malformed file."), + file, nrow(object), ncol(object), n_var, n_obs))) + } + if (!identical(colnames(object), cell.names)) { + stop(.scconvert_data_error(sprintf( + paste0("Post-read verification failed for '%s': cell names/order in the ", + "loaded object do not match the file's obs index."), + file))) + } + actual_features <- rownames(object) + mangled <- gsub("_", "-", feature.names, fixed = TRUE) + features_ok <- identical(actual_features, feature.names) || + identical(actual_features, mangled) || + identical(actual_features, make.unique(mangled)) + if (!features_ok) { + stop(.scconvert_data_error(sprintf( + paste0("Post-read verification failed for '%s': feature names/order in ", + "the loaded object do not match the file's var index."), + file))) + } + if (anyDuplicated(colnames(object)) > 0L) { + stop(.scconvert_data_error(sprintf( + "Post-read verification failed for '%s': duplicate cell barcodes in loaded object.", + file))) + } + if (anyDuplicated(actual_features) > 0L) { + stop(.scconvert_data_error(sprintf( + "Post-read verification failed for '%s': duplicate feature names in loaded object.", + file))) + } + invisible(TRUE) +} + +#' Record what the reader actually did in misc$scConvert_read +#' +#' So callers (and any read-back verification downstream) can branch on facts +#' rather than trusting \code{readH5AD()}'s success return: which slot became +#' the counts layer, where X went, whether names had to be de-duplicated, and +#' which scConvert version performed the read. +#' +#' @keywords internal +#' @noRd +.h5ad_record_provenance <- function(object, file, counts_source, x_mapped_to, + dedup_cells = FALSE, dedup_features = FALSE) { + object@misc[["scConvert_read"]] <- list( + package_version = as.character(utils::packageVersion("scConvert")), + source_file = tryCatch(normalizePath(file), error = function(e) file), + source_format = "h5ad", + counts_source = counts_source, + x_mapped_to = x_mapped_to, + n_obs = ncol(object), + n_var = nrow(object), + duplicate_cell_names_renamed = isTRUE(dedup_cells), + duplicate_feature_names_renamed = isTRUE(dedup_features), + verified = TRUE + ) + object +} + +#' Stamp /uns/scConvert provenance into a freshly-written h5ad file +#' +#' Records the writing package version and -- by inspecting the file itself, +#' not the writer's intent -- where the raw counts landed +#' (\code{layers/counts}, \code{raw/X}, \code{X}, or \code{none}). Readers use +#' the stamp as the authoritative, version-aware branch for counts resolution, +#' eliminating the layout guessing that made older files ambiguous. +#' +#' Best-effort: returns invisible(FALSE) on any failure rather than failing +#' the write. +#' +#' @param target An open hdf5r H5File in write mode, or a file path +#' @return invisible(TRUE) on success +#' @keywords internal +#' @noRd +.h5ad_stamp_provenance <- function(target) { + own_handle <- is.character(target) + h5 <- if (own_handle) { + tryCatch(hdf5r::H5File$new(target, mode = "r+"), error = function(e) NULL) + } else { + target + } + if (is.null(h5)) return(invisible(FALSE)) + + ok <- tryCatch({ + counts_location <- if (.h5ad_slot_exists(h5, "layers/counts")) { + "layers/counts" + } else if (.h5ad_slot_exists(h5, "raw/X")) { + "raw/X" + } else if (isTRUE(h5$exists("X"))) { + "X" + } else { + "none" + } + + if (!h5$exists("uns")) { + uns <- h5$create_group("uns") + AddAnndataEncoding(uns, encoding_type = "dict", encoding_version = "0.1.0") + } + uns <- h5[["uns"]] + if (uns$exists("scConvert")) uns$link_delete("scConvert") + sg <- uns$create_group("scConvert") + AddAnndataEncoding(sg, encoding_type = "dict", encoding_version = "0.1.0") + # Scalar strings are written the same way WriteUnsItem writes them (no + # element encoding attr): anndata's legacy fallback reads them cleanly, + # which the existing python-validation CI already exercises. + version_str <- as.character(utils::packageVersion("scConvert")) + sg$create_dataset("version", robj = version_str, + dtype = CachedUtf8Type(), chunk_dims = 1L) + sg$create_dataset("counts_location", robj = counts_location, + dtype = CachedUtf8Type(), chunk_dims = 1L) + TRUE + }, error = function(e) FALSE) + + if (own_handle) tryCatch(h5$close_all(), error = function(e) NULL) + invisible(isTRUE(ok)) +} diff --git a/R/WriteH5AD.R b/R/WriteH5AD.R index e10510e..af76139 100644 --- a/R/WriteH5AD.R +++ b/R/WriteH5AD.R @@ -515,8 +515,10 @@ DirectSeuratToH5AD <- function( if (verbose) message(" Writing uns...") uns_grp <- dfile$create_group("uns") - # Internal keys managed separately (varp, lazy-load bookkeeping) - skip_keys <- c("__varp__", ".__h5ad_path__", ".__h5ad_loaded__") + # Internal keys managed separately (varp, lazy-load bookkeeping, and the + # reader-side provenance record, which describes a past read, not this file) + skip_keys <- c("__varp__", ".__h5ad_path__", ".__h5ad_loaded__", + "scConvert_read") skip_keys <- c(skip_keys, grep("^__varp__\\.", names(misc), value = TRUE)) gzip <- GetCompressionLevel() @@ -560,6 +562,13 @@ DirectSeuratToH5AD <- function( dfile$create_attr(attr_name = 'encoding-version', robj = '0.1.0', dtype = CachedGuessDType('0.1.0'), space = ScalarSpace()) + # ========== Provenance stamp ========== + # /uns/scConvert/{version, counts_location}: records -- by inspecting the + # file just written, not the writer's intent -- where the raw counts + # landed, so readers resolve the counts layer by version-aware fact + # instead of layout heuristics. + .h5ad_stamp_provenance(dfile) + dfile$flush() dfile$close_all() diff --git a/R/zzz.R b/R/zzz.R index 3b2942c..07132a6 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -1133,12 +1133,16 @@ SafeSetLayerData <- function(object, layer, value) { # Append misc and spatial data via R if present. misc <- tryCatch(Misc(object), error = function(e) list()) images <- tryCatch(Images(object), error = function(e) character(0)) - skip_keys <- c("__varp__", ".__h5ad_path__", ".__h5ad_loaded__") + skip_keys <- c("__varp__", ".__h5ad_path__", ".__h5ad_loaded__", + "scConvert_read") skip_keys <- c(skip_keys, grep("^__varp__\\.", names(misc), value = TRUE)) misc_to_write <- misc[!names(misc) %in% skip_keys] misc_to_write <- misc_to_write[!vapply(misc_to_write, is.null, logical(1))] has_varp <- !is.null(tryCatch(Misc(object)[["__varp__"]], error = function(e) NULL)) - if (length(misc_to_write) > 0 || length(images) > 0 || has_varp) { + if (!(length(misc_to_write) > 0 || length(images) > 0 || has_varp)) { + # Nothing to append; stamp provenance by path. + .h5ad_stamp_provenance(filename) + } else { h5 <- hdf5r::H5File$new(filename, mode = "r+") on.exit(h5$close_all(), add = TRUE) if (!h5$exists("uns")) h5$create_group("uns") @@ -1177,6 +1181,9 @@ SafeSetLayerData <- function(object, layer, value) { }, error = function(e) NULL) } } + # Stamp last so /uns/scConvert reflects the final file (a stale + # misc$scConvert echoed back through WriteUnsItem is replaced). + .h5ad_stamp_provenance(h5) } return(invisible(filename)) } diff --git a/man/DecodeCategorical.Rd b/man/DecodeCategorical.Rd index 759c6f1..2b9551a 100644 --- a/man/DecodeCategorical.Rd +++ b/man/DecodeCategorical.Rd @@ -4,17 +4,24 @@ \alias{DecodeCategorical} \title{Decode AnnData categorical encoding to R factor} \usage{ -DecodeCategorical(codes, categories) +DecodeCategorical(codes, categories, ordered = FALSE) } \arguments{ \item{codes}{Integer vector of 0-based category codes (-1 = NA)} \item{categories}{Character vector of category labels} + +\item{ordered}{Logical; produce an ordered factor (AnnData's +\code{ordered} categorical flag). Default \code{FALSE}.} } \value{ A factor vector } \description{ -Decode AnnData categorical encoding to R factor +The factor levels follow the stored category order verbatim (pandas +preserves an explicit category order; re-sorting it alphabetically would +silently reorder positional palettes and any order-dependent downstream +logic). When the AnnData categorical carries \code{ordered = TRUE}, pass +it here to get an ordered factor back. } \keyword{internal} diff --git a/man/readH5AD.Rd b/man/readH5AD.Rd index 125c9a7..9016462 100644 --- a/man/readH5AD.Rd +++ b/man/readH5AD.Rd @@ -10,7 +10,8 @@ readH5AD( use.bpcells = NULL, components = NULL, use.c = TRUE, - verbose = TRUE + verbose = TRUE, + reductions = NULL ) } \arguments{ @@ -33,6 +34,16 @@ loading via \code{\link{scLoadMeta}}.} FALSE to force the pure-R hdf5r path.} \item{verbose}{Show progress messages} + +\item{reductions}{Which \code{obsm} entries to load as dimensional +reductions. \code{NULL} (default) loads all of them under cleaned names +(leading \code{X_} stripped). Pass a character vector of obsm keys to +load a subset, optionally named to control the Seurat reduction names: +\code{reductions = c(scvi = "X_scVI", umap = "X_umap")} loads only those +two keys as reductions \code{scvi} and \code{umap}. Name collisions +(e.g. \code{X_pca} and \code{pca} both present) are resolved by keeping +the raw obsm key for later claimants, with a warning, instead of the +previous silent overwrite.} } \value{ A \code{Seurat} object. If \code{use.bpcells} is set, the count matrix @@ -44,3 +55,17 @@ Supports optional BPCells on-disk matrix loading for large datasets that exceed available memory. When compiled C routines are available, uses a fast native reader (typically 2-3x faster than the pure-R path). } +\section{Post-read verification}{ + +Before returning, the loaded object is verified against the file: dims +must match the file's obs/var counts (orientation check), cell names must +equal the obs index in order, feature names must equal the var index +modulo Seurat's documented underscore-to-dash replacement, and no +duplicate barcodes/features may survive. Violations raise +\code{scConvert_data_error}. Duplicate names in the file are made unique +with a \code{scConvert_names_warning}; a counts layer left holding +non-integer values raises a \code{scConvert_counts_warning}. What the +reader did (which slot became the counts layer, where X went, the +scConvert version) is recorded in \code{misc$scConvert_read}. +} + diff --git a/tests/testthat/test-reverse-conversion.R b/tests/testthat/test-reverse-conversion.R new file mode 100644 index 0000000..9e0275d --- /dev/null +++ b/tests/testthat/test-reverse-conversion.R @@ -0,0 +1,422 @@ +# Reverse-conversion (h5ad -> Seurat) integrity fixes: +# * version/layout-aware counts resolution (layers/counts, raw/X, X, +# /uns/scConvert stamp) +# * non-integer counts warning +# * post-read shape/orientation + name-identity verification +# * duplicate barcode/gene handling +# * reduction key collision handling + remapping +# * categorical level order + ordered flag +# * gene identity preservation (orig_var_index) +# * read/write provenance records + +library(scConvert) + +# ---- helpers ----------------------------------------------------------------- + +scalar_str_attr <- function(obj, name, val) { + invisible(obj$create_attr(name, robj = val, + dtype = hdf5r::H5T_STRING$new(size = Inf), + space = hdf5r::H5S$new(type = "scalar"))) +} + +write_csr <- function(parent, name, m) { + # m is a dense (cells x genes) matrix; store as h5ad CSR + csr <- as(Matrix::t(as(m, "CsparseMatrix")), "CsparseMatrix") + g <- parent$create_group(name) + g$create_dataset("data", robj = as.numeric(csr@x)) + g$create_dataset("indices", robj = as.integer(csr@i)) + g$create_dataset("indptr", robj = as.integer(csr@p)) + invisible(g$create_attr("shape", robj = as.integer(c(nrow(m), ncol(m))))) + scalar_str_attr(g, "encoding-type", "csr_matrix") + invisible(g) +} + +# Build a synthetic h5ad; returns the ground-truth matrices/names +make_h5ad <- function(tmp, n_cells = 12, n_genes = 6, lognorm_x = TRUE, + layers_counts = TRUE, obsm = NULL, obs_cat = NULL, + stamp = NULL, gene_names = NULL, cell_names = NULL) { + set.seed(7) + counts <- matrix(rpois(n_cells * n_genes, 4), nrow = n_cells) + xmat <- if (lognorm_x) log1p(counts / rowSums(counts) * 1e4) else counts + if (is.null(cell_names)) cell_names <- paste0("Cell-", seq_len(n_cells)) + if (is.null(gene_names)) gene_names <- paste0("Gene-", seq_len(n_genes)) + h5 <- hdf5r::H5File$new(tmp, mode = "w") + write_csr(h5, "X", xmat) + if (layers_counts) { + write_csr(h5$create_group("layers"), "counts", counts) + } + obs <- h5$create_group("obs") + obs$create_dataset("_index", robj = cell_names) + scalar_str_attr(obs, "encoding-type", "dataframe") + if (!is.null(obs_cat)) { + cg <- obs$create_group("celltype") + cg$create_dataset("codes", robj = obs_cat$codes) + cg$create_dataset("categories", robj = obs_cat$categories) + scalar_str_attr(cg, "encoding-type", "categorical") + invisible(cg$create_attr("ordered", robj = isTRUE(obs_cat$ordered))) + # anndata always records column-order; the C reader enumerates via it + invisible(obs$create_attr("column-order", robj = "celltype")) + } + var <- h5$create_group("var") + var$create_dataset("_index", robj = gene_names) + scalar_str_attr(var, "encoding-type", "dataframe") + var$create_dataset("gene_ids", robj = paste0("ENSG", seq_len(n_genes))) + invisible(var$create_attr("column-order", robj = "gene_ids")) + if (!is.null(obsm)) { + og <- h5$create_group("obsm") + for (nm in names(obsm)) og$create_dataset(nm, robj = obsm[[nm]]) + } + if (!is.null(stamp)) { + ug <- h5$create_group("uns") + sg <- ug$create_group("scConvert") + sg$create_dataset("counts_location", robj = stamp) + sg$create_dataset("version", robj = "0.3.0") + } + h5$close_all() + invisible(list(counts = counts, x = xmat, genes = gene_names, + cells = cell_names)) +} + +# ---- counts-layer resolution (item: counts-layer ambiguity) ------------------ + +test_that("layers/counts becomes the counts layer and X becomes data (both readers)", { + skip_if_not_installed("hdf5r") + skip_if_not_installed("Seurat") + + tmp <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp), add = TRUE) + info <- make_h5ad(tmp, lognorm_x = TRUE, layers_counts = TRUE) + + for (usec in c(FALSE, TRUE)) { + obj <- suppressWarnings(readH5AD(tmp, verbose = FALSE, use.c = usec)) + cts <- as.matrix(Seurat::GetAssayData(obj, layer = "counts")) + dat <- as.matrix(Seurat::GetAssayData(obj, layer = "data")) + # counts layer holds the integer counts from /layers/counts + expect_true(all(cts == round(cts)), info = paste("use.c =", usec)) + expect_true(all(abs(t(cts) - info$counts) < 1e-8), info = paste("use.c =", usec)) + # data layer holds the (log-normalized) /X values + expect_true(all(abs(t(dat) - info$x) < 1e-6), info = paste("use.c =", usec)) + # provenance records what happened + prov <- obj@misc$scConvert_read + expect_identical(prov$counts_source, "layers/counts") + expect_identical(prov$x_mapped_to, "data") + expect_true(isTRUE(prov$verified)) + } +}) + +test_that("layers/counts is honored even when `components` excludes layers", { + skip_if_not_installed("hdf5r") + skip_if_not_installed("Seurat") + + tmp <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp), add = TRUE) + info <- make_h5ad(tmp, lognorm_x = TRUE, layers_counts = TRUE) + + for (usec in c(FALSE, TRUE)) { + obj <- suppressWarnings( + readH5AD(tmp, verbose = FALSE, use.c = usec, components = c("X", "obs")) + ) + cts <- as.matrix(Seurat::GetAssayData(obj, layer = "counts")) + expect_true(all(cts == round(cts)), info = paste("use.c =", usec)) + expect_true(all(abs(t(cts) - info$counts) < 1e-8), info = paste("use.c =", usec)) + } +}) + +test_that("non-integer counts raise a classed scConvert_counts_warning", { + skip_if_not_installed("hdf5r") + skip_if_not_installed("Seurat") + + tmp <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp), add = TRUE) + make_h5ad(tmp, lognorm_x = TRUE, layers_counts = FALSE) + + expect_warning( + readH5AD(tmp, verbose = FALSE, use.c = FALSE), + class = "scConvert_counts_warning" + ) + + # Integer-valued X without a dedicated counts slot stays quiet + tmp2 <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp2), add = TRUE) + make_h5ad(tmp2, lognorm_x = FALSE, layers_counts = FALSE) + got <- character(0) + withCallingHandlers( + readH5AD(tmp2, verbose = FALSE, use.c = FALSE), + warning = function(w) { + got <<- c(got, class(w)[1]) + invokeRestart("muffleWarning") + } + ) + expect_false("scConvert_counts_warning" %in% got) +}) + +test_that("the /uns/scConvert counts_location stamp overrides layout heuristics", { + skip_if_not_installed("hdf5r") + skip_if_not_installed("Seurat") + + tmp <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp), add = TRUE) + # File has layers/counts, but the stamp pins counts to X: the stamp wins. + make_h5ad(tmp, lognorm_x = TRUE, layers_counts = TRUE, stamp = "X") + + obj <- suppressWarnings(readH5AD(tmp, verbose = FALSE, use.c = FALSE)) + expect_identical(obj@misc$scConvert_read$counts_source, "X") +}) + +# ---- post-read verification (items: transpose correctness, read-back) -------- + +test_that(".h5ad_verify_read catches orientation, identity, and duplicate faults", { + skip_if_not_installed("Seurat") + + m <- matrix(rpois(50, 3), nrow = 5, ncol = 10, + dimnames = list(paste0("g", 1:5), paste0("c", 1:10))) + obj <- suppressWarnings(Seurat::CreateSeuratObject(as(m, "CsparseMatrix"))) + + # Correct expectation passes + expect_true(scConvert:::.h5ad_verify_read(obj, paste0("c", 1:10), + paste0("g", 1:5), "f.h5ad")) + # Swapped dims = transpose/orientation fault + err <- tryCatch( + scConvert:::.h5ad_verify_read(obj, paste0("c", 1:5), paste0("g", 1:10), "f.h5ad"), + error = function(e) e + ) + expect_s3_class(err, "scConvert_data_error") + expect_match(conditionMessage(err), "transpose/orientation") + + # Same dims, different cell order = identity fault + err2 <- tryCatch( + scConvert:::.h5ad_verify_read(obj, rev(paste0("c", 1:10)), + paste0("g", 1:5), "f.h5ad"), + error = function(e) e + ) + expect_s3_class(err2, "scConvert_data_error") + expect_match(conditionMessage(err2), "cell names") + + # Underscore mangling is tolerated (documented Seurat behavior) + expect_true(scConvert:::.h5ad_verify_read(obj, paste0("c", 1:10), + paste0("g", 1:5), "f.h5ad")) + m2 <- m + rownames(m2) <- paste0("g_", 1:5) + obj2 <- suppressWarnings(Seurat::CreateSeuratObject(as(m2, "CsparseMatrix"))) + expect_true(scConvert:::.h5ad_verify_read(obj2, paste0("c", 1:10), + paste0("g_", 1:5), "f.h5ad")) +}) + +# ---- duplicate names (items: cell order/identity, read-back verification) ---- + +test_that("duplicate feature names warn, dedupe, and preserve original identity", { + skip_if_not_installed("hdf5r") + skip_if_not_installed("Seurat") + + tmp <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp), add = TRUE) + info <- make_h5ad(tmp, gene_names = c("G_1", "G_1", "G_2", "G_3", "G_4", "G_5")) + + for (usec in c(FALSE, TRUE)) { + got <- character(0) + obj <- withCallingHandlers( + readH5AD(tmp, verbose = FALSE, use.c = usec), + warning = function(w) { + got <<- c(got, class(w)[1]) + invokeRestart("muffleWarning") + } + ) + expect_true("scConvert_names_warning" %in% got, info = paste("use.c =", usec)) + expect_equal(anyDuplicated(rownames(obj)), 0L) + # counts relocation still lands despite underscore->dash mangling + cts <- as.matrix(Seurat::GetAssayData(obj, layer = "counts")) + expect_true(all(abs(t(cts) - info$counts) < 1e-8), info = paste("use.c =", usec)) + # original var index preserved as feature metadata + fmeta <- obj[["RNA"]][[]] + expect_true("orig_var_index" %in% colnames(fmeta)) + expect_identical(as.character(fmeta$orig_var_index), info$genes) + # var columns still positionally aligned after mangling + expect_identical(as.character(fmeta$gene_ids), paste0("ENSG", 1:6)) + expect_true(isTRUE(obj@misc$scConvert_read$duplicate_feature_names_renamed)) + } +}) + +test_that("duplicate cell barcodes warn, dedupe, and are recorded", { + skip_if_not_installed("hdf5r") + skip_if_not_installed("Seurat") + + tmp <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp), add = TRUE) + cells <- c("AAA", "AAA", paste0("C", 3:12)) + make_h5ad(tmp, cell_names = cells) + + got <- character(0) + obj <- withCallingHandlers( + readH5AD(tmp, verbose = FALSE, use.c = FALSE), + warning = function(w) { + got <<- c(got, class(w)[1]) + invokeRestart("muffleWarning") + } + ) + expect_true("scConvert_names_warning" %in% got) + expect_equal(anyDuplicated(colnames(obj)), 0L) + expect_identical(colnames(obj), make.unique(cells)) + expect_true(isTRUE(obj@misc$scConvert_read$duplicate_cell_names_renamed)) +}) + +# ---- reduction keys (item: reduction key collisions) ------------------------- + +test_that("colliding obsm keys are kept under distinct names with a warning", { + skip_if_not_installed("hdf5r") + skip_if_not_installed("Seurat") + + tmp <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp), add = TRUE) + make_h5ad(tmp, obsm = list(X_pca = matrix(rnorm(36), 12, 3), + pca = matrix(rnorm(24), 12, 2))) + + for (usec in c(FALSE, TRUE)) { + got <- character(0) + obj <- withCallingHandlers( + readH5AD(tmp, verbose = FALSE, use.c = usec), + warning = function(w) { + got <<- c(got, class(w)[1]) + invokeRestart("muffleWarning") + } + ) + expect_true("scConvert_reduction_warning" %in% got, info = paste("use.c =", usec)) + # Both embeddings survive (previously the second silently overwrote the first) + expect_length(names(obj@reductions), 2L) + } +}) + +test_that("`reductions` selects and renames obsm keys", { + skip_if_not_installed("hdf5r") + skip_if_not_installed("Seurat") + + tmp <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp), add = TRUE) + make_h5ad(tmp, obsm = list(X_pca = matrix(rnorm(36), 12, 3), + X_umap = matrix(rnorm(24), 12, 2))) + + for (usec in c(FALSE, TRUE)) { + obj <- suppressWarnings( + readH5AD(tmp, verbose = FALSE, use.c = usec, + reductions = c(scvi = "X_pca")) + ) + expect_identical(names(obj@reductions), "scvi") + expect_equal(ncol(Seurat::Embeddings(obj, "scvi")), 3L) + } + + # Requesting a missing key warns with the reduction class + expect_warning( + suppressWarnings( + readH5AD(tmp, verbose = FALSE, use.c = FALSE, + reductions = c("X_pca", "X_missing")), + classes = "updatedKeyWarning" + ), + class = "scConvert_reduction_warning" + ) +}) + +test_that(".h5ad_plan_reductions maps, excludes spatial, and resolves collisions", { + plan <- scConvert:::.h5ad_plan_reductions(c("X_pca", "X_umap", "spatial")) + expect_identical(unname(plan), c("X_pca", "X_umap")) + expect_identical(names(plan), c("pca", "umap")) + + plan2 <- scConvert:::.h5ad_plan_reductions(c("X_pca", "X_umap", "spatial"), + exclude_spatial = FALSE) + expect_true("spatial" %in% names(plan2)) + + expect_warning( + plan3 <- scConvert:::.h5ad_plan_reductions(c("X_pca", "pca")), + class = "scConvert_reduction_warning" + ) + expect_equal(anyDuplicated(names(plan3)), 0L) + expect_length(plan3, 2L) + + plan4 <- scConvert:::.h5ad_plan_reductions(c("X_pca", "X_umap"), + reductions = c(lat = "X_pca")) + expect_identical(plan4, stats::setNames("X_pca", "lat")) +}) + +# ---- categorical order + ordered flag (item: factor/category level order) ---- + +test_that("categorical level order and the ordered flag survive the read", { + skip_if_not_installed("hdf5r") + skip_if_not_installed("Seurat") + + tmp <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp), add = TRUE) + # Deliberately non-alphabetical category order + cats <- c("zeta", "alpha", "mid") + make_h5ad(tmp, obs_cat = list(codes = as.integer(rep(c(0, 1, 2), 4)), + categories = cats, ordered = TRUE)) + + for (usec in c(FALSE, TRUE)) { + obj <- suppressWarnings(readH5AD(tmp, verbose = FALSE, use.c = usec)) + ct <- obj[["celltype"]][, 1] + expect_s3_class(ct, "factor") + expect_identical(levels(ct), cats, info = paste("use.c =", usec)) + expect_true(is.ordered(ct), info = paste("use.c =", usec)) + } + + # ordered = FALSE round-trips as an unordered factor + tmp2 <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp2), add = TRUE) + make_h5ad(tmp2, obs_cat = list(codes = as.integer(rep(c(0, 1, 2), 4)), + categories = cats, ordered = FALSE)) + obj2 <- suppressWarnings(readH5AD(tmp2, verbose = FALSE, use.c = FALSE)) + ct2 <- obj2[["celltype"]][, 1] + expect_identical(levels(ct2), cats) + expect_false(is.ordered(ct2)) +}) + +test_that("DecodeCategorical preserves category order and honors ordered=", { + f <- scConvert:::DecodeCategorical(c(0L, 2L, 1L), c("z", "a", "m")) + expect_identical(levels(f), c("z", "a", "m")) + expect_false(is.ordered(f)) + f2 <- scConvert:::DecodeCategorical(c(0L, 2L, 1L), c("z", "a", "m"), + ordered = TRUE) + expect_true(is.ordered(f2)) + expect_identical(as.character(f2), c("z", "m", "a")) +}) + +# ---- writer provenance stamp (items: counts ambiguity, version pinning) ------ + +test_that("writeH5AD stamps /uns/scConvert and the round trip restores layers", { + skip_if_not_installed("hdf5r") + skip_if_not_installed("Seurat") + + set.seed(11) + m <- matrix(rpois(60, 3), 10, 6, + dimnames = list(paste0("G", 1:10), paste0("C", 1:6))) + sobj <- suppressWarnings( + Seurat::CreateSeuratObject(counts = as(m, "CsparseMatrix")) + ) + sobj <- suppressWarnings(Seurat::NormalizeData(sobj, verbose = FALSE)) + + tmp <- tempfile(fileext = ".h5ad") + on.exit(unlink(tmp), add = TRUE) + suppressWarnings(suppressMessages(writeH5AD(sobj, tmp, verbose = FALSE))) + + # Stamp present and correct: counts live in layers/counts + h5 <- hdf5r::H5File$new(tmp, mode = "r") + expect_true(h5$exists("uns")) + expect_true(h5[["uns"]]$exists("scConvert")) + expect_identical( + as.character(h5[["uns"]][["scConvert"]][["counts_location"]]$read()), + "layers/counts" + ) + expect_identical( + as.character(h5[["uns"]][["scConvert"]][["version"]]$read()), + as.character(utils::packageVersion("scConvert")) + ) + h5$close_all() + + # Round trip: integer counts back in counts, log-normalized values in data + back <- suppressWarnings(readH5AD(tmp, verbose = FALSE)) + cts <- as.matrix(Seurat::GetAssayData(back, layer = "counts")) + dat <- as.matrix(Seurat::GetAssayData(back, layer = "data")) + expect_true(all(cts == round(cts))) + expect_equal(cts, as.matrix(Seurat::GetAssayData(sobj, layer = "counts")), + ignore_attr = TRUE, tolerance = 1e-8) + expect_equal(dat, as.matrix(Seurat::GetAssayData(sobj, layer = "data")), + ignore_attr = TRUE, tolerance = 1e-6) + expect_identical(back@misc$scConvert_read$counts_source, "layers/counts") +})