This notebook is aimed showcasing the iter_subplots function. Some nice dual-color super-resolution data on a large number of cells would be perfect for this. unfortunately, I don’t have any, so we’ll first generate some.
%matplotlib notebook
import matplotlib.pyplot as plt
from colicoords import load, save, CellPlot, CellListPlot, SynthCell, IterCellPlot, iter_subplots, SynthCellList, AutoIterCellPlot, Data, Cell, CellList
from tqdm.auto import tqdm
from itertools import combinations
import os
import numpy as np
np.random.seed(43)We start out by loading some cell-objects, these were measured and analyzed by Yichen Li
cells = load('starting_cells.hdf5')
len(cells), cells.data.names(19, ['binary', 'brightfield', 'storm'])fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(6, 2))
cp = CellPlot(cells[0])
cp.imshow('binary', ax=ax1)
cp.plot_outline(ax=ax1)
cp.imshow('brightfield', ax=ax2)
cp.plot_storm(data_name='storm', ax=ax3, method='gauss', upscale=5)
for ax in [ax1, ax2, ax3]:
ax.set_xticks([])
ax.set_yticks([])<IPython.core.display.Javascript object>
HBox(children=(IntProgress(value=0, max=2104), HTML(value='')))As we can see in the cells above, we have a total of 19 cells, which have a 1-color super resolution dataset, outlining the membrane. We’d like to a) add more cells and b) add a second channel dataset.
To create more cells, we can define two functions to mirror the existing cells either horizontally or vertically:
def mirror_vertical(cell):
"""
Mirrors a Cell object verically and returns the mirrored cell
"""
c = cell.copy()
#mirroring the images by indexing
bf = c.data.data_dict['brightfield'][:, ::-1]
binary = c.data.data_dict['binary'][:, ::-1]
data = Data()
data.add_data(binary, 'binary')
data.add_data(bf, 'brightfield')
#Mirroring STORM data elements by inverting the x-axis
for k, v in c.data.storm_dict.items():
st_m = v
st_m['x'] = c.data.shape[1] - st_m['x']
data.add_data(st_m, 'storm', k)
c_out = Cell(data)
c_out.name = c.name + '_v'
return c_out
def mirror_horizontal(cell):
"""
Mirrors a Cell object horizontally and returns the mirrored cell
"""
c = cell.copy()
#mirroring the images by indexing
bf = c.data.data_dict['brightfield'][::-1, :]
binary = c.data.data_dict['binary'][::-1, :]
data = Data()
data.add_data(binary, 'binary')
data.add_data(bf, 'brightfield')
#Mirroring STORM data elements by inverting the y-axis
for k, v in c.data.storm_dict.items():
st_m = v
st_m['y'] = c.data.shape[0] - st_m['y']
data.add_data(st_m, 'storm', k)
c_out = Cell(data)
c_out.name = c.name + '_h'
return c_outNext, we want to add something that could represent the second channel. For this, we’ll generate some random locations within the cell, lets say two ‘blobs’ left and right which could represent nucleoids of the dividing cell.
We define some functions below to help us do this. The general process is to generate random values along the longitinudinal and radial coordinates of the cells, and then transforming them back to Cartesian coordinates and add them to the cell Object. We check for AssertionErrors because sometimes out-of-bounds values are generated, and we simply skip those and try again.
The gen_values function returns both negative and positive values, indicating postitions above and below the cell midline. In ColiCoords only positive values are used and the parameter phi is used to indicate the direction. In gen_nucleoid_storm this information is transferred to phi.
def gen_values(cell, l_pos, num):
"""
Generate a number of random positions in the cell centered around ``l_pos``.
"""
l_mu = l_pos*np.random.normal(1, 0.2)*cell.length
l_std = 0.1*np.random.normal(1, 0.2)*cell.length
r_mu = 0+np.random.normal(0, 0.1)*cell.radius
r_std = 0.15*np.random.normal(1, 0.2)*cell.radius
pos = np.array([[np.random.normal(l_mu, l_std), np.random.normal(r_mu, r_std)] for i in range(num)]).T
l, r = pos
return l, r
def gen_nucleoid_storm(cell, l, r):
"""
Given a Cell object and some random cellular coordinates `l` and `r`, create a STORM data element.
"""
phi = np.zeros_like(r)
phi[np.sign(r) == 1] = 180
x_res, y_res = cell.coords.rev_transform(l, np.abs(r), phi, l_norm=False)
storm = np.recarray((len(x_res, )), dtype=[('x', float), ('y', float), ('frame', int), ('intensity', int)])#
storm['x'] = x_res
storm['y'] = y_res
storm['frame'] = np.zeros_like(x_res)
storm['intensity'] = np.ones_like(x_res)
return storm
def add_nucleoids(c):
"""
Takes a Cell object and adds 'nucleoid' data elements at 0.25 l and 0.75 l.
"""
c_new = c.copy()
while True:
try:
l, r = gen_values(c, 0.25, int(np.random.normal(75, 15)))
storm1 = gen_nucleoid_storm(c_new, l, r)
break
except AssertionError:
pass
while True:
try:
l, r = gen_values(c, 0.75, int(np.random.normal(75, 15)))
storm2 = gen_nucleoid_storm(c_new, l, r)
break
except AssertionError:
pass
c_new.data.add_data(storm1, 'storm', 'nuc1')
c_new.data.add_data(storm2, 'storm', 'nuc2')
return c_newCombining all functions, we can write a for loop which multiplies the number of cells 4-fold and adds two nucleod datasets
functions = [lambda x: x, mirror_horizontal, mirror_vertical, lambda x: mirror_horizontal(mirror_vertical(x))]
total_cells = []
for f in functions:
for c in cells:
c_n = add_nucleoids(f(c))
total_cells.append(c_n)
cl = CellList(total_cells)
len(cl), cl[0].data.names(76, ['binary', 'brightfield', 'storm', 'nuc1', 'nuc2'])We now have 76 cells with two data elements added. We can save the resulting Cell objects and display one to have a look at the result.
save('output_cells.hdf5', cl)fig, ax = plt.subplots( figsize=(6, 4))
cp = CellPlot(cl[0])
cp.imshow(np.zeros(c.data.shape), cmap='gray')
cp.plot_storm(data_name='nuc1', method='gauss', upscale=5, alpha_cutoff=0.25, cmap='inferno')
cp.plot_storm(data_name='nuc2', method='gauss', upscale=5, alpha_cutoff=0.25, cmap='inferno')
cp.plot_storm(data_name='storm', method='gauss', upscale=5, alpha_cutoff=0.25, cmap='viridis')
t = ax.set_xticks([])
t = ax.set_yticks([])<IPython.core.display.Javascript object>
HBox(children=(IntProgress(value=0, max=2104), HTML(value='')))Nice! The upscale parameter controls the number of pixels in the recontructed image. It scales the pixel size down, so setting it to 5 means that both pixel dimensions are reduced by 5 so the amount of pixels is increased by a factor 25. When setting this parameter too high the computation time for the reconstruction can get very long. The alpha_cutoff is the fraction of maximum intensity below which the values are set to maximum alpha. This is to allow for overlap of multiple colors.
fig, (ax1, ax2) = iter_subplots(1, 2, figsize=(6.1, 2))
cp = IterCellPlot(cl[:20][::-1])
cp.imshow('brightfield', ax=ax1)
cp.plot_storm(ax=ax1, data_name='storm', color='g', alpha=0.1)
cp.plot_outline(ax=ax1)
cp.imshow(np.stack([np.zeros(cp.cell_list.data.shape) for i in range(len(cp.cell_list))]), cmap='gray', ax=ax2)
cp.plot_storm(data_name='nuc1', method='gauss', upscale=5, alpha_cutoff=0.15, cmap='inferno', ax=ax2)
cp.plot_storm(data_name='nuc2', method='gauss', upscale=5, alpha_cutoff=0.15, cmap='inferno', ax=ax2)
cp.plot_storm(data_name='storm', method='gauss', upscale=5, alpha_cutoff=0.15, cmap='viridis', ax=ax2)
fig.display()
ax1.set_axis_off()
ax2.set_axis_off()
plt.tight_layout()<IPython.core.display.Javascript object>
VBox(children=(HBox(children=(Button(description='First', style=ButtonStyle()), Button(description='Prev', sty…And done! I’ve limited the amount of cells plotted to 20 because rendering all storm images takes a long time (yes I know why bother generating more). The rendering of STORM images is to be improved in a future performance update. Currently, the experimental numba branch already provides some speedup.