Other helpers¶

The small functions that do not belong to any one analysis: relabelling, tidying, batch correction, and a couple of building blocks.

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

tl.rename — merge or relabel categories¶

The mapping is keyed by the new label, with one old label or a list of them as the value. That is scimap's convention, and it reads naturally when several categories collapse into one.

A new column is written, so the original labels stay available.

In [3]:

adata = sp.tl.rename(
    adata,
    rename={"Tumour": "ECAD+",
            "Immune": ["Treg", "NK cells", "Dendritic cells", "Other myeloid cells"],
            "Stroma": "SMA+"},
    fromColumn="phenotype",
    toColumn="broad",
    verbose=False,
)
adata.obs["broad"].value_counts()

adata = sp.tl.rename( adata, rename={"Tumour": "ECAD+", "Immune": ["Treg", "NK cells", "Dendritic cells", "Other myeloid cells"], "Stroma": "SMA+"}, fromColumn="phenotype", toColumn="broad", verbose=False, ) adata.obs["broad"].value_counts()

Out[3]:

broad
Tumour     7112
Immune     3534
Stroma      509
Unknown      46
Name: count, dtype: int64

In [4]:

sp.pl.spatialScatterPlot(adata, colorBy=["phenotype", "broad"], s=2, ncols=2,
                         figsize=(4.5, 4.5), fontSize=6)

sp.pl.spatialScatterPlot(adata, colorBy=["phenotype", "broad"], s=2, ncols=2, figsize=(4.5, 4.5), fontSize=6)

tl.classify — label cells by a marker rule¶

A cell passes when it is at or above threshold for every marker in pos and below it for every marker in neg.

In [5]:

adata = sp.tl.classify(
    adata,
    pos=["CD45"],
    neg=["ECAD"],
    classifyLabel="CD45+ECAD-",
    failedLabel="Other",
    label="cd45_only",
    verbose=False,
)
adata.obs["cd45_only"].value_counts()

adata = sp.tl.classify( adata, pos=["CD45"], neg=["ECAD"], classifyLabel="CD45+ECAD-", failedLabel="Other", label="cd45_only", verbose=False, ) adata.obs["cd45_only"].value_counts()

Out[5]:

cd45_only
Other         11135
CD45+ECAD-       66
Name: count, dtype: int64

subclassifyPhenotype restricts the rule to cells that already carry particular labels, and collapseFailed=True gives everything else its existing phenotype back — so you split one population without disturbing the rest.

In [6]:

adata = sp.tl.classify(
    adata,
    pos=["CD16"],
    phenotype="phenotype",
    subclassifyPhenotype=["Other myeloid cells"],
    classifyLabel="CD16+",
    collapseFailed=True,
    showPhenotypeLabel=True,
    label="phenotype_refined",
    verbose=False,
)
adata.obs["phenotype_refined"].value_counts()

adata = sp.tl.classify( adata, pos=["CD16"], phenotype="phenotype", subclassifyPhenotype=["Other myeloid cells"], classifyLabel="CD16+", collapseFailed=True, showPhenotypeLabel=True, label="phenotype_refined", verbose=False, ) adata.obs["phenotype_refined"].value_counts()

Out[6]:

phenotype_refined
ECAD+                        7112
Other myeloid cells-CD16+    2419
Dendritic cells               863
SMA+                          509
Treg                          216
Unknown                        46
NK cells                       35
Immune                          1
Name: count, dtype: int64

pp.dropFeatures — remove markers, cells, columns, groups¶

Everything in one pass: groups first, then cells, then markers, then obs columns.

In [7]:

trimmed = sp.pp.dropFeatures(
    adata.copy(),
    dropMarkers=["NCAM"],
    dropGroups=["Unknown"],
    groupsColumn="phenotype",
    dropMetaColumns=["Solidity", "Extent"],
    verbose=False,
)
print(adata.shape, "->", trimmed.shape)
print("NCAM gone from raw:", "NCAM" not in trimmed.raw.var_names)

trimmed = sp.pp.dropFeatures( adata.copy(), dropMarkers=["NCAM"], dropGroups=["Unknown"], groupsColumn="phenotype", dropMetaColumns=["Solidity", "Extent"], verbose=False, ) print(adata.shape, "->", trimmed.shape) print("NCAM gone from raw:", "NCAM" not in trimmed.raw.var_names)

(11201, 9) -> (11155, 8)
NCAM gone from raw: True

subsetRaw=True (the default) prunes .raw alongside .X. Forget it and later calls with layer='raw' will disagree with .X.

!!! warning "SpatialData tables need their annotation columns" For SpatialData input, keep the table's region_key and instance_key obs columns — dropping them makes the table fail validation on write-back.

pp.log1p — log-transform into a layer¶

The original matrix is untouched; the transform lands in a layer.

In [8]:

logged = sp.pp.log1p(adata.copy(), targetLayer="log_manual", verbose=False)
list(logged.layers)

logged = sp.pp.log1p(adata.copy(), targetLayer="log_manual", verbose=False) list(logged.layers)

Out[8]:

['log', None, 'log_manual']

pp.combat — batch correction¶

ComBat removes additive and multiplicative batch effects while preserving biological variation. The demo data is one image, so we split it in two to have something to correct between.

!!! note These are synthetic groups, not biological ones — the point is to show the call, not a real batch effect.

In [9]:

batched = adata.copy()
median_x = batched.obs["X_centroid"].median()
batched.obs["batch"] = ["a" if x < median_x else "b" for x in batched.obs["X_centroid"]]
batched.obs["batch"].value_counts()

batched = adata.copy() median_x = batched.obs["X_centroid"].median() batched.obs["batch"] = ["a" if x < median_x else "b" for x in batched.obs["X_centroid"]] batched.obs["batch"].value_counts()

Out[9]:

batch
b    5601
a    5600
Name: count, dtype: int64

In [10]:

corrected = sp.pp.combat(batched, batch="batch", layer="raw", log=True,
                         label="combat", verbose=False)
list(corrected.layers)

corrected = sp.pp.combat(batched, batch="batch", layer="raw", log=True, label="combat", verbose=False) list(corrected.layers)

Found 2 batches.
Adjusting for 0 covariate(s) or covariate level(s).
Standardizing Data across genes.
Fitting L/S model and finding priors.
Finding parametric adjustments.
Adjusting the Data

Out[10]:

['log', None, 'combat']

In [11]:

sp.pl.heatmap(corrected, groupBy="batch", layer="combat", standardScale=None,
              clusterRows=False, clusterColumns=False, figsize=(5, 2))

sp.pl.heatmap(corrected, groupBy="batch", layer="combat", standardScale=None, clusterRows=False, clusterColumns=False, figsize=(5, 2))

replaceOriginal=True writes the corrected matrix into .X instead of a layer. A single-batch object raises an Exception — there is nothing to correct between.

pp.nGraph — a k-NN graph you can use yourself¶

Returns a bare igraph.Graph rather than modifying the object. [sp.tl.cluster][scimappro.tl.cluster] builds its own graph internally for Leiden; this is for when you want igraph's own algorithms.

In [12]:

graph = sp.pp.nGraph(adata, layer="raw", standardScale=True, runPCA=True,
                     kNeighbors=15, nPcs=8)
print(graph.vcount(), "vertices,", graph.ecount(), "edges")

graph = sp.pp.nGraph(adata, layer="raw", standardScale=True, runPCA=True, kNeighbors=15, nPcs=8) print(graph.vcount(), "vertices,", graph.ecount(), "edges")

11201 vertices, 168015 edges

In [13]:

communities = graph.community_multilevel()
adata.obs["igraph_community"] = [str(c) for c in communities.membership]
adata.obs["igraph_community"].value_counts().head()

communities = graph.community_multilevel() adata.obs["igraph_community"] = [str(c) for c in communities.membership] adata.obs["igraph_community"].value_counts().head()

Out[13]:

igraph_community
10    1123
2      977
9      852
14     809
5      680
Name: count, dtype: int64

pp.mergeAdataObs — stack metadata across objects¶

Reads only obs, so it stays cheap across many large files.

In [14]:

outputDir = Path("tutorial_output")
outputDir.mkdir(exist_ok=True)
adata.write(outputDir / "helper_a.h5ad")
adata.write(outputDir / "helper_b.h5ad")

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

outputDir = Path("tutorial_output") outputDir.mkdir(exist_ok=True) adata.write(outputDir / "helper_a.h5ad") adata.write(outputDir / "helper_b.h5ad") obs = sp.pp.mergeAdataObs( [str(outputDir / "helper_a.h5ad"), str(outputDir / "helper_b.h5ad")], verbose=False, ) obs.shape

Out[14]:

(22402, 16)

pp.scimapToCsv — a flat table¶

Covered in Export data.

Not ported from scimap¶

scimap	Status
sm.hl.animate	Not ported.
sm.pl.gate_finder, sm.pl.napariGater	Not ported. Use [sp.pp.rescale][scimappro.pp.rescale]'s automatic GMM gates, or a manual_gates.csv — see Prior-knowledge phenotyping.
sm.tl.spatial_cluster	Folded into tl.cluster(mode="spatial", layer=<uns key>) — see Latent motifs.

The full mapping is in Migrating from scimap.

Next¶

• Export data — getting results out.
• Command line — running these steps from a shell script.