diff --git a/notebooks/Xenium_pre-process.ipynb b/notebooks/Xenium_pre-process.ipynb index 7eb5fe954..7f13bcbfd 100644 --- a/notebooks/Xenium_pre-process.ipynb +++ b/notebooks/Xenium_pre-process.ipynb @@ -20,6 +20,18 @@ "text": [ "env: ANYWIDGET_HMR=1\n" ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/jishar/Documents/celldega/dega/lib/python3.13/site-packages/dask/dataframe/__init__.py:31: FutureWarning: The legacy Dask DataFrame implementation is deprecated and will be removed in a future version. Set the configuration option `dataframe.query-planning` to `True` or None to enable the new Dask Dataframe implementation and silence this warning.\n", + " warnings.warn(\n", + "/Users/jishar/Documents/celldega/dega/lib/python3.13/site-packages/xarray_schema/__init__.py:1: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.\n", + " from pkg_resources import DistributionNotFound, get_distribution\n", + "/Users/jishar/Documents/celldega/dega/lib/python3.13/site-packages/anndata/__init__.py:70: FutureWarning: Importing read_text from `anndata` is deprecated. Import anndata.io.read_text instead.\n", + " return module_get_attr_redirect(attr_name, deprecated_mapping=_DEPRECATED)\n" + ] } ], "source": [ @@ -49,13 +61,13 @@ "output_type": "stream", "text": [ "Starting preprocessing for sample: Xenium_V1_human_Pancreas_FFPE_outs\n", - "Created directory: data/landscape_files/Xenium_V1_human_Pancreas_FFPE_outs_test\n", + "Created directory: data/landscape_files/Xenium_V1_human_Pancreas_FFPE_outs_test_img_intensity_fix\n", "\n", "========Unzip and extract Xenium-related files========\n", "All files have been successfully extracted or skipped.\n", "\n", "========Write xenium transform file from the Zarr folder========\n", - "Transformation matrix saved to 'data/landscape_files/Xenium_V1_human_Pancreas_FFPE_outs_test/micron_to_image_transform.csv'.\n", + "Transformation matrix saved to 'data/landscape_files/Xenium_V1_human_Pancreas_FFPE_outs_test_img_intensity_fix/micron_to_image_transform.csv'.\n", "\n", "========Check if all required files or directories exist========\n", "All required files or directories for technology 'Xenium' are present in 'data/xenium_data/Xenium_V1_human_Pancreas_FFPE_outs'.\n", @@ -78,7 +90,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "/Users/jishar/Documents/celldega/src/celldega/pre/__init__.py:218: PerformanceWarning: Concatenating sparse arrays with multiple fill values: '[True, False]'. Picking the first and converting the rest.\n", + "/Users/jishar/Documents/celldega/src/celldega/pre/__init__.py:225: PerformanceWarning: Concatenating sparse arrays with multiple fill values: '[True, False]'. Picking the first and converting the rest.\n", " df_sig = df_sig.dropna(axis=1, how=\"all\")\n" ] }, @@ -93,11 +105,6 @@ "Calculating mean expression\n", "Calculating variance\n", "All meta gene files are succesfully saved.\n", - "data/landscape_files/Xenium_V1_human_Pancreas_FFPE_outs_test/cbg\n", - "cbg.index before mapping:\n", - "Index(['aaaadnje-1', 'aaacalai-1', 'aaacjgil-1', 'aaacpcil-1', 'aaadhocp-1'], dtype='object', name=0)\n", - "cbg.index after mapping:\n", - "Index([0, 1, 2, 3, 4], dtype='int64', name=0)\n", "Processing gene 0: ABCC11\n", "All gene-specific parquet files are succesfully saved.\n", "\n", @@ -121,15 +128,22 @@ "\n", "========Generating image tiles========\n", "------ xenium\n", - "generating dapi image tiles ...\n" + "Using morphology image: data/xenium_data/Xenium_V1_human_Pancreas_FFPE_outs/morphology_focus/morphology_focus_0000.ome.tif\n", + "OME shape: (4, 13770, 34155)\n", + "OME axes: CYX\n", + "OME dtype: uint16\n", + "Reading channel 'dapi' (index 0)\n", + "generating dapi image tiles ...\n", + "dapi: lo=0.00, p95=1596.00, gamma=1\n", + "Reading channel 'bound' (index 1)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ - "/Users/jishar/Documents/celldega/dega/lib/python3.13/site-packages/skimage/_shared/utils.py:328: UserWarning: /Users/jishar/Documents/celldega/notebooks/data/landscape_files/Xenium_V1_human_Pancreas_FFPE_outs_test/dapi_output_regular.tif is a low contrast image\n", - " return func(*args, **kwargs)\n" + "/Users/jishar/Documents/celldega/dega/lib/python3.13/site-packages/tifffile/tifffile.py:9310: UserWarning: reading array from closed file\n", + " warnings.warn(\n" ] }, { @@ -137,8 +151,41 @@ "output_type": "stream", "text": [ "generating bound image tiles ...\n", + "bound: lo=0.00, p95=2544.00, gamma=1\n", + "Reading channel 'rna' (index 2)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/jishar/Documents/celldega/dega/lib/python3.13/site-packages/tifffile/tifffile.py:9310: UserWarning: reading array from closed file\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ "generating rna image tiles ...\n", + "rna: lo=0.00, p95=4301.00, gamma=1\n", + "Reading channel 'prot' (index 3)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/jishar/Documents/celldega/dega/lib/python3.13/site-packages/tifffile/tifffile.py:9310: UserWarning: reading array from closed file\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ "generating prot image tiles ...\n", + "prot: lo=0.00, p95=1074.00, gamma=1\n", "Image tiles created successfully.\n", "\n", "======== Transcript Tiles========\n" @@ -148,8 +195,8 @@ "name": "stderr", "output_type": "stream", "text": [ - "Processing chunks: 100%|███████████████████████| 81/81 [00:00<00:00, 900.17it/s]\n", - "Processing coarse tiles: 84tile [00:21, 3.97tile/s]\n" + "Processing chunks: 100%|███████████████████████| 81/81 [00:00<00:00, 611.97it/s]\n", + "Processing coarse tiles: 84tile [00:20, 4.11tile/s]\n" ] }, { @@ -161,45 +208,14 @@ "======== Cell Boundary Tiles ========\n", "\n", "========Create cell boundary spatial tiles========\n", - "technology Xenium\n", - " geometry_micron \\\n", - "cell_id \n", - "0 POLYGON ((445.613 1697.663, 444.763 1698.300, ... \n", - "1 POLYGON ((442.850 1730.812, 440.938 1731.875, ... \n", - "2 POLYGON ((470.688 1706.163, 470.475 1706.375, ... \n", - "3 POLYGON ((429.888 1703.400, 429.038 1704.038, ... \n", - "4 POLYGON ((478.125 1702.125, 476.213 1703.188, ... \n", - "\n", - " GEOMETRY \\\n", - "cell_id \n", - "0 [[[2096.9999288922727, 7988.9998603837885], [2... \n", - "1 [[[2083.9998724224242, 8144.999389125], [2074.... \n", - "2 [[[2214.999833875, 8028.999857383789], [2213.9... \n", - "3 [[[2022.9999057198486, 8015.999513689697], [20... \n", - "4 [[[2249.99983125, 8009.999399249999], [2240.99... \n", - "\n", - " geometry center_x \\\n", - "cell_id \n", - "0 POLYGON ((2097.000 7989.000, 2093.000 7992.000... 2099.898777 \n", - "1 POLYGON ((2084.000 8144.999, 2075.000 8149.999... 2076.246077 \n", - "2 POLYGON ((2215.000 8029.000, 2214.000 8029.999... 2192.645181 \n", - "3 POLYGON ((2023.000 8016.000, 2019.000 8019.000... 2027.103432 \n", - "4 POLYGON ((2250.000 8009.999, 2241.000 8014.999... 2240.146290 \n", - "\n", - " center_y \n", - "cell_id \n", - "0 8005.963877 \n", - "1 8168.366465 \n", - "2 8057.422615 \n", - "3 8034.700824 \n", - "4 8051.613206 \n" + "technology Xenium\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ - "Processing coarse tiles: 100%|██████████████████| 14/14 [00:19<00:00, 1.41s/it]\n" + "Processing coarse tiles: 100%|██████████████████| 14/14 [00:17<00:00, 1.22s/it]\n" ] }, { @@ -217,7 +233,7 @@ "source": [ "sample = 'Xenium_V1_human_Pancreas_FFPE_outs'\n", "data_dir = f'data/xenium_data/'\n", - "path_landscape_files=f'data/landscape_files/{sample}_test'\n", + "path_landscape_files=f'data/landscape_files/{sample}_test_img_intensity_fix'\n", "\n", "tile_size=250\n", "\n", @@ -227,6 +243,9 @@ " tile_size=tile_size,\n", " path_landscape_files=path_landscape_files,\n", " use_int_index=True,\n", + " # use_row_groups=True,\n", + " upper_percentile=95,\n", + " white_level=40\n", " )" ] }, @@ -240,22 +259,22 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 4, "id": "4a1320eb", "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "80f3435bda344cbf890fcbe3d56da16b", + "model_id": "156e2217606e41f4be8f5de155e04781", "version_major": 2, "version_minor": 1 }, "text/plain": [ - "Landscape(base_url='http://localhost:53080/data/landscape_files/Xenium_V1_human_Pancreas_FFPE_outs_test', cell…" + "Landscape(base_url='http://localhost:64759/data/landscape_files/Xenium_V1_human_Pancreas_FFPE_outs_test_img_in…" ] }, - "execution_count": 3, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } diff --git a/src/celldega/pre/__init__.py b/src/celldega/pre/__init__.py index 5f76671f3..3d8cfe77d 100644 --- a/src/celldega/pre/__init__.py +++ b/src/celldega/pre/__init__.py @@ -334,20 +334,23 @@ def create_cluster_and_meta_cluster( return clusters -def _process_image_channel(path_dega_files, channel_info, img): +def _process_image_channel( + path_dega_files, + channel_info, + img, + upper_percentile=99, + scale_non_dapi=1.0, + white_level=100, +): """ - Process a single image channel for tiling. - - Parameters: - - path_dega_files: Landscape files path - - channel_info: Dictionary with channel information (name, index) - - img: Optional pre-loaded image array - - Returns: - - None + Process a single image channel for tiling using simple per-channel windowing, + similar to Xenium Explorer: + - choose a per-channel upper bound from a high percentile + - clip to [0, upper_bound] + - normalize to 0-1 + - convert to 8-bit for display """ channel_name = channel_info["name"] - channel_index = channel_info.get("index", 0) print(f"generating {channel_name} image tiles ...") @@ -355,22 +358,51 @@ def _process_image_channel(path_dega_files, channel_info, img): if pyramid_path.exists(): return - # Extract and process the channel - scale = 1 if channel_name.lower() == "dapi" else 2 # Adjust intensity for better visualization - if img.ndim == 3: - image_data = img[..., channel_index] * scale - elif img.ndim == 2: - image_data = img * scale + if img.ndim != 2: + raise ValueError(f"Expected a 2D channel image, got shape {img.shape}") + + channel = img.astype(np.float32) + + # Keep scaling neutral unless you intentionally want a boost + scale = 1.0 if channel_name.lower() == "dapi" else scale_non_dapi + channel *= scale + + # Per-channel display window (high-percentile upper bound) + sample = channel[::10, ::10] + + if not (0 <= upper_percentile <= 100): + raise ValueError( + f"upper_percentile must be between 0 and 100 (inclusive); got {upper_percentile!r}" + ) + + hi = np.percentile(sample, upper_percentile) + + print(f"{channel_name}: p{upper_percentile}={hi:.2f}") + + if hi > 0: + # Clip to display range and normalize from zero. + channel = np.clip(channel, 0.0, hi) + channel = channel / hi + + # Clamp white_level to the valid 8-bit display range [0, 255] + white_level_safe = float(white_level) + if white_level_safe < 0.0 or white_level_safe > 255.0: + warnings.warn( + f"white_level ({white_level_safe}) is outside [0, 255]; " + "clamping to this range for display.", + stacklevel=2, + ) + white_level_safe = min(255.0, max(0.0, white_level_safe)) + + image_data = (channel * white_level_safe).astype(np.uint8) else: - raise ValueError(f"Unsupported image dimensions: {img.ndim}. Expected 2D or 3D image.") + image_data = np.zeros_like(channel, dtype=np.uint8) output_path = Path(path_dega_files) / f"{channel_name}_output_regular.tif" - imsave(output_path, image_data) + imsave(output_path, image_data, check_contrast=False) - # Convert the image to PNG format image_png = _convert_to_png(str(output_path)) - # Create a DeepZoom pyramid for the channel make_deepzoom_pyramid( image_png, str(Path(path_dega_files) / "pyramid_images"), @@ -379,7 +411,14 @@ def _process_image_channel(path_dega_files, channel_info, img): ) -def create_image_tiles(technology, data_dir, path_dega_files, image_tile_layer="dapi"): +def create_image_tiles( + technology, + data_dir, + path_dega_files, + image_tile_layer="dapi", + upper_percentile=99, + white_level=100, +): """ Creates image tiles for visualization from the Xenium morphology image. @@ -389,6 +428,12 @@ def create_image_tiles(technology, data_dir, path_dega_files, image_tile_layer=" path_dega_files (str): Path to the directory where the image tiles and pyramid will be saved. image_tile_layer (str, optional): Specifies which image layers to process. Options for Xenium are 'dapi' (default) or 'all'. Use the filename of the .scn file for h&e Landscapes. + upper_percentile (float, optional): Upper intensity percentile used to rescale/clip the image + before tile generation. Must be between 0 and 100 (inclusive). Values close to 100 (e.g. 95-99) + reduce the influence of very bright outliers while preserving most detail. + white_level (float, optional): Factor controlling how intensities are mapped toward white during + normalization. Must be non-negative. Higher values produce brighter tiles; typical values are + in the range 40-255, with 100 as a balanced default. Raises: ValueError: If the specified technology is not supported or if the image_tile_layer is invalid. @@ -397,7 +442,13 @@ def create_image_tiles(technology, data_dir, path_dega_files, image_tile_layer=" print("\n========Generating image tiles========") if technology == "Xenium": print("------ xenium") - create_image_tiles_xenium(data_dir, path_dega_files, image_tile_layer=image_tile_layer) + create_image_tiles_xenium( + data_dir, + path_dega_files, + image_tile_layer=image_tile_layer, + upper_percentile=upper_percentile, + white_level=white_level, + ) elif technology == "MERSCOPE": print("------ merscope") create_image_tiles_merscope(data_dir, path_dega_files, image_tile_layer=image_tile_layer) @@ -457,40 +508,65 @@ def remove_intermediate_files(path_dega_files): file.unlink() -def create_image_tiles_xenium(data_dir, path_dega_files, image_tile_layer="dapi"): +def create_image_tiles_xenium( + data_dir, path_dega_files, image_tile_layer="dapi", upper_percentile=99, white_level=100 +): """ Creates image tiles for visualization from the Xenium morphology image. - Args: - data_dir (str): Path to the directory containing the data (e.g., morphology_focus_0000.ome.tif). - path_dega_files (str): Path to the directory where the image tiles and pyramid will be saved. - image_tile_layer (str, optional): Specifies which image layers to process. Options are 'dapi' (default) or 'all'. - Raises: - FileNotFoundError: If the required input image file is not found. + This version: + - resolves the morphology OME-TIFF via resolve_xenium_morphology_ome_path + - avoids loading the full OME-TIFF into memory + - reads one channel at a time from CYX images + - applies per-channel intensity windowing (upper_percentile / white_level) """ + if image_tile_layer not in ["dapi", "all"]: raise ValueError(f"Invalid image_tile_layer: {image_tile_layer}. Must be 'dapi' or 'all'.") file_path = resolve_xenium_morphology_ome_path(data_dir) + print(f"Using morphology image: {file_path}") - # Load the morphology image once if processing multiple channels - img = imread(file_path) - - if image_tile_layer == "all" and file_path.name == "morphology.ome.tif": - raise ValueError( - "image_tile_layer='all' needs a multi-channel morphology_focus OME-TIFF; " - "this bundle only has morphology.ome.tif. Use image_tile_layer='dapi' or " - "supply morphology_focus/*.ome.tif from the instrument output." + channel_map = [{"name": "dapi", "index": 0}] + if image_tile_layer == "all": + channel_map.extend( + [ + {"name": "bound", "index": 1}, + {"name": "rna", "index": 2}, + {"name": "prot", "index": 3}, + ] ) - # Process the DAPI channel - if image_tile_layer in ["dapi", "all"]: - _process_image_channel(path_dega_files, {"name": "dapi", "index": 0}, img) + # Use tifffile to safely read OME-TIFF without loading the full image + with tifffile.TiffFile(file_path) as tif: + series = tif.series[0] - # Process additional channels if image_tile_layer is 'all' - if image_tile_layer == "all": - for idx, channel in enumerate(["bound", "rna", "prot"]): - _process_image_channel(path_dega_files, {"name": channel, "index": idx + 1}, img) + print(f"OME shape: {series.shape}") + print(f"OME axes: {series.axes}") + print(f"OME dtype: {series.dtype}") + + # Ensure expected Xenium morphology layout + if series.axes != "CYX": + raise ValueError( + f"Expected Xenium morphology image axes to be 'CYX', got '{series.axes}'" + ) + + for channel_info in channel_map: + channel_name = channel_info["name"] + channel_index = channel_info["index"] + + print(f"Reading channel '{channel_name}' (index {channel_index})") + + # Read one channel at a time to avoid large memory usage + channel_2d = series.asarray(key=channel_index) + + _process_image_channel( + path_dega_files=path_dega_files, + channel_info={"name": channel_name, "index": 0}, + img=channel_2d, + upper_percentile=upper_percentile, + white_level=white_level, + ) remove_intermediate_files(path_dega_files) @@ -1678,6 +1754,7 @@ def _check_required_files(technology, data_dir): f"The following required files or directories are missing in directory '{data_dir}' " f"for technology '{technology}': {', '.join(missing_files_or_dir)}" ) + print( f"All required files or directories for technology '{technology}' are present in '{data_dir}'." ) diff --git a/src/celldega/pre/run_pre_processing.py b/src/celldega/pre/run_pre_processing.py index f44ddd24d..4cafc9550 100644 --- a/src/celldega/pre/run_pre_processing.py +++ b/src/celldega/pre/run_pre_processing.py @@ -155,6 +155,8 @@ def main( max_workers=1, use_row_groups=False, max_row_groups_per_file=400, + upper_percentile=99, + white_level=100, ): """ Main function to preprocess Xenium or MERSCOPE data and generate landscape files. @@ -168,11 +170,17 @@ def main( image_tile_layer (str): Image layers to be tiled. 'dapi' or 'all'. path_dega_files (str): Directory to save the landscape files. use_int_index (bool): Use integer index for smaller files and faster rendering. + max_workers (int): Maximum number of worker processes to use when generating + tiles. Defaults to 1 (no parallelism). use_row_groups (bool): If True, save tiles as row groups in chunked parquet files instead of individual tile files. Defaults to False. max_row_groups_per_file (int): Maximum row groups per parquet file when using row groups mode. Lower values create more files but avoid parquet-wasm memory issues with dense datasets. Defaults to 400. + upper_percentile (float): Upper percentile used for image intensity clipping when + generating image tiles. Defaults to 99. + white_level (float): White level (intensity scale) used when normalizing images + for tiling and visualization. Defaults to 100. Example: change directory to celldega, and run: @@ -311,7 +319,12 @@ def make_column_names_unique_fast(df): if technology in ["MERSCOPE", "Xenium"]: print("\n======== Image Tiles========") dega.pre.create_image_tiles( - technology, str(data_dir), path_dega_files, image_tile_layer=image_tile_layer + technology, + str(data_dir), + path_dega_files, + image_tile_layer=image_tile_layer, + upper_percentile=upper_percentile, + white_level=white_level, ) # Optionally pack image tiles into parquet row groups diff --git a/tests/unit/test_pre/test_row_groups.py b/tests/unit/test_pre/test_row_groups.py index fe6a9420d..bc5ac56e9 100644 --- a/tests/unit/test_pre/test_row_groups.py +++ b/tests/unit/test_pre/test_row_groups.py @@ -6,10 +6,10 @@ """ import json +from pathlib import Path import pyarrow as pa import pyarrow.parquet as pq -import pytest class TestRowGroupMetadata: @@ -141,7 +141,7 @@ def test_cbg_gene_metadata(self, tmp_path): with pq.ParquetWriter( str(output_path), schema_with_metadata, write_statistics=False ) as writer: - for gene in sorted(genes): + for _gene in sorted(genes): table = pa.Table.from_pydict( { "cell_id": [1, 2, 3], @@ -255,11 +255,11 @@ def test_landscape_parameters_structure(self, tmp_path): } output_path = tmp_path / "landscape_parameters.json" - with open(output_path, "w") as f: + with Path.open(output_path, "w") as f: json.dump(params, f, indent=2) # Read back and verify - with open(output_path) as f: + with Path.open(output_path) as f: read_params = json.load(f) # Verify key fields