Yeast Dosage-Response Curve

This is data from a real experiment, measuring the dose response of an engineered yeast line to isopentyladenine, or IP. The yeast is engineered with a basic GFP reporter, and GFP fluorescence is measured after 12 hours, at which time we expect the cells to be at steady-state.

(Reference: Chen et al, Nature Biotech 2005)

Set up the Jupyter matplotlib integration, and import the cytoflow module.

# in some buggy versions of the Jupyter notebook, this needs to be in its OWN CELL.
%matplotlib inline
import cytoflow as flow

# if your figures are too big or too small, you can scale them by changing matplotlib's DPI
import matplotlib
matplotlib.rc('figure', dpi = 160)

First we set up the Experiment. This is primarily the mapping between the flow data (wells on a microtitre plate, here) and the metadata (experimental conditions, in this case the IP concentration.)

For this experiment, I did a serial dilution with three dilutions per log. Note that I specify that this metadata is of type log, which is represented internally as a numpy type float, but is plotted on a log scale when we go to visualize it.

Also note that, in order to keep the size of the source distribution manageable, I’ve only included 1000 events from each tube in the example files. The actual data set has 30,000 events from each, which gives a much prettier dose-response curve!

inputs = {
    "Yeast_B1_B01.fcs" : 5.0,
    "Yeast_B2_B02.fcs" : 3.75,
    "Yeast_B3_B03.fcs" : 2.8125,
    "Yeast_B4_B04.fcs" : 2.109,
    "Yeast_B5_B05.fcs" : 1.5820,
    "Yeast_B6_B06.fcs" : 1.1865,
    "Yeast_B7_B07.fcs" : 0.8899,
    "Yeast_B8_B08.fcs" : 0.6674,
    "Yeast_B9_B09.fcs" : 0.5,
    "Yeast_B10_B10.fcs" : 0.3754,
    "Yeast_B11_B11.fcs" : 0.2816,
    "Yeast_B12_B12.fcs" : 0.2112,
    "Yeast_C1_C01.fcs" : 0.1584,
    "Yeast_C2_C02.fcs" : 0.1188,
    "Yeast_C3_C03.fcs" : 0.0892,
    "Yeast_C4_C04.fcs" : 0.0668,
    "Yeast_C5_C05.fcs" : 0.05,
    "Yeast_C6_C06.fcs" : 0.0376,
    "Yeast_C7_C07.fcs" : 0.0282,
    "Yeast_C8_C08.fcs" : 0.0211,
    "Yeast_C9_C09.fcs" : 0.0159

tubes = []
for filename, ip in inputs.items():
    tubes.append(flow.Tube(file = "data/" + filename, conditions = {'IP' : ip}))

ex = flow.ImportOp(conditions = {'IP' : "float"},
                   tubes = tubes).apply()

Set logicle as the default plot scale:


The first step in most cytometry workflows is to filter out debris and aggregates based on morphological parameters (the forward- and side-scatter measurements.) So let’s have a look at the FSC-A and SSC-A channels.

Remember, because we’re not specifying otherwise, this is a plot of all of the data in the experiment, not a single tube!

flow.ScatterplotView(xchannel = "FSC-A",
                     ychannel = "SSC-A").plot(ex, alpha = 0.05)
/home/brian/src/cytoflow/cytoflow/utility/ CytoflowWarning: Channel FSC-A doesn't have any negative data. Try a log transform instead.

Oops! The biexponential scale is really intended for channels that have events around 0, with some negative events; let’s use a log scale on both FSC-A and SSC-A in the future. (The fluorescence channels should be fine on a logicle scale, though.)

Because this is yeast growing in log phase on a drum roller, the distributions in FSC-A and SSC-A are pretty tight. Let’s fit a 2D gaussian to it, and gate out any data that is more than three standard deviations from the mean.

g = flow.GaussianMixtureOp(name = "Debris_Filter",
                           channels = ["FSC-A", "SSC-A"],
                           scale = {'FSC-A' : 'log',
                                    'SSC-A' : 'log'},
                           num_components = 1,
                           sigma = 3)
g.default_view().plot(ex, alpha = 0.05)

ex2 = g.apply(ex)
/home/brian/src/cytoflow/cytoflow/operations/ RuntimeWarning: invalid value encountered in less_equal

One thing to note: if you’ve read the Machine Learning notebook, you know that a GaussianMixtureOp adds categorical metadata as well. In this case, we have a new boolean column named Debris_Filter_1, which is True if the event is less than sigma from the centroid and False otherwise.
FITC-A FSC-A FSC-H FSC-W IP PerCP-Cy5-5-A SSC-A SSC-H SSC-W Time Debris_Filter_1
0 2734.459961 3720.110107 3771.0 64651.589844 5.0 -306.440002 4036.360107 3915.0 67567.531250 2.3 True
1 3316.320068 6256.010254 6119.0 67003.414062 5.0 -185.179993 2631.060059 2519.0 68451.437500 3.2 True
2 3551.320068 5292.209961 4708.0 73668.281250 5.0 -234.059998 3729.919922 3376.0 72406.406250 9.2 True
3 3110.459961 5934.479980 5855.0 66425.640625 5.0 395.739990 3499.619873 3419.0 67081.335938 12.7 True
4 2560.560059 4605.700195 4631.0 65177.964844 5.0 -165.440002 3196.939941 3139.0 66745.671875 13.0 True

Now, let’s look at the FITC-A channel, which is the one we’re interested in. First of all, just check histograms of the FITC-A distribution in each tube. Make sure to set the subset trait to only include the events that aren’t debris!

flow.HistogramView(channel = "FITC-A",
                   xfacet = "IP",
                   subset = "Debris_Filter_1 == True").plot(ex2, col_wrap = 3)

Sure enough, looking at the plots, it’s clear that the histogram peak moves to the right as the IP concentration increases. We can get a more compact representation of the same data with one of cytoflow’s fancier plots, a violin plot:

flow.ViolinPlotView(variable = "IP",
                    channel = "FITC-A",
                    scale = "logicle",
                    subset = "Debris_Filter_1 == True").plot(ex2)

A quick aside - sometimes you get ugly axes because of overlapping labels. In this case, we want a wider plot. While we can directly specify the height of a plot, we can’t directly specify its width, only its aspect ratio (the ratio of the width to the height.) cytoflow defaults to 1.5; in this case, let’s widen it out to 4. If this results in a plot that’s wider than your browser, the Jupyter notebook will scale it down for you.

flow.ViolinPlotView(variable = "IP",
                    channel = "FITC-A",
                    scale = "logicle",
                    subset = "Debris_Filter_1 == True").plot(ex2, aspect = 4)

At least on my screen that’s better, but a violin plot is really best for categorical data (with a modest number of categories.) Instead, let’s make a true does-response curve. The ChannelStatisticsOp divides the data into subsets by different values of the IP variable and applies flow.geom_mean to each subset. Then, the Stats1DView plots that statistic.

(Confused? Have a look at the Statistics example notebook for a more in-depth discussion.)

ex3 = flow.ChannelStatisticOp(name = "ByIP",
                              channel = "FITC-A",
                              by = ["IP"],
                              function = flow.geom_mean,
                              subset = "Debris_Filter_1 == True").apply(ex2)

flow.Stats1DView(statistic = ("ByIP", "geom_mean"),
                 scale = "log",
                 variable = "IP",
                 variable_scale = "log").plot(ex3)

And of course, because the underlying data is just a pandas.DataFrame, we can manipulate it with the rest of the Scientific Python ecosystem. Here, we get the compute the geometric means using pandas.DataFrame.groupby.

("Debris_Filter_1 == True")[['IP', 'FITC-A']]
0.0159 110.209301
0.0211 143.507986
0.0282 157.332377
0.0376 184.959273
0.0500 250.297256
0.0668 413.781993
0.0892 678.979138
0.1188 955.664252
0.1584 1201.557123
0.2112 1164.530618
0.2816 1503.467572
0.3754 1691.126815
0.5000 1812.333210
0.6674 2219.928952
0.8899 2330.923961
1.1865 2423.731503
1.5820 2362.649599
2.1090 2341.478585
2.8125 2422.780181
3.7500 2549.603737
5.0000 2258.993753