PICASSO 5-color example

This tutorial documents the interactive script user_scripts/unmix_picasso_5color_example.py.

It is the main public example for genuine multi-channel blind unmixing beyond the 3-channel case.

The PICASSO-family blind-unmixing logic used here is motivated by the original PICASSO paper:

Seo, J., Sim, Y., Kim, J. et al. PICASSO allows ultra-multiplexed fluorescence imaging of spatially overlapping proteins without reference spectra measurements. Nature Communications 13, 2475 (2022). https://doi.org/10.1038/s41467-022-30168-z

How to use this tutorial

The script is intended for interactive execution in a cell-based editor.

The recommended workflow is:

1. open user_scripts/unmix_picasso_5color_example.py,

2. inspect the measured channels first,

3. compare the generalized MATLAB-style path with the source-sink path.

The subsections below follow the same order as the script.

What this tutorial is good for

This is the best public example in the repository for understanding how the package behaves when there are many channels and when the cross-talk graph is not obvious in advance.

Imports

from __future__ import annotations

from pathlib import Path

PROJECT_ROOT = Path(__file__).resolve().parents[1]

from spectral_unmixing import report_path_from_output_path, unmix_picasso
from spectral_unmixing.viewer import show_all_channels_in_napari

Define input and output paths

INPUT_PATH = (PROJECT_ROOT / "example_data" / "PICASSO_examples" / "5_color_unmixing_simulation.tif")
INPUT_NAME = INPUT_PATH.stem

OUTPUT_DIR = INPUT_PATH.parent / "unmixed"
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

Inspect the measured channels

show_all_channels_in_napari(INPUT_PATH, layer_prefix="Measured 5-color simulation")

As in the other PICASSO tutorials, it is worth inspecting the raw measured channels first before deciding which blind-unmixing strategy is the better fit for the dataset.

Composite view of the raw five-channel stack (top). Channel 0 is shown in cyan (top center), Channel 1 in magenta (center), Channel 2 in yellow (bottom center), Channel 3 in green (bottom), and Channel 4 in red (top right). The channels clearly exhibit cross-talk.

matlab_n blind unmixing on five channels

OUTPUT_PICASSO_MATLAB_N = OUTPUT_DIR / f"{INPUT_NAME}_picasso_matlab_n.tif"
picasso_matlab_n_output = unmix_picasso(
    input_path=INPUT_PATH,
    output_path=OUTPUT_PICASSO_MATLAB_N,
    channels=[0, 1, 2, 3, 4],
    # method="picasso",  # default
    implementation="matlab_n",  # "matlab_3c" or "matlab_n" or "source_sink_n"
    # alpha_mode="reference_t",
    # alpha_reference_t=0,
    background_percentile=1.0,
    # preprocess_alpha_inputs=True,  # recorded for compatibility
    mi_bins=64,
    alpha_max=1.0,
    max_iter=50,
    tolerance=1e-4,
    max_alpha_voxels=250_000,
    step_size=0.2,
    qn=100,
    pixel_bin_size=16,
    alpha_clip=0.5,
    # negativity_threshold=0.0009,
    # clip_every_n_iterations=50,
    random_state=42,
    clip_negative=True,
    output_dtype="float32",
    verbose=True,
)
print(picasso_matlab_n_output)
print(report_path_from_output_path(picasso_matlab_n_output).read_text(encoding="utf-8"))
show_all_channels_in_napari(
    picasso_matlab_n_output,
    layer_prefix="PICASSO MATLAB-N unmixed 5-color simulation")

This method applies the explicit N-channel generalization of the MATLAB-style PICASSO iteration to all five selected channels.

The settings that matter most are:

• channels: selects which measured channels are included in the blind-unmixing run.

• implementation="matlab_n": chooses the generalized MATLAB-style iteration for the selected channels.

• background_percentile: low-percentile background estimate used during preprocessing before the update sequence is estimated.

• preprocess_alpha_inputs: enables or disables that shared preprocessing step.

• mi_bins: retained in the shared API and JSON report for compatibility, but not used by the MATLAB-like implementation itself.

• alpha_max: likewise retained in the shared API and JSON report, but not used directly by the MATLAB-like implementation.

• max_iter: number of iteration steps. More iterations can unmix more strongly but may also increase instability.

• tolerance: convergence threshold for the iterative update sequence.

• max_alpha_voxels: optional voxel cap for large stacks. Lower values speed up the run; higher values use more image information.

• step_size: strength of each update step. Larger values make the update more aggressive; smaller values make it gentler.

• qN: quantization parameter for the MATLAB-style mutual-information estimate.

• pixel_bin_size: spatial binning factor before the mutual-information calculation. Larger values smooth more strongly; smaller values preserve more detail.

• alpha_clip: clipping bound for pairwise coefficients. Larger values allow stronger pairwise subtraction.

• optional negativity_threshold and clip_every_n_iterations: control negativity handling during the iterative updates.

• random_state: random seed used whenever voxel subsampling is needed.

• clip_negative: clips final negative values in the written output to zero.

• output_dtype: controls the saved data type of the unmixed stack.

• verbose: enables terminal progress output during the run.

• alpha_mode and alpha_reference_t: become relevant for real multi-time-point stacks. reference_t estimates one update sequence from one chosen reference time point; per_t estimates one sequence per time point.

This is usually the best choice when you want a broad, symmetric blind-unmixing strategy across many channels without encoding a very explicit source-sink model yourself.

Results of the MATLAB-N blind-unmixing run on the five-channel stack. The workflow did a decent job of removing cross-talk between the channels, but not for all of them. For instance, channel 1 still exhibits some cross-talk from channel 0, and channel 3 still exhibits some cross-talk from channel 2. The same is true for channel 4, which still exhibits some cross-talk from channel 3. Further tuning of the parameters may improve the results, but it is also possible that a more explicit source-sink model is needed to achieve better separation of the channels.

source_sink_n blind unmixing on five channels

OUTPUT_PICASSO_SOURCE_SINK = OUTPUT_DIR / f"{INPUT_NAME}_picasso_source_sink.tif"

sink_channels = [0, 1, 2, 3, 4]
neutral_channels = []

# Example for a more selective run:
#
# sink_channels = [1, 3]
# neutral_channels = [4]

picasso_source_sink_output = unmix_picasso(
    input_path=INPUT_PATH,
    output_path=OUTPUT_PICASSO_SOURCE_SINK,
    channels=[0, 1, 2, 3, 4],
    # method="picasso",  # default
    implementation="source_sink_n",
    # alpha_mode="reference_t",
    # alpha_reference_t=0,
    sink_channels=sink_channels,
    neutral_channels=neutral_channels,
    # source_sink_matrix=[
    #     [1, -1, 0, 0, 0],
    #     [0, 1, 0, -1, 0],
    #     [0, 0, 1, 0, 0],
    #     [0, -1, 0, 1, 0],
    #     [0, 0, 0, 0, 1],
    # ],
    background_percentile=1.0,
    # preprocess_alpha_inputs=True,  # recorded for compatibility
    mi_bins=64,
    alpha_max=1.0,
    max_iter=50,
    tolerance=1e-4,
    max_alpha_voxels=250_000,
    source_sink_optimize_background=True,
    source_sink_max_background=0.2,
    source_sink_n_restarts=6,
    source_sink_joint_optimization=True,
    random_state=42,
    clip_negative=True,
    output_dtype="float32",
    verbose=True,
)
print(picasso_source_sink_output)
print(report_path_from_output_path(picasso_source_sink_output).read_text(encoding="utf-8"))
show_all_channels_in_napari(picasso_source_sink_output, layer_prefix="PICASSO source-sink-N unmixed 5-color simulation")

This method uses the more explicit source-sink formulation.

In the current implementation, all sources contributing to one sink are optimized jointly by default, and the workflow can additionally estimate a small background offset for each modeled source-sink relation. This brings the method closer to the source-sink formulation used by the napari PICASSO plugin, while still relying on histogram-based mutual information rather than the plugin’s neural MINE estimator.

The most important settings are:

• channels: as above, selects which measured channels participate in the explicit source-sink run.

• implementation="source_sink_n": chooses the source-sink PICASSO-family workflow rather than the MATLAB-like implementations.

• sink_channels: defines which channels should actually be corrected as sinks.

• neutral_channels: defines which channels should stay untouched and be excluded from active correction roles.

• optional source_sink_matrix: gives explicit manual control over the allowed source-to-sink graph.

• background_percentile: same role as above, but now used in source-sink coefficient estimation.

• preprocess_alpha_inputs: same shared preprocessing switch as above.

• alpha_max: now actively used as the upper bound for source-to-sink coefficients.

• mi_bins: now actively used as the histogram resolution for the mutual-information objective.

• max_iter: controls the maximum number of optimizer iterations for each sink.

• tolerance: stopping tolerance for the numerical optimizer.

• max_alpha_voxels: same optional voxel cap as above.

• source_sink_optimize_background: if enabled, estimate one small background offset beta per modeled source-sink relation before subtracting the source contribution.

• source_sink_max_background: upper bound for these optimized background offsets on normalized intensities.

• source_sink_joint_optimization: if enabled, optimize all sources contributing to the same sink together instead of fitting them greedily one after another.

• source_sink_n_restarts: number of multi-start optimizer initializations used for each sink.

• random_state: same random seed used for optional voxel subsampling.

• clip_negative: same final clipping behavior as above.

• output_dtype: same output dtype control as above.

• verbose: same terminal verbosity control as above.

• alpha_mode and alpha_reference_t: same time-axis controls as above for real multi-time-point stacks.

For real five-channel data this mode is often attractive because it lets you move from a broad first-pass model to a more selective graph once you have a better idea which channels are plausible sinks and which should remain neutral. The example script therefore starts with a broad all-to-all non-neutral model, but it also shows how to restrict the sink list or declare selected channels neutral once prior knowledge about the cross-talk graph becomes available.

THis time, channel 0 and 1 are clearly separated. The same is true for channel 4. However, cross-talk remains between channel 1, 2 and 3. Further tuning of the parameters may improve the results, but it is also possible that a more explicit source-sink model is needed to achieve better separation of the channels.