Scarf Single Cell

Scarf single-cell analysis

Scarf streams Zarr-backed count matrices in bounded blocks and records every computation as an immutable, content-addressed artifact. This skill encodes how to drive core Scarf for single-cell RNA analysis: which calls to make, in what order, what evidence to check, and what not to do. It does not use or describe scarf.agent.

Read this file first. Then open only the reference modules you need (index below). Paths in this skill are relative to the skill directory. The maintained copy lives in the Scarf repository at https://github.com/NygenAnalytics/scarf/tree/master/skills/scarf-single-cell. Every recipe in the modules was run against a real store. Numbers quoted there come from the 10x 5K PBMC documentation dataset and are illustrations, not thresholds.

When to Use This Skill

• Use when a task involves a Scarf .zarr store, scarf.DataStore, ds.pipeline or a PipelineRun.
• Use when converting H5AD, 10x, MTX or Seurat data into a Scarf store, or opening and mounting a Cytebase dataset.
• Use for single-cell RNA-seq analysis with core Scarf: QC with removal audits, HVG, PCA and neighbour graphs, Leiden and Paris clustering, UMAP, markers and cautious annotation, batch correction, donor-level comparisons, headless plotting, provenance and export.
• Use when counts are too large for memory or live on object storage, so they must be streamed under an explicit memory budget.
• Do not use it for scarf.agent (the automated agent package), which this skill does not cover.

Setup

• Install into the environment you run Python from: pip install "scarf[extra]>=1.0.0rc17", plus scarf[cytebase] for Cytebase. The explicit pre-release floor matters: a bare scarf[extra] resolves to the old 0.32 series, whose API this skill does not describe.
• Set resources per process before importing Scarf. Defaults claim all detected RAM and every CPU: SCARF_MEM_BUDGET=8G SCARF_WORKERS=8 (memory specs need a unit; a bare 8 is rejected).
• Headless: MPLBACKEND=Agg, call plots with show=False, then result.save(path) and result.close().
• Quiet logs for scripts: scarf.configure_output(level="WARNING", progress=False).
• First look at a store: open it read-only and follow "Open and inspect a store" (references/data-access.md) and "Check the matrix first" (references/quality-control.md) for cell and feature counts, runs, artifact kinds and whether counts look raw. Then review every metadata column yourself: design columns (donor, sample, batch, condition), author annotations to hold out (rule 14) and author-derived columns. The maintained copy of this skill in the Scarf repository also ships scripts/inspect_store.py, a read-only command that runs these checks and flags those columns by name (annotation values stay hidden unless you pass --show-annotation-values). This catalogue copy is Markdown only and does not include the script.
• Long steps (a pipeline run on tens of thousands of cells over remote counts can take many minutes): write the step as a script that logs to a file and ends with a DONE or FAILED line, run it in the background, and wait with a bounded loop that also stops on a Traceback or a dead process. If you launch through a wrapper such as uv run, put timeout inside it (uv run timeout N python step.py; the reverse can leave stale running records), and track the process by $! rather than pgrep -f. See references/performance-and-export.md. Never poll without a time limit.

Mental model

• A DataStore is one Zarr directory: a shared cell table ds.cells, one group per assay (ds.RNA, feature table ds.RNA.feats, counts as cell-major counts plus gene-major countsT), artifacts, and pipeline/runs. QC columns are assay-prefixed: RNA_nCounts, RNA_nFeatures, RNA_percentMito, RNA_percentRibo.
• I is the live boolean cell key. Analytical filters never edit it; they return a frozen cell selection artifact.
• Every result is an ArtifactRef (scope, kind, artifact_id, assay). Producers take refs and return refs: select_hvgs -> run_normalization -> run_pca -> build_ann_index -> query_neighbors -> build_connectivity_map -> run_leiden_clustering / run_umap -> run_marker_search. An identical call returns the existing artifact without recomputing.
• ds.pipeline.run(label=...) runs that whole RNA recipe (filtering, cell cycle, HVG, normalization, PCA, graph, UMAP, Leiden at 0.5/0.75/1.0/1.25 with a silhouette pick, Paris, doublet scores, markers) and returns a durable PipelineRun: a mapping of output names to refs plus frozen run.cells / run.features views and run.report().
• Artifact payload rows follow the artifact's cell selection, not ds.cells.N.
• A mount is a local writable store whose counts stay in a remote source (matrixSource). Every pass over counts is a network read; artifacts are written locally.

Golden rules

1. Inspect read-only first. scarf.DataStore(path, zarr_mode="r") writes nothing. A writable open (the default) prepares new stores and applies min_features_per_cell (default 10) to I permanently. Reopen existing stores and mounts with scarf.DataStore(path, min_features_per_cell=-1).
2. Never filter by editing I. Keep the cell_selection ref a filter returns and pass it on, or insert a boolean column and use cell_key=. ds.cells.reset_key("I") undoes an accidental open-time filter.
3. Keep every ref you will reuse in a dict and record them. Never guess artifact IDs or read private Zarr paths. Use ds.inspect_artifact(ref), ds.load_artifact(ref) and ds.list_artifacts(...). Cell selections need scope="datastore" when listing.
4. Check the matrix, then audit QC removals before accepting a filter. Published matrices can be corrected (for example SCTransform) or already filtered, which makes QC bounds meaningless; the read-only checks in references/quality-control.md ("Check the matrix first") report this. Pooled MAD filters routinely remove low-complexity populations (platelets, erythrocytes, neutrophils) and high-RNA ones (plasma, cycling cells). Audit retention per sample and per marker-defined group without author labels, and look at the markers of removed cells (references/quality-control.md).
5. Labels are immutable and runs cannot be resumed. Use a new label for every variant. A rerun reuses every complete artifact, so it is cheap. Choose snapshot_columns (design columns you want inside run.cells and exports) on the first run, because changing them recomputes everything downstream.
6. Read pipeline outputs instead of recomputing them. run["markers"], run["cell_cycle"] and run["doublets"] already exist. A direct run_marker_search(run["clusters"], ...) creates a new artifact and rereads all counts.
7. Align arrays through the cell selection. Use run.cells.fetch(...) (run cells only) or fetch_all(...) (full length, with -1, NaN or "" outside the run). For explicit artifacts use ds.inspect_artifact(ref).input_ref("cell_selection"). Join tables on cell ids, never on row order.
8. Budget count passes on mounts. QC metrics, the HVG summary, normalization per selection, markers, gene plots and AUCell each read counts. For many passes, make a local copy once: python -m scarf.tools.repack_zarr MOUNT.zarr LOCAL.zarr --mem-budget 4G.
9. One DataStore per process; it is not thread-safe. Run parallel analyses in separate processes on separate stores, each with its own SCARF_MEM_BUDGET/SCARF_WORKERS.
10. Payload names differ by kind. Leiden labels, embeddings, doublet and strength scores live under values; Paris (cluster_cut) under labels; PCA under data. get_markers returns string group_id while Leiden labels are integers.
11. get_markers defaults hide genes. min_score=0.25, min_frac_exp=0.2; pass -1 for both to see negatives and full panels.
12. Labels are hypotheses. Name a cluster only with two or more specific positive markers plus low lineage-negative markers and a plausible QC and doublet profile. Keep unresolved for mixed, doublet-like or marker-poor clusters.
13. Cells are not replicates, and Scarf does not check your design. Condition claims need donor-level aggregation and at least two independent donors per group. A donor sampled twice is repeated measures. A donor with samples in two arms must not count in both: run_statistical_testing(sample_by="sample_id") silently does that, while sample_by="donor_id" refuses such a design. run_harmony runs silently on a column confounded with the biology you want to compare, so audit donors x batch x condition first. When inserting design columns, replace missing values explicitly: ds.cells.insert stores None as "" and pd.NA as "<NA>".
14. Hold out annotation columns when they will judge the result. Do not use author labels (cell_type, author_cell_type, cell.type.*, singler, predicted.*, ...) to choose QC, parameters or labels. Use them only in a final, clearly separated comparison.
15. Pass an explicit HVG blacklist. On stores with current gene symbols the default blacklist misses replication-dependent histones (^HIST matches only pre-2020 names), and its ^CCN family removes the CCN1-6 matricellular genes and non-cell-cycle cyclins. It also leaves Ig V/J genes in, which can split plasma cells by light chain instead of biology. Use the species string from references/gene-blacklists.md (params={"hvg": {"blacklist": ...}} in the pipeline). Add the clonotype add-on only on marker evidence. A blacklist= string replaces the default entirely, and changing it creates new HVG, PCA and graph artifacts.

The analysis loop

Work in this order. After each step, write the decision and its evidence to an analysis log.

Quick start

A complete first pass on an RNA store. Every step after the run reuses the run's artifacts.

from pathlib import Path

import matplotlib

matplotlib.use("Agg")  # headless; select the backend before importing scarf
import numpy as np
import pandas as pd
import scarf

scarf.configure_output(level="WARNING", progress=False)
store = "analysis.zarr"
out = Path("analysis_out")
(out / "figures").mkdir(parents=True, exist_ok=True)

ds = scarf.DataStore(store, min_features_per_cell=-1)

# QC: compare candidate filters and audit them before choosing (quality-control.md).
cells = ds.snapshot_cell_selection("I")
mad5 = ds.auto_filter_cells(cell_selection=cells, n_mads=5)
kept = np.asarray(ds.load_artifact(mad5)["values"][:], dtype=bool)
print(f"MAD 5 keeps {kept.sum()} of {ds.cells.N} cells")

# Baseline run with the chosen policy; freeze design columns you will need later.
run = ds.pipeline.run(
    label="baseline_mad5",
    filtering={"method": "mad", "n_mads": 5},
    snapshot_columns=[c for c in ("donor_id", "sample_id") if c in ds.cells.columns],
)
print(run.status, list(run.keys()))
(out / "run_report.md").write_text(run.report(format="markdown"))

# Figures and tables.
umap = ds.plots.embedding(run=run, color_by="clusters", show=False)
umap.save(out / "figures" / "umap_clusters.png", dpi=150)
umap.close()
markers = ds.get_markers(marker=run["markers"])
markers.to_csv(out / "markers_top.csv", index=False)
print(markers.groupby("group_id", sort=False).head(5)
      .groupby("group_id", sort=False)["feature_name"].agg(", ".join))
sizes = pd.Series(run.cells.fetch("clusters")).value_counts().sort_index()
print(sizes.to_dict())