Save segger segmentation as a SOPA-compatible SpatialData object

[Tobiaspk]

Add support for a SOPA-compatible Segger output. Multiple PRs try to tackle that, but no agreed-upon integration and output format yet.

Proposal based on a discussion with @EliHei2

Outputs: SpatialData object, containing

• Shapes: Cell Masks
• Points: Transcript
• Table: Anndata (reuse existing scripts)

Integration: Implement a new segger export interface with parameters

• Input: Path to seggers output folder
• Output: Defaults to input path
• Overwrite: Default false

Open questions

• Cell mask method: delauney, delauney + smoothing or convex_hull? (after chatting with Eli, I err on side of delauney + minimal smoothing
• Transcript: Should SpatialData include all, or only segmented transcripts?
• Are images required? Or can we skip in first iteration?
• Which other parameters are needed? Let's keep it minimal.
• Do we need to transform coordinate systems?
• Which (segger) metadata should be included? (segger version)
• Are you able to install SpatialData into the Segger environment? Or should this live in a lightweight segger-io package?

Eli, could you also share the segmentation masks that you generated? Thanks!

[EliHei2]

Hi @Tobiaspk, thanks for making this! re inclluding all transcripts --> I think so, because the object will be full and and one could apply other seg. methods to this and do other downstream anslysis or even other seg. masks form segger itself.
, these are the masks for `delauney + smoothing:

[Tobiaspk]

Got it. In what scenario does a user run another seg method on top of segger?

As a user, I would expect that segger outputs a clean, filtered dataset with assigned transcripts only per default. Similar to other pipelines. I see segger_segmentation.parquet as the raw-data that mostly advanced users work with. I think we have to assume that a typical user does not know about filtering segger thresholds and control genes.

Thoughts? Maybe add a --use-raw param to include all transcripts?

[EliHei2]

Hey Tobz, sorry I rushed through before leaving work today. These are some of more concrete reasons I think it'd be good to have full and not only assigned transcripts:

• The same way 10x outputs their transcript results, it would be nice for the user to check out which/why/how many transcripts are left out. --> moreover, a different thresholding on segger scores (which turns out to be pretty much confidence) could be applied by the user to tighten or loosen the assignment.
• In paper (new verison) we show how different receptive fields give different levels of specificity/sensitivity. A user might want to check for different levels of expansion to look at the unassigned transcripts of interest.
• segger's envisioned fragment mode (see branch) as well as tools like troutpy argue that there's a lot of mis-segmented transcripts and one would might want to apply a regrouping on the unassigned transcripts to get "cell-like" groups or "fragments".
• as mentioned before, in benchmarks like singlecell open benchmarks continuously would run on multiple segmentations to compare them, that is, having all transcripts would be essential to compare.

I, however, understand the gegen-argument :)), that is the user would just care about one segger run and let's just drop all others. But even in segger_segmentation.parquet as far as I know all tx are kept even the unassigned ones.

Would love to have some output from @andrewmoorman and @LucaMarconato and @quentinblampey.

[EliHei2]

and I missed it: I like the idea of --use-raw or --include-all-transctips or --only-assigned. maybe this is actually the idea way to go. with some sensible default (I agree only assigned should be the default!)

[quentinblampey]

Hi,

Really great to see you working on such an integration!

For other segmentation methods within Sopa, we keep the default SpatialData object obtained via spatialdata-io (+ some sopa-related spatial elements), and we just add on top of it (i) one shape element and (ii) one table.

NB: we could also decide to support segger directly from sopa (i.e., converting segger's output on sopa's side), but both work and I'll of course let you decide what's best for you!

[Tobiaspk]

Thanks @EliHei2 for the context and use cases. Yeah I'm good to move ahead with a --include-all-transcripts flag (i like verbose names).

Quick follow-up: You mentioned you'd like to take write_anndata out of the main segger function, and provide this as a segger export anndata interface instead. Want to do the same for polygons, as in segger export boundaries for example?

@quentinblampey for a direct sopa integration do you only need an AnnData h5ad and a geo-dataframe with polygons? How about coordinate transformations? I'd love to prioritise this, ideally with direct support from sopa. Lmk how we can help with this.

[Tobiaspk]

@EliHei2 Gentle bump. Which scenario would you prefer in this iteration? Assume that both still provide a segger export anndata/boundaries [--include-all-transcripts] [--<other-params>] interface:

1st User-focused - Output all these per default (names up for debate):

• transcripts_all.parquet
• transcripts_assigned.parquet
• adata.h5ad
• cell_boundaries.parquet

2nd Flexibility - Output only transcripts_all.parquet and require user-built pipelines.

@quentinblampey Also quick bump - do we need coordinate transformations or other info for the integration?

[EliHei2]

Hola @Tobiaspk I didn't quite get your Q. but if you mean intermediate resutls, I'd output only what they asked, indeed a spatialdata object for example would include all of these. but anndata, for instance, would not need boundaries or transcripts for that matter (cell area is computed using convex hull there). Did I answer your question or did I not get it?

[quentinblampey]

for a direct sopa integration do you only need an AnnData h5ad and a geo-dataframe with polygons? How about coordinate transformations? I'd love to prioritise this, ideally with direct support from sopa. Lmk how we can help with this.

Yes, an AnnData and a GeoDataFrame are enough, except if you don't use the same coordinate system as the transcripts. If you use the same coordinate system, we can just use the transformations of the transcripts for the GeoDataFrame (as we do within Sopa for baysor/proseg/comseg), else I'll need the coordinate transformations

[EliHei2]

Thanks @quentinblampey. @Tobiaspk if you're happy with the the overall plan, shall we ask @quentinblampey kindly to go through the PR or would you liek me to break it down to smaller PRs? I think for a true sopa export all of the features would be needed. what do you think?

[Tobiaspk]

@quentinblampey Sounds good, thanks! Shapes and transcripts will share the coordinate space.

@EliHei2 Yeah, my understanding is that segger only needs to output the anndata, and the geodataframe with cell boundary polygons to be SOPA compatible. Let's start with that, and create a SpatialData interface in the next PR. Does that work for you?

Needed:

• segger export anndata/boundaries interface (I'd be open to generating these per default)
• implementation for delauney + minimal smoothing
• short paragraph in README.md describing the outputs
• create an issue/PR in SOPA to support segger outputs

[EliHei2]

Hey Tobias, I'm generating a minimal pr, but I think although SOPA is one of the main outputs, many people would like to export other formats, e.g., xenium bundle/explorer output. so I suggest this (which is pretty much already implemented but I'm breaking down):

segger export anndata|sopa|xenium

what do you think? @Tobiaspk ?

[Tobiaspk]

Agreed. :)

Let's start with the new interface for anndata & boundaries and get this right. From this, it will be easy to re-bundle xenium and spatialdata, and in parallel start the SOPA integration.