Registration

Final data is in one common coordinate system. It takes quite a few steps to get that. This file describe how genes and barcodes are registered together.

Overview

The command is run as follows:

iss register --path relative/path/to/data --prefix genes_round

By default, the script will perform only missing steps. There is a --force-redo flag that will force the script to re-run all steps, even if the output files already exist.

The script will perform different steps depending on the type of acquisition. See details for sequencing and fluorescent acquisitions.

To check that everything worked, look in figures/registration/prefix for diagnostic plots. If there is no motion on the mp4 files for sequencing acquisitions, you’re good. For fluorescent acquisitions, you should see that the spots are aligned across channels on the png.

If it looks bad, see the relevant Troubleshooting section below.

Register Sequencing acquisitions

These are acquisitions that have multiple rounds, namely genes_round and barcode_round. The steps are:

        flowchart TD
start[Start] --> regref[run_register_reference_tile];
    batch_est(((register_tile)));

    subgraph register_reference_tile
        regref --> diag_ref([check_ref_tile_registration]);
    end



    regref --> batch_est;
    batch_est --> corr;
    subgraph correct_shifts
        subgraph run_correct_shifts
            corr[run_correct_shifts];
            corr --> filt[filter_ransac_shifts];
        end

        subgraph check_tile_shifts
            filt --> ctr([check_tile_registration]);
        end
        subgraph check_shift_correction
            filt --> csc([check_shift_correction]);
        end
        subgraph check_tile_registration
            ctr --> diag_tile([check_tile_registration]);
        end
    end

    style batch_est fill:#E1BEE7,stroke:#424242

    style corr stroke:#000000,fill:#E1BEE7
    style regref stroke:#000000,fill:#E1BEE7
    style filt stroke:#000000,fill:#E1BEE7

    style csc fill:#BBDEFB,stroke:#616161,color:#000000
    style diag_ref fill:#BBDEFB,stroke:#616161,color:#000000
    style diag_tile fill:#BBDEFB,stroke:#616161,color:#000000
    style ctr fill:#BBDEFB,stroke:#616161,color:#000000

    style run_correct_shifts fill:#EEEEEE, stroke:#424242
    style check_tile_shifts fill:#EEEEEE, stroke:#424242
    style check_shift_correction fill:#EEEEEE, stroke:#424242
    style check_tile_registration fill:#EEEEEE, stroke:#424242

    style correct_shifts fill:#AAAAAA, stroke:#424242
    style register_reference_tile fill:#AAAAAA, stroke:#424242

Troubleshooting

If the registration looks bad, we need to find which step failed.

  • Is the reference tile registered properly?

In the figures/registration folder, look at the files starting with registration_reference_tile. If there is no or little signal: pick a better tile (change the 'ref_tile' parameter in the ops.yml file). If there is signal and it still looks bad, double check the figures/registration/PREFIX folder, look at the affine_debug_PREFIX... png file. You might not have enough signal for affine registration. If that’s the case, you will have to try similarity transform (not supported anymore, but might still work).

  • Are most shifts estimated correctly?

Parameters you can tweak: ops["ransac_max_shift"], ops["ransac_min_tiles"], and ops["ransac_residual_threshold"].

Registering sequencing rounds

We need to register the channels and rounds together and the tiles with their neighbours.

Short version:

Final transforms for channel and round registration are saved in "reg" / f"tforms_best_{prefix}_{roi}_{tilex}_{tiley}.npz". The second part must be run each time.

Detailed explanation part 1: Estimating angle, scale and shift

For each acquisition we need to find how the channels register together. It needs to be done for each acquisition as the mirrors wobble a bit and the gain of the stage motor seem to vary a bit.

Register reference tile

We do that on a manually tile that has signal with iss register-ref-tile.

This will save f"tforms_{prefix}.npz" in the main data_folder. The npz contains:

  • angles_within_channels: rotation angles between rounds for each channel

  • shifts_within_channels: shifts between rounds for each channel

  • scales_between_channels: scaling between channels (common for all rounds)

  • angles_between_channels: rotation angles between channels (common for all rounds)

  • shifts_between_channels: shifts between channels (common for all rounds)

To estimate these values, the algorithm first align images for each channel across rounds. This is much more reliable than registering different channels for the same acquisition, as the sequencing dyes have limited bleedthrough across channels. On the other hand, when aligning between rounds, many rolonies will have the same base and therefore show up across rounds, providing a robust signal for registration.

Registration is done by iterative grid search. We first search over an initial range of rotation angles and compute phase correlation for each angle. We then determine the best angle and narrow the search range around this value. It is important that the initial spacing between angles is fine enough that we can find this peak. This will yield angles_within_channels and shifts_within_channels.

Once we have registered together rounds for each channel, we can use the resulting angles and shifts to compute mean and STD projections across rounds (we use the STD projections because rolonies show up very nicely on them). These projection should capture all rolonies and will look very similar across channels. They are provide ideal signal for registration across channels.

To register channels we need to correct for scaling as well as rotation due to chromatic aberration and small differences in alignment of the tube lenses for each camera. This is done using grid search, similar to how angles_within_channels are estimated. We search for the best angles and scales while iteratively refining the search range. This will yield scales_between_channels, angles_between_channels, and shifts_between_channels.

Estimate for all tiles

We can then use the parameters estimated for the reference tile to register all tiles with:

iss estimate-shifts

This is necessary for two reasons. First, dichroic wobble slightly during and between acquisitions resulting in different shifts between channels. Second, the gain of the microscope stage seems to vary from day to day. Therefore, the microscope does not consistently move to the same position for each tile from round to round, resulting in different shifts across rounds. Therefore, we will re-estimate shifts, both within and across channel, but will not change angles_within_channels, angles_between_channels and scales_between_channels.

Note

angles_within_channels and angles_between_channels might actually vary due to the dichroic wobble but in practice registration works well using values from the reference tile.

The output is saved in the reg subfolder as f"tforms_{prefix}_{roi}_{tilex}_{tiley}.npz"

Correct shift with ransac

The single tile estimation tends to fail sporadically if there is not enough signal. This can be corrected given that the main change of shifts from tile to tile is a linear function of X and Y (probably due to change in gain of the stage). We do that with RANSAC robust regression in:

iss correct-shifts.

Once again, this does not re-estimate angles and scale changes, just shifts. The output is saved in the reg folder as f"tforms_corrected_{prefix}_{roi}_{tilex}_{tiley}.npz"

However, this correction is not ideal for tiles that were already properly registered and can introduce bigger shifts. Therefore, we only apply this correction to tiles that have a shift above a certain threshold. This threshold is currently set in ops['ransac_residual_threshold']

The final transformation is then saved in the reg folder as f"tforms_best_{prefix}_{roi}_{tilex}_{tiley}.npz"

Detailed explanation part 2: Stitching tiles

The information computed above allows us to load all tiles in their “acquisition” coordinates (same for all tiles of one prefix, but different across prefixes).

Find tile shifts

We estimate how much overlap there is between tiles (and therefore how much we need to shift them to merge) by phase correlation. This also takes into account that the cameras may not be perfectly aligned with the stage, therefore there might be (and usually will be) a shift in both X and Y between both rows and columns.

This is done by calling:

shift_right, shift_down, tile_shape = iss.pipeline.register_adjacent_tiles(
    data_path, ref_coors=ops['ref_tile'], prefix='genes_round_1_1')

The output is currently not saved.

Merge coordinates

With these tile shift we can find the position of each tile, simply by multiplying the tile number by the shift.

This can be done with:

roi_dims = np.load(processed_path / data_path / f"{prefix}_roi_dims.npy")
ntiles = roi_dims[roi_dims[:, 0] == 1, 1:][0] + 1
tile_origins, tile_centers = iss.pipeline.calculate_tile_positions(
        shift_right, shift_down, tile_shape, ntiles)

The output is currently not saved.

Registering acquisition together

The final reference coordinate is (for now) genes_round. We can register each acquisition independently first. Then we want to merge them. To do that we generate a downsampled stitched image of the reference acquisition and the acquisition we want to register.

This is done for raw images with iss.pipeline.stitch_and_register. It returns the two registered mosaic at full resolution as well as the transformation parameter: shift and angle.

This output is not saved for now.

For spots, the same function is called by iss-reg2ref align-spots