Sergio-ote · Sep 9, 2026 · Sep 8, 2026 · Sep 8, 2026 · Sep 8, 2026 · Sep 8, 2026
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
@@ -41,7 +41,7 @@ standards.
 ```bash
 git clone https://github.com/cbib/SpaceBlocks && cd SpaceBlocks
 # Install the development environment
-conda create -f workflow/envs/dev.yaml
+conda env create -f workflow/envs/dev.yaml
 conda activate spaceblocks_dev
 snakemake -n --sdm conda      # dry-run: builds the DAG, validates the config, provisions envs
 snakemake -s workflow/Snakefile -d .test -n --workflow-profile none  # decoupled smoke test
@@ -99,20 +99,23 @@ changes. See the `xenium5k` head as a reference implementation.
 
 ## Validating a change
 
-Run the applicable checks locally before pushing. The first two mirror CI and use the committed
-`mode: decoupled` fixture under `.test/`:
+Run the applicable checks locally before pushing. The lint and dry-run commands mirror CI and use
+the committed `mode: decoupled` fixture under `.test/`:
 
 ```bash
-# 1. Workflow parses + lints (against the committed test fixture)
-snakemake -s workflow/Snakefile -d .test --lint
+# 1. Formatting and file-hygiene hooks pass
+pre-commit run --all-files
+
+# 2. Workflow parses + lints (against the committed test fixture)
+snakemake -s workflow/Snakefile -d .test --lint --workflow-profile none
 
-# 2. The decoupled CoreBlock DAG builds
+# 3. The decoupled CoreBlock DAG builds
 snakemake -s workflow/Snakefile -d .test -n --workflow-profile none
 
-# 3. Config still validates against the schema
+# 4. Config still validates against the schema
 python -c "import yaml,jsonschema; jsonschema.validate(yaml.safe_load(open('config/config.yaml')), yaml.safe_load(open('workflow/schemas/config.schema.yaml'))); print('Configuration schema validation passed')"
 
-# 4. Docs build cleanly (only if you touched docs/)
+# 5. Docs build cleanly (only if you touched docs/)
 # Activate the development environment
 conda activate spaceblocks_dev
 # Build the docs

diff --git a/config/README.md b/config/README.md
@@ -104,7 +104,8 @@ Parameters are settings that live exclusively in `config/config.yaml` and determ
 | `use_precomputed_clusters` | **Optional** (reproducibility) | If `true`, reuse Leiden clusters/metadata from `precomputed_metadata_dir` instead of recomputing them. |
 | `ingest_ref_label_key` | `ingest_ref` reference set | Column in the `ingest_ref` reference that holds the reference cell-type labels. |
 | `integration.integrate_key` | **Mandatory** | Variable Harmony corrects over during integration (e.g. `sample`). |
-| `extra_annotations.columns` | **Optional** | `core_samples.tsv` columns to carry into `obs` and surface in downstream plots (e.g. `[patient, batch]`). |
+| `extra_annotations.columns` | **Optional** | `core_samples.tsv` columns to surface in downstream plots (e.g. `[patient, batch]`). All sample-sheet columns are carried into `obs` and can be used by pseudobulk models whether or not they are plotted. |
+| `analysis.pseudobulk.analyses` | Pseudobulk | Named aggregation/model specifications containing the comparison variables, explicit contrasts, exclusions, optional pairing/covariates, and optional LRT. |
 
 ## 4. Color scale customization
 
@@ -132,8 +133,11 @@ annotation_colors:
   cell_type_external:          # declare explicitly; not inherited from cell_type_tsv
     "Tcells": "#1f77b4"
 
-# Region palette — lives under analysis, keyed by the region_levels
+# Preferred region order in QC-sweep, subcluster, and region-only pseudobulk split heatmaps.
+# Entries are shown first-to-last (typically left-to-right); unlisted levels follow.
+# region_colors sets the palette. Pseudobulk filtering uses exclude_levels.
 analysis:
+  region_levels: ["Tumor area", "Healthy area"]
   region_colors:
     "Tumor area": "#46337EFF"
     "Healthy area": "#FDE725FF"
@@ -151,7 +155,7 @@ The rest of the configuration lives in nested blocks. Files and single parameter
 | `atera` | `mode: atera` | Atera head settings: `atera_dir`, Zarr and pyramid options, plus optional registered H&E image, alignment, and keypoint patterns. |
 | `merscope` | `mode: merscope` | MERSCOPE head settings: `merscope_dir`, selected z-plane, embedded-image resolution, and image channels. |
 | `contract` | **Mandatory** | Semantic keys of the hand-off object: `sample_key`, `spatial_key`, `require_region`, `require_raw_counts`, `mito_prefix`. |
-| `analysis` | **Mandatory** | The bulk of the run: QC filters (`min_counts` / `min_genes` / `min_cells` / `max_counts` / `max_pct_mt`), the Leiden `resolution_scan_*`, the pseudobulk `analysis_levels` + thresholds, the `region_levels`, and the `run_*` toggles. Per-sample QC overrides come from `per_sample_qc`. |
+| `analysis` | **Mandatory** | The bulk of the run: QC filters (`min_counts` / `min_genes` / `min_cells` / `max_counts` / `max_pct_mt`), the Leiden `resolution_scan_*`, named pseudobulk analyses + thresholds, and the `run_*` toggles. Per-sample QC overrides come from `per_sample_qc`. |
 | `resources` | **Mandatory** | Per-rule `mem_mb` / `runtime` / `threads` (with a `default`). Memory scales with the retry attempt, so an OOM-killed job is resubmitted with more RAM. |
 | `external_annotation` | *optional* | Overlay labels from an external tool: `enabled`, `column`, `keep_unannotated` (see [section 6](#6-advanced-reproducibility-and-reusability)). |
 | `qc_sweep` | *optional* | Candidate-threshold QC diagnostics (never filters). |
@@ -162,6 +166,153 @@ The rest of the configuration lives in nested blocks. Files and single parameter
 !!! note
     A few one-key blocks are documented elsewhere for readability: `integration` (`integrate_key`) and `extra_annotations` (`columns`) are parameters in [section 3](#3-parameters); `cluster_annotations` is a file in [section 2](#2-paths-files-and-sample-sheets); and the colour blocks (`sample_colors`, `annotation_colors`, `analysis.region_colors`) are in [section 4](#4-color-scale-customization).
 
+### Pseudobulk experimental designs
+
+Pseudobulk analyses are configured under `analysis.pseudobulk.analyses`. Each entry
+has a stable `name`, one aggregation strategy, the sample metadata that define the
+comparison groups, and the statistical tests to run. Counts are summed within each
+pseudobulk observation described in the table below.
+
+The available aggregations are:
+
+| `aggregation` | One pseudobulk observation per | Separate result set per |
+| --- | --- | --- |
+| `all_cells` | sample | — |
+| `by_celltype` | sample and cell type | cell type |
+| `by_region` | sample and region | — |
+| `by_celltype_region` | sample, cell type and region | cell type |
+
+`cell_type` is a portable alias: SpaceBlocks resolves it to the cell-type column for
+the active `annotation_types` entry, so manual TSV, external, ingest-transferred, and
+refined cell-type annotations can use the same pseudobulk configuration.
+
+#### Independent phenotype comparison
+
+Add the experimental variables to `core_samples.tsv`:
+
+```tsv
+sample	phenotype	batch
+Cancer_P1	Cancer	Batch1
+Cancer_P2	Cancer	Batch2
+Cancer_P3	Cancer	Batch3
+Normal_P4	Normal	Batch1
+Normal_P5	Normal	Batch2
+Normal_P6	Normal	Batch3
+```
+
+Then compare phenotypes for each cell type while adjusting for batch:
+
+```yaml
+analysis:
+  run_pseudobulk_de: true
+  min_replicates: 3
+  pseudobulk:
+    analyses:
+      - name: phenotype_by_celltype
+        aggregation: by_celltype
+        group_by: [phenotype]
+        covariates: [batch]
+        exclude_levels:
+          cell_type: [Unannotated]
+        contrasts:
+          - name: cancer_vs_normal
+            numerator: {phenotype: Cancer}
+            denominator: {phenotype: Normal}
+        lrt:
+          enabled: false
+```
+
+The numerator determines the positive log2-fold-change direction. Levels not named
+in a pairwise contrast can remain in the fitted model and help dispersion estimation.
+Use `exclude_levels` only for values that should be removed from the analysis entirely.
+The exclusions are also applied to the LRT.
+
+In this example, P1–P6 are six independent biological samples: no individual
+contributes to both phenotypes, so `paired_by` must be omitted. The batch covariate is
+valid because every batch contains one Cancer and one Normal sample; omit `covariates`
+when there is no nuisance variable to adjust for. Use `aggregation: all_cells` instead
+if one whole-sample result is wanted rather than a separate result for each cell type.
+
+#### Phenotype and treatment
+
+Multiple `group_by` columns form a combined categorical group. In the example below,
+Cancer and Normal phenotype levels are combined with Drug and Vehicle treatment
+levels. This supports clear comparisons between observed combinations without exposing
+interaction formulas:
+
+```yaml
+- name: treatment_by_celltype
+  aggregation: by_celltype
+  group_by: [phenotype, treatment]
+  exclude_levels:
+    cell_type: [Unannotated]
+  contrasts:
+    - name: drug_vs_vehicle_in_cancer
+      numerator: {phenotype: Cancer, treatment: Drug}
+      denominator: {phenotype: Cancer, treatment: Vehicle}
+```
+
+This tests one combination against another; it is **not** a formal statistical
+interaction or difference-of-differences test.
+
+#### Paired regions or matched samples
+
+Use `paired_by: sample` when the same sample contributes both levels of a contrast,
+as in tumour and healthy regions from the same section:
+
+```yaml
+- name: regions_by_celltype
+  aggregation: by_celltype_region
+  group_by: [region_annotation]
+  paired_by: sample
+  exclude_levels:
+    region_annotation: [Unlabeled, Bubble]
+    cell_type: [Unannotated]
+  contrasts:
+    - name: tumor_vs_healthy
+      numerator: {region_annotation: Tumor area}
+      denominator: {region_annotation: Healthy area}
+  lrt:
+    enabled: true
+```
+
+For a paired Wald contrast, SpaceBlocks retains only matched units containing both
+requested levels and interprets `min_replicates` as the minimum number of complete pairs.
+
+When Cancer and Normal are separate samples matched from the same patients, add a
+shared patient identifier to `core_samples.tsv` and pair by that identifier:
+
+```tsv
+sample	phenotype	patient
+Cancer_P1	Cancer	P1
+Normal_P1	Normal	P1
+Cancer_P2	Cancer	P2
+Normal_P2	Normal	P2
+Cancer_P3	Cancer	P3
+Normal_P3	Normal	P3
+```
+
+```yaml
+- name: matched_phenotype_by_celltype
+  aggregation: by_celltype
+  group_by: [phenotype]
+  paired_by: patient
+  contrasts:
+    - name: cancer_vs_normal
+      numerator: {phenotype: Cancer}
+      denominator: {phenotype: Normal}
+```
+
+Thus, use `paired_by: sample` for multiple region levels from the same spatial sample,
+`paired_by: patient` for separate samples matched by patient, and omit `paired_by` for
+independent Cancer-versus-Normal samples.
+
+When `lrt.enabled` is true, SpaceBlocks runs an omnibus test across every non-excluded
+level with sufficient replication. With `paired_by`, the full model includes the paired
+unit term and the reduced model removes only the comparison group. DEGpatterns is run
+when the LRT finds enough significant genes; no order or trend direction is imposed by
+the configuration.
+
 ## 6. Advanced: Reproducibility and reusability
 
 Even though we have ensured the highest reproducibility standards when creating SpaceBlocks, some steps are never 100% reproducible between systems (e.g. UMAP calculation, Leiden clustering).

diff --git a/config/config.yaml b/config/config.yaml
@@ -110,15 +110,58 @@ analysis:
   # pseudobulk
   run_pseudobulk_de: true
   run_leiden_analysis: true
-  analysis_levels:                  # by_niche_region also available; off by default (too many groups)
-    - by_region
-    - by_celltype_region
   min_cells_per_pseudobulk: 10
   min_replicates: 3
   min_counts_per_pseudobulk: 1000
   padj_threshold: 0.05
   lfc_threshold: 0.25
-  # regions to consider in pseudobulk (omitted names get removed from pseudobulk DE)
+  # Each named analysis defines one aggregation and one DE model. `cell_type` is
+  # resolved to the annotation selected by annotation_types. Multiple group_by
+  # columns form one combined comparison group; interactions are not fitted.
+  pseudobulk:
+    analyses:
+      - name: regions_all_cells
+        aggregation: by_region
+        group_by:
+          - region_annotation
+        paired_by: sample
+        exclude_levels:
+          region_annotation:
+            - Unlabeled
+            - Bubble
+        contrasts:
+          - name: tumor_vs_healthy
+            numerator:
+              region_annotation: Tumor area
+            denominator:
+              region_annotation: Healthy area
+        # Omnibus test across every sufficiently replicated, non-excluded region.
+        lrt:
+          enabled: true
+
+      - name: regions_by_celltype
+        aggregation: by_celltype_region
+        group_by:
+          - region_annotation
+        paired_by: sample
+        exclude_levels:
+          region_annotation:
+            - Unlabeled
+            - Bubble
+          cell_type:
+            - Unannotated
+        contrasts:
+          - name: tumor_vs_healthy
+            numerator:
+              region_annotation: Tumor area
+            denominator:
+              region_annotation: Healthy area
+        lrt:
+          enabled: true
+
+  # Preferred region order for QC-sweep, subcluster, and region-only pseudobulk split
+  # heatmaps, first-to-last (typically left-to-right). This does not filter pseudobulk
+  # levels; use each analysis' exclude_levels for filtering.
   region_levels:
     - Healthy area
     - Healthy over tumor
@@ -128,7 +171,7 @@ analysis:
     - In situ area
     - Tumor area
     - Invasive area
-  # color palette: keys must match region_levels above
+  # Palette used when region_annotation is the single comparison variable.
   region_colors:
     Healthy area: '#FDE725FF'
     Healthy over tumor: '#9FDA3AFF'

diff --git a/docs/design.md b/docs/design.md
@@ -91,7 +91,7 @@ For what every file is and why it is useful, see [Outputs](outputs.md); the tree
 │   ├── annotation/, adata_{sample}_annotated.h5ad   annotate_cells
 │   └── neighbourhood_analysis/{annot_type}/      neighbourhood_analysis
 ├── integrated_samples/                           concatenated / harmony / sketched / samples_report.pdf
-├── pseudobulk/{annot_type}/{analysis_level}/     aggregated/ + de_results/
+├── pseudobulk/{annot_type}/{analysis_name}/      aggregated/ + de_results/
 ├── Subcompartments/{name}/{Harmony,NoHarmony}/   subcluster
 ├── spatial_niches/                               tsv/, plots/
 ├── gene_exploration/                             expression_ranges.tsv, {entry}/…

diff --git a/docs/getting-started.md b/docs/getting-started.md
@@ -115,7 +115,14 @@ snakemake run_postprocessing --sdm conda  # annotate → integrate → DE → ne
 
 This command will produce the integrated object, a per-sample report `sample_report.pdf` to share with collaborators, pseudobulk DE (*optional*, but recommended), neighbourhood/co-occurrence and optional subcompartment re-clustering (for cell subtyping exploration).
 
-## 7. Explore genes and signatures
+## 7. Differential expression
+
+SpaceBlocks performs differential expression on raw counts summed into pseudobulk observations. Comparison variables come from the `core_samples.tsv` sample sheet and may define a single contrast (for example, treated versus untreated), combined categorical groups (for example, phenotype and treatment), or a multilevel variable (for example, cancer stage).
+Importantly, combined groups are compared directly, since SpaceBlocks does not currently fit statistical interaction terms. An optional omnibus LRT can test all retained levels of a multilevel variable, and the R package DEGreport provides summarization of gene expression trends across levels.
+
+Independent samples should omit `paired_by`. When comparison levels are matched within the same patient or spatial sample, `paired_by` identifies that repeated unit. For each explicit pairwise contrast, only complete pairs are tested. See the [pseudobulk experimental-design configuration](configuration.md#pseudobulk-experimental-designs) for independent, combined-group, multilevel, and paired examples.
+
+## 8. Explore genes and signatures
 
 Once the integrated file is available, you may explore any gene or gene set of interest (generated in the results, or completely external).
 
@@ -125,18 +132,18 @@ List the genes / gene sets in the `gene_exploration.queries` TSV; integrated and
 snakemake explore_genes --sdm conda       # AUCell scoring + spatial/dotplot composites
 ```
 
-## 8. Reproducibility
+## 9. Reproducibility
 
-### 8.1 Recommendations at a glance
+### 9.1 Recommendations at a glance
 
 - **Colours** are consistent across every plot — set them once in `config.yaml`
   ([color scale configuration](configuration.md#4-color-scale-customization)).
 - **Reproducibility**: to close the UMAP/Leiden non-determinism gap across machines, share
   pre-computed clusters/annotations or a pre-built contract h5ad (`mode: decoupled`) — see
   [reproducibility configuration](configuration.md#6-advanced-reproducibility-and-reusability).
-- Make sure to back up the intermediate files (see section 7.2) and to [save a locked copy your environment versions](environments.md).
+- Make sure to back up the intermediate files (see section 9.2) and to [save a locked copy of your environment versions](environments.md).
 
-### 8.2 Files to ensure reproducibility
+### 9.2 Files to ensure reproducibility
 
 To reproduce your full run, make sure to back up:
 1. The design files (`config/core_samples.tsv`; `visiumhd_samples.csv`).
@@ -145,7 +152,7 @@ To reproduce your full run, make sure to back up:
 4. The annotated metadata, generated during the postprocessing. Optionally, if using manual annotation, the TSV mapping cluster-to-cell type equivalences.
 5. The environment versions (recommended for publication).
 
-## 9. Example runs
+## 10. Example runs
 
 You may check our [public data end-to-end example runs](demos.md) to see how a full run looks.
 

diff --git a/docs/outputs.md b/docs/outputs.md
@@ -23,7 +23,7 @@ The tag after each folder marks which CoreBlock step produces it — **(pre)** p
 ├── spatial_niches/                                  (pre)   optional
 ├── cluster_annotations_template.tsv                 (pre)   fill this in, then annotate
 ├── integrated_samples/                              (post)  concatenated + integrated objects, report
-├── pseudobulk/<annot_type>/<analysis_level>/        (post)  aggregated/ + de_results/
+├── pseudobulk/<annot_type>/<analysis_name>/         (post)  aggregated/ + de_results/
 ├── Subcompartments/<name>/                          (post)  optional
 └── gene_exploration/                                (expl)
 ```
@@ -50,8 +50,8 @@ The tag after each folder marks which CoreBlock step produces it — **(pre)** p
 | `annotate_cells` | `annotation/` | Per-sample annotation and marker plots |
 | `integrate_samples` | `integrated_samples/harmony_integrated.h5ad` | All samples concatenated and batch-integrated |
 | `integrate_samples` | `integrated_samples/sketched.h5ad` | Geometric-sketch subsample for fast plotting |
-| `pseudobulk_aggregate` | `pseudobulk/<annot_type>/<analysis_level>/aggregated/` | Pseudobulk count matrices and their PCA QC |
-| `pseudobulk_de` | `pseudobulk/<annot_type>/<analysis_level>/de_results/` | Differential-expression tables and plots |
+| `pseudobulk_aggregate` | `pseudobulk/<annot_type>/<analysis_name>/aggregated/` | Pseudobulk count matrices, model metadata, manifest, and PCA QC for one named analysis |
+| `pseudobulk_de` | `pseudobulk/<annot_type>/<analysis_name>/de_results/` | Explicit pairwise differential-expression tables/plots and optional `LRT/` results |
 | `neighbourhood_analysis` | `neighbourhood_analysis/<annot_type>/` | Cell-type co-occurrence and enrichment per sample |
 | `subcluster` | `Subcompartments/<name>/` | Optional re-clustering of a chosen compartment |
 | `sample_report` | `integrated_samples/samples_report.pdf` | One page per sample with UMAP, spatial, dotplot and barplot |