differential_expression

Files (2)

• PROTOCOL.md
• SKILL.md

SKILL.md (first lines)

## Read this before planning step 3

A condition contrast is a statement about the CONDITION, so the unit of replication is the
SAMPLE (donor / animal / library), never the cell. Cells from one donor are not independent
observations of that donor's condition. A Wilcoxon test over cells treats them as if they
were, and the p-values are then anti-conservative by orders of magnitude.

This is measured, not theoretical. On a synthetic 4-donor design where **exactly 1 of 400
genes** was made different between the arms:

| test | called significant (padj<0.05) |
|---|---|
| `run_de`, Wilcoxon over cells | **310 of 400 (78%)** |
| `run_pseudobulk_de`, over donors | 0 (and it ranked the true gene #1) |

So: **use `run_pseudobulk_de` for a condition contrast whenever the data has ≥2 samples per
arm.** `run_de` at its default (`reference="rest"`) is for markers — one cluster versus the rest
within a sample — which is what it is valid for. `run_de` can also run a real cell-level contrast
(`reference=<control level>`, `stratify_by=<cell-type column>`), but that inherits the
pseudoreplication above: use it only when there are no replicates, and label it as an exploratory
ranking rather than differential expression.

## Read this before planning the contrast, too: SEQUENCING DEPTH

Replication is not the only way a condition contrast goes wrong. If one arm was sequenced deeper
per cell, EVERY comparison inherits it: a gene detected in more cells because more molecules were
sampled scores as up-regulated, and whole pathways of abundant transcripts — translation,
ribosome, RNA metabolism — move together in the deeper arm. Log-normalisation does NOT remove
this; it rescales totals, while detection is a function of how many molecules were drawn.

The dataset profile gives you the median counts per cell for each arm. **Compare them before you
plan.** If they differ by more than ~1.3x, the plan MUST contain a depth-matched check, because
with one library per arm there is no replication that could separate depth from biology.

`run_depth_matched_de` is that check as a deterministic tool: within each cell type it down-samples
the deeper arm quantile by quantile until the per-cell UMI distributions match (never up-sampling),
re-runs the SAME contrast, and Spearman-correlates the original ranking against the matched one,
separately per direction. High rho = not an artefact. Rho near zero = depth decided the ranking.
NEGATIVE rho = the ranking inverts once depth is equal, which is the strongest evidence that the
finding must not be reported as biology.

Name that tool. Do NOT plan this as free-form `run_code`: it existed only as a step brief before,
and three different models wrote three different wrong versions — one plan asked the executor to
"correlate per-gene logFC against the between-arm median nCount_RNA", which is not computable (a
per-gene vector against a single scalar), and the step failed on every attempt.

Condition-comparison protocol. Use when the goal is to find the genes that change between two
experimental GROUPS of cells — disease vs control, knockout vs wild-type, treated vs untreated —
and interpret them. NOT for assigning cell-type labels (use `celltype_annotation` for that).

## Figuring out WHAT to compare from the data (do this even with a vague/empty question)

The user may not know the biology, or may just say "analyze this dataset". You do NOT need them to
spell it out — the dataset's own metadata defines the experiment. Read the DATASET PROFILE (the
obs columns + their category values are given to you at planning time) and infer:

1. **The condition/group column.** Look in obs for a low-cardinality (usually 2-level) column that
   names an experimental variable: `sampleid`, `condition`, `genotype`, `treatment`, `group`,
   `orig.ident`, `disease`, etc. Its two values ARE the comparison (e.g. `sampleid=[DDX41, WT]`
   → compare DDX41 vs WT). The dataset / file name often confirms intent (e.g. `*_DEG` = they
   want differential expression).
2. **The reference (control) group.** Pick the level whose name looks like the baseline —
   `WT`, `wild-type`, `control`, `ctrl`, `untreated`, `vehicle`, `DMSO`, `normal`, `sham`. The
   OTHER level is the condition of interest. State which you chose; if it is genuinely unclear,
   ask ONE clarify question in plan mode rather than guessing.
3. **The cell-type stratification.** If obs already has a cell-type / cluster label column
   (`majorclass`, `celltype`, a predicted-label column), REUSE it — run the comparison
   SEPARATELY WITHIN EACH cell type. Do NOT re-cluster or re-annotate from scratch. Only if no
   label column exists do you `run_clustering` first and compare across clusters.

So the default study for a dataset with a 2-group condition column + existing cell-type labels is:
**"condition vs control differential expression within each major cell type, then pathway
enrichment on the changed genes"** — plan that without waiting for the user to describe it.

## Ordered plan

0. **Establish what the matrix IS, and the rules, before touching it** (`run_code`) — read-only,
   and BEFORE QC. Two outputs, written as tables the later steps read back:
   * *Provenance.* **Not for a 10x Cell Ranger input whose data profile already rules on it**
     (its `.h5` matrices were read at upload: integer UMI counts, checked against Cell Ranger's
     `metrics_summary.csv`). Cite that ruling; do NOT plan a `run_code` step to re-derive it, and
     write the decision rule below into the plan's first step instead. Otherwise: inspect `X`,
     `.raw`, `layers`, value ranges, sparsity and integer-ness, and check whether matrix-derived
     totals reproduce any stored `nCount`/`nFeature` fields. A
     published object often arrives already normalized and log1p'd: normalizing it again has NO
     symptom — every tool succeeds, every figure renders, and every number after it is wrong. So
     rule from the NUMBERS, and say which way you ruled:
     - **Counts.** Every value non-negative and integer-valued to tolerance, the row sums
       reproduce a stored total-count field (`nCount_RNA`, `total_counts`), the per-row non-zero
       counts reproduce the stored detected-feature field, and no `layers`/`.raw` entry carries a
       different scale. **That is sufficient — proceed.** An empty `uns` and a missing
       transformation history are a LIMITATION TO REPORT, not a stop condition: a matrix that
       reproduces its own stored QC fields exactly is not made unusable by the absence of a
       written record, and stopping there forfeits the whole quantitative analysis over a
       provenance note.
     - **Not counts.** Any negative or non-integer value, row sums that do not reproduce the
       stored totals, or a differently-scaled representation present. Then say so and drop the
       count-dependent steps — never recover counts by exponentiating, rounding, or renaming.
     Count-based tests (`run_pseudobulk_de`, `run_depth_matched_de`) need the first ruling; say
     which one you applied and which representation supports which operation.
   * *The decision rule.* Name the primary estimand, the unit of replication, the minimum cells
     per arm per cell type, and the effect-size gate — and write them down NOW, while no DE
     result is visible. A cutoff chosen after seeing which one flatters the answer is not a
     cutoff. Say what happens when a gate FAILS: that stratum is reported unsupported, not
     relaxed until it passes. Most of these numbers already have defaults in the Parameters
     table below — this step is about COMMITTING to them (and to any change, with its reason)
     before the results can influence the choice.
1. **QC** (`run_scanpy_qc` for one matrix; `run_cellqc` INSTEAD when the data profile says the input is a folder of 10x Cell Ranger outputs (raw + filtered matrices): it corrects ambient RNA, removes doublets and writes the same normalised checkpoint): per-cell metrics, filter, normalize + log1p (+ HVG). Report counts.
   Honor any QC columns already in the data (e.g. `percent.mt`, doublet calls).
2. **Define the comparison** from the profile (above): name the condition column, the two groups,
   the reference group, and the cell-type column you will stratify by. If labels already exist,
   skip clustering.
3. **Per-cell-type differential expression — choose the test from the DESIGN, not from habit.**
   First find the SAMPLE column (donor / animal / library / `orig.ident`), which is usually a
   different column from the condition. Then:

   - **≥2 samples per arm → `run_pseudobulk_de`** with `sample_key`, `condition_key`, and
     `group_key` = the cell-type column. It aggregates counts per sample, tests across
     samples, and returns `skipped_groups` for any cell type without enough samples. Report

From RCHENLAB/AiScientist, MIT.

More from RCHENLAB/AiScientist