<jupyter_start><jupyter_text>The purpose of this notebook is to demonstrate the use of the diff_register package for cell morphology analyses including the entire process from raw microscopy images to quantification and statistical analyses on cell morphological properties. By the end of the demonstration, users will be able to* Use filters to **clean up** raw images* Create **binary images** from input images* Perform **morphological operations** on binary images such as dilation and erosion* **Segment and label** individual cells from binary images* **Output dataframes** of cell morphological quantities Before jumping into the body of the analysis, I would like to help conceptualize the process of morphology analysis into two sets of functions: * Input raw image file, output binary image* Input binary image, output cell morphological quantitiesIt's helpful in data science to think in terms of inputs and outputs. Diff_classifier provides two powerhouse functions that performs these two steps for you, but users can easily create their own functions using the principles demo'd in the previous training notebook to refine their analyses. Mix and match as you see fit.  Importing and Prepping Images We're going to use Mike's CD11b images, since they're so pretty. We'll read them in from AWS using my diff_classifier `aws` module.<jupyter_code>from diff_classifier import aws as aws
filename = 'ipsilateral_s1_t3.tif'
folder = '12_11_18_CD11b_microglia_images'
bucket = 'mckenna.data'

#aws.download_s3('{}/{}'.format(folder, filename), filename, bucket_name=bucket)<jupyter_output><empty_output><jupyter_text>Let's read in the file using `scikit-image`'s `imread` command.<jupyter_code>import skimage.io as sio
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
image = sio.imread(filename)<jupyter_output><empty_output><jupyter_text>Let's check some of the image properties before going further to make sure we know what we're dealing with.<jupyter_code>image.shape<jupyter_output><empty_output><jupyter_text>Ah, ok, so we have a 2D image with two channels. Let's see what's in each one.<jupyter_code>sio.imshow(image[0, :, :], cmap='gray', vmin=0, vmax=100)<jupyter_output><empty_output><jupyter_text>The DAPI channel doesn't seem too exciting, but let's check out the GFAP channel:<jupyter_code>sio.imshow(image[1, :, :], cmap='gray', vmin=80, vmax=200)<jupyter_output><empty_output><jupyter_text>Now we're getting something exciting! Before we do anything else, let's create a variable that only has the channel of interest for easier manipulation.<jupyter_code>mglia = image[1, :, :]<jupyter_output><empty_output><jupyter_text>Binarizing Input Images The final flow of diff_register isn't set in stone yet. In fact, I significantly edited the flow this weekend to completely leave the binary step open to the user, because my self-built functions were a bit clunky, and its easier for the user to directly manipulate binarization steps in the notebook itself. I want you to toy with the binarization steps, maybe even build your own inspired from a few methods below. All these use the commands found in the training notebook, including:* Filters (e.g. mean and median)* Edge detectors* Background subtraction* Simple thresholding* Opening and closing operations* Fill homes* Remove small objects<jupyter_code>from skimage.filters import sobel, prewitt, scharr, gaussian, median, rank
from skimage.morphology import square, opening, closing, skeletonize, star, disk, dilation
from skimage.morphology import remove_small_objects, reconstruction, white_tophat
from skimage.feature import canny
from scipy.ndimage.morphology import binary_fill_holes
from skimage.measure import regionprops, label
from skan import csr, draw
import numpy as np
from skimage.util import pad

from diff_register import im_process as imp

from scipy import ndimage as ndi
from skimage.morphology import watershed
from skimage.feature import peak_local_max<jupyter_output><empty_output><jupyter_text>The first binarization function I built looked something like this (but a lot messier-- now it looks nice and tidy, with easy one-liners of code with a defined task). The basic idea is to filter the image to get rid of noise, subtract background, threshold such that anything greater than *x* is defined as cell, fill holes, and remove small objects.<jupyter_code># Old binarization
filt_image = median(mglia, selem=star(2))
backg = median(mglia, selem=star(15))
filt_min = filt_image - 0.5*backg
biim = filt_min > 51
fullim = closing(biim, selem=star(1))
filled_im = binary_fill_holes(fullim)
clean_im1 = remove_small_objects(filled_im, min_size=350)

sio.imshow(clean_im1, cmap='gray')<jupyter_output><empty_output><jupyter_text>You may or may not like this second workflow better. Instead of using a simple thresholding technique, I use an edge detector. Edge detectors seem to catch more detail at the cell edges, but it could potentially mess with the skeletonization process downstream. Each individial "bump" will be counted as a process.<jupyter_code>#New binarization

edges = canny(mglia/255., sigma=0.001)
fullim = closing(edges, selem=star(1))
filled_im = binary_fill_holes(fullim)
clean_im2 = remove_small_objects(filled_im, min_size=350)

sio.imshow(clean_im2)
plt.axis('off')<jupyter_output><empty_output><jupyter_text>The final binarization workflow is based off a design I used for studying diffusion in and out of cells in *in vitro* cell culture. This particular implementation doesn't capture a lot of the winding at cell edges with may or may not be good for your downstream analysis.<jupyter_code># Alternate binarization
filt_image = rank.mean(mglia, selem=star(2))
edge = prewitt(filt_image)
#closeim = closing(edge, selem=star(1))
openim = opening(edge, selem=star(2))
backg = 1*np.mean(openim)
binarim = openim > backg

seed = np.copy(binarim)
seed[1:-1, 1:-1] = binarim.max()
mask = binarim
filledim = reconstruction(seed, mask, method='erosion')
filled_im = binary_fill_holes(filledim)
clean_im3 = remove_small_objects(filled_im, min_size=350)

sio.imshow(clean_im3)
plt.axis('off')
# Implementation of watershed

padw = 150
cleanpad = pad(clean_im1, ((padw, padw),), mode='constant')
distance = ndi.distance_transform_edt(cleanpad)
local_maxi = peak_local_max(distance, indices=False, footprint=np.ones((90, 90)),
                            labels=cleanpad)
markers = ndi.label(local_maxi)[0]
labels = watershed(-distance, markers, mask=cleanpad)
labels = labels[padw:-padw, padw:-padw]

label_edge = scharr(labels)
label_edge = label_edge < 0.0000000000001
split_cells = np.logical_and(label_edge, clean_im1)
split_cells = remove_small_objects(split_cells, min_size=350)

sio.imshow(split_cells)
plt.axis('off')
from skimage.util import pad
cleanpad = pad(clean_im1, ((70, 70),), mode='constant')
test = cleanpad[70:-70, 70:-70]
sio.imshow(cleanpad)
plt.axis('off')<jupyter_output><empty_output><jupyter_text>Skeletonization In order to calculate branch features associated with each cell, we need to create a skeletonized image from the binarized image. Here's what the raw output of the `skimage.morphology.skeletonize` function looks likeafter a dilation step to make them clearer:<jupyter_code>skelim = dilation(skeletonize(clean_im1), selem=disk(1))
#skelim = skeletonize(clean_im1)
sio.imshow(skelim, cmap='gray')
plt.axis('off')<jupyter_output><empty_output><jupyter_text>Make sure to compare how the skeletons differ depending on the binary image you use.<jupyter_code>skelim = dilation(skeletonize(clean_im2), selem=disk(1))
#skelim = skeletonize(clean_im1)
sio.imshow(skelim, cmap='gray')
plt.axis('off')
skelim = dilation(skeletonize(clean_im3), selem=disk(1))
#skelim = skeletonize(clean_im1)
sio.imshow(skelim, cmap='gray')
plt.axis('off')<jupyter_output><empty_output><jupyter_text>I have wrapped the skeletonization process up with the calculation of branch properties in a single function called `imp.labelandskel`. It creates a `class` object with the skeletonized image and the output properties.<jupyter_code>skel = imp.labelandskel(clean_im)<jupyter_output><empty_output><jupyter_text>Morphological Analysis The final function `imp.mglia_features` unpacks the skeletonized object into a human-readable array using `pandas`. And wallah! Data!<jupyter_code>Dfeatures = imp.mglia_features(skel, umppx=1.24)<jupyter_output><empty_output><jupyter_text>Let's take a look at the data using the `head` function.<jupyter_code>Dfeatures.head()<jupyter_output><empty_output><jupyter_text>We can get a quick summary of the stats using the `describe` function:<jupyter_code>Dfeatures.describe()<jupyter_output><empty_output><jupyter_text>And for a quick look at the distribution we can use the `plot` or `hist` functions:<jupyter_code>Dfeatures['total_branches'].plot('hist', bins=np.linspace(0, 50, 11))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08686369254972469, 'fraction_numerical': 0.013709405551185527, 'mean_word_length': 6.759372275501308, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 71, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '295107', 'license': 'MIT', 'n_tokens_mistral': 2870, 'n_tokens_neox': 2706, 'n_words': 941, 'path': 'notebooks/cell_morphology_demo.ipynb', 'repo_name': 'ccurtis7/diff_register', 'size': '8899'}