Expand 10x renal dataset smoke test workflow

Author: hutaobo

README.md

@@ -59,6 +59,22 @@ python examples/smoke_test_10x_renal_ffpe_protein.py \
   "Y:/long/10X_datasets/Xenium/Xenium_Renal/Xenium_V1_Human_Kidney_FFPE_Protein"
 ```
 
+To also write a compact Markdown/JSON/CSV report bundle:
+
+```bash
+python examples/smoke_test_10x_renal_ffpe_protein.py \
+  "Y:/long/10X_datasets/Xenium/Xenium_Renal/Xenium_V1_Human_Kidney_FFPE_Protein" \
+  --output-dir ./smoke_test_outputs
+```
+
+To export the loaded object for downstream analysis:
+
+```bash
+python examples/smoke_test_10x_renal_ffpe_protein.py \
+  "Y:/long/10X_datasets/Xenium/Xenium_Renal/Xenium_V1_Human_Kidney_FFPE_Protein" \
+  --write-h5ad ./renal_ffpe_protein.h5ad
+```
+
 Quick Start
 -----------
 

examples/smoke_test_10x_renal_ffpe_protein.py

@@ -6,6 +6,7 @@
 import sys
 from pathlib import Path
 
+import numpy as np
 import pandas as pd
 
 try:
@@ -27,7 +28,54 @@
 EXPECTED_PROTEIN_MARKERS = 27
 
 
-def build_summary(base_path: str, prefer: str) -> dict:
+def _top_rna_features(adata, top_n: int) -> list[dict]:
+    feature_names = (
+        adata.var["name"].astype(str).tolist()
+        if "name" in adata.var.columns
+        else adata.var_names.astype(str).tolist()
+    )
+    nnz = np.asarray(adata.X.getnnz(axis=0)).ravel()
+    total = np.asarray(adata.X.sum(axis=0)).ravel()
+
+    order = np.argsort(-total)[:top_n]
+    rows = []
+    for idx in order:
+        rows.append(
+            {
+                "feature": feature_names[idx],
+                "detected_cells": int(nnz[idx]),
+                "total_counts": float(total[idx]),
+            }
+        )
+    return rows
+
+
+def _top_protein_markers(adata, top_n: int) -> list[dict]:
+    protein = adata.obsm["protein"]
+    protein_df = protein if isinstance(protein, pd.DataFrame) else pd.DataFrame(protein, index=adata.obs_names)
+
+    mean_signal = protein_df.mean(axis=0).sort_values(ascending=False)
+    rows = []
+    for marker, value in mean_signal.head(top_n).items():
+        rows.append(
+            {
+                "marker": str(marker),
+                "mean_signal": float(value),
+                "positive_cells": int((protein_df[marker] > 0).sum()),
+            }
+        )
+    return rows
+
+
+def _top_clusters(adata, top_n: int) -> list[dict]:
+    if "cluster" not in adata.obs.columns:
+        return []
+
+    counts = adata.obs["cluster"].astype(str).value_counts().head(top_n)
+    return [{"cluster": str(cluster), "n_cells": int(count)} for cluster, count in counts.items()]
+
+
+def build_summary(base_path: str, prefer: str, top_n: int = 10) -> tuple[dict, object]:
     adata = load_xenium_gene_protein(base_path=base_path, prefer=prefer)
 
     protein = adata.obsm.get("protein")
@@ -47,6 +95,9 @@ def build_summary(base_path: str, prefer: str) -> dict:
         "has_cluster": "cluster" in adata.obs.columns,
         "obsm_keys": sorted(adata.obsm.keys()),
         "metrics_summary_num_cells_detected": None,
+        "top_rna_features_by_total_counts": _top_rna_features(adata, top_n=top_n),
+        "top_protein_markers_by_mean_signal": _top_protein_markers(adata, top_n=top_n),
+        "largest_clusters": _top_clusters(adata, top_n=top_n),
     }
 
     metrics_path = Path(base_path) / "metrics_summary.csv"
@@ -55,7 +106,83 @@ def build_summary(base_path: str, prefer: str) -> dict:
         if "num_cells_detected" in metrics.columns and not metrics.empty:
             summary["metrics_summary_num_cells_detected"] = int(metrics.loc[0, "num_cells_detected"])
 
-    return summary
+    return summary, adata
+
+
+def render_markdown_report(payload: dict) -> str:
+    summary = payload["summary"]
+    validated = payload["validated_reference"]
+    issues = payload["issues"]
+
+    lines = [
+        "# pyXenium Smoke Test Report",
+        "",
+        f"Dataset: {summary['dataset_title']}",
+        f"Source: {summary['dataset_url']}",
+        f"Local path: `{summary['base_path']}`",
+        f"Backend preference: `{summary['prefer']}`",
+        "",
+        "## Core Results",
+        "",
+        f"- Cells: `{summary['n_cells']}`",
+        f"- RNA features: `{summary['n_rna_features']}`",
+        f"- Protein markers: `{summary['n_protein_markers']}`",
+        f"- Sparse matrix nnz: `{summary['x_nnz']}`",
+        f"- Spatial coordinates present: `{summary['has_spatial']}`",
+        f"- Cluster labels present: `{summary['has_cluster']}`",
+        f"- metrics_summary.csv detected cells: `{summary['metrics_summary_num_cells_detected']}`",
+        "",
+        "## Validated Reference",
+        "",
+        f"- Expected cells: `{validated['expected_cells']}`",
+        f"- Expected RNA features: `{validated['expected_rna_features']}`",
+        f"- Expected protein markers: `{validated['expected_protein_markers']}`",
+        "",
+        "## Largest Clusters",
+        "",
+    ]
+
+    for row in summary["largest_clusters"]:
+        lines.append(f"- `{row['cluster']}`: `{row['n_cells']}` cells")
+
+    lines.extend(["", "## Top RNA Features by Total Counts", ""])
+    for row in summary["top_rna_features_by_total_counts"]:
+        lines.append(
+            f"- `{row['feature']}`: total counts `{row['total_counts']:.0f}`, detected cells `{row['detected_cells']}`"
+        )
+
+    lines.extend(["", "## Top Protein Markers by Mean Signal", ""])
+    for row in summary["top_protein_markers_by_mean_signal"]:
+        lines.append(
+            f"- `{row['marker']}`: mean signal `{row['mean_signal']:.4f}`, positive cells `{row['positive_cells']}`"
+        )
+
+    lines.extend(["", "## Issues", ""])
+    if issues:
+        lines.extend(f"- {issue}" for issue in issues)
+    else:
+        lines.append("- No issues detected.")
+
+    lines.append("")
+    return "\n".join(lines)
+
+
+def write_output_artifacts(payload: dict, output_dir: str | None) -> None:
+    if not output_dir:
+        return
+
+    out = Path(output_dir)
+    out.mkdir(parents=True, exist_ok=True)
+
+    (out / "summary.json").write_text(json.dumps(payload, indent=2) + "\n", encoding="utf-8")
+    (out / "report.md").write_text(render_markdown_report(payload), encoding="utf-8")
+    pd.DataFrame(payload["summary"]["top_rna_features_by_total_counts"]).to_csv(
+        out / "top_rna_features.csv", index=False
+    )
+    pd.DataFrame(payload["summary"]["top_protein_markers_by_mean_signal"]).to_csv(
+        out / "top_protein_markers.csv", index=False
+    )
+    pd.DataFrame(payload["summary"]["largest_clusters"]).to_csv(out / "largest_clusters.csv", index=False)
 
 
 def validate_summary(summary: dict) -> list[str]:
@@ -108,6 +235,12 @@ def parse_args() -> argparse.Namespace:
         default="auto",
         help="Preferred matrix backend passed to load_xenium_gene_protein().",
     )
+    parser.add_argument(
+        "--top-n",
+        type=int,
+        default=10,
+        help="Number of top RNA features, protein markers, and clusters to report.",
+    )
     parser.add_argument(
         "--allow-mismatch",
         action="store_true",
@@ -118,12 +251,22 @@ def parse_args() -> argparse.Namespace:
         default=None,
         help="Optional path to write the summary JSON.",
     )
+    parser.add_argument(
+        "--output-dir",
+        default=None,
+        help="Optional directory for report.md, summary.json, and CSV summaries.",
+    )
+    parser.add_argument(
+        "--write-h5ad",
+        default=None,
+        help="Optional path to export the loaded AnnData object as an .h5ad file.",
+    )
     return parser.parse_args()
 
 
 def main() -> int:
     args = parse_args()
-    summary = build_summary(base_path=args.base_path, prefer=args.prefer)
+    summary, adata = build_summary(base_path=args.base_path, prefer=args.prefer, top_n=args.top_n)
     issues = validate_summary(summary)
 
     payload = {
@@ -141,6 +284,12 @@ def main() -> int:
 
     if args.output_json:
         Path(args.output_json).write_text(rendered + "\n", encoding="utf-8")
+    write_output_artifacts(payload, args.output_dir)
+
+    if args.write_h5ad:
+        h5ad_path = Path(args.write_h5ad)
+        h5ad_path.parent.mkdir(parents=True, exist_ok=True)
+        adata.write_h5ad(h5ad_path)
 
     if issues and not args.allow_mismatch:
         return 1