Skip to content

[main] Add spatialdata element to segger export - #80

Open
Tobiaspk wants to merge 19 commits into
mainfrom
feature/spatialdata_export
Open

[main] Add spatialdata element to segger export#80
Tobiaspk wants to merge 19 commits into
mainfrom
feature/spatialdata_export

Conversation

@Tobiaspk

@Tobiaspk Tobiaspk commented Aug 25, 2026

Copy link
Copy Markdown
Collaborator

Added better support for spatialdata in response to #25 #48 #67 #73.

Adds a spatialdata element to segger export. Requires that sdata.zarr already exists (e.g. created with spatialdata_io.<technology>). Run like this:

pixi run -e cuda121 segger export spatialdata \
  -s "$PATH_OUTPUT/segger_segmentation.parquet" \
  --sdata "$PATH_INPUT/sdata.zarr"

Output — this edits sdata.zarr in place:

  • Appends segger's per-transcript columns (segger_cell_id, segger_similarity, segger_similarity_threshold, segger_converged, segger_filtered, segger_seen) to the existing transcripts points element.
  • Adds a cell_boundaries_segger shapes element (cell polygons).
  • Adds a table_segger table element (cell × gene AnnData).

Element names are configurable:

  • --sdata-transcripts-name (default transcripts) — existing points element to append segger's columns to.
  • --sdata-cell-boundaries-name (default cell_boundaries_segger)
  • --sdata-table-name (default table_segger)

Notes:

  • Fails fast with FileExistsError if the target boundaries/table name already exists, or KeyError if the transcripts element is missing.
  • The transcripts element keeps all of its original transcripts (including unassigned and negative-control ones); segger's columns are null for transcripts it did not assign, and segger_seen marks the ones it did.
  • --min-counts (default 10) drops cells with fewer assigned transcripts from the table. Must be >= 3 for spatialdata, since boundaries need >= 3 transcripts to form a polygon. The table is aligned to the boundaries, so it never annotates a missing shape.

Implementation notes:

Other changes:

@Tobiaspk Tobiaspk changed the title Add spatialdata element to segger export [main] \Add spatialdata element to segger export Aug 25, 2026
@Tobiaspk Tobiaspk changed the title [main] \Add spatialdata element to segger export [main] Add spatialdata element to segger export Aug 25, 2026
@Tobiaspk

Copy link
Copy Markdown
Collaborator Author

Example user workflow:

# define inputs and outputs
PATH_INPUT=../../data/inputs/xenium_sample_a
PATH_OUTPUT=../../data/outputs/xenium_sample_a

# save xenium as a spatialdata object (this can take a bit)
conda run -n sc python -c "
import spatialdata_io
sdata = spatialdata_io.xenium('$PATH_INPUT')
sdata.write('$PATH_INPUT/sdata.zarr')
"

# segment with segger
pixi run -e cuda121 segger segment \
  --input-directory "$PATH_INPUT" \
  --output-directory "$PATH_OUTPUT"

# export spatialdata
pixi run -e cuda121 segger export spatialdata \
  -s "$PATH_OUTPUT/segger_segmentation.parquet" \
  --sdata "$PATH_INPUT/sdata.zarr"

@Tobiaspk

Tobiaspk commented Aug 27, 2026

Copy link
Copy Markdown
Collaborator Author

Asking for your feedback on this workflow @MeyerBender, @quentinblampey. This adds segger elements to an existing spatialdata object. Would this be useful for you?

@quentinblampey

Copy link
Copy Markdown

Tagging @alihamraoui since he started working on integrating segger to sopa

@MeyerBender

Copy link
Copy Markdown

Thanks a lot, looks really cool and will be super useful! One thing I don't fully understand is why the transcripts_segger points would only contain the assigned transcripts. Since we are working on QC for ST segmentation, it is often useful to have all points in the transcripts, even the unassigned ones. Wouldn't this be what the filtered column in the points df would be for anyway? Or alternatively, could there be a switch allowing users to in-/exclude those transcripts in transcripts_segger?

Also tagging @LazDaria.

@EliHei2

EliHei2 commented Aug 28, 2026

Copy link
Copy Markdown
Collaborator

Agree with @MeyerBender I belive keeping all transcripts is essential for downstream QC, so default should contain all, while user can select only to keep assigned transcripts.

@Tobiaspk

Copy link
Copy Markdown
Collaborator Author

Appreciate your inputs. How about negative controls and antisense probes? And low quality probes (qv < 20)?

If we keep all the same transcripts, it doesn't necessarily make sense to store a new element. Perhaps just add a column with seggers segmentation to the existing transcripts element?

@EliHei2 with the new filtered flag, we could default to all transcripts. Less chance of error for the user. Same question as above though.

@MeyerBender

Copy link
Copy Markdown

I think that would be a good idea, having just one points element should also make things less ambiguous. I would even keep quality probes in there, and just have them as "Unassigned". I can ofc only speak for SegTraQ, but we usually filter out quality probes beforehand, and then only look at the "regular" gene probes to compute the number of unassigned transcripts. What would be important is that the number of transcripts assigned in the points element matches the number in the table exactly (in the old ProSeg version, this wasn't enforced, which led to quite some inconsistencies, see here).

@Tobiaspk

Copy link
Copy Markdown
Collaborator Author

I see, this is useful. We could join the following info with configurable prefix to the existing transcripts then:

  • segger_similarity
  • segger_cell_id
  • segger_similarity_threshold
  • segger_converged
  • segger_filtered

How about the cell_boundaries and table? Should we override, or add new elements?

-- we don't use spatialdata, so really reffering to your inputs here. thanks @MeyerBender

@Tobiaspk

Copy link
Copy Markdown
Collaborator Author

Are there naming conventions or other aspects to look out for for SOPA integration? @alihamraoui

@MeyerBender

Copy link
Copy Markdown

Sounds like a resonable plan, from my POV the most important ones would be segger_cell_id and segger_filtered. I would set segger_cell_id to None for all cells that are unassigned, regardless of the reason (low QV, control probe, or no convergence). And segger_filtered would then simply be a boolean that indicates if segger even attempted to assign the transcript or filtered it out beforehand, correct?

1 similar comment
@MeyerBender

Copy link
Copy Markdown

Sounds like a resonable plan, from my POV the most important ones would be segger_cell_id and segger_filtered. I would set segger_cell_id to None for all cells that are unassigned, regardless of the reason (low QV, control probe, or no convergence). And segger_filtered would then simply be a boolean that indicates if segger even attempted to assign the transcript or filtered it out beforehand, correct?

tobiaspk and others added 14 commits August 31, 2026 15:58
Runs cell_boundary for each cell in a ProcessPoolExecutor sized to the
SLURM-allocated CPU count instead of a single-threaded loop.
spatialdata's xarray_schema dependency imports pkg_resources, which
comes from setuptools; not otherwise pulled in by the resolved env.
zip(*results) produced tuples, which geopandas' set_geometry mistakes
for a column-name lookup instead of geometry values, dumping every
polygon into the error message. Convert to lists. Also leave one CPU
free instead of using every allocated core.
- Disable dask's query-planning backend at the top of segger/__init__.py:
  cudf is imported unconditionally by every CLI subcommand and pulls in
  dask.dataframe before export code ever runs, so setting the config
  lazily inside the export path was too late.
- _write_to_sdata now adds elements directly to the given --sdata store
  instead of copytree-ing it into the output directory first; avoids
  the 'cannot overwrite' error from the raw Xenium sdata's transcripts/
  cell_boundaries/table already occupying those element names.
- Fail fast: check the target element names don't already exist before
  doing any of the (expensive) segmentation/boundary work.
- Add --spatialdata-element-prefix (default '_segger') to name segger's
  elements distinctly from the source sdata's own.
- -o/--output-directory is now optional, only required for the
  anndata/transcripts/boundaries elements.
Points and table centroids include z when present; no-op until a reader emits it.
Spatialdata writes now fail on a name collision instead of overwriting.

Co-authored-by: Elyas Heidari <55977725+EliHei2@users.noreply.github.com>
Co-authored-by: enric-bazz <enrico.bazzacco02@outlook.it>
…w one

Per PR #80 discussion: transcripts_segger duplicated the sdata transcripts
table and dropped unassigned/control transcripts. Now left-joins segger's
columns (segger_cell_id, segger_similarity, segger_similarity_threshold,
segger_converged, segger_filtered) onto the existing points element in
place, with unassigned transcripts left null and a new segger_unused flag
for rows segger never saw at all.
Rework transcripts loading around a single load_transcripts() that trusts
segger's own 'filtered' flag instead of recomputing similarity thresholds
client-side, and rename the sdata element options to explicit
--sdata-transcripts-name/--sdata-cell-boundaries-name/--sdata-table-name.
Fix two gaps that would have broken the spatialdata path: load_transcripts
now returns the full per-transcript frame alongside the assigned subset
(needed to append segger's columns onto the sdata transcripts table), and
_merge_sdata_transcripts keeps row_index through its select so the join
onto the sdata table doesn't fail.
@Tobiaspk
Tobiaspk force-pushed the feature/spatialdata_export branch from aee0b10 to a1e223b Compare August 31, 2026 20:00
Matches _write_to_sdata's signature style instead of a dict.
…azily in dask

PointsModel/ShapesModel.parse default to an Identity transform on 'global' when
none is passed, silently dropping whatever scale/affine the base transcripts
element carried. _write_to_sdata now reads it via get_transformation() and
passes it through to both the rewritten transcripts and the new boundaries.

_merge_sdata_transcripts no longer round-trips the base transcripts through
pandas/polars/pandas; it merges lazily in dask (broadcast join, cumsum-built
row_index) and is only materialized once, right before write_element()
overwrites the on-disk store in place (required, since write_element refuses
to overwrite a path a live dask graph still reads from).

include_all_transcripts now defaults to True.
write_element(overwrite=True) always fails for an element inside its own
store (discussion #520). Back up transcripts, detach in-memory handles,
then delete-then-write. Also drop the row_index join key from output.
Drop cells below min_counts from the table/anndata, and keep only cells
that produced a boundary polygon so the table annotates only real shapes
(no dangling SpatialData annotations). Require min_counts >= 3 for
spatialdata, since boundaries need >= 3 transcripts.
@Tobiaspk

Tobiaspk commented Sep 1, 2026

Copy link
Copy Markdown
Collaborator Author

Current usage:

pixi run -e cuda121 segger export spatialdata \
  -s "$PATH_OUTPUT/segger_segmentation.parquet" \
  --sdata "$PATH_INPUT/sdata.zarr" \

  # below: optional element names, showing defaults:
  --sdata-transcripts-name transcripts \
  --sdata-cell-boundaries-name cell_boundaries_segger \
  --sdata-table-name table_segger

This

  1. asserts that transcripts exists, but boundaries and table doesn't
  2. append segger columns like above with None for unassigned, and write the other two elements.

Other changes

  • Make load_transcripts more human readable, using seggers new output structure + backwards compatibility
  • Change include-all-transcripts to true by default
  • Filter at least 3 transcripts per cell (default --mincounts 10), otherwise cells could be in adata but not in cell boundaries.

@MeyerBender @EliHei2 do you want to test this? Shouldn't need to rerun segger, just the export should be fine.

Tobiaspk added a commit that referenced this pull request Sep 9, 2026
…> export spatialdata)

Anticipates #80: segger export spatialdata edits an existing SpatialData
Zarr store in place rather than copying into the output directory, so
update the flag table and outputs page to match, and add the recommended
Xenium workflow to the quickstart.
@MeyerBender

Copy link
Copy Markdown

@Tobiaspk it took me a while to get the installation to work, I think the issue is that the current installation instructions install two competing versions of cupy. I did manage to manually make it work in the end, but now I keep running into RuntimeError: cuSpatial returned an invalid quadtree after 6 attempts (see segger issue #40).. As per issue #81, I will just try a couple more times and see if that fixes it, but if you have any better suggestions please let me know.

@alihamraoui

Copy link
Copy Markdown

@MeyerBender,
It took me a while too, turned out to be numpy getting pulled to 2.x (breaks cupy/RAPIDS), plus a few undeclared deps (opencv, scikit-image, geopandas, shapely...). Switched to uv and pinned everything, works reliably now: https://github.com/alihamraoui/segger/blob/main/pyproject.toml

Python 3.11, uv venv, then uv pip install -e .

@MeyerBender

Copy link
Copy Markdown

Cool, will have a look, thanks!

I fixed the quadtree bug by changing max_size += 10000 → max_size = int(max_size * 1.5) in quadtree.py and adding max_retries=20 in get_quadtree_index() in query.py.

It seems like segger is running now, but I cannot install spatialdata into the same environment. Specifically, I get the following error: spatialdata 0.7.3 requires dask<2026.1.2,>=2025.12.0, but you have dask 2024.9.0 which is incompatible.

I will try to run the segmentation and export in two separate environments now to see if it works that way.

@alihamraoui

Copy link
Copy Markdown

@MeyerBender that's expected, RAPIDS pins dask to an old version, and spatialdata needs dask ≥2025.12. Worse, they can't be reconciled anytime soon, because RAPIDS stopped publishing cuSpatial (v25.06, repo's archived), so segger is stuck on old RAPIDS/dask as long as it depends on cuspatial @Tobiaspk.

I'm running segger (uv) and spatialdata (conda) in separate environments for now.

@quentinblampey

Copy link
Copy Markdown

Thanks @alihamraoui for these details! This will be an issue for containerization, no?

@alihamraoui

Copy link
Copy Markdown

I think so, two images/two venvs, with segger writing to disk and a second container picking it up for the spatialdata step. But I think it'll work fine for nf-core/sopa. Hope the maintainers find a better alternative.

@Tobiaspk

Tobiaspk commented Sep 11, 2026

Copy link
Copy Markdown
Collaborator Author

Thanks for looking into this @MeyerBender. Did you try installing it with pixi? I set it up so that it should work with both. I just reset my environment, then ran this which installs the environment (segger + spatialdata) and it worked fine:

pixi run -e cuda121 segger export spatialdata \
  -s "$PATH_OUTPUT/segger_segmentation.parquet" \
  --sdata "$PATH_INPUT/sdata.zarr"

Unfortunately, cuspatial's quadtree is terribly inconsistent; it's random if a tree is "valid" or not, hence the repetitions are required. But, a new gpu accelerated package is expected to be released early next week that will allow us to move away from this dependency. We'll implement that immediately once available. cuspatial also relies on older packages, which complicates the installation and containerization as @alihamraoui pointed out. That should also be better once we moved away from cuspatial

Re-running export against transcripts that already have segger_*
columns caused merge() to silently suffix them with _x/_y instead of
overwriting, corrupting the segger_seen/segger_filtered computation.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants