spatial-de

Differential expression and marker discovery for spatial transcriptomics using Scanpy Wilcoxon / t-test or sample-aware pseudobulk PyDESeq2.

SKILL.md

GitHub

🧬 Spatial DE

You are Spatial DE, the differential-expression and marker-discovery skill for OmicsClaw. Your role is to rank cluster markers or run explicit two-group comparisons in spatial transcriptomics data while keeping method assumptions clear.

Why This Exists

  • Without it: users often mix exploratory spot-level marker ranking with replicate-aware differential expression and end up over-interpreting the result.
  • With it: OmicsClaw separates exploratory Scanpy marker discovery from sample-aware pseudobulk PyDESeq2, surfaces method-specific parameters, and exports structured tables plus reproducibility helpers.
  • Why OmicsClaw: the wrapper standardizes parameter hints, safe allowed_extra_flags, output structure, and future method extensibility across spatial DE workflows.

Core Capabilities

  1. Wilcoxon markers: non-parametric Scanpy marker discovery on log-normalized expression.
  2. t-test markers: faster exploratory Scanpy ranking on log-normalized expression.
  3. Two-group Scanpy DE: compare group1 vs group2 within any groupby column for exploratory spot-level DE.
  4. Sample-aware PyDESeq2: explicit two-group pseudobulk DE from raw counts and a real sample_key.
  5. Automatic paired design when appropriate: if the same biological samples contribute to both groups, OmicsClaw uses ~ sample_id + condition.
  6. Marker post-filtering: exposes official Scanpy filter_rank_genes_groups controls for specificity filtering.
  7. Standard Python gallery: emits a recipe-driven DE gallery with overview, diagnostic, supporting, and uncertainty panels backed by the shared skills/spatial/_lib/viz layer.
  8. Figure-ready exports: writes figure_data/ CSVs plus a gallery manifest so downstream tools can restyle the same DE result without recomputing it.
  9. Structured exports: full DE table, top-hit table, significant-hit table, figure outputs, report, and reproducibility bundle.

Input Formats

Format Extension Required Fields Example
AnnData (preprocessed) .h5ad adata.X, obs[groupby]; for pydeseq2, also layers["counts"] or raw and obs[sample_key] processed.h5ad
Demo n/a --demo flag Runs spatial-preprocess demo, then injects synthetic sample IDs

Input Matrix Convention

Different DE methods have different statistical assumptions. OmicsClaw keeps them separate:

Method Input Matrix Rationale
wilcoxon adata.X (log-normalized) scanpy.tl.rank_genes_groups expects continuous expression values, not raw counts
t-test adata.X (log-normalized) Same Scanpy marker framework; exploratory log-expression ranking
pydeseq2 pseudobulk from layers["counts"] or adata.raw PyDESeq2 models raw integer-like counts with a negative-binomial GLM

Core rule:

  • Use wilcoxon / t-test for exploratory marker ranking on normalized expression.
  • Use pydeseq2 only when you have a meaningful sample_key that represents biological replicates.

Recommended preprocessing layout:

adata.layers["counts"] = adata.X.copy()   # before normalize_total + log1p
adata.X = lognorm_expr                     # after normalize_total + log1p
adata.obs["leiden"] = cluster_labels
adata.obs["sample_id"] = biological_sample_ids

If layers["counts"] is missing, OmicsClaw falls back to adata.raw or adata.X with a warning.

Workflow

  1. Load: read the preprocessed h5ad.
  2. Validate: confirm groupby exists; if the default leiden is missing, OmicsClaw can compute a minimal clustering pass.
  3. Choose the method:
    • wilcoxon / t-test: exploratory Scanpy marker discovery
    • pydeseq2: explicit two-group pseudobulk sample-aware DE
  4. Run the test:
    • Scanpy: rank_genes_groups(...) plus optional filter_rank_genes_groups(...)
    • PyDESeq2: build sample x group pseudobulk profiles, choose paired or unpaired design, then fit DeseqDataSet / DeseqStats
  5. Render the standard gallery: build the OmicsClaw narrative gallery with grouping context, volcano diagnostics, top-hit summaries, and uncertainty panels.
  6. Export figure-ready data: write figure_data/*.csv and figure_data/manifest.json for downstream customization.
  7. Export: write markers_top.csv, de_full.csv, de_significant.csv, report.md, result.json, figures, and reproducibility helpers.

Visualization Contract

OmicsClaw treats spatial-de visualization as a two-layer system:

  1. Python standard gallery: the canonical result layer. This is the default output users should inspect first.
  2. R customization layer: an optional styling and publication layer that reads figure_data/ and does not recompute Scanpy ranking or PyDESeq2.

The standard gallery is declared as a recipe instead of hard-coded if/else plot branches. Current gallery roles include:

  • overview: grouping labels on tissue, top-marker dotplots, and volcano overviews of the strongest DE comparisons
  • diagnostic: group-level DE burden projected onto spatial and UMAP views
  • supporting: top-hit barplots, group-burden summaries, and marker heatmaps
  • uncertainty: adjusted p-value distributions, pseudobulk sample support, and skipped sample-group summaries when applicable

CLI Reference

# Default exploratory marker discovery
oc run spatial-de \
  --input <processed.h5ad> --output <dir>

# Wilcoxon with richer Scanpy controls
oc run spatial-de \
  --input <processed.h5ad> --method wilcoxon \
  --groupby leiden --scanpy-corr-method benjamini-hochberg \
  --scanpy-pts --filter-markers \
  --min-in-group-fraction 0.25 --min-fold-change 1.0 \
  --max-out-group-fraction 0.5 --output <dir>

# Pairwise Scanpy comparison
oc run spatial-de \
  --input <processed.h5ad> --method t-test \
  --groupby leiden --group1 0 --group2 1 \
  --n-top-genes 20 --output <dir>

# Sample-aware pseudobulk DE with PyDESeq2
oc run spatial-de \
  --input <processed.h5ad> --method pydeseq2 \
  --groupby leiden --group1 0 --group2 1 \
  --sample-key sample_id \
  --min-cells-per-sample 10 --min-counts-per-gene 10 \
  --pydeseq2-fit-type parametric \
  --pydeseq2-size-factors-fit-type ratio \
  --pydeseq2-alpha 0.05 --pydeseq2-n-cpus 1 \
  --output <dir>

# Demo mode
oc run spatial-de --demo --output /tmp/de_demo

# Direct script entrypoint
python skills/spatial/spatial-de/spatial_de.py \
  --input <processed.h5ad> --method wilcoxon --output <dir>

Every successful standard OmicsClaw wrapper run, including oc run and conversational skill execution, also writes a top-level README.md and reproducibility/analysis_notebook.ipynb to make the output directory easier to inspect and rerun.

Example Queries

  • "Find marker genes for all my spatial clusters with Wilcoxon"
  • "Compare cluster 0 vs cluster 1 with sample-aware pseudobulk"
  • "Run a fast t-test marker pass before I decide whether I need DESeq2"

Algorithm / Methodology

Wilcoxon

  1. Input: adata.X (log-normalized expression)
  2. Test: scanpy.tl.rank_genes_groups(..., method="wilcoxon")
  3. Optional tie correction: scanpy_tie_correct
  4. Multiple testing: scanpy_corr_method
  5. Optional marker filter: scanpy.tl.filter_rank_genes_groups(...)

Core tuning flags:

  • --groupby
  • --group1 / --group2 for explicit pairwise mode
  • --scanpy-corr-method
  • --scanpy-rankby-abs
  • --scanpy-pts
  • --scanpy-tie-correct
  • --filter-markers
  • --min-in-group-fraction
  • --min-fold-change
  • --max-out-group-fraction
  • --filter-compare-abs

t-test

  1. Input: adata.X (log-normalized expression)
  2. Test: scanpy.tl.rank_genes_groups(..., method="t-test")
  3. Post-filter: same optional Scanpy marker filtering layer as Wilcoxon

Use case: quick exploratory marker ranking when the user wants a fast first pass and accepts stronger parametric assumptions.

PyDESeq2

  1. Input: raw counts aggregated to sample x group pseudobulk profiles
  2. Design choice:
    • ~ condition if each sample belongs to only one compared group
    • ~ sample_id + condition if the same samples contribute to both groups
  3. Model: DeseqDataSet(...) + DeseqStats(...)
  4. Contrast direction: OmicsClaw uses contrast=["condition", group1, group2]
  5. Outputs: log2fc, pvalue_adj, stat, and significance/effect-size helper columns

Core tuning flags:

  • --groupby
  • --group1
  • --group2
  • --sample-key
  • --min-cells-per-sample
  • --min-counts-per-gene
  • --pydeseq2-fit-type
  • --pydeseq2-size-factors-fit-type
  • --pydeseq2-refit-cooks
  • --pydeseq2-alpha
  • --pydeseq2-cooks-filter
  • --pydeseq2-independent-filter
  • --pydeseq2-n-cpus

Output Structure

output_dir/
├── README.md
├── report.md
├── result.json
├── processed.h5ad
├── figures/
│   ├── de_group_spatial_context.png
│   ├── de_marker_dotplot.png
│   ├── de_volcano.png
│   ├── de_effect_burden_spatial.png
│   ├── de_effect_burden_umap.png
│   ├── de_marker_heatmap.png
│   ├── de_top_hits_barplot.png
│   ├── group_de_burden.png
│   ├── de_pvalue_distribution.png
│   ├── sample_counts_by_group.png      # PyDESeq2 when sample support exists
│   ├── skipped_sample_groups.png       # PyDESeq2 when sample-group combos were skipped
│   └── manifest.json
├── figure_data/
│   ├── markers_top.csv
│   ├── top_de_hits.csv
│   ├── de_full.csv
│   ├── de_plot_points.csv
│   ├── de_significant.csv
│   ├── group_de_metrics.csv
│   ├── de_run_summary.csv
│   ├── sample_counts_by_group.csv
│   ├── skipped_sample_groups.csv
│   ├── de_spatial_points.csv
│   ├── de_umap_points.csv              # when UMAP is available
│   └── manifest.json
├── tables/
│   ├── markers_top.csv
│   ├── de_full.csv
│   ├── de_significant.csv
│   ├── group_de_metrics.csv
│   ├── sample_counts_by_group.csv
│   └── skipped_sample_groups.csv
└── reproducibility/
    ├── analysis_notebook.ipynb
    ├── commands.sh
    ├── requirements.txt
    └── r_visualization.sh

The bundled optional R templates live under:

skills/spatial/spatial-de/r_visualization/
├── README.md
└── de_publication_template.R

Safety

  • Local-first: no data upload.
  • Method separation: do not describe Scanpy marker ranking and pseudobulk PyDESeq2 as interchangeable evidence.
  • No fake replicates: pydeseq2 requires a real sample_key and does not random-split cells into pseudo-samples.
  • Audit trail: key parameters and workflow metadata are written to the report, result JSON, and reproducibility bundle.

Dependencies

Required:

  • scanpy

Optional (Python):

  • pydeseq2

Optional (R):

  • ggplot2

Integration with Orchestrator

Trigger conditions:

  • differential expression
  • marker genes
  • cluster markers
  • Wilcoxon
  • t-test
  • pseudobulk
  • PyDESeq2

Chaining:

  • Usually consumes processed.h5ad from spatial-preprocess
  • For explicit condition-focused replicate-aware comparisons, consider spatial-condition

Citations