Example Stitching Notebook

Contents

Example Stitching Notebook#

scPortrait uses Ashlar for stitching images. When stitching from .tif files, Ashlar reads channel and tile position information from filenames according to a predefined pattern. Hence, filenames matter when stitching from .tif files.

[1]:

import os

os.environ["JAVA_HOME"] = "/Users/sophia/mambaforge/envs/scPortrait/lib/jvm"
import matplotlib.pyplot as plt

from scportrait.data._datasets import dataset_stitching_example
from scportrait.tools.stitch import ParallelStitcher, Stitcher

/Users/sophia/mambaforge/envs/scportrait_new/lib/python3.10/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.
  warnings.warn(
/Users/sophia/mambaforge/envs/scportrait_new/lib/python3.10/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Single-threaded Stitching#

Initializing the `Stitcher` object#

[ ]:

input_dir = dataset_stitching_example()
slidename = "stitching_test"
outdir = os.path.join(str(input_dir).replace("stitching_example", "example_projects/stitching"), slidename)

row = str(2).zfill(2)  # specify the row of the well you want to stitch, here = 2
well = str(4).zfill(2)  # specifc the well number you wish to stitch, here = 4
zstack_value = str(1).zfill(
    3
)  # specify the zstack you want to stitch. for multiple zstacks please make a loop and iterate through each of them.
timepoint = str(1).zfill(3)  # specify the timepoint you wish to stitch

pattern = f"Timepoint{timepoint}_Row{row}_Well{well}_{{channel}}_zstack{zstack_value}_r{{row:03}}_c{{col:03}}.tif"

# initialize stitcher
stitcher = Stitcher(
    os.path.abspath(input_dir),
    slidename,
    outdir,
    pattern=pattern,
    stitching_channel="Alexa488",
    channel_order=[
        "DAPI",
        "Alexa488",
        "mCherry",
    ],  # this can be used to override the order of channels in the final stitched image for image types like ome.zarr, if not specified alphabetical order is used
    overlap=0.1,
    max_shift=30,
    filter_sigma=0,
    rescale_range={"Alexa488": (1, 99), "DAPI": (1, 99), "mCherry": (1, 99)},
    overwrite=True,
)

Output directory at /Users/sophia/Documents/GitHub/scPortrait/scportrait_data/example_projects/stitching/stitching_test already exists, overwriting.

You can access information on which channels are present in your dataset and the specific values of the parameters relevant to stitching by executing:

[22]:

stitcher.get_stitching_information()

Tile positions will be calculated based on channel: Alexa488
Channel Names: ['Alexa488', 'DAPI', 'mCherry']
Overlap of image tiles: 0.1
Max Shift value: 30
Filter Sigma value: 0
Output will be written to: /Users/sophia/Documents/GitHub/scPortrait/scportrait_data/example_projects/stitching/stitching_test

Generating thumbnails#

[3]:

stitcher.generate_thumbnail()

    assembling thumbnail 9/9

[4]:

# thumbnail is saved in the stitcher object and can be accessed via stitcher.thumbnail
plt.imshow(stitcher.thumbnail)

[4]:

<matplotlib.image.AxesImage at 0xb08b0b580>

../../../_images/pages_tools_stitching_example_stitching_notebook_10_1.png

[5]:

# alterantively it can be saved to a tif file
stitcher.write_thumbnail()

Generating full-scale stitched image#

[6]:

stitcher.stitch()

performing stitching on channel Alexa488 with id number 0
    quantifying alignment error 1000/1000
    aligning edge 12/12

WARNING:root:Folder /var/folders/35/p4c58_4n3bb0bxnzgns1t7kh0000gn/T/temp_mmap_rajh1lzh with temp mmap arrays is being deleted. All existing temp mmapp arrays will be unusable!
WARNING:root:New temp folder location. Temp mmap arrays are written to /Users/sophia/Documents/GitHub/scPortrait/scportrait_data/example_projects/stitching/stitching_test/temp_mmap_zhxh9_b2. Cleanup of this folder is OS dependant, and might need to be triggered manually!

Alignment complete.
assembling mosaic with shape (3, 3040, 3038)
created tempmmap array for assembled mosaic at /Users/sophia/Documents/GitHub/scPortrait/scportrait_data/example_projects/stitching/stitching_test/temp_mmap_zhxh9_b2/temp_mmap_267690355536480897.hdf

  0%|          | 0/3 [00:00<?, ?it/s]

        merging tile 9/9

 33%|███▎      | 1/3 [00:01<00:02,  1.06s/it]


        merging tile 9/9

 67%|██████▋   | 2/3 [00:02<00:01,  1.06s/it]


        merging tile 9/9

100%|██████████| 3/3 [00:03<00:00,  1.06s/it]

../../../_images/pages_tools_stitching_example_stitching_notebook_13_12.png

../../../_images/pages_tools_stitching_example_stitching_notebook_13_13.png

[7]:

# the stitched image is saved in the stitcher object and can be accessed via stitcher.assembled_mosaic
stitcher.assembled_mosaic

[7]:

	Array	Chunk
Bytes	52.85 MiB	52.85 MiB
Shape	(3, 3040, 3038)	(3, 3040, 3038)
Dask graph	1 chunks in 2 graph layers
Data type	uint16 numpy.ndarray

[8]:

# the stitched image can then be written to a variety of output formats

stitcher.write_tif(export_xml=True)
stitcher.write_ome_zarr()

[9]:

import matplotlib.pyplot as plt
from tifffile import imread

fig, axs = plt.subplots(1, 3, figsize=(30, 10))
axs[0].imshow(imread(f"{outdir}/stitching_test_Alexa488.tif"))
axs[0].axis("off")

axs[1].imshow(imread(f"{outdir}/stitching_test_DAPI.tif"))
axs[1].axis("off")

axs[2].imshow(imread(f"{outdir}/stitching_test_mCherry.tif"))
axs[2].axis("off")

[9]:

(-0.5, 3037.5, 3039.5, -0.5)

../../../_images/pages_tools_stitching_example_stitching_notebook_16_1.png

[10]:

del stitcher

Multi-threaded Stitching#

The ParallelStitcher class can speed up stitching by using multiple threads. The code to start stitching remains the same, but ParallelStitcher takes an additional argument threads, specifying the number of parallel threads to use.

Initializing the `ParallelStitcher` object#

[11]:

input_dir = dataset_stitching_example()
slidename = "stitching_test_parallel"
outdir_parallel = os.path.join(str(input_dir).replace("stitching_example", "example_projects/stitching"), slidename)


row = str(2).zfill(2)  # specify the row of the well you want to stitch, here = 2
well = str(4).zfill(2)  # specifc the well number you wish to stitch, here = 4
zstack_value = str(1).zfill(
    3
)  # specify the zstack you want to stitch. for multiple zstacks please make a loop and iterate through each of them.
timepoint = str(1).zfill(3)  # specify the timepoint you wish to stitch

pattern = f"Timepoint{timepoint}_Row{row}_Well{well}_{{channel}}_zstack{zstack_value}_r{{row:03}}_c{{col:03}}.tif"

# initialize stitcher
stitcher = ParallelStitcher(
    input_dir,
    slidename,
    outdir_parallel,
    pattern=pattern,
    stitching_channel="Alexa488",
    channel_order=["DAPI", "Alexa488", "mCherry"],
    overlap=0.1,
    max_shift=30,
    filter_sigma=0,
    rescale_range={"Alexa488": (1, 99), "DAPI": (1, 99), "mCherry": (1, 99)},
    overwrite=True,
    threads=12,
)

Output directory at /Users/sophia/Documents/GitHub/scPortrait/scportrait_data/example_projects/stitching/stitching_test_parallel already exists, overwriting.

Generating thumbnails#

[12]:

stitcher.generate_thumbnail()

    assembling thumbnail 9/9

[13]:

# thumbnail is saved in the stitcher object and can be accessed via stitcher.thumbnail
plt.imshow(stitcher.thumbnail)

[13]:

<matplotlib.image.AxesImage at 0xb61216c80>

../../../_images/pages_tools_stitching_example_stitching_notebook_23_1.png

[14]:

# alterantively it can be saved to a tif file
stitcher.write_thumbnail()

Generating full-scale stitched image#

[15]:

stitcher.stitch()

performing stitching on channel Alexa488 with id number 0
    quantifying alignment error: 100%|██████████| 1000/1000 [00:01<00:00, 935.36it/s]
                  aligning edge: 100%|██████████| 12/12 [00:00<00:00, 451.91it/s]
using graph-tool to build spanning tree

WARNING:root:Folder /Users/sophia/Documents/GitHub/scPortrait/scportrait_data/example_projects/stitching/stitching_test/temp_mmap_zhxh9_b2 with temp mmap arrays is being deleted. All existing temp mmapp arrays will be unusable!
WARNING:root:New temp folder location. Temp mmap arrays are written to /Users/sophia/Documents/GitHub/scPortrait/scportrait_data/example_projects/stitching/stitching_test_parallel/temp_mmap_yq7657gu. Cleanup of this folder is OS dependant, and might need to be triggered manually!

Alignment complete.
current channel order:  [0, 1, 2]
new channel order [1, 0, 2]
assembling mosaic with shape (3, 3040, 3038)
created tempmmap array for assembled mosaic at /Users/sophia/Documents/GitHub/scPortrait/scportrait_data/example_projects/stitching/stitching_test_parallel/temp_mmap_yq7657gu/temp_mmap_2248175831853828840.hdf
assembling channels with 3 workers
assembling mosaic:   0%|          | 0/3 [00:00<?, ?it/s]



























assembling channel 1: 100%|██████████| 9/9 [00:01<00:00,  7.58it/s]
assembling channel 0: 100%|██████████| 9/9 [00:01<00:00,  7.54it/s]
assembling mosaic:  33%|███▎      | 1/3 [00:01<00:02,  1.20s/it]

assembling channel 2: 100%|██████████| 9/9 [00:01<00:00,  7.53it/s]
assembling mosaic: 100%|██████████| 3/3 [00:01<00:00,  2.50it/s]

../../../_images/pages_tools_stitching_example_stitching_notebook_26_3.png

../../../_images/pages_tools_stitching_example_stitching_notebook_26_4.png

[16]:

# the stitched image is saved in the stitcher object and can be accessed via stitcher.assembled_mosaic
stitcher.assembled_mosaic

[16]:

	Array	Chunk
Bytes	52.85 MiB	52.85 MiB
Shape	(3, 3040, 3038)	(3, 3040, 3038)
Dask graph	1 chunks in 2 graph layers
Data type	uint16 numpy.ndarray

[17]:

# the stitched image can then be written to a variety of output formats

stitcher.write_tif(export_xml=True)

# alternatively the stitched images can also be written out to tifs in a multi-threaded format
# speedups here are limited by write speed to disk, the number of threads is limited by the number of channels available
stitcher.write_tif_parallel(export_xml=True)

stitcher.write_ome_zarr()

writing tif files: 100%|██████████| 3/3 [00:00<00:00, 221.21it/s]

Visualize Stitching Output#

[18]:

import matplotlib.pyplot as plt
from tifffile import imread

fig, axs = plt.subplots(1, 3, figsize=(30, 10))
axs[0].imshow(imread(f"{outdir_parallel}/stitching_test_parallel_Alexa488.tif"))
axs[0].axis("off")

axs[1].imshow(imread(f"{outdir_parallel}/stitching_test_parallel_DAPI.tif"))
axs[1].axis("off")

axs[2].imshow(imread(f"{outdir_parallel}/stitching_test_parallel_mCherry.tif"))
axs[2].axis("off")

[18]:

(-0.5, 3037.5, 3039.5, -0.5)

../../../_images/pages_tools_stitching_example_stitching_notebook_30_1.png

[19]:

# compare parallel and non-parallel stitching

imread(f"{outdir_parallel}/stitching_test_parallel_Alexa488.tif") == imread(f"{outdir}/stitching_test_Alexa488.tif")

[19]:

array([[ True,  True,  True, ...,  True,  True,  True],
       [ True, False, False, ..., False, False,  True],
       [ True, False, False, ..., False, False,  True],
       ...,
       [ True, False, False, ..., False, False,  True],
       [ True, False, False, ..., False, False,  True],
       [ True,  True,  True, ...,  True,  True,  True]])

[20]:

del stitcher