Skip to content

Extract

Viewer

view_filter

Function to view filtering of raw data in napari. There will be 2 scrollbars in 3D. One to change between raw/difference_of_hanning/difference_of_hanning+smoothed and one to change z-plane.

There are also sliders to change the parameters for the filtering/smoothing. When the sliders are changed, the time taken for the new filtering/smoothing will be printed to the console.

Note

When r_smooth is set to [1, 1, 1], no smoothing will be performed. When this is the case, changing the filtering radius using the slider will be quicker because it will only do filtering and not any smoothing.

If r == anchor_round and c == dapi_channel, the filtering will be tophat filtering and no smoothing will be allowed. Otherwise, the filtering will be convolution with a difference of hanning kernel.

The current difference of hanning kernel can be viewed by pressing the 'h' key.

Requires access to nb.file_names.input_dir

Parameters:

Name Type Description Default
nb Optional[Notebook]

Notebook for experiment. If no Notebook exists, pass config_file instead.

 None
t int

npy (as opposed to nd2 fov) tile index to view. For an experiment where the tiles are arranged in a 4 x 3 (ny x nx) grid, tile indices are indicated as below:

| 2 | 1 | 0 |

| 5 | 4 | 3 |

| 8 | 7 | 6 |

| 11 | 10 | 9 |

 0
r int

round to view

 0
c int

Channel to view.

 0
use_z Optional[Union[int, List[int]]]

Which z-planes to load in from raw data. If None, will use load all z-planes (except from first one if config['basic_info']['ignore_first_z_plane'] == True).

 None
config_file Optional[str]

path to config file for experiment.

 None
Source code in coppafish/plot/extract/viewer.py 
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
def __init__(self, nb: Optional[Notebook] = None, t: int = 0,
             r: int = 0,
             c: int = 0,
             use_z: Optional[Union[int, List[int]]] = None, config_file: Optional[str] = None):
    """
    Function to view filtering of raw data in napari.
    There will be 2 scrollbars in 3D.
    One to change between *raw/difference_of_hanning/difference_of_hanning+smoothed* and one to change z-plane.

    There are also sliders to change the parameters for the filtering/smoothing.
    When the sliders are changed, the time taken for the new filtering/smoothing
    will be printed to the console.

    !!! note
        When `r_smooth` is set to `[1, 1, 1]`, no smoothing will be performed.
        When this is the case, changing the filtering radius using the slider will
        be quicker because it will only do filtering and not any smoothing.

    If `r == anchor_round` and `c == dapi_channel`, the filtering will be tophat filtering and no smoothing
    will be allowed. Otherwise, the filtering will be convolution with a difference of hanning kernel.

    The current difference of hanning kernel can be viewed by pressing the 'h' key.

    !!! warning "Requires access to `nb.file_names.input_dir`"

    Args:
        nb: *Notebook* for experiment. If no *Notebook* exists, pass `config_file` instead.
        t: npy (as opposed to nd2 fov) tile index to view.
            For an experiment where the tiles are arranged in a 4 x 3 (ny x nx) grid, tile indices are indicated as
            below:

            | 2  | 1  | 0  |

            | 5  | 4  | 3  |

            | 8  | 7  | 6  |

            | 11 | 10 | 9  |
        r: round to view
        c: Channel to view.
        use_z: Which z-planes to load in from raw data. If `None`, will use load all z-planes (except from first
            one if `config['basic_info']['ignore_first_z_plane'] == True`).
        config_file: path to config file for experiment.
    """
    if nb is None:
        nb = Notebook(config_file=config_file)
    add_basic_info_no_save(nb)  # deal with case where there is no notebook yet
    if use_z is None:
        use_z = nb.basic_info.use_z
    t, r, c, use_z = number_to_list([t, r, c, use_z])
    self.image_raw = get_raw_images(nb, t, r, c, use_z)[0, 0, 0]

    self.is_3d = nb.basic_info.is_3d
    if not self.is_3d:
        self.image_raw = extract.focus_stack(self.image_raw)
        self.image_plot = np.zeros((3, nb.basic_info.tile_sz, nb.basic_info.tile_sz), dtype=np.int32)
        self.image_plot[0] = self.image_raw
    else:
        self.image_plot = np.zeros((3, len(use_z), nb.basic_info.tile_sz, nb.basic_info.tile_sz), dtype=np.int32)
        self.image_plot[0] = np.moveaxis(self.image_raw, 2, 0)  # put z axis first for plotting
    self.raw_max = self.image_raw.max()

    # find faulty columns and update self.image_raw so don't have to run strip_hack each time we filter
    self.image_raw, self.bad_columns = extract.strip_hack(self.image_raw)

    # Get default filter info
    # label for each image in image_plot
    self.ax0_labels = ['Raw', 'Difference of Hanning', 'Difference of Hanning and Smoothed']
    if not self.is_3d:
        self.ax0_labels[0] = 'Raw (Focus Stacked)'
    config = nb.get_config()['extract']
    if r[0] == nb.basic_info.anchor_round and c[0] == nb.basic_info.dapi_channel:
        self.dapi = True
        if config['r_dapi'] is None:
            if config['r_dapi_auto_microns'] is not None:
                config['r_dapi'] = extract.get_pixel_length(config['r_dapi_auto_microns'],
                                                            nb.basic_info.pixel_size_xy)
            else:
                config['r_dapi'] = 48  # good starting value
        self.r_filter = config['r_dapi']
        self.image_plot = self.image_plot[:2]  # no smoothing if dapi
        self.ax0_labels = self.ax0_labels[:2]
        self.update_filter_image()
        r_filter_lims = [10, 70]
    else:
        self.dapi = False
        self.r_filter = config['r1']
        r2 = config['r2']
        if self.r_filter is None:
            self.r_filter = extract.get_pixel_length(config['r1_auto_microns'], nb.basic_info.pixel_size_xy)
        if r2 is None:
            r2 = self.r_filter * 2
        r_filter_lims = [2, 10]
        self.update_filter_image(r2)
        self.r_filter2 = r2

        # Get default smoothing info
        if config['r_smooth'] is None:
            # start with no smoothing. Quicker to change filter params as no need to update smoothing too.
            config['r_smooth'] = [1, 1, 1]
        if not nb.basic_info.is_3d:
            config['r_smooth'] = config['r_smooth'][:2]
        self.r_smooth = config['r_smooth']
        self.update_smooth_image()

    self.viewer = napari.Viewer()
    self.viewer.add_image(self.image_plot, name=f"Tile {t[0]}, Round {r[0]}, Channel{c[0]}")
    # Set min image contrast to 0 for better comparison between images
    self.viewer.layers[0].contrast_limits = [0, 0.9 * self.viewer.layers[0].contrast_limits[1]]
    self.ax0_ind = 0
    self.viewer.dims.set_point(0, self.ax0_ind)  # set filter type to raw initially
    self.viewer.dims.events.current_step.connect(self.filter_type_status)

    self.filter_slider = QSlider(Qt.Orientation.Horizontal)
    self.filter_slider.setRange(r_filter_lims[0], r_filter_lims[1])
    self.filter_slider.setValue(self.r_filter)
    # When dragging, status will show r_filter value
    self.filter_slider.valueChanged.connect(lambda x: self.show_filter_radius(x))
    # On release of slider, filtered / smoothed images updated
    self.filter_slider.sliderReleased.connect(self.filter_slider_func)
    if self.dapi:
        filter_slider_name = 'Tophat kernel radius'
    else:
        filter_slider_name = 'Difference of Hanning Radius'
    self.viewer.window.add_dock_widget(self.filter_slider, area="left", name=filter_slider_name)

    if not self.dapi:
        self.smooth_yx_slider = QSlider(Qt.Orientation.Horizontal)
        self.smooth_yx_slider.setRange(1, 5)  # gets very slow with large values
        self.smooth_yx_slider.setValue(self.r_smooth[0])
        # When dragging, status will show r_smooth value
        self.smooth_yx_slider.valueChanged.connect(lambda x: self.show_smooth_radius_yx(x))
        # On release of slider, smoothed image updated
        self.smooth_yx_slider.sliderReleased.connect(self.smooth_slider_func)
        smooth_title = "Smooth Radius"
        if self.is_3d:
            smooth_title = smooth_title + " YX"
        self.viewer.window.add_dock_widget(self.smooth_yx_slider, area="left", name=smooth_title)

    if self.is_3d and not self.dapi:
        self.smooth_z_slider = QSlider(Qt.Orientation.Horizontal)
        self.smooth_z_slider.setRange(1, 5)  # gets very slow with large values
        self.smooth_z_slider.setValue(self.r_smooth[2])
        # When dragging, status will show r_smooth value
        self.smooth_z_slider.valueChanged.connect(lambda x: self.show_smooth_radius_z(x))
        # On release of slider, smoothed image updated
        self.smooth_z_slider.sliderReleased.connect(self.smooth_slider_func)
        self.viewer.window.add_dock_widget(self.smooth_z_slider, area="left", name="Smooth Radius Z")

    if self.is_3d:
        self.viewer.dims.axis_labels = ['Filter Method', 'z', 'y', 'x']
    else:
        self.viewer.dims.axis_labels = ['Filter Method', 'y', 'x']

    self.key_call_functions()
    napari.run()

Diagnostics

thresh_box_plots

Function to plot distribution of auto_threshold values amongst tiles for each round and channel.

Parameters:

Name Type Description Default
nb Notebook

Notebook containing the extract NotebookPage.

 required
Source code in coppafish/plot/extract/diagnostics.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
def thresh_box_plots(nb: Notebook):
    """
    Function to plot distribution of auto_threshold values amongst tiles for each round and channel.

    Args:
        nb: Notebook containing the extract NotebookPage.
    """
    box_data = [nb.extract.auto_thresh[:, r, c] for c in nb.basic_info.use_channels for r in nb.basic_info.use_rounds]
    if nb.basic_info.use_anchor:
        box_data = box_data + [nb.extract.auto_thresh[:, nb.basic_info.anchor_round, nb.basic_info.anchor_channel]]
    fig, ax1 = plt.subplots(figsize=(10, 6))
    fig.subplots_adjust(left=0.075, right=0.95, bottom=0.15)
    n_use_channels = len(nb.basic_info.use_channels)
    # different boxplot color for each channel (+ different color for anchor channel)
    # Must be distinct from black and white
    channel_colors = distinctipy.get_colors(n_use_channels + int(nb.basic_info.use_anchor), [(0, 0, 0), (1, 1, 1)])
    bp = ax1.boxplot(box_data, notch=0, sym='+', patch_artist=True)

    c = -1
    tick_labels = np.tile(nb.basic_info.use_rounds, n_use_channels).tolist()
    leg_labels = nb.basic_info.use_channels
    if nb.basic_info.use_anchor:
        tick_labels = tick_labels + ['Anchor']
        leg_labels = leg_labels + [f'Anchor ({nb.basic_info.anchor_channel})']
    ax1.set_xticklabels(tick_labels)
    if nb.basic_info.use_anchor:
        ticks = ax1.get_xticklabels()
        ticks[-1].set_rotation(90)

    leg_markers = []
    for i in range(len(box_data)):
        if i % n_use_channels == 0:
            c += 1
            leg_markers = leg_markers + [bp['boxes'][i]]
        bp['boxes'][i].set_facecolor(channel_colors[c])
    ax1.legend(leg_markers, leg_labels, title='Channel')
    ax1.set_xlabel('Round')
    ax1.set_ylabel('Auto Threshold')
    ax1.set_title('Boxplots showing distribution of Auto Threshold amongst tiles for each round and channel')
    plt.show()

histogram_plots

Plots histograms showing distribution of intensity values combined from all tiles for each round and channel. There is also a Norm button which equalises color channels so all color channels should have most intensity values between -1 and 1.

In the normalised histograms, a good channel will have a sharp peak near 0 accounting for non-spot pixels and a long tail from around 0.1 to just beyond 1 accounting for spot pixels.

Parameters:

Name Type Description Default
nb Notebook

Notebook containing the extract NotebookPage.

 required
Source code in coppafish/plot/extract/diagnostics.py 
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
def __init__(self, nb: Notebook):
    """
    Plots histograms showing distribution of intensity values combined from all tiles for each round and channel.
    There is also a Norm button which equalises color channels so all color channels should have most intensity
    values between -1 and 1.

    In the normalised histograms, a good channel will have a sharp peak near 0 accounting for non-spot pixels
    and a long tail from around 0.1 to just beyond 1 accounting for spot pixels.

    Args:
        nb: Notebook containing the extract NotebookPage.
    """
    self.use_rounds = nb.basic_info.use_rounds
    self.use_channels = nb.basic_info.use_channels
    self.hist_values = nb.extract.hist_values
    self.hist_values_norm = np.arange(-10, 10.004, 0.005)
    n_use_rounds = len(self.use_rounds)
    n_use_channels = len(self.use_channels)
    self.fig, self.ax1 = plt.subplots(n_use_channels, n_use_rounds, figsize=(10, 6), sharey=True, sharex=True)
    self.ax1 = self.ax1.flatten()

    # Compute color_norm_factor as it will be computed at call_spots step of pipeline
    config = nb.get_config()['call_spots']
    rc_ind = np.ix_(self.use_rounds, self.use_channels)
    hist_counts_use = np.moveaxis(np.moveaxis(nb.extract.hist_counts, 0, -1)[rc_ind], -1, 0)
    self.color_norm_factor = color_normalisation(self.hist_values, hist_counts_use,
                                                 config['color_norm_intensities'],
                                                 config['color_norm_probs'], config['bleed_matrix_method'])
    self.color_norm_factor = self.color_norm_factor.T.flatten()

    i = 0
    min_value = 3  # clip hist counts to this so don't get log(0) error
    self.plot_lines = []
    self.hist_counts_norm = []
    self.hist_counts = []
    for c in self.use_channels:
        for r in self.use_rounds:
            # Clip histogram to stop log(0) error.
            self.hist_counts = self.hist_counts + [np.clip(nb.extract.hist_counts[:, r, c], min_value, np.inf)]
            self.hist_counts_norm = self.hist_counts_norm + \
                                    [resample_histogram(self.hist_values / self.color_norm_factor[i],
                                                        self.hist_counts[i], self.hist_values_norm)]

            # Normalise histograms to give probabilities
            self.hist_counts[i] = self.hist_counts[i] / np.sum(nb.extract.hist_counts[:, r, c])
            self.hist_counts_norm[i] = self.hist_counts_norm[i] / np.sum(nb.extract.hist_counts[:, r, c])
            self.plot_lines = self.plot_lines + self.ax1[i].plot(self.hist_values, self.hist_counts[i])
            if r == nb.basic_info.use_rounds[0]:
                self.ax1[i].set_ylabel(c)
            if c == nb.basic_info.use_channels[-1]:
                self.ax1[i].set_xlabel(r)
            i += 1
    self.ax1[0].set_yscale('log')
    self.fig.supylabel('Channel')
    self.fig.supxlabel('Round')
    plt.suptitle('Histograms showing distribution of intensity values combined from '
                 'all tiles for each round and channel')

    self.norm = False
    self.xlims_norm = [-1, 1]
    self.xlims = [-300, 300]
    self.ax1[0].set_xlim(self.xlims[0], self.xlims[1])
    self.norm_button_ax = self.fig.add_axes([0.85, 0.02, 0.1, 0.05])
    self.norm_button = Button(self.norm_button_ax, 'Norm', hovercolor='0.275')
    self.norm_button.on_clicked(self.change_norm)

    plt.show()