Streaming large files¶
A whole-slide dataset can be several gigabytes on disk and rather more once
AnnData has materialised X, raw.X, and every layer. Most spatial
statistics do not need all of that at once:
[sp.tl.spatialDistance][scimappro.tl.spatialDistance] needs four obs
columns, [sp.tl.phenotype][scimappro.tl.phenotype] needs the marker matrix a
group at a time.
streamData=True runs those functions against the file on disk through
CAP-AnnData, reading only what
a function touches and writing back only what it changed.
The demo data is small, so this tutorial demonstrates the contract rather than the speed-up — but the contract is what trips people up.
from pathlib import Path
import matplotlib.pyplot as plt
import pandas as pd
import scimappro as sp
# The demo data is not shipped with the package (it is 247 MB). Point DATA at
# wherever you unpacked it; this looks in the repository root first, then in the
# path the docs are built from.
DATA = next(
path for path in (Path("example_data"), Path("../../../example_data"))
if path.exists()
)
DATA
PosixPath('../../../example_data')
import anndata as ad
adata = ad.read_h5ad(DATA / "adata_scimap.h5ad")
adata = sp.pp.rescale(adata, gate=str(DATA / "manual_gates.csv"), verbose=False)
adata = sp.tl.phenotype(adata, phenotype=str(DATA / "phenotype_workflow.csv"), verbose=False)
adata.obs["phenotype"].value_counts()
phenotype ECAD+ 7112 Other myeloid cells 2419 Dendritic cells 863 SMA+ 509 Treg 216 Unknown 46 NK cells 35 Immune 1 Name: count, dtype: int64
Rule 1: data must be a path¶
Streaming an in-memory object makes no sense — it is already in memory. Passing one warns and runs normally.
import warnings
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always")
sp.tl.spatialDistance(adata.copy(), phenotype="phenotype",
streamData=True, verbose=False)
print(caught[0].message)
streamData=True needs a path to an .h5ad file; running in memory on the AnnData that was passed.
Rule 2: the file is modified in place, and None comes back¶
This is the one place scimappro mutates your input. Copy the file first if you want to keep the original.
outputDir = Path("tutorial_output")
outputDir.mkdir(exist_ok=True)
streamPath = outputDir / "streamed.h5ad"
adata.write(streamPath)
returned = sp.tl.spatialDistance(str(streamPath), phenotype="phenotype",
streamData=True, verbose=False)
print("returned:", returned)
[Streaming Mode] Writing 'spatial_distance' to the 'uns' section... returned: None
onDisk = ad.read_h5ad(streamPath)
onDisk.uns["spatial_distance"].head()
| Dendritic cells | ECAD+ | Immune | NK cells | Other myeloid cells | SMA+ | Treg | Unknown | |
|---|---|---|---|---|---|---|---|---|
| obs | ||||||||
| quant_1 | 554.265017 | 508.856310 | 0.000000 | 2053.512851 | 532.311179 | 505.979757 | 575.266490 | 2086.992073 |
| quant_2 | 114.340629 | 10.602318 | 776.673240 | 1287.763826 | 25.455242 | 27.967658 | 225.297772 | 1311.182079 |
| quant_3 | 122.911429 | 10.602318 | 771.606273 | 1292.597225 | 24.200678 | 34.714545 | 226.977790 | 1315.778661 |
| quant_4 | 11.936186 | 13.731018 | 890.139838 | 1178.924606 | 39.819156 | 36.240050 | 168.802786 | 1204.626638 |
| quant_5 | 135.256707 | 15.790519 | 755.938553 | 1308.403525 | 40.006555 | 58.347658 | 215.734114 | 1331.056574 |
Rule 3: the answer is the same¶
Streaming is an execution strategy, not a different algorithm. Both paths give the same numbers.
import numpy as np
inMemory = sp.tl.spatialDistance(adata.copy(), phenotype="phenotype", verbose=False)
a = inMemory.uns["spatial_distance"]
b = onDisk.uns["spatial_distance"].reindex(columns=a.columns)
np.allclose(a.to_numpy(dtype=float), b.to_numpy(dtype=float))
True
Stochastic methods — UMAP, LDA, NMF, Leiden — agree to floating-point noise
rather than exactly, which the stream-vs-memory tests in tests/ cover.
A whole pipeline on disk¶
Each step edits the same file, so nothing round-trips through memory and there are no intermediate files to manage.
pipelinePath = outputDir / "pipeline.h5ad"
ad.read_h5ad(DATA / "adata_scimap.h5ad").write(pipelinePath)
sp.pp.rescale(str(pipelinePath), gate=str(DATA / "manual_gates.csv"),
streamData=True, verbose=False)
sp.tl.phenotype(str(pipelinePath), phenotype=str(DATA / "phenotype_workflow.csv"),
streamData=True, verbose=False)
sp.tl.spatialDistance(str(pipelinePath), phenotype="phenotype",
streamData=True, verbose=False)
sp.tl.spatialCooccurrence(str(pipelinePath), phenotype="phenotype",
permutation=100, streamData=True, verbose=False)
final = ad.read_h5ad(pipelinePath)
print(final.obs["phenotype"].value_counts().to_dict())
print(sorted(final.uns))
[Streaming Mode] Writing 'spatial_distance' to the 'uns' section...
/Users/aj/Partners HealthCare Dropbox/Ajit Nirmal/nirmal lab/softwares/dev/scimappro/.venv/lib/python3.12/site-packages/joblib/externals/loky/process_executor.py:787: UserWarning: A worker stopped while some jobs were given to the executor. This can be caused by a too short worker timeout or by a memory leak. warnings.warn(
{'ECAD+': 7112, 'Other myeloid cells': 2419, 'Dendritic cells': 863, 'SMA+': 509, 'Treg': 216, 'Unknown': 46, 'NK cells': 35, 'Immune': 1}
['all_markers', 'gates', 'spatialCooccurrence', 'spatial_distance']
outputDir in streaming mode¶
Streaming edits the input file, and most functions ignore outputDir
in this mode. The exceptions are [pp.log1p][scimappro.pp.log1p],
[pp.combat][scimappro.pp.combat], [tl.umap][scimappro.tl.umap], and
[tl.spatialSimilarityLookup][scimappro.tl.spatialSimilarityLookup], which
copy the file there after modifying it in place — so the original is
still changed either way.
To keep an untouched original, copy the file yourself first.
keepPath = outputDir / "keep_original.h5ad"
workPath = outputDir / "work_copy.h5ad"
ad.read_h5ad(DATA / "adata_scimap.h5ad").write(keepPath)
import shutil
shutil.copyfile(keepPath, workPath)
sp.tl.cluster(str(workPath), method="kmeans", k=5, label="km",
streamData=True, verbose=False)
print("original untouched:", "km" not in ad.read_h5ad(keepPath).obs)
print("copy updated: ", "km" in ad.read_h5ad(workPath).obs)
original untouched: True copy updated: True
Controlling memory¶
maxWorkers bounds the parallelism, and each worker holds its own slice — so
lowering it lowers peak memory.
sp.tl.spatialCooccurrence(str(pipelinePath), phenotype="phenotype",
permutation=100, maxWorkers=2, streamData=True,
verbose=False)
print("done")
done
Which functions stream¶
pp: rescale, log1p, combat.
tl: phenotype, cluster, umap, foldChange, spatialDistance,
spatialCooccurrence, spatialProximityScore, spatialAggregate,
spatialSimilarityLookup, neighCount, neighExp, neighLDA, neighNMF.
pl: barplot, which reads obs only.
Anything not listed reads little enough that streaming would not help, or needs the whole matrix anyway.
Plotting streamed results¶
Plotting functions do not stream, but backed="r" gets you a light handle for
anything that only reads uns or obs.
backed = ad.read_h5ad(pipelinePath, backed="r")
sp.pl.spatialCooccurrence(backed)
Caveats¶
!!! warning "The file is open for writing"
While a streaming call runs, the .h5ad is held open in r+ mode. Do not
read it from another process at the same time, and do not interrupt a
streaming call mid-write — the file may be left partially updated.
!!! warning "SpatialData is not streamable"
streamData=True applies to .h5ad paths only. With a SpatialData or a
.zarr store it warns and continues in memory, since the tables are already
loaded.
Next¶
- The streaming guide — what each function reads and writes.
- Command line — the same pipeline as a shell script.