diff --git a/.github/workflows/build_container.yml b/.github/workflows/build_container.yml
index 1b97a6a..7c81d4a 100644
--- a/.github/workflows/build_container.yml
+++ b/.github/workflows/build_container.yml
@@ -1,13 +1,16 @@
-name: Build R container
+name: Build R containers
 
-# The image tag is read from conf/containers.config, so that single line is both what the
-# pipeline pulls and what this workflow publishes -- the two cannot drift.
-# Runs on pushes to a feature branch, whenever the container changes.
+# The image tag of the container is read from conf/containers.config, so that the tag is
+# what the pipeline pulls and what this workflow publishes. Each container is built independently,
+# so updating one tag in conf/containers.config rebuilds only that image:
+#   feature branch - rebuilds its tag on every push, unless main or devel is pinned to that tag
+#   main or devel  - builds only a tag that has never been published, so most merges are a no-op
+#   manual run     - 'force' rebuilds and overwrites the tag on any branch
 
 on:
   push:
     paths:
-      - 'containers/r/**'
+      - 'containers/**'
       - 'conf/containers.config'
       - '.github/workflows/build_container.yml'
   workflow_dispatch:
@@ -22,81 +25,130 @@ concurrency:
   cancel-in-progress: true
 
 jobs:
-  build:
-    name: Build and publish
+  discover:
+    name: Select containers to build
     runs-on: ubuntu-latest
     permissions:
       contents: read
-      packages: write
+      packages: read
+    outputs:
+      containers: ${{ steps.select.outputs.containers }}
 
     steps:
       - uses: actions/checkout@v4
 
-      - name: Resolve image from conf/containers.config
-        id: image
+      - uses: docker/login-action@v3
+        with:
+          registry: ghcr.io
+          username: ${{ github.actor }}
+          password: ${{ secrets.GITHUB_TOKEN }}
+
+      - name: Select containers to build
+        id: select
         run: |
-          IMAGE=$(sed -n "s|.*withLabel: *'r'.*container *= *\"\([^\"]*\)\".*|\1|p" conf/containers.config)
-          if [ -z "$IMAGE" ]; then
-            echo "::error::could not parse the 'r' container from conf/containers.config"
-            exit 1
+          RELEASED=false
+          case "${{ github.ref_name }}" in main|devel) RELEASED=true ;; esac
+
+          # Each base branch is fetched into its own ref so both stay readable in the loop
+          BASES=""
+          for BASE in devel main; do
+            if git fetch --no-tags --depth=1 origin "$BASE:refs/base/$BASE" 2>/dev/null; then
+              BASES="$BASES refs/base/$BASE"
+            fi
+          done
+
+          # devel is at or ahead of main, so it is the baseline for "did this branch edit a container"
+          DEVEL=refs/base/devel
+
+          # a container is only rebuilt by a push that touches it, so the previous push is the
+          # baseline for that. devel is the fallback when that commit cannot be fetched, which is
+          # the case on the first push of a branch, on a force push and on a manual run
+          PUSH_BASE=$DEVEL
+          BEFORE="${{ github.event.before }}"
+          if [ -n "$BEFORE" ] && git fetch --no-tags --depth=1 origin "$BEFORE" 2>/dev/null; then
+            PUSH_BASE=$BEFORE
           fi
-          case "$IMAGE" in
-            ghcr.io/goekelab/*) ;;
-            *) echo "::error::refusing to push outside ghcr.io/goekelab: $IMAGE"; exit 1 ;;
-          esac
-          echo "image=$IMAGE" >> "$GITHUB_OUTPUT"
-          echo "Resolved container: $IMAGE"
-
-      - name: Refuse to overwrite a tag main or devel is pinned to
-        if: github.ref_name != 'main' && github.ref_name != 'devel'
-        run: |
-          IMAGE="${{ steps.image.outputs.image }}"
-          for BASE in main devel; do
-            git fetch --no-tags --depth=1 origin "$BASE" 2>/dev/null || continue
-            PINNED=$(git show "FETCH_HEAD:conf/containers.config" 2>/dev/null \
-              | sed -n "s|.*withLabel: *'r'.*container *= *\"\([^\"]*\)\".*|\1|p")
-            if [ "$IMAGE" = "$PINNED" ]; then
-              echo "::error::$BASE is pinned to $IMAGE -- bump the tag in conf/containers.config before rebuilding it"
+
+          # Every containers/<label>/Dockerfile is an image whose <label> is its process label
+          : > selected.txt
+          for DOCKERFILE in containers/*/Dockerfile; do
+            LABEL=$(basename "$(dirname "$DOCKERFILE")")
+            IMAGE=$(sed -n "s|.*withLabel: *'$LABEL'.*container *= *\"\([^\"]*\)\".*|\1|p" conf/containers.config)
+
+            if [ -z "$IMAGE" ]; then
+              echo "::error::could not parse the '$LABEL' container from conf/containers.config"
               exit 1
             fi
+            case "$IMAGE" in
+              ghcr.io/goekelab/*) ;;
+              *) echo "::error::refusing to push outside ghcr.io/goekelab: $IMAGE"; exit 1 ;;
+            esac
+
+            EDITED_ON_BRANCH=false
+            if ! git diff --quiet "$DEVEL" HEAD -- "containers/$LABEL"; then
+              EDITED_ON_BRANCH=true
+            fi
+
+            EDITED_BY_PUSH=false
+            if ! git diff --quiet "$PUSH_BASE" HEAD -- "containers/$LABEL"; then
+              EDITED_BY_PUSH=true
+            fi
+
+            # A tag main or devel is pinned to is immutable, so an edit needs a tag bump first
+            if [ "$RELEASED" = false ] && [ "$EDITED_ON_BRANCH" = true ]; then
+              for REF in $BASES; do
+                PINNED=$(git show "$REF:conf/containers.config" 2>/dev/null \
+                  | sed -n "s|.*withLabel: *'$LABEL'.*container *= *\"\([^\"]*\)\".*|\1|p")
+                if [ "$IMAGE" = "$PINNED" ]; then
+                  echo "::error::${REF#refs/base/} is pinned to $IMAGE -- bump the tag in conf/containers.config before rebuilding it"
+                  exit 1
+                fi
+              done
+            fi
+
+            # An unpublished tag is always built, since the pipeline cannot pull what is not there
+            if [ "${{ inputs.force }}" != "true" ] && docker manifest inspect "$IMAGE" >/dev/null 2>&1; then
+              if [ "$RELEASED" = true ] || [ "$EDITED_BY_PUSH" = false ]; then
+                echo "::notice::$IMAGE is already published and unchanged -- skipping $LABEL"
+                continue
+              fi
+            fi
+
+            echo "$LABEL $IMAGE" >> selected.txt
           done
 
+          CONTAINERS=$(jq -R -s -c 'split("\n") | map(select(length > 0) | split(" ")) | map({label: .[0], image: .[1]})' < selected.txt)
+          echo "containers=$CONTAINERS" >> "$GITHUB_OUTPUT"
+          echo "Building: $CONTAINERS"
+
+  build:
+    name: Build and publish (${{ matrix.container.label }})
+    needs: discover
+    if: needs.discover.outputs.containers != '[]'
+    runs-on: ubuntu-latest
+    permissions:
+      contents: read
+      packages: write
+    strategy:
+      fail-fast: false
+      matrix:
+        container: ${{ fromJson(needs.discover.outputs.containers) }}
+
+    steps:
+      - uses: actions/checkout@v4
+
       - uses: docker/login-action@v3
         with:
           registry: ghcr.io
           username: ${{ github.actor }}
           password: ${{ secrets.GITHUB_TOKEN }}
 
-      - name: Decide whether to build
-        id: decide
-        run: |
-          IMAGE="${{ steps.image.outputs.image }}"
-          if docker manifest inspect "$IMAGE" >/dev/null 2>&1; then EXISTS=true; else EXISTS=false; fi
-
-          if [ "${{ inputs.force }}" = "true" ]; then
-            BUILD=true
-          elif [ "${{ github.ref_name }}" = "main" ] || [ "${{ github.ref_name }}" = "devel" ]; then
-            # Released tags are immutable: the merge that brought the bump in has already
-            # published it from the feature branch, so there is nothing left to do
-            [ "$EXISTS" = "false" ] && BUILD=true || BUILD=false
-          else
-            # Feature branch: the tag is unreleased (guarded above), so rebuild freely
-            BUILD=true
-          fi
-
-          echo "build=$BUILD" >> "$GITHUB_OUTPUT"
-          if [ "$BUILD" = "false" ]; then
-            echo "::notice::$IMAGE is already published -- skipping build"
-          fi
-
       - uses: docker/setup-buildx-action@v3
-        if: steps.decide.outputs.build == 'true'
 
       - uses: docker/build-push-action@v6
-        if: steps.decide.outputs.build == 'true'
         with:
-          context: containers/r
+          context: containers/${{ matrix.container.label }}
           push: true
-          tags: ${{ steps.image.outputs.image }}
-          cache-from: type=gha
-          cache-to: type=gha,mode=max
+          tags: ${{ matrix.container.image }}
+          cache-from: type=gha,scope=${{ matrix.container.label }}
+          cache-to: type=gha,mode=max,scope=${{ matrix.container.label }}
diff --git a/.github/workflows/smoke_test.yml b/.github/workflows/smoke_test.yml
index cfea629..27baa44 100644
--- a/.github/workflows/smoke_test.yml
+++ b/.github/workflows/smoke_test.yml
@@ -11,7 +11,7 @@ concurrency:
 
 jobs:
   smoke_test:
-    name: "${{ matrix.profile }} (NF ${{ matrix.nextflow_version }})${{ matrix.quantification_mode && format(' [{0}]', matrix.quantification_mode) || '' }}${{ matrix.bam_only && ' [bam_only]' || '' }}"
+    name: "${{ matrix.profile }} (NF ${{ matrix.nextflow_version }})${{ matrix.quantification_mode && format(' [{0}]', matrix.quantification_mode) || '' }}${{ matrix.bam_only && ' [bam_only]' || '' }}${{ matrix.extra_args && format(' [{0}]', matrix.extra_args) || '' }}"
     runs-on: ubuntu-latest
     strategy:
       matrix:
@@ -21,16 +21,24 @@ jobs:
           - {nextflow_version: "26.04.0", profile: test_sc_bam}
           - {nextflow_version: "26.04.0", profile: test_sc_multi}
           - {nextflow_version: "26.04.0", profile: test_visium}
+          - {nextflow_version: "26.04.0", profile: test_visium_hd}
           - {nextflow_version: "26.04.0", profile: test_custom}
-          # latest-stable - all profiles (default EM_clusters)
+          - {nextflow_version: "26.04.0", profile: test_sc_quant_data}
+          - {nextflow_version: "26.04.0", profile: test_visium_hd_quant_data}
+          # latest-stable - all profiles (default clusteredEM)
           - {nextflow_version: "latest-stable", profile: test_sc_fastq}
           - {nextflow_version: "latest-stable", profile: test_sc_bam}
           - {nextflow_version: "latest-stable", profile: test_sc_multi}
           - {nextflow_version: "latest-stable", profile: test_visium}
+          - {nextflow_version: "latest-stable", profile: test_visium_hd}
           - {nextflow_version: "latest-stable", profile: test_custom}
+          - {nextflow_version: "latest-stable", profile: test_sc_quant_data}
+          - {nextflow_version: "latest-stable", profile: test_visium_hd_quant_data}
           # latest-stable - extra parameter variants
           - {nextflow_version: "latest-stable", profile: test_sc_fastq, quantification_mode: EM}
           - {nextflow_version: "latest-stable", profile: test_sc_multi, quantification_mode: EM}
+          - {nextflow_version: "latest-stable", profile: test_visium_hd, quantification_mode: EM}
+          - {nextflow_version: "latest-stable", profile: test_visium_hd, extra_args: "--banksy false"}
           - {nextflow_version: "latest-stable", profile: test_sc_fastq, quantification_mode: no_quant}
           - {nextflow_version: "latest-stable", profile: test_sc_fastq, bam_only: true}
     steps:
@@ -42,4 +50,5 @@ jobs:
         run: |
           nextflow run . -profile test_base,${{ matrix.profile }},docker \
             ${{ matrix.quantification_mode && format('--quantification_mode {0}', matrix.quantification_mode) || '' }} \
-            ${{ matrix.bam_only && '--bam_only' || '' }}
+            ${{ matrix.bam_only && '--bam_only' || '' }} \
+            ${{ matrix.extra_args || '' }}
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7425217..4b471d0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,32 @@ This file contains all notable changes to Bambu-Pipe.
 
 ---
 
+## [v0.10.0] - 2026-08-05
+
+### Added
+- Visium HD workflow (`--visium_hd`), run as a single sample from a Spaceranger-aligned, barcode-tagged BAM
+  - Transcript discovery and read-to-transcript assignment at the native 2 µm resolution, with counts aggregated to every bin listed in `--bins`
+  - `--barcode_mappings` for the Spaceranger `barcode_mappings.parquet`, used to assign 2 µm spots to bins
+  - Out-of-tissue reads filtered from the BAM using the 2 µm `tissue_positions.parquet`
+  - Spatially aware clustering with Banksy (`--banksy`, `--banksy_lambda`, `--banksy_k_geom`), or gene expression alone
+  - `--clustering_bin` to select the resolution to cluster at; cluster labels are expanded back to 2 µm spots for quantification
+  - Spot-level quantification at every resolution under `--quantification_mode EM`
+  - `test_visium_hd` smoke test profile with synthetic example data
+- `--manual_clustering` to restart the pipeline from cluster assignments generated outside the pipeline, for both standard and Visium HD runs
+- `test_sc_quant_data` and `test_visium_hd_quant_data` smoke test profiles covering the manual clustering restart
+- Self-hosted `bambu` and `seurat` container images published to `ghcr.io/goekelab`, built by the `build_container.yml` GitHub Actions workflow
+- Shared R helpers in `bin/` for transcript discovery, Seurat object creation, count saving, and cluster mapping
+
+### Changed
+- Renamed `--resolution` to `--seurat_resolution`
+- Renamed the `--quantification_mode` option `EM_clusters` to `clusteredEM`, matching the `bambu.singlecell` API
+- `quant_data.rds` and `extended_annotations.rds` are now always published to `intermediate_R/`, so a manual clustering run can restart from them
+- Seurat objects are built from the published count directories and Bambu's `colData` instead of the `SummarizedExperiment`
+- `clusters.rds` is now a named vector of `id -> cluster` label, replacing the per-sample list of `CompressedCharacterList`
+- Cluster-level quantification moved into a single module shared by the standard and Visium HD workflows
+- Restructured modules into `standard/`, `visium_hd/`, and `shared/` directories
+- Smoke tests now run on pull request and manual dispatch only, with in-progress runs cancelled on a new push
+
 ## [v0.9.1] - 2026-05-20
 
 ### Added
diff --git a/README.md b/README.md
index 25e7d0a..1341e4a 100644
--- a/README.md
+++ b/README.md
@@ -4,8 +4,10 @@
 - [Installation](#installation)
 - [General Usage](#general-usage)
 - [Parameters](#parameters)
-- [Output](#output)
 - [Spatial Analysis](#spatial-analysis)
+  - [Visium Spatial Gene Expression](#visium-spatial-gene-expression)
+  - [Visium HD](#visium-hd)
+- [Output](#output)
 - [Advanced Usage](#advanced-usage)
 - [Additional Information](#additional-information)
 - [Release History](#release-history)
@@ -32,7 +34,7 @@ This pipeline performs context-aware transcript discovery and quantification fro
 5. Read class construction and transcript discovery with [Bambu](https://github.com/GoekeLab/bambu/tree/BambuDev) (performed jointly across all samples)
 6. (Optional) Transcript quantification with Bambu using one of two modes:
    - Cluster-level EM: Gene expression-based cell clustering with [Seurat](https://github.com/satijalab/seurat) across all samples, followed by per-sample cluster-level transcript quantification
-   - Single-cell EM: Per-cell transcript quantification 
+   - Single-cell EM: Transcript quantification for each cell, or for each spot in spatial samples
 
 ### **Installation** 
 Install the following dependencies before running the pipeline:
@@ -106,12 +108,16 @@ For the following chemistries, the pipeline handles the full workflow — FASTQ
 - `visium-v4` (Visium CytAssist Spatial Gene Expression Slide 6.5 mm; serial prefix V4)
 - `visium-v5` (Visium CytAssist Spatial Gene Expression Slide 11mm; serial prefix V5)
 
-> **Note:** Visium samples must be run one sample at a time. Multi-sample runs are not supported for Visium chemistries.
+> **Note:** Visium Spatial Gene Expression samples (`visium-v*`) must be run one at a time. Multi-sample runs are not supported for these chemistries.
 
-*Custom Chemistry*
+*Custom Chemistries*
 
 If your dataset uses a chemistry not listed above, or if you prefer to handle FASTQ preprocessing and genome alignment manually, provide a pre-processed, demultiplexed BAM file as input. See the [Advanced Usage](#advanced-usage) section for details.
 
+*Visium HD*
+
+The pipeline also supports Visium HD samples. See the [Visium HD](#visium-hd) section for more details.
+
 **Pipeline Configuration**
 
 *Nextflow Profiles*
@@ -131,7 +137,7 @@ To configure the executor and container, pass profile types via the `-profile` a
 **Mandatory**
 - `--input` [string]: Path to the samplesheet `.csv` file 
 - `--genome` [string]: Path to the reference genome `.fa`, `.fasta`, or `.fa.gz` file 
-- `--annotation` [string]: Path to the reference annotation `.gtf`, `.gff`, `.gtf.gz`, or `.gff.gz` file 
+- `--annotation` [string]: Path to the reference annotation `.gtf`, `.gff`, `.gtf.gz`, or `.gff.gz` file. When restarting the pipeline with `--manual_clustering`, provide the `extended_annotations.rds` file instead 
 
 **Optional**
 - `--output_dir` [string, default: 'output']: Path to the output directory
@@ -141,11 +147,100 @@ To configure the executor and container, pass profile types via the `-profile` a
 - `--qscore_filtering` [boolean, default: true]: Enable or disable quality score filtering of reads
 - `--ndr` [float, default: null]: NDR threshold for Bambu transcript discovery. If not set, Bambu will recommend a suitable value
 - `--deduplicate_umis` [boolean, default: true]: If true, Bambu will perform UMI deduplication 
-- `--quantification_mode` [string, default: "EM_clusters"]: Quantification mode for transcript counts. Available options are:
+- `--quantification_mode` [string, default: "clusteredEM"]: Quantification mode for transcript counts. Available options are:
   - "no_quant": Transcript quantification is not performed
   - "EM": Performs transcript quantification for each cell/spatial coordinate
-  - "EM_clusters": Performs gene expression-based cell clustering using [Seurat](https://satijalab.org/seurat/), followed by transcript quantification at the cluster level
-- `--resolution` [float, default: 0.8]: Seurat clustering resolution
+  - "clusteredEM": Performs gene expression-based cell clustering using [Seurat](https://satijalab.org/seurat/), followed by transcript quantification at the cluster level
+- `--seurat_resolution` [float, default: 0.8]: Seurat clustering resolution
+- `--manual_clustering` [boolean, default: false]: If true, skips clustering and quantifies from cluster assignments generated outside the pipeline (see Advanced Usage section)
+
+**Visium HD**
+
+See the [Visium HD](#visium-hd) section for the required input files and samplesheet format.
+- `--visium_hd` [boolean, default: false]: Enable the Visium HD workflow
+- `--bins` [string, default: null]: Path to a `.csv` listing the resolutions to quantify. Required when `--visium_hd` is set
+- `--barcode_mappings` [string, default: null]: Path to the Spaceranger `barcode_mappings.parquet`. Required when `--visium_hd` is set
+- `--clustering_bin` [integer, default: 8]: Bin size (µm) to cluster at. Must be listed in `--bins`
+- `--banksy` [boolean, default: true]: If true, clusters with [Banksy](https://github.com/prabhakarlab/Banksy) using both expression and spatial position. If false, clusters on gene expression alone
+- `--banksy_lambda` [float, default: 0.8]: Weight given to spatial information
+- `--banksy_k_geom` [integer, default: 50]: Number of spatial neighbours per bin
+
+### **Spatial Analysis**
+
+#### **Visium Spatial Gene Expression**
+The pipeline applies the same processing steps to both 10x Single Cell and Visium Spatial Gene Expression (`visium-v*`) samples. The only difference is that the generated `SummarizedExperiment` objects are attached with spatial mapping information, which is stored in `colData`.
+
+*Example: Spatial Mapping Information*
+
+In the `SummarizedExperiment` object, each row in `colData` contains a spatial barcode together with its X and Y coordinates on the slide. 
+
+| barcode            | x_coordinate | y_coordinate | 
+|:---|:---|:---|
+| AAACAACGAATAGTTC | 17 | 1 |
+| AAACAAGTATCTCCCA  | 103 | 51 |
+| AAACAATCTACTAGCA | 44 | 4 |
+
+
+#### **Visium HD**
+
+This pipeline does not demultiplex Visium HD reads. Barcode assignment and alignment must be performed beforehand using one of the following tools, and the resulting BAM file supplied in the samplesheet:
+
+- [Percula](https://epi2me.nanoporetech.com/epi2me-docs/tools/percula/) for ONT reads
+- [visium-hd-long-reads](https://github.com/10XGenomics/visium-hd-long-reads) for PacBio reads
+
+*Required Input Files*
+
+Alongside the BAM file, Visium HD runs require two Spaceranger outputs:
+
+| Spaceranger output | Description |
+|---|---|
+| `tissue_positions.parquet` | One file per resolution, from `binned_outputs/square_XXXum/spatial/`|
+| `barcode_mappings.parquet` | Maps each 2 µm barcode to its bin at every resolution |
+
+For more information, please refer to the [Space Ranger Spatial Outputs](https://www.10xgenomics.com/support/software/space-ranger/latest/analysis/outputs/spatial-outputs) documentation.
+
+*Samplesheet (CSV)*
+
+Visium HD runs take exactly one sample and only needs the `sample` and `path` columns in the input samplesheet. 
+```csv
+sample,path
+visium_hd_example,examples/visium_hd_example.bam
+```
+
+*Bins Samplesheet (CSV)*
+
+Visium HD captures data at 2 µm resolution, which can be aggregated into larger square bins (e.g. 8 µm or 16 µm). In this pipeline, count matrices are generated at every resolution specified. To specify the resolutions for aggregation, supply a `.csv` file to `--bins`:
+
+```csv
+resolution,tissue_positions
+2,binned_outputs/square_002um/spatial/tissue_positions.parquet
+8,binned_outputs/square_008um/spatial/tissue_positions.parquet
+16,binned_outputs/square_016um/spatial/tissue_positions.parquet
+```
+
+The samplesheet must contain two columns: `resolution`, the bin size in µm, and `tissue_positions`, the path to that resolution's `tissue_positions.parquet` file. A row for the 2 µm resolution is required.
+
+*Running the pipeline*
+
+```bash
+nextflow run main.nf \
+  --input samplesheet_visium_hd.csv \
+  --genome examples/GRCh38.primary_assembly.genome.chr21.fa.gz \
+  --annotation examples/gencode.v49.primary_assembly.annotation.chr21.gtf.gz \
+  --visium_hd \
+  --bins bins.csv \
+  --barcode_mappings barcode_mappings.parquet \
+  --clustering_bin 8 \
+  -profile singularity,hpc
+```
+
+For more information on the Visium HD parameters, please refer to the [Parameters](#parameters) section.
+
+*Clustering*
+
+Clustering is performed at a single resolution only, set by `--clustering_bin`, when `--quantification_mode` is set to `clusteredEM`. By default it uses Seurat with [Banksy](https://github.com/prabhakarlab/Banksy), which combines gene expression with spatial position. Set `--banksy false` to cluster on gene expression alone.
+
+When clustering is skipped under `--quantification_mode EM` (i.e., spot-level EM), transcript counts are produced at every resolution listed in `--bins`.
 
 ### **Output**
 All outputs from the pipeline are written to the directory specified by the `--output_dir` parameter. The pipeline produces per-sample alignment files and the combined transcript discovery and quantification results. 
@@ -153,6 +248,14 @@ All outputs from the pipeline are written to the directory specified by the `--o
 *Output Structure*
 ```
 output/
+├── software_versions.yml
+│
+├── pipeline_info/
+│   ├── execution_timeline.html
+│   ├── execution_report.html
+│   ├── execution_trace.txt
+│   └── pipeline_dag.svg
+│
 ├── bam/                                
 │   ├── <sample>_demultiplexed.bam
 │   └── <sample>_demultiplexed.bam.bai
@@ -173,7 +276,11 @@ output/
 │   ├── barcodes.tsv.gz
 │   └── features.tsv.gz
 │
-│   # single-cell EM:
+├── intermediate_R/
+│   ├── quant_data.rds
+│   └── extended_annotations.rds
+│
+│   # --quantification_mode EM:
 ├── transcript_counts_singlecell/
 │   ├── se_transcript_counts_singlecell.rds
 │   ├── counts.mtx.gz
@@ -183,8 +290,9 @@ output/
 │   ├── barcodes.tsv.gz
 │   └── features.tsv.gz
 │
-│   # cluster-level EM:
+│   # --quantification_mode clusteredEM:
 ├── seurat_obj.rds
+│
 ├── transcript_counts_clusters/
 │   ├── se_transcript_counts_clusters.rds
 │   ├── counts.mtx.gz
@@ -193,20 +301,48 @@ output/
 │   ├── uniqueCounts.mtx.gz
 │   ├── barcodes.tsv.gz
 │   └── features.tsv.gz
-├── gene_counts_clusters/
-│   ├── se_gene_counts_clusters.rds
-│   ├── counts.mtx.gz
-│   ├── CPM.mtx.gz
-│   ├── barcodes.tsv.gz
-│   └── features.tsv.gz
+│
+└── gene_counts_clusters/
+    ├── se_gene_counts_clusters.rds
+    ├── counts.mtx.gz
+    ├── CPM.mtx.gz
+    ├── barcodes.tsv.gz
+    └── features.tsv.gz
+```
+
+*Output Structure (Visium HD)*
+
+Visium HD runs write the count matrices for each resolution specified in the `.csv` provided to `--bins`. Each count directory holds the same files as above (`se_*.rds`, `counts.mtx.gz`, `barcodes.tsv.gz`, `features.tsv.gz`, and any additional assays).
+
+```
+output/
+├── software_versions.yml
 │
 ├── pipeline_info/
-│   ├── execution_timeline.html
-│   ├── execution_report.html
-│   ├── execution_trace.txt
-│   └── pipeline_dag.svg
 │
-└── software_versions.yml
+├── extended_annotations.gtf
+│
+├── unique_counts/
+│   ├── unique_counts_002um/
+│   └── unique_counts_<res>/
+│
+├── gene_counts/
+│   ├── gene_counts_002um/
+│   └── gene_counts_<res>/
+│
+├── intermediate_R/
+│   ├── quant_data.rds
+│   └── extended_annotations.rds
+│
+│   # --quantification_mode EM:
+├── transcript_counts/
+│   ├── transcript_counts_002um/
+│   └── transcript_counts_<res>/
+│
+│   # --quantification_mode clusteredEM:
+├── seurat_obj.rds
+├── transcript_counts_clusters/
+└── gene_counts_clusters/
 ```
 
 **Description of the Output Files**
@@ -215,12 +351,17 @@ output/
 | `<sample>_demultiplexed.bam` | BAM file containing demultiplexed, trimmed and aligned reads
 | `<sample>_demultiplexed.bam.bai` | BAM index for the corresponding BAM file
 | `extended_annotations.gtf` | A `.gtf` file containing the novel transcripts discovered by Bambu as well as the reference annotations provided by the user.
-| `seurat_obj.rds` | A [SeuratObject](https://satijalab.github.io/seurat-object/reference/Seurat-class.html) containing normalised counts, PCA embeddings, and cluster assignments. For multi-sample runs, also contains Harmony-integrated embeddings corrected for sequencing technology and capture chemistry. UMAP has not been computed. Only produced when `--quantification_mode` is set to `EM_clusters`.
+| `seurat_obj.rds` | A [SeuratObject](https://satijalab.github.io/seurat-object/reference/Seurat-class.html) containing normalised counts, PCA embeddings, and cluster assignments. For multi-sample runs, also contains Harmony-integrated embeddings corrected for sequencing technology and capture chemistry. For Visium HD runs, holds the bins of the `--clustering_bin` resolution and, with `--banksy`, the BANKSY assay and embeddings. UMAP has not been computed. Only produced when `--quantification_mode` is set to `clusteredEM`.
 | `unique_counts/se_unique_counts.rds` | A [RangedSummarizedExperiment](https://www.rdocumentation.org/packages/SummarizedExperiment/versions/1.2.3/topics/RangedSummarizedExperiment-class) object containing transcript-level unique counts at single-cell resolution, produced prior to EM quantification. Columns follow the `sampleName_barcode` naming convention.
 | `gene_counts/se_gene_counts.rds` | A RangedSummarizedExperiment object containing gene-level counts at single-cell resolution. Columns follow the `sampleName_barcode` naming convention.
+| `intermediate_R/quant_data.rds` | Bambu's read-to-transcript assignments for the run. Required to restart the pipeline with `--manual_clustering`.
+| `intermediate_R/extended_annotations.rds` | The extended annotations as a GRangesList. Required to restart the pipeline with `--manual_clustering`.
 | `transcript_counts_singlecell/se_transcript_counts_singlecell.rds` | A RangedSummarizedExperiment object containing per-cell transcript counts after EM quantification. Columns follow the `sampleName_barcode` naming convention. Only produced when `--quantification_mode` is set to `EM`.
-| `transcript_counts_clusters/se_transcript_counts_clusters.rds` | A RangedSummarizedExperiment object containing cluster-level transcript counts after EM quantification. Columns follow the `clusterId` naming convention for single-sample runs, and `sampleName_clusterId` for multi-sample runs. Only produced when `--quantification_mode` is set to `EM_clusters`.
-| `gene_counts_clusters/se_gene_counts_clusters.rds` | A RangedSummarizedExperiment object containing cluster-level gene counts. Columns follow the `clusterId` naming convention for single-sample runs, and `sampleName_clusterId` for multi-sample runs. Only produced when `--quantification_mode` is set to `EM_clusters`.
+| `transcript_counts_clusters/se_transcript_counts_clusters.rds` | A RangedSummarizedExperiment object containing cluster-level transcript counts after EM quantification. Columns follow the `clusterId` naming convention for single-sample runs, and `sampleName_clusterId` for multi-sample runs. Only produced when `--quantification_mode` is set to `clusteredEM`.
+| `gene_counts_clusters/se_gene_counts_clusters.rds` | A RangedSummarizedExperiment object containing cluster-level gene counts. Columns follow the `clusterId` naming convention for single-sample runs, and `sampleName_clusterId` for multi-sample runs. Only produced when `--quantification_mode` is set to `clusteredEM`.
+| `unique_counts/unique_counts_<res>/` | Visium HD only. Transcript-level unique counts at resolution `<res>`, produced prior to EM quantification. Columns follow the `sampleName_barcode` naming convention.
+| `gene_counts/gene_counts_<res>/` | Visium HD only. Gene-level counts at resolution `<res>`.
+| `transcript_counts/transcript_counts_<res>/` | Visium HD only. Transcript counts after EM quantification, one directory per resolution. Only produced when `--quantification_mode` is set to `EM`.
 | `<assay>.mtx.gz` | Gzip-compressed sparse count matrix in [Matrix Market Exchange](https://math.nist.gov/MatrixMarket/formats.html) format (features × barcodes), one per assay present in the `SummarizedExperiment` (`counts.mtx.gz`, `CPM.mtx.gz`, `fullLengthCounts.mtx.gz`, `uniqueCounts.mtx.gz`). See the **Count Matrices** section below for what each assay contains.
 | `barcodes.tsv.gz` | Gzip-compressed list of cell barcodes (one per line) corresponding to the columns of each `.mtx.gz`. Barcodes follow the `sampleName_barcode` naming convention for multi-sample runs.
 | `features.tsv.gz` | Gzip-compressed three-column TSV corresponding to the rows of each `.mtx.gz`. Column 1 is the feature ID (transcript ID for transcript-level matrices, gene ID for gene-level matrices). Column 2 is the feature name (same as column 1). Column 3 is the feature type (`Transcript Expression` for transcript-level matrices, `Gene Expression` for gene-level matrices).
@@ -241,62 +382,12 @@ The [RangedSummarizedExperiment](https://www.rdocumentation.org/packages/Summari
 > **Note:** In `se_unique_counts.rds`, unique counts are stored under the `counts` assay, not `uniqueCounts`.
 
 
-### **Spatial Analysis**
-The pipeline applies the same processing steps to both single-cell and spatial samples. However, for spatial data, the generated `SummarizedExperiment` object is appended with spatial mapping information, which is stored in `colData`.  
-
-**Example - Spatial Mapping Information (`visium-v*`)**:
-
-For `visium-v*` samples, `colData` contains the spatial barcode and the corresponding X and Y spatial coordinates. 
-
-| barcode            | x_coordinate | y_coordinate | 
-|:---|:---|:---|
-| AAACAACGAATAGTTC | 17 | 1 |
-| AAACAAGTATCTCCCA  | 103 | 51 |
-| AAACAATCTACTAGCA | 44 | 4 |
-
-
-### **Visium HD Spatial Analysis (Under Development)**
-
-This feature is still under development and will be released in a future update.
-
-
 ### **Fusion Transcript Analysis (Under Development)**
 This feature is still under development and will be released in a future update.
 
 
 ### **Advanced Usage**
 
-**Minimal End-to-End Smoke Test**
-
-Example data and pre-configured profiles are provided in `examples/` to run the pipeline end-to-end automatically without preparing your own data. The commands below must be run from the project's root directory. Combine the profile `test_base` with one of the profiles below and a container profile (`singularity` or `docker`).
-
-| Profile | Description |
-|---|---|
-| `test_sc_fastq` | Single-cell, single-sample ONT run from raw reads |
-| `test_sc_bam` | Single-cell, single-sample ONT run from demultiplexed BAM |
-| `test_sc_multi` | Single-cell, multiple-sample run with different chemistries and technologies |
-| `test_visium` | Spatial (Visium), single-sample ONT run from raw reads |
-| `test_custom` | Custom chemistry, single-sample ONT run from demultiplexed BAM |
-
-```bash
-# Single-cell: test from FASTQ input
-nextflow run . -profile test_base,test_sc_fastq,singularity
-
-# Single-cell: test from BAM input
-nextflow run . -profile test_base,test_sc_bam,singularity
-
-# Single-cell: test with multiple samples (ONT + PacBio)
-nextflow run . -profile test_base,test_sc_multi,singularity
-
-# Spatial: test Visium from FASTQ input
-nextflow run . -profile test_base,test_visium,singularity
-
-# Custom chemistry: test from demultiplexed BAM
-nextflow run . -profile test_base,test_custom,singularity
-```
-
-The output files from the smoke tests are written to `.smoke_test/<profile>/output/`.
-
 **Running Pipeline with a Custom Chemistry or Pre-aligned BAM**
 
 If your dataset uses a chemistry not listed under Supported 10x Library Chemistries, or if you prefer to perform FASTQ preprocessing and genome alignment manually, start the pipeline directly from a pre-processed, demultiplexed BAM file. The BAM file must have the barcode and UMI information encoded either in the `CB`/`UB` column, or in the read name using the format `<barcode>_<umi>#<read_id>`.
@@ -373,18 +464,102 @@ DimPlot(obj, reduction = "umap", group.by = "harmony_clusters", label = TRUE)
 DimPlot(obj, reduction = "umap", group.by = "sample")
 ```
 
-**Manual Clustering (Under Development)**
+**Clustering Outside the Pipeline**
 
-Currently, cell clustering is performed automatically as part of the pipeline. In a future release, a tutorial will be provided that allows users to stop the pipeline after transcript discovery, perform their own custom clustering, and then resume the pipeline to run Bambu transcript quantification using their cluster assignments.
+Set `--manual_clustering` to quantify from cluster assignments you generated yourself instead of the pipeline's Seurat clustering. This is done in two stages.
 
-### **Additional Information**
-UMI correction is done at the barcode level. The longest read for each unique barcode-UMI combination is kept for analysis.
+First, stop the pipeline after transcript discovery:
+
+```bash
+nextflow run . -profile singularity \
+  --input samplesheet.csv --genome genome.fa --annotation annotation.gtf \
+  --quantification_mode no_quant
+```
 
-### **Release History** 
+Cluster the published `gene_counts/` matrix however you like, then write a `.csv` with `id` and `cluster` columns, where each `id` is a line of `gene_counts/barcodes.tsv.gz`:
 
-- v0.1-beta: 2025-May-19
-- v0.9-beta: 2026-May-11
+```csv
+id,cluster
+sample1_AAACCCAAGAAACACT-1,cluster_0
+sample1_AAACCCAAGAAACCAT-1,cluster_1
+```
 
+Next, restart the pipeline. The samplesheet supplied to `--input` must contain two columns: `clusters_path`, the cluster assignment file, and `quant_data_path`, the `quant_data.rds` written to `intermediate_R/` by the first stage.
+
+```csv
+clusters_path,quant_data_path
+clusters.csv,output/intermediate_R/quant_data.rds
+```
+
+```bash
+nextflow run . -profile singularity \
+  --input manual_clusters.csv --genome genome.fa \
+  --annotation output/intermediate_R/extended_annotations.rds \
+  --manual_clustering --output_dir output_manual
+```
+
+`--annotation` takes the `extended_annotations.rds` written to `intermediate_R/` by the first stage.
+
+*Visium HD*
+
+The same two stages described above applies. Cluster the `gene_counts/gene_counts_<res>/` matrix of whichever resolution you want, then restart with `--visium_hd` and set `--clustering_bin` to that resolution:
+
+```bash
+nextflow run . -profile singularity \
+  --input manual_clusters.csv --genome genome.fa \
+  --annotation output/intermediate_R/extended_annotations.rds \
+  --visium_hd --clustering_bin 8 \
+  --barcode_mappings barcode_mappings.parquet \
+  --manual_clustering --output_dir output_manual
+```
+
+`--bins` is not needed on the restart, and `--barcode_mappings` is only needed when `--clustering_bin` is greater than 2.
+
+**Minimal End-to-End Smoke Test**
+
+Example data and pre-configured profiles are provided in `examples/` to run the pipeline end-to-end automatically without preparing your own data. The commands below must be run from the project's root directory. Combine the profile `test_base` with one of the profiles below and a container profile (`singularity` or `docker`).
+
+| Profile | Description |
+|---|---|
+| `test_sc_fastq` | Single-cell, single-sample ONT run from raw reads |
+| `test_sc_bam` | Single-cell, single-sample ONT run from demultiplexed BAM |
+| `test_sc_multi` | Single-cell, multiple-sample run with different chemistries and technologies |
+| `test_visium` | Spatial (Visium), single-sample ONT run from raw reads |
+| `test_visium_hd` | Spatial (Visium HD), single-sample run from demultiplexed BAM |
+| `test_custom` | Custom chemistry, single-sample ONT run from demultiplexed BAM |
+| `test_sc_quant_data` | Single-cell manual clustering run |
+| `test_visium_hd_quant_data` | Visium HD manual clustering run, with clusters made at 8 µm |
+
+```bash
+# Single-cell: test from FASTQ input
+nextflow run . -profile test_base,test_sc_fastq,singularity
+
+# Single-cell: test from BAM input
+nextflow run . -profile test_base,test_sc_bam,singularity
+
+# Single-cell: test with multiple samples (ONT + PacBio)
+nextflow run . -profile test_base,test_sc_multi,singularity
+
+# Spatial: test Visium from FASTQ input
+nextflow run . -profile test_base,test_visium,singularity
+
+# Spatial: test Visium HD from a synthetic barcode-tagged BAM
+nextflow run . -profile test_base,test_visium_hd,singularity
+
+# Custom chemistry: test from demultiplexed BAM
+nextflow run . -profile test_base,test_custom,singularity
+
+# Single-cell: test clustered EM from manually generated clusters
+nextflow run . -profile test_base,test_sc_quant_data,singularity
+
+# Spatial: test Visium HD clustered EM from manually generated clusters
+nextflow run . -profile test_base,test_visium_hd_quant_data,singularity
+```
+
+The output files from the smoke tests are written to `.smoke_test/<profile>/output/`.
+
+### **Additional Information**
+UMI correction is done at the barcode level. The longest read for each unique barcode-UMI combination is kept for analysis.
 
 ### **Citation**
 If you use this pipeline, please cite our paper:
@@ -393,6 +568,9 @@ Sim, A., Ling, M. H., Chen, Y., Lu, H., See, Y. X., Perrin, A., Leng Agnes, O. B
 
 The following are citations for the other tools used in this pipeline:
 
+#### Banksy
+Singhal, V., Chou, N., Lee, J., Yue, Y., Liu, J., Chock, W. K., Lin, L., Chang, Y.-C., Teo, E. M. L., Aow, J., Lee, H. K., Chen, K. H., & Prabhakar, S. (2024, February 27). Banksy unifies cell typing and tissue domain segmentation for scalable spatial omics data analysis. Nature News. https://www.nature.com/articles/s41588-024-01664-3 
+
 #### Chopper
 De Coster Wouter, & Rademakers, R. (2023). NanoPack2: Population scale evaluation of long-read sequencing data. Bioinformatics, 39(5). https://doi.org/10.1093/bioinformatics/btad311
 
diff --git a/bin/bambu_discovery.R b/bin/bambu_discovery.R
new file mode 100755
index 0000000..37cdae8
--- /dev/null
+++ b/bin/bambu_discovery.R
@@ -0,0 +1,31 @@
+# Shared discovery sequence for the standard and Visium HD transcript-discovery
+# processes. Assumes bambu is loaded by the caller.
+# bambuDiscovery() arguments:
+#   reads       named character vector of readClass .rds paths (names = sample
+#               names), or a single path for a one-sample run.
+#   annotation  the bambu transcript annotation object (a GRangesList).
+#   genome      path to the reference genome FASTA.
+#   ncore       integer, number of cores passed to bambu.singlecell.
+#   ndr         numeric, the novel discovery rate (NDR) threshold for discovery.
+#   sampleData  optional path(s) to spatial / bin metadata joined into colData by
+#               the quantData stage; NULL (default) attaches none.
+#
+# Returns a named list with extendedAnno, quantData and seDiscovery (the
+# transcript-level unique-counts SE). No file I/O.
+bambuDiscovery <- function(reads, annotation, genome, ncore, ndr, sampleData = NULL) {
+    # Transcript discovery
+    extendedAnno <- bambu.singlecell(reads = reads, output = "extendedAnnotations",
+        annotations = annotation, genome = genome, ncore = ncore,
+        verbose = FALSE, NDR = ndr)
+
+    # Read to transcript assignment
+    quantData <- bambu.singlecell(reads = reads, output = "quantData",
+        annotations = extendedAnno, genome = genome, ncore = ncore,
+        verbose = FALSE, sampleData = sampleData)
+
+    # Quantification without EM
+    seDiscovery <- bambu.singlecell(reads = quantData, output = "uniqueCounts",
+        annotations = extendedAnno)
+
+    list(extendedAnno = extendedAnno, quantData = quantData, seDiscovery = seDiscovery)
+}
diff --git a/bin/create_seurat_object.R b/bin/create_seurat_object.R
new file mode 100755
index 0000000..7fe3361
--- /dev/null
+++ b/bin/create_seurat_object.R
@@ -0,0 +1,14 @@
+# createSeuratObject() - build a Seurat object from a bambu count directory, so the
+# clustering never has to read a SummarizedExperiment.
+#   countsDir - path to a count directory written by saveCounts()
+#   metadata  - data.frame of cell metadata, one row per barcode with matching row names
+#   assay, minCells, project - passed through to CreateSeuratObject()
+# Assumes Seurat is loaded by the caller
+
+createSeuratObject <- function(countsDir, assay = "RNA", metadata = NULL, minCells = 0, project = "SeuratProject") {
+    counts <- ReadMtx(mtx      = file.path(countsDir, "counts.mtx.gz"),
+                      cells    = file.path(countsDir, "barcodes.tsv.gz"),
+                      features = file.path(countsDir, "features.tsv.gz"))
+
+    CreateSeuratObject(counts, assay = assay, meta.data = metadata, min.cells = minCells, project = project)
+}
diff --git a/bin/map_spots_to_clusters.R b/bin/map_spots_to_clusters.R
new file mode 100755
index 0000000..988620c
--- /dev/null
+++ b/bin/map_spots_to_clusters.R
@@ -0,0 +1,20 @@
+# Cluster assignments arrive either from Seurat or from the user, and bambu quantifies
+# at 2um, so coarser bin labels have to be expanded down to the spots it quantified.
+# Assumes data.table is available to the caller.
+
+# readClusters() - load a cluster assignment as a named vector of id -> cluster label
+#   path - a .rds holding a named vector, or a .csv/.tsv/.txt with id and cluster columns
+readClusters <- function(path) {
+    if (endsWith(path, ".rds")) return(readRDS(path))
+
+    clusterDf <- data.table::fread(path, data.table = FALSE)
+    setNames(as.character(clusterDf$cluster), clusterDf$id)
+}
+
+# Seurat stores the map of bin-level barcode (e.g., 8um) -> cluster. To allow Bambu,
+# to run quantification, we have to generate the 2um barocde -> cluster map
+# as quantData is generated at the 2um resolution
+mapSpotsToClusters <- function(binToClustersMap, spotMappings) {
+    spotToClustersMap <- setNames(unname(binToClustersMap[spotMappings$sample_bin]), spotMappings$sample_barcode)
+    spotToClustersMap[!is.na(spotToClustersMap)]
+}
diff --git a/bin/save_counts.R b/bin/save_counts.R
index eb01f8a..72b12a7 100755
--- a/bin/save_counts.R
+++ b/bin/save_counts.R
@@ -1,17 +1,17 @@
-# save_counts() - write a SummarizedExperiment as a count directory: one
+# saveCounts() - write a SummarizedExperiment as a count directory: one
 # MTX file per assay, plus a single shared barcodes.tsv.gz / features.tsv.gz,
-# and the SummarisedExperiment object. 
+# and the SummarisedExperiment object.
 # Assumes DropletUtils and bambu have been loaded
 
-save_counts <- function(se, dir, gene.type = "Gene Expression") {
-    write10xCounts(dir, assays(se)$counts, version = "3", gene.type = gene.type)
+saveCounts <- function(se, dir, geneType = "Gene Expression") {
+    write10xCounts(dir, assays(se)$counts, version = "3", gene.type = geneType)
     file.rename(file.path(dir, "matrix.mtx.gz"), file.path(dir, "counts.mtx.gz"))
-    extra_assays <- setdiff(assayNames(se), "counts")
-    for (assay_name in extra_assays) {
-        mat      <- as(assays(se)[[assay_name]], "CsparseMatrix")
-        mtx_path <- file.path(dir, paste0(assay_name, ".mtx"))
-        Matrix::writeMM(mat, mtx_path)
-        R.utils::gzip(mtx_path, overwrite = TRUE)
+    extraAssays <- setdiff(assayNames(se), "counts")
+    for (assayName in extraAssays) {
+        mat     <- as(assays(se)[[assayName]], "CsparseMatrix")
+        mtxPath <- file.path(dir, paste0(assayName, ".mtx"))
+        Matrix::writeMM(mat, mtxPath)
+        R.utils::gzip(mtxPath, overwrite = TRUE)
     }
     saveRDS(se, file.path(dir, paste0("se_", dir, ".rds")))
 }
diff --git a/bin/visium_hd_aggregate_resolution.R b/bin/visium_hd_aggregate_resolution.R
new file mode 100755
index 0000000..753d846
--- /dev/null
+++ b/bin/visium_hd_aggregate_resolution.R
@@ -0,0 +1,55 @@
+# Assumes SummarizedExperiment, Matrix and dplyr are loaded by the caller.
+# A spot is a 2um barcode, i.e. one column of se; a bin is the coarser square
+# (e.g. 8um) that several spots fall into.
+# Build a bin-resolution SummarizedExperiment from a 2um-resolution SE by summing
+# counts across the spots that share each bin. spotMappings gives the spot -> bin
+# assignment for one resolution; tissuePositions supplies the new colData,
+# alongside the id and sampleName carried over from the 2um SE.
+# rowRanges/rowData and metadata are carried over.
+aggregateResolution <- function(se, tissuePositions, spotMappings) {
+    # Every column of se is one spot, so look up the bin each spot falls in. The
+    # prefixed forms are used so the aggregated columns are named like bambu ids.
+    matchIdx     <- match(colData(se)$id, spotMappings$sample_barcode)
+    spotToBinMap <- factor(spotMappings$sample_bin[matchIdx])
+
+    # The bins become the columns of the aggregated matrix, in level order
+    binIds      <- levels(spotToBinMap) # bambu ids, e.g. <sample>_s_008um_<array_row>_<array_col>-1
+    spotsPerBin <- t(fac2sparse(spotToBinMap)) # spots x bins, 1 where a spot falls in a bin
+
+    # Sum each bin's spots: (features x spots) %*% (spots x bins) = features x bins
+    spotCounts <- assays(se)$counts
+    binCounts  <- as(spotCounts %*% spotsPerBin, "CsparseMatrix")
+
+    ## build new colData for the aggregated SE object
+    
+    # bin id -> barcode (without sampleName prefix), e.g. sample_s_008um_00247_00090-1 -> s_008um_00247_00090-1
+    binIdToBarcodeMap <- distinct(spotMappings, sample_bin, barcode = bin)
+
+    # bin id -> the 2um barcodes the bin contains
+    binToSpotsMap <- data.frame(sample_bin = as.character(spotToBinMap), barcode = colData(se)$barcode) %>%
+        group_by(sample_bin) %>%
+        summarise(barcodes = list(barcode), .groups = "drop")
+
+    # creates dataframe containing one row per bin, each row contains metadata information (e.g., barcode, in_tissue, array and pixel coordinates)
+    metadataPerBin <- data.frame(sample_bin = binIds) %>%
+        left_join(binIdToBarcodeMap, by = "sample_bin") %>%
+        left_join(binToSpotsMap,     by = "sample_bin") %>%
+        left_join(tissuePositions,   by = "barcode") %>%
+        select(-sample_bin)
+
+    sampleName <- unique(colData(se)$sampleName)
+    binColData <- DataFrame(id = binIds, sampleName = sampleName, metadataPerBin, row.names = binIds)
+
+    binSe <- SummarizedExperiment(assays = SimpleList(counts = binCounts),
+                                  rowRanges = rowRanges(se),
+                                  colData = binColData)
+
+    # Carry the bambu metadata matrices over, summed the same way
+    incompatiblePerSpot <- metadata(se)$incompatibleCounts
+    nonuniquePerSpot    <- metadata(se)$nonuniqueCounts
+
+    metadata(binSe)$incompatibleCounts <- as(incompatiblePerSpot %*% spotsPerBin, "CsparseMatrix")
+    metadata(binSe)$nonuniqueCounts    <- as(nonuniquePerSpot %*% spotsPerBin, "CsparseMatrix")
+    metadata(binSe)$seType             <- metadata(se)$seType
+    binSe
+}
diff --git a/bin/visium_hd_cluster.R b/bin/visium_hd_cluster.R
new file mode 100755
index 0000000..e2fff5e
--- /dev/null
+++ b/bin/visium_hd_cluster.R
@@ -0,0 +1,51 @@
+# Seurat clustering of a single Visium HD bin resolution, adapted from
+# https://satijalab.org/seurat/articles/visiumhd_analysis_vignette
+# Assumes Seurat is loaded by the caller, plus SeuratWrappers and Banksy
+# when the spatially aware path is used.
+
+# Spatially aware clustering using Banksy (Refer to vignette for more information)
+clusterBanksy <- function(object, assay, lambda, kGeom, clusterResolution, npcs = 30) {
+    # dimx/dimy name the meta.data columns carrying each bin's coordinates
+    object <- RunBanksy(object, lambda = lambda, assay = assay, slot = "data",
+                        features = "variable", k_geom = kGeom,
+                        dimx = "pxl_col_in_fullres", dimy = "pxl_row_in_fullres", verbose = FALSE)
+
+    DefaultAssay(object) <- "BANKSY"
+    npcs   <- min(npcs, ncol(object) - 1)
+    object <- RunPCA(object, assay = "BANKSY", reduction.name = "pca.banksy",
+                     features = rownames(object), npcs = npcs, verbose = FALSE)
+    object <- FindNeighbors(object, reduction = "pca.banksy", dims = 1:npcs, verbose = FALSE)
+    FindClusters(object, cluster.name = "clusters", resolution = clusterResolution, verbose = FALSE)
+}
+
+# Expression-only clustering (Refer to vignette for more information)
+clusterExpression <- function(object, assay, clusterResolution, dims = 15, ncells = 50000) {
+    object <- FindVariableFeatures(object, verbose = FALSE)
+    object <- ScaleData(object, verbose = FALSE)
+
+    sketched <- ncol(object) > ncells
+    if (sketched) {
+        object <- SketchData(object, ncells = ncells, method = "LeverageScore", sketched.assay = "sketch")
+        DefaultAssay(object) <- "sketch"
+        object <- FindVariableFeatures(object, verbose = FALSE)
+        object <- ScaleData(object, verbose = FALSE)
+    }
+
+    reduction <- if (sketched) "pca.sketch" else "pca"
+    npcs      <- min(dims, ncol(object) - 1)
+    object    <- RunPCA(object, reduction.name = reduction, npcs = npcs, verbose = FALSE)
+    object    <- FindNeighbors(object, reduction = reduction, dims = 1:npcs, verbose = FALSE)
+    object    <- FindClusters(object, cluster.name = "clusters", resolution = clusterResolution, verbose = FALSE)
+
+    if (sketched) {
+        # project the sketched labels onto the full set of bins, then adopt them as the
+        # cluster assignment so both paths emit the same 'clusters' column
+        object <- ProjectData(object, assay = assay, full.reduction = "full.pca.sketch",
+                              sketched.assay = "sketch", sketched.reduction = "pca.sketch",
+                              dims = 1:npcs, refdata = list(clusters.projected = "clusters"))
+        DefaultAssay(object) <- assay
+        object$clusters      <- object$clusters.projected
+    }
+
+    object
+}
diff --git a/bin/visium_hd_convert_barcode_mappings.py b/bin/visium_hd_convert_barcode_mappings.py
new file mode 100755
index 0000000..42c58d4
--- /dev/null
+++ b/bin/visium_hd_convert_barcode_mappings.py
@@ -0,0 +1,36 @@
+#!/usr/bin/env python3
+import argparse
+import sys
+
+import pandas as pd
+import pyarrow.parquet as pq
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Subset a Visium HD barcode_mappings.parquet to a single bin resolution, "
+        "writing a CSV of the 2um-barcode -> bin assignment with columns 'barcode,bin'"
+    )
+    parser.add_argument("barcode_mappings", help="Path to barcode_mappings.parquet")
+    parser.add_argument("resolution", help="Bin resolution token, e.g. 008um, 016um")
+    parser.add_argument("output", help="Output path for the converted CSV")
+    args = parser.parse_args()
+
+    # read the schema first so only the two needed columns are loaded from the parquet;
+    # the file has one row per 2um barcode and one column per resolution
+    schema = pq.read_schema(args.barcode_mappings)
+    barcode_column = schema.names[0]
+    bin_column = f"square_{args.resolution}"
+    if bin_column not in schema.names:
+        sys.exit(
+            f"error: column '{bin_column}' not found in {args.barcode_mappings} "
+            f"(available: {', '.join(schema.names)})"
+        )
+
+    barcode_mappings = pd.read_parquet(args.barcode_mappings, columns=[barcode_column, bin_column])
+    barcode_mappings = barcode_mappings.rename(columns={barcode_column: "barcode", bin_column: "bin"})
+    barcode_mappings.to_csv(args.output, index=False)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/bin/visium_hd_convert_tissue_positions.py b/bin/visium_hd_convert_tissue_positions.py
new file mode 100755
index 0000000..8bfe917
--- /dev/null
+++ b/bin/visium_hd_convert_tissue_positions.py
@@ -0,0 +1,25 @@
+#!/usr/bin/env python3
+import argparse
+
+import pandas as pd
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Convert a Visium HD tissue_positions.parquet to CSV; for the 2um resolution also write the in-tissue barcode list used to filter the BAM")
+    parser.add_argument("tissue_positions", help="Path to tissue_positions.parquet")
+    parser.add_argument("resolution", help="Bin resolution token, e.g. 002um, 008um, 016um")
+    parser.add_argument("csv", help="Output path for the converted CSV")
+    parser.add_argument("barcodes", help="Output path for the in-tissue barcode list (written only for the 002um resolution)")
+    args = parser.parse_args()
+
+    tissue_positions = pd.read_parquet(args.tissue_positions)
+    tissue_positions.to_csv(args.csv, index=False)
+
+    # the BAM carries 2um barcodes, so only the 2um resolution feeds the samtools tissue filter
+    if args.resolution == "002um":
+        in_tissue = tissue_positions.loc[tissue_positions["in_tissue"] == 1, "barcode"]
+        in_tissue.to_csv(args.barcodes, index=False, header=False)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/bin/visium_hd_spot_bin_mappings.py b/bin/visium_hd_spot_bin_mappings.py
new file mode 100755
index 0000000..486363d
--- /dev/null
+++ b/bin/visium_hd_spot_bin_mappings.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env python3
+import argparse
+
+import pandas as pd
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Map every 2um spot bambu quantified to its bin at one resolution")
+    parser.add_argument("barcode_mappings", help="Path to the resolution's barcode mappings CSV (columns: barcode,bin)")
+    parser.add_argument("barcodes", help="Path to barcodes.tsv.gz from the 2um count directory; one bambu id per column")
+    parser.add_argument("sample", help="Sample name that bambu prefixes onto each barcode to form the id")
+    parser.add_argument("output", help="Output path for the spot to bin mapping CSV")
+    args = parser.parse_args()
+
+    # bambu ids are sampleName_barcode, so strip the prefix to get back to the mapping's barcodes
+    prefix = f"{args.sample}_"
+    spot_barcodes = pd.read_csv(args.barcodes, header=None)[0].str.removeprefix(prefix)
+
+    # keep the bins of the spots bambu quantified, then add the prefixed forms so
+    # downstream modules never have to rebuild a bambu id
+    barcode_mappings = pd.read_csv(args.barcode_mappings)
+    spot_mappings = barcode_mappings[barcode_mappings["barcode"].isin(spot_barcodes)].copy()
+    spot_mappings["sample_barcode"] = prefix + spot_mappings["barcode"]
+    spot_mappings["sample_bin"] = prefix + spot_mappings["bin"]
+
+    spot_mappings[["sample_barcode", "barcode", "sample_bin", "bin"]].to_csv(args.output, index=False)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/conf/base.config b/conf/base.config
index ae2574c..76247ad 100644
--- a/conf/base.config
+++ b/conf/base.config
@@ -19,5 +19,5 @@ process {
     withLabel: 'long'           { time = { 12.h * task.attempt } }
 
     // R processes exit with code 1 on BiocParallel/OOM failures (not signal-based codes)
-    withLabel: 'r'              { errorStrategy = { task.exitStatus in [1, 130, 137, 139, 143] ? 'retry' : 'finish' } }
+    withLabel: 'bambu|seurat'   { errorStrategy = { task.exitStatus in [1, 130, 137, 139, 143] ? 'retry' : 'finish' } }
 }
\ No newline at end of file
diff --git a/conf/containers.config b/conf/containers.config
index 3cac435..f5fd818 100644
--- a/conf/containers.config
+++ b/conf/containers.config
@@ -1,6 +1,8 @@
 process {
-    withLabel: 'r'                  { container = "ghcr.io/goekelab/bambu-pipe-r:1.0.0" }
+    withLabel: 'bambu'              { container = "ghcr.io/goekelab/bambu-pipe-bambu:1.0.0" }
+    withLabel: 'seurat'             { container = "ghcr.io/goekelab/bambu-pipe-seurat:1.0.0" }
     withLabel: 'spaceranger'        { container = "quay.io/nf-core/spaceranger:9c5e7dc93c32448e" }
     withLabel: 'minimap2_samtools'  { container = "community.wave.seqera.io/library/minimap2_samtools:b09096fc890429ce" }
     withLabel: 'preprocess'         { container = "community.wave.seqera.io/library/chopper_cutadapt_flexiplex_pigz:077c3bc67452482c" }
-}
\ No newline at end of file
+    withLabel: 'pyarrow_pandas'     { container = "community.wave.seqera.io/library/pip_pandas_pyarrow:e60a5c578ef189b2" }
+}
diff --git a/conf/dev.config b/conf/dev.config
index d91ddae..ba53fc2 100644
--- a/conf/dev.config
+++ b/conf/dev.config
@@ -3,7 +3,7 @@ params {
     bambu_path                  = null
     valid_chemistries           = ['10x3v2', '10x3v3', '10x3v4', '10x5v2', '10x5v3', 'visium-v1', 'visium-v2', 'visium-v3', 'visium-v4', 'visium-v5']
     valid_technologies          = ['ONT', 'PacBio']
-    valid_quantification_modes  = ['no_quant', 'EM', 'EM_clusters']
+    valid_quantification_modes  = ['no_quant', 'EM', 'clusteredEM']
     save_intermediates          = false
     qfilter_threshold           = 10
     flexiplex_f_5prime          = 8
diff --git a/conf/smoke_test.config b/conf/smoke_test.config
index 8576796..4d87e50 100644
--- a/conf/smoke_test.config
+++ b/conf/smoke_test.config
@@ -46,6 +46,15 @@ profiles {
         }
     }
 
+    test_sc_quant_data {
+        params {
+            input             = "${projectDir}/examples/samplesheet_test_sc_quant_data.csv"
+            annotation        = "${projectDir}/examples/extended_annotations_sc_example.rds"
+            output_dir        = "${projectDir}/.smoke_test/test_sc_quant_data/output"
+            manual_clustering = true
+        }
+    }
+
     test_visium {
         params {
             input      = "${projectDir}/examples/samplesheet_test_visium.csv"
@@ -53,6 +62,30 @@ profiles {
         }
     }
 
+    test_visium_hd {
+        params {
+            input            = "${projectDir}/examples/samplesheet_visium_hd_synthetic.csv"
+            bins             = "${projectDir}/examples/bins_visium_hd_synthetic.csv"
+            barcode_mappings = "${projectDir}/examples/barcode_mappings_synthetic.parquet"
+            output_dir       = "${projectDir}/.smoke_test/test_visium_hd/output"
+            visium_hd        = true
+            clustering_bin   = 8
+            banksy           = true
+        }
+    }
+
+    test_visium_hd_quant_data {
+        params {
+            input             = "${projectDir}/examples/samplesheet_visium_hd_quant_data.csv"
+            annotation        = "${projectDir}/examples/extended_annotations_visium_hd_synthetic.rds"
+            barcode_mappings  = "${projectDir}/examples/barcode_mappings_synthetic.parquet"
+            output_dir        = "${projectDir}/.smoke_test/test_visium_hd_quant_data/output"
+            visium_hd         = true
+            manual_clustering = true
+            clustering_bin    = 8
+        }
+    }
+
     test_custom {
         params {
             input      = "${projectDir}/examples/samplesheet_test_custom.csv"
diff --git a/containers/r/Dockerfile b/containers/bambu/Dockerfile
similarity index 84%
rename from containers/r/Dockerfile
rename to containers/bambu/Dockerfile
index fcdd656..c98b8b1 100644
--- a/containers/r/Dockerfile
+++ b/containers/bambu/Dockerfile
@@ -7,22 +7,20 @@ RUN apt-get update && apt-get install -y libcurl4-openssl-dev procps git \
 ARG MAMBA_DOCKERFILE_ACTIVATE=1
 
 # Install bioconductor version of bambu for its dependencies
-# TODO: data.table 1.17.8 installation is a temp fix for clustered_EM 
+# TODO: data.table 1.17.8 installation is a temp fix for clustered_EM
 RUN micromamba install -y -n base -c conda-forge -c bioconda \
     r-base=4.5.3 \
     r-devtools=2.5.2 \
     r-biocmanager=1.30.27 \
     bioconductor-bambu=3.12.1 \
     bioconductor-dropletutils=1.30.0 \
-    r-seurat=5.4.0 \
-    r-harmony=2.0.2 \
     r-data.table=1.17.8 \
     gxx=15.2.0 \
     && micromamba clean --all --yes
 
-# Clone bambu single cell feature branch 
+# Clone bambu single cell feature branch
 RUN git clone --branch devel_pre_v4 \
     https://github.com/GoekeLab/bambu.git /opt/bambu \
     && R -e 'devtools::install("/opt/bambu")'
 
-ENV PATH="$MAMBA_ROOT_PREFIX/bin:$PATH"
\ No newline at end of file
+ENV PATH="$MAMBA_ROOT_PREFIX/bin:$PATH"
diff --git a/containers/seurat/Dockerfile b/containers/seurat/Dockerfile
new file mode 100644
index 0000000..ce6d557
--- /dev/null
+++ b/containers/seurat/Dockerfile
@@ -0,0 +1,13 @@
+FROM rocker/r-ver:4.5.2
+
+RUN apt-get update && apt-get install -y --no-install-recommends \
+        procps \
+    && rm -rf /var/lib/apt/lists/*
+
+RUN R -e 'install.packages("pak")'
+
+RUN R -e 'options(timeout = 1800); \
+          pak::pkg_install(c("Seurat", "harmony", "R.utils", "bioc::Banksy", \
+                             "github::satijalab/seurat-wrappers@9586c29"))'
+
+RUN R -e 'library(Seurat); library(SeuratWrappers); library(Banksy); library(harmony)'
\ No newline at end of file
diff --git a/examples/barcode_mappings_synthetic.parquet b/examples/barcode_mappings_synthetic.parquet
new file mode 100644
index 0000000..851e293
Binary files /dev/null and b/examples/barcode_mappings_synthetic.parquet differ
diff --git a/examples/bins_visium_hd_synthetic.csv b/examples/bins_visium_hd_synthetic.csv
new file mode 100644
index 0000000..6cfbd98
--- /dev/null
+++ b/examples/bins_visium_hd_synthetic.csv
@@ -0,0 +1,4 @@
+resolution,tissue_positions
+2,examples/tissue_positions_002um_synthetic.parquet
+8,examples/tissue_positions_008um_synthetic.parquet
+16,examples/tissue_positions_016um_synthetic.parquet
diff --git a/examples/clusters_sc_example.rds b/examples/clusters_sc_example.rds
new file mode 100644
index 0000000..4dd32d7
Binary files /dev/null and b/examples/clusters_sc_example.rds differ
diff --git a/examples/clusters_visium_hd_synthetic.rds b/examples/clusters_visium_hd_synthetic.rds
new file mode 100644
index 0000000..fd5b788
Binary files /dev/null and b/examples/clusters_visium_hd_synthetic.rds differ
diff --git a/examples/extended_annotations_sc_example.rds b/examples/extended_annotations_sc_example.rds
new file mode 100644
index 0000000..a7791d2
Binary files /dev/null and b/examples/extended_annotations_sc_example.rds differ
diff --git a/examples/extended_annotations_visium_hd_synthetic.rds b/examples/extended_annotations_visium_hd_synthetic.rds
new file mode 100644
index 0000000..35581c7
Binary files /dev/null and b/examples/extended_annotations_visium_hd_synthetic.rds differ
diff --git a/examples/quant_data_sc_example.rds b/examples/quant_data_sc_example.rds
new file mode 100644
index 0000000..f05ae04
Binary files /dev/null and b/examples/quant_data_sc_example.rds differ
diff --git a/examples/quant_data_visium_hd_synthetic.rds b/examples/quant_data_visium_hd_synthetic.rds
new file mode 100644
index 0000000..d60edf5
Binary files /dev/null and b/examples/quant_data_visium_hd_synthetic.rds differ
diff --git a/examples/samplesheet_test_sc_quant_data.csv b/examples/samplesheet_test_sc_quant_data.csv
new file mode 100644
index 0000000..2a8a9d2
--- /dev/null
+++ b/examples/samplesheet_test_sc_quant_data.csv
@@ -0,0 +1,2 @@
+clusters_path,quant_data_path
+examples/clusters_sc_example.rds,examples/quant_data_sc_example.rds
\ No newline at end of file
diff --git a/examples/samplesheet_visium_hd_quant_data.csv b/examples/samplesheet_visium_hd_quant_data.csv
new file mode 100644
index 0000000..d99ac36
--- /dev/null
+++ b/examples/samplesheet_visium_hd_quant_data.csv
@@ -0,0 +1,2 @@
+clusters_path,quant_data_path
+examples/clusters_visium_hd_synthetic.rds,examples/quant_data_visium_hd_synthetic.rds
diff --git a/examples/samplesheet_visium_hd_synthetic.csv b/examples/samplesheet_visium_hd_synthetic.csv
new file mode 100644
index 0000000..d3290bd
--- /dev/null
+++ b/examples/samplesheet_visium_hd_synthetic.csv
@@ -0,0 +1,2 @@
+sample,path
+visium_hd_synthetic,examples/visium_hd_synthetic.bam
diff --git a/examples/tissue_positions_002um_synthetic.parquet b/examples/tissue_positions_002um_synthetic.parquet
new file mode 100644
index 0000000..1e8727b
Binary files /dev/null and b/examples/tissue_positions_002um_synthetic.parquet differ
diff --git a/examples/tissue_positions_008um_synthetic.parquet b/examples/tissue_positions_008um_synthetic.parquet
new file mode 100644
index 0000000..e9c78ef
Binary files /dev/null and b/examples/tissue_positions_008um_synthetic.parquet differ
diff --git a/examples/tissue_positions_016um_synthetic.parquet b/examples/tissue_positions_016um_synthetic.parquet
new file mode 100644
index 0000000..897e1df
Binary files /dev/null and b/examples/tissue_positions_016um_synthetic.parquet differ
diff --git a/examples/visium_hd_synthetic.bam b/examples/visium_hd_synthetic.bam
new file mode 100644
index 0000000..992a7a3
Binary files /dev/null and b/examples/visium_hd_synthetic.bam differ
diff --git a/lib/Validation.groovy b/lib/Validation.groovy
index bc7ae47..e40f304 100644
--- a/lib/Validation.groovy
+++ b/lib/Validation.groovy
@@ -23,11 +23,90 @@ class Validation {
             throw new Exception("Invalid params.quantification_mode '${params.quantification_mode}' — must be one of: ${params.valid_quantification_modes.join(', ')}")
 
         // Numeric range checks
-        if (params.resolution <= 0)
-            throw new Exception("Invalid params.resolution '${params.resolution}' — must be a positive number")
+        if (params.seurat_resolution <= 0)
+            throw new Exception("Invalid params.seurat_resolution '${params.seurat_resolution}' — must be a positive number")
 
         if (params.ndr != null && (params.ndr < 0 || params.ndr > 1))
             throw new Exception("Invalid params.ndr '${params.ndr}' — must be a float between 0 and 1")
+
+        // Visium HD checks
+        if (params.visium_hd) {
+            // the manual path only needs the mappings when it has to expand bins down to 2um
+            def needsBarcodeMappings = !params.manual_clustering || params.clustering_bin != 2
+
+            if (needsBarcodeMappings) {
+                if (params.barcode_mappings == null)
+                    throw new Exception("params.barcode_mappings is required when params.visium_hd is true")
+
+                if (!params.barcode_mappings.exists())
+                    throw new Exception("params.barcode_mappings '${params.barcode_mappings}' does not exist")
+
+                if (params.barcode_mappings.extension != 'parquet')
+                    throw new Exception("params.barcode_mappings '${params.barcode_mappings}' must be a .parquet file")
+            }
+
+            // the bins samplesheet only drives the automatic path
+            if (!params.manual_clustering) {
+                if (params.bins == null)
+                    throw new Exception("params.bins is required when params.visium_hd is true")
+
+                if (!params.bins.exists())
+                    throw new Exception("params.bins '${params.bins}' does not exist")
+
+                if (params.bins.extension != 'csv')
+                    throw new Exception("params.bins '${params.bins}' must be a CSV file")
+
+                validateVisiumHDBins(params.bins)
+
+                // only the clustering mode reads the clustering bin, so only it has to resolve
+                if (params.quantification_mode == 'clusteredEM')
+                    validateClusteringBin(params.bins, params.clustering_bin)
+            }
+        }
+    }
+
+    static def validateVisiumHDRows(rows) {
+        if (rows.size() != 1)
+            throw new Exception("Visium HD requires exactly 1 sample in the samplesheet, but found ${rows.size()}")
+
+        def row = rows[0]
+        ["sample", "path"].each { col ->
+            if (!row.containsKey(col))
+                throw new Exception("Samplesheet is missing a required '${col}' column")
+            if (!row[col])
+                throw new Exception("A row in the samplesheet has an empty '${col}' value")
+        }
+
+        if (!row.path.endsWith('.bam'))
+            throw new Exception("Visium HD sample '${row.sample}' must point to a pre-aligned, barcode-tagged BAM file")
+    }
+
+    static def validateVisiumHDBins(bins) {
+        def lines = bins.text.readLines().findAll { line -> line.trim() }
+        if (lines.size() < 2)
+            throw new Exception("params.bins '${bins}' must have a header row and at least one resolution row")
+
+        def header = lines[0].split(',').collect { col -> col.trim() }
+        ["resolution", "tissue_positions"].each { col ->
+            if (!header.contains(col))
+                throw new Exception("params.bins '${bins}' is missing a required '${col}' column")
+        }
+
+        if (!readBinResolutions(bins).contains("2"))
+            throw new Exception("params.bins '${bins}' must include a row for the native 2 um resolution (resolution = 2)")
+    }
+
+    static def readBinResolutions(bins) {
+        def lines = bins.text.readLines().findAll { line -> line.trim() }
+        def header = lines[0].split(',').collect { col -> col.trim() }
+        def resIdx = header.indexOf("resolution")
+        lines[1..-1].collect { line -> line.split(',')[resIdx].trim() }
+    }
+
+    static def validateClusteringBin(bins, clusteringBin) {
+        def resolutions = readBinResolutions(bins)
+        if (!resolutions.contains(clusteringBin.toString()))
+            throw new Exception("params.clustering_bin '${clusteringBin}' is not listed in params.bins '${bins}' — available bins: ${resolutions.join(', ')}")
     }
 
     static def validateVisiumSampleCount(samples) {
diff --git a/main.nf b/main.nf
index 2b45d93..9237623 100644
--- a/main.nf
+++ b/main.nf
@@ -2,16 +2,29 @@
 
 nextflow.enable.types = true
 
-include { DECOMPRESS as DECOMPRESS_GENOME }     from './modules/decompress.nf'
-include { DECOMPRESS as DECOMPRESS_ANNOTATION } from './modules/decompress.nf'
-include { PREPARE_INPUT_STANDARD } from './subworkflows/prepare_input_standard.nf'
-include { PREPROCESS_FASTQ } from './modules/preprocess_fastq.nf'
-include { ALIGNMENT } from './subworkflows/alignment.nf'
-include { BAMBU_CONSTRUCT_READ_CLASS } from './modules/bambu/construct_read_class.nf'
-include { BAMBU_PREPARE_ANNOTATION } from './modules/bambu/prepare_annotation.nf'
-include { BAMBU_TRANSCRIPT_DISCOVERY } from './modules/bambu/transcript_discovery.nf'
-include { CLUSTERING } from './subworkflows/clustering.nf'
-include { BAMBU_EM } from './modules/bambu/EM_quant.nf'
+// subworkflows
+include { PREPARE_INPUT_STANDARD }               from './subworkflows/prepare_input_standard.nf'
+include { PREPARE_INPUT_VISIUM_HD }              from './subworkflows/prepare_input_visium_hd.nf'
+include { ALIGNMENT }                            from './subworkflows/alignment.nf'
+include { CLUSTERING }                           from './subworkflows/clustering_standard.nf'
+
+// modules
+include { DECOMPRESS as DECOMPRESS_GENOME }      from './modules/decompress.nf'
+include { DECOMPRESS as DECOMPRESS_ANNOTATION }  from './modules/decompress.nf'
+include { PREPROCESS_FASTQ }                     from './modules/preprocess_fastq.nf'
+include { BAMBU_PREPARE_ANNOTATION }             from './modules/bambu/shared/prepare_annotation.nf'
+include { BAMBU_CONSTRUCT_READ_CLASS }           from './modules/bambu/shared/construct_read_class.nf'
+include { BAMBU_CLUSTER_LEVEL_QUANTIFICATION }   from './modules/bambu/shared/cluster_level_quantification.nf'
+include { BAMBU_TRANSCRIPT_DISCOVERY }           from './modules/bambu/standard/transcript_discovery.nf'
+include { BAMBU_SINGLE_CELL_QUANTIFICATION }     from './modules/bambu/standard/single_cell_quantification.nf'
+include { BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD } from './modules/bambu/visium_hd/transcript_discovery.nf'
+include { BAMBU_SPOT_LEVEL_QUANTIFICATION }      from './modules/bambu/visium_hd/spot_level_quantification.nf'
+include { AGGREGATE_BINS_VISIUM_HD }             from './modules/bambu/visium_hd/aggregate_bins.nf'
+include { EXTRACT_SPOT_BARCODES }                from './modules/bambu/visium_hd/extract_spot_barcodes.nf'
+include { MAP_CLUSTERS_TO_2UM_SPOTS }            from './modules/bambu/visium_hd/map_clusters_to_2um_spots.nf'
+include { CONVERT_BARCODE_MAPPINGS }             from './modules/prepare_input/visium_hd/convert_barcode_mappings.nf'
+include { SPOT_BIN_MAPPINGS }                    from './modules/prepare_input/visium_hd/spot_bin_mappings.nf'
+include { SEURAT_VISIUM_HD }                     from './modules/seurat/visium_hd/clustering.nf'
 
 params {
     input: Path
@@ -25,14 +38,178 @@ params {
     ndr: Float?
     deduplicate_umis: Boolean
     quantification_mode: String
-    resolution: Float
+    seurat_resolution: Float
+    visium_hd: Boolean
+    barcode_mappings: Path?
+    bins: Path?
+    clustering_bin: Integer
+    banksy: Boolean
+    banksy_lambda: Float
+    banksy_k_geom: Integer
+    manual_clustering: Boolean
+}
+
+workflow STANDARD {
+    take:
+    ch_rows: Channel<Map>
+    ch_genome: Path
+    ch_annotation: Path
+    ndr: Float?
+    manual_clustering: Boolean
+
+    main:
+    if (!manual_clustering) {
+        def ndrArg = ndr != null ? ndr : 'NULL'
+
+        // load config files
+        ch_barcode_coordinate_config = file("${projectDir}/assets/10x_config/barcode_coordinate_config.csv", checkIfExists: true)
+        ch_adapter_seq_config = file("${projectDir}/assets/10x_config/adapter_seq_config.csv", checkIfExists: true)
+        ch_flank_seq_config = file("${projectDir}/assets/10x_config/flank_seq_config.csv", checkIfExists: true)
+
+        ch_n_samples = ch_rows.count()
+
+        PREPARE_INPUT_STANDARD(ch_rows, ch_barcode_coordinate_config)
+
+        // input files are split by type (fastq, bam)
+        ch_input_fastq = PREPARE_INPUT_STANDARD.out.fastq
+        ch_input_bam = PREPARE_INPUT_STANDARD.out.bam
+
+        // process fastq samples
+        ch_preprocess_fastq_in = ch_input_fastq.map { sample, path, meta -> [sample, path, meta, meta.barcode] } // add whitelist path to fastq input tuple
+        PREPROCESS_FASTQ(ch_preprocess_fastq_in, ch_flank_seq_config, ch_adapter_seq_config)
+        ALIGNMENT(PREPROCESS_FASTQ.out.fastq, ch_genome, ch_annotation)
+
+        if (!params.bam_only) {
+            // process bam samples
+            ch_bam_files = ALIGNMENT.out.bam.mix(ch_input_bam)
+            BAMBU_PREPARE_ANNOTATION(ch_annotation)
+            BAMBU_CONSTRUCT_READ_CLASS(ch_bam_files, ch_genome, BAMBU_PREPARE_ANNOTATION.out.annotation)
+
+            // collect the read class files for joint clustering
+            ch_rds_files_collect = BAMBU_CONSTRUCT_READ_CLASS.out.rds
+                .map { sample, path, meta -> [sample, path, meta, meta.spatial_metadata] }
+                .collect(flat:false)
+                .map { collected_tup ->
+                    def (samples, paths, metas, spatial_metadatas) = collected_tup.transpose()
+                    def has_spatial = metas.any { meta -> meta.chemistry.startsWith('visium') } // for non-visium samples set the spatial metadata to an empty list (for staging)
+                    [samples, paths, metas, has_spatial ? spatial_metadatas : []]
+                }
+            BAMBU_TRANSCRIPT_DISCOVERY(ch_rds_files_collect, ch_genome, BAMBU_PREPARE_ANNOTATION.out.annotation, ndrArg)
+
+            // cluster the cells first, then pool each cluster's cells for the EM
+            if (params.quantification_mode == 'clusteredEM') {
+                CLUSTERING(BAMBU_TRANSCRIPT_DISCOVERY.out.gene_counts, BAMBU_TRANSCRIPT_DISCOVERY.out.col_data, ch_n_samples)
+                BAMBU_CLUSTER_LEVEL_QUANTIFICATION(CLUSTERING.out.clusters, BAMBU_TRANSCRIPT_DISCOVERY.out.quant_data, BAMBU_TRANSCRIPT_DISCOVERY.out.extended_annotations, ch_genome)
+            } else if (params.quantification_mode == 'EM') {
+                BAMBU_SINGLE_CELL_QUANTIFICATION(BAMBU_TRANSCRIPT_DISCOVERY.out.quant_data, BAMBU_TRANSCRIPT_DISCOVERY.out.extended_annotations, ch_genome)
+            }
+        }
+
+    } else {
+        // restarting pipeline from quantData and manual cluster map
+        ch_clusters   = ch_rows.map { row -> file(row.clusters_path, checkIfExists: true) }
+        ch_quant_data = ch_rows.map { row -> file(row.quant_data_path, checkIfExists: true) }
+
+        BAMBU_CLUSTER_LEVEL_QUANTIFICATION(ch_clusters, ch_quant_data, ch_annotation, ch_genome)
+    }
+}
+
+workflow VISIUM_HD {
+    take:
+    ch_rows: Channel<Map>
+    ch_genome: Path
+    ch_annotation: Path
+    ndr: Float?
+    manual_clustering: Boolean
+
+    main:
+    def requested_bin = String.format('%03dum', params.clustering_bin) // convert clustering_bin specified as an integer into Spaceranger format
+
+    if (!manual_clustering) {
+        def ndrArg = ndr != null ? ndr : 'NULL'
+
+        PREPARE_INPUT_VISIUM_HD(ch_rows)
+        BAMBU_PREPARE_ANNOTATION(ch_annotation)
+        BAMBU_CONSTRUCT_READ_CLASS(PREPARE_INPUT_VISIUM_HD.out.bam, ch_genome, BAMBU_PREPARE_ANNOTATION.out.annotation)
+        BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD(BAMBU_CONSTRUCT_READ_CLASS.out.rds, ch_genome, BAMBU_PREPARE_ANNOTATION.out.annotation, ndrArg, PREPARE_INPUT_VISIUM_HD.out.tissue_positions_002um)
+
+        // perform transcript discovery at 2um first
+        ch_quant_data     = BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD.out.quant_data.first()
+        ch_extended_anno  = BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD.out.extended_annotations.first()
+        ch_unique_002um   = BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD.out.se_unique_002um.first()   // unique counts at 2um resolution
+        ch_barcodes_002um = BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD.out.barcodes_002um.first()    // list of 2um barcodes (same as the column names)
+
+        // map every 2um spot in the SE to its bin, once per resolution; every module below reads this file
+        SPOT_BIN_MAPPINGS(PREPARE_INPUT_VISIUM_HD.out.barcode_mappings, ch_barcodes_002um, PREPARE_INPUT_VISIUM_HD.out.sample_name)
+
+        // pair each bin's tissue positions with its spot mappings, keyed on resolution
+        ch_bins = PREPARE_INPUT_VISIUM_HD.out.tissue_positions_bins.join(SPOT_BIN_MAPPINGS.out.csv) // [resolution, tissue_positions, spot_mappings]
+
+        // aggregate the 2um SEs into each requested bin resolution (e.g., 8um/16um)
+        AGGREGATE_BINS_VISIUM_HD(ch_bins, ch_unique_002um)
+
+        if (params.quantification_mode == 'clusteredEM') {
+            // the 2um counts come straight from transcript discovery, every coarser bin from the aggregated SEs
+            ch_counts_002um = BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD.out.gene_counts_002um
+                .combine(BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD.out.col_data_002um)
+                .map { gene_counts, col_data -> ['002um', gene_counts, col_data] }
+            ch_counts_bins = AGGREGATE_BINS_VISIUM_HD.out.gene_counts.join(AGGREGATE_BINS_VISIUM_HD.out.col_data)
+
+            // perform clustering at the requested resolution only
+            ch_clustering = ch_counts_002um.concat(ch_counts_bins)
+                .filter { resolution, _gene_counts, _col_data -> resolution == requested_bin } // [resolution, gene_counts, col_data]
+            SEURAT_VISIUM_HD(ch_clustering)
+
+            // since quantData is built at 2um, if clustering is performed at the bin-level,
+            // we need to convert clusters from bin-level barcode -> cluster label
+            // to 2um spot/barcode -> cluster label
+            if (requested_bin == '002um') {
+                ch_2um_spot_clusters = SEURAT_VISIUM_HD.out.clusters
+            } else {
+                ch_spot_mappings = SPOT_BIN_MAPPINGS.out.csv
+                    .filter { resolution, _csv -> resolution == requested_bin }
+                    .map { _resolution, csv -> csv }
+                MAP_CLUSTERS_TO_2UM_SPOTS(SEURAT_VISIUM_HD.out.clusters, ch_spot_mappings)
+                ch_2um_spot_clusters = MAP_CLUSTERS_TO_2UM_SPOTS.out.clusters
+            }
+
+            BAMBU_CLUSTER_LEVEL_QUANTIFICATION(ch_2um_spot_clusters, ch_quant_data, ch_extended_anno, ch_genome)
+
+        } else if (params.quantification_mode == 'EM') {
+            // run spot level quantification on all resolutions
+            // at 2um resolution, tissue_positions and spot_mappings are not required
+            ch_resolution_002um = channel.of(['002um', [], []])
+            ch_resolutions      = ch_resolution_002um.concat(ch_bins) // [resolution, tissue_positions, spot_mappings]
+            BAMBU_SPOT_LEVEL_QUANTIFICATION(ch_resolutions, ch_quant_data, ch_extended_anno, ch_genome)
+        }
+
+    } else {
+        // restaring pipeline from quantData and manual cluster mapping
+        ch_clusters   = ch_rows.map { row -> file(row.clusters_path, checkIfExists: true) }
+        ch_quant_data = ch_rows.map { row -> file(row.quant_data_path, checkIfExists: true) }
+
+        if (requested_bin == '002um') {
+            // clusters made at 2um already name the spots bambu quantified
+            ch_2um_spot_clusters = ch_clusters
+        } else {
+            // rebuild the spot mapping file
+            EXTRACT_SPOT_BARCODES(ch_quant_data)
+            ch_sample = EXTRACT_SPOT_BARCODES.out.sample_name.map { name -> name.text.trim() }
+            CONVERT_BARCODE_MAPPINGS(requested_bin, channel.value(params.barcode_mappings))
+            SPOT_BIN_MAPPINGS(CONVERT_BARCODE_MAPPINGS.out.csv, EXTRACT_SPOT_BARCODES.out.barcodes, ch_sample)
+
+            ch_spot_mappings = SPOT_BIN_MAPPINGS.out.csv.map { _resolution, csv -> csv }
+            MAP_CLUSTERS_TO_2UM_SPOTS(ch_clusters, ch_spot_mappings)
+            ch_2um_spot_clusters = MAP_CLUSTERS_TO_2UM_SPOTS.out.clusters
+        }
+
+        BAMBU_CLUSTER_LEVEL_QUANTIFICATION(ch_2um_spot_clusters, ch_quant_data, ch_annotation, ch_genome)
+    }
 }
 
 workflow {
     Validation.validateParams(params, workflow)
 
-    def ndr = params.ndr ?: 'NULL'
-
     // load reference files
     ch_genome     = channel.value(params.genome)
     ch_annotation = channel.value(params.annotation)
@@ -47,55 +224,15 @@ workflow {
         ch_annotation = DECOMPRESS_ANNOTATION.out
     }
 
-    // load config files
-    ch_barcode_coordinate_config = file("${projectDir}/assets/10x_config/barcode_coordinate_config.csv", checkIfExists: true)
-    ch_adapter_seq_config = file("${projectDir}/assets/10x_config/adapter_seq_config.csv", checkIfExists: true)
-    ch_flank_seq_config = file("${projectDir}/assets/10x_config/flank_seq_config.csv", checkIfExists: true)
-
     // parsing samplesheet csv file
     ch_input = channel.of(params.input)
+    ch_rows  = ch_input.splitCsv(header:true, sep:',')
 
-    ch_standard  = ch_input.splitCsv(header:true, sep:',')
-    ch_n_samples = ch_standard.count()
-
-    PREPARE_INPUT_STANDARD(ch_standard, ch_barcode_coordinate_config)
-
-    // input files are split by type (fastq, bam)
-    ch_input_fastq = PREPARE_INPUT_STANDARD.out.fastq
-    ch_input_bam = PREPARE_INPUT_STANDARD.out.bam
-
-    // process fastq samples
-    ch_preprocess_fastq_in = ch_input_fastq.map { sample, path, meta -> [sample, path, meta, meta.barcode] } // add whitelist path to fastq input tuple
-    PREPROCESS_FASTQ(ch_preprocess_fastq_in, ch_flank_seq_config, ch_adapter_seq_config)
-    ALIGNMENT(PREPROCESS_FASTQ.out.fastq, ch_genome, ch_annotation)
-
-    if (!params.bam_only) {
-        // process bam samples
-        ch_bam_files = ALIGNMENT.out.bam.mix(ch_input_bam)
-        BAMBU_PREPARE_ANNOTATION(ch_annotation)
-        BAMBU_CONSTRUCT_READ_CLASS(ch_bam_files, ch_genome, BAMBU_PREPARE_ANNOTATION.out.annotation)
-
-        // collect the read class files for joint clustering
-        ch_rds_files_collect = BAMBU_CONSTRUCT_READ_CLASS.out.rds
-            .map { sample, path, meta -> [sample, path, meta, meta.spatial_metadata] }
-            .collect(flat:false)
-            .map { collected_tup ->
-                def (samples, paths, metas, spatial_metadatas) = collected_tup.transpose()
-                def has_spatial = metas.any { meta -> meta.chemistry.startsWith('visium') } // for non-visium samples set the spatial metadata to an empty list (for staging)
-                [samples, paths, metas, has_spatial ? spatial_metadatas : []]
-            }
-        BAMBU_TRANSCRIPT_DISCOVERY(ch_rds_files_collect, ch_genome, BAMBU_PREPARE_ANNOTATION.out.annotation, ndr)
-
-        if (params.quantification_mode != 'no_quant') {
-            if (params.quantification_mode == 'EM_clusters') {
-                CLUSTERING(BAMBU_TRANSCRIPT_DISCOVERY.out.se_gene_counts, ch_n_samples)
-                ch_clusters = CLUSTERING.out.clusters.map { clusters -> [true, clusters] } // flag to indicate that clustering was performed
-            } else {
-                ch_clusters = channel.value([false, []]) // flag to indicate that clustering was not performed
-            }
-            BAMBU_EM(BAMBU_TRANSCRIPT_DISCOVERY.out.quant_data, BAMBU_TRANSCRIPT_DISCOVERY.out.extended_annotations, ch_clusters, ch_genome)
-        }
+    if (params.visium_hd) {
+        VISIUM_HD(ch_rows, ch_genome, ch_annotation, params.ndr, params.manual_clustering)
+    } else {
+        STANDARD(ch_rows, ch_genome, ch_annotation, params.ndr, params.manual_clustering)
     }
 
     channel.topic('versions').collectFile(name: 'software_versions.yml', storeDir: "${params.output_dir}")
-}
\ No newline at end of file
+}
diff --git a/modules/bambu/EM_quant.nf b/modules/bambu/shared/cluster_level_quantification.nf
similarity index 55%
rename from modules/bambu/EM_quant.nf
rename to modules/bambu/shared/cluster_level_quantification.nf
index b10b28e..408f060 100644
--- a/modules/bambu/EM_quant.nf
+++ b/modules/bambu/shared/cluster_level_quantification.nf
@@ -1,22 +1,20 @@
-process BAMBU_EM{
-    publishDir "$params.output_dir", mode: 'copy', pattern: 'transcript_counts_singlecell'
+process BAMBU_CLUSTER_LEVEL_QUANTIFICATION {
     publishDir "$params.output_dir", mode: 'copy', pattern: 'transcript_counts_clusters'
     publishDir "$params.output_dir", mode: 'copy', pattern: 'gene_counts_clusters'
-    label "r"
+    label "bambu"
     label "low_cpu"
     label "high_mem"
     label "long"
 
     input:
+    path(clusters)
     path(quant_data)
     path(extended_annotation)
-    tuple val(has_clusters), path(clusters)
     path(genome)
 
     output:
-    path ('transcript_counts_singlecell'), optional: true
-    path ('transcript_counts_clusters'), optional: true
-    path ('gene_counts_clusters'), optional: true
+    path ('transcript_counts_clusters')
+    path ('gene_counts_clusters')
     path "versions.yml", topic: 'versions'
 
     script:
@@ -28,29 +26,26 @@ process BAMBU_EM{
 
     extendedAnno <- readRDS("$extended_annotation")
     quantData    <- readRDS("$quant_data")
-    clusters     <- if ("$has_clusters" == "true") readRDS("$clusters") else NULL
-    degBias      <- !is.null(clusters)
+
+    # load clusters if it is a .rds file, else provide path to bambu
+    clusters <- if (endsWith("$clusters", ".rds")) readRDS("$clusters") else "$clusters"
 
     se <- bambu.singlecell(
         reads = quantData,
-        output = if (is.null(clusters)) "EM" else "clusteredEM",
+        output = "clusteredEM",
         annotations = extendedAnno,
         genome = "$genome",
         ncore = $task.cpus,
         verbose = FALSE,
-        opt.em = list(degradationBias = degBias),
+        opt.em = list(degradationBias = TRUE),
         clusters = clusters
     )
 
-    if (is.null(clusters)) {
-        save_counts(se, "transcript_counts_singlecell", "Transcript Expression")
-    } else {
-        save_counts(se, "transcript_counts_clusters", "Transcript Expression")
+    saveCounts(se, "transcript_counts_clusters", "Transcript Expression")
 
-        se_gene <- transcriptToGeneExpression(se)
-        save_counts(se_gene, "gene_counts_clusters")
-    }
+    seGene <- transcriptToGeneExpression(se)
+    saveCounts(seGene, "gene_counts_clusters")
 
     writeLines(c('"${task.process}":', paste0('    R: ', R.Version()\$version.string), paste0('    bambu: ', as.character(packageVersion("bambu")))), "versions.yml")
     """
-}
\ No newline at end of file
+}
diff --git a/modules/bambu/construct_read_class.nf b/modules/bambu/shared/construct_read_class.nf
similarity index 98%
rename from modules/bambu/construct_read_class.nf
rename to modules/bambu/shared/construct_read_class.nf
index 8a7270c..28a3aac 100644
--- a/modules/bambu/construct_read_class.nf
+++ b/modules/bambu/shared/construct_read_class.nf
@@ -1,6 +1,6 @@
 process BAMBU_CONSTRUCT_READ_CLASS{
     publishDir "$params.output_dir/intermediate_R/read_class", mode: 'copy', pattern: '*_read_class.rds', enabled: params.save_intermediates
-    label "r"
+    label "bambu"
     label "low_cpu"
     label "high_mem"
     label "medium"
diff --git a/modules/bambu/prepare_annotation.nf b/modules/bambu/shared/prepare_annotation.nf
similarity index 97%
rename from modules/bambu/prepare_annotation.nf
rename to modules/bambu/shared/prepare_annotation.nf
index 942fe88..ee6fbc3 100644
--- a/modules/bambu/prepare_annotation.nf
+++ b/modules/bambu/shared/prepare_annotation.nf
@@ -1,5 +1,5 @@
 process BAMBU_PREPARE_ANNOTATION{
-    label "r"
+    label "bambu"
     label "low_cpu"
     label "low_mem"
     label "short"
diff --git a/modules/bambu/standard/single_cell_quantification.nf b/modules/bambu/standard/single_cell_quantification.nf
new file mode 100644
index 0000000..d5d212f
--- /dev/null
+++ b/modules/bambu/standard/single_cell_quantification.nf
@@ -0,0 +1,41 @@
+process BAMBU_SINGLE_CELL_QUANTIFICATION {
+    publishDir "$params.output_dir", mode: 'copy', pattern: 'transcript_counts_singlecell'
+    label "bambu"
+    label "low_cpu"
+    label "high_mem"
+    label "long"
+
+    input:
+    path(quant_data)
+    path(extended_annotation)
+    path(genome)
+
+    output:
+    path ('transcript_counts_singlecell')
+    path "versions.yml", topic: 'versions'
+
+    script:
+    """
+    #!/usr/bin/env Rscript
+    if ("$params.bambu_path" == "null") { library("bambu") } else { library("devtools"); load_all("$params.bambu_path") }
+    library(DropletUtils)
+    source(Sys.which("save_counts.R"))
+
+    extendedAnno <- readRDS("$extended_annotation")
+    quantData    <- readRDS("$quant_data")
+
+    se <- bambu.singlecell(
+        reads = quantData,
+        output = "EM",
+        annotations = extendedAnno,
+        genome = "$genome",
+        ncore = $task.cpus,
+        verbose = FALSE,
+        opt.em = list(degradationBias = FALSE)
+    )
+
+    saveCounts(se, "transcript_counts_singlecell", "Transcript Expression")
+
+    writeLines(c('"${task.process}":', paste0('    R: ', R.Version()\$version.string), paste0('    bambu: ', as.character(packageVersion("bambu")))), "versions.yml")
+    """
+}
diff --git a/modules/bambu/transcript_discovery.nf b/modules/bambu/standard/transcript_discovery.nf
similarity index 64%
rename from modules/bambu/transcript_discovery.nf
rename to modules/bambu/standard/transcript_discovery.nf
index 8214740..2a3a34a 100644
--- a/modules/bambu/transcript_discovery.nf
+++ b/modules/bambu/standard/transcript_discovery.nf
@@ -2,8 +2,11 @@ process BAMBU_TRANSCRIPT_DISCOVERY{
     publishDir "$params.output_dir", mode: 'copy', pattern: 'extended_annotations.gtf'
     publishDir "$params.output_dir", mode: 'copy', pattern: 'unique_counts'
     publishDir "$params.output_dir", mode: 'copy', pattern: 'gene_counts'
-    publishDir "$params.output_dir/intermediate_R", mode: 'copy', pattern: '*.rds', enabled: params.save_intermediates
-    label "r"
+    // published unconditionally so a manual clustering run can restart from them
+    publishDir "$params.output_dir/intermediate_R", mode: 'copy', pattern: 'quant_data.rds'
+    publishDir "$params.output_dir/intermediate_R", mode: 'copy', pattern: 'extended_annotations.rds'
+    publishDir "$params.output_dir/intermediate_R", mode: 'copy', pattern: 'col_data.rds', enabled: params.save_intermediates
+    label "bambu"
     label "medium_cpu"
     label "high_mem"
     label "medium"
@@ -17,11 +20,11 @@ process BAMBU_TRANSCRIPT_DISCOVERY{
 	output:
     path ('quant_data.rds'), emit: quant_data
     path ('unique_counts/se_unique_counts.rds')
-    path ('gene_counts/se_gene_counts.rds'), emit: se_gene_counts
 	path ('extended_annotations.rds'), emit: extended_annotations
     path ('extended_annotations.gtf')
     path ('unique_counts')
-    path ('gene_counts')
+    path ('gene_counts'), emit: gene_counts
+    path ('col_data.rds'), emit: col_data
     path "versions.yml", topic: 'versions'
 
 	script:
@@ -29,6 +32,7 @@ process BAMBU_TRANSCRIPT_DISCOVERY{
 	#!/usr/bin/env Rscript
     if ("$params.bambu_path" == "null") { library("bambu") } else { library("devtools"); load_all("$params.bambu_path") }
     library(DropletUtils)
+    source(Sys.which("bambu_discovery.R"))
     source(Sys.which("save_counts.R"))
 
     annotation <- readRDS("$bambu_annotation")
@@ -37,28 +41,26 @@ process BAMBU_TRANSCRIPT_DISCOVERY{
     sampleData <- strsplit("${spatial_metadata_files.join(',')}", ",")[[1]]
     chemistry  <- setNames(strsplit("${meta.collect { m -> m.chemistry }.join(',')}", ",")[[1]], sampleNames)
     technology <- setNames(strsplit("${meta.collect { m -> m.technology }.join(',')}", ",")[[1]], sampleNames)
-
-    # Transcript discovery
-    extendedAnno <- bambu.singlecell(reads = readClassFile, output = "extendedAnnotations",
-    annotations = annotation, genome = "$genome", ncore = $task.cpus, verbose = FALSE, NDR = $ndr)
-    saveRDS(extendedAnno, "extended_annotations.rds")
-    writeToGTF(extendedAnno, "extended_annotations.gtf")
-
-    # Read to transcript assignment
     sampleData <- if (any(startsWith(chemistry, "visium-v"))) sampleData else NULL # Add spatial metadata for visium samples
-    quantData <- bambu.singlecell(reads = readClassFile, output = "quantData",
-    annotations = extendedAnno, genome = "$genome", ncore = $task.cpus, verbose = FALSE, sampleData = sampleData)
-    saveRDS(quantData, "quant_data.rds")
 
-    # Quantification without EM
-    seDiscovery <- bambu.singlecell(reads = quantData, output = "uniqueCounts", annotations = extendedAnno)
+    result <- bambuDiscovery(reads = readClassFile, annotation = annotation, genome = "$genome",
+        ncore = $task.cpus, ndr = $ndr, sampleData = sampleData)
+    saveRDS(result\$extendedAnno, "extended_annotations.rds")
+    writeToGTF(result\$extendedAnno, "extended_annotations.gtf")
+    saveRDS(result\$quantData, "quant_data.rds")
+
+    seDiscovery <- result\$seDiscovery
     colData(seDiscovery)\$chemistry  <- unname(chemistry[colData(seDiscovery)\$sampleName]) # Add chemistry into colData (for subsequent batch correction)
     colData(seDiscovery)\$technology <- unname(technology[colData(seDiscovery)\$sampleName]) # Add technology into colData (for subsequent batch correction)
-    save_counts(seDiscovery, "unique_counts", "Transcript Expression")
+    saveCounts(seDiscovery, "unique_counts", "Transcript Expression")
 
     # Generate gene counts SE from unique counts SE
-    seDiscovery.gene <- transcriptToGeneExpression(seDiscovery)
-    save_counts(seDiscovery.gene, "gene_counts")
+    seDiscoveryGene <- transcriptToGeneExpression(seDiscovery)
+    saveCounts(seDiscoveryGene, "gene_counts")
+
+    # colData is saved to create the Seurat object's metadata
+    geneColData <- as.data.frame(colData(seDiscoveryGene))
+    saveRDS(geneColData, "col_data.rds")
 
     writeLines(c('"${task.process}":', paste0('    R: ', R.Version()\$version.string), paste0('    bambu: ', as.character(packageVersion("bambu")))), "versions.yml")
 	"""
diff --git a/modules/bambu/visium_hd/aggregate_bins.nf b/modules/bambu/visium_hd/aggregate_bins.nf
new file mode 100644
index 0000000..e65c1ff
--- /dev/null
+++ b/modules/bambu/visium_hd/aggregate_bins.nf
@@ -0,0 +1,46 @@
+process AGGREGATE_BINS_VISIUM_HD {
+    publishDir "$params.output_dir/unique_counts", mode: 'copy', pattern: 'unique_counts_*'
+    publishDir "$params.output_dir/gene_counts", mode: 'copy', pattern: 'gene_counts_*'
+    label "bambu"
+    label "low_cpu"
+    label "medium_mem"
+    label "short"
+
+    input:
+    tuple val(resolution), path(tissue_positions), path(spot_mappings)
+    path(se_unique_002um)
+
+    output:
+    tuple val(resolution), path("unique_counts_${resolution}"), emit: unique_counts
+    tuple val(resolution), path("gene_counts_${resolution}"),   emit: gene_counts
+    tuple val(resolution), path("col_data_${resolution}.rds"),  emit: col_data
+    path "versions.yml", topic: 'versions'
+
+    script:
+    """
+    #!/usr/bin/env Rscript
+    if ("$params.bambu_path" == "null") { library("bambu") } else { library("devtools"); load_all("$params.bambu_path") }
+    library(DropletUtils)
+    library(Matrix)
+    library(dplyr)
+    source(Sys.which("save_counts.R"))
+    source(Sys.which("visium_hd_aggregate_resolution.R"))
+
+    unique002um     <- readRDS("$se_unique_002um")
+    tissuePositions <- read.csv("$tissue_positions", stringsAsFactors = FALSE)
+    spotMappings    <- read.csv("$spot_mappings", stringsAsFactors = FALSE)
+
+    # aggregate transcript and gene counts
+    uniqueAgg <- aggregateResolution(unique002um, tissuePositions, spotMappings)
+    geneAgg   <- transcriptToGeneExpression(uniqueAgg)
+
+    saveCounts(uniqueAgg, "unique_counts_$resolution", "Transcript Expression")
+    saveCounts(geneAgg,   "gene_counts_$resolution")
+
+    # colData is saved to create the Seurat object's metadata
+    binColData <- as.data.frame(colData(geneAgg))
+    saveRDS(binColData, "col_data_${resolution}.rds")
+
+    writeLines(c('"${task.process}":', paste0('    R: ', R.Version()\$version.string), paste0('    bambu: ', as.character(packageVersion("bambu")))), "versions.yml")
+    """
+}
diff --git a/modules/bambu/visium_hd/extract_spot_barcodes.nf b/modules/bambu/visium_hd/extract_spot_barcodes.nf
new file mode 100644
index 0000000..2cb953d
--- /dev/null
+++ b/modules/bambu/visium_hd/extract_spot_barcodes.nf
@@ -0,0 +1,29 @@
+process EXTRACT_SPOT_BARCODES {
+    label "bambu"
+    label "low_cpu"
+    label "medium_mem"
+    label "short"
+
+    input:
+    path(quant_data)
+
+    output:
+    path ('barcodes.tsv.gz'), emit: barcodes
+    path ('sample_name.txt'), emit: sample_name
+    path "versions.yml", topic: 'versions'
+
+    script:
+    """
+    #!/usr/bin/env Rscript
+    if ("$params.bambu_path" == "null") { library("bambu") } else { library("devtools"); load_all("$params.bambu_path") }
+
+    # a Visium HD run holds a single sample, whose ids name the 2um spots bambu quantified
+    quantData  <- readRDS("$quant_data")
+    sampleData <- quantData[[1]]@sampleData
+
+    write.table(sampleData\$id, gzfile("barcodes.tsv.gz"), quote = FALSE, row.names = FALSE, col.names = FALSE)
+    writeLines(unique(sampleData\$sampleName), "sample_name.txt")
+
+    writeLines(c('"${task.process}":', paste0('    R: ', R.Version()\$version.string), paste0('    bambu: ', as.character(packageVersion("bambu")))), "versions.yml")
+    """
+}
diff --git a/modules/bambu/visium_hd/map_clusters_to_2um_spots.nf b/modules/bambu/visium_hd/map_clusters_to_2um_spots.nf
new file mode 100644
index 0000000..efd3ab9
--- /dev/null
+++ b/modules/bambu/visium_hd/map_clusters_to_2um_spots.nf
@@ -0,0 +1,33 @@
+process MAP_CLUSTERS_TO_2UM_SPOTS {
+    label "bambu"
+    label "low_cpu"
+    label "medium_mem"
+    label "short"
+
+    input:
+    path(clusters)
+    path(spot_mappings)
+
+    output:
+    path ('spot_clusters.rds'), emit: clusters
+    path "versions.yml", topic: 'versions'
+
+    script:
+    """
+    #!/usr/bin/env Rscript
+    source(Sys.which("map_spots_to_clusters.R"))
+
+    binToClustersMap <- readClusters("$clusters")
+    spotMappings     <- read.csv("$spot_mappings", stringsAsFactors = FALSE)
+
+    # bambu quantifies at 2um, so map the bin level labels down to the 2um spots
+    spotToClustersMap <- mapSpotsToClusters(binToClustersMap, spotMappings)
+    if (length(spotToClustersMap) == 0) {
+        stop("no cluster id matches a bin in $spot_mappings -- check that the clusters were made at the ${params.clustering_bin}um resolution")
+    }
+
+    saveRDS(spotToClustersMap, "spot_clusters.rds")
+
+    writeLines(c('"${task.process}":', paste0('    R: ', R.Version()\$version.string)), "versions.yml")
+    """
+}
diff --git a/modules/bambu/visium_hd/spot_level_quantification.nf b/modules/bambu/visium_hd/spot_level_quantification.nf
new file mode 100644
index 0000000..af19419
--- /dev/null
+++ b/modules/bambu/visium_hd/spot_level_quantification.nf
@@ -0,0 +1,63 @@
+process BAMBU_SPOT_LEVEL_QUANTIFICATION {
+    publishDir "$params.output_dir/transcript_counts", mode: 'copy', pattern: 'transcript_counts_*'
+    label "bambu"
+    label "low_cpu"
+    label "high_mem"
+    label "long"
+
+    input:
+    tuple val(resolution), path(tissue_positions), path(spot_mappings)
+    path(quant_data)
+    path(extended_annotation)
+    path(genome)
+
+    output:
+    path ("transcript_counts_${resolution}")
+    path "versions.yml", topic: 'versions'
+
+    script:
+    """
+    #!/usr/bin/env Rscript
+    if ("$params.bambu_path" == "null") { library("bambu") } else { library("devtools"); load_all("$params.bambu_path") }
+    library(DropletUtils)
+    library(dplyr)
+    source(Sys.which("save_counts.R"))
+
+    extendedAnno <- readRDS("$extended_annotation")
+    quantData    <- readRDS("$quant_data")
+
+    # For bin-level quantification, we need to use the clusters argument to collapse the quantData object
+    # to the desired resolution, since quantData is generated at 2um resolution.
+    if ("$resolution" == "002um") {
+        clusters <- NULL
+    } else {
+        # bambu prefixes the sample name onto each cluster label, so pass the bare bin
+        spotMappings <- read.csv("$spot_mappings", stringsAsFactors = FALSE)
+        clusters     <- setNames(spotMappings\$bin, spotMappings\$sample_barcode)
+    }
+
+    se <- bambu.singlecell(
+        reads = quantData,
+        output = if (is.null(clusters)) "EM" else "clusteredEM",
+        annotations = extendedAnno,
+        genome = "$genome",
+        ncore = $task.cpus,
+        verbose = FALSE,
+        opt.em = list(degradationBias = FALSE),
+        clusters = clusters
+    )
+
+    # At the bin-level resolutions bambu builds a fresh colData, so we need to reattach the bin metadata.
+    if (!is.null(clusters)) {
+        tissuePositions <- read.csv("$tissue_positions", stringsAsFactors = FALSE)
+        binColData      <- as.data.frame(colData(se)) %>%
+            left_join(tissuePositions, by = c("cluster" = "barcode")) %>%
+            rename(barcode = cluster)
+        colData(se) <- DataFrame(binColData, row.names = colnames(se))
+    }
+
+    saveCounts(se, "transcript_counts_$resolution", "Transcript Expression")
+
+    writeLines(c('"${task.process}":', paste0('    R: ', R.Version()\$version.string), paste0('    bambu: ', as.character(packageVersion("bambu")))), "versions.yml")
+    """
+}
diff --git a/modules/bambu/visium_hd/transcript_discovery.nf b/modules/bambu/visium_hd/transcript_discovery.nf
new file mode 100644
index 0000000..ffe8e3b
--- /dev/null
+++ b/modules/bambu/visium_hd/transcript_discovery.nf
@@ -0,0 +1,62 @@
+process BAMBU_TRANSCRIPT_DISCOVERY_VISIUM_HD {
+    publishDir "$params.output_dir", mode: 'copy', pattern: 'extended_annotations.gtf'
+    publishDir "$params.output_dir/unique_counts", mode: 'copy', pattern: 'unique_counts_002um'
+    publishDir "$params.output_dir/gene_counts", mode: 'copy', pattern: 'gene_counts_002um'
+    // published unconditionally so a manual clustering run can restart from them
+    publishDir "$params.output_dir/intermediate_R", mode: 'copy', pattern: 'quant_data.rds'
+    publishDir "$params.output_dir/intermediate_R", mode: 'copy', pattern: 'extended_annotations.rds'
+    publishDir "$params.output_dir/intermediate_R", mode: 'copy', pattern: 'col_data_002um.rds', enabled: params.save_intermediates
+    label "bambu"
+    label "medium_cpu"
+    label "high_mem"
+    label "medium"
+
+    input:
+    tuple val(sample), path(rds_files), val(meta)
+    path(genome)
+    path(bambu_annotation)
+    val(ndr)
+    path(tissue_positions)
+
+    output:
+    path ('quant_data.rds'), emit: quant_data
+    path ('extended_annotations.rds'), emit: extended_annotations
+    path ('extended_annotations.gtf')
+    path ('unique_counts_002um')
+    path ('gene_counts_002um'), emit: gene_counts_002um
+    path ('unique_counts_002um/se_unique_counts_002um.rds'), emit: se_unique_002um
+    path ('unique_counts_002um/barcodes.tsv.gz'), emit: barcodes_002um
+    path ('col_data_002um.rds'), emit: col_data_002um
+    path "versions.yml", topic: 'versions'
+
+    script:
+    """
+    #!/usr/bin/env Rscript
+    if ("$params.bambu_path" == "null") { library("bambu") } else { library("devtools"); load_all("$params.bambu_path") }
+    library(DropletUtils)
+    source(Sys.which("bambu_discovery.R"))
+    source(Sys.which("save_counts.R"))
+
+    annotation <- readRDS("$bambu_annotation")
+    readClassFile <- setNames("$rds_files", "$sample")
+
+    # the 2um tissue_positions is passed as sampleData so bambu joins the per-spot spatial metadata into colData
+    result <- bambuDiscovery(reads = readClassFile, annotation = annotation, genome = "$genome",
+        ncore = $task.cpus, ndr = $ndr, sampleData = "$tissue_positions")
+    saveRDS(result\$extendedAnno, "extended_annotations.rds")
+    writeToGTF(result\$extendedAnno, "extended_annotations.gtf")
+    saveRDS(result\$quantData, "quant_data.rds")
+
+    unique002um <- result\$seDiscovery
+    gene002um   <- transcriptToGeneExpression(unique002um)
+
+    saveCounts(unique002um, "unique_counts_002um", "Transcript Expression")
+    saveCounts(gene002um,   "gene_counts_002um")
+
+    # colData is saved to create the Seurat object's metadata
+    spotColData <- as.data.frame(colData(gene002um))
+    saveRDS(spotColData, "col_data_002um.rds")
+
+    writeLines(c('"${task.process}":', paste0('    R: ', R.Version()\$version.string), paste0('    bambu: ', as.character(packageVersion("bambu")))), "versions.yml")
+    """
+}
diff --git a/modules/prepare_input_standard/extract_barcodes.nf b/modules/prepare_input/standard/extract_barcodes.nf
similarity index 100%
rename from modules/prepare_input_standard/extract_barcodes.nf
rename to modules/prepare_input/standard/extract_barcodes.nf
diff --git a/modules/prepare_input_standard/extract_spatial_coordinates.nf b/modules/prepare_input/standard/extract_spatial_coordinates.nf
similarity index 100%
rename from modules/prepare_input_standard/extract_spatial_coordinates.nf
rename to modules/prepare_input/standard/extract_spatial_coordinates.nf
diff --git a/modules/prepare_input/visium_hd/convert_barcode_mappings.nf b/modules/prepare_input/visium_hd/convert_barcode_mappings.nf
new file mode 100644
index 0000000..422ba04
--- /dev/null
+++ b/modules/prepare_input/visium_hd/convert_barcode_mappings.nf
@@ -0,0 +1,26 @@
+process CONVERT_BARCODE_MAPPINGS {
+    publishDir "$params.output_dir/intermediate_visium_hd", mode: 'copy', pattern: '*.csv', enabled: params.save_intermediates
+    label "pyarrow_pandas"
+    label "low_cpu"
+    label "low_mem"
+    label "short"
+
+    input:
+    val(resolution)
+    path(barcode_mappings)
+
+    output:
+    tuple val(resolution), path("barcode_mappings_${resolution}.csv"), emit: csv
+    path "versions.yml", topic: 'versions'
+
+    script:
+    """
+    visium_hd_convert_barcode_mappings.py $barcode_mappings ${resolution} barcode_mappings_${resolution}.csv
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        python: \$(python3 --version 2>&1)
+        pandas: \$(python3 -c 'import pandas; print(pandas.__version__)')
+    END_VERSIONS
+    """
+}
diff --git a/modules/prepare_input/visium_hd/convert_tissue_positions.nf b/modules/prepare_input/visium_hd/convert_tissue_positions.nf
new file mode 100644
index 0000000..580df75
--- /dev/null
+++ b/modules/prepare_input/visium_hd/convert_tissue_positions.nf
@@ -0,0 +1,27 @@
+process CONVERT_TISSUE_POSITIONS {
+    publishDir "$params.output_dir/intermediate_visium_hd", mode: 'copy', pattern: '*.csv', enabled: params.save_intermediates
+    publishDir "$params.output_dir/intermediate_visium_hd", mode: 'copy', pattern: 'barcodes_in_tissue.txt', enabled: params.save_intermediates
+    label "pyarrow_pandas"
+    label "low_cpu"
+    label "low_mem"
+    label "short"
+
+    input:
+    tuple val(resolution), path(tissue_positions)
+
+    output:
+    tuple val(resolution), path("tissue_positions_${resolution}.csv"), emit: csv
+    path("barcodes_in_tissue.txt"), optional: true, emit: barcodes
+    path "versions.yml", topic: 'versions'
+
+    script:
+    """
+    visium_hd_convert_tissue_positions.py $tissue_positions ${resolution} tissue_positions_${resolution}.csv barcodes_in_tissue.txt
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        python: \$(python3 --version 2>&1)
+        pandas: \$(python3 -c 'import pandas; print(pandas.__version__)')
+    END_VERSIONS
+    """
+}
diff --git a/modules/prepare_input/visium_hd/filter_barcoded_bam.nf b/modules/prepare_input/visium_hd/filter_barcoded_bam.nf
new file mode 100644
index 0000000..9b9720b
--- /dev/null
+++ b/modules/prepare_input/visium_hd/filter_barcoded_bam.nf
@@ -0,0 +1,27 @@
+process FILTER_BARCODED_BAM {
+    publishDir "$params.output_dir/intermediate_bam", mode: 'copy', pattern: '*_filtered.bam*', enabled: params.save_intermediates
+    label "minimap2_samtools"
+    label "medium_cpu"
+    label "medium_mem"
+    label "medium"
+
+    input:
+    tuple val(sample), path(bam), val(meta)
+    path(barcodes)
+
+    output:
+    tuple val(sample), path("${sample}_filtered.bam"), val(meta), emit: bam
+    path("${sample}_filtered.bam.bai")
+    path "versions.yml", topic: 'versions'
+
+    script:
+    """
+    samtools view -@ $task.cpus -D CB:$barcodes -o ${sample}_filtered.bam $bam
+    samtools index -@ $task.cpus ${sample}_filtered.bam
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        samtools: \$(samtools --version 2>&1 | head -1)
+    END_VERSIONS
+    """
+}
diff --git a/modules/prepare_input/visium_hd/spot_bin_mappings.nf b/modules/prepare_input/visium_hd/spot_bin_mappings.nf
new file mode 100644
index 0000000..8b0a9ac
--- /dev/null
+++ b/modules/prepare_input/visium_hd/spot_bin_mappings.nf
@@ -0,0 +1,27 @@
+process SPOT_BIN_MAPPINGS {
+    publishDir "$params.output_dir/intermediate_visium_hd", mode: 'copy', pattern: '*.csv.gz', enabled: params.save_intermediates
+    label "pyarrow_pandas"
+    label "low_cpu"
+    label "low_mem"
+    label "short"
+
+    input:
+    tuple val(resolution), path(barcode_mappings)
+    path(barcodes_002um)
+    val(sample)
+
+    output:
+    tuple val(resolution), path("spot_bin_mappings_${resolution}.csv.gz"), emit: csv
+    path "versions.yml", topic: 'versions'
+
+    script:
+    """
+    visium_hd_spot_bin_mappings.py $barcode_mappings $barcodes_002um $sample spot_bin_mappings_${resolution}.csv.gz
+
+    cat <<-END_VERSIONS > versions.yml
+    "${task.process}":
+        python: \$(python3 --version 2>&1)
+        pandas: \$(python3 -c 'import pandas; print(pandas.__version__)')
+    END_VERSIONS
+    """
+}
diff --git a/modules/seurat/multi_sample.nf b/modules/seurat/standard/multi_sample_clustering.nf
similarity index 60%
rename from modules/seurat/multi_sample.nf
rename to modules/seurat/standard/multi_sample_clustering.nf
index 2316cfb..b0883f9 100644
--- a/modules/seurat/multi_sample.nf
+++ b/modules/seurat/standard/multi_sample_clustering.nf
@@ -1,13 +1,13 @@
 process SEURAT_MULTI_SAMPLE {
     publishDir "$params.output_dir", mode: 'copy', pattern: 'seurat_obj.rds'
     publishDir "$params.output_dir/intermediate_R", mode: 'copy', pattern: 'clusters.rds', enabled: params.save_intermediates
-    label "r"
+    label "seurat"
     label "medium_cpu"
     label "high_mem"
     label "medium"
 
     input:
-    path(se)
+    tuple path(gene_counts), path(col_data)
 
     output:
     path ('clusters.rds'), emit: clusters
@@ -17,23 +17,18 @@ process SEURAT_MULTI_SAMPLE {
     script:
     """
     #!/usr/bin/env Rscript
-    library(SummarizedExperiment)
     library(Seurat)
+    source(Sys.which("create_seurat_object.R"))
 
-    # Extract gene count matrix and colData metadata
-    se     <- readRDS("$se")
-    counts <- assays(se)\$counts
-    dim    <- $params.seurat_dim_multi
+    metadata <- readRDS("$col_data")
+    dim      <- $params.seurat_dim_multi
 
-    # Create Seurat object and append metadata
-    cellMix <- CreateSeuratObject(counts = counts, project = "cellMix", min.cells = 1)
-    cellMix\$sample     <- setNames(colData(se)\$sampleName,  colnames(se))
-    cellMix\$chemistry  <- setNames(colData(se)\$chemistry,   colnames(se))
-    cellMix\$technology <- setNames(colData(se)\$technology,  colnames(se))
-    cellMix\$orig.ident <- cellMix\$sample
+    # Create Seurat object with bambu's colData as its metadata
+    cellMix <- createSeuratObject("$gene_counts", metadata = metadata, minCells = 1, project = "cellMix")
+    cellMix\$orig.ident <- cellMix\$sampleName
 
     # scRNA-seq multi-sample integration using Harmony adapted from https://satijalab.org/seurat/articles/seurat5_integration
-    cellMix[["RNA"]] <- split(cellMix[["RNA"]], f = cellMix\$sample)
+    cellMix[["RNA"]] <- split(cellMix[["RNA"]], f = cellMix\$sampleName)
     cellMix <- NormalizeData(cellMix)
     cellMix <- FindVariableFeatures(cellMix)
     cellMix <- ScaleData(cellMix)
@@ -50,7 +45,7 @@ process SEURAT_MULTI_SAMPLE {
 
     dim <- min(dim, ncol(cellMix[["harmony"]]))
     cellMix <- FindNeighbors(cellMix, reduction = "harmony", dims = 1:dim)
-    cellMix <- FindClusters(cellMix, resolution = $params.resolution, cluster.name = "harmony_clusters")
+    cellMix <- FindClusters(cellMix, resolution = $params.seurat_resolution, cluster.name = "harmony_clusters")
     saveRDS(cellMix, "seurat_obj.rds")
 
     clusters <- setNames(paste0("cluster_", cellMix\$harmony_clusters), names(cellMix\$harmony_clusters))
@@ -58,9 +53,9 @@ process SEURAT_MULTI_SAMPLE {
 
     writeLines(c(
         '"${task.process}":',
-        paste0('    R: ',                   R.Version()\$version.string),
-        paste0('    seurat: ',              as.character(packageVersion("Seurat"))),
-        paste0('    SummarizedExperiment: ', as.character(packageVersion("SummarizedExperiment")))
+        paste0('    R: ',       R.Version()\$version.string),
+        paste0('    seurat: ',  as.character(packageVersion("Seurat"))),
+        paste0('    harmony: ', as.character(packageVersion("harmony")))
     ), "versions.yml")
     """
 }
diff --git a/modules/seurat/single_sample.nf b/modules/seurat/standard/single_sample_clustering.nf
similarity index 67%
rename from modules/seurat/single_sample.nf
rename to modules/seurat/standard/single_sample_clustering.nf
index 283374c..406bc82 100644
--- a/modules/seurat/single_sample.nf
+++ b/modules/seurat/standard/single_sample_clustering.nf
@@ -1,13 +1,13 @@
 process SEURAT_SINGLE_SAMPLE {
     publishDir "$params.output_dir", mode: 'copy', pattern: 'seurat_obj.rds'
     publishDir "$params.output_dir/intermediate_R", mode: 'copy', pattern: 'clusters.rds', enabled: params.save_intermediates
-    label "r"
+    label "seurat"
     label "medium_cpu"
     label "medium_mem"
     label "medium"
 
     input:
-    path(se)
+    tuple path(gene_counts), path(col_data)
 
     output:
     path ('clusters.rds'), emit: clusters
@@ -17,23 +17,24 @@ process SEURAT_SINGLE_SAMPLE {
     script:
     """
     #!/usr/bin/env Rscript
-    library(SummarizedExperiment)
     library(Seurat)
+    source(Sys.which("create_seurat_object.R"))
 
-    se     <- readRDS("$se")
-    counts <- assays(se)\$counts
-    dim    <- $params.seurat_dim_single
+    metadata <- readRDS("$col_data")
+    dim      <- $params.seurat_dim_single
+
+    # Create Seurat object with bambu's colData as its metadata
+    cellMix <- createSeuratObject("$gene_counts", metadata = metadata, minCells = 1, project = "cellMix")
 
     # Single sample scRNA-seq clustering adapted from https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
-    cellMix <- CreateSeuratObject(counts = counts, project = "cellMix", min.cells = 1)
     cellMix <- NormalizeData(cellMix, normalization.method = "LogNormalize", scale.factor = 10000)
     cellMix <- FindVariableFeatures(cellMix, selection.method = "vst", nfeatures = 2500)
     cellMix <- ScaleData(cellMix)
-    npcs    <- ifelse(ncol(counts) > 50, 50, ncol(counts) - 1)
+    npcs    <- ifelse(ncol(cellMix) > 50, 50, ncol(cellMix) - 1)
     cellMix <- RunPCA(cellMix, features = VariableFeatures(object = cellMix), npcs = npcs)
     dim     <- ifelse(dim >= dim(cellMix@reductions\$pca)[2], dim(cellMix@reductions\$pca)[2], dim)
     cellMix <- FindNeighbors(cellMix, dims = 1:dim)
-    cellMix <- FindClusters(cellMix, resolution = $params.resolution, cluster.name = "clusters")
+    cellMix <- FindClusters(cellMix, resolution = $params.seurat_resolution, cluster.name = "clusters")
     
     saveRDS(cellMix, "seurat_obj.rds")
 
@@ -42,9 +43,8 @@ process SEURAT_SINGLE_SAMPLE {
     
     writeLines(c(
         '"${task.process}":',
-        paste0('    R: ',                   R.Version()\$version.string),
-        paste0('    seurat: ',              as.character(packageVersion("Seurat"))),
-        paste0('    SummarizedExperiment: ', as.character(packageVersion("SummarizedExperiment")))
+        paste0('    R: ',      R.Version()\$version.string),
+        paste0('    seurat: ', as.character(packageVersion("Seurat")))
     ), "versions.yml")
     """
 }
diff --git a/modules/seurat/visium_hd/clustering.nf b/modules/seurat/visium_hd/clustering.nf
new file mode 100644
index 0000000..6427759
--- /dev/null
+++ b/modules/seurat/visium_hd/clustering.nf
@@ -0,0 +1,55 @@
+process SEURAT_VISIUM_HD {
+    publishDir "$params.output_dir", mode: 'copy', pattern: 'seurat_obj.rds'
+    label "seurat"
+    label "medium_cpu"
+    label "high_mem"
+    label "long"
+
+    input:
+    tuple val(resolution), path(gene_counts), path(col_data)
+
+    output:
+    path ("clusters.rds"), emit: clusters
+    path ("seurat_obj.rds")
+    path "versions.yml", topic: 'versions'
+
+    script:
+    """
+    #!/usr/bin/env Rscript
+    library(Seurat)
+    source(Sys.which("create_seurat_object.R"))
+    source(Sys.which("visium_hd_cluster.R"))
+
+    banksy <- as.logical("$params.banksy")
+    assay  <- "Spatial.$resolution"
+
+    # Create Seurat object with bambu's colData as its metadata
+    colData      <- readRDS("$col_data")
+    binMetadata  <- colData[, colnames(colData) != "barcodes"] # drop barcodes column (list) to prevent downstream issues with Seurat
+    object       <- createSeuratObject("$gene_counts", assay = assay, metadata = binMetadata)
+
+    DefaultAssay(object) <- assay
+    object <- NormalizeData(object, verbose = FALSE)
+
+    # Cluster spatially with Banksy, or on gene expression alone
+    if (banksy) {
+        library(SeuratWrappers)
+        library(Banksy)
+        object <- clusterBanksy(object, assay, $params.banksy_lambda, $params.banksy_k_geom, $params.seurat_resolution)
+    } else {
+        object <- clusterExpression(object, assay, $params.seurat_resolution)
+    }
+
+    # Save the cluster assignment at the clustered resolution, plus the Seurat object
+    clusters <- setNames(paste0("cluster_", object\$clusters), names(object\$clusters))
+    saveRDS(clusters, "clusters.rds")
+    saveRDS(object, "seurat_obj.rds")
+
+    writeLines(c(
+        '"${task.process}":',
+        paste0('    R: ',      R.Version()\$version.string),
+        paste0('    seurat: ', as.character(packageVersion("Seurat"))),
+        paste0('    banksy: ', as.character(packageVersion("Banksy")))
+    ), "versions.yml")
+    """
+}
diff --git a/nextflow.config b/nextflow.config
index 5c8bf0a..e6a0a32 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -12,33 +12,57 @@ params {
     input                   = null                          // Path to samplesheet .csv file
     genome                  = null                          // Path to .fa or .fasta file
     annotation              = null                          // Path to .gtf or .gff file
-    
+
     // Optional: Output directory
     output_dir              = "output"                      // Path to output directory
 
     /*
     Optional: Samplesheet settings (Non Visium HD samples only)
     Note: Use this if all samples share the same chemistry/technology
-    */ 
-    chemistry               = null                          // Examples: "10x3v2", "10x3v3", "10x5v2", "visium-v1" 
-    technology              = null                          // Options: "ONT", "PacBio" 
-    
+    */
+    chemistry               = null                          // Examples: "10x3v2", "10x3v3", "10x5v2", "visium-v1"
+    technology              = null                          // Options: "ONT", "PacBio"
+
     // Optional: Stop after alignment and save BAM files only
     bam_only                = false                         // boolean
 
     // Optional: Q-score filtering
-    qscore_filtering        = true                          // boolean  
+    qscore_filtering        = true                          // boolean
 
     // Optional: Bambu parameters
-    ndr                     = null                          // null or float 
+    ndr                     = null                          // null or float
     deduplicate_umis        = true                          // boolean
 
     // Optional: Quantification mode
-    quantification_mode     = "EM_clusters"                 // Options: "no_quant", "EM", "EM_clusters"
+    quantification_mode     = "clusteredEM"                 // Options: "no_quant", "EM", "clusteredEM"
 
     // Optional: Seurat clustering
-    resolution              = 0.8                           // float
-  
+    seurat_resolution       = 0.8                           // float
+
+    /*
+    Optional: Manual clustering
+    Note: Set --manual_clustering true to quantify from clusters generated outside the pipeline.
+    The samplesheet takes 'clusters_path' and 'quant_data_path' columns instead of the usual ones,
+    and --annotation takes the extended_annotations.rds written by the earlier run.
+    */
+    manual_clustering       = false                         // boolean
+
+    /*
+    Optional: Visium HD
+    Note: Set --visium_hd true for a single-sample run starting from a pre-aligned, barcode-tagged BAM file.
+    The samplesheet only needs 'sample' and 'path' columns — 'chemistry' and 'technology' are not used.
+    Counts are produced at 2um and at every resolution listed in the bins samplesheet.
+    */
+    visium_hd                = false                        // boolean
+    bins                     = null                         // CSV samplesheet of resolutions (columns: resolution,tissue_positions); must include a 2 row (native base) plus any bins; required when --visium_hd true
+    barcode_mappings         = null                         // Path to barcode_mappings.parquet (required when --visium_hd true)
+
+    // Optional: Visium HD clustering
+    clustering_bin           = 8                            // Integer bin size to cluster at; must be listed in the bins samplesheet
+    banksy                   = true                         // boolean; spatially aware clustering
+    banksy_lambda            = 0.8                          // float; 0.8 segments tissue domains, 0.2 gives spatially informed cell types
+    banksy_k_geom            = 50                           // integer; spatial neighbours per bin
+
 }
 
 includeConfig 'conf/base.config'
@@ -50,7 +74,7 @@ profiles {
 	singularity {
         singularity.enabled           = true
         singularity.autoMounts        = true
-        docker.enabled                = false  
+        docker.enabled                = false
 	    singularity.runOptions        = "--bind ${env('PWD')}"
         singularity.envWhitelist      = "JAVA_HOME"
     }
@@ -88,7 +112,7 @@ profiles {
 
 includeConfig 'conf/smoke_test.config'
 
-// Output an html timeline report 
+// Output an html timeline report
 timeline { enabled = true; file = "${params.output_dir}/pipeline_info/execution_timeline.html" }
 
 // Output resource and runtime reports for a workflow run
@@ -98,4 +122,4 @@ report   { enabled = true; file = "${params.output_dir}/pipeline_info/execution_
 trace    { enabled = true; file = "${params.output_dir}/pipeline_info/execution_trace.txt" }
 
 // Produce a workflow diagram
-dag      { enabled = true; file = "${params.output_dir}/pipeline_info/pipeline_dag.svg" }
\ No newline at end of file
+dag      { enabled = true; file = "${params.output_dir}/pipeline_info/pipeline_dag.svg" }
diff --git a/subworkflows/clustering.nf b/subworkflows/clustering.nf
deleted file mode 100644
index 08a3979..0000000
--- a/subworkflows/clustering.nf
+++ /dev/null
@@ -1,22 +0,0 @@
-include { SEURAT_SINGLE_SAMPLE } from '../modules/seurat/single_sample.nf'
-include { SEURAT_MULTI_SAMPLE  } from '../modules/seurat/multi_sample.nf'
-
-workflow CLUSTERING {
-    take:
-    ch_se_gene_counts
-    ch_n_samples
-
-    main:
-    ch_se_branched = ch_se_gene_counts
-        .combine(ch_n_samples)
-        .branch { _se, n ->
-            single: n == 1
-            multi:  n > 1
-        }
-
-    SEURAT_SINGLE_SAMPLE(ch_se_branched.single.map { se, _n -> se })
-    SEURAT_MULTI_SAMPLE(ch_se_branched.multi.map { se, _n -> se })
-
-    emit:
-    clusters = SEURAT_SINGLE_SAMPLE.out.clusters.mix(SEURAT_MULTI_SAMPLE.out.clusters)
-}
diff --git a/subworkflows/clustering_standard.nf b/subworkflows/clustering_standard.nf
new file mode 100644
index 0000000..4d9168f
--- /dev/null
+++ b/subworkflows/clustering_standard.nf
@@ -0,0 +1,24 @@
+include { SEURAT_SINGLE_SAMPLE } from '../modules/seurat/standard/single_sample_clustering.nf'
+include { SEURAT_MULTI_SAMPLE  } from '../modules/seurat/standard/multi_sample_clustering.nf'
+
+workflow CLUSTERING {
+    take:
+    ch_gene_counts
+    ch_col_data
+    ch_n_samples
+
+    main:
+    ch_branched = ch_gene_counts
+        .combine(ch_col_data)
+        .combine(ch_n_samples)
+        .branch { _gene_counts, _col_data, n ->
+            single: n == 1
+            multi:  n > 1
+        }
+
+    SEURAT_SINGLE_SAMPLE(ch_branched.single.map { gene_counts, col_data, _n -> [gene_counts, col_data] })
+    SEURAT_MULTI_SAMPLE(ch_branched.multi.map { gene_counts, col_data, _n -> [gene_counts, col_data] })
+
+    emit:
+    clusters = SEURAT_SINGLE_SAMPLE.out.clusters.mix(SEURAT_MULTI_SAMPLE.out.clusters)
+}
diff --git a/subworkflows/prepare_input_standard.nf b/subworkflows/prepare_input_standard.nf
index 4275965..8e01e07 100644
--- a/subworkflows/prepare_input_standard.nf
+++ b/subworkflows/prepare_input_standard.nf
@@ -1,5 +1,5 @@
-include { EXTRACT_10X_BARCODES } from '../modules/prepare_input_standard/extract_barcodes.nf'
-include { EXTRACT_10X_SPATIAL_COORDINATES } from '../modules/prepare_input_standard/extract_spatial_coordinates.nf'
+include { EXTRACT_10X_BARCODES } from '../modules/prepare_input/standard/extract_barcodes.nf'
+include { EXTRACT_10X_SPATIAL_COORDINATES } from '../modules/prepare_input/standard/extract_spatial_coordinates.nf'
 
 workflow PREPARE_INPUT_STANDARD {
     take:
diff --git a/subworkflows/prepare_input_visium_hd.nf b/subworkflows/prepare_input_visium_hd.nf
new file mode 100644
index 0000000..4866c7c
--- /dev/null
+++ b/subworkflows/prepare_input_visium_hd.nf
@@ -0,0 +1,53 @@
+include { FILTER_BARCODED_BAM }      from '../modules/prepare_input/visium_hd/filter_barcoded_bam.nf'
+include { CONVERT_BARCODE_MAPPINGS } from '../modules/prepare_input/visium_hd/convert_barcode_mappings.nf'
+include { CONVERT_TISSUE_POSITIONS } from '../modules/prepare_input/visium_hd/convert_tissue_positions.nf'
+
+workflow PREPARE_INPUT_VISIUM_HD {
+    take:
+    ch_rows  // raw samplesheet rows
+
+    main:
+    // Visium HD: single sample, starting from a pre-aligned, barcode-tagged BAM file
+    ch_rows.collect(flat: false).map { rows -> Validation.validateVisiumHDRows(rows) }
+
+    ch_sample = ch_rows.map { row ->
+        def sample_path = file(row.path, checkIfExists: true)
+        def meta = [chemistry: 'visium-hd', technology: 'NA']
+        [row.sample, sample_path, meta]
+    }
+
+    // parse the bins samplesheet into one [resolution, tissue_positions] tuple per row (includes mandatory 2um base)
+    ch_resolutions = channel.of(params.bins)
+        .splitCsv(header: true, sep: ',')
+        .map { row ->
+            // Convert integer values in the resolution column of the CSV into Spaceranger resolution, e.g. 8 -> "008um"
+            def resolution = String.format('%03dum', row.resolution as Integer)
+            def tissue_positions = file(row.tissue_positions, checkIfExists: true)
+            [resolution, tissue_positions]
+        }
+    ch_barcode_mappings = channel.value(params.barcode_mappings)
+
+    // convert every resolution's tissue positions parquet to CSV; the 2um task also extracts the in-tissue barcode list
+    CONVERT_TISSUE_POSITIONS(ch_resolutions)
+
+    // the 2um in-tissue barcodes are used to filter out-of-tissue reads from the BAM file
+    FILTER_BARCODED_BAM(ch_sample, CONVERT_TISSUE_POSITIONS.out.barcodes)
+
+    // convert the barcode mappings parquet to one 'barcode,bin' CSV per bin resolution level (e.g., 8um/16um)
+    ch_bin_resolutions = ch_resolutions
+        .map { resolution, _tissue_positions -> resolution }
+        .filter { resolution -> resolution != '002um' } 
+    CONVERT_BARCODE_MAPPINGS(ch_bin_resolutions, ch_barcode_mappings)
+
+    // extract the 2um tissue position separately since transcript discovery is performed on 2um resolution
+    // first, before aggregating the SE object to the other lower resolutions
+    ch_tissue_positions_002um = CONVERT_TISSUE_POSITIONS.out.csv.filter { resolution, _csv -> resolution == '002um' }.map { _resolution, csv -> csv }
+    ch_tissue_positions_bins  = CONVERT_TISSUE_POSITIONS.out.csv.filter { resolution, _csv -> resolution != '002um' }
+
+    emit:
+    bam                    = FILTER_BARCODED_BAM.out.bam
+    sample_name            = ch_sample.map { sample, _path, _meta -> sample }.first()
+    tissue_positions_002um = ch_tissue_positions_002um
+    tissue_positions_bins  = ch_tissue_positions_bins
+    barcode_mappings       = CONVERT_BARCODE_MAPPINGS.out.csv
+}