Quantifying wound healing/cell migration in 3D

Timo Friedrich and Jan Kaslin showed me some 3D fluorescent labelled cells invading a wound that they would like to quantify. Initially, quantification along a single axis is acceptable. Let's load a sample dataset:



In [3]:

    
%matplotlib inline
import os
import numpy as np
from matplotlib import pyplot as plt

datadir = '/Users/nuneziglesiasj/Data/zebrafish-lesion'
os.chdir(datadir)



In [7]:

    
from gala import imio
im = imio.read_image_stack('sci-quant.lzf.h5')
im.shape, im.dtype, im.min(), im.max()









    



WARNING:root: progressbar package not installed. Progress cannot be shown. See http://pypi.python.org/simple/progressbar or type "sudo easy_install progressbar" to fix.






    



pylibtiff not available: http://www.lfd.uci.edu/~gohlke/pythonlibs/#pylibtiff






    Out[7]:





((104, 18, 512, 512), dtype('uint8'), 0, 255)

Now, let's look at a slice to make sure everything looks all right.



In [8]:

    
from matplotlib import cm



In [9]:

    
import vis
cshow = vis.cshow



In [10]:

    
cshow(im)









    Out[10]:





<matplotlib.image.AxesImage at 0x10a30df10>

The idea Timo and I discussed was to use the centres of the cell masses at the top and bottom to define an axis along which to measure cell migration. Below, we try to reliably find the centres:



In [11]:

    
sample = im[52]
top = sample[:, 0, :]
bottom = sample[:, -1, :]
plt.subplot(121); cshow(top)
plt.subplot(122); cshow(bottom)









    Out[11]:





<matplotlib.image.AxesImage at 0x10a3463d0>

Whoops! The anisotropy is killing the visualisation here! let's "expand" the images to see something a bit more reasonable:



In [12]:

    
top2 = vis.expandx(top)
cshow(top2)









    Out[12]:





<matplotlib.image.AxesImage at 0x10b833150>

Much better. Now, we want the centre of that blob of cells. My suggestion is a gaussian blur, thresholding, and center-of-mass computation. Let's try that.



In [13]:

    
from scipy import ndimage as nd
topg = nd.gaussian_filter(top, sigma=[5, 25])
expandx = vis.expandx
cshow(expandx(topg))









    Out[13]:





<matplotlib.image.AxesImage at 0x10b8f9410>



In [14]:

    
np.unravel_index(np.argmax(topg), topg.shape)









    Out[14]:





(15, 138)

That looks like a pretty bloody good estimate. Let's try the same thing with the bottom:



In [15]:

    
bottomg = nd.gaussian_filter(bottom, sigma=[5, 25])
plt.subplot(121); cshow(expandx(bottom))
plt.subplot(122); cshow(expandx(bottomg))









    Out[15]:





<matplotlib.image.AxesImage at 0x10bb12650>

The "reflect" mode padding the array might be biasing us away from the actual center. Not sure how to deal with this yet but it's something to think about.

For profiling, we might want to simplify the problem to 2D. Let's see what this looks like:



In [16]:

    
plt.plot(bottomg.sum(axis=0))









    Out[16]:





[<matplotlib.lines.Line2D at 0x10bbffdd0>]

Looks pretty Gaussian! We extract the mode, and find the limits (which will define our acquisition width) using the width-at-half-maximum (WHM):



In [30]:

    
bottom_distrib = bottomg.sum(axis=0)
top_distrib = topg.sum(axis=0)
flattened_image = sample.sum(axis=0)
cshow(flattened_image)
bottom_loc = bottom_distrib.argmax()
bottom_max = bottom_distrib[bottom_loc]
bottom_whm = (bottom_distrib > (bottom_max / 2)).astype(np.uint8).sum()
top_loc = top_distrib.argmax()
top_max = top_distrib[top_loc]
top_whm = (top_distrib > (top_max / 2)).astype(np.uint8).sum()
print(top_whm, bottom_whm)

angle = np.arctan(np.abs(float(bottom_loc - top_loc)) / flattened_image.shape[0])
width = max(top_whm, bottom_whm) * np.cos(angle)

from lineprofile import profile_line
prof = profile_line(flattened_image, ((top_loc, 0), (bottom_loc, flattened_image.shape[0] - 1)), int(width))
plt.figure(); plt.plot(prof)









    



(83, 85)






    Out[30]:





[<matplotlib.lines.Line2D at 0x10c5a7410>]