Complete analysis with Skan#

In this notebook, we reproduce the analysis leading to Figure 2 of our paper, submitted to PeerJ.

0. Getting the data#

First, we need to download the input data, which is posted at the Open Science Framework (OSF) as a zip file. Many thanks to the OSF for providing hosting for scientific data. The data be downloaded at this URL:

Download and extract the zipfile, and make sure the schizonts folder contained within it is placed in this notebook’s working directory.

1. Using skan.pipe to get a branch dataset for analysis#

from glob import glob

schizont_files = glob('schizonts/*.tif')
# remove files ending with '01.tif', which are low-res overview images
schizont_files = list(filter(lambda x: not x.endswith('01.tif'),

Now we process all 96 images. This should take about a minute on a modern laptop, but could be faster on a multicore machine, as it processes images using as many CPUs as are avaliable.

from skan import pipe
import toolz as tz

data, per_image_data = tz.last(pipe.process_images(

Next, we use the filename to infer the infection status of each image.

def infection_status(filename):
    if 'Uninf' in filename:
        return 'normal'
        return 'infected'

data['infection'] = data['filename'].apply(infection_status)

Also encoded in the filename are the cell number (unique by infection status) and the field number (each cell membrane is imaged at several non-overlapping locations).

import re  # regular expressions, for text matching

def cell_number(filename):
    regex = r'.*RBC(\d+)_\d\d.tif'
    regex = re.compile(regex)
    match = regex.match(filename)
    if match is not None:
        return int(
        return None
def field_number(filename):
    regex = r'.*RBC\d+_(\d\d).tif'
    regex = re.compile(regex)
    match = regex.match(filename)
    if match is not None:
        return int(
        return None
data['cell number'] = data['filename'].apply(cell_number)
data['field number'] = data['filename'].apply(field_number)

2. Cleaning up the data#

Next, we filter the branches by using the shape index. We have used a very simple method to extract skeletons (see Getting started), which does an acceptable job but creates a lot of false branches. Since the goal of Skan is to analyse skeletons, rather than generate them, we attempt to filter the branches, and measure only those that look like ridges according to the shape index.

ridges = ((data['mean-shape-index'] < 0.625) &
          (data['mean-shape-index'] > 0.125))

For the same reason, we only look at junction-to-junction branches, which are more accurately identified by our method than junction-to-endpoint branches.

j2j = data['branch-type'] == 2
datar = data.loc[ridges & j2j].copy()

Finally, we make a new column of measurements in a more natural scale for our purpose.

datar['branch distance (nm)'] = datar['branch-distance'] * 1e9

3. Making the figure#

%matplotlib inline
%config InlineBackend.figure_format='retina'

import numpy as np
import imageio as iio
from skimage import morphology
import matplotlib.pyplot as plt
import seaborn as sns

from skan.pre import threshold
fig, axes = plt.subplots(2, 2)
ax = axes.ravel()

# display an arbitrary image
crop = (slice(20, -20),) * 2
image_raw = iio.v2.imread('schizonts/schizont4_UninfRBC7_06.tif',
image = image_raw[crop]
ax[0].imshow(image, cmap='gray')

# add a 300nm scale bar
height, width = image.shape
scalenm = image.meta['Scan']['PixelHeight'] * 10**9
ax[0].plot((1096, 1096 + 300 / scalenm), (940, 940),
           c='k', lw=2);

# Next, show the smoothed binary image and corresponding skeleton
radius = int(np.ceil(50/scalenm))
sigma = radius * 0.1
thresholded = threshold(image, sigma=sigma,
                        radius=radius, offset=0.075,
viz = np.zeros(image.shape + (3,), dtype=float)
viz[thresholded] = [1, 1, 1]
skeleton = morphology.skeletonize(thresholded)
fat_skeleton = morphology.binary_dilation(skeleton)
viz[fat_skeleton] = [1, 0, 0]
ax[1].plot((1096, 1096 + 300 / scalenm), (940, 940),
           c='g', lw=2)

# Add a panel with a histogram of all branch lengths by infection status
_, bins = np.histogram(datar['branch distance (nm)'], bins='auto')
for inf, df in (datar.sort_values(by='infection', ascending=False)
                .groupby('infection', sort=False)):
    ax[2].hist(df['branch distance (nm)'], bins=bins,
               density=True, alpha=0.5, label=inf)

ax[2].set_xlabel('branch distance (nm)')

# Finally, a panel grouping the data by cell, showing the difference
# between infected and uninfected cells
cellmeans = (datar.groupby(['infection', 'cell number'])
sns.stripplot(x='infection', y='branch distance (nm)', data=cellmeans,
              jitter=True, order=('normal', 'infected'), ax=ax[3])
ax[3].set_xlabel('infection status')
ax[3].set_ylabel('mean branch distance\nby cell (nm)')

# Use matplotlib's automatic layout algorithm
/tmp/ipykernel_3656/ FutureWarning: The default value of numeric_only in DataFrameGroupBy.mean is deprecated. In a future version, numeric_only will default to False. Either specify numeric_only or select only columns which should be valid for the function.