How to register multiple .ome.tif at the same time & issue registering .ome.tif

[Boehmin]

Hi @Yu-AnChen !

I just tried palom on two cycles of FISH and the registration so far seems really really impressive (in particular since one of the cycles is half out of focus). First of all, thank you so much for making this tool available!

I was wondering how it would be possible to register multiple .ome.tif images simultaneously. Additionally, the registration worked fine with one of my datasets (converted from .h5 to .ome.tif which was a big pain already), however I have issues with a smaller .ome.tif dataset. I get the below error message:

(palom) [boehm0002@sgi27 ]$python3 palom_dan_ome_R1toR2.py
2023-09-04 14:12:42.717 | WARNING  | palom.reader:auto_format_pyramid:64 - Unable to detect pyramid levels, it may take a while to compute thumbnails during coarse alignment
2023-09-04 14:12:42.729 | WARNING  | palom.reader:auto_format_pyramid:64 - Unable to detect pyramid levels, it may take a while to compute thumbnails during coarse alignment
Traceback (most recent call last):
  File "/scicore/home/rueegg/boehm0002/spatial_transcriptomics/dans_test1/HiPlex_Example_DJH/palom_dan_ome_R1toR2.py", line 27, in <module>
    ref_thumbnail=c1r.read_level_channels(THUMBNAIL_LEVEL, 2).compute(),
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/dask/threaded.py", line 89, in get
    results = get_async(
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/dask/local.py", line 511, in get_async
    raise_exception(exc, tb)
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/dask/local.py", line 319, in reraise
    raise exc
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/dask/local.py", line 224, in execute_task
    result = _execute_task(task, data)
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/zarr/core.py", line 821, in __getitem__
    result = self.get_basic_selection(pure_selection, fields=fields)
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/zarr/core.py", line 947, in get_basic_selection
    return self._get_basic_selection_nd(selection=selection, out=out,
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/zarr/core.py", line 990, in _get_basic_selection_nd
    return self._get_selection(indexer=indexer, out=out, fields=fields)
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/zarr/core.py", line 1285, in _get_selection
    self._chunk_getitem(chunk_coords, chunk_selection, out, out_selection,
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/zarr/core.py", line 1994, in _chunk_getitem
    cdata = self.chunk_store[ckey]
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/zarr/storage.py", line 738, in __getitem__
    return self._mutable_mapping[key]
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/tifffile/tifffile.py", line 12274, in __getitem__
    return self._getitem(key)
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/tifffile/tifffile.py", line 12951, in _getitem
    keyframe, page, chunkindex, offset, bytecount = self._parse_key(key)
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/tifffile/tifffile.py", line 13050, in _parse_key
    page = series[pageindex]
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/tifffile/tifffile.py", line 12135, in __getitem__
    return self._getitem(int(key))
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/tifffile/tifffile.py", line 12115, in _getitem
    return self.parent.pages._getitem(page.index + key)
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/tifffile/tifffile.py", line 7568, in _getitem
    self._seek(key)
  File "/scicore/home/rueegg/boehm0002/miniconda3/envs/palom/lib/python3.10/site-packages/tifffile/tifffile.py", line 7402, in _seek
    raise ValueError('seek of closed file')
ValueError: seek of closed file

Thanks ever so much for your help!

[Yu-AnChen]

Hi @Boehmin thanks for the kind words :)

Palom might be able to read h5 files depending on the h5 file structure. If you have a small dataset to share I'll find time to take a look.

For the .ome.tif dataset that gave you the error, it looks to me that the file might have some issues. Would you be able to share that?

To register multiple files, our recommended approach is to register everything (cycles 2, 3, 4, and so on) to the first cycle. Please checkout the example script and lmk for any questions.

[Boehmin]

Hi @Yu-AnChen

Thanks for pointing me to the example script. I'll try this!

I have a smaller fused example h5 file (only one timepoint) and the ome tif dataset here. The .h5 consists of z-stacks so I do not know if that would be an issue for palom? If not that would be even better. Channel 4 is DAPI and channel 3 is a fiducial channel so in theory I could register over the fiducial channel as well if that helps?

Again, thanks for looking into this.

[Yu-AnChen]

Thanks for sharing the images, I was able to reproduce the error when using your ome-tiff files as input. Will look into a solution there.

On the h5 dataset (seems to be generated with BigDataViewer), I was able to convert it into a pyramidal ome-tiff using palom. Note that installing h5py in the environment that you run palom is required.

Since there are Z-stacks, I did a maximum intensity projection on the Z-axis in the following script (the output file size is ~50 MB; it can be opened in QuPath)

import palom
import h5py
import dask.array as da
import pathlib

img_path = pathlib.Path('h5_small/small_fuse_fused.h5')
h5_img = h5py.File(img_path)

# X-Y pixel size at full resolution; unit is µm/pixel
PIXEL_SIZE = 1
CHANNEL_NAMES = ['marker-1', 'marker-2', 'marker-fiducial', 'DNA']

mosaics = [
    # max projection on the Z axis
    da.from_array(h5_img[f"t00000/{channel}/0/cells"]).max(axis=0)
    for channel in h5_img['t00000']
]

palom.pyramid.write_pyramid(
    mosaics=mosaics,
    output_path=img_path.parent / f"{img_path.stem}-pyramid.ome.tif",
    pixel_size=PIXEL_SIZE,
    channel_names=CHANNEL_NAMES,
    downscale_factor=2,
    compression='zlib',
    tile_size=1024,
    save_RAM=True
)

With the output ome-tiffs from the above script, I think you should be able to run palom without running into the ValueError: seek of closed file error in your first comment.

[Yu-AnChen]

As for the issue of ValueError: seek of closed file, I've found a workaround. Please try it out on your end - I don't have a complete solution to it yet but hopefully that can get you up and running.

[lguerard]

Since there are Z-stacks, I did a maximum intensity projection on the Z-axis in the following script (the output file size is ~50 MB; it can be opened in QuPath)

Sorry to hijack this, is palom working on 2D files only or did you do a projection to speed it up ?
Thanks !

[Boehmin]

@Yu-AnChen, I tried both solutions and they both worked! Thank you.
Also curious whether this could work in 3D. :)

[Yu-AnChen]

The current image file reader and writer only works with images that have 2 or 3 dimensions. In our own application, it's (Channel/Band, Y, X). Palom could compute and find out the transformation/mapping between two (C, Z, Y, X) datasets and apply the transformation plane by plane, but the image file reader and writing will need some custom formatting.

[moataz-youssef]

Hi @Boehmin thanks for the kind words :)

Palom might be able to read h5 files depending on the h5 file structure. If you have a small dataset to share I'll find time to take a look.

For the .ome.tif dataset that gave you the error, it looks to me that the file might have some issues. Would you be able to share that?

To register multiple files, our recommended approach is to register everything (cycles 2, 3, 4, and so on) to the first cycle. Please checkout the example script and lmk for any questions.

Hello Yu-AnChen,

Sorry to hijack the thread and thanks for this great tool. I discovered your tool after hitting a wall with Ashlar, as my data is already stitched. I used your sample script in this thread and just passed in the location of 3 of my images (CyCIF runs, each DAPI + 3 other channels) as well as an output destination. My files were originally .VSI images, which I converted to OME-TIFF files using QuPath.

I get the following error:

2024-12-15 23:21:24.168 | DEBUG    | palom.register:cv2_feature_detect_and_match:288 - keypts L:5000, keypts R:5000
2024-12-15 23:21:24.261 | DEBUG    | palom.register:cv2_feature_detect_and_match:288 - keypts L:5000, keypts R:5000
2024-12-15 23:21:24.355 | DEBUG    | palom.register:cv2_feature_detect_and_match:288 - keypts L:4165, keypts R:4165
2024-12-15 23:21:31.402 | INFO     | palom.align:compute_shifts:203 - Computing block-wise shifts
Computing shifts: 100%|##############################################################################################################################################################################| 34021/34021 [01:49<00:00, 309.32it/s]
Traceback (most recent call last):
  File "D:\valis_container\palom_v1.py", line 74, in <module>
    fig = aligner.plot_shifts()
  File "C:\Users\yousseah\AppData\Local\anaconda3\envs\palom\lib\site-packages\palom\align.py", line 276, in plot_shifts
    viz_shifts(shifts, self.grid_shape, ax=axs[1])
  File "C:\Users\yousseah\AppData\Local\anaconda3\envs\palom\lib\site-packages\palom\align.py", line 90, in viz_shifts
    dcenter = skimage.filters.threshold_triangle(distances[is_finite])
  File "C:\Users\yousseah\AppData\Local\anaconda3\envs\palom\lib\site-packages\skimage\filters\thresholding.py", line 960, in threshold_triangle
    arg_level = np.argmax(length) + arg_low_level
  File "C:\Users\yousseah\AppData\Local\anaconda3\envs\palom\lib\site-packages\numpy\core\fromnumeric.py", line 1229, in argmax
    return _wrapfunc(a, 'argmax', axis=axis, out=out, **kwds)
  File "C:\Users\yousseah\AppData\Local\anaconda3\envs\palom\lib\site-packages\numpy\core\fromnumeric.py", line 59, in _wrapfunc
    return bound(*args, **kwds)
ValueError: attempt to get argmax of an empty sequence

[Yu-AnChen]

Hi @moataz-youssef , no worries and thanks for giving it a try. Apologies that I should have included this detail in my recent updates - VSI files are now supported through slideio (github repo) in palom v2024.10.2. You can pass the path to the VSI file to the palom.reader.VsiReader (ref).

From the above message, I'm guessing when converting to ome-tiff via QuPath, the tile-size wasn't specified and the default is 256x256 pixels. I found that small tile size will break how it's currently filtering out tiles that failed the alignment.

If you could start from the VSI files, the tile size will be set to 1024x1024 pixels internally, and the error you are running into would likely be gone.

It would be helpful if you could share a thumbnail of your scan. And here's a jupyter notebook that might give you some ideas.

[moataz-youssef]

Thanks a lot @Yu-AnChen.

I actually specified in QuPath that the tile size is 1024 pixels, but I will double check, what the final tile size in the converted OME-TIF files was. I can also share the scans/converted files with you privately, if you have time for this.

Thanks a lot for this great tool, I will update you on my progress.

[moataz-youssef]

@Yu-AnChen your solution worked! Using palom.reader.VsiReader worked and the script finished. The output had only 3 channels from the 3 original files (each with DAPI + 3 channels), so I should expect 12 channels. I can't seem to find how to do this, other than to specify that the registration should be done on channel 0 (DAPI). Do you have any suggestions for me? Thanks a lot in advance.

[moataz-youssef]

And just to add, the jupyter notebook also worked great and gave me a composite OME-TIFF with the DAPI channels from the 3 VSI images. How can I add the other channels through this method as well?

I am really impressed with the results, they are comparable, if not better, than manual adjustment and elastic deformation using Warpy in QuPath and Fiji.

[Yu-AnChen]

Thanks, and glad to hear that the alignment seemed to work well!

Would you be able to share the image or images from some test samples? I suspect that the num_channels property of the scene in the VSI file wasn't detected properly. If that's the case, I think it's also valuable to the slideio author.

[Yu-AnChen]

Re: the jupyter notebook that only gives you the first channel of all the "cycles", try changing the reader.pyramid[0][0] in block In [11]: to reader.pyramid[0], this should give you all the channels detected in the output.

[moataz-youssef]

Re: the jupyter notebook that only gives you the first channel of all the "cycles", try changing the reader.pyramid[0][0] in block In [11]: to reader.pyramid[0], this should give you all the channels detected in the output.

Thanks for this suggestion. Removing one [0] from reader.pyramid, gave me 7 channels. Removing one [0] from r1.pyramid in the same block. The output of r1.pyramid is:

 dask.array<transpose, shape=(3, 30757, 39900), dtype=uint16, chunksize=(3, 1024, 1024), chunktype=numpy.ndarray>,
 dask.array<transpose, shape=(3, 15378, 19950), dtype=uint16, chunksize=(3, 1024, 1024), chunktype=numpy.ndarray>,
 dask.array<transpose, shape=(3, 7689, 9975), dtype=uint16, chunksize=(3, 1024, 1024), chunktype=numpy.ndarray>,
 dask.array<transpose, shape=(3, 3844, 4987), dtype=uint16, chunksize=(3, 1024, 1024), chunktype=numpy.ndarray>,
 dask.array<transpose, shape=(3, 1922, 2493), dtype=uint16, chunksize=(3, 1024, 1024), chunktype=numpy.ndarray>,
 dask.array<transpose, shape=(3, 961, 1246), dtype=uint16, chunksize=(3, 961, 1024), chunktype=numpy.ndarray>,
 dask.array<transpose, shape=(3, 480, 623), dtype=uint16, chunksize=(3, 480, 623), chunktype=numpy.ndarray>,
 dask.array<transpose, shape=(3, 240, 311), dtype=uint16, chunksize=(3, 240, 311), chunktype=numpy.ndarray>]

So I guess somehow, my images are seen as 3-channel images.

Assembling mosaic  1/ 3 (channel  1/ 3): 100%|###################################| 14275/14275 [04:47<00:00, 49.59it/s]
Assembling mosaic  1/ 3 (channel  2/ 3): 100%|###################################| 14275/14275 [04:56<00:00, 48.22it/s]
Assembling mosaic  1/ 3 (channel  3/ 3): 100%|###################################| 14275/14275 [04:56<00:00, 48.11it/s]
Assembling mosaic  2/ 3 (channel  1/ 3): 100%|###################################| 24230/24230 [05:32<00:00, 72.81it/s]
Assembling mosaic  2/ 3 (channel  2/ 3): 100%|###################################| 24230/24230 [05:25<00:00, 74.54it/s]
Assembling mosaic  2/ 3 (channel  3/ 3): 100%|###################################| 24230/24230 [05:24<00:00, 74.72it/s]
Assembling mosaic  3/ 3 (channel  1/ 3): 100%|###################################| 23834/23834 [05:47<00:00, 68.53it/s]
Assembling mosaic  3/ 3 (channel  2/ 3): 100%|###################################| 23834/23834 [05:46<00:00, 68.76it/s]
Assembling mosaic  3/ 3 (channel  3/ 3): 100%|###################################| 23834/23834 [05:47<00:00, 68.65it/s]
2024-12-17 22:31:30.615 | INFO     | palom.pyramid:write_pyramid:182 - Generating pyramid
2024-12-17 22:31:30.615 | INFO     | palom.pyramid:write_pyramid:185 -     Level 1 (30758 x 39901)
Processing channel: 100%|################################################################| 9/9 [01:57<00:00, 13.04s/it]
2024-12-17 22:33:28.146 | INFO     | palom.pyramid:write_pyramid:185 -     Level 2 (15379 x 19951)
Processing channel: 100%|################################################################| 9/9 [00:36<00:00,  4.02s/it]
2024-12-17 22:34:04.411 | INFO     | palom.pyramid:write_pyramid:185 -     Level 3 (7690 x 9976)
Processing channel: 100%|################################################################| 9/9 [00:15<00:00,  1.74s/it]
2024-12-17 22:34:20.349 | INFO     | palom.pyramid:write_pyramid:185 -     Level 4 (3845 x 4988)
Processing channel: 100%|################################################################| 9/9 [00:10<00:00,  1.22s/it]
2024-12-17 22:34:31.365 | INFO     | palom.pyramid:write_pyramid:185 -     Level 5 (1923 x 2494)
Processing channel: 100%|################################################################| 9/9 [00:09<00:00,  1.09s/it]
2024-12-17 22:34:41.161 | INFO     | palom.pyramid:write_pyramid:185 -     Level 6 (962 x 1247)
Processing channel: 100%|################################################################| 9/9 [00:09<00:00,  1.06s/it]
2024-12-17 22:34:50.687 | INFO     | palom.pyramid:write_pyramid:185 -     Level 7 (481 x 624)
Processing channel: 100%|################################################################| 9/9 [00:09<00:00,  1.04s/it]

The registration of the 9 channels is just fantastic. So I have the 3 DAPIs + 6 other channels. Unfortunately, I can't distinguish which is which. Is there a way to add the channel names to the final mosaic OME-TIFF, either by supplying a csv file with the channel names and order or by querying the file name? This will be just awesome!

I will send you my images (6 in total from 6 staining cycles) as well as the jupyter notebook with the output of all blocks on the email linked to your profile.