sc-de

When to use

The user has a preprocessed scRNA-seq AnnData and wants to know either (a) which genes mark each cluster (Wilcoxon / t-test / logreg ranking) or (b) which genes change between conditions in a replicate-aware way (deseq2_r pseudobulk). Five backends are exposed; the wrapper enforces the matrix contract (normalized expression vs raw counts) per backend because mixing them is the most common silent-wrong-answer failure mode.

Inputs & Outputs

<!-- AUTO-GENERATED from skill.yaml (interface) — do not edit by hand. Regenerate: python scripts/generate_skill_md.py <skill_dir> -->

Inputs

• Modalities: scrna
• File types: .h5ad
• Requires a preprocessed AnnData (X normalised, PCA/neighbours present)

Outputs

• tables/counts.csv
• tables/de_full.csv
• tables/de_group_summary.csv
• tables/de_top_markers.csv
• tables/deseq2_results.csv
• tables/gene_expression.csv
• tables/markers_top.csv
• tables/mast_results.csv
• tables/metadata.csv
• tables/pseudobulk_summary.csv
• figures/marker_dotplot.png
• figures/pseudobulk_group_summary.png
• figures/r_de_heatmap.png
• figures/r_de_manhattan.png
• figures/r_de_volcano.png
• figures/r_feature_cor.png
• figures/r_feature_violin.png
• figures/rank_genes_groups.png
• input.h5ad
• processed.h5ad
• report.md
• result.json
• Processed AnnData (saves_h5ad)

Flow

1. Load the AnnData and inspect which matrix is in X (normalized vs raw).
2. Validate the requested method's matrix contract — fail fast on mismatch.
3. For exploratory paths: run Scanpy rank_genes_groups and export marker dotplot + rank summary.
4. For mast: hand the log-normalized matrix to the R MAST bridge.
5. For deseq2_r: pseudobulk-aggregate by sample_key × celltype_key, drop bins under --pseudobulk-min-cells / --pseudobulk-min-counts, run R DESeq2.
6. Render direct figures (and R-enhanced ones if --r-enhanced).
7. Write processed.h5ad, tables, figure_data/manifest.json, report.md, result.json, and reproducibility script.

Gotchas

• Exploratory ranking paths are NOT replicate-aware. Wilcoxon / t-test / logreg / mast treat each cell as an independent observation, which inflates Type I error for treated-vs-control comparisons across samples. Use deseq2_r whenever biological replicates exist and the question is condition-level inference — not cluster markers.
• deseq2_r requires raw counts. It looks for layers["counts"] first, then adata.raw, then adata.X, gated by a matrix_looks_count_like heuristic on each candidate. After every pseudobulk run, check result.json["summary"]["expression_source"] — it should read layers.counts or adata.raw, not adata.X. If the heuristic mis-classifies a normalized matrix as count-like, the run will fall through to adata.X silently.
• --group1 and --group2 are required for deseq2_r. sc_de.py raises ValueError when either is missing on the pseudobulk path. Other methods treat them as optional (cluster-vs-rest if absent).
• mast and deseq2_r need a working R stack. mast needs MAST + SingleCellExperiment + zellkonverter; deseq2_r needs DESeq2 + same companions. The wrapper raises if the R bridge is unavailable — install via the project's R bootstrap, not pip.
• Pseudobulk bin filtering silently drops entire cell types. Per skills/singlecell/_lib/pseudobulk.py:103, any sample × cell type bin with fewer than --pseudobulk-min-cells (default 10) cells is continue'd with no diagnostic written to result.json. When DESeq2 reports "no DEGs in cell type X," sanity-check the per-sample cell counts manually before assuming biological null.
• sample_key is statistical design, not a label. DESeq2 fits a per-sample dispersion; using a non-replicate column (e.g. cell_type itself) gives nonsense. The wrapper does not currently catch this — sanity-check the column has >=2 distinct values per condition.

Key CLI

# Demo: PBMC3k Wilcoxon cluster markers
python omicsclaw.py run sc-de --demo

# Exploratory: cluster markers via Wilcoxon
python omicsclaw.py run sc-de \
  --input processed.h5ad --output results/ \
  --groupby leiden --method wilcoxon --n-top-genes 20

# Replicate-aware: treated vs control via DESeq2 pseudobulk
python omicsclaw.py run sc-de \
  --input processed.h5ad --output results/ \
  --method deseq2_r --groupby condition \
  --group1 treated --group2 control \
  --sample-key sample_id --celltype-key cell_type \
  --pseudobulk-min-cells 10 --pseudobulk-min-counts 1000

See also

• references/parameters.md — every CLI flag and per-method tuning hint
• references/methodology.md — Five DE paths, scope boundary, input expectations, workflow
• references/output_contract.md — exact output directory layout + visualization contract
• references/r_visualization.md — five R-enhanced renderers
• Adjacent skills: sc-clustering (upstream cluster discovery), sc-cell-annotation (upstream cell type labels for celltype_key), sc-markers (lighter cluster-marker-only path), sc-enrichment (downstream pathway enrichment of DEG lists), bulkrna-de / spatial-de (sibling DE skills for the other two data modalities)