From c87e91b9cb6007bcbeb09cced4dc21a7f0c24b97 Mon Sep 17 00:00:00 2001 From: Sergio-ote Date: Wed, 9 Sep 2026 13:12:22 +0200 Subject: [PATCH] fix: handle annotation availability and external labels, reporting only available cell annotations --- config/README.md | 11 +- config/config.yaml | 4 +- docs/rules.md | 4 +- workflow/Snakefile | 3 +- workflow/rules/common.smk | 16 ++- workflow/schemas/config.schema.yaml | 10 +- workflow/scripts/annotate_cells.py | 181 +++++++++++++++----------- workflow/scripts/annotation_utils.py | 93 +++++++++++++ workflow/scripts/integrate_samples.py | 27 ++-- workflow/scripts/preprocess_umap.py | 38 ++++-- workflow/scripts/sample_report.py | 22 ++-- workflow/scripts/validate_input.py | 58 ++++++--- 12 files changed, 324 insertions(+), 143 deletions(-) create mode 100644 workflow/scripts/annotation_utils.py diff --git a/config/README.md b/config/README.md index 7dc080a..89fbb4f 100644 --- a/config/README.md +++ b/config/README.md @@ -325,7 +325,16 @@ SpaceBlocks allows you to input: - **Externally assembled h5ad AnnData objects** — run `mode: decoupled` and point `contract_dir` at the directory of pre-built contract h5ads. The heads are skipped; the core validates and analyses them directly. - **Pre-computed clusters / annotations** — set `use_precomputed_clusters: true` and `precomputed_metadata_dir` to reuse Leiden clusters and metadata; for niches, set `spatial_niches.use_precomputed: true` with `spatial_niches.niche_dir`. -- **Externally annotated data** — set `external_annotation.enabled: true` with `external_annotation.column`, and choose whether to keep or discard unannotated barcodes downstream via `external_annotation.keep_unannotated` (see [section 5](#5-key-sections)). +- **Externally annotated data** — set `external_annotation.enabled: true` and use + `external_annotation.column` for the source column. Provide it in each + `metadata_{sample}.tsv`, or, in `mode: decoupled`, directly in each contract + h5ad's `obs`. The source column may have any name; SpaceBlocks copies it to the + canonical `cell_type_external` column used by downstream rules. Choose + whether to keep or discard unannotated barcodes via + `external_annotation.keep_unannotated` (see [section 5](#5-key-sections)). With + `keep_unannotated: true`, a mixture of real labels and `Unannotated` remains an + active annotation and both are shown in overview and composition plots. Only a + column containing no real labels is omitted from annotation plots and reports. These files are generated during the run, and can be shared with minimum effort to reproduce downstream results from raw data. diff --git a/config/config.yaml b/config/config.yaml index 8123b70..97818dd 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -190,10 +190,10 @@ analysis: snakemake_cell_markers: config/snakemake_cell_markers.tsv # (default) cluster → cell-type TSV (fill template from generate_annotation_template) cluster_annotations: config/cluster_annotations.tsv -# (optional) use external tsv with annotations +# (optional) external labels from metadata TSVs, or from contract h5ad obs in decoupled mode external_annotation: enabled: false - column: manual_celltype + column: manual_celltype # arbitrary source name; normalized to cell_type_external keep_unannotated: true ingest_ref: "/path/to/ingest_ref_data/bcc_and_normal_Garnier_PNAS_2024.h5ad" # optional annotated scRNA-seq reference (ingest_ref) ingest_ref_label_key: 04_celltypes diff --git a/docs/rules.md b/docs/rules.md index dcefff6..309277f 100644 --- a/docs/rules.md +++ b/docs/rules.md @@ -155,7 +155,7 @@ This second CoreBlock includes rules from sample annotation to integration and d #### `annotate_cells` Adds cell-type annotations from a manually curated TSV (cluster → cell type), or from an external annotation TSV (barcode → cell type). -Writes `adata_{sample}_annotated.h5ad` plus annotation plots and composition barplots. +Writes `adata_{sample}_annotated.h5ad` plus annotation plots and composition barplots. The output records which annotation columns contain real labels; placeholder-only columns such as an all-`Unannotated` `cell_type_tsv` are retained for compatibility but are not treated as available annotations. The external source column configured by `external_annotation.column` may have any name and is copied to the canonical `cell_type_external` column. It can come from `metadata_{sample}.tsv` or, in decoupled mode without external metadata, directly from the contract h5ad's `obs`. When `keep_unannotated: true`, mixed real and `Unannotated` labels remain active and the unannotated cells are retained in overview and composition plots. #### `integrate_samples` Concatenates all annotated samples, runs Harmony batch correction on the specified variable (`integration.integrate_key`), geosketch subsampling, and composition barplots. @@ -189,7 +189,7 @@ Generates integrated (Harmony) and unintegrated results, with UMAPs coloured by #### `sample_report` Assembles a compact, multi-page PDF with graphical outputs per sample, ideal to share with internal and external collaborators. -Each page contains the most important data-descriptive information: cluster / annotation / region / niche UMAPs and matching spatial maps, bar plots for absolute and relative composition, and a dot plot with the top-10 expressed markers per cell type. +Each page contains the most important data-descriptive information: cluster / annotation / region / niche UMAPs and matching spatial maps, bar plots for absolute and relative composition, and a dot plot with the top-10 expressed markers per cell type. Only annotation sources containing real labels are shown; in runs without a manual TSV annotation, the external or ingest annotation also drives cell-type composition. --- diff --git a/workflow/Snakefile b/workflow/Snakefile index f3df267..1b2627f 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -239,7 +239,8 @@ RUN_DE = bool(ANALYSIS.get("run_pseudobulk_de", False)) RUN_LEIDEN_ANALYSIS = bool(ANALYSIS.get("run_leiden_analysis", True)) -# External annotation: fail fast if enabled but the column is not present in every sample's metadata +# External annotation: validate metadata headers now, or decoupled h5ad obs per sample +# in validate_input when no external metadata directory is configured. check_external_annotation() # Consistency: analysing the external annotation downstream requires it to be produced. diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 64e7c0e..34494f8 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -58,9 +58,9 @@ def _expand_base_dir(node, root="config"): # ── Small derivations called from the Snakefile globals block ──────────────── def check_external_annotation(): """Fail fast (at parse time, before any job) when external_annotation is enabled - but not fully in place: the annotation column must be present in the metadata for - EVERY sample. External labels can only come from precomputed_metadata_dir (the - pipeline's own metadata does not carry them), so that dir is required.""" + but not fully in place. Metadata TSV columns are checked here; in decoupled mode, + an absent metadata directory means the configured column will instead be checked + in each contract h5ad by validate_input.""" cfg = config.get("external_annotation", {}) or {} if not cfg.get("enabled", False): return @@ -71,10 +71,13 @@ def check_external_annotation(): ) meta_dir = config.get("precomputed_metadata_dir", "") or "" if not meta_dir: + if IS_DECOUPLED: + return sys.exit( "[config error] external_annotation.enabled requires 'precomputed_metadata_dir' " f"to point at a directory of metadata_{{sample}}.tsv files carrying the '{col}' " - "column (the pipeline's own metadata does not contain external labels)." + "column. In decoupled mode only, the column may instead be stored directly " + "in each contract h5ad's obs." ) missing_file, missing_col = [], [] for s in SAMPLE_IDS: @@ -470,10 +473,11 @@ def _geojson_validation_report(sample): def _external_metadata_input(wildcards): """Tracked per-sample external metadata, or no input when the feature is off.""" cfg = config.get("external_annotation", {}) or {} - if not cfg.get("enabled", False): + metadata_dir = config.get("precomputed_metadata_dir", "") or "" + if not cfg.get("enabled", False) or not metadata_dir: return [] return os.path.join( - config.get("precomputed_metadata_dir", "") or "", + metadata_dir, f"metadata_{wildcards.sample}.tsv", ) diff --git a/workflow/schemas/config.schema.yaml b/workflow/schemas/config.schema.yaml index 0dd450d..c34d0a9 100644 --- a/workflow/schemas/config.schema.yaml +++ b/workflow/schemas/config.schema.yaml @@ -501,13 +501,19 @@ properties: default: false column: type: string - description: obs/metadata column holding the external label. + description: >- + Source column holding the external label, either in + metadata_{sample}.tsv or, in decoupled mode without external metadata, + directly in the contract h5ad obs. The name is arbitrary; SpaceBlocks + copies it to the canonical cell_type_external obs column downstream. keep_unannotated: type: boolean default: true description: >- false → keep only externally-annotated cells and skip pipeline QC (external - labels drive QC); true → normal QC, labels overlaid at annotation. + labels drive QC); true → normal QC, labels overlaid at annotation and + cells labelled Unannotated retained. A mixed labelled/Unannotated column + remains available for plots; only placeholder-only columns are omitted. contract: type: ["object", "null"] diff --git a/workflow/scripts/annotate_cells.py b/workflow/scripts/annotate_cells.py index 38ffbe6..3a394b0 100644 --- a/workflow/scripts/annotate_cells.py +++ b/workflow/scripts/annotate_cells.py @@ -32,7 +32,11 @@ _here = os.getcwd() sys.path.insert(0, _here) from composition_barplots import composition_pair, find_niche_column -from external_metadata import normalized_labels, optional_input_path, read_cell_metadata +from annotation_utils import ( + normalize_annotation_labels, + record_active_annotation_columns, +) +from external_metadata import optional_input_path, read_cell_metadata # ── Logging ────────────────────────────────────────────────────────────────── @@ -256,77 +260,98 @@ def generate_annotation_plots(adata, annot_key, label, plots_dir, sample_id, annot_col, resolution = find_sample_column(annot_df, sample_id) if annot_col is None: - log.warning("Sample '%s' not found in TSV. All cells → 'Unannotated'.", sample_id) - adata.obs["cell_type_tsv"] = "Unannotated" - Path(out_adata_path).parent.mkdir(parents=True, exist_ok=True) - adata.write(out_adata_path) - sys.exit(0) - - log.info("Using column '%s' (resolution: %s)", annot_col, - resolution if resolution else "not specified") - - if resolution is not None: - leiden_key = f"leiden_{resolution.replace('.', '_')}" + log.warning( + "Sample '%s' not found in cluster_annotations TSV. " + "cell_type_tsv set to 'Unannotated'; continuing with any " + "external/ingest annotation.", + sample_id, + ) + adata.obs["cell_type_tsv"] = pd.Categorical( + ["Unannotated"] * adata.n_obs + ) else: - leiden_cols = [c for c in adata.obs.columns if c.startswith("leiden_")] - if leiden_cols: - leiden_key = leiden_cols[0] - else: - raise ValueError("No leiden columns found in adata.") + log.info("Using column '%s' (resolution: %s)", annot_col, + resolution if resolution else "not specified") - if leiden_key not in adata.obs.columns: - available = [c for c in adata.obs.columns if c.startswith("leiden_")] - raise ValueError( - f"Leiden column '{leiden_key}' not found. Available: {available}." + if resolution is not None: + leiden_key = f"leiden_{resolution.replace('.', '_')}" + else: + leiden_cols = [c for c in adata.obs.columns if c.startswith("leiden_")] + if leiden_cols: + leiden_key = leiden_cols[0] + else: + raise ValueError("No leiden columns found in adata.") + + if leiden_key not in adata.obs.columns: + available = [c for c in adata.obs.columns if c.startswith("leiden_")] + raise ValueError( + f"Leiden column '{leiden_key}' not found. Available: {available}." + ) + + cluster_to_celltype = annot_df[annot_col].dropna().to_dict() + cluster_to_celltype = {str(k): str(v) for k, v in cluster_to_celltype.items()} + log.info("Cluster → cell type mapping (%d entries): %s", + len(cluster_to_celltype), cluster_to_celltype) + + adata.obs["cell_type_tsv"] = ( + adata.obs[leiden_key].astype(str) + .map(cluster_to_celltype) + .fillna("Unannotated") + .astype("category") ) - cluster_to_celltype = annot_df[annot_col].dropna().to_dict() - cluster_to_celltype = {str(k): str(v) for k, v in cluster_to_celltype.items()} - log.info("Cluster → cell type mapping (%d entries): %s", - len(cluster_to_celltype), cluster_to_celltype) - - adata.obs["cell_type_tsv"] = ( - adata.obs[leiden_key].astype(str) - .map(cluster_to_celltype) - .fillna("Unannotated") - .astype("category") - ) - - tsv_counts = adata.obs["cell_type_tsv"].value_counts() - log.info("TSV annotation distribution:\n%s", tsv_counts.to_string()) + tsv_counts = adata.obs["cell_type_tsv"].value_counts() + log.info("TSV annotation distribution:\n%s", tsv_counts.to_string()) # ── 2. External annotation (from metadata TSV column) ──────────────── # When enabled the pipeline MUST find the column for this sample — a missing # source or column is a hard error (fail loudly rather than silently skipping), # so a run only proceeds when the external labels are in place for every sample. - ext_enabled = False if isinstance(EXT_ANNOT_CFG, dict) and EXT_ANNOT_CFG.get("enabled", False): ext_col = EXT_ANNOT_CFG.get("column", "") if not ext_col: raise ValueError("external_annotation.enabled is true but no 'column' is set in config.") - if not external_metadata_path: - raise FileNotFoundError( - f"external_annotation enabled but no metadata found for sample " - f"'{sample_id}'.") + # Prefer the explicitly supplied metadata TSV. In decoupled mode, a custom + # source column may instead already be present in the contract h5ad and is + # preserved by preprocessing. + if external_metadata_path: + _ext_source = external_metadata_path + log.info( + "Loading external annotation from '%s' in %s …", + ext_col, + _ext_source, + ) + saved = read_cell_metadata( + _ext_source, required_columns=(ext_col,) + ) + external_labels = saved[ext_col].reindex( + adata.obs_names.astype(str) + ) + elif ext_col in adata.obs.columns: + _ext_source = f"processed contract h5ad obs['{ext_col}']" + log.info("Loading external annotation from %s …", _ext_source) + external_labels = adata.obs[ext_col] + else: + raise ValueError( + f"external_annotation column '{ext_col}' not found for sample " + f"'{sample_id}': no matching precomputed metadata column or h5ad " + "obs column is available" + ) + adata.obs["cell_type_external"] = normalize_annotation_labels( + external_labels + ).astype("category") + adata.uns["external_annotation_source_column"] = ext_col + n_external_unannotated = int( + (adata.obs["cell_type_external"] == "Unannotated").sum() + ) log.info( - "Loading external annotation from '%s' in %s …", + " External annotation '%s' copied to 'cell_type_external': %d types " + "(%d cells unmatched → Unannotated)", ext_col, - external_metadata_path, - ) - saved = read_cell_metadata( - external_metadata_path, required_columns=(ext_col,) - ) - labels = normalized_labels(saved, ext_col) - adata.obs["cell_type_external"] = ( - labels.reindex(adata.obs_names.astype(str)) - .fillna("Unannotated") - .astype("category") + adata.obs["cell_type_external"].nunique(), + n_external_unannotated, ) - ext_enabled = True - log.info(" External annotation: %d types (%d cells unmatched → Unannotated)", - adata.obs["cell_type_external"].nunique(), - int((adata.obs["cell_type_external"] == "Unannotated").sum())) # ── 2b. Spatial niche labels (from spatial_niches rule / external) ─── _ni = getattr(snakemake.input, "spatial_niche", "") @@ -363,11 +388,25 @@ def generate_annotation_plots(adata, annot_key, label, plots_dir, sample_id, # ── 3. Plots ───────────────────────────────────────────────────────── log.info("Generating annotation plots …") + active_annotations = record_active_annotation_columns(adata) + log.info("Active annotation columns: %s", active_annotations or "none") + if ( + isinstance(EXT_ANNOT_CFG, dict) + and EXT_ANNOT_CFG.get("enabled", False) + and "cell_type_external" not in active_annotations + ): + log.warning( + "External annotation is enabled, but '%s' contains no non-placeholder " + "labels matching this sample. No external-annotation plots will be made.", + EXT_ANNOT_CFG.get("column", ""), + ) + # Apply custom palettes if configured if isinstance(ANNOTATION_COLORS, dict): - for obs_key in ["cell_type_tsv", "cell_type_ingest", "cell_type_external"]: + for obs_key in active_annotations: cd = ANNOTATION_COLORS.get(obs_key, {}) if cd and obs_key in adata.obs.columns: + adata.obs[obs_key] = adata.obs[obs_key].astype("category") cats = adata.obs[obs_key].cat.categories adata.uns[f"{obs_key}_colors"] = [cd.get(str(c), "#cccccc") for c in cats] @@ -378,22 +417,20 @@ def generate_annotation_plots(adata, annot_key, label, plots_dir, sample_id, REGION_COLORS.get(str(c), "#cccccc") for c in cats ] - generate_annotation_plots(adata, "cell_type_tsv", "tsv", plots_dir, - sample_id, DE_N_GENES, library_id) - - # If ingest labels exist, plot those too - if has_ingest: - generate_annotation_plots(adata, "cell_type_ingest", "ingest", plots_dir, - sample_id, DE_N_GENES, library_id) - - if ext_enabled: - generate_annotation_plots(adata, "cell_type_external", "external", plots_dir, - sample_id, DE_N_GENES, library_id) + annotation_labels = { + "cell_type_tsv": "tsv", + "cell_type_ingest": "ingest", + "cell_type_external": "external", + "cell_type_refined": "refined", + } + for annot_col in active_annotations: + generate_annotation_plots( + adata, annot_col, annotation_labels.get(annot_col, annot_col), + plots_dir, sample_id, DE_N_GENES, library_id, + ) # Side-by-side comparison (dynamic number of panels) - annot_cols_present = [c for c in ["cell_type_tsv", "cell_type_ingest", - "cell_type_external"] - if c in adata.obs.columns] + annot_cols_present = active_annotations n_panels = len(annot_cols_present) if n_panels >= 2: fig, axes = plt.subplots(1, n_panels, figsize=(8 * n_panels, 6), @@ -423,9 +460,7 @@ def generate_annotation_plots(adata, annot_key, label, plots_dir, sample_id, log.warning("Spatial comparison plot failed: %s", e) # ── Per-sample composition barplots (sample / region / niche) ──────── - annot_cols_present = [c for c in ["cell_type_tsv", "cell_type_ingest", - "cell_type_external"] - if c in adata.obs.columns] + annot_cols_present = active_annotations # sample column (one bar per sample; here a single sample). Per-sample adatas # carry "sample"; "sample_batch" is an integration-time concat label and does diff --git a/workflow/scripts/annotation_utils.py b/workflow/scripts/annotation_utils.py new file mode 100644 index 0000000..34e4b94 --- /dev/null +++ b/workflow/scripts/annotation_utils.py @@ -0,0 +1,93 @@ +"""Shared helpers for deciding which cell-type annotations are meaningful. + +The names below are SpaceBlocks' canonical ``adata.obs`` columns. In particular, +``external_annotation.column`` may name any column in the source metadata TSV; +``annotate_cells.py`` copies that source column to ``cell_type_external`` before +these helpers are called. +""" + +from __future__ import annotations + + +ACTIVE_ANNOTATION_COLUMNS_KEY = "active_annotation_columns" +CANONICAL_ANNOTATION_COLUMNS = ( + "cell_type_tsv", + "cell_type_ingest", + "cell_type_external", + "cell_type_refined", +) +_PLACEHOLDER_LABELS = { + "", "", "n/a", "na", "nan", "none", "null", "unannotated", +} + + +def non_placeholder_annotation_mask(labels): + """Return a boolean mask selecting real labels from a pandas Series.""" + normalized = labels.astype("string").str.strip().str.casefold() + return labels.notna() & ~normalized.isin(_PLACEHOLDER_LABELS) + + +def normalize_annotation_labels(labels, placeholder="Unannotated"): + """Strip real labels and collapse missing/placeholder values for display.""" + normalized = labels.astype("string").str.strip() + return normalized.where(non_placeholder_annotation_mask(labels), placeholder) + + +def meaningful_annotation_labels(adata, column): + """Return distinct, non-placeholder labels from an ``obs`` column.""" + if not column or column not in adata.obs.columns: + return [] + labels = adata.obs[column] + return ( + labels[non_placeholder_annotation_mask(labels)] + .astype(str) + .str.strip() + .unique() + .tolist() + ) + + +def has_meaningful_annotation(adata, column, min_labels=1): + """Whether ``column`` contains at least ``min_labels`` real labels.""" + return len(meaningful_annotation_labels(adata, column)) >= min_labels + + +def record_active_annotation_columns(adata, candidates=CANONICAL_ANNOTATION_COLUMNS): + """Store and return canonical annotation columns carrying at least one real label. + + Placeholder labels are ignored for this availability decision only. Thus, a + column containing both ``Unannotated`` and real cell types remains active and + retains its unannotated cells in downstream plots. + """ + active = [column for column in candidates if has_meaningful_annotation(adata, column)] + adata.uns[ACTIVE_ANNOTATION_COLUMNS_KEY] = active + return active + + +def active_annotation_columns(adata, candidates=CANONICAL_ANNOTATION_COLUMNS): + """Return usable annotation columns, with a fallback for older h5ad files. + + New outputs record the active columns explicitly. Objects written before that + metadata existed are inspected semantically, so an all-``Unannotated`` column + never becomes a report panel merely because it is present in ``obs``. + """ + candidates = list(candidates) + if ACTIVE_ANNOTATION_COLUMNS_KEY not in adata.uns: + return [ + column + for column in candidates + if has_meaningful_annotation(adata, column) + ] + + recorded = adata.uns[ACTIVE_ANNOTATION_COLUMNS_KEY] + if recorded is None: + recorded = [] + elif isinstance(recorded, str): + recorded = [recorded] + else: + recorded = [str(column) for column in recorded] + return [ + column + for column in candidates + if column in recorded and has_meaningful_annotation(adata, column) + ] diff --git a/workflow/scripts/integrate_samples.py b/workflow/scripts/integrate_samples.py index d2d6e97..36e93db 100644 --- a/workflow/scripts/integrate_samples.py +++ b/workflow/scripts/integrate_samples.py @@ -32,6 +32,10 @@ composition_pair, find_niche_column, composition_grouped, build_niche_palette, ) +from annotation_utils import ( + has_meaningful_annotation, + record_active_annotation_columns, +) try: import geosketch as sketch @@ -73,16 +77,6 @@ INTEGRATE_KEY = str(getattr(snakemake.params, "integrate_key", "sample") or "sample") -def _has_real_annotation(ad, col): - """True when an annotation column carries at least one meaningful label. External-only - runs leave cell_type_tsv all-'Unannotated', so colouring by it yields a single grey, - meaningless UMAP/barplot — skip those rather than emit them.""" - if not col or col not in ad.obs.columns: - return False - vals = set(ad.obs[col].astype(str).str.strip().unique()) - return len(vals - {"Unannotated", "", "nan", "NA", "None"}) > 0 - - def _apply_design_palette(ad, col): """Colour an obs design column from SAMPLE_COLORS. Values absent from a CONFIGURED palette fall back to grey (#cccccc); when NO palette is configured for the column, @@ -144,6 +138,8 @@ def _niche_umap(ad, title, out_path): log.info("Concatenating …") adata = sc.concat(adatas, join="inner", label="sample_batch", keys=sample_ids) log.info(" Concatenated: %d cells, %d genes", adata.n_obs, adata.n_vars) + active_annotations = record_active_annotation_columns(adata) + log.info(" Active annotation columns: %s", active_annotations or "none") # Carry over raw_counts if available if "raw_counts" in adatas[0].layers: @@ -168,8 +164,7 @@ def _niche_umap(ad, title, out_path): # Apply custom palettes if configured if isinstance(ANNOTATION_COLORS, dict): - for obs_key in ["sample_batch", "cell_type_tsv", - "cell_type_ingest", "cell_type_external", "leiden_uncorrected"]: + for obs_key in ["sample_batch", *active_annotations, "leiden_uncorrected"]: cd = ANNOTATION_COLORS.get(obs_key, {}) if cd and obs_key in adata.obs.columns: cats = adata.obs[obs_key].cat.categories if hasattr(adata.obs[obs_key], "cat") else [] @@ -183,7 +178,7 @@ def _niche_umap(ad, title, out_path): dpi=DPI, bbox_inches="tight") plt.close() - if _has_real_annotation(adata, "cell_type_tsv"): + if has_meaningful_annotation(adata, "cell_type_tsv"): sc.pl.umap(adata, color=["cell_type_tsv"], size=UMAP_POINT_SIZE, frameon=False, title="Uncorrected – by cell type (TSV)") plt.savefig(os.path.join(output_dir, "UMAP_uncorrected_by_celltype.png"), @@ -240,7 +235,7 @@ def _niche_umap(ad, title, out_path): dpi=DPI, bbox_inches="tight") plt.close() - if _has_real_annotation(adata_harmony, "cell_type_tsv"): + if has_meaningful_annotation(adata_harmony, "cell_type_tsv"): sc.pl.umap(adata_harmony, color=["cell_type_tsv"], size=UMAP_POINT_SIZE, frameon=False, title="Harmony – by cell type (TSV)") plt.savefig(os.path.join(output_dir, "UMAP_harmony_by_celltype.png"), @@ -266,9 +261,7 @@ def _niche_umap(ad, title, out_path): log.info("Composition barplots …") bar_dir = os.path.join(output_dir, "barplots") os.makedirs(bar_dir, exist_ok=True) - annot_cols = [c for c in ["cell_type_tsv", "cell_type_ingest", - "cell_type_external"] - if _has_real_annotation(adata, c)] + annot_cols = active_annotations has_regions = ("region_annotation" in adata.obs.columns and adata.obs["region_annotation"].nunique() > 1 and not all(adata.obs["region_annotation"] == "Unlabeled")) diff --git a/workflow/scripts/preprocess_umap.py b/workflow/scripts/preprocess_umap.py index 2a579cb..7934403 100644 --- a/workflow/scripts/preprocess_umap.py +++ b/workflow/scripts/preprocess_umap.py @@ -25,10 +25,11 @@ try: _here = os.path.dirname(os.path.abspath(__file__)) -except NameError: # very old Snakemake +except NameError: # very old Snakemake _here = os.getcwd() sys.path.insert(0, _here) -from external_metadata import annotated_cell_ids, optional_input_path, read_cell_metadata +from annotation_utils import non_placeholder_annotation_mask +from external_metadata import optional_input_path, read_cell_metadata # ── Logging ────────────────────────────────────────────────────────────────── log_handlers = [logging.StreamHandler(sys.stderr)] @@ -71,7 +72,7 @@ KEEP_UNANNOTATED = bool(getattr(snakemake.params, "keep_unannotated", True)) _ext_meta_in = getattr(snakemake.input, "external_meta", None) _precomputed_meta_in = getattr(snakemake.input, "precomputed_meta", None) -EXTERNAL_MASK = bool(EXTERNAL_ENABLED and not KEEP_UNANNOTATED and _ext_meta_in) +EXTERNAL_MASK = bool(EXTERNAL_ENABLED and not KEEP_UNANNOTATED) out_adata = str(snakemake.output.adata) out_metadata = str(snakemake.output.metadata) @@ -200,21 +201,38 @@ def _range(series): # and SKIP the pipeline QC thresholds (the user's external labels ARE the QC # decision). Genes expressed in no kept cell are dropped (data hygiene, not a # cell-QC threshold, so PCA/normalisation stay well-behaved). - ext_path = optional_input_path(_ext_meta_in) - _saved = read_cell_metadata(ext_path, required_columns=(EXTERNAL_COLUMN,)) - _annot_bc = annotated_cell_ids(_saved, EXTERNAL_COLUMN) + if _ext_meta_in: + ext_path = optional_input_path(_ext_meta_in) + _saved = read_cell_metadata( + ext_path, required_columns=(EXTERNAL_COLUMN,) + ) + _lab = _saved[EXTERNAL_COLUMN] + external_source = ext_path + elif EXTERNAL_COLUMN in adata.obs.columns: + _lab = adata.obs[EXTERNAL_COLUMN] + external_source = f"contract h5ad obs['{EXTERNAL_COLUMN}']" + else: + raise ValueError( + f"external_annotation column '{EXTERNAL_COLUMN}' is absent from " + "the contract h5ad obs and no external metadata TSV was supplied" + ) + _annot_bc = set( + _lab[non_placeholder_annotation_mask(_lab)].index.astype(str) + ) _keep = adata.obs_names.astype(str).isin(_annot_bc) adata = adata[_keep].copy() sc.pp.filter_genes(adata, min_cells=1) THRESHOLDS_SOURCE = "external_annotation (keep_unannotated=false; pipeline QC skipped)" MIN_COUNTS = MAX_COUNTS = MIN_GENES = MIN_CELLS = MAX_PCT_MT = None - log.info("External-driven cell set for %s: kept %d / %d cells (%d externally " - "annotated); pipeline QC thresholds SKIPPED.", - sample_id, adata.n_obs, n_cells_before, len(_annot_bc)) + log.info("External-driven cell set for %s from %s: kept %d / %d cells " + "(%d externally annotated); pipeline QC thresholds SKIPPED.", + sample_id, external_source, adata.n_obs, n_cells_before, + len(_annot_bc)) if adata.n_obs == 0: raise ValueError( f"No externally-annotated cells matched the contract barcodes for " - f"'{sample_id}'. Check that {ext_path} barcodes match adata.obs_names.") + f"'{sample_id}'. Check that {external_source} contains real labels " + "whose barcodes match adata.obs_names.") else: log.info("QC (source: %s): min_counts=%s min_cells=%s min_genes=%s " "max_counts=%s max_pct_mt=%s", THRESHOLDS_SOURCE, diff --git a/workflow/scripts/sample_report.py b/workflow/scripts/sample_report.py index d12fabe..d8ec6ff 100644 --- a/workflow/scripts/sample_report.py +++ b/workflow/scripts/sample_report.py @@ -32,6 +32,7 @@ sys.path.insert(0, _here) from composition_barplots import (draw_stacked_composition, find_niche_column, build_niche_palette) +from annotation_utils import active_annotation_columns, has_meaningful_annotation # ── Logging ────────────────────────────────────────────────────────────────── @@ -109,9 +110,10 @@ def _pick_keys(adata, niche_column=None): if manual not in adata.obs.columns: leiden_cols = sorted(c for c in adata.obs.columns if c.startswith("leiden_")) manual = leiden_cols[0] if leiden_cols else None - tsv = "cell_type_tsv" if "cell_type_tsv" in adata.obs.columns else None + active = set(active_annotation_columns(adata)) + tsv = "cell_type_tsv" if "cell_type_tsv" in active else None auto = next((c for c in ["cell_type_ingest", "cell_type_external"] - if c in adata.obs.columns), None) + if c in active), None) niche = find_niche_column(adata, niche_column) return manual, tsv, auto, niche @@ -179,8 +181,9 @@ def _build_page2(adata, sample_id, keys, annotation_colors, sample_col, bar and an accurately-sized dot-size legend. """ manual, tsv, auto, niche = keys - cmap = (annotation_colors.get(tsv, {}) - if (tsv and isinstance(annotation_colors, dict)) else {}) + composition_key = tsv or auto + cmap = (annotation_colors.get(composition_key, {}) + if (composition_key and isinstance(annotation_colors, dict)) else {}) figs = [] # Group the composition by region/niche when available; only fall back to @@ -193,7 +196,7 @@ def _build_page2(adata, sample_id, keys, annotation_colors, sample_col, if not bases: bases.append((sample_col, "Sample")) - if tsv: + if composition_key: try: ncol = 2 * len(bases) figb = plt.figure(figsize=(6.5 * ncol, 5.5)) @@ -201,7 +204,7 @@ def _build_page2(adata, sample_id, keys, annotation_colors, sample_col, fontsize=16, fontweight="bold", y=1.02) gs = gridspec.GridSpec(1, ncol, figure=figb, wspace=0.5) for i, (gkey, glabel) in enumerate(bases): - ct = pd.crosstab(adata.obs[gkey], adata.obs[tsv]) + ct = pd.crosstab(adata.obs[gkey], adata.obs[composition_key]) ax_abs = figb.add_subplot(gs[0, 2 * i]) draw_stacked_composition(ax_abs, ct, cmap, normalize=False, ylabel="Number of cells", xlabel=glabel, @@ -260,9 +263,7 @@ def _build_page2(adata, sample_id, keys, annotation_colors, sample_col, "cell_type_ingest": "cell type (ingest)"} marker_key = next( (k for k in (tsv, auto) - if k and k in adata.obs.columns - and len(set(adata.obs[k].astype(str).str.strip()) - - {"Unannotated", "", "nan", "NA", "None"}) >= 2), + if k and has_meaningful_annotation(adata, k, min_labels=2)), None) genes = _top_markers_by(adata, marker_key, n=10) if marker_key else [] if genes and marker_key: @@ -332,8 +333,7 @@ def _build_page2(adata, sample_id, keys, annotation_colors, sample_col, library_id = list(adata.uns["spatial"].keys())[0] # palettes - for obs_key in ["cell_type_tsv", "cell_type_ingest", - "cell_type_external"]: + for obs_key in active_annotation_columns(adata): apply_palette(adata, obs_key, ANNOTATION_COLORS) apply_region_palette(adata, REGION_COLORS) diff --git a/workflow/scripts/validate_input.py b/workflow/scripts/validate_input.py index acdd6b9..0e4f457 100644 --- a/workflow/scripts/validate_input.py +++ b/workflow/scripts/validate_input.py @@ -39,10 +39,11 @@ try: _here = os.path.dirname(os.path.abspath(__file__)) -except NameError: # very old Snakemake +except NameError: # very old Snakemake _here = os.getcwd() sys.path.insert(0, _here) -from external_metadata import annotated_cell_ids, optional_input_path, read_cell_metadata +from annotation_utils import non_placeholder_annotation_mask +from external_metadata import optional_input_path, read_cell_metadata # ── Logging ────────────────────────────────────────────────────────────────── log_handlers = [logging.StreamHandler(sys.stderr)] @@ -195,25 +196,46 @@ def _valid_image_entry(entry): ext_overlap = None if ext_enabled and ext_column: try: - _sv = read_cell_metadata(ext_path, required_columns=(ext_column,)) - _bc = annotated_cell_ids(_sv, ext_column) - _match = len(_bc & set(adata.obs_names.astype(str))) - frac = round(_match / len(_bc), 4) if _bc else 0.0 - ext_overlap = { - "annotated_barcodes": len(_bc), - "matched_obs_names": _match, - "fraction_matched": frac, - } - if _bc and _match == 0: + if ext_path: + _sv = read_cell_metadata( + ext_path, required_columns=(ext_column,) + ) + _lab = _sv[ext_column] + elif ext_column in adata.obs.columns: + _lab = adata.obs[ext_column] + else: errors.append( - f"external_annotation: 0 of {len(_bc)} annotated " - "barcodes match obs_names — check barcode formatting" + f"external_annotation column '{ext_column}' is absent from the " + "contract h5ad obs and no external metadata TSV was supplied" ) - elif _bc and frac < 0.5: - warnings.append( - f"external_annotation: only {_match}/{len(_bc)} annotated " - f"barcodes match obs_names ({frac:.0%})" + _lab = None + + if _lab is not None: + _bc = set( + _lab[non_placeholder_annotation_mask(_lab)].index.astype(str) ) + _match = len(_bc & set(adata.obs_names.astype(str))) + frac = round(_match / len(_bc), 4) if _bc else 0.0 + ext_overlap = { + "annotated_barcodes": len(_bc), + "matched_obs_names": _match, + "fraction_matched": frac, + } + if _bc and _match == 0: + errors.append( + f"external_annotation: 0 of {len(_bc)} annotated barcodes " + "match obs_names — check barcode formatting" + ) + elif _bc and frac < 0.5: + warnings.append( + f"external_annotation: only {_match}/{len(_bc)} annotated " + f"barcodes match obs_names ({frac:.0%})" + ) + elif not _bc: + warnings.append( + "external_annotation column contains no non-placeholder " + "labels; external-annotation plots will be omitted" + ) except Exception as e: errors.append(f"external_annotation metadata is invalid: {e}")