diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 0b5b9ae..cdb3578 100644 --- 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 index a279db1..f5988a1 100644 --- 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 index c383732..d5a97c8 100644 --- 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 index c68222d..949da68 100644 --- 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 index 7789f0f..c4119c5 100644 --- 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 index 07ab397..6360f3e 100644 --- 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/// (post) aggregated/ + de_results/ +├── pseudobulk/// (post) aggregated/ + de_results/ ├── Subcompartments// (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///aggregated/` | Pseudobulk count matrices and their PCA QC | -| `pseudobulk_de` | `pseudobulk///de_results/` | Differential-expression tables and plots | +| `pseudobulk_aggregate` | `pseudobulk///aggregated/` | Pseudobulk count matrices, model metadata, manifest, and PCA QC for one named analysis | +| `pseudobulk_de` | `pseudobulk///de_results/` | Explicit pairwise differential-expression tables/plots and optional `LRT/` results | | `neighbourhood_analysis` | `neighbourhood_analysis//` | Cell-type co-occurrence and enrichment per sample | | `subcluster` | `Subcompartments//` | Optional re-clustering of a chosen compartment | | `sample_report` | `integrated_samples/samples_report.pdf` | One page per sample with UMAP, spatial, dotplot and barplot | diff --git a/docs/rules.md b/docs/rules.md index db07ff1..dcefff6 100644 --- a/docs/rules.md +++ b/docs/rules.md @@ -163,19 +163,20 @@ Concatenates all annotated samples, runs Harmony batch correction on the specifi Writes `concatenated.h5ad`, `harmony_integrated.h5ad`, and `sketched.h5ad`. This is the input to pseudobulk, subclustering (CoreBlock 2), and gene exploration (CoreBlock 3). #### `pseudobulk_aggregate` -Builds pseudobulk count matrices via [decoupler](https://decoupler.readthedocs.io/en/latest/) per `region`, `cell type` and (optionally) `niche` from the integrated object. +Builds raw-count pseudobulk matrices via [decoupler](https://decoupler.readthedocs.io/en/latest/) from the integrated object. Named analyses in `analysis.pseudobulk.analyses` choose one of four aggregations: all cells per sample, per cell type, per region, or per cell type and region. Sample metadata from `core_samples.tsv` are preserved for downstream models, and configured levels are excluded consistently before aggregation and testing. Produces diagnostic QC and PCA plots coloured by sample and by each extra-annotation column. #### `pseudobulk_de` *(optional; `analysis.run_pseudobulk_de: true`)* -Identifies deregulated genes between conditions and annotated regions via differential expression with DESeq2 (R package) on the pseudobulk matrices. +Identifies deregulated genes with DESeq2 on the pseudobulk matrices. The comparison may use region annotations or one or more sample-level variables such as phenotype and treatment. -Differential expression runs in three ways: -- pairwise Wald contrasts between level pairs (region A vs region B), -- each level against the rest, -- and, if there are more than 2 levels, a likelihood-ratio test (LRT) with DEGpatterns for gene-group clustering. +Differential expression runs only for tests requested in the named analysis: -Produces TSV tables, volcano plots, and metadata-annotated heatmaps between conditions. +- explicit pairwise Wald contrasts, with optional categorical covariates; +- paired Wald contrasts using `paired_by: sample`, retaining complete pairs; +- and an optional omnibus likelihood-ratio test (LRT) across sufficiently replicated, non-excluded levels, with DEGpatterns gene-group clustering when enough genes are significant. + +Multiple `group_by` columns create a combined categorical group for comparisons such as drug versus vehicle within one phenotype. Produces TSV tables, volcano plots, and metadata-annotated heatmaps between configured groups. Note that SpaceBlocks does not currently support fitting statistical interactions. #### `neighbourhood_analysis` Runs squidpy neighbourhood enrichment per sample and per `region`, plus per-region cell-type co-occurrence (a full co-occurrence figure per region and compact pairwise heatmaps). diff --git a/workflow/Snakefile b/workflow/Snakefile index 43e9600..f3df267 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -9,7 +9,9 @@ directly below; helper functions live in rules/common.smk. import glob import os +import re import sys +from itertools import combinations import pandas as pd from snakemake.utils import validate, min_version @@ -208,7 +210,29 @@ _ANNOT_TYPE_TO_COL = { DEFAULT_ANNOT_COL = _ANNOT_TYPE_TO_COL.get(PRIMARY_ANNOT_TYPE, "cell_type_tsv") ANNOT_TYPES = list(config.get("annotation_types") or [PRIMARY_ANNOT_TYPE]) -ANALYSIS_LEVELS = ANALYSIS.get("analysis_levels", ["by_region", "by_celltype_region"]) +_PSEUDOBULK_AGGREGATIONS = { + "all_cells", + "by_celltype", + "by_region", + "by_celltype_region", +} + + +PSEUDOBULK = ANALYSIS.get("pseudobulk", {}) or {} +_configured_pseudobulk_analyses = PSEUDOBULK.get("analyses") +if _configured_pseudobulk_analyses: + PSEUDOBULK_ANALYSES = [ + _validate_pseudobulk_analysis(spec) for spec in _configured_pseudobulk_analyses + ] +else: + PSEUDOBULK_ANALYSES = [ + _validate_pseudobulk_analysis(spec) for spec in _legacy_pseudobulk_analyses() + ] + +PSEUDOBULK_ANALYSIS_BY_NAME = {spec["name"]: spec for spec in PSEUDOBULK_ANALYSES} +if len(PSEUDOBULK_ANALYSIS_BY_NAME) != len(PSEUDOBULK_ANALYSES): + sys.exit("[config error] analysis.pseudobulk analysis names must be unique.") +PSEUDOBULK_ANALYSIS_NAMES = list(PSEUDOBULK_ANALYSIS_BY_NAME) # Pseudobulk DESeq2 DE is OPTIONAL and on by default RUN_DE = bool(ANALYSIS.get("run_pseudobulk_de", False)) @@ -357,7 +381,7 @@ rule run_postprocessing: expand( rules.pseudobulk_de.output.results_dir, annot_type=ANNOT_TYPES, - analysis_level=ANALYSIS_LEVELS, + analysis_name=PSEUDOBULK_ANALYSIS_NAMES, ) if RUN_DE else [] diff --git a/workflow/envs/dev.yaml b/workflow/envs/dev.yaml index 25df53f..620f59b 100644 --- a/workflow/envs/dev.yaml +++ b/workflow/envs/dev.yaml @@ -10,6 +10,7 @@ dependencies: # Formatting and linting tools - pre_commit - bioconda::snakefmt + - conda-forge::go-shfmt # Documentation tools - mkdocs diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index f501c5d..64e7c0e 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -134,6 +134,209 @@ def _extra_annot_color_triples(): return cols, vals, colors +def _validate_pseudobulk_analysis(spec): + """Validate pseudobulk settings that JSON Schema cannot check across fields.""" + if not isinstance(spec, dict): + sys.exit( + "[config error] each analysis.pseudobulk.analyses entry must be a mapping." + ) + + name = str(spec.get("name", "")) + if not re.fullmatch(r"[A-Za-z0-9_-]+", name): + sys.exit( + "[config error] pseudobulk analysis names may contain only letters, " + f"numbers, '_' and '-'; got '{name}'." + ) + + aggregation = str(spec.get("aggregation", "")) + if aggregation not in _PSEUDOBULK_AGGREGATIONS: + sys.exit( + f"[config error] pseudobulk analysis '{name}' has unsupported aggregation " + f"'{aggregation}'. Choose one of {sorted(_PSEUDOBULK_AGGREGATIONS)}." + ) + + group_by = list(spec.get("group_by") or []) + if not group_by or any(not isinstance(c, str) or not c for c in group_by): + sys.exit( + f"[config error] pseudobulk analysis '{name}' needs at least one " + "non-empty group_by column." + ) + if len(group_by) != len(set(group_by)): + sys.exit( + f"[config error] pseudobulk analysis '{name}' repeats a group_by column." + ) + + paired_by = spec.get("paired_by", "") or "" + if not isinstance(paired_by, str): + sys.exit( + f"[config error] pseudobulk analysis '{name}'.paired_by must be a string." + ) + + covariates = list(spec.get("covariates") or []) + if any(not isinstance(c, str) or not c for c in covariates): + sys.exit( + f"[config error] pseudobulk analysis '{name}' has an invalid covariate." + ) + if len(covariates) != len(set(covariates)): + sys.exit(f"[config error] pseudobulk analysis '{name}' repeats a covariate.") + overlap = sorted(set(group_by) & set(covariates)) + if overlap: + sys.exit( + f"[config error] pseudobulk analysis '{name}' lists {overlap} as both " + "group_by and covariates." + ) + if paired_by and paired_by in group_by: + sys.exit( + f"[config error] pseudobulk analysis '{name}' lists '{paired_by}' as both " + "paired_by and group_by." + ) + if paired_by and paired_by in covariates: + sys.exit( + f"[config error] pseudobulk analysis '{name}' lists '{paired_by}' as both " + "paired_by and a covariate." + ) + + exclude_levels = spec.get("exclude_levels") or {} + if not isinstance(exclude_levels, dict): + sys.exit( + f"[config error] pseudobulk analysis '{name}'.exclude_levels must be a mapping." + ) + + contrast_names = set() + for contrast in list(spec.get("contrasts") or []): + if not isinstance(contrast, dict): + sys.exit( + f"[config error] contrasts in pseudobulk analysis '{name}' must be mappings." + ) + contrast_name = str(contrast.get("name", "")) + if not re.fullmatch(r"[A-Za-z0-9_-]+", contrast_name): + sys.exit( + f"[config error] contrast names in pseudobulk analysis '{name}' may " + "contain only letters, numbers, '_' and '-'." + ) + if contrast_name in contrast_names: + sys.exit( + f"[config error] duplicate contrast '{contrast_name}' in '{name}'." + ) + contrast_names.add(contrast_name) + for side in ("numerator", "denominator"): + values = contrast.get(side) + if not isinstance(values, dict) or set(values) != set(group_by): + sys.exit( + f"[config error] contrast '{contrast_name}' {side} must specify " + f"exactly the group_by columns {group_by}." + ) + excluded_values = { + column: {str(value) for value in exclude_levels.get(column, [])} + for column in group_by + } + conflicts = [ + f"{column}={value}" + for column, value in values.items() + if str(value) in excluded_values[column] + ] + if conflicts: + sys.exit( + f"[config error] contrast '{contrast_name}' {side} in '{name}' " + f"uses excluded level(s): {conflicts}." + ) + if contrast["numerator"] == contrast["denominator"]: + sys.exit( + f"[config error] contrast '{contrast_name}' in '{name}' has identical " + "numerator and denominator groups." + ) + + return { + **spec, + "name": name, + "aggregation": aggregation, + "group_by": group_by, + "paired_by": paired_by, + "covariates": covariates, + "exclude_levels": dict(exclude_levels), + "contrasts": list(spec.get("contrasts") or []), + "lrt": dict(spec.get("lrt") or {}), + } + + +def _legacy_pseudobulk_analyses(): + """Translate the former region-only config for existing demo/user configs. + + New configurations should use analysis.pseudobulk.analyses. The unimplemented + by_niche_region value is deliberately rejected rather than carried forward. + """ + levels = list(ANALYSIS.get("analysis_levels", ["by_region", "by_celltype_region"])) + unsupported = sorted(set(levels) - {"by_region", "by_celltype_region"}) + if unsupported: + sys.exit( + "[config error] unsupported legacy pseudobulk analysis_levels: " + f"{unsupported}. 'by_niche_region' has been removed because it was never implemented." + ) + contrasts = [ + { + "name": re.sub(r"[^A-Za-z0-9_-]+", "_", f"{a}_vs_{b}").strip("_"), + "numerator": {"region_annotation": str(a)}, + "denominator": {"region_annotation": str(b)}, + } + for a, b in combinations(REGION_LEVELS, 2) + ] + return [ + { + "name": level, + "aggregation": level, + "group_by": ["region_annotation"], + "paired_by": "sample", + "covariates": [], + "exclude_levels": {"region_annotation": ["Unlabeled", "Bubble"]}, + "contrasts": contrasts, + "lrt": {"enabled": len(REGION_LEVELS) > 2}, + } + for level in levels + ] + + +def pseudobulk_analysis(name): + """Return one validated named pseudobulk analysis from the Snakefile globals.""" + return PSEUDOBULK_ANALYSIS_BY_NAME[str(name)] + + +def pseudobulk_group_label(spec, values): + """Stable display/factor label for one configured group combination.""" + columns = spec["group_by"] + if len(columns) == 1: + return str(values[columns[0]]) + return " | ".join(f"{column}={values[column]}" for column in columns) + + +def pseudobulk_contrast_param(analysis_name, field): + """Flatten named contrast mappings for safe Snakemake Python-to-R transfer.""" + spec = pseudobulk_analysis(analysis_name) + if field == "names": + return [str(c["name"]) for c in spec["contrasts"]] + if field not in {"numerator", "denominator"}: + raise ValueError(f"Unknown pseudobulk contrast field: {field}") + return [pseudobulk_group_label(spec, c[field]) for c in spec["contrasts"]] + + +def pseudobulk_group_palette(analysis_name): + """Palette for a one-column comparison; combined groups use R fallbacks.""" + spec = pseudobulk_analysis(analysis_name) + if len(spec["group_by"]) != 1: + return {} + column = spec["group_by"][0] + if column == "region_annotation": + return ANALYSIS.get("region_colors", {}) or {} + return SAMPLE_COLORS.get(column, {}) if isinstance(SAMPLE_COLORS, dict) else {} + + +def pseudobulk_condition_order(analysis_name): + """Preferred heatmap order for a region-only comparison.""" + spec = pseudobulk_analysis(analysis_name) + if spec["group_by"] == ["region_annotation"]: + return [str(level) for level in REGION_LEVELS] + return [] + + def core_sample_meta(sample): """Experimental-design metadata for `sample` (design columns only, excludes the `sample` key). Returns {} when no core sheet is configured.""" @@ -403,13 +606,13 @@ def get_all_targets(wildcards): targets += expand( rules.pseudobulk_aggregate.output.agg_dir, annot_type=ANNOT_TYPES, - analysis_level=ANALYSIS_LEVELS, + analysis_name=PSEUDOBULK_ANALYSIS_NAMES, ) if RUN_DE: targets += expand( rules.pseudobulk_de.output.results_dir, annot_type=ANNOT_TYPES, - analysis_level=ANALYSIS_LEVELS, + analysis_name=PSEUDOBULK_ANALYSIS_NAMES, ) targets.append(rules.integrate_samples.output.concatenated) targets.append(rules.integrate_samples.output.harmony) diff --git a/workflow/rules/pseudobulk_aggregate.smk b/workflow/rules/pseudobulk_aggregate.smk index 312f2d0..5bd1287 100644 --- a/workflow/rules/pseudobulk_aggregate.smk +++ b/workflow/rules/pseudobulk_aggregate.smk @@ -1,28 +1,28 @@ rule pseudobulk_aggregate: """ - Aggregate single cells into pseudobulk count matrices using decoupler. + Aggregate single cells into pseudobulk count matrices using a named analysis. Output structure: - {annot_type}/{analysis_level}/ + {annot_type}/{analysis_name}/aggregated/ ├── matrices/ count TSVs ├── metadata/ metadata TSVs - ├── plots/ sample-level QC (pseudobulk bar, PCA) + ├── plots/ sample-level QC └── manifest.tsv """ input: adata=rules.integrate_samples.output.concatenated, output: agg_dir=directory( - f"{OUTDIR_PP}/pseudobulk/{{annot_type}}/{{analysis_level}}/aggregated" + f"{OUTDIR_PP}/pseudobulk/{{annot_type}}/{{analysis_name}}/aggregated" ), log: - out=f"{LOGDIR}/pseudobulk_aggregate/{{annot_type}}_{{analysis_level}}.out", - err=f"{LOGDIR}/pseudobulk_aggregate/{{annot_type}}_{{analysis_level}}.err", + out=f"{LOGDIR}/pseudobulk_aggregate/{{annot_type}}_{{analysis_name}}.out", + err=f"{LOGDIR}/pseudobulk_aggregate/{{annot_type}}_{{analysis_name}}.err", benchmark: - f"{LOGDIR}/benchmarks/pseudobulk_aggregate/{{annot_type}}_{{analysis_level}}.tsv" + f"{LOGDIR}/benchmarks/pseudobulk_aggregate/{{annot_type}}_{{analysis_name}}.tsv" wildcard_constraints: - annot_type="tsv_annotation|ingest_annotation", - analysis_level="by_region|by_celltype_region|by_niche_region", + annot_type="tsv_annotation|external_annotation|ingest_annotation|refined_annotation", + analysis_name="[A-Za-z0-9_-]+", conda: "../envs/pseudobulk_aggregate.yaml" threads: get_resource("pseudobulk_aggregate", "threads") @@ -31,7 +31,7 @@ rule pseudobulk_aggregate: runtime=get_resource("pseudobulk_aggregate", "runtime"), params: annot_type=lambda wc: wc.annot_type, - analysis_level=lambda wc: wc.analysis_level, + analysis_spec=lambda wc: pseudobulk_analysis(wc.analysis_name), min_cells_per_pseudobulk=ANALYSIS.get("min_cells_per_pseudobulk", 10), min_counts_per_pseudobulk=ANALYSIS.get("min_counts_per_pseudobulk", 1000), annotation_colors=config.get("annotation_colors", {}), diff --git a/workflow/rules/pseudobulk_de.smk b/workflow/rules/pseudobulk_de.smk index 15c5536..0a2ca20 100644 --- a/workflow/rules/pseudobulk_de.smk +++ b/workflow/rules/pseudobulk_de.smk @@ -1,22 +1,19 @@ rule pseudobulk_de: - """ - DE with R DESeq2 (Wald + LRT), EnhancedVolcano, ComplexHeatmap, DEGpatterns. - Region levels and colors from config control ordering and visualization. - """ + """Explicit pairwise Wald contrasts and optional omnibus LRT with R/DESeq2.""" input: agg_dir=rules.pseudobulk_aggregate.output.agg_dir, output: results_dir=directory( - f"{OUTDIR_PP}/pseudobulk/{{annot_type}}/{{analysis_level}}/de_results" + f"{OUTDIR_PP}/pseudobulk/{{annot_type}}/{{analysis_name}}/de_results" ), log: - out=f"{LOGDIR}/pseudobulk_de/{{annot_type}}_{{analysis_level}}.out", - err=f"{LOGDIR}/pseudobulk_de/{{annot_type}}_{{analysis_level}}.err", + out=f"{LOGDIR}/pseudobulk_de/{{annot_type}}_{{analysis_name}}.out", + err=f"{LOGDIR}/pseudobulk_de/{{annot_type}}_{{analysis_name}}.err", benchmark: - f"{LOGDIR}/benchmarks/pseudobulk_de/{{annot_type}}_{{analysis_level}}.tsv" + f"{LOGDIR}/benchmarks/pseudobulk_de/{{annot_type}}_{{analysis_name}}.tsv" wildcard_constraints: - annot_type="tsv_annotation|refined_annotation|ingest_annotation", - analysis_level="by_region|by_celltype_region|by_niche_region", + annot_type="tsv_annotation|external_annotation|ingest_annotation|refined_annotation", + analysis_name="[A-Za-z0-9_-]+", conda: "../envs/pseudobulk_de.yaml" threads: get_resource("pseudobulk_de", "threads") @@ -24,18 +21,32 @@ rule pseudobulk_de: mem_mb=mem_mb_attempt("pseudobulk_de"), runtime=get_resource("pseudobulk_de", "runtime"), params: - annot_type=lambda wc: wc.annot_type, - analysis_level=lambda wc: wc.analysis_level, + analysis_name=lambda wc: wc.analysis_name, min_replicates=ANALYSIS.get("min_replicates", 3), de_n_genes=ANALYSIS.get("de_n_genes", 10), padj_threshold=ANALYSIS.get("padj_threshold", 0.05), lfc_threshold=ANALYSIS.get("lfc_threshold", 0.5), - region_levels=REGION_LEVELS, - # Pass colors as two parallel lists for safe Python→R conversion - region_color_names=list(ANALYSIS.get("region_colors", {}).keys()), - region_color_values=list(ANALYSIS.get("region_colors", {}).values()), - # Design columns to annotate heatmaps with (chosen via config design.columns); - # palettes flattened into parallel lists (grey fallback applied in R). + dpi=ANALYSIS.get("plot_dpi", 300), + group_by_columns=lambda wc: pseudobulk_analysis(wc.analysis_name)["group_by"], + paired_by=lambda wc: pseudobulk_analysis(wc.analysis_name)["paired_by"], + covariates=lambda wc: pseudobulk_analysis(wc.analysis_name)["covariates"], + contrast_names=lambda wc: pseudobulk_contrast_param(wc.analysis_name, "names"), + contrast_numerators=lambda wc: pseudobulk_contrast_param( + wc.analysis_name, "numerator" + ), + contrast_denominators=lambda wc: pseudobulk_contrast_param( + wc.analysis_name, "denominator" + ), + lrt_enabled=lambda wc: bool( + pseudobulk_analysis(wc.analysis_name)["lrt"].get("enabled", False) + ), + group_color_names=lambda wc: list( + pseudobulk_group_palette(wc.analysis_name).keys() + ), + group_color_values=lambda wc: list( + pseudobulk_group_palette(wc.analysis_name).values() + ), + condition_level_order=lambda wc: pseudobulk_condition_order(wc.analysis_name), extra_annot_columns=EXTRA_ANNOT_COLUMNS, extra_anno_col_names=EXTRA_ANNO_COLS, extra_anno_values=EXTRA_ANNO_VALS, diff --git a/workflow/schemas/config.schema.yaml b/workflow/schemas/config.schema.yaml index 6e5bbba..c4f95dc 100644 --- a/workflow/schemas/config.schema.yaml +++ b/workflow/schemas/config.schema.yaml @@ -210,7 +210,8 @@ properties: type: ["object", "null"] description: >- core_samples.tsv columns to surface downstream as UMAP colourings, pseudobulk - QC panels, and DE heatmap annotation tracks. Annotation-only (no DESeq2 covariate). + QC panels, and DE heatmap annotation tracks. This controls plotting only; + sample-sheet columns may independently be selected as pseudobulk covariates. properties: columns: type: array @@ -302,6 +303,32 @@ properties: type: boolean default: false description: Gate the optional R/DESeq2 pseudobulk_de step. + min_cells_per_pseudobulk: + type: integer + minimum: 1 + default: 10 + description: Minimum cells contributing to one retained pseudobulk observation. + min_counts_per_pseudobulk: + type: integer + minimum: 1 + default: 1000 + description: Minimum summed raw counts in one retained pseudobulk observation. + min_replicates: + type: integer + minimum: 2 + default: 3 + description: >- + Minimum independent samples per unpaired level, or complete pairs when + paired_by is configured, required to run a contrast or include an LRT level. + padj_threshold: + type: number + exclusiveMinimum: 0 + maximum: 1 + default: 0.05 + lfc_threshold: + type: number + minimum: 0 + default: 0.5 run_leiden_analysis: type: boolean default: true @@ -313,17 +340,110 @@ properties: type: array items: type: string - description: Pseudobulk analysis levels (e.g. by_region, by_celltype_region). + enum: ["by_region", "by_celltype_region"] + description: >- + Deprecated region-only pseudobulk configuration retained for older + configs. New configs should use analysis.pseudobulk.analyses. region_levels: type: array items: type: string - description: Ordered region names (must match region_annotation values in the data). + description: >- + Optional preferred first-to-last display order for region names in QC-sweep, + subcluster, and region-only pseudobulk split heatmaps. This does not filter + pseudobulk observations; use pseudobulk analyses' exclude_levels instead. region_colors: type: object additionalProperties: type: string description: region name -> hex colour. + pseudobulk: + type: object + additionalProperties: false + description: >- + Named pseudobulk aggregations and DE models. Each analysis defines how + cells are pooled, which sample metadata form the comparison group, + exclusions, explicit Wald contrasts, and an optional omnibus LRT. + required: [analyses] + properties: + analyses: + type: array + minItems: 1 + items: + type: object + additionalProperties: false + required: [name, aggregation, group_by] + properties: + name: + type: string + pattern: "^[A-Za-z0-9_-]+$" + description: Stable analysis identifier used in output paths. + aggregation: + type: string + enum: ["all_cells", "by_celltype", "by_region", "by_celltype_region"] + group_by: + type: array + minItems: 1 + uniqueItems: true + items: + type: string + description: >- + One or more obs/core_samples columns whose values define the + comparison groups. Multiple columns form a combined group. + paired_by: + type: string + default: "" + description: >- + Optional matched-unit column, normally sample for comparisons + between regions from the same sample. Empty means unpaired. + covariates: + type: array + uniqueItems: true + default: [] + items: + type: string + description: Categorical sample metadata adjusted for before group_by. + exclude_levels: + type: object + default: {} + additionalProperties: + type: array + uniqueItems: true + items: + type: string + description: >- + Column -> levels removed before aggregation. The same exclusion + applies to explicit Wald contrasts and the optional LRT. + contrasts: + type: array + default: [] + items: + type: object + additionalProperties: false + required: [name, numerator, denominator] + properties: + name: + type: string + pattern: "^[A-Za-z0-9_-]+$" + numerator: + type: object + additionalProperties: + type: string + denominator: + type: object + additionalProperties: + type: string + lrt: + type: object + additionalProperties: false + default: {} + properties: + enabled: + type: boolean + default: false + description: >- + Run an omnibus LRT across all sufficiently replicated levels + remaining after exclude_levels. plot_dpi: type: integer minimum: 1 diff --git a/workflow/scripts/pseudobulk_aggregate.py b/workflow/scripts/pseudobulk_aggregate.py index 9a61ef4..6ae41c7 100644 --- a/workflow/scripts/pseudobulk_aggregate.py +++ b/workflow/scripts/pseudobulk_aggregate.py @@ -11,6 +11,7 @@ import logging import os +import re import sys import traceback from pathlib import Path @@ -48,6 +49,7 @@ "tsv_annotation": "cell_type_tsv", "ingest_annotation": "cell_type_ingest", "external_annotation": "cell_type_external", + "refined_annotation": "cell_type_refined", } @@ -216,134 +218,249 @@ def _plot_pseudobulk_qc(pdata, condition_col, sample_col, prefix, plots_dir, plt.close() +def _safe_name(value): + safe = re.sub(r"[^A-Za-z0-9_.-]+", "_", str(value)).strip("_") + return safe or "group" + + +def _group_label(row, columns): + if len(columns) == 1: + return str(row[columns[0]]) + return " | ".join(f"{column}={row[column]}" for column in columns) + + +class DesignMetadataError(ValueError): + """Requested model metadata is not constant within a count unit.""" + + +def _attach_design_metadata(pdata, source, sample_col, groups_col, columns): + """Attach metadata that is constant within every pseudobulk count unit. + + decoupler retains grouping columns, but preservation of unrelated obs columns + is version-dependent. Reconstructing them explicitly makes the DE contract + stable and detects invalid designs where a requested value varies inside one + pseudobulk observation. + """ + columns = list(dict.fromkeys(c for c in columns if c not in {sample_col, groups_col})) + if not columns: + return pdata + + keys = [sample_col, groups_col] + frame = source.obs[keys + columns].copy() + for column in keys + columns: + frame[column] = frame[column].astype(str) + + grouped = frame.groupby(keys, observed=True, dropna=False) + for column in columns: + variable = grouped[column].nunique(dropna=False) + if (variable > 1).any(): + examples = [" / ".join(map(str, idx if isinstance(idx, tuple) else [idx])) + for idx in variable[variable > 1].index[:5]] + raise DesignMetadataError( + f"Design column '{column}' varies within pseudobulk unit(s): {examples}. " + "Use an aggregation that keeps this variable separate." + ) + + lookup = grouped[columns].first() + pdata_keys = pdata.obs[keys].astype(str) + lookup_index = pd.MultiIndex.from_frame(pdata_keys) + for column in columns: + pdata.obs[column] = lookup[column].reindex(lookup_index).to_numpy() + return pdata + + +def _aggregate(source, groups_col, design_columns): + pdata = dc.get_pseudobulk( + source, + sample_col="sample", + groups_col=groups_col, + layer="raw_counts", + mode="sum", + min_cells=MIN_CELLS, + min_counts=MIN_COUNTS, + ) + return _attach_design_metadata( + pdata, source, "sample", groups_col, design_columns + ) + + # ── Parameters ─────────────────────────────────────────────────────────────── -annot_type = snakemake.params.annot_type -analysis_level = snakemake.params.analysis_level -MIN_CELLS = int(snakemake.params.min_cells_per_pseudobulk) -MIN_COUNTS = int(snakemake.params.min_counts_per_pseudobulk) +annot_type = str(snakemake.params.annot_type) +analysis_spec = dict(snakemake.params.analysis_spec) +analysis_name = str(analysis_spec["name"]) +aggregation = str(analysis_spec["aggregation"]) +group_by = list(analysis_spec["group_by"]) +paired_by = str(analysis_spec.get("paired_by", "") or "") +covariates = list(analysis_spec.get("covariates") or []) +exclude_levels = dict(analysis_spec.get("exclude_levels") or {}) +MIN_CELLS = int(snakemake.params.min_cells_per_pseudobulk) +MIN_COUNTS = int(snakemake.params.min_counts_per_pseudobulk) ANNOTATION_COLORS = snakemake.params.annotation_colors -REGION_COLORS = snakemake.params.region_colors -DPI = int(getattr(snakemake.params, "dpi", 300)) +REGION_COLORS = snakemake.params.region_colors +DPI = int(getattr(snakemake.params, "dpi", 300)) EXTRA_ANNOT_COLUMNS = list(getattr(snakemake.params, "extra_annot_columns", []) or []) -SAMPLE_COLORS = getattr(snakemake.params, "sample_colors", {}) or {} -agg_dir = str(snakemake.output.agg_dir) +SAMPLE_COLORS = getattr(snakemake.params, "sample_colors", {}) or {} +agg_dir = str(snakemake.output.agg_dir) try: log.info("=" * 70) - log.info("Pseudobulk aggregation: %s / %s", annot_type, analysis_level) + log.info("Pseudobulk aggregation: %s / %s", annot_type, analysis_name) + log.info(" aggregation=%s, group_by=%s", aggregation, group_by) log.info(" min_cells=%d, min_counts=%d", MIN_CELLS, MIN_COUNTS) log.info("=" * 70) matrices_dir = os.path.join(agg_dir, "matrices") metadata_dir = os.path.join(agg_dir, "metadata") - plots_dir = os.path.join(agg_dir, "plots") - for d in [agg_dir, matrices_dir, metadata_dir, plots_dir]: - os.makedirs(d, exist_ok=True) + plots_dir = os.path.join(agg_dir, "plots") + for directory in [agg_dir, matrices_dir, metadata_dir, plots_dir]: + os.makedirs(directory, exist_ok=True) adata = sc.read_h5ad(str(snakemake.input.adata)) log.info("Loaded: %d cells, %d genes", adata.n_obs, adata.n_vars) - cell_type_col = ANNOT_COL_MAP.get(annot_type) - if cell_type_col and cell_type_col not in adata.obs.columns: - log.warning("Column '%s' not found. Skipping.", cell_type_col) - Path(os.path.join(agg_dir, f"SKIPPED_no_{cell_type_col}.txt")).write_text( - f"Column {cell_type_col} not found.\n" - ) - sys.exit(0) + # sample_batch is created by integrate_samples from the authoritative core + # sample sheet. Prefer it over a potentially renamed contract obs['sample']. + if "sample_batch" in adata.obs.columns: + adata.obs["sample"] = adata.obs["sample_batch"].astype(str).to_numpy() + elif "sample" not in adata.obs.columns: + raise ValueError("No canonical sample column found.") - sample_col = None - for c in ["sample", "sample_batch"]: - if c in adata.obs.columns: - sample_col = c - break - if sample_col is None: - raise ValueError("No sample column found.") - - region_col = "region_annotation" - if region_col not in adata.obs.columns: - log.warning("No region_annotation. Skipping.") + cell_type_col = ANNOT_COL_MAP.get(annot_type) + referenced_columns = set(group_by + covariates + list(exclude_levels)) + if paired_by: + referenced_columns.add(paired_by) + needs_cell_type = ( + aggregation in {"by_celltype", "by_celltype_region"} + or "cell_type" in referenced_columns + ) + if needs_cell_type: + if not cell_type_col or cell_type_col not in adata.obs.columns: + log.warning("Column '%s' not found. Skipping.", cell_type_col) + Path(os.path.join(agg_dir, f"SKIPPED_no_{cell_type_col}.txt")).write_text( + f"Column {cell_type_col} not found.\n" + ) + sys.exit(0) + adata.obs["cell_type"] = adata.obs[cell_type_col].astype(str).to_numpy() + + needs_region = aggregation in {"by_region", "by_celltype_region"} + if needs_region and "region_annotation" not in adata.obs.columns: + log.warning("No region_annotation for region-based aggregation. Skipping.") Path(os.path.join(agg_dir, "SKIPPED_no_regions.txt")).write_text("") sys.exit(0) - for col in [region_col, sample_col]: - adata = _clean_nan(adata, col) - if cell_type_col: - adata = _clean_nan(adata, cell_type_col) - - valid_regions = [r for r in adata.obs[region_col].unique() - if r not in ("Unlabeled", "Bubble")] - adata = adata[adata.obs[region_col].isin(valid_regions)].copy() - log.info("Regions: %s (%d cells)", valid_regions, adata.n_obs) + required_columns = list(dict.fromkeys( + ["sample"] + group_by + covariates + ([paired_by] if paired_by else []) + )) + missing = [column for column in required_columns if column not in adata.obs.columns] + if missing: + raise ValueError( + f"Pseudobulk analysis '{analysis_name}' requires missing obs column(s): {missing}. " + "Add them to core_samples.tsv or correct group_by/covariates/paired_by." + ) - if len(valid_regions) < 2: - Path(os.path.join(agg_dir, "SKIPPED_fewer_than_2_regions.txt")).write_text("") + for column in required_columns: + adata = _clean_nan(adata, column) + + # Exclusions are shared by aggregation, explicit Wald contrasts, and LRT. + for column, excluded in exclude_levels.items(): + if column not in adata.obs.columns: + raise ValueError( + f"exclude_levels references missing column '{column}' in analysis '{analysis_name}'." + ) + excluded = {str(value) for value in (excluded or [])} + if excluded: + values = adata.obs[column].astype(str) + keep = ~values.isin(excluded) + log.info(" Excluding %d cells from %s levels: %s", + int((~keep).sum()), column, sorted(excluded)) + adata = adata[keep].copy() + + if adata.n_obs == 0: + Path(os.path.join(agg_dir, "SKIPPED_no_cells_after_exclusion.txt")).write_text("") sys.exit(0) - if "raw_counts" not in adata.layers: raise ValueError("raw_counts layer missing.") + group_palette = {} + if len(group_by) == 1: + group_column = group_by[0] + if group_column == "region_annotation": + group_palette = REGION_COLORS if isinstance(REGION_COLORS, dict) else {} + elif isinstance(SAMPLE_COLORS, dict): + group_palette = SAMPLE_COLORS.get(group_column, {}) or {} + + # Plot-only annotations are optional. Model columns were validated above and + # remain mandatory, while absent plotting annotations should not abort DE. + present_extra_annotations = [ + column for column in EXTRA_ANNOT_COLUMNS if column in adata.obs.columns + ] + design_columns = list(dict.fromkeys( + group_by + covariates + ([paired_by] if paired_by else []) + present_extra_annotations + )) manifest_rows = [] + used_prefixes = set() + + if aggregation in {"all_cells", "by_celltype"}: + adata.obs["__all__"] = "all" + groups_col = "__all__" + else: + groups_col = "region_annotation" - if analysis_level == "by_region": - pdata = dc.get_pseudobulk( - adata, sample_col=sample_col, groups_col=region_col, - layer="raw_counts", mode="sum", - min_cells=MIN_CELLS, min_counts=MIN_COUNTS, + if aggregation in {"all_cells", "by_region"}: + subgroups = [("pooled", "all_cells", adata)] + else: + subgroups = [] + for cell_type in sorted(adata.obs["cell_type"].astype(str).unique()): + prefix = _safe_name(cell_type) + if prefix in used_prefixes: + raise ValueError( + f"Cell-type names collide after filename sanitization: '{cell_type}' -> '{prefix}'." + ) + used_prefixes.add(prefix) + subset = adata[adata.obs["cell_type"].astype(str) == cell_type].copy() + subgroups.append((prefix, cell_type, subset)) + + for prefix, grouping, subset in subgroups: + if subset.n_obs < MIN_CELLS: + log.warning(" '%s': too few cells (%d). Skipping.", grouping, subset.n_obs) + continue + try: + pdata = _aggregate(subset, groups_col, design_columns) + except DesignMetadataError: + raise + except Exception as exc: + log.warning(" '%s': aggregation failed: %s", grouping, exc) + continue + if pdata.n_obs == 0: + log.warning(" '%s': no pseudobulk observations passed QC. Skipping.", grouping) + continue + + pdata.obs[".comparison_group"] = [ + _group_label(row, group_by) for _, row in pdata.obs.iterrows() + ] + _save_pseudobulk(pdata, matrices_dir, metadata_dir, prefix) + _plot_pseudobulk_qc( + pdata, ".comparison_group", "sample", prefix, plots_dir, + group_palette, ANNOTATION_COLORS, + EXTRA_ANNOT_COLUMNS, SAMPLE_COLORS, ) - _save_pseudobulk(pdata, matrices_dir, metadata_dir, "pooled") - _plot_pseudobulk_qc(pdata, region_col, sample_col, "pooled", plots_dir, - REGION_COLORS, ANNOTATION_COLORS, - EXTRA_ANNOT_COLUMNS, SAMPLE_COLORS) manifest_rows.append({ - "prefix": "pooled", "grouping": "all_cells", - "n_samples": pdata.n_obs, "n_genes": pdata.n_vars, - "condition_col": region_col, "sample_col": sample_col, + "prefix": prefix, + "grouping": grouping, + "n_samples": pdata.n_obs, + "n_genes": pdata.n_vars, + "condition_col": ".comparison_group", + "sample_col": "sample", + "analysis_name": analysis_name, + "aggregation": aggregation, }) - elif analysis_level == "by_celltype_region": - adata = adata[adata.obs[cell_type_col] != "Unannotated"].copy() - cell_types = sorted(adata.obs[cell_type_col].unique()) - - for ct in cell_types: - safe_ct = ct.replace("/", "_").replace(" ", "_") - ct_adata = adata[adata.obs[cell_type_col] == ct].copy() - - if ct_adata.n_obs < MIN_CELLS: - log.warning(" '%s': too few cells (%d). Skipping.", ct, ct_adata.n_obs) - continue - - try: - pdata = dc.get_pseudobulk( - ct_adata, sample_col=sample_col, groups_col=region_col, - layer="raw_counts", mode="sum", - min_cells=MIN_CELLS, min_counts=MIN_COUNTS, - ) - except Exception as e: - log.warning(" '%s': aggregation failed: %s", ct, e) - continue - - if pdata.n_obs < 4: - log.warning(" '%s': too few pseudobulk samples (%d). Skipping.", ct, pdata.n_obs) - continue - - _save_pseudobulk(pdata, matrices_dir, metadata_dir, safe_ct) - _plot_pseudobulk_qc(pdata, region_col, sample_col, safe_ct, plots_dir, - REGION_COLORS, ANNOTATION_COLORS, - EXTRA_ANNOT_COLUMNS, SAMPLE_COLORS) - manifest_rows.append({ - "prefix": safe_ct, "grouping": ct, - "n_samples": pdata.n_obs, "n_genes": pdata.n_vars, - "condition_col": region_col, "sample_col": sample_col, - }) - - elif analysis_level == "by_niche_region": - log.warning("by_niche_region not yet implemented.") - Path(os.path.join(agg_dir, "SKIPPED_niche_not_implemented.txt")).write_text("") - if manifest_rows: pd.DataFrame(manifest_rows).to_csv( os.path.join(agg_dir, "manifest.tsv"), sep="\t", index=False ) + else: + Path(os.path.join(agg_dir, "SKIPPED_no_pseudobulks.txt")).write_text("") log.info("Aggregation complete.") diff --git a/workflow/scripts/pseudobulk_de.R b/workflow/scripts/pseudobulk_de.R index 6316fcd..7805e08 100644 --- a/workflow/scripts/pseudobulk_de.R +++ b/workflow/scripts/pseudobulk_de.R @@ -1,13 +1,10 @@ ################################################################################ -## pseudobulk_de.R – DESeq2 Wald + LRT, ComplexHeatmap, DEGpatterns, -## DEGpatterns +## pseudobulk_de.R – explicit DESeq2 Wald contrasts and optional omnibus LRT ################################################################################ -# ── Redirect logs FIRST — before anything else ────────────────────────────── -# This ensures any crash during param parsing is captured in the log files. +# Redirect logs before parsing parameters so early failures are captured. log_out <- tryCatch(snakemake@log[["out"]], error = function(e) "") log_err <- tryCatch(snakemake@log[["err"]], error = function(e) "") - if (!is.null(log_out) && nzchar(log_out)) { dir.create(dirname(log_out), recursive = TRUE, showWarnings = FALSE) log_out_con <- file(log_out, open = "wt") @@ -20,422 +17,514 @@ if (!is.null(log_err) && nzchar(log_err)) { } message("=== pseudobulk_de.R starting ===") -message("Loading libraries …") - suppressPackageStartupMessages({ library(DESeq2) library(ComplexHeatmap) library(circlize) - library(RColorBrewer) - library(tidyverse) + library(dplyr) + library(tibble) + library(ggplot2) library(ggrepel) - library(DEGreport) }) -message("Libraries loaded OK") - -# ── Snakemake params (with diagnostics) ────────────────────────────────────── -message("Parsing snakemake params …") - -agg_dir <- snakemake@input[["agg_dir"]] -results_dir <- snakemake@output[["results_dir"]] -message(" agg_dir: ", agg_dir) -message(" results_dir: ", results_dir) - +agg_dir <- snakemake@input[["agg_dir"]] +results_dir <- snakemake@output[["results_dir"]] +analysis_name <- as.character(snakemake@params[["analysis_name"]]) min_replicates <- as.integer(snakemake@params[["min_replicates"]]) -de_n_genes <- as.integer(snakemake@params[["de_n_genes"]]) -padj_thr <- as.numeric(snakemake@params[["padj_threshold"]]) -lfc_thr <- as.numeric(snakemake@params[["lfc_threshold"]]) -# Plot resolution (config param); defaults to 300 if the rule doesn't pass it. -.dpi_param <- snakemake@params[["dpi"]] -res_dpi <- if (is.null(.dpi_param)) 300L else as.integer(.dpi_param) -px_scale <- res_dpi / 150 # keep physical size, scale pixels -message(" min_reps: ", min_replicates, ", de_n_genes: ", de_n_genes, - ", padj: ", padj_thr, ", lfc: ", lfc_thr) - -# Region levels -message(" Parsing region_levels …") -region_levels <- tryCatch({ - rl <- snakemake@params[["region_levels"]] - if (is.null(rl)) character(0) else as.character(rl) -}, error = function(e) { - message(" WARNING: Failed to parse region_levels: ", conditionMessage(e)) - character(0) -}) -message(" region_levels: [", paste(region_levels, collapse = ", "), "]") - -# Region colors — reconstruct named vector from parallel lists -message(" Parsing region_colors …") -region_colors <- tryCatch({ - rc_names <- snakemake@params[["region_color_names"]] - rc_values <- snakemake@params[["region_color_values"]] - message(" rc_names type: ", class(rc_names), " length: ", length(rc_names)) - message(" rc_values type: ", class(rc_values), " length: ", length(rc_values)) - if (length(rc_names) > 0 && length(rc_names) == length(rc_values)) { - setNames(as.character(rc_values), as.character(rc_names)) - } else { - character(0) - } -}, error = function(e) { - message(" WARNING: Failed to parse region_colors: ", conditionMessage(e)) - character(0) -}) -message(" region_colors: ", paste(names(region_colors), region_colors, sep = "=", collapse = ", ")) +de_n_genes <- as.integer(snakemake@params[["de_n_genes"]]) +padj_thr <- as.numeric(snakemake@params[["padj_threshold"]]) +lfc_thr <- as.numeric(snakemake@params[["lfc_threshold"]]) +res_dpi <- as.integer(snakemake@params[["dpi"]]) +group_by_columns <- as.character(snakemake@params[["group_by_columns"]]) +paired_by <- as.character(snakemake@params[["paired_by"]]) +covariates <- as.character(snakemake@params[["covariates"]]) +contrast_names <- as.character(snakemake@params[["contrast_names"]]) +contrast_numerators <- as.character(snakemake@params[["contrast_numerators"]]) +contrast_denominators <- as.character(snakemake@params[["contrast_denominators"]]) +lrt_enabled <- isTRUE(as.logical(snakemake@params[["lrt_enabled"]])) +paired_by <- paired_by[nzchar(paired_by)] +covariates <- covariates[nzchar(covariates)] +condition_level_order <- tryCatch( + as.character(snakemake@params[["condition_level_order"]]), + error = function(e) character(0) +) +condition_level_order <- unique(condition_level_order[nzchar(condition_level_order)]) +px_scale <- res_dpi / 150 + +if (!(length(contrast_names) == length(contrast_numerators) && + length(contrast_names) == length(contrast_denominators))) { + stop("Contrast parameter lists have inconsistent lengths.") +} -# ── Design annotation columns + palettes (parallel lists → named vectors) ───── -extra_annot_columns <- tryCatch(as.character(snakemake@params[["extra_annot_columns"]]), - error = function(e) character(0)) +# Reconstruct configured palettes from parallel lists. +group_colors <- tryCatch({ + nms <- as.character(snakemake@params[["group_color_names"]]) + vals <- as.character(snakemake@params[["group_color_values"]]) + if (length(nms) > 0 && length(nms) == length(vals)) setNames(vals, nms) + else character(0) +}, error = function(e) character(0)) + +extra_annot_columns <- tryCatch( + as.character(snakemake@params[["extra_annot_columns"]]), + error = function(e) character(0) +) extra_annot_columns <- extra_annot_columns[nzchar(extra_annot_columns)] extra_palettes <- tryCatch({ - dc <- as.character(snakemake@params[["extra_anno_col_names"]]) - dv <- as.character(snakemake@params[["extra_anno_values"]]) - dh <- as.character(snakemake@params[["extra_anno_colors"]]) - pals <- list() - if (length(dc) > 0 && length(dc) == length(dv) && length(dv) == length(dh)) { - for (u in unique(dc)) { - idx <- which(dc == u) - pals[[u]] <- setNames(dh[idx], dv[idx]) - } - } - pals -}, error = function(e) { - message(" WARNING: Failed to parse design palettes: ", conditionMessage(e)) - list() -}) -message(" extra_annot_columns: [", paste(extra_annot_columns, collapse = ", "), "]") - -# Build a ComplexHeatmap top annotation: Region (existing behaviour) plus one -# track per design column present in `meta_df`, coloured from extra_palettes with -# a grey (#cccccc) fallback for values absent from the config palette. Design -# columns are annotation-only — they are NOT added to the DESeq2 model. -build_top_annotation <- function(meta_df, condition_col, region_cols) { - anno_args <- list(Region = meta_df[[condition_col]]) - anno_cols <- list() - avail <- region_cols[names(region_cols) %in% unique(as.character(meta_df[[condition_col]]))] - if (length(avail) > 0) anno_cols[["Region"]] <- avail - for (dc in extra_annot_columns) { - if (dc %in% colnames(meta_df)) { - vals <- as.character(meta_df[[dc]]) - anno_args[[dc]] <- vals - lv <- unique(vals) - full <- setNames(rep("#cccccc", length(lv)), lv) - pal <- extra_palettes[[dc]] - if (!is.null(pal)) { - common <- intersect(names(pal), lv) - if (length(common) > 0) full[common] <- pal[common] - } - anno_cols[[dc]] <- full + cols <- as.character(snakemake@params[["extra_anno_col_names"]]) + vals <- as.character(snakemake@params[["extra_anno_values"]]) + hex <- as.character(snakemake@params[["extra_anno_colors"]]) + palettes <- list() + if (length(cols) > 0 && length(cols) == length(vals) && length(vals) == length(hex)) { + for (column in unique(cols)) { + idx <- which(cols == column) + palettes[[column]] <- setNames(hex[idx], vals[idx]) } } - if (length(anno_cols) > 0) anno_args[["col"]] <- anno_cols - anno_args[["show_legend"]] <- TRUE - do.call(HeatmapAnnotation, anno_args) + palettes +}, error = function(e) list()) + +message(" analysis: ", analysis_name) +message(" group_by: [", paste(group_by_columns, collapse = ", "), "]") +message(" paired_by: ", ifelse(length(paired_by), paired_by, "")) +message(" covariates: [", paste(covariates, collapse = ", "), "]") +message(" contrasts: [", paste(contrast_names, collapse = ", "), "]") +message(" LRT: ", lrt_enabled) +dir.create(results_dir, recursive = TRUE, showWarnings = FALSE) + + +safe_name <- function(value) { + safe <- gsub("[^A-Za-z0-9_.-]+", "_", as.character(value)) + safe <- gsub("^_+|_+$", "", safe) + ifelse(nzchar(safe), safe, "contrast") } -message("Params parsed OK") -dir.create(results_dir, recursive = TRUE, showWarnings = FALSE) +heatmap_condition_factor <- function(values) { + observed <- unique(as.character(values)) + preferred <- condition_level_order[condition_level_order %in% observed] + factor(as.character(values), levels = unique(c(preferred, observed))) +} -# ── ComplexHeatmap helper ──────────────────────────────────────────────────── -draw_complex_heatmaps <- function(mat_scaled, meta, condition_col, top_genes, - contrast_name, out_dir, region_lvls, region_cols) { - # Reorder columns by region level - if (length(region_lvls) > 0) { - col_order_factor <- factor(meta[[condition_col]], levels = region_lvls, ordered = TRUE) - } else { - col_order_factor <- factor(meta[[condition_col]]) - } - ord <- order(col_order_factor) - mat_ordered <- mat_scaled[, ord, drop = FALSE] - meta_ordered <- meta[ord, , drop = FALSE] - split_factor <- factor(meta_ordered[[condition_col]], - levels = if (length(region_lvls) > 0) region_lvls else unique(meta_ordered[[condition_col]]), - ordered = TRUE) - - # Drop genes that became all-NaN after scaling (zero variance) — they - # collapse the layout and break row clustering. - finite_rows <- apply(mat_ordered, 1, function(z) all(is.finite(z))) - if (any(!finite_rows)) { - message(" Dropping ", sum(!finite_rows), - " zero-variance gene(s) from heatmap") - mat_ordered <- mat_ordered[finite_rows, , drop = FALSE] +formula_for <- function(condition_col, paired = character(0), covars = character(0), + include_condition = TRUE) { + terms <- c(covars, paired) + if (include_condition) terms <- c(terms, condition_col) + terms <- unique(terms[nzchar(terms)]) + if (length(terms) == 0) return(~1) + as.formula(paste("~", paste(sprintf("`%s`", terms), collapse = " + "))) +} + + +assert_full_rank <- function(meta, design_formula, label) { + matrix <- model.matrix(design_formula, data = meta) + rank <- qr(matrix)$rank + if (rank < ncol(matrix)) { + stop( + label, " model matrix is not full rank (rank ", rank, "/", ncol(matrix), + "). Check for confounded covariates or pairing variables." + ) } - if (nrow(mat_ordered) < 2) { - message(" <2 finite genes left — skipping heatmap") - return(invisible(NULL)) + invisible(matrix) +} + + +build_top_annotation <- function(meta, condition_col) { + annotation_args <- list(Group = as.character(meta[[condition_col]])) + annotation_colors <- list() + present_groups <- unique(as.character(meta[[condition_col]])) + available <- group_colors[names(group_colors) %in% present_groups] + if (length(available) > 0) annotation_colors[["Group"]] <- available + + for (column in extra_annot_columns) { + if (column %in% colnames(meta)) { + values <- as.character(meta[[column]]) + annotation_args[[column]] <- values + levels <- unique(values) + palette <- setNames(rep("#cccccc", length(levels)), levels) + configured <- extra_palettes[[column]] + if (!is.null(configured)) { + common <- intersect(names(configured), levels) + palette[common] <- configured[common] + } + annotation_colors[[column]] <- palette + } } + if (length(annotation_colors) > 0) annotation_args[["col"]] <- annotation_colors + annotation_args[["show_legend"]] <- TRUE + do.call(HeatmapAnnotation, annotation_args) +} - # Pin the heatmap BODY so each gene row has a fixed physical height. - n_row <- nrow(mat_ordered); n_col <- ncol(mat_ordered) - row_mm <- 5 - body_mm <- max(n_row * row_mm, 40) - body_h <- unit(body_mm, "mm") - dev_h_px <- round((body_mm / 25.4 + 2.6) * res_dpi) # +2.6in for title/anno/legend - # Region + design-column annotation tracks (grey fallback for missing values) - col_anno <- build_top_annotation(meta_ordered, condition_col, region_cols) +draw_complex_heatmaps <- function(mat_scaled, meta, condition_col, contrast_name, out_dir) { + finite_rows <- apply(mat_scaled, 1, function(values) all(is.finite(values))) + mat_scaled <- mat_scaled[finite_rows, , drop = FALSE] + if (nrow(mat_scaled) < 2) return(invisible(NULL)) - col_fun <- colorRamp2(c(-2, 0, 2), c("blue", "white", "red")) + condition <- heatmap_condition_factor(meta[[condition_col]]) + order_idx <- order(condition) + mat_ordered <- mat_scaled[, order_idx, drop = FALSE] + meta_ordered <- meta[order_idx, , drop = FALSE] + split_factor <- droplevels(condition[order_idx]) + annotation <- build_top_annotation(meta_ordered, condition_col) + color_function <- colorRamp2(c(-2, 0, 2), c("blue", "white", "red")) + + n_row <- nrow(mat_ordered) + n_col <- ncol(mat_ordered) + body_mm <- max(n_row * 5, 40) + body_height <- grid::unit(body_mm, "mm") + device_height <- round((body_mm / 25.4 + 2.6) * res_dpi) - # 1. Unsplit heatmap tryCatch({ - ht_unsplit <- Heatmap( - mat_ordered, - name = "Z-score", - col = col_fun, - top_annotation = col_anno, - height = body_h, - cluster_columns = TRUE, - cluster_rows = TRUE, - show_row_names = TRUE, - show_column_names = TRUE, - row_names_gp = gpar(fontsize = 7), - column_names_gp = gpar(fontsize = 7), - column_title = paste0("Top DE genes — ", contrast_name), - heatmap_legend_param = list(title = "Z-score") + heatmap <- Heatmap( + mat_ordered, name = "Z-score", col = color_function, + top_annotation = annotation, height = body_height, + cluster_columns = TRUE, cluster_rows = TRUE, + show_row_names = TRUE, show_column_names = TRUE, + row_names_gp = grid::gpar(fontsize = 7), column_names_gp = grid::gpar(fontsize = 7), + column_title = paste0("Top DE genes — ", contrast_name) + ) + png( + file.path(out_dir, paste0("heatmap_unsplit_", safe_name(contrast_name), ".png")), + width = round(max(800, n_col * 50) * px_scale), + height = device_height, res = res_dpi ) - png(file.path(out_dir, paste0("heatmap_unsplit_", contrast_name, ".png")), - width = round(max(800, n_col * 50) * px_scale), - height = dev_h_px, res = res_dpi) - draw(ht_unsplit, merge_legend = TRUE) + draw(heatmap, merge_legend = TRUE) dev.off() }, error = function(e) message(" Unsplit heatmap failed: ", conditionMessage(e))) - # 2. Split heatmap by region level tryCatch({ - ht_split <- Heatmap( - mat_ordered, - name = "Z-score", - col = col_fun, - top_annotation = col_anno, - column_split = split_factor, - height = body_h, - cluster_columns = FALSE, - cluster_column_slices = TRUE, - cluster_rows = TRUE, - show_row_names = TRUE, - show_column_names = TRUE, - row_names_gp = gpar(fontsize = 7), - column_names_gp = gpar(fontsize = 7), - column_title = paste0("Top DE genes — ", contrast_name, " (split)"), - heatmap_legend_param = list(title = "Z-score") + heatmap <- Heatmap( + mat_ordered, name = "Z-score", col = color_function, + top_annotation = annotation, column_split = split_factor, + height = body_height, cluster_columns = FALSE, + cluster_column_slices = FALSE, cluster_rows = TRUE, + show_row_names = TRUE, show_column_names = TRUE, + row_names_gp = grid::gpar(fontsize = 7), column_names_gp = grid::gpar(fontsize = 7), + column_title = paste0("Top DE genes — ", contrast_name, " (split)") + ) + png( + file.path(out_dir, paste0("heatmap_split_", safe_name(contrast_name), ".png")), + width = round(max(900, n_col * 55) * px_scale), + height = device_height, res = res_dpi ) - png(file.path(out_dir, paste0("heatmap_split_", contrast_name, ".png")), - width = round(max(900, n_col * 55) * px_scale), - height = dev_h_px, res = res_dpi) - draw(ht_split, merge_legend = TRUE) + draw(heatmap, merge_legend = TRUE) dev.off() }, error = function(e) message(" Split heatmap failed: ", conditionMessage(e))) } -# ── Helper: run DE on one pseudobulk subgroup ──────────────────────────────── -run_de_for_prefix <- function(prefix, condition_col, sample_col, out_base) { +write_volcano <- function(results_df, contrast_name, prefix, out_dir) { + volcano <- results_df %>% + filter(!is.na(padj) & !is.na(log2FoldChange)) %>% + mutate( + sig_group = case_when( + padj < padj_thr & log2FoldChange > lfc_thr ~ "Up", + padj < padj_thr & log2FoldChange < -lfc_thr ~ "Down", + TRUE ~ "NS" + ), + neg_log10p = -log10(pmax(padj, 1e-300)) + ) + labels <- bind_rows( + volcano %>% filter(sig_group == "Up") %>% arrange(desc(log2FoldChange)) %>% head(10), + volcano %>% filter(sig_group == "Down") %>% arrange(log2FoldChange) %>% head(10) + ) + plot <- ggplot(volcano, aes(x = log2FoldChange, y = neg_log10p, color = sig_group)) + + geom_point(size = 0.8, alpha = 0.6) + + scale_color_manual(values = c("NS" = "grey70", "Up" = "#e41a1c", "Down" = "#377eb8"), name = "") + + geom_hline(yintercept = -log10(padj_thr), linetype = "dashed", color = "grey40", linewidth = 0.3) + + geom_vline(xintercept = c(-lfc_thr, lfc_thr), linetype = "dashed", color = "grey40", linewidth = 0.3) + + ggrepel::geom_text_repel( + data = labels, aes(label = gene), size = 2.5, max.overlaps = 20, + color = "black", segment.size = 0.2, segment.color = "grey50" + ) + + labs(title = contrast_name, subtitle = paste0(prefix, " — Wald test"), + x = "log2 Fold Change", y = "-log10(padj)") + + theme_bw() + ggsave( + file.path(out_dir, paste0("volcano_", safe_name(contrast_name), ".png")), + plot = plot, width = 9, height = 7, dpi = res_dpi + ) +} - counts_file <- file.path(agg_dir, "matrices", paste0("counts_", prefix, ".tsv")) - metadata_file <- file.path(agg_dir, "metadata", paste0("metadata_", prefix, ".tsv")) - if (!file.exists(counts_file) || !file.exists(metadata_file)) { - message(" Skipping '", prefix, "': files not found.") - return(NULL) +write_de_heatmap <- function(results_df, vsd, meta, condition_col, contrast_name, out_dir) { + significant <- results_df %>% filter(padj < padj_thr) + top_up <- significant %>% + filter(log2FoldChange > 0) %>% arrange(desc(log2FoldChange)) %>% + head(de_n_genes) %>% pull(gene) + top_down <- significant %>% + filter(log2FoldChange < 0) %>% arrange(log2FoldChange) %>% + head(de_n_genes) %>% pull(gene) + genes <- unique(c(top_up, top_down)) + genes <- genes[genes %in% rownames(assay(vsd))] + if (length(genes) >= 2) { + scaled <- t(scale(t(assay(vsd)[genes, , drop = FALSE]))) + draw_complex_heatmaps(scaled, meta, condition_col, contrast_name, out_dir) } +} - dir.create(out_base, recursive = TRUE, showWarnings = FALSE) - counts <- read.delim(counts_file, row.names = 1, check.names = FALSE) - meta <- read.delim(metadata_file, row.names = 1, check.names = FALSE) +fit_wald <- function(counts, meta, condition_col, paired = character(0)) { + design_formula <- formula_for(condition_col, paired, covariates, TRUE) + assert_full_rank(meta, design_formula, "Wald") + counts_t <- t(counts) + counts_t <- counts_t[, rownames(meta), drop = FALSE] + dds <- DESeqDataSetFromMatrix(countData = counts_t, colData = meta, design = design_formula) + dds <- DESeq(dds, test = "Wald") + vsd <- tryCatch( + vst(dds, blind = FALSE), + error = function(e) { + message(" vst failed, using normTransform: ", conditionMessage(e)) + normTransform(dds) + } + ) + list(dds = dds, vsd = vsd) +} - counts <- round(counts) - counts <- as.matrix(counts) - storage.mode(counts) <- "integer" - common <- intersect(rownames(counts), rownames(meta)) - if (length(common) == 0) { - message(" '", prefix, "': no common samples. Skipping.") - return(NULL) - } - counts <- counts[common, , drop = FALSE] - meta <- meta[common, , drop = FALSE] - - # Filter to configured region levels if provided - if (length(region_levels) > 0) { - keep <- meta[[condition_col]] %in% region_levels - if (sum(!keep) > 0) { - dropped <- unique(meta[[condition_col]][!keep]) - message(" Filtering out region levels not in config: ", paste(dropped, collapse = ", ")) - } - counts <- counts[keep, , drop = FALSE] - meta <- meta[keep, , drop = FALSE] - meta[[condition_col]] <- factor(meta[[condition_col]], levels = region_levels) - } else { - meta[[condition_col]] <- factor(meta[[condition_col]]) +run_lrt <- function(counts, meta, condition_col, sample_col, prefix, out_base) { + if (!lrt_enabled) return(invisible(NULL)) + message(" Explicit omnibus LRT …") + + unit_col <- if (length(paired_by)) paired_by else sample_col + replicate_counts <- meta %>% + mutate(.unit = as.character(.data[[unit_col]])) %>% + group_by(.data[[condition_col]]) %>% + summarise(n = n_distinct(.unit), .groups = "drop") + valid <- as.character(replicate_counts[[condition_col]][replicate_counts$n >= min_replicates]) + lrt_meta <- meta[as.character(meta[[condition_col]]) %in% valid, , drop = FALSE] + lrt_counts <- counts[rownames(lrt_meta), , drop = FALSE] + lrt_meta[[condition_col]] <- droplevels(factor(lrt_meta[[condition_col]])) + + if (length(paired_by)) { + informative <- lrt_meta %>% + mutate(.row = rownames(lrt_meta)) %>% + group_by(.data[[paired_by]]) %>% + filter(n_distinct(.data[[condition_col]]) >= 2) %>% + pull(.row) + lrt_meta <- lrt_meta[informative, , drop = FALSE] + lrt_counts <- lrt_counts[informative, , drop = FALSE] + lrt_meta[[paired_by]] <- droplevels(factor(lrt_meta[[paired_by]])) } - meta[[condition_col]] <- droplevels(meta[[condition_col]]) - conditions <- levels(meta[[condition_col]]) - - cond_counts <- table(meta[[condition_col]]) - valid_conds <- names(cond_counts[cond_counts >= min_replicates]) - if (length(valid_conds) < 2) { - message(" '", prefix, "': < 2 conditions with >= ", min_replicates, " replicates. Skipping.") - writeLines("Insufficient replicates.", file.path(out_base, "SKIPPED_insufficient.txt")) - return(NULL) + lrt_meta[[condition_col]] <- droplevels(factor(lrt_meta[[condition_col]])) + if (nlevels(lrt_meta[[condition_col]]) < 2) { + message(" LRT skipped: fewer than two non-excluded levels have sufficient replicates.") + lrt_dir <- file.path(out_base, "LRT") + dir.create(lrt_dir, recursive = TRUE, showWarnings = FALSE) + writeLines("Insufficient replicated levels for LRT.", file.path(lrt_dir, "SKIPPED_insufficient.txt")) + return(invisible(NULL)) } - keep_samples <- meta[[condition_col]] %in% valid_conds - counts <- counts[keep_samples, , drop = FALSE] - meta <- meta[keep_samples, , drop = FALSE] - meta[[condition_col]] <- droplevels(meta[[condition_col]]) - conditions <- levels(meta[[condition_col]]) - gene_sums <- colSums(counts) - counts <- counts[, gene_sums > 0, drop = FALSE] + full_formula <- formula_for(condition_col, paired_by, covariates, TRUE) + reduced_formula <- formula_for(condition_col, paired_by, covariates, FALSE) + assert_full_rank(lrt_meta, full_formula, "LRT full") + assert_full_rank(lrt_meta, reduced_formula, "LRT reduced") + counts_t <- t(lrt_counts) + counts_t <- counts_t[, rownames(lrt_meta), drop = FALSE] + dds <- DESeqDataSetFromMatrix(countData = counts_t, colData = lrt_meta, design = full_formula) + dds <- DESeq(dds, test = "LRT", reduced = reduced_formula) + result <- results(dds) + # An LRT p-value tests the group term as a whole. DESeq2 also returns one + # coefficient's fold change, which is not an omnibus effect size and is easy to + # misinterpret when there are more than two levels, so omit it here. + result_df <- as.data.frame(result) %>% + rownames_to_column("gene") %>% + select(gene, baseMean, stat, pvalue, padj) %>% + arrange(padj) + lrt_dir <- file.path(out_base, "LRT") + dir.create(lrt_dir, recursive = TRUE, showWarnings = FALSE) + write.table(result_df, file.path(lrt_dir, "LRT_results.tsv"), + sep = "\t", quote = FALSE, row.names = FALSE) + + significant <- rownames(result)[which(result$padj < padj_thr)] + message(" LRT: ", length(significant), " significant genes") + if (length(significant) < 10) return(invisible(NULL)) - message(" '", prefix, "': ", nrow(counts), " samples, ", ncol(counts), " genes, ", - length(conditions), " conditions: ", paste(conditions, collapse = " -> ")) + tryCatch({ + transformed <- tryCatch(vst(dds, blind = FALSE), error = function(e) normTransform(dds)) + matrix <- assay(transformed)[significant, , drop = FALSE] + clusters <- DEGreport::degPatterns( + matrix, metadata = lrt_meta, time = condition_col, plot = FALSE + ) + if (!is.null(clusters$df)) { + normalized <- clusters$normalized + old_levels <- sort(unique(normalized$cluster)) + remap <- setNames(seq_along(old_levels), as.character(old_levels)) + normalized$cluster <- as.integer(remap[as.character(normalized$cluster)]) + cluster_df <- clusters$df + cluster_df$cluster <- as.integer(remap[as.character(cluster_df$cluster)]) + cluster_df <- cluster_df %>% arrange(cluster, genes) + + plot <- DEGreport::degPlotCluster( + normalized, time = condition_col, color = condition_col, points = TRUE + ) + + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + if (length(group_colors) > 0) { + plot <- plot + scale_color_manual(values = group_colors) + } + ggsave(file.path(lrt_dir, "DEGpatterns_groups.png"), plot = plot, + width = 14, height = 10, dpi = max(res_dpi, 500L)) + + grouped_genes <- split(cluster_df$genes, paste0("group_", cluster_df$cluster)) + max_length <- max(lengths(grouped_genes)) + gene_table <- as.data.frame(lapply(grouped_genes, function(genes) { + c(genes, rep("", max_length - length(genes))) + })) + write.table(gene_table, file.path(lrt_dir, "gene_groups.tsv"), + sep = "\t", quote = FALSE, row.names = FALSE) + } + }, error = function(e) message(" DEGpatterns failed: ", conditionMessage(e))) +} - counts_t <- t(counts) - if (!identical(colnames(counts_t), rownames(meta))) { - meta <- meta[colnames(counts_t), , drop = FALSE] - } +run_de_for_prefix <- function(prefix, condition_col, sample_col, out_base) { + counts_file <- file.path(agg_dir, "matrices", paste0("counts_", prefix, ".tsv")) + metadata_file <- file.path(agg_dir, "metadata", paste0("metadata_", prefix, ".tsv")) + if (!file.exists(counts_file) || !file.exists(metadata_file)) return(NULL) + dir.create(out_base, recursive = TRUE, showWarnings = FALSE) - design_formula <- as.formula(paste("~", condition_col)) - dds <- DESeqDataSetFromMatrix( - countData = counts_t, colData = meta, design = design_formula - ) + counts <- read.delim(counts_file, row.names = 1, check.names = FALSE) + meta <- read.delim(metadata_file, row.names = 1, check.names = FALSE) + counts <- round(as.matrix(counts)) + storage.mode(counts) <- "integer" + common <- intersect(rownames(counts), rownames(meta)) + counts <- counts[common, , drop = FALSE] + meta <- meta[common, , drop = FALSE] + + required <- unique(c(condition_col, sample_col, paired_by, covariates)) + missing <- setdiff(required, colnames(meta)) + if (length(missing) > 0) stop("Missing pseudobulk metadata columns: ", paste(missing, collapse = ", ")) + complete <- complete.cases(meta[, required, drop = FALSE]) + counts <- counts[complete, , drop = FALSE] + meta <- meta[complete, , drop = FALSE] + meta[[condition_col]] <- factor(meta[[condition_col]]) + for (column in unique(c(paired_by, covariates))) meta[[column]] <- factor(meta[[column]]) + + unit_col <- if (length(paired_by)) paired_by else sample_col + unit_groups <- meta[, c(unit_col, condition_col), drop = FALSE] + duplicate_units <- duplicated(unit_groups) | duplicated(unit_groups, fromLast = TRUE) + if (any(duplicate_units)) { + examples <- unique(apply(unit_groups[duplicate_units, , drop = FALSE], 1, paste, collapse = " / ")) + stop( + "More than one pseudobulk observation exists for an independent unit/group: ", + paste(head(examples, 5), collapse = ", "), + ". Choose an aggregation aligned with group_by; repeated observations are not independent." + ) + } + if (!length(paired_by)) { + groups_per_sample <- tapply( + as.character(meta[[condition_col]]), as.character(meta[[sample_col]]), + function(values) length(unique(values)) + ) + if (any(groups_per_sample > 1)) { + stop( + "At least one sample contributes multiple comparison levels. Configure paired_by: ", + sample_col, " for a paired design." + ) + } + } - # ── LRT ──────────────────────────────────────────────────────────────── - message(" LRT …") - dds_lrt <- NULL - tryCatch({ - dds_lrt <- DESeq(dds, test = "LRT", reduced = ~ 1) - res_lrt <- results(dds_lrt) - res_lrt_df <- as.data.frame(res_lrt) %>% - rownames_to_column("gene") %>% - arrange(desc(log2FoldChange)) - # LRT results saved inside DEGpatterns dir (created later) - lrt_save_dir <- file.path(out_base, "DEGpatterns") - dir.create(lrt_save_dir, recursive = TRUE, showWarnings = FALSE) - write.table(res_lrt_df, file.path(lrt_save_dir, "LRT_results.tsv"), - sep = "\t", quote = FALSE, row.names = FALSE) - n_sig <- sum(res_lrt_df$padj < padj_thr, na.rm = TRUE) - message(" LRT: ", n_sig, " significant genes (padj < ", padj_thr, ")") - }, error = function(e) { - message(" LRT failed: ", conditionMessage(e)) - }) - - # ── Pairwise Wald ────────────────────────────────────────────────────── - message(" Wald pairwise tests …") - dds_wald <- DESeq(dds, test = "Wald") + counts <- counts[, colSums(counts) > 0, drop = FALSE] + if (nrow(counts) == 0 || ncol(counts) == 0) { + writeLines("No usable counts.", file.path(out_base, "SKIPPED_no_counts.txt")) + return(NULL) + } - vsd <- tryCatch( - vst(dds_wald, blind = FALSE), + message(" '", prefix, "': ", nrow(counts), " observations, ", ncol(counts), " genes") + tryCatch( + run_lrt(counts, meta, condition_col, sample_col, prefix, out_base), error = function(e) { - message(" vst failed, using normTransform: ", conditionMessage(e)) - normTransform(dds_wald) + message(" LRT failed: ", conditionMessage(e)) + lrt_dir <- file.path(out_base, "LRT") + dir.create(lrt_dir, recursive = TRUE, showWarnings = FALSE) + writeLines(conditionMessage(e), file.path(lrt_dir, "ERROR.txt")) } ) pairwise_dir <- file.path(out_base, "pairwise") dir.create(pairwise_dir, recursive = TRUE, showWarnings = FALSE) summary_rows <- list() - pairs <- combn(conditions, 2, simplify = FALSE) - for (pair in pairs) { - cond_a <- pair[1] - cond_b <- pair[2] - contrast_name <- paste0(cond_a, "_vs_", cond_b) - message(" Contrast: ", contrast_name) + # Unpaired analyses share a single fit across all sufficiently replicated levels. + shared_fit <- NULL + shared_meta <- NULL + if (!length(paired_by) && length(contrast_names) > 0) { + replicate_counts <- meta %>% + mutate(.sample = as.character(.data[[sample_col]])) %>% + group_by(.data[[condition_col]]) %>% + summarise(n = n_distinct(.sample), .groups = "drop") + valid <- as.character(replicate_counts[[condition_col]][replicate_counts$n >= min_replicates]) + shared_meta <- meta[as.character(meta[[condition_col]]) %in% valid, , drop = FALSE] + shared_meta[[condition_col]] <- droplevels(factor(shared_meta[[condition_col]])) + if (nlevels(shared_meta[[condition_col]]) >= 2) { + shared_counts <- counts[rownames(shared_meta), , drop = FALSE] + shared_fit <- fit_wald(shared_counts, shared_meta, condition_col) + } + } + + for (index in seq_along(contrast_names)) { + contrast_name <- contrast_names[index] + numerator <- contrast_numerators[index] + denominator <- contrast_denominators[index] + message(" Contrast: ", contrast_name, " (", numerator, " vs ", denominator, ")") tryCatch({ - res <- results(dds_wald, contrast = c(condition_col, cond_a, cond_b)) - res_df <- as.data.frame(res) %>% - rownames_to_column("gene") %>% - arrange(desc(log2FoldChange)) + if (length(paired_by)) { + pair_meta <- meta[as.character(meta[[condition_col]]) %in% c(numerator, denominator), , drop = FALSE] + pair_meta[[condition_col]] <- droplevels(factor(pair_meta[[condition_col]])) + complete_units <- pair_meta %>% + mutate(.row = rownames(pair_meta)) %>% + group_by(.data[[paired_by]]) %>% + filter(all(c(numerator, denominator) %in% as.character(.data[[condition_col]]))) %>% + pull(.row) + pair_meta <- pair_meta[complete_units, , drop = FALSE] + pair_meta[[condition_col]] <- droplevels(factor(pair_meta[[condition_col]])) + pair_meta[[paired_by]] <- droplevels(factor(pair_meta[[paired_by]])) + n_replicates <- n_distinct(pair_meta[[paired_by]]) + if (n_replicates < min_replicates) { + stop("only ", n_replicates, " complete pairs; need ", min_replicates) + } + pair_counts <- counts[rownames(pair_meta), , drop = FALSE] + fit <- fit_wald(pair_counts, pair_meta, condition_col, paired_by) + fit_meta <- pair_meta + } else { + if (is.null(shared_fit) || + !all(c(numerator, denominator) %in% levels(shared_meta[[condition_col]]))) { + stop("one or both levels have fewer than ", min_replicates, " independent samples") + } + fit <- shared_fit + fit_meta <- shared_meta + } - write.table(res_df, file.path(pairwise_dir, paste0(contrast_name, ".tsv")), + result <- results( + fit$dds, + contrast = c(condition_col, numerator, denominator), + alpha = padj_thr + ) + result_df <- as.data.frame(result) %>% + rownames_to_column("gene") %>% arrange(desc(log2FoldChange)) + output_name <- safe_name(contrast_name) + write.table(result_df, file.path(pairwise_dir, paste0(output_name, ".tsv")), sep = "\t", quote = FALSE, row.names = FALSE) - n_sig <- sum(res_df$padj < padj_thr, na.rm = TRUE) - n_up <- sum(res_df$padj < padj_thr & res_df$log2FoldChange > 0, na.rm = TRUE) - n_down <- sum(res_df$padj < padj_thr & res_df$log2FoldChange < 0, na.rm = TRUE) - message(" ", n_sig, " sig (up:", n_up, ", down:", n_down, ")") - + n_sig <- sum(result_df$padj < padj_thr, na.rm = TRUE) + n_up <- sum(result_df$padj < padj_thr & result_df$log2FoldChange > 0, na.rm = TRUE) + n_down <- sum(result_df$padj < padj_thr & result_df$log2FoldChange < 0, na.rm = TRUE) summary_rows[[contrast_name]] <- tibble( - contrast = contrast_name, n_tested = nrow(res_df), - n_significant = n_sig, n_up = n_up, n_down = n_down + contrast = contrast_name, numerator = numerator, denominator = denominator, + n_tested = nrow(result_df), n_significant = n_sig, n_up = n_up, n_down = n_down, + error = NA_character_ ) - - # Custom volcano plot (grey=NS, red=up, blue=down, top genes labelled) - tryCatch({ - volc_df <- res_df %>% - filter(!is.na(padj) & !is.na(log2FoldChange)) %>% - mutate( - sig_group = case_when( - padj >= padj_thr ~ "NS", - log2FoldChange > lfc_thr ~ "Up", - log2FoldChange < -lfc_thr ~ "Down", - TRUE ~ "NS" - ), - neg_log10p = -log10(pmax(padj, 1e-300)) - ) - - # Top genes to label - top_up <- volc_df %>% filter(sig_group == "Up") %>% - arrange(desc(log2FoldChange)) %>% head(10) - top_down <- volc_df %>% filter(sig_group == "Down") %>% - arrange(log2FoldChange) %>% head(10) - label_genes <- bind_rows(top_up, top_down) - - p <- ggplot(volc_df, aes(x = log2FoldChange, y = neg_log10p, color = sig_group)) + - geom_point(size = 0.8, alpha = 0.6) + - scale_color_manual(values = c("NS" = "grey70", "Up" = "#e41a1c", "Down" = "#377eb8"), - name = "") + - geom_hline(yintercept = -log10(padj_thr), linetype = "dashed", color = "grey40", linewidth = 0.3) + - geom_vline(xintercept = c(-lfc_thr, lfc_thr), linetype = "dashed", color = "grey40", linewidth = 0.3) + - ggrepel::geom_text_repel( - data = label_genes, aes(label = gene), - size = 2.5, max.overlaps = 20, color = "black", - segment.size = 0.2, segment.color = "grey50" - ) + - labs(title = contrast_name, subtitle = paste0(prefix, " — Wald test"), - x = "log2 Fold Change", y = "-log10(padj)") + - theme_bw() + - theme(legend.position = "right") - - ggsave(file.path(pairwise_dir, paste0("volcano_", contrast_name, ".png")), - plot = p, width = 9, height = 7, dpi = res_dpi) - }, error = function(e) { - message(" Volcano failed: ", conditionMessage(e)) - }) - - # ComplexHeatmap - tryCatch({ - sig_genes <- res_df %>% filter(padj < padj_thr) - top_up <- sig_genes %>% filter(log2FoldChange > 0) %>% head(de_n_genes) %>% pull(gene) - top_down <- sig_genes %>% filter(log2FoldChange < 0) %>% head(de_n_genes) %>% pull(gene) - top_genes <- c(top_up, top_down) - top_genes <- top_genes[top_genes %in% rownames(assay(vsd))] - - if (length(top_genes) >= 2) { - mat <- assay(vsd)[top_genes, , drop = FALSE] - mat_scaled <- t(scale(t(mat))) - - draw_complex_heatmaps( - mat_scaled, meta, condition_col, top_genes, - contrast_name, pairwise_dir, conditions, region_colors - ) - } - }, error = function(e) { - message(" Heatmap failed: ", conditionMessage(e)) - }) - + write_volcano(result_df, contrast_name, prefix, pairwise_dir) + write_de_heatmap(result_df, fit$vsd, fit_meta, condition_col, contrast_name, pairwise_dir) }, error = function(e) { - message(" Wald test failed: ", conditionMessage(e)) - summary_rows[[contrast_name]] <- tibble( - contrast = contrast_name, n_tested = 0, - n_significant = 0, n_up = 0, n_down = 0, error = conditionMessage(e) + message(" Contrast skipped: ", conditionMessage(e)) + summary_rows[[contrast_name]] <<- tibble( + contrast = contrast_name, numerator = numerator, denominator = denominator, + n_tested = 0L, n_significant = 0L, n_up = 0L, n_down = 0L, + error = conditionMessage(e) ) }) } @@ -445,389 +534,31 @@ run_de_for_prefix <- function(prefix, condition_col, sample_col, out_base) { write.table(file.path(out_base, "pairwise_summary.tsv"), sep = "\t", quote = FALSE, row.names = FALSE) } - - # ── DEGpatterns ──────────────────────────────────────────────────────── - if (!is.null(dds_lrt)) { - tryCatch({ - res_lrt_all <- results(dds_lrt) - sig_genes <- rownames(res_lrt_all)[which(res_lrt_all$padj < padj_thr)] - - if (length(sig_genes) >= 10) { - message(" DEGpatterns on ", length(sig_genes), " LRT-significant genes …") - sig_mat <- assay(vsd)[sig_genes, , drop = FALSE] - - deg_dir <- file.path(out_base, "DEGpatterns") - dir.create(deg_dir, recursive = TRUE, showWarnings = FALSE) - - clusters <- degPatterns( - sig_mat, metadata = meta, time = condition_col, plot = FALSE - ) - - if (!is.null(clusters$df)) { - # ── Remap cluster ids to a contiguous 1..N (degPatterns skips ids) - norm <- clusters$normalized - old_levels <- sort(unique(norm$cluster)) - remap <- setNames(seq_along(old_levels), as.character(old_levels)) - # keep cluster as INTEGER — degPlotCluster joins on it internally and a - # factor here triggers "Can't join with ". - norm$cluster <- as.integer(remap[as.character(norm$cluster)]) - - cluster_df <- clusters$df - cluster_df$cluster <- as.integer(remap[as.character(cluster_df$cluster)]) - cluster_df <- cluster_df %>% arrange(cluster, genes) - - # Base plot from degPlotCluster (now facetted by 1..N) - p <- degPlotCluster( - norm, time = condition_col, color = condition_col, points = TRUE - ) - - # Apply custom colors if available - if (length(region_colors) > 0) { - p <- p + - ggplot2::aes(col = .data[[condition_col]]) + - ggplot2::scale_color_manual(values = region_colors) - } - - # Loess trend + rotated x labels - p <- p + - ggplot2::geom_smooth( - mapping = ggplot2::aes(x = .data[[condition_col]], y = value, group = 1), - method = "loess", color = "black", se = FALSE, linewidth = 1.2 - ) + - ggplot2::theme( - axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5) - ) - - ggsave(file.path(deg_dir, "DEGpatterns_groups.png"), # renamed - plot = p, width = 14, height = 10, dpi = max(res_dpi, 500L)) - message(" DEGpatterns: ", length(old_levels), " groups (renumbered 1..", - length(old_levels), ")") - - # Gene groups TSV — use the remapped, contiguous group numbers - tryCatch({ - grouped_genes <- split(cluster_df$genes, - paste0("group_", cluster_df$cluster)) - max_len <- max(sapply(grouped_genes, length)) - gene_df <- as.data.frame(do.call(cbind, lapply(grouped_genes, function(x) { - c(x, rep("", max_len - length(x))) - }))) - colnames(gene_df) <- gsub("[^A-Za-z0-9_]", "_", colnames(gene_df)) - colnames(gene_df) <- gsub("_+", "_", colnames(gene_df)) - colnames(gene_df) <- gsub("_$", "", colnames(gene_df)) - write.table(gene_df, file.path(deg_dir, "gene_groups.tsv"), - sep = "\t", quote = FALSE, row.names = FALSE) - message(" Gene groups saved: ", ncol(gene_df), " groups") - }, error = function(e) { - message(" Gene groups export failed: ", conditionMessage(e)) - }) - - # ── Per-group heatmaps: the genes of every DEGpatterns group ───── - # Reuses the DE heatmap helper (unsplit + region-split, z-scored vsd), - # written alongside the group plots in DEGpatterns/. - tryCatch({ - for (g in sort(unique(cluster_df$cluster))) { - genes_g <- cluster_df$genes[cluster_df$cluster == g] - genes_g <- genes_g[genes_g %in% rownames(assay(vsd))] - if (length(genes_g) >= 2) { - mat_g_scaled <- t(scale(t(assay(vsd)[genes_g, , drop = FALSE]))) - draw_complex_heatmaps( - mat_g_scaled, meta, condition_col, genes_g, - paste0("group_", g), deg_dir, conditions, region_colors - ) - } - } - message(" DEGpatterns per-group heatmaps written") - }, error = function(e) { - message(" DEGpatterns per-group heatmaps failed: ", conditionMessage(e)) - }) - } - } else { - message(" Too few LRT-significant genes for DEGpatterns (", length(sig_genes), ")") - } - }, error = function(e) { - message(" DEGpatterns failed: ", conditionMessage(e)) - }) - } - - # ── One-vs-rest comparisons + unique marker identification ────────── - message(" One-vs-rest comparisons …") - ovr_dir <- file.path(out_base, "one_vs_rest") - dir.create(ovr_dir, recursive = TRUE, showWarnings = FALSE) - - ovr_results <- list() # store results for cross-comparison - - for (focal in conditions) { - focal_safe <- gsub("[^A-Za-z0-9_]", "_", focal) - message(" ", focal, " vs rest") - - tryCatch({ - # Create binary factor: focal vs rest - meta_ovr <- meta - meta_ovr[["ovr_group"]] <- ifelse(meta_ovr[[condition_col]] == focal, focal, "rest") - meta_ovr[["ovr_group"]] <- factor(meta_ovr[["ovr_group"]], levels = c("rest", focal)) - - dds_ovr <- DESeqDataSetFromMatrix( - countData = counts_t, colData = meta_ovr, design = ~ ovr_group - ) - dds_ovr <- DESeq(dds_ovr, test = "Wald") - res_ovr <- results(dds_ovr, contrast = c("ovr_group", focal, "rest")) - - res_ovr_df <- as.data.frame(res_ovr) %>% - rownames_to_column("gene") %>% - arrange(desc(log2FoldChange)) - - write.table(res_ovr_df, file.path(ovr_dir, paste0(focal_safe, "_vs_rest.tsv")), - sep = "\t", quote = FALSE, row.names = FALSE) - - n_sig <- sum(res_ovr_df$padj < padj_thr, na.rm = TRUE) - n_up <- sum(res_ovr_df$padj < padj_thr & res_ovr_df$log2FoldChange > 0, na.rm = TRUE) - message(" ", n_sig, " sig (", n_up, " up)") - - # Store for cross-comparison - ovr_results[[focal]] <- res_ovr_df - - # Volcano plot (same style as pairwise) - tryCatch({ - volc_df <- res_ovr_df %>% - filter(!is.na(padj) & !is.na(log2FoldChange)) %>% - mutate( - sig_group = case_when( - padj >= padj_thr ~ "NS", - log2FoldChange > lfc_thr ~ "Up", - log2FoldChange < -lfc_thr ~ "Down", - TRUE ~ "NS" - ), - neg_log10p = -log10(pmax(padj, 1e-300)) - ) - top_up <- volc_df %>% filter(sig_group == "Up") %>% arrange(desc(log2FoldChange)) %>% head(10) - top_down <- volc_df %>% filter(sig_group == "Down") %>% arrange(log2FoldChange) %>% head(10) - label_genes <- bind_rows(top_up, top_down) - - p <- ggplot(volc_df, aes(x = log2FoldChange, y = neg_log10p, color = sig_group)) + - geom_point(size = 0.8, alpha = 0.6) + - scale_color_manual(values = c("NS" = "grey70", "Up" = "#e41a1c", "Down" = "#377eb8"), name = "") + - geom_hline(yintercept = -log10(padj_thr), linetype = "dashed", color = "grey40", linewidth = 0.3) + - geom_vline(xintercept = c(-lfc_thr, lfc_thr), linetype = "dashed", color = "grey40", linewidth = 0.3) + - ggrepel::geom_text_repel(data = label_genes, aes(label = gene), size = 2.5, - max.overlaps = 20, color = "black", - segment.size = 0.2, segment.color = "grey50") + - labs(title = paste0(focal, " vs rest"), x = "log2 Fold Change", y = "-log10(padj)") + - theme_bw() + theme(legend.position = "right") - - ggsave(file.path(ovr_dir, paste0("volcano_", focal_safe, "_vs_rest.png")), - plot = p, width = 9, height = 7, dpi = res_dpi) - }, error = function(e) { - message(" OVR volcano failed: ", conditionMessage(e)) - }) - - # Top DE gene heatmap - tryCatch({ - sig_genes_ovr <- res_ovr_df %>% filter(padj < padj_thr) - top_up_genes <- sig_genes_ovr %>% filter(log2FoldChange > 0) %>% head(de_n_genes) %>% pull(gene) - top_down_genes <- sig_genes_ovr %>% filter(log2FoldChange < 0) %>% head(de_n_genes) %>% pull(gene) - top_genes_ovr <- c(top_up_genes, top_down_genes) - top_genes_ovr <- top_genes_ovr[top_genes_ovr %in% rownames(assay(vsd))] - - if (length(top_genes_ovr) >= 2) { - draw_complex_heatmaps( - t(scale(t(assay(vsd)[top_genes_ovr, , drop = FALSE]))), - meta, condition_col, top_genes_ovr, - paste0(focal_safe, "_vs_rest"), ovr_dir, conditions, region_colors - ) - } - }, error = function(e) { - message(" OVR heatmap failed: ", conditionMessage(e)) - }) - - }, error = function(e) { - message(" OVR test failed for '", focal, "': ", conditionMessage(e)) - }) - } - - # ── Unique marker identification ─────────────────────────────────── - # A gene is a unique marker for level X if it is: - # - significantly UP in X vs rest (padj < thr, log2FC > lfc_thr) - # - NOT significantly UP in any other level vs rest - if (length(ovr_results) >= 2) { - message(" Identifying unique markers …") - tryCatch({ - # For each level, get the set of significantly upregulated genes - up_gene_sets <- lapply(ovr_results, function(df) { - df %>% filter(padj < padj_thr & log2FoldChange > lfc_thr) %>% pull(gene) - }) - - unique_markers <- list() - for (lvl in names(up_gene_sets)) { - others <- setdiff(names(up_gene_sets), lvl) - other_genes <- unique(unlist(up_gene_sets[others])) - unique_to_lvl <- setdiff(up_gene_sets[[lvl]], other_genes) - - if (length(unique_to_lvl) > 0) { - # Get the DE stats for these genes - lvl_df <- ovr_results[[lvl]] %>% - filter(gene %in% unique_to_lvl) %>% - arrange(desc(log2FoldChange)) - unique_markers[[lvl]] <- lvl_df - message(" ", lvl, ": ", nrow(lvl_df), " unique markers") - } else { - message(" ", lvl, ": 0 unique markers") - } - } - - # Save as combined TSV - if (length(unique_markers) > 0) { - all_markers <- bind_rows(unique_markers, .id = "level") %>% - arrange(level, desc(log2FoldChange)) - write.table(all_markers, file.path(ovr_dir, "unique_markers.tsv"), - sep = "\t", quote = FALSE, row.names = FALSE) - - # Also save as column-format (one column per level, padded) - max_len <- max(sapply(unique_markers, nrow)) - marker_cols <- lapply(unique_markers, function(df) { - c(df$gene, rep("", max_len - nrow(df))) - }) - marker_df <- as.data.frame(marker_cols) - colnames(marker_df) <- gsub("[^A-Za-z0-9_]", "_", colnames(marker_df)) - write.table(marker_df, file.path(ovr_dir, "unique_markers_by_level.tsv"), - sep = "\t", quote = FALSE, row.names = FALSE) - } - - # Summary barplot: number of unique markers per level - if (length(unique_markers) > 0) { - counts_bar <- tibble( - level = factor(names(unique_markers), levels = conditions), - n_markers = sapply(unique_markers, nrow) - ) - p <- ggplot(counts_bar, aes(x = level, y = n_markers, fill = level)) + - geom_col(show.legend = FALSE) + - geom_text(aes(label = n_markers), vjust = -0.3, size = 3) + - labs(title = "Unique upregulated markers per level", - x = "", y = "Number of unique markers") + - theme_bw() + - theme(axis.text.x = element_text(angle = 45, hjust = 1)) - ggsave(file.path(ovr_dir, "unique_markers_barplot.png"), - plot = p, width = max(6, length(unique_markers) * 1.2), height = 5, dpi = res_dpi) - } - - # Heatmap: ALL unique marker genes (rows), samples (cols), row-split by - # the level each gene is unique to. May be tall — that is intended. - if (length(unique_markers) > 0) { - tryCatch({ - # gene -> level map, preserving level order in `conditions` - gene_level <- unlist(lapply(names(unique_markers), function(lvl) - setNames(rep(lvl, nrow(unique_markers[[lvl]])), - unique_markers[[lvl]]$gene))) - genes_all <- names(gene_level) - genes_all <- genes_all[genes_all %in% rownames(assay(vsd))] - - if (length(genes_all) >= 2) { - mat <- assay(vsd)[genes_all, , drop = FALSE] - mat_scaled <- t(scale(t(mat))) - - # order columns by region level - col_ord <- order(factor(meta[[condition_col]], - levels = if (length(conditions) > 0) conditions - else unique(meta[[condition_col]]))) - mat_scaled <- mat_scaled[, col_ord, drop = FALSE] - meta_ord <- meta[col_ord, , drop = FALSE] - - # drop zero-variance genes (all-NaN rows after scaling) - finite_rows <- apply(mat_scaled, 1, function(z) all(is.finite(z))) - mat_scaled <- mat_scaled[finite_rows, , drop = FALSE] - - if (nrow(mat_scaled) >= 2) { - row_split <- factor(gene_level[rownames(mat_scaled)], - levels = if (length(conditions) > 0) conditions - else unique(gene_level)) - - col_anno <- build_top_annotation(meta_ord, condition_col, region_colors) - - col_fun <- colorRamp2(c(-2, 0, 2), c("blue", "white", "red")) - n_row <- nrow(mat_scaled); n_col <- ncol(mat_scaled) - res_png <- res_dpi - extra_in <- 3.0 # title + anno + legend - max_px <- 30000 # png/cairo device ceiling - row_mm <- 3 # desired mm per gene row - body_mm <- max(n_row * row_mm, 40) - dev_h_px <- round((body_mm / 25.4 + extra_in) * res_png) - if (dev_h_px > max_px) { # too tall for the device: - dev_h_px <- max_px # cap and shrink rows to fit - body_mm <- (max_px / res_png - extra_in) * 25.4 - } - row_fs <- if (n_row > 400) 3 else 5 # smaller font when dense - - ht_uniq <- Heatmap( - mat_scaled, name = "Z-score", col = col_fun, - top_annotation = col_anno, - height = unit(body_mm, "mm"), - row_split = row_split, - cluster_row_slices = FALSE, - cluster_columns = FALSE, - cluster_rows = TRUE, - show_row_names = (n_row <= 1500), # labels unreadable beyond this - show_column_names = TRUE, - row_names_gp = gpar(fontsize = row_fs), - column_names_gp = gpar(fontsize = 7), - row_title_gp = gpar(fontsize = 8), - column_title = "Unique markers per level (all genes)", - heatmap_legend_param = list(title = "Z-score") - ) - png(file.path(ovr_dir, "unique_markers_heatmap.png"), - width = round(max(900, n_col * 55) * px_scale), - height = dev_h_px, res = res_png) - draw(ht_uniq, merge_legend = TRUE) - dev.off() - message(" Unique markers heatmap: ", n_row, " genes") - } - } - }, error = function(e) { - message(" Unique markers heatmap failed: ", conditionMessage(e)) - }) - } - - }, error = function(e) { - message(" Unique marker identification failed: ", conditionMessage(e)) - }) - } - - return(summary_rows) + invisible(summary_rows) } -# ── Main ───────────────────────────────────────────────────────────────────── manifest_file <- file.path(agg_dir, "manifest.tsv") - if (!file.exists(manifest_file)) { - message("No manifest.tsv found.") writeLines("No pseudobulk matrices found.", file.path(results_dir, "SKIPPED_no_manifest.txt")) - if (exists("log_err_con")) { sink(type = "message"); close(log_err_con) } - if (exists("log_out_con")) { sink(); close(log_out_con) } - quit(save = "no", status = 0) -} - -manifest <- read.delim(manifest_file) -message("Manifest: ", nrow(manifest), " subgroups") - -for (i in seq_len(nrow(manifest))) { - row <- manifest[i, ] - prefix <- row$prefix - condition_col <- row$condition_col - sample_col <- row$sample_col - - message("\n--- Processing: ", prefix, " ---") - - out_base <- if (prefix == "pooled") results_dir else file.path(results_dir, prefix) - tryCatch( - run_de_for_prefix(prefix, condition_col, sample_col, out_base), - error = function(e) { - message(" FAILED for '", prefix, "': ", conditionMessage(e)) - dir.create(out_base, recursive = TRUE, showWarnings = FALSE) - writeLines(conditionMessage(e), file.path(out_base, "ERROR.txt")) - } - ) +} else { + manifest <- read.delim(manifest_file) + message("Manifest: ", nrow(manifest), " subgroup(s)") + for (index in seq_len(nrow(manifest))) { + row <- manifest[index, ] + prefix <- row$prefix + out_base <- if (prefix == "pooled") results_dir else file.path(results_dir, prefix) + tryCatch( + run_de_for_prefix(prefix, row$condition_col, row$sample_col, out_base), + error = function(e) { + message(" FAILED for '", prefix, "': ", conditionMessage(e)) + dir.create(out_base, recursive = TRUE, showWarnings = FALSE) + writeLines(conditionMessage(e), file.path(out_base, "ERROR.txt")) + } + ) + } } -message("\n=== Pseudobulk DE complete ===") - +message("=== Pseudobulk DE complete ===") if (exists("log_err_con")) { sink(type = "message"); close(log_err_con) } if (exists("log_out_con")) { sink(); close(log_out_con) }