Repository navigation
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #427 +/- ##
==========================================
+ Coverage 65.79% 66.60% +0.80%
==========================================
Files 26 26
Lines 3263 3345 +82
==========================================
+ Hits 2147 2228 +81
- Misses 1116 1117 +1
🚀 New features to boost your workflow:
|
b741f27 to
8150719
Compare
XOA < 1.3.0 exports, and GEO deposits reduced to the CSV outputs, ship no parquet and no cells.zarr.zip. xenium() opened cells.zarr.zip (via _XeniumCells.open) and read the parquet files unconditionally, so it raised FileNotFoundError on these bundles before any CSV-based path could run. When cells.zarr.zip and the parquet files are absent: - read the table metadata, cell/nucleus boundaries, and transcripts from the .csv[.gz] outputs (columns identical to the parquet ones); - reconstruct the raster cell/nucleus labels by filling each boundary polygon with its integer label into a uint32 raster (as the zarr masks are) -- keyed by cell_id for pre-1.3.0 bundles, or by the label_index recovered from the boundary CSVs' label_id column for v2/v3 (hex cell_id, multinucleate nuclei supported). The table's cell_labels column and the raster share these exact values, honoring labels_models_kwargs (chunking/multiscale); - read the cell feature matrix from a cell_feature_matrix/ MatrixMarket directory or cell_feature_matrix.tar.gz when cell_feature_matrix.h5 is missing. The zarr path is gated on cells.zarr.zip existing, not on the XOA version. Mapping the table to cell_labels degrades gracefully (warn, keep circles) when a cell lacks a boundary, rather than raising. Add committed CSV-only fixtures -- 30 cells of the 10x Mouse Brain v1.0.2 dataset (integer ids) and 12 cells of the Xenium Prime Mouse Brain v3 dataset (hex ids + label_id) -- with tests asserting the reconstructed labels match the table's cell_labels, plus the circles, uncompressed-.csv, and mtx variants.
8150719 to
642f014
Compare
|
|
||
| def _get_labels_from_boundaries( | ||
| shapes: GeoDataFrame, | ||
| specs: dict[str, Any], | ||
| indices_mapping: pd.DataFrame | None = None, | ||
| labels_models_kwargs: Mapping[str, Any] = MappingProxyType({}), | ||
| ) -> DataArray: | ||
| """Reconstruct a raster labels element from boundary polygons. | ||
|
|
||
| CSV-only bundles have no cells.zarr.zip mask arrays, so the labels are rebuilt by filling each | ||
| boundary polygon with its integer label into a ``uint32`` raster (as the zarr masks are). The | ||
| label is the integer ``cell_id`` (the GeoDataFrame index) for pre-1.3.0 bundles, or the | ||
| ``label_index`` that ``indices_mapping`` maps the string ``cell_id`` to for v2/v3 bundles. The | ||
| raster spans the polygons' bounding box with a transform that keeps it aligned in ``global``. | ||
| """ | ||
| if indices_mapping is not None and not pd.api.types.is_integer_dtype(shapes.index): | ||
| label = _cell_id_to_label_index(indices_mapping).loc[shapes.index].to_numpy() | ||
| else: | ||
| label = shapes.index.to_numpy() | ||
| # geometry is in microns; global == pixels == microns / pixel_size | ||
| inv = 1.0 / specs["pixel_size"] | ||
| minx, miny, maxx, maxy = shapes.total_bounds * inv | ||
| x0, y0 = math.floor(minx), math.floor(miny) | ||
| raster = np.zeros((math.ceil(maxy) - y0, math.ceil(maxx) - x0), dtype=np.uint32) | ||
| for geom, value in zip(shapes.geometry.to_numpy(), label.astype(np.uint32), strict=True): | ||
| xs, ys = geom.exterior.coords.xy | ||
| rr, cc = polygon(np.asarray(ys) * inv - y0, np.asarray(xs) * inv - x0, shape=raster.shape) | ||
| raster[rr, cc] = value | ||
| transform = Translation([x0, y0], axes=("x", "y")) | ||
| return Labels2DModel.parse(raster, dims=("y", "x"), transformations={"global": transform}, **labels_models_kwargs) | ||
|
|
There was a problem hiding this comment.
This could be replaced by a call of spatialdata.rasterize()
There was a problem hiding this comment.
But the code is short and clear because it makes use of the knowledge of the Xenium data looks like and we are rasterizing something specific. So not a huge deal to leave as this.
There was a problem hiding this comment.
To be fair, spatialdata.rasterize() would have to address this: scverse/spatialdata#987 (easy). An advantage would be that rasterize() should be faster as it removes the for loop. But if you are not concerned about performance, we can put a comment here about this possible improvement, and leave as it.
There was a problem hiding this comment.
Yeah, then let's open this as a follow-up maybe issue. It's for legacy data so we probably don't have to overoptimize this, newer data won't use this path. Added a comment so that we can find it later again
| if cells_zarr_ctx is not None: | ||
| metadata = cells_zarr_ctx.get_cell_metadata(path) | ||
| else: | ||
| metadata = _read_cell_metadata(path) |
There was a problem hiding this comment.
Do we even need this if else? Can't we just call metadata = _read_cell metadata(path), which is what ultimately cells_zarr_ctx.get_cell_metadata(path) does?
| """Build the cell_id <-> label_index mapping from a boundary CSV's ``label_id`` column. | ||
|
|
||
| v2/v3 boundary CSVs carry a ``label_id`` column (the same integer the zarr stores as | ||
| ``label_index``); pre-1.3.0 CSVs do not, in which case there is nothing to map and ``None`` is | ||
| returned so labels fall back to the integer ``cell_id``. |
There was a problem hiding this comment.
Good, that was a critical point that I'm happy to see to be covered.
|
Thanks Tim, I'll leave some comments directly in the code. One here:
|
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
| if cells_zarr_ctx is not None: | ||
| cell_im = cells_zarr_ctx.cell_indices_mapping | ||
| nucleus_im = cells_zarr_ctx.nucleus_indices_mapping | ||
| else: | ||
| cell_im = ( | ||
| _csv_indices_mapping(path, XeniumKeys.CELL_BOUNDARIES_FILE_CSV) | ||
| if (cells_boundaries or cells_labels) | ||
| else None | ||
| ) | ||
| nucleus_im = ( | ||
| _csv_indices_mapping(path, XeniumKeys.NUCLEUS_BOUNDARIES_FILE_CSV) | ||
| if (nucleus_boundaries or nucleus_labels) | ||
| else None | ||
| ) |
There was a problem hiding this comment.
This assumes no zarr available -> .csv case. In theory one could find a store without Zarr, with Parquet, without .csv. This would crash the reader.
But it would have also crashed the reader before this PR. So we can add a comment for the future, but leave as is.
| minx, miny, maxx, maxy = shapes.total_bounds * inv | ||
| x0, y0 = math.floor(minx), math.floor(miny) | ||
| raster = np.zeros((math.ceil(maxy) - y0, math.ceil(maxx) - x0), dtype=np.uint32) | ||
| for geom, value in zip(shapes.geometry.to_numpy(), label.astype(np.uint32), strict=True): |
There was a problem hiding this comment.
This assumes that labels is int. I think that for XOA versions >= 1.3.0 and < 2.0.0, we are in the else case (line 728) because indices_mapping is None, while shapes has non-integer dtype.
If this happens, we could get an exception when label.astype(np.uint32) is called.
I suggest to:
- if you can test on real data, please do, and add a small test dataset
- if you can't, please leave a comment on the code (on the "else" above), but we could leave the code as is as we are not introducing a new bug: the previous code, in this rare edge case, would have been failing due to the lack of the Zarr file. It could be that such bug will never occur in practice.
There was a problem hiding this comment.
now raising, added comment
| if h5.is_file(): | ||
| adata = sc.read_10x_h5(h5, gex_only=gex_only) | ||
| # Undo fixed-point scaling factor applied to Xenium Protein data stored in HDF5. | ||
| with h5py.File(h5, "r") as f: | ||
| if "protein_scaling_factor" in f.attrs: | ||
| protein_feats = np.flatnonzero(adata.var["feature_types"] == "Protein Expression") | ||
| if len(protein_feats) > 0: | ||
| adata.X[:, protein_feats] /= f.attrs["protein_scaling_factor"] | ||
| return adata |
There was a problem hiding this comment.
We don't know if the protein data needs descaling as we do for when we read it from HDF5. I'd suggest to verify this on real data.
Or if this is not possible, conservatively we could set gex_only=False for when we read the matrix outside hdf5.
There was a problem hiding this comment.
Yeah, we can't get that info if it's not hdf5, added a warning
- attribute the two CSV test fixtures (10x datasets, CC BY 4.0): source, license, and the edits made (tests/fixtures/ATTRIBUTION.md) - warn when protein counts are read from a non-HDF5 matrix: they stay in raw scaled units because protein_scaling_factor lives only in the .h5 (verified against the 10x Protein Kidney bundle: mtx protein == scaled h5, factor 10) - raise a clear error instead of a cryptic astype() crash when boundary polygons have non-integer cell_ids but no label_index mapping - comment the rasterize() alternative (scverse/spatialdata#987) and the parquet-without-csv assumption
3 participants
GEO Xenium deposits are commonly stripped to the flat outputs: no
cells.zarr.zip, no parquet, and the counts matrix as MatrixMarket instead ofcell_feature_matrix.h5.xenium()openedcells.zarr.zipand read the parquet unconditionally, so these raisedFileNotFoundErrorbefore any fallback could run.When those are absent, read the equivalents:
.csv[.gz](columns identical to the parquet)cell_id, or bylabel_indexfrom the boundary CSVs'label_idcolumn for hex-id v2/v3 bundles (multinucleate nuclei included)cell_feature_matrix/MatrixMarket dir orcell_feature_matrix.tar.gzwhen.h5is missingThe zarr path is gated on
cells.zarr.zipexisting, not on the XOA version — hex-id ≥1.3.0 runs are deposited without it.GEO examples this targets:
.h5: GSE271693Tests: two committed CSV-only fixtures (10x Mouse Brain v1.0.2, integer id; Xenium Prime v3, hex id +
label_id). Also checked against the full 10x v1.0.2 bundle (reconstructed label id-sets match the zarr masks, 36,602 cells) and a real GEO deposit (GSE283843, 13,178 cells).Note: GEO strips
experiment.xenium; a reconstructed manifest must carrypixel_size(the reader needs it for the micron↔pixel scale and does not guess it).