PR: Backbone subsampling — stable sequence set across reruns

CLAUDE.md

@@ -60,6 +60,35 @@ for forward-compatibility but does not use it). Omitting it defaults to today wi
 runs with the same `--run-date` and the same config produce the same `config_hash` and `run_id`
 in the manifest. When not passed, `augur subsample` has no upper date bound.
 
+**`--backbone-from` semantics:** opt-in stable backbone — pass the workdir of a previous run to
+force-include its subsampled strain list in the new subsample.  The new subsample becomes the
+union of *(stable backbone strains) + (freshly-selected new sequences)*, making results comparable
+across runs (e.g. a 2-year-later rerun of an RSV build).
+
+```bash
+# Run A (initial build — e.g. June 2024)
+flexpipe-run \
+    --config   builds/rsv-a-brazil/config.yaml \
+    --workdir  /path/to/workdir/rsv-A \
+    --run-date 2024-06-01
+
+# Run B (2 years later — backbone anchors the 2024 selection)
+flexpipe-run \
+    --config        builds/rsv-a-brazil/config.yaml \
+    --workdir       /path/to/workdir/rsv-B \
+    --run-date      2026-06-01 \
+    --backbone-from /path/to/workdir/rsv-A
+```
+
+Key limitations of `--backbone-from`:
+- **Best-effort retention:** strains dropped by upstream QC (`augur filter`) or
+  `clade_filter` cannot be force-kept — `include` only applies within `augur subsample`.
+  Backbone retention is bounded by the quality contract.
+- **Runtime-only:** never put the backbone path in `builds/<name>/config.yaml` (it bakes a
+  machine-specific absolute path into version control).  Always pass it via the CLI flag.
+- **Self-reference guard:** pointing `--backbone-from` at the current workdir exits with code 2.
+- **Missing previous run:** a warning is logged and the run proceeds without a backbone.
+
 ### Workflow Control
 ```bash
 # Run ingest only
@@ -243,6 +272,13 @@ path, not the workdir-local resolved snapshot. The Snakefile reads
 - Reads `builds/<name>/subsample.yaml`
 - For YFV Brazil: subsamples by division (state) and year
 - Generates subsampled metadata and sequences in `<workdir>/results/subsampled/`
+- **Backbone support** (`--backbone-from`): when enabled, `resolve_subsample_config` injects a
+  synthetic `samples.__backbone__: {include: <workdir>/config/backbone_strains.txt}` entry so
+  `augur subsample` force-keeps the previous run's strains regardless of group caps.
+  The `backbone_strains.txt` file is written by `_materialize_backbone()` in `run.py` before
+  Snakemake is invoked; `SubsamplingConfig.backbone_strains` carries the path through the
+  resolved config to the Snakefile.  Feature is a complete no-op when `backbone_strains` is
+  `None` (the default).
 
 **Colors and coordinates** (can be run in parallel):
 - `flexpipe-name2hue`: deterministic hue assignment from subsampled metadata

flexpipe/config.py

@@ -306,6 +306,15 @@ def validate_rare_state_label(cls, value: str) -> str:
 class SubsamplingConfig(BaseModel):
     model_config = ConfigDict(extra="forbid")
     random_seed: int = 42
+    backbone_strains: str | None = None
+    """Absolute path to a one-strain-per-line include-list from a previous run.
+
+    Populated at runtime by ``flexpipe-run --backbone-from``; never authored in build
+    ``config.yaml`` (it would bake a machine-specific absolute path into version control).
+    When set, ``resolve_subsample_config`` injects a synthetic ``samples.__backbone__``
+    sample set so ``augur subsample`` force-keeps the listed strains regardless of group
+    caps.  ``None`` (the default) disables the feature entirely.
+    """
 
 
 class LineageColumnsConfig(BaseModel):
@@ -712,34 +721,49 @@ def resolve_subsample_config(
     raw: dict,
     run_date: str | None,
     subsample_config_path: str | Path | None = None,
+    backbone_strains: str | None = None,
 ) -> dict:
-    """Return a copy of the subsample config dict with ``defaults.max_date`` injected.
+    """Return a copy of the subsample config dict with runtime overrides injected.
 
     When *run_date* is provided, it is written into the ``defaults`` section of the
     subsample config as the ``max_date`` upper bound for ``augur subsample``.  This
     ensures a scheduled rerun with ``--run-date 2026-01-01`` is bounded by that date
     rather than anchored to the system clock.
 
-    When *run_date* is empty or ``None`` the dict is returned unchanged — preserving
-    current behaviour for direct ``snakemake`` invocations that do not pass
-    ``--config run_date=``.
+    When *backbone_strains* is provided (an absolute path to a one-strain-per-line
+    include-list from a previous run), a synthetic sample set ``__backbone__`` is
+    injected.  ``augur subsample`` force-keeps every listed strain that is still
+    present in the input regardless of ``group_by`` / ``sequences_per_group`` caps,
+    producing a stable union of (previous selection) + (newly-subsampled sequences).
+
+    Both overrides default to ``None`` / ``""`` — i.e. no-op — preserving current
+    behaviour for direct ``snakemake`` invocations that do not pass the extra flags.
 
     Args:
         raw: Parsed subsample config dict (e.g. from ``builds/<name>/subsample.yaml``).
         run_date: Reference date in ``YYYY-MM-DD`` format, or ``None`` / ``""`` to skip.
+        subsample_config_path: Path to the source YAML; used to resolve relative
+            include/exclude paths.  ``None`` skips path resolution.
+        backbone_strains: Absolute path to the backbone include-list written by
+            ``flexpipe-run --backbone-from``.  ``None`` disables backbone injection.
 
     Returns:
-        A shallow copy of *raw* with the ``defaults`` section updated.
+        A deep copy of *raw* with the injected overrides applied.
     """
     out = (
         resolve_subsample_paths(raw, subsample_config_path)
         if subsample_config_path is not None
         else copy.deepcopy(raw)
     )
-    if not run_date:
-        return out
-    out.setdefault("defaults", {})["max_date"] = run_date
-    logger.debug("resolve_subsample_config: set defaults.max_date=%s", run_date)
+    if run_date:
+        out.setdefault("defaults", {})["max_date"] = run_date
+        logger.debug("resolve_subsample_config: set defaults.max_date=%s", run_date)
+    if backbone_strains:
+        out.setdefault("samples", {})["__backbone__"] = {"include": backbone_strains}
+        logger.debug(
+            "resolve_subsample_config: injected __backbone__ include-list from %s",
+            backbone_strains,
+        )
     return out
 
 

flexpipe/paths.py

@@ -40,6 +40,7 @@
         latlongs.tsv              # generated by coordinates step
         colour_scheme.tsv         # generated by colours step
         name2hue.tsv              # generated by generate_name2hue step
+        backbone_strains.txt      # optional; written by flexpipe-run --backbone-from
       logs/
         ingest.log
         phylo.log
@@ -129,6 +130,17 @@ def subsample_config_resolved(self) -> Path:
         """Workdir-local subsample config with ``defaults.max_date`` injected from run_date."""
         return self.generated_config_dir / "subsample_resolved.yaml"
 
+    @property
+    def backbone_strains(self) -> Path:
+        """One-strain-per-line include-list materialized from a previous run's subsampled metadata.
+
+        Written by ``flexpipe-run --backbone-from`` before Snakemake is invoked; referenced as a
+        synthetic ``samples.__backbone__.include`` entry in the resolved subsample config so
+        ``augur subsample`` force-keeps the listed strains regardless of group caps.
+        The parent ``config/`` directory is created by :meth:`ensure_dirs`.
+        """
+        return self.generated_config_dir / "backbone_strains.txt"
+
     # ── mutable cache ────────────────────────────────────────────────────────
     @property
     def cache_dir(self) -> Path:

flexpipe/run.py

@@ -17,7 +17,8 @@
         --workdir /data/runs/yfv-brazil/2025-06-01 \\
         --run-date 2025-06-01 \\
         [--stage ingest|phylo|all] \\
-        [--cores 4]
+        [--cores 4] \\
+        [--backbone-from /data/runs/yfv-brazil/2023-06-01]
 
 Exit codes:
     0 — success
@@ -26,6 +27,7 @@
 """
 
 import argparse
+import csv
 import logging
 import subprocess
 import sys
@@ -97,6 +99,88 @@ def _seed_coordinate_cache_with_shared(
     )
 
 
+def _materialize_backbone(
+    backbone_from: Path | None,
+    paths: WorkdirPaths,
+) -> Path | None:
+    """Extract the previous run's subsampled strain list and write it to the workdir.
+
+    Reads ``<backbone_from>/results/subsampled/metadata.tsv``, pulls the ``strain``
+    column, and writes one strain per line to ``paths.backbone_strains``.  Returns the
+    output path so the caller can set ``cfg.subsampling.backbone_strains``.
+
+    Returns ``None`` (no-op) when:
+
+    * *backbone_from* is ``None`` (feature disabled).
+    * *backbone_from* resolves to the current workdir (self-reference guard — exits 2).
+    * The previous metadata file is missing (warns, continues without backbone).
+    * The previous metadata has no strains (warns, continues without backbone).
+
+    Args:
+        backbone_from: Path to the previous run's workdir, or ``None``.
+        paths: :class:`~flexpipe.paths.WorkdirPaths` for the current run.
+
+    Returns:
+        Absolute :class:`~pathlib.Path` to the written include-list, or ``None``.
+    """
+    if backbone_from is None:
+        return None
+
+    prev_root = Path(backbone_from).resolve()
+    if prev_root == paths.root:
+        raise SystemExit(
+            f"--backbone-from points at the current workdir ({paths.root}).\n"
+            "Pass a *previous* run's workdir, not the one being built."
+        )
+
+    prev_metadata = WorkdirPaths.from_root(prev_root).subsampled_metadata
+    if not prev_metadata.exists():
+        logger.warning(
+            "backbone: no subsampled metadata found at %s — proceeding without backbone.",
+            prev_metadata,
+        )
+        return None
+
+    strains: list[str] = []
+    try:
+        with open(prev_metadata, newline="", encoding="utf-8") as fh:
+            reader = csv.DictReader(fh, delimiter="\t")
+            if reader.fieldnames is None or "strain" not in reader.fieldnames:
+                logger.warning(
+                    "backbone: previous metadata at %s has no 'strain' column — "
+                    "proceeding without backbone.",
+                    prev_metadata,
+                )
+                return None
+            for row in reader:
+                s = row.get("strain", "").strip()
+                if s:
+                    strains.append(s)
+    except Exception as exc:
+        logger.warning(
+            "backbone: could not read %s (%s) — proceeding without backbone.", prev_metadata, exc
+        )
+        return None
+
+    if not strains:
+        logger.warning(
+            "backbone: previous metadata at %s contains no strains — proceeding without backbone.",
+            prev_metadata,
+        )
+        return None
+
+    out = paths.backbone_strains
+    out.parent.mkdir(parents=True, exist_ok=True)
+    out.write_text("\n".join(strains) + "\n", encoding="utf-8")
+    logger.info(
+        "backbone: materialized %d strains from %s → %s",
+        len(strains),
+        prev_metadata,
+        out,
+    )
+    return out
+
+
 def _run_snakemake(
     snakefile: Path,
     config_path: Path,
@@ -149,6 +233,7 @@ def run_pipeline(
     run_date: str,
     stage: str = "all",
     cores: int = 4,
+    backbone_from: Path | None = None,
 ) -> int:
     """Run the pipeline for one build.
 
@@ -158,6 +243,10 @@ def run_pipeline(
         run_date: Reference date for this run (``YYYY-MM-DD``).
         stage: ``"ingest"``, ``"phylo"``, or ``"all"`` (default).
         cores: Number of CPU cores to pass to Snakemake.
+        backbone_from: Path to a previous run's workdir.  When provided, the previous
+            run's subsampled strain list is force-included in the new subsample so the
+            sequence SET stays stable across reruns.  ``None`` (default) disables the
+            feature — behaviour is identical to the current pipeline.
 
     Returns:
         Exit code (0 = success).
@@ -207,6 +296,7 @@ def run_pipeline(
             run_date=run_date,
             stage=stage,
             cores=cores,
+            backbone_from=backbone_from,
         )
     finally:
         lock.release()
@@ -221,8 +311,15 @@ def _run_pipeline_locked(
     run_date: str,
     stage: str,
     cores: int,
+    backbone_from: Path | None = None,
 ) -> int:
     """Inner pipeline body — called inside the workdir lock."""
+    # Materialize the backbone strain list before writing the resolved config so
+    # cfg.subsampling.backbone_strains propagates into the single Snakemake --configfile.
+    backbone_path = _materialize_backbone(backbone_from, paths)
+    if backbone_path is not None:
+        cfg.subsampling.backbone_strains = str(backbone_path)
+
     snakemake_overrides = write_snakemake_config_overrides(
         cfg, paths.snakemake_config_overrides, config_path
     )
@@ -238,6 +335,9 @@ def _run_pipeline_locked(
 
     manifest = Manifest(run_date=run_date, build_name=build_name, config_path=config_path)
     manifest.record_provenance(cfg, snakemake_overrides)
+    if backbone_path is not None:
+        manifest.record("backbone_from", str(backbone_from))
+        manifest.record("backbone_strain_count", len(backbone_path.read_text().splitlines()))
 
     min_sequences = cfg.qc.min_sequences
 
@@ -329,6 +429,20 @@ def main() -> None:
         default=4,
         help="Number of CPU cores for Snakemake (default: 4)",
     )
+    parser.add_argument(
+        "--backbone-from",
+        default=None,
+        type=Path,
+        metavar="PREV_WORKDIR",
+        help=(
+            "Path to a previous run's workdir.  When provided, the strain list from that "
+            "run's results/subsampled/metadata.tsv is force-included in the new subsample "
+            "via augur subsample's per-sample include mechanism.  The new subsample becomes "
+            "the union of (stable backbone strains) + (freshly-selected new sequences), "
+            "making results comparable across runs.  Omit for a fully-fresh subsample "
+            "(default)."
+        ),
+    )
     parser.add_argument(
         "--log-level",
         default="INFO",
@@ -356,5 +470,6 @@ def main() -> None:
         run_date=run_date,
         stage=args.stage,
         cores=args.cores,
+        backbone_from=args.backbone_from.resolve() if args.backbone_from else None,
     )
     sys.exit(rc)

ingest/Snakefile

@@ -421,27 +421,37 @@ rule qc_summary:
 # ─── Subsampling ──────────────────────────────────────────────────────────────
 
 rule resolve_subsample_config:
-    """Write a workdir-local copy of the subsample config with defaults.max_date injected.
+    """Write a workdir-local copy of the subsample config with runtime overrides injected.
 
     When flexpipe-run passes --config run_date=<YYYY-MM-DD>, this rule injects that date
     as the upper bound for augur subsample so the analysis window is reproducibly anchored
     to the declared reference date rather than the system clock.
 
+    When flexpipe-run passes --backbone-from <prev_workdir>, a synthetic sample set
+    __backbone__ is injected with an include path pointing at the materialized strain list
+    so augur subsample force-keeps those strains regardless of group caps.
+
     When run_date is empty (direct snakemake invocation), the source subsample config is
     copied unchanged — preserving current behaviour.
     """
     input:
         subsample_cfg = files.subsample_config
     params:
-        run_date = _run_date
+        run_date        = _run_date,
+        backbone_strains = _sub.get("backbone_strains") or ""
     output:
         resolved = f"{_wd}/config/subsample_resolved.yaml"
     run:
         import yaml
         from flexpipe.config import resolve_subsample_config
         with open(input.subsample_cfg, encoding="utf-8") as fh:
             raw = yaml.safe_load(fh) or {}
-        resolved = resolve_subsample_config(raw, params.run_date, input.subsample_cfg)
+        resolved = resolve_subsample_config(
+            raw,
+            params.run_date,
+            input.subsample_cfg,
+            backbone_strains=params.backbone_strains or None,
+        )
         import os; os.makedirs(os.path.dirname(output.resolved), exist_ok=True)
         with open(output.resolved, "w", encoding="utf-8") as fh:
             yaml.safe_dump(resolved, fh, default_flow_style=False, sort_keys=False)