Skills › Scientific Agent Skills › scanpy › _common.py
_common.py
251 lines · skills/scanpy/scripts/_common.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Shared helpers for the scanpy script toolkit.
5Every CLI script in this directory imports from this module so that data
6loading, saving, figure configuration, and logging behave consistently.
7This file is NOT a CLI itself; import it:
9 from _common import load_anndata, save_anndata, configure_scanpy, info
10"""
12import os
13import sys
16def info(msg):
17 """Print a progress message to stderr-friendly stdout with a marker."""
18 print(f"[scanpy] {msg}", flush=True)
21def die(msg, code=1):
22 """Print an error and exit."""
23 print(f"Error: {msg}", file=sys.stderr, flush=True)
24 sys.exit(code)
27def _import_scanpy():
28 try:
29 import scanpy as sc # noqa: F401
30 return sc
31 except ImportError:
32 die('scanpy not installed. Install with: uv pip install "scanpy[leiden]"')
35def configure_scanpy(figdir="figures", dpi=120, verbosity=1, autosave=False,
36 file_format="png"):
37 """Apply consistent scanpy settings and return the scanpy module.
39 Note: autosave is left off by default because the toolkit scripts pass
40 explicit ``save=`` suffixes to plotting calls for predictable filenames.
41 """
42 sc = _import_scanpy()
43 sc.settings.verbosity = verbosity
44 sc.set_figure_params(dpi=dpi, facecolor="white")
45 sc.settings.figdir = figdir
46 sc.settings.file_format_figs = file_format
47 sc.settings.autosave = autosave
48 sc.settings.autoshow = False
49 os.makedirs(figdir, exist_ok=True)
50 return sc
53def load_anndata(path, var_names="gene_symbols"):
54 """Load an AnnData object, dispatching on the file extension / layout.
56 Supported inputs:
57 * ``.h5ad`` -> sc.read_h5ad
58 * ``.h5`` (10x CellRanger) -> sc.read_10x_h5
59 * ``.csv`` / ``.tsv`` / ``.txt`` -> sc.read_csv / read_text
60 * ``.loom`` -> sc.read_loom
61 * ``.mtx`` -> sc.read (matrix market)
62 * a directory -> sc.read_10x_mtx (10x mtx folder)
63 """
64 sc = _import_scanpy()
65 if not os.path.exists(path):
66 die(f"input not found: {path}")
68 path = os.fspath(path)
69 lower = path.lower()
70 if os.path.isdir(path):
71 info(f"Reading 10x mtx directory: {path}")
72 return sc.read_10x_mtx(path, var_names=var_names)
73 if lower.endswith(".h5ad"):
74 return sc.read_h5ad(path)
75 if lower.endswith(".h5"):
76 info("Reading 10x HDF5 (.h5)")
77 return sc.read_10x_h5(path)
78 if lower.endswith(".loom"):
79 return sc.read_loom(path)
80 if lower.endswith(".csv"):
81 return sc.read_csv(path)
82 if lower.endswith((".tsv", ".txt")):
83 return sc.read_text(path)
84 if lower.endswith(".mtx") or lower.endswith(".mtx.gz"):
85 return sc.read_mtx(path)
86 die(f"unrecognized input format: {path}")
89def save_anndata(adata, path):
90 """Write an AnnData object to .h5ad, creating parent dirs as needed."""
91 parent = os.path.dirname(os.path.abspath(path))
92 os.makedirs(parent, exist_ok=True)
93 adata.write_h5ad(path)
94 info(f"Wrote {path} ({adata.n_obs} cells x {adata.n_vars} genes)")
97def add_io_args(parser, default_output=None):
98 """Attach the standard input/output/figdir arguments to an argparse parser."""
99 parser.add_argument("input", help="Input file (.h5ad, .h5, .csv, .loom, or 10x mtx dir)")
100 parser.add_argument("-o", "--output", default=default_output,
101 help="Output .h5ad path" +
102 (f" (default: {default_output})" if default_output else ""))
103 parser.add_argument("--figdir", default="figures",
104 help="Directory for saved figures (default: figures)")
105 return parser
108def _named_keys(mapping):
109 """Named keys of an AnnData mapping, in order.
111 anndata >= 0.13 reports an unnamed `None` key on `.layers` standing for X
112 itself. Joining that into a string raises TypeError, so filter it out.
113 """
114 return [key for key in mapping.keys() if isinstance(key, str)]
117def summarize(adata):
118 """Return a short human-readable summary string of an AnnData object."""
119 lines = [f"{adata.n_obs} cells x {adata.n_vars} genes"]
120 if len(adata.obs.columns):
121 lines.append("obs: " + ", ".join(adata.obs.columns[:20]))
122 obsm = _named_keys(adata.obsm)
123 if obsm:
124 lines.append("obsm: " + ", ".join(obsm))
125 layers = _named_keys(adata.layers)
126 if layers:
127 lines.append("layers: " + ", ".join(layers))
128 return "\n".join(lines)
131def validate_counts(matrix, label="counts", require_nonzero_cells=True):
132 """Validate finite, nonnegative integer counts without densifying sparse input.
134 Values alone cannot establish provenance: the caller must identify the actual
135 unnormalized assay. Floating dtypes are accepted when every value is integral.
136 """
137 import numpy as np
138 from scipy import sparse
140 if matrix is None or len(matrix.shape) != 2 or min(matrix.shape) == 0:
141 die(f"{label}: expected a nonempty cells x genes count matrix")
142 values = matrix.data if sparse.issparse(matrix) else np.asarray(matrix)
143 if (not np.isfinite(values).all() or (values < 0).any()
144 or not np.equal(values, np.floor(values)).all()):
145 die(f"{label}: expected finite, nonnegative integer counts; verify assay provenance")
146 if require_nonzero_cells and (np.asarray(matrix.sum(axis=1)).ravel() <= 0).any():
147 die(f"{label}: contains zero-total cells; filter them before normalization")
150def prepare_counts(adata, layer=None):
151 """Select explicitly identified counts and preserve an independent copy."""
152 import numpy as np
153 from scipy import sparse
155 if layer is not None and layer not in adata.layers:
156 die(f"counts layer '{layer}' is missing")
157 matrix = adata.X if layer is None else adata.layers[layer]
158 validate_counts(matrix, "X" if layer is None else f"layers[{layer!r}]")
159 if layer is None and "log1p" in adata.uns:
160 die("X has log1p provenance; select the original --counts-layer explicitly")
161 if layer is None and "counts" in adata.layers:
162 old = adata.layers["counts"]
163 delta = old - matrix
164 differs = delta.nnz != 0 if sparse.issparse(delta) else np.any(delta != 0)
165 if differs:
166 die("X differs from the preserved counts layer; select --counts-layer counts")
167 adata.X = matrix.copy()
168 adata.layers["counts"] = matrix.copy()
169 adata.uns.pop("log1p", None)
172def normalize_hvg(sc, adata, *, target_sum=1e4, n_top_genes=2000,
173 flavor="seurat", batch_key=None):
174 """Keep all genes/counts, select HVGs on the flavor's required representation."""
175 if target_sum <= 0 or n_top_genes < 1:
176 die("target_sum and n_top_genes must be positive")
177 if batch_key is not None:
178 ensure_categories(adata, batch_key)
179 count_flavor = flavor in {"seurat_v3", "seurat_v3_paper"}
180 if count_flavor:
181 sc.pp.highly_variable_genes(adata, layer="counts", n_top_genes=n_top_genes,
182 flavor=flavor, batch_key=batch_key)
183 sc.pp.normalize_total(adata, target_sum=target_sum)
184 sc.pp.log1p(adata)
185 if not count_flavor:
186 sc.pp.highly_variable_genes(adata, n_top_genes=n_top_genes,
187 flavor=flavor, batch_key=batch_key)
188 if not adata.var["highly_variable"].any():
189 die("HVG selection returned no genes; inspect expression/filtering parameters")
190 adata.raw = adata.copy()
193def ensure_categories(adata, key):
194 """Use observed, nonmissing categories; numeric labels become lexical strings."""
195 if key not in adata.obs or adata.obs[key].isna().any():
196 die(f"obs column '{key}' must exist and contain no missing labels")
197 adata.obs[key] = adata.obs[key].astype(str).astype("category")
200def compute_pca(sc, adata, n_comps=50):
201 """ARPACK requires fewer components than both cells and selected features."""
202 n_features = int(adata.var["highly_variable"].sum()) if "highly_variable" in adata.var else adata.n_vars
203 actual = min(n_comps, adata.n_obs - 1, n_features - 1)
204 if actual < 1:
205 die("PCA requires at least two cells and two selected genes")
206 sc.pp.pca(adata, n_comps=actual, svd_solver="arpack", random_state=0)
207 return actual
210def build_neighbors(sc, adata, *, use_rep="X_pca", n_pcs=40, n_neighbors=15):
211 if use_rep not in adata.obsm:
212 die(f"embedding '{use_rep}' is missing")
213 if adata.n_obs < 3 or n_pcs < 1 or n_neighbors < 2:
214 die("neighbors requires at least three cells, n_pcs >= 1 and n_neighbors >= 2")
215 actual = min(n_pcs, adata.obsm[use_rep].shape[1])
216 sc.pp.neighbors(adata, use_rep=use_rep, n_pcs=actual,
217 n_neighbors=min(n_neighbors, adata.n_obs - 1), random_state=0)
220def clear_graph(adata, clear_pca=False):
221 """Remove standard derived results after changing expression or cell selection."""
222 for key in ("neighbors", "umap", "paga", "rank_genes_groups", "diffmap_evals", "iroot"):
223 adata.uns.pop(key, None)
224 for key in list(adata.uns):
225 if key.startswith("dendrogram_"):
226 adata.uns.pop(key)
227 for key in ("connectivities", "distances"):
228 adata.obsp.pop(key, None)
229 for key in ("X_umap", "X_tsne", "X_diffmap", "X_draw_graph_fa"):
230 adata.obsm.pop(key, None)
231 if clear_pca:
232 adata.uns.pop("pca", None)
233 adata.varm.pop("PCs", None)
234 for key in ("X_pca", "X_pca_harmony"):
235 adata.obsm.pop(key, None)
238def integrate_harmony(adata, key):
239 """Use harmonypy 2's cells x PCs output; Scanpy 1.12.4 transposes it incorrectly."""
240 import numpy as np
241 import harmonypy
242 from importlib.metadata import version
244 if int(version("harmonypy").split(".", 1)[0]) != 2:
245 die("this Harmony path targets harmonypy 2.x; install harmonypy==2.0.2")
246 pcs = np.asarray(adata.obsm["X_pca"], dtype=np.float64)
247 result = harmonypy.run_harmony(pcs, adata.obs, key, random_state=0, ncores=1)
248 corrected = np.asarray(result.Z_corr)
249 if corrected.shape != pcs.shape or not np.isfinite(corrected).all():
250 die("Harmony returned invalid cells x PCs coordinates")
251 adata.obsm["X_pca_harmony"] = corrected