Examples

Dataset overview plot

This example demonstrates basic data visualization with dclab and matplotlib. To run this script, download the reference dataset calibration_beads.rtdc [RHMG19] and place it in the same directory.

You will find more examples in the advanced usage section of this documentation.

overview_plot.py

import matplotlib.pylab as plt
import numpy as np

import dclab

# Dataset to display
DATASET_PATH = "calibration_beads.rtdc"
# Features for scatter plot
SCATTER_X = "area_um"
SCATTER_Y = "deform"
# Event index to display
EVENT_INDEX = 100

xlabel = dclab.dfn.get_feature_label(SCATTER_X)
ylabel = dclab.dfn.get_feature_label(SCATTER_Y)

ds = dclab.new_dataset(DATASET_PATH)

fig = plt.figure(figsize=(8, 7))


ax1 = plt.subplot(221, title="Simple scatter plot")
ax1.plot(ds[SCATTER_X], ds[SCATTER_Y], "o", color="k", alpha=.2, ms=1)
ax1.set_xlabel(xlabel)
ax1.set_ylabel(ylabel)
ax1.set_xlim(19, 40)
ax1.set_ylim(0.005, 0.03)

ax2 = plt.subplot(222, title="KDE scatter plot")
sc = ax2.scatter(ds[SCATTER_X], ds[SCATTER_Y],
                 c=ds.get_kde_scatter(xax=SCATTER_X,
                                      yax=SCATTER_Y,
                                      kde_type="multivariate"),
                 s=3)
plt.colorbar(sc, label="kernel density [a.u]", ax=ax2)
ax2.set_xlabel(xlabel)
ax2.set_ylabel(ylabel)
ax2.set_xlim(19, 40)
ax2.set_ylim(0.005, 0.03)

ax3 = plt.subplot(425, title="Event image with contour")
ax3.imshow(ds["image"][EVENT_INDEX], cmap="gray")
ax3.plot(ds["contour"][EVENT_INDEX][:, 0],
         ds["contour"][EVENT_INDEX][:, 1],
         c="r")
ax3.set_xlabel("Detector X [px]")
ax3.set_ylabel("Detector Y [px]")

ax4 = plt.subplot(427, title="Event mask with µm-scale")
pxsize = ds.config["imaging"]["pixel size"]
ax4.imshow(ds["mask"][EVENT_INDEX],
           extent=[0, ds["mask"].shape[2] * pxsize,
                   0, ds["mask"].shape[1] * pxsize],
           cmap="gray")
ax4.set_xlabel("Detector X [µm]")
ax4.set_ylabel("Detector Y [µm]")

ax5 = plt.subplot(224, title="Fluorescence traces")
flsamples = ds.config["fluorescence"]["samples per event"]
flrate = ds.config["fluorescence"]["sample rate"]
fltime = np.arange(flsamples) / flrate * 1e6
# here we plot "fl?_raw"; you may also plot "fl?_med"
ax5.plot(fltime, ds["trace"]["fl1_raw"][EVENT_INDEX],
         c="#15BF00", label="fl1_raw")
ax5.plot(fltime, ds["trace"]["fl2_raw"][EVENT_INDEX],
         c="#BF8A00", label="fl2_raw")
ax5.plot(fltime, ds["trace"]["fl3_raw"][EVENT_INDEX],
         c="#BF0C00", label="fl3_raw")
ax5.legend()
ax5.set_xlim(ds["fl1_pos"][EVENT_INDEX] - 2*ds["fl1_width"][EVENT_INDEX],
             ds["fl1_pos"][EVENT_INDEX] + 2*ds["fl1_width"][EVENT_INDEX])
ax5.set_xlabel("Event time [µs]")
ax5.set_ylabel("Fluorescence [a.u.]")

plt.tight_layout()

plt.show()

Young’s modulus computation from data on DCOR

This example reproduces the lower right subplot of figure 10 in [Her17]. It illustrates how the Young’s modulus of elastic beads can be retrieved correctly (independent of the flow rate, with correction for pixelation and shear-thinning) using the area-deformation look-up table implemented in dclab (right plot). For comparison, the flow-rate-dependent deformation is also shown (left plot).

Note that this example uses the ‘buyukurganci-2022’ model for computing the viscosity, which was introduced in dclab 0.48.0.

The dataset is loaded directly from DCOR and thus an active internet connection is required for this example.

emodulus_dcor.py

import dclab
import matplotlib.pylab as plt

# The dataset is also available on figshare
# (https://doi.org/10.6084/m9.figshare.12721436.v1), but we
# are accessing it through the DCOR API, because we do not
# have the time to download the entire dataset. The dataset
# name is figshare-12721436-v1. These are the resource IDs:
ds_loc = ["e4d59480-fa5b-c34e-0001-46a944afc8ea",
          "2cea205f-2d9d-26d0-b44c-0a11d5379152",
          "2cd67437-a145-82b3-d420-45390f977a90",
          ]
ds_list = []  # list of opened datasets
labels = []  # list of flow rate labels

# load the data
for loc in ds_loc:
    ds = dclab.new_dataset(loc)
    labels.append("{:.2f}".format(ds.config["setup"]["flow rate"]))
    # emodulus computation
    ds.config["calculation"]["emodulus lut"] = "LE-2D-FEM-19"
    ds.config["calculation"]["emodulus medium"] = ds.config["setup"]["medium"]
    ds.config["calculation"]["emodulus temperature"] = \
        ds.config["setup"]["temperature"]
    ds.config["calculation"]["emodulus viscosity model"] = 'buyukurganci-2022'
    # filtering
    ds.config["filtering"]["area_ratio min"] = 1.0
    ds.config["filtering"]["area_ratio max"] = 1.1
    ds.config["filtering"]["deform min"] = 0
    ds.config["filtering"]["deform max"] = 0.035
    # This option will remove "nan" events that appear in the "emodulus"
    # feature. If you are not working with DCOR, this might lead to a
    # longer computation time, because all available features are
    # computed locally. For data on DCOR, this computation already has
    # been done.
    ds.config["filtering"]["remove invalid events"] = True
    ds.apply_filter()
    # Create a hierarchy child for convenience reasons
    # (Otherwise we would have to do e.g. ds["deform"][ds.filter.all]
    # everytime we need to access a feature)
    ds_list.append(dclab.new_dataset(ds))

# plot
fig = plt.figure(figsize=(8, 4))

# box plot for deformation
ax1 = plt.subplot(121)
ax1.set_ylabel(dclab.dfn.get_feature_label("deform"))
data_deform = [di["deform"] for di in ds_list]
# Uncomment this line if you are not filtering invalid events (above)
# data_deform = [d[~np.isnan(d)] for d in data_deform]
bplot1 = ax1.boxplot(data_deform,
                     vert=True,
                     patch_artist=True,
                     labels=labels,
                     )

# box plot for Young's modulus
ax2 = plt.subplot(122)
ax2.set_ylabel(dclab.dfn.get_feature_label("emodulus"))
data_emodulus = [di["emodulus"] for di in ds_list]
# Uncomment this line if you are not filtering invalid events (above)
# data_emodulus = [d[~np.isnan(d)] for d in data_emodulus]
bplot2 = ax2.boxplot(data_emodulus,
                     vert=True,
                     patch_artist=True,
                     labels=labels,
                     )

# colors
colors = ["#0008A5", "#A5008D", "#A50100"]
for bplot in (bplot1, bplot2):
    for patch, color in zip(bplot['boxes'], colors):
        patch.set_facecolor(color)

# axes
for ax in [ax1, ax2]:
    ax.grid()
    ax.set_xlabel("flow rate [µL/s]")

plt.tight_layout()
plt.show()

Viscosity models for Young’s modulus estimation

This example visualizes the different viscosity models for the MC-PBS media implemented in dclab. We reproduce the lower left part of figure 3 in [RB23] (channel width is 20 µm).

viscosity_models.py

import matplotlib.pylab as plt
import matplotlib.lines as mlines
from matplotlib import cm
import numpy as np

from dclab.features.emodulus import viscosity


visc_res = {}

for medium in ["0.49% MC-PBS", "0.59% MC-PBS"]:
    visc_her = {}
    visc_buy = {}

    kwargs = {
        "medium": medium,
        "channel_width": 20.0,
        "temperature": np.linspace(19, 37, 100, endpoint=True),
    }

    flow_rate = np.arange(0.02, 0.13, 0.02)

    for fr in flow_rate:
        visc_her[fr] = viscosity.get_viscosity_mc_pbs_herold_2017(
            flow_rate=fr, **kwargs)
        visc_buy[fr] = viscosity.get_viscosity_mc_pbs_buyukurganci_2022(
            flow_rate=fr, **kwargs)

    visc_res[medium] = [visc_her, visc_buy]


fig, axes = plt.subplots(1, 2, figsize=(8, 4), sharey="all", sharex="all")
colors = [cm.get_cmap("viridis")(x) for x in np.linspace(.8, 0,
                                                         len(flow_rate))]

for ii, medium in enumerate(visc_res):
    visc_her, visc_buy = visc_res[medium]
    ax = axes.flatten()[ii]
    ax.set_title(medium)

    for jj, fr in enumerate(flow_rate):
        ax.plot(kwargs["temperature"], visc_her[fr], color=colors[jj], ls="--")
        ax.plot(kwargs["temperature"], visc_buy[fr], color=colors[jj], ls="-")

    ax.set_xlabel("Temperature [°C]")
    ax.set_ylabel("Viscosity [mPa·s]")
    ax.grid()
    ax.set_ylim(2, 12)

handles = []
for jj, fr in enumerate(flow_rate):
    handles.append(
        mlines.Line2D([], [], color=colors[jj], label=f'{fr:.4g} µL/s'))
handles.append(
    mlines.Line2D([], [], color='gray', label='Büyükurgancı 2022'))
handles.append(
    mlines.Line2D([], [], color='gray', ls="--", label='Herold 2017'))
axes[0].legend(handles=handles)


plt.tight_layout()
plt.show()

lme4: Linear mixed-effects models

We would like to quantify the difference between human skeletal stem cells (SSC) and the human osteosarcoma cell line MG-63 (which is often used as a model system for SSCs) using a likelihood ratio test based on LMM.

This example illustrates a basic LMM analysis. The data are loaded from DCOR ([XRM+20], DCOR:figshare-11662773-v2). We treat SSC as our “treatment” and MG-63 as our “control” group. These are just names that remind us that we are comparing one type of sample against another type.

We are interested in the p-value, which is 0.01256 for deformation. We repeat the analysis with area (0.0002183) and Young’s modulus (0.0002771). The p-values indicate that MG-63 (mean elastic modulus 1.26 kPa) cells are softer than SSCs (mean elastic modulus 1.54 kPa). The figure reproduces the last subplot of figure 6b im [HMMO18].

lme4_lmer.py

import dclab
from dclab import lme4

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt


# https://dcor.mpl.mpg.de/dataset/figshare-11662773-v2
# SSC_16uls_rep1_20150611.rtdc
ds_ssc_rep1 = dclab.new_dataset("86cc5a47-364b-cf58-f9e3-cc114dd38e55")
# SSC_16uls_rep2_20150611.rtdc
ds_ssc_rep2 = dclab.new_dataset("ab95c914-0311-6a46-4eba-8fabca7d27d6")
# MG63_pure_16uls_rep1_20150421.rtdc
ds_mg63_rep1 = dclab.new_dataset("42cb33d4-2f7c-3c22-88e1-b9102d64d7e9")
# MG63_pure_16uls_rep2_20150422.rtdc
ds_mg63_rep2 = dclab.new_dataset("a4a98fcb-1de1-1048-0efc-b0a84d4ab32e")
# MG63_pure_16uls_rep3_20150422.rtdc
ds_mg63_rep3 = dclab.new_dataset("0a8096ce-ea7a-e36d-1df3-42c7885cd71c")

datasets = [ds_ssc_rep1, ds_ssc_rep2, ds_mg63_rep1, ds_mg63_rep2, ds_mg63_rep3]
for ds in datasets:
    # perform filtering
    ds.config["filtering"]["area_ratio min"] = 0
    ds.config["filtering"]["area_ratio max"] = 1.05
    ds.config["filtering"]["area_um min"] = 120
    ds.config["filtering"]["area_um max"] = 550
    ds.config["filtering"]["deform min"] = 0
    ds.config["filtering"]["deform max"] = 0.1
    ds.apply_filter()
    # enable computation of Young's modulus
    ds.config["calculation"]["emodulus lut"] = "LE-2D-FEM-19"
    ds.config["calculation"]["emodulus medium"] = "CellCarrier"
    ds.config["calculation"]["emodulus temperature"] = 23.0
    ds.config["calculation"]["emodulus viscosity model"] = 'buyukurganci-2022'

# setup lme4 analysis
rlme4 = lme4.Rlme4(model="lmer")
rlme4.add_dataset(ds_ssc_rep1, group="treatment", repetition=1)
rlme4.add_dataset(ds_ssc_rep2, group="treatment", repetition=2)
rlme4.add_dataset(ds_mg63_rep1, group="control", repetition=1)
rlme4.add_dataset(ds_mg63_rep2, group="control", repetition=2)
rlme4.add_dataset(ds_mg63_rep3, group="control", repetition=3)

# perform analysis for deformation
for feat in ["area_um", "deform", "emodulus"]:
    res = rlme4.fit(feature=feat)
    print("Results for {}:".format(feat))
    print("  p-value", res["anova p-value"])
    print("  mean of MG-63", res["fixed effects intercept"])
    print("  fixed effect size", res["fixed effects treatment"])

# prepare for plotting
df = pd.DataFrame()
for ds in datasets:
    group = ds.config["experiment"]["sample"].split()[0]
    rep = ds.config["experiment"]["sample"].split()[-1]
    dfi = pd.DataFrame.from_dict(
        {"area_m": ds["area_um"][ds.filter.all],
         "deform": ds["deform"][ds.filter.all],
         "emodulus": ds["emodulus"][ds.filter.all],
         "group and repetition": [group + " " + rep] * ds.filter.all.sum(),
         "group": [group] * ds.filter.all.sum(),
         })
    df = df.append(dfi)

# plot
fig = plt.figure(figsize=(8, 5))
ax = sns.boxplot(x="group and repetition", y="emodulus", data=df, hue="group")
# note that `res` is still the result for "emodulus"
numstars = sum([res["anova p-value"] < .05,
                res["anova p-value"] < .01,
                res["anova p-value"] < .001,
                res["anova p-value"] < .0001])
# significance bars
h = .1
y1 = 6
y2 = 4.2
y3 = 6.2
ax.plot([-.5, -.5, 1, 1], [y1, y1+h, y1+h, y1], lw=1, c="k")
ax.plot([2, 2, 4.5, 4.5], [y2, y2+h, y2+h, y2], lw=1, c="k")
ax.plot([.25, .25, 3.25, 3.25], [y1+h, y1+2*h, y1+2*h, y2+h], lw=1, c="k")
ax.text(2, y3, "*"*numstars, ha='center', va='bottom', color="k")
ax.set_ylim(0, 7)

plt.tight_layout()
plt.show()

lme4: Generalized linear mixed-effects models with differential deformation

This example illustrates how to perform a differential feature (including reservoir data) GLMM analysis. The example data are taken from DCOR ([XRM+20], DCOR:figshare-11662773-v2). As in the previous example, we treat SSC as our “treatment” and MG-63 as our “control” group.

The p-value for the differential deformation is magnitudes lower than the p-value for the (non-differential) deformation in the previous example. This indicates that there is a non-negligible initial deformation of the cells in the reservoir.

lme4_glmer_diff.py

from dclab import lme4, new_dataset

# https://dcor.mpl.mpg.de/dataset/figshare-11662773-v2
datasets = [
    # SSC channel
    [new_dataset("86cc5a47-364b-cf58-f9e3-cc114dd38e55"), "treatment", 1],
    [new_dataset("ab95c914-0311-6a46-4eba-8fabca7d27d6"), "treatment", 2],
    # SSC reservoir
    [new_dataset("761ab515-0416-ede8-5137-135c1682580c"), "treatment", 1],
    [new_dataset("3b83d47b-d860-4558-51d6-dcc524f5f90d"), "treatment", 2],
    # MG-63 channel
    [new_dataset("42cb33d4-2f7c-3c22-88e1-b9102d64d7e9"), "control", 1],
    [new_dataset("a4a98fcb-1de1-1048-0efc-b0a84d4ab32e"), "control", 2],
    [new_dataset("0a8096ce-ea7a-e36d-1df3-42c7885cd71c"), "control", 3],
    # MG-63 reservoir
    [new_dataset("56c449bb-b6c9-6df7-6f70-6744b9960980"), "control", 1],
    [new_dataset("387b5ac9-1cc6-6cac-83d1-98df7d687d2f"), "control", 2],
    [new_dataset("7ae49cd7-10d7-ef35-a704-72443bb32da7"), "control", 3],
]

# perform filtering
for ds, _, _ in datasets:
    ds.config["filtering"]["area_ratio min"] = 0
    ds.config["filtering"]["area_ratio max"] = 1.05
    ds.config["filtering"]["area_um min"] = 120
    ds.config["filtering"]["area_um max"] = 550
    ds.config["filtering"]["deform min"] = 0
    ds.config["filtering"]["deform max"] = 0.1
    ds.apply_filter()

# perform LMM analysis for differential deformation
# setup lme4 analysis
rlme4 = lme4.Rlme4(feature="deform")
for ds, group, repetition in datasets:
    rlme4.add_dataset(ds, group=group, repetition=repetition)

# LMM
lmer_result = rlme4.fit(model="lmer")
print("LMM p-value", lmer_result["anova p-value"])  # 0.00000351

# GLMM with log link function
glmer_result = rlme4.fit(model="glmer+loglink")
print("GLMM p-value", glmer_result["anova p-value"])  # 0.000868

Plotting isoelastics

This example illustrates how to plot dclab isoelastics by reproducing figure 3 (lower left) of [MMM+17].

Warning

This example does not work anymore since dclab 0.46.0, because the isoelasticity lines of the analytical model have different Young’s moduli than the ones of the revised LE-2D-FEM-19 model. For the sake of completeness, we keep this example here. If you would like to extract lines at specific Young’s moduli, please take a look at the next example.

isoelastics.py

import matplotlib.pylab as plt
import matplotlib.lines as mlines
from matplotlib import cm
import numpy as np

import dclab

# parameters for isoelastics
kwargs = {"col1": "area_um",  # x-axis
          "col2": "deform",  # y-axis
          "channel_width": 20,  # [um]
          "flow_rate": 0.04,  # [ul/s]
          "viscosity": 15,  # [mPa s]
          "add_px_err": False  # no pixelation error
          }

isos = dclab.isoelastics.get_default()
analy = isos.get(lut_identifier="LE-2D-ana-18", **kwargs)
numer = isos.get(lut_identifier="LE-2D-FEM-19", **kwargs)

plt.figure(figsize=(8, 4))
ax = plt.subplot(111, title="elastic sphere isoelasticity lines")
colors = [cm.get_cmap("jet")(x) for x in np.linspace(0, 1, len(analy))]
for aa, nn, cc in zip(analy, numer, colors):
    ax.plot(aa[:, 0], aa[:, 1], color=cc)
    ax.plot(nn[:, 0], nn[:, 1], color=cc, ls=":")

line = mlines.Line2D([], [], color='k', label='analytical')
dotted = mlines.Line2D([], [], color='k', ls=":", label='numerical')
ax.legend(handles=[line, dotted])

ax.set_xlim(50, 240)
ax.set_ylim(0, 0.02)
ax.set_xlabel(dclab.dfn.get_feature_label("area_um"))
ax.set_ylabel(dclab.dfn.get_feature_label("deform"))

plt.tight_layout()
plt.show()

Plotting custom isoelastics

This example illustrates how to extract custom isoelasticity lines from the dclab look-up tables by reproducing figure 3 (right) of [WRM+22].

Note that at the boundary of the support of a look-up table, the isoelasticity lines may break away in perpendicular directions. The underlying reason is that the look-up table is first mapped onto a grid from which the constant isoelasticity lines are extracted. Since the Young’s modulus values are linearly interpolated from the LUT onto that grid, there can be inaccuracies for pixels that are at the LUT boundary.

An elaborate way of getting rid of these inaccuracies (and this is how the isoelasticity lines for dclab are extracted), is to extend the LUT by fitting a polynomial to isoelasticity lines which are well-defined within the LUT and extrapolating these lines beyond the boundary of the LUT. This technique is documented in the scripts directory of the dclab repository.

A quicker and much less elaborate way of getting around this issue is to simply crop the individual isoelasticity lines where necessary.

isoelastics_custom.py

import matplotlib.pylab as plt
import numpy as np
import skimage

import dclab
from dclab.features import emodulus


colors = ["r", "b"]
linestyles = [":", "-"]

plt.figure(figsize=(8, 4))
ax = plt.subplot(111,
                 title="Comparison of the isoelasticity lines of two LUTs")

grid_sie = 250

for ii, lut_name in enumerate(["LE-2D-FEM-19", "HE-3D-FEM-22"]):
    area_um = np.linspace(0, 350, grid_sie, endpoint=True)
    deform = np.linspace(0, 0.2, grid_sie, endpoint=True)
    area_um_grid, deform_grid = np.meshgrid(area_um, deform, indexing="ij")

    emod = emodulus.get_emodulus(area_um=area_um_grid,
                                 deform=deform_grid,
                                 medium=6.0,
                                 channel_width=20,
                                 flow_rate=0.04,
                                 px_um=0,
                                 temperature=None,
                                 visc_model=None,
                                 lut_data=lut_name)

    levels = [0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 3.0, 6.0]
    for level in levels:
        conts = skimage.measure.find_contours(emod, level=level)
        if not conts:
            continue
        # get the longest one
        idx = np.argmax([len(cc) for cc in conts])
        cc = conts[idx]
        # remove nan values
        cc = cc[~np.isnan(np.sum(cc, axis=1))]
        # scale isoelastics back
        cc_sc = np.copy(cc)
        cc_sc[:, 0] = cc[:, 0] / grid_sie * 350
        cc_sc[:, 1] = cc[:, 1] / grid_sie * 0.2
        plt.plot(cc_sc[:, 0], cc_sc[:, 1],
                 color=colors[ii],
                 ls=linestyles[ii],
                 label=lut_name if level == levels[0] else None)

ax.set_ylim(-0.005, 0.1)
ax.set_xlabel(dclab.dfn.get_feature_label("area_um"))
ax.set_ylabel(dclab.dfn.get_feature_label("deform"))
plt.legend()
plt.tight_layout()
plt.show()

Working with plugin features

This example shows how to load a user-defined plugin feature recipe in dclab and use it in a scatter plot.

Please also download the plugin_example.py file for this example.

plugin_usage.py

import pathlib

import matplotlib.pyplot as plt

import dclab


plugin_path = pathlib.Path(__file__).parent

# load a single plugin feature
dclab.load_plugin_feature(plugin_path / "plugin_example.py")

# load some data from DCOR
ds = dclab.new_dataset("fb719fb2-bd9f-817a-7d70-f4002af916f0")

# access the features
circ_per_area = ds["circ_per_area"]
circ_times_area = ds["circ_times_area"]

# create a plot with a plugin feature
plt.figure(figsize=(8, 4))
xlabel = dclab.dfn.get_feature_label("circ_times_area")
ylabel = dclab.dfn.get_feature_label("deform")

ax1 = plt.subplot(title="Plot with a plugin feature")
ax1.plot(ds["circ_times_area"], ds["deform"],
         "o", color="k", alpha=.2, ms=1)
ax1.set_xlabel(xlabel)
ax1.set_ylabel(ylabel)
ax1.set_xlim(20, 40)
ax1.set_ylim(0.0025, 0.025)

plt.tight_layout()
plt.show()