Export data¶

Getting results out of scimappro, for collaborators, for a paper, or for a tool that does not read .h5ad.

In [1]:

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

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

Out[1]:

PosixPath('../../../example_data')

In [2]:

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()

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()

Out[2]:

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

In [3]:

adata = sp.tl.spatialDistance(adata, phenotype="phenotype", verbose=False)
adata = sp.tl.spatialCooccurrence(adata, phenotype="phenotype", permutation=200,
                                  verbose=False)
outputDir = Path("tutorial_output")
outputDir.mkdir(exist_ok=True)
sorted(adata.uns)

adata = sp.tl.spatialDistance(adata, phenotype="phenotype", verbose=False) adata = sp.tl.spatialCooccurrence(adata, phenotype="phenotype", permutation=200, verbose=False) outputDir = Path("tutorial_output") outputDir.mkdir(exist_ok=True) sorted(adata.uns)

Out[3]:

['all_markers', 'gates', 'spatialCooccurrence', 'spatial_distance']

The whole object¶

.h5ad is the format to keep. It carries everything — expression, metadata, every result in uns, obsm, and layers — and every scimappro function accepts the path directly, so you never have to load it to continue.

In [4]:

adata.write(outputDir / "analysis.h5ad")
(outputDir / "analysis.h5ad").stat().st_size / 1e6

adata.write(outputDir / "analysis.h5ad") (outputDir / "analysis.h5ad").stat().st_size / 1e6

Out[4]:

5.213348

A flat CSV¶

[sp.pp.scimapToCsv][scimappro.pp.scimapToCsv] gives one row per cell: a cell id, one column per marker, then every obs column. That is the shape most non-Python tools expect.

In [5]:

flat = sp.pp.scimapToCsv(adata, layer="raw", outputDir=str(outputDir),
                         fileName="cells.csv", verbose=False)
flat.shape

flat = sp.pp.scimapToCsv(adata, layer="raw", outputDir=str(outputDir), fileName="cells.csv", verbose=False) flat.shape

Out[5]:

(11201, 22)

In [6]:

flat.head()

flat.head()

Out[6]:

	CellID	ELANE	CD57	CD45	CD11B	SMA	CD16	ECAD	FOXP3	NCAM	...	Area	MajorAxisLength	MinorAxisLength	Eccentricity	Solidity	Extent	Orientation	CellID	imageid	phenotype
0	quant_1	1179.529915	161.410256	708.111111	378.692308	246.256410	332.769231	707.367521	205.401709	853.205128	...	117	12.402944	12.006487	0.250814	0.959016	0.812500	-1.146733	1	quant	Immune
1	quant_2	1124.728261	177.217391	790.271739	910.434783	423.445652	649.989130	1027.836957	387.500000	902.489130	...	92	11.874070	9.982065	0.541562	0.948454	0.696970	-0.435290	2	quant	ECAD+
2	quant_3	1249.500000	156.637931	791.189655	1034.827586	477.362069	599.482759	1334.327586	287.810345	936.741379	...	58	10.113305	7.629922	0.656364	0.878788	0.585859	1.221658	3	quant	ECAD+
3	quant_4	1159.676580	187.431227	853.814126	816.914498	2221.639405	375.884758	1030.483271	310.550186	891.988848	...	269	25.433196	15.183300	0.802251	0.835404	0.531621	-0.705293	4	quant	SMA+
4	quant_5	1568.526961	162.904412	318.497549	827.286765	362.482843	366.365196	1325.740196	269.232843	896.955882	...	408	26.604670	19.759781	0.669604	0.937931	0.739130	-0.711002	5	quant	ECAD+

5 rows × 22 columns

layer=None exports the rescaled values in .X instead of the raw intensities:

In [7]:

rescaled = sp.pp.scimapToCsv(adata, layer=None, verbose=False)
rescaled[["CellID", "CD45", "ECAD", "phenotype"]].head()

rescaled = sp.pp.scimapToCsv(adata, layer=None, verbose=False) rescaled[["CellID", "CD45", "ECAD", "phenotype"]].head()

Out[7]:

	CellID	CellID	CD45	ECAD	phenotype
0	quant_1	1	0.614936	0.604100	Immune
1	quant_2	2	0.637269	0.663324	ECAD+
2	quant_3	3	0.637505	0.704700	ECAD+
3	quant_4	4	0.653004	0.663731	SMA+
4	quant_5	5	0.292861	0.703676	ECAD+

Individual results¶

Everything in uns is a DataFrame, so exporting one is a to_csv call.

In [8]:

adata.uns["spatial_distance"].to_csv(outputDir / "spatial_distance.csv")
adata.uns["spatialCooccurrence"].to_csv(outputDir / "cooccurrence.csv", index=False)
sorted(p.name for p in outputDir.iterdir())

adata.uns["spatial_distance"].to_csv(outputDir / "spatial_distance.csv") adata.uns["spatialCooccurrence"].to_csv(outputDir / "cooccurrence.csv", index=False) sorted(p.name for p in outputDir.iterdir())

Out[8]:

['adata_scimap.zarr',
 'analysis.h5ad',
 'cells.csv',
 'cooccurrence.csv',
 'keep_original.h5ad',
 'pipeline.h5ad',
 'spatial_distance.csv',
 'streamed.h5ad',
 'work_copy.h5ad']

Plotting functions hand back the exact table behind a figure with returnData=True, which is usually what you want for a supplementary table — the numbers as plotted, not as stored.

In [9]:

plotted = sp.pl.spatialDistanceHeatmap(adata, phenotype="phenotype",
                                       returnData=True, show=False)
plotted.round(1)

plotted = sp.pl.spatialDistanceHeatmap(adata, phenotype="phenotype", returnData=True, show=False) plotted.round(1)

Out[9]:

	Dendritic cells	ECAD+	Immune	NK cells	Other myeloid cells	SMA+	Treg	Unknown
phenotype
Dendritic cells	17.9	23.7	1594.1	805.2	24.8	74.4	55.4	1246.9
ECAD+	146.8	15.2	1889.4	562.3	37.7	65.7	94.2	1214.2
Immune	554.3	508.9	0.0	2053.5	532.3	506.0	575.3	2087.0
NK cells	127.8	33.0	2336.4	51.0	85.0	93.9	152.8	311.4
Other myeloid cells	88.4	18.4	1766.7	642.5	17.6	58.6	63.8	1194.6
SMA+	141.1	18.4	1969.2	544.9	30.1	36.1	90.0	1219.1
Treg	101.1	19.0	1980.2	534.8	20.2	54.9	62.8	1193.7
Unknown	127.9	30.6	2398.3	63.6	89.5	102.1	195.6	19.6

Just the metadata¶

[sp.pp.mergeAdataObs][scimappro.pp.mergeAdataObs] stacks the obs tables of several objects without touching their expression matrices, which is cheap even across many large files.

In [10]:

adata.write(outputDir / "sample_a.h5ad")
adata.write(outputDir / "sample_b.h5ad")

obs = sp.pp.mergeAdataObs(
    [str(outputDir / "sample_a.h5ad"), str(outputDir / "sample_b.h5ad")],
    verbose=False,
)
obs.shape

adata.write(outputDir / "sample_a.h5ad") adata.write(outputDir / "sample_b.h5ad") obs = sp.pp.mergeAdataObs( [str(outputDir / "sample_a.h5ad"), str(outputDir / "sample_b.h5ad")], verbose=False, ) obs.shape

Out[10]:

(22402, 12)

In [11]:

obs["phenotype"].value_counts()

obs["phenotype"].value_counts()

Out[11]:

phenotype
ECAD+                  14224
Other myeloid cells     4838
Dendritic cells         1726
SMA+                    1018
Treg                     432
Unknown                   92
NK cells                  70
Immune                     2
Name: count, dtype: int64

Figures¶

Every plotting function takes outputDir plus fileName; the extension decides the format. Use show=False in scripts so nothing tries to open a window.

In [12]:

figures = outputDir / "figures"
sp.pl.heatmap(adata, groupBy="phenotype", outputDir=str(figures),
              fileName="heatmap.pdf", show=False)
sp.pl.spatialScatterPlot(adata, colorBy="phenotype", s=2, outputDir=str(figures),
                         fileName="spatial.png", show=False)
sorted(p.name for p in figures.iterdir())

figures = outputDir / "figures" sp.pl.heatmap(adata, groupBy="phenotype", outputDir=str(figures), fileName="heatmap.pdf", show=False) sp.pl.spatialScatterPlot(adata, colorBy="phenotype", s=2, outputDir=str(figures), fileName="spatial.png", show=False) sorted(p.name for p in figures.iterdir())

Out[12]:

['heatmap.pdf', 'spatial.png']

transparent=True and dpi are there for figures going into a manuscript.

!!! note "saveDir is deprecated" Some plotting functions still accept saveDir alongside outputDir, for scripts carried over from scimap. outputDir wins when both are given. Use outputDir.

To SpatialData¶

For the scverse ecosystem, [sp.pp.toSpatialData][scimappro.pp.toSpatialData] writes a .zarr store — see the SpatialData tutorial.

In [13]:

sp.pp.toSpatialData(adata, outputDir=str(outputDir / "stores"), verbose=False)
sorted(p.name for p in (outputDir / "stores").iterdir())

sp.pp.toSpatialData(adata, outputDir=str(outputDir / "stores"), verbose=False) sorted(p.name for p in (outputDir / "stores").iterdir())

Out[13]:

['adata_scimap.zarr']

Next¶

• Other helpers — tidying an object before export.
• SpatialData workflow — the scverse route.