Is it possible to use different --align-channel for different cycles?

[ciszew]

In some of our experiments we have different number of channels for some of the cycles. This results in mismatched numbering in the order of the channels, for example:
cycle 1 has DAPI only -> DAPI channel=0
cycle 2 has 4 colors with DAPI as last channel -> DAPI channel=3
cycle 3 has 5 colors with DAPI as last channel-> DAPI channel=4
...more cycles like cycle 2 and 3
when setting up ashlar run we specify align channel option, in above scenario if selecting channel 0 for alignment some of the cycles (small number 0-2 cycles out of 15 cycles in recent experiment) have a large shift in alignment in range of hundreds to few thousands um. Image for this cycle appears cut or shifted by large distance.
Im guessing this is due to the different channels being used for alignment for different cycles (in above example channel 0 is DAPI for cycle 1 and Cy7 for cycle 2 and 3). If my assumption is correct is there a way to specify different alignment channels for different cycles (in example above 0, 3, 4 to use DAPI signal each time).
I tested alignment of selected matching cycles with the same true DAPI channel and i didnt get the large shift described above.
When searching for some answerers i noticed that issue #153 has similar large shift mentioned, possibly its related?

[jmuhlich]

I understand the problem! The command line interface doesn't support this but the internal Python code does. You could write a short custom script to call the Python API directly and choose separate channel numbers for each cycle. If you could provide a sample ashlar command line that you use, I can convert that to a Python script where you can customize the channel numbers per cycle.

[ciszew]

Thank you for creating this great tool and thank you for continuous support.
It may not be the most elegant solution but what we are using right now is a bash script that ties ims to ome.tif conversion using bfconvert tool and then stitching and registration with ashlar (this 2 are most time consuming steps and its nice to run them together in automated way). Example of relevant ashlar part of the script looks like this:

WIDTH=(3)
HEIGHT=(4)
PIXEL_SIZE=0.157
OVERLAP=0.1
LAYOUT=snake
DIRECTION=vertical
export MAX_SHIFT=30
export SIGMA=2
export ALIGN_CHANNEL=4
export FILENAME_HEADER="mcmicro_test_40x
 c+=( "fileseries|${OUTPUTDIR}|pattern=${FILENAME_HEADER}${k}_F{series}.ome.tiff|width=${WIDTH[$(($(($i))-1))]}|height=${HEIGHT[$(($(($i))-1))]}|pixel_size=${PIXEL_SIZE}|overlap=${OVERLAP}|layout=${LAYOUT}|direction=${DIRECTION}" )
ashlar "${c[@]}" --flip-y --pyramid --maximum-shift $MAX_SHIFT --filter-sigma $SIGMA --align-channel $ALIGN_CHANNEL -o $OUTPUT_DIR"image_"$i"".ome.tiff

[ciszew]

I have similar question regarding --output-channels argument so decided to piggyback in this thread. Is it possibly to specify to output different channels per cycle? or even easier to specify channels that should be omitted from some of the cycles (for example if we dont have 4 antibodies in one of the cycles we often run the same acquisition protocol which generates one of the channels as intermediate background, it simplifies the setup, keeps all the channels in order and give us a little bit of QC window for signal quenching between actual antibodies) but for final stack image file for analysis it would be nice to have an option to remove this "empty" channels if possible.

[josenimo]

hello @ciszew, did you figure out how to use different channel numbers, I am in a similar situation as you were.
if you could share the solution, I would be very grateful!
Thanks,
Jose

[Yu-AnChen]

This currently cannot be achieved using the ashlar script/command. For a more complex configuration @ciszew described, I think constructing a configuration file will be needed.

Channel used for alignment is set when instantiating reg.EdgeAligner (example) and reg.LayerAligner (example) by passing channel=N (N is an integer >= 0; 0-based indexing) and output channels in the output ome-tiff are specified when instantiating reg.Mosaic (example), channels=[1, 2, 3] to select second, third, and fourth channels, for example.

[ciszew]

Hi @josenimo
I only have an easy work-arounds right now. We either acquire images so DAPI channel that we use for ashlar alignment is always 1st in each cycle (align channel=0 for every cycle regardless number of acquired channels in each cycle) or with some investigators we always acquire the same number of channels (5) for all cycles with DAPI being last (align channel=4) which creates some "empty channels"
I have a collaborator who is finalizing data processing pipeline and they are using several python packages to convert data from our microscope (ims format to tiff) and in the process they can also rearrange and and remove some of the channels. Im guessing this is what we will use when finalized.
Or possibly, with new version of ashlar we will be able to go directly from individual ims files to stack and with an option to choose different channels for alignment :)

[josenimo]

Thank you both,

@Yu-AnChen I will try again to take a closer look at this file, I believe I can manage to understand it with my poor python, but right now I am under a lot of pressure to present some results soon :`). Thank you for the help.

@ciszew this is what I would have done have every channel start at 0, or image empty. However this was our first run at CYCIF and they always imaged DAPI last, and cycles have different number of channels.

My current hack is the following, read the image in python with skimage.io.imread() and then reverse the numpy array with numpy.flip(array, axis=0), and then exporting with imsave(). This works well in theory, always switching DAPI from the last channel to the first channel. However I have few concerns:

1. This effectively converts my .ome.tif images into .tif images, their sizes fell down from ~85gb to ~60gb. I am not 100% sure what their difference is, I assume it is the tiling metadata present in the .ome.tif. I am afraid ashlar cannot run properly with these 60,000x80,000px images without tile information.
2. I am afraid my workstation cannot handle the final result, I am stacking 4 cycles, each averaging 80gb, this sums up to more than my 256gb of RAM. Could I even visualize the output in a sensible way, with napari or Qupath?

Any ideas and feedback are very welcome!
Best,
Jose

[Yu-AnChen]

@josenimo thanks for elaborating on your use case. I added a branch on my fork that allows passing channel -1 when calling ashlar and it's indexing from the end of the sequence, i.e. -1 indicates the last channel.

Not sure if you are running ashlar within mcmicro pipeline, if so you likely will need to build a custom docker container. Otherwise I would pull the code from the branch and pip install it in editable mode

git clone https://github.com/labsyspharm/ashlar.git
cd ashlar
git remote add yu-anchen https://github.com/yu-anchen/ashlar.git
git fetch yu-anchen
git checkout -b last-channel yu-anchen/last-channel

# install in editable mode
pip install -e .

Hopefully, with this you don't have to re-write the files.

On you two other concerns

1. It seems to me that you are using fileseries/filepattern reader, i.e. the position of each FOV is denoted in the filename and you specify the scanning grid spec. If so, I wouldn't worry too much about the metadata in the original metadata other than image pixel size. If the stage positions are recorded in the ome-tiffs, then maybe using a companion file and bioformats reader is better.
2. ashlar writes tiled multi-resolution pyramidal ome-tiff, and it can be load lazily, meaning you wouldn't need to load the whole image into RAM. Both QuPath and Napari works, avivator is another easy solution too.

[josenimo]

Dear @Yu-AnChen , thank you for your time, explanation, and sharing your ashlar fork. I have run into another issue with ashlar and I will try to explain, erring on the side of oversharing some details, thank you for bearing with me, I have learnt lots in the past few weeks.

Context:
My colleague imaged 4 cycles of a lung cancer with strong autofluorescence, she imaged DAPI always as the last channel.

Preprocessing:
In order to remove the autofluorescence, we imaged before staining with antibodies for every cycle, with the same microscope settings. To use the RAM intensive Schapiro's background subtraction I decided to register and stitch each autofluorescent cycle with the antibody staining of each cycle, then I ran the background subtraction for every cycle individually.

Right now I have 4 stacks, one for each cycle. I tried your ashlar fork unsuccessfully, see stack_trace below:

(ashlar_yuanchen) E:\Jose\bg_sub>ashlar --output ./Slide1.ome.tif -c -1 --filter-sigma 1 -m 30 bg_sub_Slide1_Cycle1.ome.tif bg_sub_Slide1_Cycle2.ome.tif bg_sub_Slide1_Cycle3.ome.tif bg_sub_Slide1_Cycle4.ome.tif
Stitching and registering input images
Cycle 0:
    reading bg_sub_Slide1_Cycle1.ome.tif
WARNING: Stage coordinates undefined; falling back to (0, 0).
Traceback (most recent call last):
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar_yuanchen\Scripts\ashlar-script.py", line 33, in <module>
    sys.exit(load_entry_point('ashlar', 'console_scripts', 'ashlar')())
  File "c:\users\administrator\jose_bi\ashlar\ashlar\scripts\ashlar.py", line 212, in main
    return process_single(
  File "c:\users\administrator\jose_bi\ashlar\ashlar\scripts\ashlar.py", line 243, in process_single
    edge_aligner.run()
  File "c:\users\administrator\jose_bi\ashlar\ashlar\reg.py", line 484, in run
    self.make_thumbnail()
  File "c:\users\administrator\jose_bi\ashlar\ashlar\reg.py", line 495, in make_thumbnail
    self.reader.thumbnail = thumbnail.make_thumbnail(
  File "c:\users\administrator\jose_bi\ashlar\ashlar\thumbnail.py", line 13, in make_thumbnail
    coordinate_max = (positions + metadata.size).max(axis=0)
  File "c:\users\administrator\jose_bi\ashlar\ashlar\reg.py", line 107, in size
    raise ValueError("Image series must all have the same dimensions")
ValueError: Image series must all have the same dimensions

I checked the dimensions of the four files, with imread and they looked like this:
Cycle1: (5, 62905, 86776)
Cycle2: (62906, 86782, 4)
Cycle3: (62904, 86780, 3)
Cycle4: (62904, 86780, 3)

I first tried transposing all of them to xyc using numpy.transpose, then cropping all of them to the same number of x,y pixels by slicing the array, and flipping the channel order with numpy.flip(array,axis=0) , writing with imsave() and running ashlar again, another error:

(ashlar) E:\Jose\bg_sub>ashlar --output ./Slide1.ome.tif -c 0 --filter-sigma 1 -m 30 crop_flip_Slide1_Cycle1.ome.tif crop_flip_Slide1_Cycle2_m30.ome.tif crop_flip_Slide1_Cycle3.ome.tif crop_flip_Slide1_Cycle4.ome.tif
Stitching and registering input images
Cycle 0:
    reading crop_flip_Slide1_Cycle1.ome.tif
WARNING: Stage coordinates undefined; falling back to (0, 0).
WARNING: Pixel size undefined; falling back to 1.0 μm.
    assembling thumbnail 1/1Traceback (most recent call last):
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar\lib\runpy.py", line 196, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar\lib\runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar\Scripts\ashlar.exe\__main__.py", line 7, in <module>
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar\lib\site-packages\ashlar\scripts\ashlar.py", line 212, in main
    return process_single(
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar\lib\site-packages\ashlar\scripts\ashlar.py", line 243, in process_single
    edge_aligner.run()
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar\lib\site-packages\ashlar\reg.py", line 482, in run
    self.make_thumbnail()
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar\lib\site-packages\ashlar\reg.py", line 493, in make_thumbnail
    self.reader.thumbnail = thumbnail.make_thumbnail(
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar\lib\site-packages\ashlar\thumbnail.py", line 20, in make_thumbnail
    img = reader.read(c=channel, series=i)
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar\lib\site-packages\ashlar\reg.py", line 427, in read
    img = self.reader.read(series, c)
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar\lib\site-packages\ashlar\reg.py", line 404, in read
    byte_array = self.metadata._reader.openBytes(index)
  File "jnius\jnius_export_class.pxi", line 1177, in jnius.JavaMultipleMethod.__call__
  File "jnius\jnius_export_class.pxi", line 885, in jnius.JavaMethod.__call__
  File "jnius\jnius_export_class.pxi", line 975, in jnius.JavaMethod.call_method
  File "jnius\jnius_utils.pxi", line 91, in jnius.check_exception
jnius.JavaException: JVM exception occurred: Array size too large: 85000 x 60000 x 2 java.lang.IllegalArgumentException

I don't understand what array too large means I did run ashlar on bigger images, so I thought ASHLAR uses yxc instead of cyx, based of Issue 104 and Issue 95.

So I then transposed cycle one to yxc, (without channel reorder, and without pixel slicing) and tried running ashlar_yuanchen again, still an error:

(ashlar_yuanchen) E:\Jose\bg_sub>ashlar --output ./Slide1.ome.tif -c -1 --filter-sigma 1 -m 30 cycle1_xyc.ome.tif bg_sub_Slide1_Cycle2.ome.tif bg_sub_Slide1_Cycle3.ome.tif bg_sub_Slide1_Cycle4.ome.tif
Stitching and registering input images
Cycle 0:
    reading cycle1_xyc.ome.tif
WARNING: Stage coordinates undefined; falling back to (0, 0).
WARNING: Pixel size undefined; falling back to 1.0 μm.
    assembling thumbnail 1/1

Cycle 1:
    reading bg_sub_Slide1_Cycle2.ome.tif
WARNING: Stage coordinates undefined; falling back to (0, 0).
Traceback (most recent call last):
  File "C:\Users\Administrator\Jose_BI\Anaconda\envs\ashlar_yuanchen\Scripts\ashlar-script.py", line 33, in <module>
    sys.exit(load_entry_point('ashlar', 'console_scripts', 'ashlar')())
  File "c:\users\administrator\jose_bi\ashlar\ashlar\scripts\ashlar.py", line 212, in main
    return process_single(
  File "c:\users\administrator\jose_bi\ashlar\ashlar\scripts\ashlar.py", line 259, in process_single
    layer_aligner.run()
  File "c:\users\administrator\jose_bi\ashlar\ashlar\reg.py", line 813, in run
    self.make_thumbnail()
  File "c:\users\administrator\jose_bi\ashlar\ashlar\reg.py", line 819, in make_thumbnail
    self.reader.thumbnail = thumbnail.make_thumbnail(
  File "c:\users\administrator\jose_bi\ashlar\ashlar\thumbnail.py", line 13, in make_thumbnail
    coordinate_max = (positions + metadata.size).max(axis=0)
  File "c:\users\administrator\jose_bi\ashlar\ashlar\reg.py", line 107, in size
    raise ValueError("Image series must all have the same dimensions")
ValueError: Image series must all have the same dimensions

When I tried slicing the pixels for the yxc arrays I keep running out of RAM:

Output exceeds the [size limit](command:workbench.action.openSettings?[). Open the full output data [in a text editor](command:workbench.action.openLargeOutput?c8ddcfc1-e556-4b33-9a82-3b080455dff2)
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
Cell In[4], line 1
----> 1 import_crop_export('bg_sub_Slide1_Cycle2.ome.tif','yxc_DAPIlast_Cycle2.ome.tif')

Cell In[3], line 20, in import_crop_export(import_path, export_path)
     14 qc_image(im_crop)
     16 #im_crop_flip = np.flip(im_crop, axis=2)
     17 #print('crop_flipped')
     18 #qc_image(im_crop_flip)
---> 20 imsave(export_path, im_crop)

File ~\Jose_BI\Anaconda\envs\Data_v2\lib\site-packages\skimage\io\_io.py:141, in imsave(fname, arr, plugin, check_contrast, **plugin_args)
    137     warn('%s is a boolean image: setting True to 255 and False to 0. '
    138          'To silence this warning, please convert the image using '
    139          'img_as_ubyte.' % fname, stacklevel=2)
    140     arr = arr.astype('uint8') * 255
--> 141 if check_contrast and is_low_contrast(arr):
    142     warn('%s is a low contrast image' % fname)
    143 return call_plugin('imsave', fname, arr, plugin=plugin, **plugin_args)

File ~\Jose_BI\Anaconda\envs\Data_v2\lib\site-packages\skimage\exposure\exposure.py:831, in is_low_contrast(image, fraction_threshold, lower_percentile, upper_percentile, method)
    828 from ..color import rgb2gray, rgba2rgb  # avoid circular import
    830 if image.shape[2] == 4:
--> 831     image = rgba2rgb(image)
...
--> 223 out = np.clip((1 - alpha) * background + alpha * channels,
    224               a_min=0, a_max=1)
    225 return out

MemoryError: Unable to allocate 114. GiB for an array with shape (60000, 85000, 3) and data type float64

Right now I am trying to recreate the background subtraction preprocessing again for cycle1, assuming something is off. but the output is still (c,y,x).

Regarding my previous concerns:

1. I am runnning the entire stack directly into ashlar, no series were used, I though ashlar used the metadata to figure out pixel size and microscope tile coordenates.
2. This is great to know, so technically I could open and explore the stacks.

If there is any simple way to slice my ome.tif files and keep the metadata ? (Right now I just know how to use skimage.io)
Or is there a way for ashlar to be more lenient with small shape differences?
Any feedback and comments are greatly appreciated!
Thank you for your help

[Yu-AnChen]

Apologies I think we are deviating from the topic of the thread. Maybe we could talk about the bg subtraction issues on image.sc.

On the ashlar side, based on what you've described, I would run ashlar on all the scans collected (w/ both Ab stains and bg-only), assuming your imager provides access to the raw FOVs. After running ashlar, it outputs a single file that contains all the channels, and the bg subtraction module uses that as input to generate the bg-subtracted image.

[josenimo]

I agree,
thank you both again for the help!

[joaomamede]

Try this in a notebook, worked for me, in my case I have 0 overlap between tiles in this set:

from ashlar import fileseries, thumbnail,reg
import matplotlib.pyplot as plt
from ashlar.scripts.ashlar import process_axis_flip
import numpy as np
# import pims

fnames = {}
n_cycles = 2



folder = '/ogmios/Janet/20230221_microglia_n1/fixed/'
fnames[0] = 'Cycle0-PQBP1647_caru_IN-cgas405_000_v{series}_PRJ.ome.tif'
fnames[1] = 'Cycle1-DAPI-IRF3488_IBA555_cgas647_v{series}_PRJ.ome.tif'
# fnames[2] = 'Cycle2-dapi-NA-CD4_PE-PRY19_APC_v{series}_PRJ.ome.tif'


readers = []

#8,9
for i in range(n_cycles):
    readers.append(fileseries.FileSeriesReader(
        folder,
        pattern=fnames[i],
        overlap=0.001,
        width=9,#7,
        height=4,#6,
        layout='snake',
        direction='horizontal',
        pixel_size=0.1072005382,
#         0.16258744,
    )
                  )
    readers[i].thumbnail = thumbnail.make_thumbnail(readers[i], channel=0)





plt.figure()
# doing log just for visualization
plt.imshow(np.log(readers[0].thumbnail))
# process_axis_flip(c1r, flip_x=False, flip_y=True)
edge_aligners = []
channel_dict = {0:3,1:0,2:3}
for i in [0]:
#     range(n_cycles):
    edge_aligners.append(reg.EdgeAligner(readers[0], max_shift=75, channel=channel_dict, filter_sigma=0.05, verbose=True,do_make_thumbnail=True))
    edge_aligners[i].run()
reg.plot_edge_quality(edge_aligners[0], img=edge_aligners[0].reader.thumbnail)
mosaics = []
mosaics.append( reg.Mosaic(
        aligner = edge_aligners[0],
        shape = edge_aligners[0].mosaic_shape,
#         channels=range(4),
    #'/home/jmamede/Data/mglCBO/2/C0-ch{channel}-DAPI-GFP-cGAMP-p241D_mosaic.ome.tif',
#     **mosaic_args
    )
             )

for j in range(1,n_cycles):
# for j in range(1,1):
    layer_aligner = reg.LayerAligner(readers[j], edge_aligners[0], max_shift=350, channel=channel_dict[j], filter_sigma=0.05, verbose=True)

    layer_aligner.run()
    mosaics.append( reg.Mosaic(
        layer_aligner, edge_aligners[0].mosaic_shape,
#         **mosaic_args
    )
                 )

output_path_format = folder+'/AllCycles.ome.tif'
quiet = False
pyramid = True
writer_args = {}
if pyramid: 
#if args.output_channels:
    writer_args['scale'] = 2
#pyramid
    writer_args['tile_size'] = 4096
if not quiet:
    print()
    print(f"Merging tiles and writing to {output_path_format}")
writer_class = reg.PyramidWriter if pyramid else reg.TiffListWriter
writer = writer_class(
    mosaics, output_path_format, verbose=not quiet, **writer_args
)
writer.run()

[joaomamede]

I have another, more complicated question.
Is it possible to align to different patterns between cycles?

Cycle 0, I don't have dapi but I have a Protein A stain
Cycle 1, I have dapi and another cycle of Protein A stain
Cycle 2, I have dapi and all different proteins?

Would I have to create a 0+1 and then align to the dapi channel? Problem is, I need to split the image in at least two (to then mosaic again since Ashlar doesn't do pre-stiched images)? I've done a fake mosaic before with Ashlar, that wouldn't be the biggest problem.