Write single-cell data as a Cellucid-ready AnnData (or MuData) store

cellucid_prepare() writes a Seurat, SingleCellExperiment or SummarizedExperiment object, or embeddings, cell and gene metadata, gene expression, a neighbour graph, and vector fields as a Cellucid-ready AnnData: a standard AnnData store, written as AnnData Zarr (format 2, or 3 with zarr_format = 3), that Python's anndata reads without a warning and the Cellucid web app opens directly, with the dataset's name, citation, and launch defaults in uns['cellucid'].

Usage

cellucid_prepare(
  out_path,
  format,
  data = NULL,
  embeddings = NULL,
  dataset_id = NULL,
  dataset_name = NULL,
  obs = NULL,
  var = NULL,
  gene_expression = NULL,
  connectivities = NULL,
  vector_fields = NULL,
  modalities = NULL,
  gene_identifiers = NULL,
  dataset_description = NULL,
  citation_text = NULL,
  citation_url = NULL,
  source_name = NULL,
  source_url = NULL,
  embedding_defaults = NULL,
  latent_key = NULL,
  obs_keys = NULL,
  var_gene_id_column = NULL,
  vector_field_default = NULL,
  normalize_embeddings = NULL,
  centroid_min_points = NULL,
  modality = NULL,
  priority_seed = 0,
  zarr_format = NULL,
  force = FALSE
)

Arguments

out_path: Required path of the store to write, ending in the suffix of format: a folder whose name ends in ".zarr" for format = "zarr", a file whose name ends in ".zarr.zip" for format = "zarr.zip". A path that does not match format is refused with the path to write instead. Missing parent folders are created.
format: Required kind of output, with no default, because each suits a different place and only you know where the store goes: "zarr" writes a .zarr folder of small files, the way to publish on GitHub (no file over 100 MB, no Git LFS); "zarr.zip" writes the same store as one .zarr.zip file, for everywhere else: single-file hosts (Zenodo, Figshare, Hugging Face, S3, a lab server), local files, cellucid serve and Jupyter. See the section Which format?. Matched as match.arg() matches. An .h5ad file is written by Python's cellucid.prepare(format = "h5ad"); this package writes no HDF5, so it depends on no HDF5 library.
data: Optional single-cell object to write: a Seurat object made by Seurat 3.0 or later (read with SeuratObject 4 or 5; Assay and Assay5 assays), a SingleCellExperiment, or a SummarizedExperiment. Every other argument then defaults to what the object holds, and an argument given replaces that one component; see the section Writing an object. With several assays (Seurat) or alternative experiments (SingleCellExperiment) the store is a Cellucid-ready MuData. SeuratObject and SingleCellExperiment are needed only for their own objects.
embeddings: Named list mapping each obsm key to a numeric matrix with one row per cell, such as list(X_umap = umap, X_pca = pca). A numeric data.frame, a Matrix matrix, or a plain numeric vector (one column) is accepted too. Each one is stored whole, every column, as a float32 array at its own coordinates; the viewer chooses the columns it plots and scales them to its view. Every value must be within the float32 range. A row holding NA or NaN (or an infinity) is a cell the embedding does not place (one computed on some of the cells): it is written as given, the viewer leaves it out of that embedding's views, and in the latent space it gets no outlier quantile; at least one cell must have finite coordinates. Each key becomes a path in the store, see the note on names below.
dataset_id: Portable identifier of the dataset, stored as uns['cellucid']['id']: 1–180 ASCII letters, numbers, ., _, or -, beginning with a letter or number, not ending with ., and not a Windows device name. Sessions and community annotations are keyed on it.
dataset_name: Non-empty name of the dataset, shown to the reader verbatim and stored as uns['cellucid']['name'] and uns['title']; see the note on displayed text below.
obs: Optional data.frame of cell metadata, one row per cell. Every column is stored, see the note on column encodings below.
var: Optional data.frame of gene metadata, one row per column of gene_expression. Every column is stored. Supplying it without gene_expression is an error.
gene_expression: Optional cells x genes matrix: a Matrix sparse matrix (stored as a gene-major csc_matrix), a dense numeric matrix (stored as a dense array in tall, gene-major chunks), or a matrix held outside R's memory – a BPCells matrix or a DelayedArray backend such as HDF5Array – read block by block (see Details) and stored as a csc_matrix. Values must be finite and within the float32 range, and are stored as float32 – or the whole matrix as float64 when every value is a whole number and one exceeds 2^24 in magnitude: integer counts float32 cannot hold exactly. That is decided on the genes stored, after gene_identifiers selects them. A stored zero of a sparse matrix is no value and is dropped; the matrix stays sparse. Rows are always cells; Cellucid never transposes the matrix, so transpose a genes x cells matrix with t() first.
connectivities: Optional cells x cells neighbour graph, a Matrix sparse matrix or a numeric matrix, stored as obsp['connectivities'] (a float32 csr_matrix). It must be exactly symmetric, with finite, non-negative weights; a stored zero is not an edge and is dropped, and so is the diagonal (a self-loop is no edge: Seurat's RNA_snn has 1 there). A directed graph (Seurat's RNA_nn) is refused.
vector_fields: Optional per-cell vectors drawn over an embedding, keyed by field and then by the embedding each array belongs to: list(velocity = list(X_umap = v2, X_umap_3d = v3, X_branches = vb)). Each array has one row per cell and one column per column of its embedding, and is stored under scVelo's name for it, obsm['<field>_<basis>'] for the embedding X_<basis> (velocity_umap, velocity_umap_3d, velocity_branches), which the viewer lists as the field '<field>' and draws over that embedding. The embedding's key must therefore start with X_, and each stored key must read back as its own field and embedding: a field named like an embedding (X_...), or a key another embedding's longer basis claims (field T_fwd on X_umap beside an embedding X_fwd_umap), is refused.
modalities: Optional character vector naming the modalities of data to write, in this order: Seurat assays or SingleCellExperiment experiments ("main", or the main experiment's name, and each altExpNames()). NULL writes every one, the Seurat active assay or the main experiment first. One modality is written as a plain AnnData.
gene_identifiers: Optional character vector of the genes to store, in the order given. NULL stores every gene. The identifiers are matched against the gene names the viewer shows (var_gene_id_column, or the var row names) after repeated names were made unique, and must be distinct.
dataset_description: Optional description of the dataset, shown verbatim; stored as uns['cellucid']['description'].
citation_text: Optional text that says how to cite the dataset, stored as uns['cellucid']['citation']['text']. Required whenever citation_url is supplied.
citation_url: Optional absolute http:// or https:// address of the publication, stored as uns['cellucid']['citation']['url'].
source_name: Optional name of where the data comes from, stored as uns['cellucid']['source']['name']. Required whenever source_url is supplied.
source_url: Optional absolute http:// or https:// address of the source, stored as uns['cellucid']['source']['url'].
embedding_defaults: Optional launch default: the embedding each view dimension shows first, by "1d", "2d", or "3d". A value is an obsm key, plotted by its first columns, such as c("2d" = "X_umap"), or list(key = , columns = ) naming the 1-based columns to plot, such as list("3d" = list(key = "X_pca", columns = 1:3)). Stored as uns['cellucid']['defaults']['embedding<d>d'], columns 0-based. Without it the viewer chooses.
latent_key: Optional embeddings key of the latent space the viewer's outlier filter measures in, such as "X_pca" or "X_scVI": a latent space such as scVI or PCA, never a 2D or 3D layout. The outlier quantiles of every categorical obs column (each cell's distance to its category's centre, as the share of its category at that distance or closer) are computed here and stored in uns['cellucid']['outlierQuantiles'], so the viewer filters outliers at once instead of reading the whole latent space to compute them. See the section Outlier quantiles. Also stored as the launch default defaults['latentKey'].
obs_keys: Optional launch default: the obs columns the viewer shows, in this order. Every column is stored either way. Stored as defaults['obsKeys'].
var_gene_id_column: Optional launch default: the var column whose values name the genes in the viewer, such as gene symbols beside Ensembl row names. NULL: the var row names name them; Data mapping in the viewer can choose a column later. A column of text or a factor names genes as it is, an integer column by its digits; a double or logical column cannot name genes and is refused. Stored as defaults['geneIdColumn'].
vector_field_default: Optional launch default: the field of vector_fields selected first, by the name used there ("velocity"), which is the id the viewer lists. One of its embeddings must be shown when the file opens: one named in embedding_defaults or found by the viewer's default rule, or any other 1 to 3 column embedding that is not latent_key. Stored as defaults['vectorFieldDefault'].
normalize_embeddings: Optional launch default: FALSE shows embeddings at their own coordinates instead of scaled to the view. Stored as defaults['normalizeEmbeddings'].
centroid_min_points: Optional positive integer: the fewest cells a category needs to get outlier quantiles (10 when NULL). Stored as the launch default defaults['centroidMinPoints'], and with latent_key as uns['cellucid']['outlierQuantiles']['minPoints'].
modality: Optional launch default for a multimodal store: the modality the viewer opens with. Stored as defaults['modality']. Without it the viewer opens the first modality named like gene expression (rna, gex, gene expression, gene_expression, sct, any letter case), else the first. var_gene_id_column then names a column of that modality's var. gene_expression, var and gene_identifiers describe one feature space and are refused for an object of several: change its assay first, or write one with modalities.
priority_seed: The seed of the cells' priority order, a whole number in [0, 2^32 - 1]. The cells are stored by ascending priority – the "uniform" rank of cellucid_subsample_cells() with this seed – so every subsample Cellucid opens with this seed is the first rows of every array and reads nothing of the cells it leaves out. Each row's original position is kept in obs[['cellucid_original_position']]; cellucid_restore_original_order() puts the cells back in the order given. NULL stores the cells in the order given, without a priority. An object in a Cellucid-ready file's stored order (its cellucid$priority record describes its cells) is put back in its original order first, whatever the seed, so every cell keeps its original position; a cellucid_original_position column with no record that describes the cells names no order and is not written (with a warning). Python's cellucid.prepare(priority_seed = ) takes the same values.
zarr_format: Optional Zarr format of the store: 2 (the default, when NULL), read by every anndata release since 0.8 and by zarr-python 2 and 3, or 3, read by anndata 0.11 and later (zarr-python 3). The web app reads both equally fast. Python's cellucid.prepare(zarr_format = ) takes the same values.
force: TRUE replaces an existing store at out_path, once the new one is complete. An existing folder that is not a Zarr store is never replaced.

Value

Invisibly, the absolute path of the written store.

Details

The store is an AnnData Zarr group (Zarr format 2 unless zarr_format = 3):

X: gene expression, gene-major float32 (float64 for integer counts above 2^24), a csc_matrix with int32 indices and indptr and no stored zeros, or a dense array;
obs and var: data frames, the index in _index;
obsm: every embedding and vector field, float32;
obsp/connectivities: the neighbour graph, float32 csr_matrix;
uns/title and uns/cellucid: schema (1), id, name, and, when supplied, description, citation (text, url), source (name, url), defaults holding only the launch defaults that were set, with latent_key, outlierQuantiles, and in priority order priority (seed, specVersion, nObs – the cell count it was written with, which a reader verifies it by – levels);
in priority order (the default, priority_seed): every per-cell element – obs, obsm, X's rows, both axes of obsp – by ascending priority, obs[['cellucid_original_position']] (uint32) holding each row's original position, and beside a sparse X varm[['cellucid_prefix_counts']] (uint32, one column per level of levels): each gene's values in the rows below that level;
the consolidated metadata (.zmetadata; format 3: the root zarr.json), so a store on a web host opens with one request.

Every element carries the encoding-type and encoding-version attributes anndata 0.8 and later write. Every chunk of values wider than one byte is byte-shuffled, so the similar high bytes of neighbouring values compress together – except the expression values themselves, which compress best as they are – and every chunk is compressed with DEFLATE level 6, which browsers decompress natively: in format 2 numcodecs' Shuffle filter and the compressor {"id": "zlib", "level": 6}; in format 3 the blosc codec (byte shuffle, zlib level 6) for a shuffled array and the gzip codec for every other, only codecs of the Zarr v3 specification (numcodecs.shuffle is not one). A chunk is the unit of every read: every array with one row per cell (data frame columns and index, embeddings as a chunk of rows with every column, masks, cell maps, outlier quantiles) holds 65,536 (2^16) values per chunk, so a subsample – the first rows of a store in priority order – reads little beyond its own cells; arrays read whole (var, the expression matrix's pointers, the neighbour graph) up to 262,144 (2^18); and the data and indices of the expression matrix, which the viewer reads one gene (or its first values) at a time, 16,384 (2^14). Python's cellucid.prepare() writes the same store, element for element, from the same input.

A .zarr.zip holds the same files, byte for byte, as entries of one ZIP archive at its top level, each stored without further compression, so the web app reads any chunk with one range request. The chunks come first and the metadata documents last, the consolidated metadata the very last entry: the web app reads the end of an archive when it opens it, and with it everything that describes the store.

The store is written beside out_path under a temporary name and moved into place only once it is complete: a failed or interrupted call leaves nothing at out_path, and with force = TRUE the store it replaces stays untouched until then.

Note on column encodings: the columns of obs and var are stored with anndata's own rules. A factor is a categorical with its levels in their order (ordered for an ordered factor). A character column is a categorical, its levels in code-point order, when a value repeats or is missing, and a string array otherwise. Category codes take the integer type pandas gives them: int8 for fewer than 127 categories, int16 for fewer than 32,767, else int32. Every other column keeps its own type, as in Python: a logical column is a boolean array, or a nullable boolean when it holds NA; an integer column is int32 (R's integer), or a nullable integer when it holds NA; a double column is float64 whether or not its values are whole numbers, with NA stored as NaN, and must hold no infinite value. Any other column, a Date for example, is refused with the conversion that makes it storable when it is given in obs; a column of an object given as data is left out instead, with a message naming it, as cellucid-python leaves it out.

Note on names: cell names are the row names obs, gene_expression, and the embeddings agree on; where several are set they must be the same names in the same order, and with none the cells are named "0" to "n - 1". Gene names are the var row names, else the column names of gene_expression. Repeated names are made unique the way anndata's obs_names_make_unique() and var_names_make_unique() do it – the second "MATR3" becomes "MATR3-1" – and what changed is reported as a message. Every obs and var column name and every obsm key becomes a path in the store, so it must be a portable file name on every platform: no slash, backslash, colon, asterisk, question mark, quotation mark, <, >, |, or control character, no trailing space or dot, not a Windows device name, not _index for a column, and distinct from its neighbours by more than letter case ("Sex" and "sex" are one folder on macOS and Windows).

Note on displayed text: dataset_name, dataset_description, citation_text, and source_name are drawn verbatim, so each must be free of control characters (U+0001–U+001F, U+007F–U+009F), of the zero-width characters U+200B, U+2060, and U+FEFF, and of leading or trailing whitespace of any kind. Nothing is trimmed: trimming would rewrite text the caller never asked to change.

Note on required arguments: out_path and format carry no default, and without data neither do embeddings, dataset_id, and dataset_name, because no value of any of them is a reasonable guess about someone else's dataset. Leaving one out is reported as what it is – cellucid_prepare() is missing 1 required argument: – followed by the argument's name and what a valid value of it is.

Writing an object

With data, the object's contents are read as the Cellucid web app reads the same object saved with saveRDS(), so a session saved on either applies to the other:

cells and obs: Seurat meta.data (and active.ident as active_ident); colData (and the size factors); columns of a type AnnData does not store are left out and named;
obsm: every Seurat reduction and every reducedDim as X_<name>, the name's letter case kept (X_umap, X_UMAP); one that does not cover every cell is left out and named;
expression: an assay's log-normalised data (else counts) as X; logcounts (else counts, else the first assay) as X; Seurat v5 split layers joined. A Cellucid-ready AnnData holds one matrix per modality, so counts beside data, the other assays, and scale.data are not written;
var: meta.features (Assay5 meta.data) with the variable features as highly_variable; rowData;
obsp/connectivities: the first symmetric neighbour graph the viewer's default rule would open with, its diagonal removed: one named connectivities, then Seurat's <active assay>_snn, then any other shared-nearest-neighbour graph (<assay>_snn, or SNN, in any case; a colPair the same), in code-point order; a graph passed over (not symmetric, or not over every cell) or left out (RNA_nn) is named;
identity and launch defaults: the object's own record, misc$cellucid or metadata(x)$cellucid, shaped as uns['cellucid']: list(id = , name = , description = , citation = list(text = , url = ), source = list(name = , url = ), defaults = list(embedding2d = "X_umap", latentKey = "X_pca", ...)); a Seurat project.name other than "SeuratProject" names it too.

Matrices kept outside R's memory (BPCells, HDF5Array and other DelayedArray backends) stay there and are read a block at a time (about $2^{24}$ positions: whole cells, or whole genes when the backend keeps genes together), never whole and never as one dgCMatrix: one pass counts each gene's values (and refuses a value that cannot be stored), then each value is placed at its gene's position, its cell's row beside it, and each gene sorted by row (on every core; its prefix counts read from its sorted rows). The genes one pass places hold at most a quarter of the machine's memory (8 bytes a stored value, 12 for float64 counts), and the pass reads only their rows; a larger matrix takes one pass per quarter, twice when it is written, X's values being written before its cells. Nothing else of the matrix is held, and there is no limit on its stored values (more than $2^{31} - 1$ are written with 64-bit pointers). The bytes are those the same values in a dgCMatrix give. Seurat objects made before 3.0, and SingleCellExperiment objects before 1.8, are refused with the line that updates them.

Several assays or experiments are written as a Cellucid-ready MuData: MuData's own Zarr layout (mudata.read_zarr() reads it), the cells' obs, obsm, obsp and uns['cellucid'] at the root, each modality an AnnData under mod/<name> (gene-major X and var), mod-order listing them, and obsmap/varmap placing each modality's cells and features. R writes MuData as Zarr only; Python's cellucid.prepare(data = mdata) writes the same store, and .h5mu.

Which format?

Where the store goes decides format:

a GitHub repository or its GitHub Pages site: "zarr". GitHub refuses files over 100 MB, and a folder of small files publishes any size without Git LFS; GitHub answers over HTTP/2, 100 requests at a time, so a folder's one request per chunk costs nothing there.
everything else – a single-file host (Zenodo, Figshare, Hugging Face, S3, a lab web server), your own disk (the web app's file picker), cellucid serve, Jupyter: "zarr.zip". One file is read in ranges that the web app merges into few requests; most single-file hosts answer over HTTP/1.1, six requests at a time, where a folder's one request per chunk makes many-gene and whole-matrix reads several times slower.

Measured on a 561,947-cell atlas, every read the web app makes, in total: over GitHub's HTTP/2 (80 ms round trip, 50 MB/s) 8.0 s from the folder and 8.2 s from the .zarr.zip; over HTTP/1.1 (80 ms, 50 MB/s) 46.8 s and 9.0 s; from a local file in a browser 2.6 and 2.7 s. In a browser, a local Zarr folder reads as fast as a single file. Python's cellucid.prepare(format = "h5ad") writes an .h5ad, which goes wherever a .zarr.zip goes; this package writes Zarr only, so it needs no HDF5 library.

Outlier quantiles

With latent_key, uns['cellucid']['outlierQuantiles'] holds latentKey, minPoints (centroid_min_points, or 10) and fields: one float32 array of one value per cell for every obs column the viewer shows as categorical – a factor, a character column, a logical column with or without NA – in obs column order. Its categories are its values; NA belongs to none. For each category of at least minPoints cells, over the latent values as stored (float32): the centroid is the column sums, accumulated cell by cell in cell order, divided by the cell count; a cell's distance is the square root of its squared differences from the centroid, accumulated column by column; its quantile is the number of the category's cells at that distance or closer divided by the cell count, both as float32. Every other cell is NaN. A column with no category of minPoints cells gets no entry, and a numeric column never does. The viewer's outlier filter hides the cells above a threshold of these quantiles. Python's cellucid.prepare() stores the same bytes, and the web app computes the same bytes for a file that stores none.

Examples

umap <- matrix(c(0, 0, 1, 0, 0, 1, 1, 1), ncol = 2, byrow = TRUE)
pca <- matrix(c(0, 3, 1, 2, 5, 1, 4, 0, 2, 1, 1, 1), ncol = 3, byrow = TRUE)
obs <- data.frame(
  cluster = factor(c("A", "A", "B", "B")),
  score = c(0.1, 0.2, 0.3, 0.4),
  row.names = paste0("cell", 1:4)
)
expr <- matrix(
  c(0, 1, 2, 3, 4, 5, 0, 0),
  nrow = 4,
  dimnames = list(rownames(obs), c("ENSG01", "ENSG02"))
)
var <- data.frame(symbol = c("CD3E", "MS4A1"), row.names = colnames(expr))

store <- cellucid_prepare(
  out_path = file.path(tempdir(), "example.zarr"),
  format = "zarr",
  embeddings = list(X_umap = umap, X_pca = pca),
  dataset_id = "cellucid-example",
  dataset_name = "Cellucid example",
  obs = obs,
  var = var,
  gene_expression = expr,
  embedding_defaults = c("2d" = "X_umap"),
  latent_key = "X_pca",
  var_gene_id_column = "symbol",
  force = TRUE
)
list.files(store, all.files = TRUE, no.. = TRUE)
#>  [1] ".zattrs"    ".zgroup"    ".zmetadata" "X"          "layers"    
#>  [6] "obs"        "obsm"       "obsp"       "uns"        "var"       
#> [11] "varm"       "varp"      

# A Seurat object, read as the web app reads it
if (requireNamespace("SeuratObject", quietly = TRUE)) {
  data("pbmc_small", package = "SeuratObject")
  cellucid_prepare(
    file.path(tempdir(), "pbmc_small.zarr"), "zarr",
    data = SeuratObject::UpdateSeuratObject(pbmc_small),
    dataset_id = "pbmc-small", dataset_name = "PBMC (small)",
    force = TRUE
  )
}
#> Validating object structure
#> Updating object slots
#> Ensuring keys are in the proper structure
#> Updating matrix keys for DimReduc ‘pca’
#> Updating matrix keys for DimReduc ‘tsne’
#> Ensuring keys are in the proper structure
#> Ensuring feature names don't have underscores or pipes
#> Updating slots in RNA
#> Updating slots in RNA_snn
#> Setting default assay of RNA_snn to RNA
#> Updating slots in pca
#> Updating slots in tsne
#> Setting tsne DimReduc to global
#> Setting assay used for NormalizeData.RNA to RNA
#> Setting assay used for ScaleData.RNA to RNA
#> Setting assay used for RunPCA.RNA to RNA
#> Setting assay used for BuildSNN.RNA.pca to RNA
#> No assay information could be found for FindClusters
#> Warning: Adding a command log without an assay associated with it
#> Setting assay used for RunTSNE.pca to RNA
#> Setting assay used for JackStraw.RNA.pca to RNA
#> Setting assay used for ScoreJackStraw.pca to RNA
#> Setting assay used for ProjectDim.RNA.pca to RNA
#> Setting assay used for FindVariableFeatures.RNA to RNA
#> Validating object structure for Assay ‘RNA’
#> Validating object structure for Graph ‘RNA_snn’
#> Validating object structure for DimReduc ‘pca’
#> Validating object structure for DimReduc ‘tsne’
#> Object representation is consistent with the most current Seurat version

# The same store as one file
cellucid_prepare(
  out_path = file.path(tempdir(), "example.zarr.zip"),
  format = "zarr.zip",
  embeddings = list(X_umap = umap),
  dataset_id = "cellucid-example",
  dataset_name = "Cellucid example",
  force = TRUE
)