01_复用前准备/validate_spatial_input.py
python · 9675 bytes
#!/usr/bin/env python3
"""Read-only structural checks for a cell-level CellBin AnnData object."""
from __future__ import annotations
import argparse
import json
import sys
from pathlib import Path
from typing import Any
def _split_names(value: str) -> list[str]:
return [item.strip() for item in value.split(",") if item.strip()]
def _serializable(value: Any) -> Any:
"""Convert common numpy/pandas scalars without importing those packages here."""
try:
return value.item()
except AttributeError:
return value
def _print_human(result: dict[str, Any]) -> None:
summary = result.get("summary", {})
print(f"H5AD: {result['path']}")
if summary:
print(
"shape: "
f"n_obs={summary.get('n_obs', '?')}, "
f"n_vars={summary.get('n_vars', '?')}"
)
print(f"obsm keys: {', '.join(summary.get('obsm_keys', [])) or '(none)'}")
print(f"layers: {', '.join(summary.get('layers', [])) or '(none)'}")
print(f"raw present: {summary.get('raw_present', False)}")
for message in result["warnings"]:
print(f"[WARN] {message}")
for message in result["errors"]:
print(f"[ERROR] {message}")
if result["errors"]:
print("Result: FAIL")
else:
print("Result: PASS (structure only; biological validity is not assessed)")
def main() -> int:
parser = argparse.ArgumentParser(
description="Validate the structural contract used by the lupus skin spatial notebooks."
)
parser.add_argument("--h5ad", type=Path, required=True, help="Path to an AnnData .h5ad file")
parser.add_argument(
"--spatial-key",
default="spatial",
help="Key in adata.obsm containing n_obs x 2 spatial coordinates (default: spatial)",
)
parser.add_argument(
"--required-obs",
default="sample",
help="Comma-separated obs columns that must exist (default: sample)",
)
parser.add_argument(
"--optional-obs",
default="donor,group,celltype,subtype,Region,TLS_region,Epi_Depth,Der_Depth",
help="Comma-separated obs columns to report when present",
)
parser.add_argument(
"--check-genes",
default="",
help="Comma-separated gene names that should exist in adata.var_names",
)
parser.add_argument(
"--duplicate-coordinate-warning",
type=float,
default=0.05,
metavar="FRACTION",
help="Warn when rounded duplicate coordinate fraction exceeds this value (default: 0.05)",
)
parser.add_argument("--json", action="store_true", help="Print machine-readable JSON")
args = parser.parse_args()
result: dict[str, Any] = {
"path": str(args.h5ad.expanduser().resolve()),
"summary": {},
"warnings": [],
"errors": [],
}
h5ad_path = args.h5ad.expanduser()
if not h5ad_path.is_file():
result["errors"].append(f"file does not exist: {h5ad_path}")
if args.json:
print(json.dumps(result, ensure_ascii=False, indent=2))
else:
_print_human(result)
return 2
try:
import anndata as ad
import numpy as np
except ImportError as exc:
result["errors"].append(
f"cannot import anndata/numpy ({exc}); install requirements_python.txt first"
)
if args.json:
print(json.dumps(result, ensure_ascii=False, indent=2))
else:
_print_human(result)
return 2
adata = None
try:
# Backed mode avoids reading the full expression matrix into memory.
adata = ad.read_h5ad(h5ad_path, backed="r")
n_obs, n_vars = adata.n_obs, adata.n_vars
obsm_keys = list(adata.obsm_keys())
layer_keys = list(adata.layers.keys())
result["summary"].update(
{
"n_obs": int(n_obs),
"n_vars": int(n_vars),
"obsm_keys": obsm_keys,
"layers": layer_keys,
"raw_present": adata.raw is not None,
}
)
if n_obs == 0 or n_vars == 0:
result["errors"].append("AnnData has zero observations or zero genes")
if not adata.obs_names.is_unique:
result["errors"].append("obs_names are not unique")
if not adata.var_names.is_unique:
result["errors"].append("var_names are not unique")
required_obs = _split_names(args.required_obs)
missing_required = [name for name in required_obs if name not in adata.obs.columns]
if missing_required:
result["errors"].append(
"missing required obs columns: " + ", ".join(missing_required)
)
present_optional = [
name for name in _split_names(args.optional_obs) if name in adata.obs.columns
]
missing_optional = [
name for name in _split_names(args.optional_obs) if name not in adata.obs.columns
]
result["summary"]["present_optional_obs"] = present_optional
if missing_optional:
result["warnings"].append(
"optional obs columns not present: " + ", ".join(missing_optional)
)
if "sample" in adata.obs.columns:
sample_series = adata.obs["sample"]
result["summary"]["n_samples"] = int(sample_series.nunique(dropna=False))
missing_sample = int(sample_series.isna().sum())
result["summary"]["missing_sample_values"] = missing_sample
if missing_sample:
result["errors"].append(f"sample has {missing_sample} missing values")
spatial_key = args.spatial_key
if spatial_key not in obsm_keys:
result["errors"].append(f"missing adata.obsm[{spatial_key!r}]")
else:
try:
coords = np.asarray(adata.obsm[spatial_key])
if coords.ndim != 2 or coords.shape[0] != n_obs or coords.shape[1] < 2:
result["errors"].append(
f"adata.obsm[{spatial_key!r}] must have shape (n_obs, >=2), "
f"got {coords.shape}"
)
else:
xy = coords[:, :2].astype(float)
if not np.isfinite(xy).all():
result["errors"].append(
f"adata.obsm[{spatial_key!r}] contains non-finite coordinates"
)
else:
rounded = np.round(xy, decimals=9)
unique_fraction = np.unique(rounded, axis=0).shape[0] / max(n_obs, 1)
duplicate_fraction = 1.0 - unique_fraction
result["summary"].update(
{
"spatial_shape": list(coords.shape),
"spatial_duplicate_fraction_rounded_1e-9": round(
float(duplicate_fraction), 6
),
"spatial_min": [float(xy[:, 0].min()), float(xy[:, 1].min())],
"spatial_max": [float(xy[:, 0].max()), float(xy[:, 1].max())],
}
)
if duplicate_fraction > args.duplicate_coordinate_warning:
result["warnings"].append(
"rounded spatial coordinates have duplicate fraction "
f"{duplicate_fraction:.3f}; verify CellBin centroid/coordinate semantics"
)
except Exception as exc: # noqa: BLE001 - report malformed obsm values clearly
result["errors"].append(
f"cannot read numeric coordinates from adata.obsm[{spatial_key!r}]: {exc}"
)
requested_genes = _split_names(args.check_genes)
if requested_genes:
gene_set = set(adata.var_names.astype(str))
missing_genes = [gene for gene in requested_genes if gene not in gene_set]
result["summary"]["checked_genes"] = requested_genes
result["summary"]["missing_genes"] = missing_genes
if missing_genes:
result["warnings"].append(
"requested genes not found in var_names: " + ", ".join(missing_genes)
)
if "donor" not in adata.obs.columns:
result["warnings"].append(
"donor is absent; add a donor/sample mapping before inferential statistics"
)
if "group" in adata.obs.columns:
groups = adata.obs["group"].dropna().astype(str).unique().tolist()
result["summary"]["group_values"] = sorted(groups)
if not groups:
result["warnings"].append("group exists but has no non-missing values")
if "celltype" in adata.obs.columns:
result["summary"]["celltype_values_n"] = int(
adata.obs["celltype"].nunique(dropna=True)
)
except Exception as exc: # noqa: BLE001 - preserve a useful CLI error
result["errors"].append(f"failed to inspect H5AD: {exc}")
finally:
if adata is not None:
try:
adata.file.close()
except (AttributeError, OSError):
pass
if args.json:
print(json.dumps(result, ensure_ascii=False, indent=2, default=_serializable))
else:
_print_human(result)
return 2 if result["errors"] else 0
if __name__ == "__main__":
sys.exit(main())
