Skip to content

OMP

Coefficients

view_omp

Diagnostic to show omp coefficients of all genes in neighbourhood of spot. Only genes for which a significant number of pixels are non-zero will be plotted.

Requires access to nb.file_names.tile_dir

Parameters:

Name Type Description Default
nb Notebook

Notebook containing experiment details. Must have run at least as far as call_reference_spots.

 required
spot_no int

Spot of interest to be plotted.

 required
method str

'anchor' or 'omp'. Which method of gene assignment used i.e. spot_no belongs to ref_spots or omp page of Notebook.

 'omp'
im_size int

Radius of image to be plotted for each gene.

 8
Source code in coppafish/plot/omp/coefs.py 
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
172
173
174
175
176
177
178
179
180
181
def __init__(self, nb: Notebook, spot_no: int, method: str = 'omp', im_size: int = 8):
    """
    Diagnostic to show omp coefficients of all genes in neighbourhood of spot.
    Only genes for which a significant number of pixels are non-zero will be plotted.

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

    Args:
        nb: Notebook containing experiment details. Must have run at least as far as `call_reference_spots`.
        spot_no: Spot of interest to be plotted.
        method: `'anchor'` or `'omp'`.
            Which method of gene assignment used i.e. `spot_no` belongs to `ref_spots` or `omp` page of Notebook.
        im_size: Radius of image to be plotted for each gene.
    """
    coef_images, min_global_yxz, max_global_yxz = get_coef_images(nb, spot_no, method, [im_size, im_size, 0])

    if method.lower() == 'omp':
        page_name = 'omp'
        config = nb.get_config()['thresholds']
        spot_score = omp_spot_score(nb.omp, config['score_omp_multiplier'], spot_no)
    else:
        page_name = 'ref_spots'
        spot_score = nb.ref_spots.score[spot_no]
    gene_no = nb.__getattribute__(page_name).gene_no[spot_no]
    t = nb.__getattribute__(page_name).tile[spot_no]
    spot_yxz = nb.__getattribute__(page_name).local_yxz[spot_no]
    gene_name = nb.call_spots.gene_names[gene_no]
    all_gene_names = list(nb.call_spots.gene_names) + [f'BG{i}' for i in range(nb.basic_info.n_channels)]
    spot_yxz_global = spot_yxz + nb.stitch.tile_origin[t]
    n_genes = nb.call_spots.bled_codes_ge.shape[0]

    n_nonzero_pixels_thresh = np.min([im_size, 5])  # If 5 pixels non-zero, plot that gene
    plot_genes = np.where(np.sum(coef_images != 0, axis=(1, 2, 3)) > n_nonzero_pixels_thresh)[0]
    coef_images = coef_images[plot_genes, :, :, 0]
    n_plot = len(plot_genes)
    # at most n_max_rows rows
    if n_plot <= 16:
        n_max_rows = 4
    else:
        n_max_rows = int(np.ceil(np.sqrt(n_plot)))
    n_cols = int(np.ceil(n_plot / n_max_rows))
    subplot_row_columns = [int(np.ceil(n_plot / n_cols)), n_cols]
    fig_size = np.clip([n_cols+5, subplot_row_columns[0]+4], 3, 12)
    subplot_adjust = [0.05, 0.775, 0.05, 0.91]
    super().__init__(coef_images, None, subplot_row_columns, subplot_adjust=subplot_adjust, fig_size=fig_size,
                     cbar_pos=[0.9, 0.05, 0.03, 0.86], slider_pos=[0.85, 0.05, 0.01, 0.86])
    # set x, y coordinates to be those of the global coordinate system
    plot_extent = [min_global_yxz[1]-0.5, max_global_yxz[1]+0.5,
                   min_global_yxz[0]-0.5, max_global_yxz[0]+0.5]
    for i in range(self.n_images):
        # Add cross-hair
        self.ax[i].axes.plot([spot_yxz_global[1], spot_yxz_global[1]], [plot_extent[2], plot_extent[3]],
                             'k', linestyle=":", lw=1)
        self.ax[i].axes.plot([plot_extent[0], plot_extent[1]], [spot_yxz_global[0], spot_yxz_global[0]],
                             'k', linestyle=":", lw=1)
        self.im[i].set_extent(plot_extent)
        self.ax[i].tick_params(labelbottom=False, labelleft=False)
        # Add title
        title_text = f'{plot_genes[i]}: {all_gene_names[plot_genes[i]]}'
        if plot_genes[i] >= n_genes:
            text_color = (0.7, 0.7, 0.7)  # If background, make grey
            title_text = all_gene_names[plot_genes[i]]
        elif plot_genes[i] == gene_no:
            text_color = 'g'
        else:
            text_color = 'w'  # TODO: maybe make color same as used in plot for each gene
        self.ax[i].set_title(title_text, color=text_color)
    plt.subplots_adjust(hspace=0.32)
    plt.suptitle(f'OMP gene coefficients for spot {spot_no} (match'
                 f' {str(np.around(spot_score, 2))} to {gene_name})',
                 x=(subplot_adjust[0] + subplot_adjust[1]) / 2, size=13)
    self.change_norm()
    plt.show()

view_omp_fit

Diagnostic to run omp on a single pixel and see which genes fitted at which iteration. Right-clicking on a particular bled code will cause coppafish.plot.call_spots.dot_product.view_score to run, indicating how the dot product calculation for that iteration was performed.

Left-clicking on background image will cause coppafish.plot.call_spots.background.view_background to run, indicating how the dot product calculation for performed.

Parameters:

Name Type Description Default
nb Notebook

Notebook containing experiment details. Must have run at least as far as call_reference_spots.

 required
spot_no int

Spot of interest to be plotted.

 required
method str

'anchor' or 'omp'. Which method of gene assignment used i.e. spot_no belongs to ref_spots or omp page of Notebook.

 'omp'
dp_thresh Optional[float]

If None, will use value in omp section of config file.

 None
max_genes Optional[int]

If None, will use value in omp section of config file.

 None
Source code in coppafish/plot/omp/coefs.py 
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
def __init__(self, nb: Notebook, spot_no: int, method: str = 'omp', dp_thresh: Optional[float] = None,
             max_genes: Optional[int] = None):
    """
    Diagnostic to run omp on a single pixel and see which genes fitted at which iteration.
    Right-clicking on a particular bled code will cause *coppafish.plot.call_spots.dot_product.view_score*
    to run, indicating how the dot product calculation for that iteration was performed.

    Left-clicking on background image will cause coppafish.plot.call_spots.background.view_background to run,
    indicating how the dot product calculation for performed.

    Args:
        nb: Notebook containing experiment details. Must have run at least as far as `call_reference_spots`.
        spot_no: Spot of interest to be plotted.
        method: `'anchor'` or `'omp'`.
            Which method of gene assignment used i.e. `spot_no` belongs to `ref_spots` or `omp` page of Notebook.
        dp_thresh: If None, will use value in omp section of config file.
        max_genes: If None, will use value in omp section of config file.
    """
    track_info, bled_codes, dp_thresh = get_track_info(nb, spot_no, method, dp_thresh, max_genes)
    # Add info so can call view_dot_product
    self.nb = nb
    self.track_info = track_info
    self.bled_codes = bled_codes
    self.dp_thresh = dp_thresh
    self.spot_no = spot_no
    self.fitting_method = method
    n_genes, n_use_rounds, n_use_channels = bled_codes.shape

    n_residual_images = track_info['residual'].shape[0]
    residual_images = [track_info['residual'][i].transpose() for i in range(n_residual_images)]
    background_image = np.zeros((n_use_rounds, n_use_channels))
    for c in range(n_use_channels):
        background_image += track_info['background_codes'][c] * track_info['background_coefs'][c]
    background_image = background_image.transpose()

    # allow for possibly adding background vector
    # TODO: Think may get index error if best gene ever was background_vector.
    bled_codes = np.append(bled_codes, track_info['background_codes'], axis=0)
    all_gene_names = list(nb.call_spots.gene_names) + [f'BG{i}' for i in nb.basic_info.use_channels]
    gene_images = [bled_codes[track_info['gene_added'][i]].transpose() *
                   track_info['coef'][i][track_info['gene_added'][i]] for i in range(2, n_residual_images)]
    all_images = residual_images + [background_image] + gene_images

    # Plot all images
    subplot_adjust = [0.06, 0.82, 0.075, 0.9]
    super().__init__(all_images, None, [2, n_residual_images], subplot_adjust=subplot_adjust, fig_size=(15, 7))

    # label axis
    self.ax[0].set_yticks(ticks=np.arange(self.im_data[0].shape[0]), labels=nb.basic_info.use_channels)
    self.ax[0].set_xticks(ticks=np.arange(self.im_data[0].shape[1]), labels=nb.basic_info.use_rounds)
    self.fig.supxlabel('Round', size=12, x=(subplot_adjust[0] + subplot_adjust[1]) / 2)
    self.fig.supylabel('Color Channel', size=12)
    plt.suptitle(f'Residual at each iteration of OMP for Spot {spot_no}. '+r'$\Delta_{thresh}$ = '+f'{dp_thresh}',
                 x=(subplot_adjust[0] + subplot_adjust[1]) / 2)

    # Add titles for each subplot
    titles = ['Initial', 'Post Background']
    for g in track_info['gene_added'][2:]:
        titles = titles + ['Post ' + all_gene_names[g]]
    for i in range(n_residual_images):
        titles[i] = titles[i] + '\nRes = {:.2f}'.format(np.linalg.norm(residual_images[i]))
    titles = titles + ['Background']
    for i in range(2, n_residual_images):
        g = track_info['gene_added'][i]
        titles = titles + [f'{g}: {all_gene_names[g]}']
        titles[-1] = titles[-1] + '\n'+r'$\Delta_{s'+f'{i-2}'+'}$ = '+'{:.2f}'.format(track_info['dot_product'][i])

    # Make title red if dot product fell below dp_thresh or if best gene background
    is_fail_thresh = False
    for i in range(self.n_images):
        if np.isin(i, np.arange(2, n_residual_images)):
            # failed if bad dot product, gene added is background or gene added has previously been added
            is_fail_thresh = np.abs(track_info['dot_product'][i]) < dp_thresh or \
                             track_info['gene_added'][i] >= n_genes or \
                             np.isin(track_info['gene_added'][i], track_info['gene_added'][2:i])
            if is_fail_thresh:
                text_color = 'r'
            else:
                text_color = 'w'
        elif i == self.n_images - 1 and is_fail_thresh:
            text_color = 'r'
        else:
            text_color = 'w'
        self.ax[i].set_title(titles[i], size=8, color=text_color)

    # Add rectangles where added gene is intense
    for i in range(len(gene_images)):
        gene_coef = track_info['coef'][i+2][track_info['gene_added'][i+2]]
        intense_gene_cr = np.where(np.abs(gene_images[i] / gene_coef) > self.intense_gene_thresh)
        for j in range(len(intense_gene_cr[0])):
            for k in [i+1, i+1+n_residual_images]:
                # can't add rectangle to multiple axes hence second for loop
                rectangle = plt.Rectangle((intense_gene_cr[1][j]-0.5, intense_gene_cr[0][j]-0.5), 1, 1,
                                          fill=False, ec="g", linestyle=':', lw=2)
                self.ax[k].add_patch(rectangle)

    self.change_norm()
    self.fig.canvas.mpl_connect('button_press_event', self.show_calc)
    self.track_info = track_info
    plt.show()

get_coef_images

Gets image of \(yxz\) dimension (2*im_size[0]+1) x (2*im_size[1]+1) x (2*im_size[2]+1) of the coefficients fitted by omp for each gene.

Parameters:

Name Type Description Default
nb Notebook

Notebook containing experiment details. Must have run at least as far as call_reference_spots.

 required
spot_no int

Spot of interest to get gene coefficient images for.

 required
method str

'anchor' or 'omp'. Which method of gene assignment used i.e. spot_no belongs to ref_spots or omp page of Notebook.

 required
im_size List[int]

\(yxz\) radius of image to get for each gene.

 required

Returns:

Type Description
np.ndarray

coef_images - float16 [n_genes x (2*im_size[0]+1) x (2*im_size[1]+1) x (2*im_size[2]+1)]. Image for each gene, axis order is \(gyxz\). coef_images[g, 0, 0, 0] refers to coefficient of gene g at global_yxz = min_global_yxz. coef_images[g, -1, -1, -1] refers to coefficient of gene g at global_yxz = max_global_yxz.
List[float]

min_global_yxz - float [3]. Min \(yxz\) coordinates of image in global coordinates.
List[float]

max_global_yxz - float [3]. Max \(yxz\) coordinates of image in global coordinates.
Source code in coppafish/plot/omp/coefs.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
def get_coef_images(nb: Notebook, spot_no: int, method: str, im_size: List[int]) -> Tuple[np.ndarray, List[float],
                                                                                          List[float]]:
    """
    Gets image of $yxz$ dimension `(2*im_size[0]+1) x (2*im_size[1]+1) x (2*im_size[2]+1)` of the coefficients
    fitted by omp for each gene.

    Args:
        nb: *Notebook* containing experiment details. Must have run at least as far as `call_reference_spots`.
        spot_no: Spot of interest to get gene coefficient images for.
        method: `'anchor'` or `'omp'`.
            Which method of gene assignment used i.e. `spot_no` belongs to `ref_spots` or `omp` page of Notebook.
        im_size: $yxz$ radius of image to get for each gene.

    Returns:
        `coef_images` - `float16 [n_genes x (2*im_size[0]+1) x (2*im_size[1]+1) x (2*im_size[2]+1)]`.
            Image for each gene, axis order is $gyxz$.
            `coef_images[g, 0, 0, 0]` refers to coefficient of gene g at `global_yxz = min_global_yxz`.
            `coef_images[g, -1, -1, -1]` refers to coefficient of gene g at `global_yxz = max_global_yxz`.
        `min_global_yxz` - `float [3]`. Min $yxz$ coordinates of image in global coordinates.
        `max_global_yxz` - `float [3]`. Max $yxz$ coordinates of image in global coordinates.
    """
    color_norm = nb.call_spots.color_norm_factor[np.ix_(nb.basic_info.use_rounds,
                                                        nb.basic_info.use_channels)]

    if method.lower() == 'omp':
        page_name = 'omp'
    else:
        page_name = 'ref_spots'
    t = nb.__getattribute__(page_name).tile[spot_no]
    spot_yxz = nb.__getattribute__(page_name).local_yxz[spot_no]

    # Subtlety here, may have y-axis flipped, but I think it is correct:
    # note im_yxz[1] refers to point at max_y, min_x+1, z. So when reshape and set plot_extent, should be correct.
    # I.e. im = np.zeros(49); im[1] = 1; im = im.reshape(7,7); plt.imshow(im, extent=[-0.5, 6.5, -0.5, 6.5])
    # will show the value 1 at max_y, min_x+1.
    im_yxz = np.array(np.meshgrid(np.arange(spot_yxz[0] - im_size[0], spot_yxz[0] + im_size[0] + 1)[::-1],
                                  np.arange(spot_yxz[1] - im_size[1], spot_yxz[1] + im_size[1] + 1),
                                  spot_yxz[2]),
                      dtype=np.int16).T.reshape(-1, 3)
    z = np.arange(-im_size[2], im_size[2]+1)
    im_yxz = np.vstack([im_yxz + [0, 0, val] for val in z])
    im_diameter_yx = [2 * im_size[0] + 1, 2 * im_size[1] + 1]
    spot_colors = get_spot_colors(im_yxz, t, nb.register.transform, nb.file_names, nb.basic_info) / color_norm

    # Only look at pixels with high enough intensity - same as in full pipeline
    spot_intensity = get_spot_intensity(np.abs(spot_colors))
    config = nb.get_config()['omp']
    if nb.has_page('omp'):
        initial_intensity_thresh = nb.omp.initial_intensity_thresh
    else:
        initial_intensity_thresh = get_initial_intensity_thresh(config, nb.call_spots)

    keep = spot_intensity > initial_intensity_thresh
    bled_codes = nb.call_spots.bled_codes_ge
    n_genes = bled_codes.shape[0]
    bled_codes = np.asarray(bled_codes[np.ix_(np.arange(n_genes),
                                              nb.basic_info.use_rounds, nb.basic_info.use_channels)])
    n_use_rounds = len(nb.basic_info.use_rounds)
    dp_norm_shift = nb.call_spots.dp_norm_shift * np.sqrt(n_use_rounds)

    dp_thresh = config['dp_thresh']
    if method.lower() == 'omp':
        alpha = config['alpha']
        beta = config['beta']
    else:
        config_call_spots = nb.get_config()['call_spots']
        alpha = config_call_spots['alpha']
        beta = config_call_spots['beta']
    max_genes = config['max_genes']
    weight_coef_fit = config['weight_coef_fit']

    all_coefs = np.zeros((spot_colors.shape[0], n_genes + nb.basic_info.n_channels))
    all_coefs[np.ix_(keep, np.arange(n_genes))], \
    all_coefs[np.ix_(keep, np.array(nb.basic_info.use_channels) + n_genes)] = \
        get_all_coefs(spot_colors[keep], bled_codes, nb.call_spots.background_weight_shift, dp_norm_shift,
                      dp_thresh, alpha, beta, max_genes, weight_coef_fit)

    n_genes = all_coefs.shape[1]
    nz = len(z)
    coef_images = np.zeros((n_genes, len(z), im_diameter_yx[0], im_diameter_yx[1]))
    for g in range(n_genes):
        ind = 0
        for z in range(nz):
            coef_images[g, z] = all_coefs[ind:ind+np.prod(im_diameter_yx), g].reshape(im_diameter_yx[0],
                                                                                      im_diameter_yx[1])
            ind += np.prod(im_diameter_yx)
    coef_images = np.moveaxis(coef_images, 1, -1)  # move z index to end
    min_global_yxz = im_yxz.min(axis=0)+nb.stitch.tile_origin[t]
    max_global_yxz = im_yxz.max(axis=0)+nb.stitch.tile_origin[t]
    return coef_images.astype(np.float16), min_global_yxz, max_global_yxz

Score Shape

view_omp_score

Diagnostic to show how score is computed in the omp method Hatched region in top plot shows pixels which contribute to the final score. Score is actually equal to the absolute sum of the top plots in the hatched regions.

Can also see how score_omp_multiplier affects the final score. The larger this is, the more the positive pixels contribute compared to the negative.

Requires access to nb.file_names.tile_dir

Parameters:

Name Type Description Default
nb Notebook

Notebook containing experiment details. Must have run at least as far as call_reference_spots.

 required
spot_no int

Spot of interest to be plotted.

 required
method str

'anchor' or 'omp'. Which method of gene assignment used i.e. spot_no belongs to ref_spots or omp page of Notebook.

 'omp'
score_multiplier Optional[float]

Initial value of score_omp_multiplier.

 None
check bool

If True, will compare score found to that saved in Notebook and raise error if they differ. Will also check that absolute sum of the top plots in the hatched regions is equal to score calculated from counting the number of pixels with the correct sign.

 False
Source code in coppafish/plot/omp/score_shape.py 
 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
 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
def __init__(self, nb: Notebook, spot_no: int, method: str = 'omp', score_multiplier: Optional[float] = None,
             check: bool = False):
    """
    Diagnostic to show how score is computed in the omp method
    Hatched region in top plot shows pixels which contribute to the final score.
    Score is actually equal to the absolute sum of the top plots in the hatched regions.

    Can also see how `score_omp_multiplier` affects the final score. The larger this is, the more
    the positive pixels contribute compared to the negative.

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

    Args:
        nb: Notebook containing experiment details. Must have run at least as far as `call_reference_spots`.
        spot_no: Spot of interest to be plotted.
        method: `'anchor'` or `'omp'`.
            Which method of gene assignment used i.e. `spot_no` belongs to `ref_spots` or `omp` page of Notebook.
        score_multiplier: Initial value of `score_omp_multiplier`.
        check: If `True`, will compare score found to that saved in *Notebook* and raise error if they differ.
            Will also check that absolute sum of the top plots in the hatched regions is equal to
            score calculated from counting the number of pixels with the correct sign.

    """
    # TODO: The textbox for this plot seems to be much less responsive than in the other diagnostics
    #  for some reason.
    if method.lower() == 'omp':
        page_name = 'omp'
    else:
        page_name = 'ref_spots'
        check = False
    self.check = check
    self.nbp_omp = nb.omp
    self.spot_no = spot_no
    self.gene_no = nb.__getattribute__(page_name).gene_no[spot_no]
    self.gene_names = nb.call_spots.gene_names
    self.expected_sign_image = nb.omp.spot_shape
    self.nz = self.expected_sign_image.shape[2]
    im_radius = ((np.asarray(self.expected_sign_image.shape) - 1) / 2).astype(int)
    self.coef_images, min_global_yxz, max_global_yxz = get_coef_images(nb, spot_no, method, im_radius)
    # Find where sign of coef image matches that of expected sign image
    self.both_positive = np.asarray([self.coef_images[self.gene_no] > 0, self.expected_sign_image > 0]).all(axis=0)
    self.both_negative = np.asarray([self.coef_images[self.gene_no] < 0, self.expected_sign_image < 0]).all(axis=0)
    self.check_neighbours()

    # Start with default multiplier
    config = nb.get_config()['thresholds']
    if score_multiplier is None:
        score_multiplier = config['score_omp_multiplier']
    self.score_multiplier = score_multiplier
    self.score_thresh = config['score_omp']

    # maximum possible value of any one pixel in expected_shape for any score_multiplier
    self.vmax_expected = 1 / np.min([np.sum(self.expected_sign_image > 0), np.sum(self.expected_sign_image < 0)])
    self.vmax_coef = np.abs(self.coef_images[self.gene_no]).max()

    # Initialize plots
    self.plot_extent = [min_global_yxz[1] - 0.5, max_global_yxz[1] + 0.5,
                        min_global_yxz[0] - 0.5, max_global_yxz[0] + 0.5]
    self.fig, self.ax = plt.subplots(2, self.nz, figsize=(14, 5), sharex=True, sharey=True)
    self.subplot_adjust = [0.05, 0.88, 0.09, 0.85]
    self.fig.subplots_adjust(left=self.subplot_adjust[0], right=self.subplot_adjust[1],
                             bottom=self.subplot_adjust[2], top=self.subplot_adjust[3])
    self.ax = self.ax.flatten()
    self.im = [None] * self.nz * 2
    expected_image_plot = self.expected_image()
    self.score = self.get_score(expected_image_plot)
    for i in range(self.nz):
        self.im[i] = self.ax[i].imshow(expected_image_plot[:, :, i],
                                       vmin=-self.vmax_expected, vmax=self.vmax_expected, cmap='bwr',
                                       extent=self.plot_extent)
        self.im[i+self.nz] = self.ax[i + self.nz].imshow(self.coef_images[self.gene_no, :, :, i],
                                                         vmin=-self.vmax_coef, vmax=self.vmax_coef, cmap='bwr',
                                                         extent=self.plot_extent)
        if i == (self.nz - 1) / 2:
            self.ax[i].set_title(f"Expected Coefficient Sign\nZ={int(np.rint(min_global_yxz[2] + i))}")
        else:
            self.ax[i].set_title(f"Z={int(np.rint(min_global_yxz[2] + i))}")
    self.set_coef_plot_title()
    self.add_hatching()

    # Set up colorbars for each plot
    mid_point = (self.subplot_adjust[2]+self.subplot_adjust[3])/2
    gap_size = 0.08
    cbar_ax = self.fig.add_axes([self.subplot_adjust[1]+0.01, mid_point+gap_size/2,
                                 0.005, self.subplot_adjust[3]-mid_point-gap_size/2])  # left, bottom, width, height
    self.fig.colorbar(self.im[0], cax=cbar_ax)
    cbar_ax = self.fig.add_axes([self.subplot_adjust[1]+0.01, self.subplot_adjust[2]+gap_size/5,
                                 0.005, mid_point-self.subplot_adjust[2]-gap_size/2])  # left, bottom, width, height
    self.fig.colorbar(self.im[self.nz], cax=cbar_ax)

    # Add titles
    self.fig.supylabel('Y')
    self.fig.supxlabel('X', size=12, x=(self.subplot_adjust[0] + self.subplot_adjust[1]) / 2)
    plt.suptitle(f"OMP Score Calculation for Spot {spot_no}, Gene {self.gene_no}: {self.gene_names[self.gene_no]}",
                 x=(self.subplot_adjust[0] + self.subplot_adjust[1]) / 2)

    # Add text box to change score multiplier
    text_ax = self.fig.add_axes([self.subplot_adjust[1] + 0.062, self.subplot_adjust[2]+gap_size/5,
                                 0.05, 0.04])
    self.text_box = TextBox(text_ax, 'Score\n'+r'Multiplier, $\rho$', str(np.around(self.score_multiplier, 2)),
                            color='k', hovercolor=[0.2, 0.2, 0.2])
    self.text_box.cursor.set_color('r')
    label = text_ax.get_children()[0]  # label is a child of the TextBox axis
    label.set_position([0.5, 2.75])  # [x,y] - change here to set the position
    # centering the text
    label.set_verticalalignment('top')
    label.set_horizontalalignment('center')
    self.text_box.on_submit(self.update)

    plt.show()

Score Histogram

histogram_score

If method is anchor, this will show the histogram of nb.ref_spots.score with the option to view the histogram of the score computed using various other configurations of background fitting and gene_efficiency. This allows one to see how the these affect the score.

If method is omp, this will show the histogram of omp score, computed with coppafish.call_spots.omp_spot_score. There will also be the option to view the histograms shown for the anchor method. I.e. we compute the dot product score for the omp spots.

Parameters:

Name Type Description Default
nb Notebook

Notebook containing at least call_spots page.

 required
method str

'anchor' or 'omp'. Which method of gene assignment used i.e. spot_no belongs to ref_spots or omp page of Notebook.

 'omp'
score_omp_multiplier Optional[float]

Can specify the value of score_omp_multiplier to use to compute omp score. If None, will use value in config file.

 None
check bool

If True, and method='anchor', will check that scores computed here match those saved to Notebook.

 False
hist_spacing float

Initial width of bin in histogram.

 0.001
show_plot bool

Whether to run plt.show() or not.

 True
Source code in coppafish/plot/omp/score_hist.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
def __init__(self, nb: Notebook, method: str = 'omp', score_omp_multiplier: Optional[float] = None,
             check: bool = False, hist_spacing: float = 0.001, show_plot: bool = True):
    """
    If method is anchor, this will show the histogram of `nb.ref_spots.score` with the option to
    view the histogram of the score computed using various other configurations of `background` fitting
    and `gene_efficiency`. This allows one to see how the these affect the score.

    If `method` is omp, this will show the histogram of omp score, computed with
    `coppafish.call_spots.omp_spot_score`.
    There will also be the option to view the histograms shown for the anchor method.
    I.e. we compute the dot product score for the omp spots.

    Args:
        nb: *Notebook* containing at least `call_spots` page.
        method: `'anchor'` or `'omp'`.
            Which method of gene assignment used i.e. `spot_no` belongs to `ref_spots` or `omp` page of Notebook.
        score_omp_multiplier: Can specify the value of score_omp_multiplier to use to compute omp score.
            If `None`, will use value in config file.
        check: If `True`, and `method='anchor'`, will check that scores computed here match those saved to Notebook.
        hist_spacing: Initial width of bin in histogram.
        show_plot: Whether to run `plt.show()` or not.
    """
    # Add data
    if score_omp_multiplier is None:
        config = nb.get_config()['thresholds']
        score_omp_multiplier = config['score_omp_multiplier']
    self.score_multiplier = score_omp_multiplier
    self.gene_names = nb.call_spots.gene_names
    self.n_genes = self.gene_names.size
    # Use all genes by default
    self.genes_use = np.arange(self.n_genes)

    # For computing score_dp
    spot_colors, spot_colors_pb, background_var = background_fitting(nb, method)
    grc_ind = np.ix_(np.arange(self.n_genes), nb.basic_info.use_rounds, nb.basic_info.use_channels)
    # Bled codes saved to Notebook should already have L2 norm = 1 over used_channels and rounds
    bled_codes = nb.call_spots.bled_codes[grc_ind]
    bled_codes_ge = nb.call_spots.bled_codes_ge[grc_ind]

    # Save score_dp for each permutation of with/without background/gene_efficiency
    self.n_plots = 5
    if method.lower() == 'omp':
        self.n_plots += 1
        self.nbp_omp = nb.omp
        self.gene_no = self.nbp_omp.gene_no
        self.score = np.zeros((self.gene_no.size, self.n_plots), dtype=np.float32)
        self.score[:, -1] = omp_spot_score(self.nbp_omp, self.score_multiplier)
        self.method = 'OMP'
    else:
        self.gene_no = nb.ref_spots.gene_no
        self.score = np.zeros((self.gene_no.size, self.n_plots), dtype=np.float32)
        self.method = 'Anchor'
    self.use = np.isin(self.gene_no, self.genes_use)  # which spots to plot
    # DP score
    self.score[:, 0] = get_dot_product_score(spot_colors_pb, bled_codes_ge, self.gene_no,
                                             nb.call_spots.dp_norm_shift, background_var)[0]
    if method.lower() != 'omp' and check:
        if np.max(np.abs(self.score[:, 0] - nb.ref_spots.score)) > self.check_tol:
            raise ValueError(f"nb.ref_spots.score differs to that computed here\n"
                             f"Set check=False to get past this error")

    # DP score no weight
    self.score[:, 1] = get_dot_product_score(spot_colors_pb, bled_codes_ge, self.gene_no,
                                             nb.call_spots.dp_norm_shift, None)[0]
    # DP score no background
    self.score[:, 2] = get_dot_product_score(spot_colors, bled_codes_ge, self.gene_no,
                                             nb.call_spots.dp_norm_shift, None)[0]
    # DP score no gene efficiency
    self.score[:, 3] = get_dot_product_score(spot_colors_pb, bled_codes, self.gene_no,
                                             nb.call_spots.dp_norm_shift, background_var)[0]
    # DP score no background or gene efficiency
    self.score[:, 4] = get_dot_product_score(spot_colors, bled_codes, self.gene_no,
                                             nb.call_spots.dp_norm_shift, None)[0]

    # Initialise plot
    self.fig, self.ax = plt.subplots(1, 1, figsize=(11, 5))
    self.subplot_adjust = [0.07, 0.85, 0.1, 0.93]
    self.fig.subplots_adjust(left=self.subplot_adjust[0], right=self.subplot_adjust[1],
                             bottom=self.subplot_adjust[2], top=self.subplot_adjust[3])
    self.ax.set_ylabel(r"Number of Spots")
    if method.lower() == 'omp':
        self.ax.set_xlabel(r"Score, $\gamma_s$ or $\Delta_s$")
    else:
        self.ax.set_xlabel(r"Score, $\Delta_s$")
    self.ax.set_title(f"Distribution of Scores for all {self.method} spots")

    # Set lower bound based on dot product score with no GE/no background as likely to be lowest
    self.hist_min = np.percentile(self.score[:, 4], 0.1)
    # Set upper bound based on dot product score with GE and background because this is likely to be highest
    self.hist_max = np.clip(np.percentile(self.score[:, 0], 99.9), 1, 2)
    self.hist_spacing = hist_spacing
    hist_bins = np.arange(self.hist_min, self.hist_max + self.hist_spacing / 2, self.hist_spacing)
    self.plots = [None] * self.n_plots
    default_colors = plt.rcParams['axes.prop_cycle']._left
    for i in range(self.n_plots):
        y, x = np.histogram(self.score[self.use, i], hist_bins)
        x = x[:-1] + self.hist_spacing / 2  # so same length as x
        self.plots[i], = self.ax.plot(x, y, color=default_colors[i]['color'])
        if method.lower() == 'omp' and i < self.n_plots - 1:
            self.plots[i].set_visible(False)
        elif i > 0 and method.lower() != 'omp':
            self.plots[i].set_visible(False)

    self.ax.set_xlim(self.hist_min, self.hist_max)
    self.ax.set_ylim(0, None)

    # Add text box to change score multiplier
    text_box_labels = ['Gene', 'Histogram\nSpacing', 'Score\n' + r'Multiplier, $\rho$']
    text_box_values = ['all', self.hist_spacing, np.around(self.score_multiplier, 2)]
    text_box_funcs = [self.update_genes, self.update_hist_spacing, self.update_score_multiplier]
    if method.lower() != 'omp':
        text_box_labels = text_box_labels[:2]
        text_box_values = text_box_values[:2]
        text_box_funcs = text_box_funcs[:2]
    self.text_boxes = [None] * len(text_box_labels)
    for i in range(len(text_box_labels)):
        text_ax = self.fig.add_axes([self.subplot_adjust[1] + 0.05,
                                     self.subplot_adjust[2] + 0.15 * (len(text_box_labels) - i - 1), 0.05, 0.04])
        self.text_boxes[i] = TextBox(text_ax, text_box_labels[i], text_box_values[i], color='k',
                                     hovercolor=[0.2, 0.2, 0.2])
        self.text_boxes[i].cursor.set_color('r')
        label = text_ax.get_children()[0]  # label is a child of the TextBox axis
        if i == 0:
            label.set_position([0.5, 1.77])  # [x,y] - change here to set the position
        else:
            label.set_position([0.5, 2.75])
            # centering the text
        label.set_verticalalignment('top')
        label.set_horizontalalignment('center')
        self.text_boxes[i].on_submit(text_box_funcs[i])

    # Add buttons to add/remove score_dp histograms
    self.buttons_ax = self.fig.add_axes([self.subplot_adjust[1] + 0.02, self.subplot_adjust[3] - 0.45, 0.15, 0.5])
    plt.axis('off')
    self.button_labels = [r"$\Delta_s$" + "\nDot Product Score",
                          r"$\Delta_s$" + "\nNo Weighting",
                          r"$\Delta_s$" + "\nNo Background",
                          r"$\Delta_s$" + "\nNo Gene Efficiency",
                          r"$\Delta_s$" + "\nNo Background\nNo Gene Efficiency"]
    label_checked = [True, False, False, False, False]
    if method.lower() == 'omp':
        self.button_labels += [r"$\gamma_s$" + "\nOMP Score"]
        label_checked += [True]
        label_checked[0] = False
    self.buttons = CheckButtons(self.buttons_ax, self.button_labels, label_checked)

    for i in range(self.n_plots):
        self.buttons.labels[i].set_fontsize(7)
        self.buttons.labels[i].set_color(default_colors[i]['color'])
        self.buttons.rectangles[i].set_color('w')
    self.buttons.on_clicked(self.choose_plots)
    if show_plot:
        plt.show()

histogram_2d_score

This plots the bivariate histogram to see the correlation between the omp spot score, \(\gamma_s\) and the dot product score \(\Delta_s\).

Parameters:

Name Type Description Default
nb Notebook

Notebook containing at least call_spots page.

 required
score_omp_multiplier Optional[float]

Can specify the value of score_omp_multiplier to use to compute omp score. If None, will use value in config file.

 None
Source code in coppafish/plot/omp/score_hist.py 
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
def __init__(self, nb: Notebook, score_omp_multiplier: Optional[float] = None):
    """
    This plots the bivariate histogram to see the correlation between the omp spot score, $\gamma_s$ and
    the dot product score $\Delta_s$.

    Args:
        nb: *Notebook* containing at least `call_spots` page.
        score_omp_multiplier: Can specify the value of score_omp_multiplier to use to compute omp score.
            If `None`, will use value in config file.
    """
    # large hist_spacing so quick as we change it anway
    super().__init__(nb, 'omp', score_omp_multiplier, False, 0.5, False)
    self.ax.clear()
    # Get rid of buttons - only use actual dot product score
    self.buttons_ax.clear()
    plt.axis('off')
    self.score = self.score[:, [0, self.n_plots - 1]]
    self.n_plots = 2
    del self.plots
    hist_bins = np.arange(self.hist_min, self.hist_max + self.hist_spacing / 2, self.hist_spacing)
    self.x_score_ind = 0
    self.hist_spacing = 0.01
    self.plot = self.ax.hist2d(self.score[:, self.x_score_ind], self.score[:, -1], hist_bins)[3]
    self.cbar = self.fig.colorbar(self.plot, ax=self.ax)
    self.ax.set_xlim(self.hist_min, self.hist_max)
    self.ax.set_ylim(self.hist_min, self.hist_max)
    self.text_boxes[1].set_val(self.hist_spacing)
    self.ax.set_xlabel(self.button_labels[0].replace('\n', ', '))
    self.ax.set_ylabel(self.button_labels[-1].replace('\n', ', '))
    plt.show()

Track Fit

get_track_info(nb, spot_no, method, dp_thresh=None, max_genes=None)

This runs omp while tracking the residual at each stage.

Parameters:

Name Type Description Default
nb Notebook

Notebook containing experiment details. Must have run at least as far as call_reference_spots.

 required
spot_no int

Spot of interest to get track_info for.

 required
method str

'anchor' or 'omp'. Which method of gene assignment used i.e. spot_no belongs to ref_spots or omp page of Notebook.

 required
dp_thresh Optional[float]

If None, will use value in omp section of config file.

 None
max_genes Optional[int]

If None, will use value in omp section of config file.

 None

Returns:

Type Description
dict

track_info - dictionary containing info about genes added at each step returned:

  • background_codes - float [n_channels x n_rounds x n_channels]. background_codes[c] is the background vector for channel c with L2 norm of 1.
  • background_coefs - float [n_channels]. background_coefs[c] is the coefficient value for background_codes[c].
  • gene_added - int [n_genes_added + 2]. gene_added[0] and gene_added[1] are -1. gene_added[2+i] is the ith gene that was added.
  • residual - float [(n_genes_added + 2) x n_rounds x n_channels]. residual[0] is the initial pixel_color. residual[1] is the post background pixel_color. residual[2+i] is the pixel_color after removing gene gene_added[2+i].
  • coef - float [(n_genes_added + 2) x n_genes]. coef[0] and coef[1] are all 0. coef[2+i] are the coefficients for all genes after the ith gene has been added.
  • dot_product - float [n_genes_added + 2]. dot_product[0] and dot_product[1] are 0. dot_product[2+i] is the dot product for the gene gene_added[2+i].
  • inverse_var - float [(n_genes_added + 2) x n_rounds x n_channels]. inverse_var[0] and inverse_var[1] are all 0. inverse_var[2+i] is the weighting used to compute dot_product[2+i], which down-weights rounds/channels for which a gene has already been fitted.
np.ndarray

bled_codes - float [n_genes x n_use_rounds x n_use_channels]. gene bled_codes used in omp with L2 norm = 1.
float

dp_thresh - threshold dot product score, above which gene is fitted.
Source code in coppafish/plot/omp/track_fit.py 
 7
 8
 9
10
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
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
def get_track_info(nb: Notebook, spot_no: int, method: str, dp_thresh: Optional[float] = None,
                   max_genes: Optional[int] = None) -> Tuple[dict, np.ndarray, float]:
    """
    This runs omp while tracking the residual at each stage.

    Args:
        nb: Notebook containing experiment details. Must have run at least as far as `call_reference_spots`.
        spot_no: Spot of interest to get track_info for.
        method: `'anchor'` or `'omp'`.
            Which method of gene assignment used i.e. `spot_no` belongs to `ref_spots` or `omp` page of Notebook.
        dp_thresh: If None, will use value in omp section of config file.
        max_genes: If None, will use value in omp section of config file.

    Returns:
        `track_info` - dictionary containing info about genes added at each step returned:

            - `background_codes` - `float [n_channels x n_rounds x n_channels]`.
                `background_codes[c]` is the background vector for channel `c` with L2 norm of 1.
            - `background_coefs` - `float [n_channels]`.
                `background_coefs[c]` is the coefficient value for `background_codes[c]`.
            - `gene_added` - `int [n_genes_added + 2]`.
                `gene_added[0]` and `gene_added[1]` are -1.
                `gene_added[2+i]` is the `ith` gene that was added.
            - `residual` - `float [(n_genes_added + 2) x n_rounds x n_channels]`.
                `residual[0]` is the initial `pixel_color`.
                `residual[1]` is the post background `pixel_color`.
                `residual[2+i]` is the `pixel_color` after removing gene `gene_added[2+i]`.
            - `coef` - `float [(n_genes_added + 2) x n_genes]`.
                `coef[0]` and `coef[1]` are all 0.
                `coef[2+i]` are the coefficients for all genes after the ith gene has been added.
            - `dot_product` - `float [n_genes_added + 2]`.
                `dot_product[0]` and `dot_product[1]` are 0.
                `dot_product[2+i]` is the dot product for the gene `gene_added[2+i]`.
            - `inverse_var` - `float [(n_genes_added + 2) x n_rounds x n_channels]`.
                `inverse_var[0]` and `inverse_var[1]` are all 0.
                `inverse_var[2+i]` is the weighting used to compute `dot_product[2+i]`,
                 which down-weights rounds/channels for which a gene has already been fitted.
        `bled_codes` - `float [n_genes x n_use_rounds x n_use_channels]`.
            gene `bled_codes` used in omp with L2 norm = 1.
        `dp_thresh` - threshold dot product score, above which gene is fitted.
    """
    color_norm = nb.call_spots.color_norm_factor[np.ix_(nb.basic_info.use_rounds,
                                                        nb.basic_info.use_channels)]
    n_use_rounds, n_use_channels = color_norm.shape
    if method.lower() == 'omp':
        page_name = 'omp'
        config_name = 'omp'
    else:
        page_name = 'ref_spots'
        config_name = 'call_spots'
    spot_color = nb.__getattribute__(page_name).colors[spot_no][
                     np.ix_(nb.basic_info.use_rounds, nb.basic_info.use_channels)] / color_norm
    n_genes = nb.call_spots.bled_codes_ge.shape[0]
    bled_codes = np.asarray(
        nb.call_spots.bled_codes_ge[np.ix_(np.arange(n_genes),
                                           nb.basic_info.use_rounds, nb.basic_info.use_channels)])
    # ensure L2 norm is 1 for bled codes
    norm_factor = np.expand_dims(np.linalg.norm(bled_codes, axis=(1, 2)), (1, 2))
    norm_factor[norm_factor == 0] = 1  # For genes with no dye in use_dye, this avoids blow up on next line
    bled_codes = bled_codes / norm_factor

    # Get info to run omp
    dp_norm_shift = nb.call_spots.dp_norm_shift * np.sqrt(n_use_rounds)
    config = nb.get_config()
    if dp_thresh is None:
        dp_thresh = config['omp']['dp_thresh']
    alpha = config[config_name]['alpha']
    beta = config[config_name]['beta']
    if max_genes is None:
        max_genes = config['omp']['max_genes']
    weight_coef_fit = config['omp']['weight_coef_fit']

    # Run omp with track to get residual at each stage
    track_info = get_all_coefs(spot_color[np.newaxis], bled_codes, nb.call_spots.background_weight_shift,
                               dp_norm_shift, dp_thresh, alpha, beta, max_genes, weight_coef_fit, True)[2]
    return track_info, bled_codes, dp_thresh