diff --git a/FlowCal/__init__.py b/FlowCal/__init__.py index 1f7100a..11341a6 100644 --- a/FlowCal/__init__.py +++ b/FlowCal/__init__.py @@ -13,5 +13,6 @@ from . import gate from . import transform from . import mef +from . import compensate from . import plot from . import stats diff --git a/FlowCal/compensate.py b/FlowCal/compensate.py new file mode 100644 index 0000000..5151a25 --- /dev/null +++ b/FlowCal/compensate.py @@ -0,0 +1,157 @@ +""" +Functions for performing multicolor compensation. + +""" + +import functools + +import numpy as np +import FlowCal.stats + +def get_transform_fxn(nfc_sample, + sfc_samples, + comp_channels, + statistic_fxn=FlowCal.stats.mean): + """ + Get a transformation function to perform multicolor compensation. + + Parameters + ---------- + nfc_sample : FCSData object or None + Data corresponding to the no-fluorophore control sample (NFC). If + None, no autofluorescence correction will be made. + sfc_samples : list of FCSData objects + Data corresponding to the single-fluorophore control samples + (SFCs). + comp_channels : list of str or int + Channels to compensate. Each channel should correspond to an + element of `sfc_samples`. + + Returns + ------- + transform_fxn : function + Transformation function to compensate flow cytometry data. + This function has the following signature:: + + data_compensated = transform_fxn(data, channels) + + Other parameters + ---------------- + statistic_fxn : function, optional + Statistical function from ``FlowCal.stats`` used to calculate the + representative fluorescence of each control sample. + + Notes + ----- + If using MEF calibration, we recommend using calibrated NFC and SFCs, + and calibrating all samples before applying compensation. Calibration + can correct for some small nonlinearities in the instrument's + fluorescence detectors, especially in older instruments. On the other + hand, compensation requires fluorescence units proportional to the + fluorescent signal. Thus, performing compensation followed by + calibration may give slightly different results than running + calibration followed by compensation as we recommend. + + The compensation method used here is based on the following analysis. + + We assume we have an instrument with :math:`n` fluorescence channels + and a sample with :math:`n` fluorophores, where we expect the signal + in channel :math:`i` to correspond to fluorophore :math:`i` only. + In reality, the signal from each channel will additionally contain + some (hopefully small) signal from all other fluorophores. In + mathematical terms: + + .. math:: s^i = a_0^i + f_i + \\sum_{j=1,j\\neq i}^{n} a^i_j \\cdot f_j + + where: + + - :math:`s^i` is the total signal observed in channel :math:`i`. + - :math:`a^i_0` is the autofluorescence signal in channel :math:`i`. + - :math:`f_i` is the signal from fluorophore :math:`i`. + - :math:`a^i_j` is a bleedthrough coefficient, which quantifies how + much signal from fluorophore :math:`j` appears in channel :math:`i`. + + In matrix notation: + + .. math:: \\mathbf{s} = \\mathbf{a_0} + \\mathbf{A} \\cdot \\mathbf{f} + + Where :math:`\\mathbf{a_0}` is the autofluorescence vector and + :math:`\\mathbf{A}` is the bleedthrough matrix, with all diagonal terms + equal to one. For an arbitrary sample, the compensation procedure + consists on solving for :math:`\\mathbf{f}` starting from the measured + signals :math:`\\mathbf{s}`: + + .. math:: \\mathbf{f} = \\mathbf{A}^{-1} (\\mathbf{s} - \\mathbf{a_0}) + + This requires knowledge of :math:`\\mathbf{A}` and + :math:`\\mathbf{a_0}`. To find these out, we use the following control + samples: + + 1. No-fluorophore control (NFC). In this case, + :math:`\\mathbf{f}_{NFC} = 0`. Therefore, + + .. math:: \\mathbf{a_0} = \\mathbf{s}_{NFC} + + Where :math:`\\mathbf{s}_{NFC}` is the vector containing the signals in + all channels when measuring the NFC. + + 2. Single-fluorophore controls (SFCs), one for each fluorophore. For + fluorophore :math:`i`, all elements in :math:`\\mathbf{f}_{SFCi}` are + zero, except for the one at position :math:`i` (:math:`f_{SFCi}`). + Therefore, + + .. math:: \\mathbf{s}_{SFCi} = \\mathbf{a_0} + \ + \\mathbf{a}_i \\cdot f_{SFCi} + + Where :math:`\\mathbf{s}_{SFCi}` is the vector containing the signals + in all channels when measuring the SFC with fluorophore :math:`i`, and + :math:`\\mathbf{a}_i` is the ith column of :math:`A`. Solving for + :math:`\\mathbf{a}_i`: + + .. math:: \\mathbf{a}_i = (\\mathbf{s}_{SFCi} - \\mathbf{a_0})/f_{SFCi} + + Finally, using the additional restriction that :math:`a_i^i=1`, we + have: + + .. math:: f_{SFCi} = s^i_{SFCi} - a^i_0 + + Therefore + + .. math:: \\mathbf{a}_i = (\\mathbf{s}_{SFCi} - \\mathbf{a_0})/ \ + (s^i_{SFCi} - a^i_0) + + """ + + # Check for appropriate number of single fluorophore controls + if len(sfc_samples) != len(comp_channels): + ValueError('number of single fluorophore controls should match' + ' the number of channels specified') + + # Calculate autofluorescence vector + if nfc_sample is None: + a0 = np.zeros(len(comp_channels)) + else: + a0 = np.array(statistic_fxn(nfc_sample[:,comp_channels])) + # Signals from the single-fluorophore controls + # Matrix S_sfc contains the signal from an SFC in each row + # S_sfc[:, i] <= s_SFCi + S_sfc = np.array( + [np.array(statistic_fxn(s[:,comp_channels])) for s in sfc_samples]).T + # Get signal minus autofluorescence + # Each column in S_sfc_noauto contains the signal from an SFC minus the + # autofluorescence vector + # S_sfc_noauto[:, i] <= s_SFCi - a_0 + S_sfc_noauto = S_sfc - a0[:, np.newaxis] + # Calculate matrix A + # The following uses broadcasting to divide column i by the (i,i) element + # of S_sfc_noauto + # A[:, i] <= (s_SFCi - a0)/(s^i_SFCi - a^i_0) + A = S_sfc_noauto / np.diag(S_sfc_noauto) + + # Make output transformation function + transform_fxn = functools.partial(FlowCal.transform.to_compensated, + comp_channels=comp_channels, + a0=a0, + A=A) + + return transform_fxn diff --git a/FlowCal/transform.py b/FlowCal/transform.py index 13075f5..fa0da30 100644 --- a/FlowCal/transform.py +++ b/FlowCal/transform.py @@ -247,7 +247,7 @@ def to_rfi(data, return data_t -def to_mef(data, channels, sc_list, sc_channels = None): +def to_mef(data, channels, sc_list, sc_channels=None): """ Transform flow cytometry data using a standard curve function. @@ -336,3 +336,111 @@ def to_mef(data, channels, sc_list, sc_channels = None): sc(data_t._range[chi][1])] return data_t + + +def to_compensated(data, channels, a0, A, comp_channels=None): + """ + Transform flow cytometry data using compensation coefficients. + + This function accepts an autofluorescence vector `a0` and a + bleedthrough matrix `A` as compensation coefficients, with rows and + columns corresponding to channels specified in `comp_channels`. + `to_compensated` automatically checks whether compensation can be + performed for every channel specified in `channels`, and throws an + error otherwise. + + This function is intended to be reduced to the following signature:: + + to_compensated_reduced(data, channels) + + by using ``functools.partial`` once compensation coefficients and + `comp_channels` are available. + + Parameters + ---------- + data : FCSData or numpy array + NxD flow cytometry data where N is the number of events and D is + the number of parameters (aka channels). + channels : int, str, list of int, list of str + Channels on which to perform the transformation. If `channels` is + None, perform transformation in all channels specified on + `comp_channels`. + a0 : array + Autofluorescence vector, with a length equal to the length of + `comp_channels`. + A : 2D array + Bleedthrough matrix, a square matrix with a size equal to the + length of `comp_channels`. + comp_channels : list of int or list of str + List of channels on which compensation can be applied. Each element + corresponds to each row and column in `a0` and `A`. + + Returns + ------- + FCSData or numpy array + NxD transformed flow cytometry data. + + Raises + ------ + ValueError + If any channel specified in `channels` is not in `comp_channels`. + + """ + # Default comp_channels + if comp_channels is None: + if data.ndim == 1: + comp_channels = range(data.shape[0]) + else: + comp_channels = range(data.shape[1]) + # Convert comp_channels to indices + if hasattr(data, '_name_to_index'): + comp_channels = data._name_to_index(comp_channels) + + # Default channels + if channels is None: + channels = comp_channels + # Convert channels to iterable + if not (hasattr(channels, '__iter__') \ + and not isinstance(channels, six.string_types)): + channels = [channels] + # Convert channels to index + if hasattr(data, '_name_to_index'): + channels_ind = data._name_to_index(channels) + else: + channels_ind = channels + # Check if every channel is in comp_channels + for chi, chs in zip(channels_ind, channels): + if chi not in comp_channels: + raise ValueError( + "no compensation coefficients for channel {}".format(chs)) + + # Check appropriate dimensions of a0 and A + if a0.shape != (len(comp_channels),): + raise ValueError('length of a0 should be equal the number of elements' + ' in comp_channels') + if A.shape != (len(comp_channels), len(comp_channels)): + raise ValueError('A should be a square matrix with size equal to the' + ' number of elements in comp_channels') + + # Copy data array + data_t = data.copy().astype(np.float64) + + # Apply compensation to data + # Compensation will be applied to all channels in `comp_channels`, but + # only data in `channels` will be copied to the output array + comp_data = np.linalg.solve(A, (data_t[:, comp_channels] - a0).T).T + comp_channels_map = [comp_channels.index(chi) for chi in channels_ind] + data_t[:, channels_ind] = comp_data[:, comp_channels_map] + + # Apply compensation to range + if hasattr(data_t, '_range'): + range_low = np.array([data_t._range[chi][0] for chi in comp_channels]) + range_high = np.array([data_t._range[chi][1] for chi in comp_channels]) + + range_low_comp = np.linalg.solve(A, range_low - a0) + range_high_comp = np.linalg.solve(A, range_high - a0) + + for chi, chmi in zip(channels_ind, comp_channels_map): + data_t._range[chi] = [range_low_comp[chmi], range_high_comp[chmi]] + + return data_t diff --git a/doc/_static/img/python_tutorial/python_tutorial_compensate_1.png b/doc/_static/img/python_tutorial/python_tutorial_compensate_1.png new file mode 100644 index 0000000..814912b Binary files /dev/null and b/doc/_static/img/python_tutorial/python_tutorial_compensate_1.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_compensate_2.png b/doc/_static/img/python_tutorial/python_tutorial_compensate_2.png new file mode 100644 index 0000000..c304f0e Binary files /dev/null and b/doc/_static/img/python_tutorial/python_tutorial_compensate_2.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_compensate_3.png b/doc/_static/img/python_tutorial/python_tutorial_compensate_3.png new file mode 100644 index 0000000..bb536f8 Binary files /dev/null and b/doc/_static/img/python_tutorial/python_tutorial_compensate_3.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_compensate_4.png b/doc/_static/img/python_tutorial/python_tutorial_compensate_4.png new file mode 100644 index 0000000..f49ce5a Binary files /dev/null and b/doc/_static/img/python_tutorial/python_tutorial_compensate_4.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_gate_density_1.png b/doc/_static/img/python_tutorial/python_tutorial_gate_density_1.png index 46c4589..a6bc6c1 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_gate_density_1.png and b/doc/_static/img/python_tutorial/python_tutorial_gate_density_1.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_gate_density_2.png b/doc/_static/img/python_tutorial/python_tutorial_gate_density_2.png index ba17198..216d7bc 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_gate_density_2.png and b/doc/_static/img/python_tutorial/python_tutorial_gate_density_2.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_gate_ellipse_1.png b/doc/_static/img/python_tutorial/python_tutorial_gate_ellipse_1.png index f3ddd8e..d07679a 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_gate_ellipse_1.png and b/doc/_static/img/python_tutorial/python_tutorial_gate_ellipse_1.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_gate_high_low_1.png b/doc/_static/img/python_tutorial/python_tutorial_gate_high_low_1.png index 28cf2d4..5d19ea1 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_gate_high_low_1.png and b/doc/_static/img/python_tutorial/python_tutorial_gate_high_low_1.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_gate_high_low_2.png b/doc/_static/img/python_tutorial/python_tutorial_gate_high_low_2.png index 7b5bb11..3a7c0b0 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_gate_high_low_2.png and b/doc/_static/img/python_tutorial/python_tutorial_gate_high_low_2.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_mef_1.png b/doc/_static/img/python_tutorial/python_tutorial_mef_1.png index b225b61..cb96de3 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_mef_1.png and b/doc/_static/img/python_tutorial/python_tutorial_mef_1.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_mef_2.png b/doc/_static/img/python_tutorial/python_tutorial_mef_2.png index f274bfb..bf43b84 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_mef_2.png and b/doc/_static/img/python_tutorial/python_tutorial_mef_2.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_mef_3.png b/doc/_static/img/python_tutorial/python_tutorial_mef_3.png index f80aa66..3ad74f2 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_mef_3.png and b/doc/_static/img/python_tutorial/python_tutorial_mef_3.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_mef_5.png b/doc/_static/img/python_tutorial/python_tutorial_mef_5.png index 6371390..458ea57 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_mef_5.png and b/doc/_static/img/python_tutorial/python_tutorial_mef_5.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_mef_6.png b/doc/_static/img/python_tutorial/python_tutorial_mef_6.png index 8c9730d..40d2923 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_mef_6.png and b/doc/_static/img/python_tutorial/python_tutorial_mef_6.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_plot_density_2d_1.png b/doc/_static/img/python_tutorial/python_tutorial_plot_density_2d_1.png index 3c7fad6..1ac6993 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_plot_density_2d_1.png and b/doc/_static/img/python_tutorial/python_tutorial_plot_density_2d_1.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_plot_density_2d_2.png b/doc/_static/img/python_tutorial/python_tutorial_plot_density_2d_2.png index 1807c39..ca21ae6 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_plot_density_2d_2.png and b/doc/_static/img/python_tutorial/python_tutorial_plot_density_2d_2.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_plot_density_and_hist_1.png b/doc/_static/img/python_tutorial/python_tutorial_plot_density_and_hist_1.png index e8b5886..1316753 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_plot_density_and_hist_1.png and b/doc/_static/img/python_tutorial/python_tutorial_plot_density_and_hist_1.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_1.png b/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_1.png index 4d5993c..d229862 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_1.png and b/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_1.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_2.png b/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_2.png index 39f8a46..832f79e 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_2.png and b/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_2.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_3.png b/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_3.png index e3dbba0..f6cc5fe 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_3.png and b/doc/_static/img/python_tutorial/python_tutorial_plot_hist1d_3.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_plot_violin_1.png b/doc/_static/img/python_tutorial/python_tutorial_plot_violin_1.png index a09622d..ed7bf0c 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_plot_violin_1.png and b/doc/_static/img/python_tutorial/python_tutorial_plot_violin_1.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_plot_violin_2.png b/doc/_static/img/python_tutorial/python_tutorial_plot_violin_2.png index c9a655a..ee1848d 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_plot_violin_2.png and b/doc/_static/img/python_tutorial/python_tutorial_plot_violin_2.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_transform_1.png b/doc/_static/img/python_tutorial/python_tutorial_transform_1.png index 8ef8d8d..df36f36 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_transform_1.png and b/doc/_static/img/python_tutorial/python_tutorial_transform_1.png differ diff --git a/doc/_static/img/python_tutorial/python_tutorial_transform_2.png b/doc/_static/img/python_tutorial/python_tutorial_transform_2.png index b011ecc..9810386 100644 Binary files a/doc/_static/img/python_tutorial/python_tutorial_transform_2.png and b/doc/_static/img/python_tutorial/python_tutorial_transform_2.png differ diff --git a/doc/python_tutorial/compensate.rst b/doc/python_tutorial/compensate.rst new file mode 100644 index 0000000..69c8728 --- /dev/null +++ b/doc/python_tutorial/compensate.rst @@ -0,0 +1,232 @@ +Compensating Flow Cytometry Data +================================ + +This tutorial focuses on how to perform multi-color compensation on flow cytometry data using ``FlowCal``, particularly by using the module :mod:`FlowCal.compensate`. + +The dataset in ``examples`` includes nine cell samples from a strain containing both sfGFP and mCherry genes, exposed to nine different inducer (aTc) levels that result in changes of sfGFP and mCherry expression. Ideally, we would be able to observe sfGFP in fluorescence channel FL1 and mCherry in FL2 without any crosstalk or interference. In reality, the signal observed in each channel is the sum of the following three components: + +* The channel's "main" fluorophore (in this case sfGFP for FL1 and mCherry for FL2). +* "Bleedthrough" from all other fluorophores (sfGFP for FL2 and mCherry for FL1, although the latter, in this particular case, is close to zero). +* Cell autofluorescence. + +Here, we show how to perform multi-color compensation to extract the true fluorophore signal from each channel. To allow this, our dataset also includes the following controls (folder ``controls``): + +* One no-fluorophore control (NFC) taken from a strain contaning no fluorophores. +* Two single-fluorophore controls (SFCs) taken from strains containing a single fluorophore each - in this case, sfGFP and mCherry. + +Multi-color compensation requires all fluorescence data to be in units proportional to the fluorophore's emission intensity. In most cytometers, fluorescence in a.u. satisfies this requirement. However, older instruments can give a signal that is slightly nonlinear with respect to intensity. MEF calibration can compensate for this, and thus it is recommended to run compensation after calibration. In addition, calibration allows samples and controls taken at different cytometer settings to be used for compensation. Please refer to the :doc:`calibration tutorial ` for more information. In order to keep this tutorial focused, compensation here will be completely performed in a.u. + +Loading data and performing gating and calibration +-------------------------------------------------- + +To start, navigate to the ``examples`` directory included with FlowCal, and open a ``python`` session therein. Then, import ``numpy``, ``matplotlib``, and ``FlowCal``. + +>>> import numpy as np +>>> import matplotlib.pyplot as plt +>>> import FlowCal + +First, we load the experimental sample files and perform density gating: + +>>> samples_filenames = ['FCFiles/sample029.fcs', +... 'FCFiles/sample030.fcs', +... 'FCFiles/sample031.fcs', +... 'FCFiles/sample032.fcs', +... 'FCFiles/sample033.fcs', +... 'FCFiles/sample034.fcs', +... 'FCFiles/sample035.fcs', +... 'FCFiles/sample036.fcs', +... 'FCFiles/sample037.fcs'] +>>> # The list ``samples`` will store processed, transformed data of cell samples +>>> samples = [] +>>> # Iterate over cell sample filenames +>>> for sample_id, sample_filename in enumerate(samples_filenames): +... # Load file +... sample = FlowCal.io.FCSData(sample_filename) +... # Transform data to RFI +... sample = FlowCal.transform.to_rfi(sample) +... # Apply density gating +... sample = FlowCal.gate.density2d( +... data=sample, +... channels=['FSC','SSC'], +... gate_fraction=0.85) +... # Save +... samples.append(sample) + +Next, we do the same for the control samples: + +>>> # File names +>>> nfc_sample_filename = 'FCFiles/controls/nfc/sample004.fcs' +>>> sfc1_sample_filename = 'FCFiles/controls/sfc1/sample007.fcs' +>>> sfc2_sample_filename = 'FCFiles/controls/sfc2/sample019.fcs' +>>> # Load files +>>> nfc_sample = FlowCal.io.FCSData(nfc_sample_filename) +>>> sfc1_sample = FlowCal.io.FCSData(sfc1_sample_filename) +>>> sfc2_sample = FlowCal.io.FCSData(sfc2_sample_filename) +>>> # Transform to RFI +>>> nfc_sample = FlowCal.transform.to_rfi(nfc_sample) +>>> sfc1_sample = FlowCal.transform.to_rfi(sfc1_sample) +>>> sfc2_sample = FlowCal.transform.to_rfi(sfc2_sample) +>>> # Perform density gating +>>> nfc_sample = FlowCal.gate.density2d( +... data=nfc_sample, +... channels=['FSC','SSC'], +... gate_fraction=0.85) +>>> sfc1_sample = FlowCal.gate.density2d( +... data=sfc1_sample, +... channels=['FSC','SSC'], +... gate_fraction=0.85) +>>> sfc2_sample = FlowCal.gate.density2d( +... data=sfc2_sample, +... channels=['FSC','SSC'], +... gate_fraction=0.85) + +Observing bleedthrough +---------------------- + +Let's now look at the FL1 and FL2 fluorescence of the SFCs. As a reminder, each control should contain only one fluorophore (sfGFP or mCherry). Ideally, the sfGFP SFC should only produce signal in FL1, and the mCherry SFC should only produce FL2 signal. + +>>> # Obtain mean autofluorescence in FL1 and FL2 to plot in all panels +>>> autofl = FlowCal.stats.mean(nfc_sample, channels=['FL1', 'FL2']) +>>> # Plot controls +>>> samples_to_plot = [ +... nfc_sample, +... sfc1_sample, +... sfc2_sample, +... ] +>>> samples_titles = [ +... "No-fluroescence control", +... "SFC, sfGFP", +... "SFC, mCherry", +... ] +>>> plt.figure(figsize=(9, 3)) +>>> for plot_id, (sample_to_plot, sample_title) in \ +... enumerate(zip(samples_to_plot, samples_titles)): +... plt.subplot(1, 3, 1 + plot_id) +... # Density plot of sample +... FlowCal.plot.density2d( +... sample_to_plot, +... channels=['FL1', 'FL2'], +... mode='scatter') +... # Plot autofluorescence lines +... plt.axvline(autofl[0], color='gray') +... plt.axhline(autofl[1], color='gray') +... # Set the axes identically accross all samples +... plt.gca().set_xscale('logicle', T=2e4) +... plt.gca().set_yscale('logicle', T=2e4) +... plt.title(sample_title) +>>> plt.tight_layout() +>>> plt.show() + +.. image:: /_static/img/python_tutorial/python_tutorial_compensate_1.png + +While the mCherry SFC produces signal above autofluorescence in FL2 only (right plot), the sfGFP SFC results in signals in both channels (middle plot). While the sfGFP-induced FL2 signal is small in this case, it can make it difficult to resolve mCherry signals that are small to begin with. To see this, let's analyze the experimental samples: + +>>> # aTc concentration of each cell sample, in ng/mL. +>>> atc = np.array([0, 0.5, 1, 1.5, 2, 3, 4, 7.5, 20]) +>>> # Plot violins of experimental samples as a function of aTc +>>> # Plot the NFC to indicate the minimum possible fluorescence +>>> plt.figure(figsize=(8, 3.5)) +>>> plt.subplot(1, 2, 1) +>>> FlowCal.plot.violin_dose_response( +... data=samples, +... channel='FL1', +... positions=atc, +... min_data=nfc_sample, +... xlabel='aTc Concentration (ng/mL)', +... xscale='log', +... yscale='log', +... ylim=(1e0,1e3), +... violin_width=0.12, +... violin_kwargs={'facecolor': 'tab:green', +... 'edgecolor':'black'}, +... ) +>>> plt.ylabel('FL1 Fluorescence (MEFL)') +>>> plt.subplot(1, 2, 2) +>>> FlowCal.plot.violin_dose_response( +... data=samples, +... channel='FL2', +... positions=atc, +... min_data=nfc_sample, +... xlabel='aTc Concentration (ng/mL)', +... xscale='log', +... yscale='log', +... ylim=(1e0,2e3), +... violin_width=0.12, +... violin_kwargs={'facecolor': 'tab:orange', +... 'edgecolor':'black'}, +... ) +>>> plt.ylabel('FL2 Fluorescence (MEPE)') +>>> plt.tight_layout() +>>> plt.show() + +.. image:: /_static/img/python_tutorial/python_tutorial_compensate_2.png + +As we can see here, at low inducer (aTc) levels, FL2 fluorescence (right plot) is small but non-zero. This may be the result of a phenomenon called "leakiness", where the output of a genetic system is not completely off in a situation where it should be. However, given that sfGFP fluorescence is high at the same inducer levels (FL1, left), it is hard to know whether the observed FL2 signal is due to leaky mCherry expression or bleedthrough from sfGFP. + +Eliminating bleedthrough via compensation +----------------------------------------- + +In ``FlowCal``, compensation is performed by creating a transformation function using ``compensate.get_transform_fxn()``, which in turn requires data from the control samples. The resulting transformation function can be used afterwards to compensate data from other samples. + +>>> # Create compensation function +>>> compensation_fxn = FlowCal.compensate.get_transform_fxn( +... nfc_sample, +... [sfc1_sample, sfc2_sample], +... ['FL1', 'FL2'], +... ) +>>> # Apply compensation to the samples and the NFC +>>> samples_compensated = [compensation_fxn(s, ['FL1', 'FL2']) for s in samples] +>>> nfc_sample_compensated = compensation_fxn(nfc_sample, ['FL1', 'FL2']) + +>>> # Plot violins with compensated data +>>> plt.figure(figsize=(8, 3.5)) +>>> plt.subplot(1, 2, 1) +>>> FlowCal.plot.violin_dose_response( +... data=samples_compensated, +... channel='FL1', +... positions=atc, +... min_data=nfc_sample_compensated, +... xlabel='aTc Concentration (ng/mL)', +... xscale='log', +... yscale='logicle', +... ylim=(-3e1, 1e3), +... violin_width=0.12, +... violin_kwargs={'facecolor': 'tab:green', +... 'edgecolor':'black'}, +... ) +>>> plt.ylabel('FL1 Fluorescence (MEFL)') +>>> plt.subplot(1, 2, 2) +>>> FlowCal.plot.violin_dose_response( +... data=samples_compensated, +... channel='FL2', +... positions=atc, +... min_data=nfc_sample_compensated, +... xlabel='aTc Concentration (ng/mL)', +... xscale='log', +... yscale='logicle', +... ylim=(-3e1, 2e3), +... violin_width=0.12, +... violin_kwargs={'facecolor': 'tab:orange', +... 'edgecolor':'black'}, +... ) +>>> plt.ylabel('FL2 Fluorescence (MEPE)') +>>> plt.tight_layout() +>>> plt.show() + +.. image:: /_static/img/python_tutorial/python_tutorial_compensate_3.png + +Here we can observe two changes. First, the NFC (black violin) is now centered around zero. This is an effect of removing the autofluorescence component, measured from the NFC itself, from the FL2 signal during compensation. Second, FL2 violins at low inducer levels are now centered around zero as well. Because both autofluorescence and bleedthrough from sfGFP were removed by the compensation process, the fact that the remaining FL2 signal is zero shows that the output of the genetic system driving mCherry is not leaky as we hypothesized above. + +A final note about compensation: most flow cytometry software packages perform compensation without taking into account autofluorescence subtraction. In fact, one can mimic this procedure in ``FlowCal`` by calling ``compensate.get_transform_fxn()`` without an NFC: + +>>> compensation_fxn = FlowCal.compensate.get_transform_fxn( +... None, +... [sfc1_sample, sfc2_sample], +... ['FL1', 'FL2'], +... ) + +Differences resulting from the usage of an NFC are negligible when sample fluorescence is much greater than autofluorescence. This may happen when the fluorescence signal is actually really large, or with modern instruments where an NFC histogram would be centered around zero (although in our experience this does not always happen perfectly). However, in cases where sample fluorescence is close to autofluorescence, ignoring the NFC can lead to nonsensical results where low fluorescence levels are brought down below autofluorescence. In fact, if we run this compensation method with our samples we obtain the following violins: + +.. image:: /_static/img/python_tutorial/python_tutorial_compensate_4.png + +We recommend using both NFCs and SFCs when possible, ideally acquired simultaneously with the experimental samples. \ No newline at end of file diff --git a/doc/python_tutorial/gate.rst b/doc/python_tutorial/gate.rst index c91b944..8b4326b 100644 --- a/doc/python_tutorial/gate.rst +++ b/doc/python_tutorial/gate.rst @@ -15,9 +15,9 @@ Also, import ``numpy`` and ``pyplot`` from ``matplotlib`` Removing Saturated Events ------------------------- -We'll start by loading the data from file ``sample006.fcs`` into an ``FCSData`` object called ``s``. Then, transform all channels into a.u. +We'll start by loading the data from file ``sample029.fcs`` into an ``FCSData`` object called ``s``. Then, transform all channels into a.u. ->>> s = FlowCal.io.FCSData('FCFiles/sample006.fcs') +>>> s = FlowCal.io.FCSData('FCFiles/sample029.fcs') >>> s = FlowCal.transform.to_rfi(s) In the :doc:`plotting tutorial ` we looked at a density plot of the forward scatter/side scatter (``FSC``/``SSC``) channels and identified several clusters of particles (events). This density plot is repeated below for convenience. @@ -56,7 +56,7 @@ Ellipse Gate >>> s_g3 = FlowCal.gate.ellipse(s_g1, ... channels=['FSC', 'SSC'], ... log=True, -... center=(2.3, 2.78), +... center=(2.2, 2.8), ... a=0.3, ... b=0.2, ... theta=30/180.*np.pi) diff --git a/doc/python_tutorial/index.rst b/doc/python_tutorial/index.rst index c6520ab..65482c1 100644 --- a/doc/python_tutorial/index.rst +++ b/doc/python_tutorial/index.rst @@ -13,4 +13,5 @@ FlowCal's Python API Tutorial plot.rst gate.rst mef.rst + compensate.rst excel_ui.rst diff --git a/doc/python_tutorial/mef.rst b/doc/python_tutorial/mef.rst index 593b5a5..35d0340 100644 --- a/doc/python_tutorial/mef.rst +++ b/doc/python_tutorial/mef.rst @@ -24,7 +24,7 @@ As mentioned in the :doc:`fundamentals` section, conv ... gate_fraction=0.3, ... full_output=True) >>> b_g = density_gate_output.gated_data ->>> c = density_gate_output.contour +>>> c = density_gate_output.contour >>> FlowCal.plot.density_and_hist(b, ... gated_data=b_g, ... gate_contour=c, @@ -62,18 +62,15 @@ The argument ``plot`` instructs :func:`FlowCal.mef.get_transform_fxn` to generat Let's now use ``to_mef`` to transform fluroescence data to MEF. >>> # Load sample ->>> s = FlowCal.io.FCSData('FCFiles/sample006.fcs') ->>> +>>> s = FlowCal.io.FCSData('FCFiles/sample029.fcs') >>> # Transform all channels to a.u., and then FL1 to MEF. >>> s = FlowCal.transform.to_rfi(s) >>> s = to_mef(s, channels='FL1') ->>> >>> # Gate >>> s_g = FlowCal.gate.high_low(s, channels=['FSC', 'SSC']) >>> s_g = FlowCal.gate.density2d(s_g, ... channels=['FSC', 'SSC'], ... gate_fraction=0.5) ->>> >>> # Plot histogram of transformed channel >>> FlowCal.plot.hist1d(s_g, channel='FL1') >>> plt.show() diff --git a/doc/python_tutorial/plot.rst b/doc/python_tutorial/plot.rst index ea2386d..f7b2d3e 100644 --- a/doc/python_tutorial/plot.rst +++ b/doc/python_tutorial/plot.rst @@ -15,9 +15,9 @@ Also, import ``numpy`` and ``pyplot`` from ``matplotlib`` Histograms ---------- -Let's load the data from file ``sample006.fcs`` into an ``FCSData`` object called ``s``, and tranform all channels to arbitrary units. +Let's load the data from file ``sample029.fcs`` into an ``FCSData`` object called ``s``, and tranform all channels to arbitrary units. ->>> s = FlowCal.io.FCSData('FCFiles/sample006.fcs') +>>> s = FlowCal.io.FCSData('FCFiles/sample029.fcs') >>> s = FlowCal.transform.to_rfi(s) One is often interested in the fluorescence distribution across a population of cells. This is represented in a histogram. Since ``FCSData`` is a numpy array, one could use the standard ``hist`` function included in matplotlib. Alternatively, ``FlowCal`` includes its own histogram function specifically tailored to work with ``FCSData`` objects. For example, one can plot the contents of the ``FL1`` channel with a single call to :func:`FlowCal.plot.hist1d`. @@ -38,7 +38,7 @@ By default, :func:`FlowCal.plot.hist1d` uses something called *logicle* scaling Finally, :func:`FlowCal.plot.hist1d` can plot several FCSData objects at the same time. Let's now load 3 FCSData objects, transform all channels to a.u., and plot the ``FL1`` channel of all three with transparency. ->>> filenames = ['FCFiles/sample{:03d}.fcs'.format(i + 9) for i in range(3)] +>>> filenames = ['FCFiles/sample{:03d}.fcs'.format(i + 33) for i in range(3)] >>> d = [FlowCal.io.FCSData(filename) for filename in filenames] >>> d = [FlowCal.transform.to_rfi(di) for di in d] >>> FlowCal.plot.hist1d(d, channel='FL1', alpha=0.7, bins=128) @@ -93,11 +93,11 @@ Violin Plots Histograms, as shown above, can be used to plot and compare data from multiple samples. However, they can easily get too crowded. A more compact way is to use a violin plot, wherein vertical, normalized, symmetrical histograms ("violins") are shown centered on corresponding x-axis values. We can do this with the :func:`FlowCal.plot.violin` function. ->>> filenames = ['FCFiles/sample{:03d}.fcs'.format(i+6) for i in range(10)] +>>> filenames = ['FCFiles/sample{:03d}.fcs'.format(i+29) for i in range(9)] >>> d = [FlowCal.io.FCSData(filename) for filename in filenames] >>> d = [FlowCal.transform.to_rfi(di) for di in d] ->>> dapg = np.array([0, 2.33, 4.36, 8.16, 15.3, 28.6, 53.5, 100, 187, 350]) ->>> FlowCal.plot.violin(data=d, channel='FL1', positions=dapg, xlabel='DAPG (uM)', xscale='log', ylim=(1e0,2e3)) +>>> atc = np.array([0, 0.5, 1, 1.5, 2, 3, 4, 7.5, 20]) +>>> FlowCal.plot.violin(data=d, channel='FL1', positions=atc, xlabel='aTc (ng/mL)', xscale='log', ylim=(1e0,2e3)) >>> plt.show() .. image:: /_static/img/python_tutorial/python_tutorial_plot_violin_1.png @@ -107,24 +107,23 @@ Note that the x axis has been plotted on a logarithmic scale using the ``xscale` "Dose response" or "transfer" functions are common in biology. These sometimes include minimum (negative) and maximum (positive) controls, and are often approximated by mathematical models. The :func:`FlowCal.plot.violin_dose_response` function can be used to plot a full dose response dataset, including min data, max data, and a mathematical model. Min and max data are illustrated to the left of the plot, and the mathematical model is correctly illustrated even when a position=0 violin is illustrated separately when ``xscale`` is ``log``. >>> # Function specifying mathematical model ->>> def dapg_sensor_model(dapg_concentration): ->>> mn = 20 ->>> mx = 250. ->>> K = 20. ->>> n = 3.57 ->>> if dapg_concentration <= 0: ->>> return mn ->>> else: ->>> return mn + ((mx-mn)/(1+((K/dapg_concentration)**n))) ->>> +>>> def atc_sensor_model(atc_concentration): +>>> mn = 16 +>>> mx = 170. +>>> K = 3. +>>> n = 5. +>>> if atc_concentration <= 0: +>>> return mx +>>> else: +>>> return mn + ((mx-mn)/(1+((atc_concentration/K)**n))) >>> # Plot >>> FlowCal.plot.violin_dose_response( >>> data=d, >>> channel='FL1', ->>> positions=dapg, +>>> positions=atc, >>> min_data=d[0], >>> max_data=d[-1], ->>> model_fxn=dapg_sensor_model, +>>> model_fxn=atc_sensor_model, >>> xscale='log', >>> yscale='log', >>> ylim=(1e0,2e3), @@ -132,7 +131,7 @@ Note that the x axis has been plotted on a logarithmic scale using the ``xscale` >>> 'linewidth':3, >>> 'zorder':-1, >>> 'solid_capstyle':'butt'}) ->>> plt.xlabel('DAPG Concentration ($\mu M$)') +>>> plt.xlabel('aTc Concentration (ng/mL)') >>> plt.ylabel('FL1 Fluorescence (a.u.)') >>> plt.show() diff --git a/doc/python_tutorial/read.rst b/doc/python_tutorial/read.rst index a485761..3acde7d 100644 --- a/doc/python_tutorial/read.rst +++ b/doc/python_tutorial/read.rst @@ -9,14 +9,14 @@ To start, navigate to the ``examples`` directory included with FlowCal, and open FCS files are standard files in which flow cytometry data is stored. Normally, one FCS file corresponds to one sample. -The object :class:`FlowCal.io.FCSData` allows a user to open an FCS file. The following instruction opens the file ``sample006.fcs`` from the ``FCFiles`` folder, loads the information into an ``FCSData`` object, and assigns it to a variable ``s``. +The object :class:`FlowCal.io.FCSData` allows a user to open an FCS file. The following instruction opens the file ``sample029.fcs`` from the ``FCFiles`` folder, loads the information into an ``FCSData`` object, and assigns it to a variable ``s``. ->>> s = FlowCal.io.FCSData('FCFiles/sample006.fcs') +>>> s = FlowCal.io.FCSData('FCFiles/sample029.fcs') An ``FCSData`` object is a 2D ``numpy`` array with a few additional features. The first dimension indexes the event number, and the second dimension indexes the flow cytometry channel (or "parameter", as called by the FCS standard). We can see the number of events and channels using the standard ``numpy``'s ``shape`` property: >>> print(s.shape) -(32224, 8) +(34440, 8) As with any ``numpy`` array, we can slice an ``FCSData`` object. For example, let's obtain the first 100 events. @@ -28,7 +28,7 @@ Note that the product of slicing an FCSData object is also an FCSData object. We >>> s_sub_ch = s[:, [3, 4, 5]] >>> print(s_sub_ch.shape) -(32224, 3) +(34440, 3) However, it is not immediately obvious what channels we are getting. Fortunately, the ``FCSData`` object contains some additional information about the acquisition settings. In particular, we can check the name of the channels with the ``channels`` property. diff --git a/doc/python_tutorial/transform.rst b/doc/python_tutorial/transform.rst index c42661e..c4e6d6d 100644 --- a/doc/python_tutorial/transform.rst +++ b/doc/python_tutorial/transform.rst @@ -10,9 +10,9 @@ To start, navigate to the ``examples`` directory included with FlowCal, and open Transforming to Arbitrary Fluorescence Units (a.u.) --------------------------------------------------- -Start by loading file ``sample006.fcs`` into an ``FCSData`` object called ``s``. +Start by loading file ``sample029.fcs`` into an ``FCSData`` object called ``s``. ->>> s = FlowCal.io.FCSData('FCFiles/sample006.fcs') +>>> s = FlowCal.io.FCSData('FCFiles/sample029.fcs') Let's now visualize the contents of the ``FL1`` channel. We will explore ``FlowCal``'s plotting functions in the :doc:`plotting tutorial `, but for now let's just use ``matplotlib``'s ``hist`` function. diff --git a/doc/reference/FlowCal.compensate.rst b/doc/reference/FlowCal.compensate.rst new file mode 100644 index 0000000..6e42b88 --- /dev/null +++ b/doc/reference/FlowCal.compensate.rst @@ -0,0 +1,7 @@ +FlowCal.compensate module +======================== + +.. automodule:: FlowCal.compensate + :members: + :undoc-members: + :show-inheritance: diff --git a/doc/reference/modules.rst b/doc/reference/modules.rst index 9395544..c27421d 100644 --- a/doc/reference/modules.rst +++ b/doc/reference/modules.rst @@ -3,6 +3,7 @@ FlowCal (Python API) Reference .. toctree:: + FlowCal.compensate FlowCal.excel_ui FlowCal.gate FlowCal.io diff --git a/examples/FCFiles/controls/nfc/sample003.fcs b/examples/FCFiles/controls/nfc/sample003.fcs new file mode 100644 index 0000000..3569191 Binary files /dev/null and b/examples/FCFiles/controls/nfc/sample003.fcs differ diff --git a/examples/FCFiles/controls/nfc/sample004.fcs b/examples/FCFiles/controls/nfc/sample004.fcs new file mode 100644 index 0000000..f6992c7 Binary files /dev/null and b/examples/FCFiles/controls/nfc/sample004.fcs differ diff --git a/examples/FCFiles/controls/sfc1/sample003.fcs b/examples/FCFiles/controls/sfc1/sample003.fcs new file mode 100644 index 0000000..f5b3566 Binary files /dev/null and b/examples/FCFiles/controls/sfc1/sample003.fcs differ diff --git a/examples/FCFiles/controls/sfc1/sample007.fcs b/examples/FCFiles/controls/sfc1/sample007.fcs new file mode 100644 index 0000000..85674ed Binary files /dev/null and b/examples/FCFiles/controls/sfc1/sample007.fcs differ diff --git a/examples/FCFiles/controls/sfc2/sample001.fcs b/examples/FCFiles/controls/sfc2/sample001.fcs new file mode 100644 index 0000000..bf22165 Binary files /dev/null and b/examples/FCFiles/controls/sfc2/sample001.fcs differ diff --git a/examples/FCFiles/controls/sfc2/sample019.fcs b/examples/FCFiles/controls/sfc2/sample019.fcs new file mode 100644 index 0000000..39f73c7 Binary files /dev/null and b/examples/FCFiles/controls/sfc2/sample019.fcs differ diff --git a/examples/FCFiles/max/sample002.fcs b/examples/FCFiles/max/sample002.fcs deleted file mode 100644 index 04501b5..0000000 Binary files a/examples/FCFiles/max/sample002.fcs and /dev/null differ diff --git a/examples/FCFiles/max/sample008.fcs b/examples/FCFiles/max/sample008.fcs deleted file mode 100644 index 5673167..0000000 Binary files a/examples/FCFiles/max/sample008.fcs and /dev/null differ diff --git a/examples/FCFiles/min/sample001.fcs b/examples/FCFiles/min/sample001.fcs deleted file mode 100644 index 158bedc..0000000 Binary files a/examples/FCFiles/min/sample001.fcs and /dev/null differ diff --git a/examples/FCFiles/min/sample004.fcs b/examples/FCFiles/min/sample004.fcs deleted file mode 100644 index 0055fa2..0000000 Binary files a/examples/FCFiles/min/sample004.fcs and /dev/null differ diff --git a/examples/FCFiles/sample001.fcs b/examples/FCFiles/sample001.fcs index 18b3942..bf22165 100644 Binary files a/examples/FCFiles/sample001.fcs and b/examples/FCFiles/sample001.fcs differ diff --git a/examples/FCFiles/sample006.fcs b/examples/FCFiles/sample006.fcs deleted file mode 100644 index 83f3db2..0000000 Binary files a/examples/FCFiles/sample006.fcs and /dev/null differ diff --git a/examples/FCFiles/sample007.fcs b/examples/FCFiles/sample007.fcs deleted file mode 100644 index 85ca13e..0000000 Binary files a/examples/FCFiles/sample007.fcs and /dev/null differ diff --git a/examples/FCFiles/sample008.fcs b/examples/FCFiles/sample008.fcs deleted file mode 100644 index a73f322..0000000 Binary files a/examples/FCFiles/sample008.fcs and /dev/null differ diff --git a/examples/FCFiles/sample009.fcs b/examples/FCFiles/sample009.fcs deleted file mode 100644 index 19a49dd..0000000 Binary files a/examples/FCFiles/sample009.fcs and /dev/null differ diff --git a/examples/FCFiles/sample010.fcs b/examples/FCFiles/sample010.fcs deleted file mode 100644 index e7f74ca..0000000 Binary files a/examples/FCFiles/sample010.fcs and /dev/null differ diff --git a/examples/FCFiles/sample011.fcs b/examples/FCFiles/sample011.fcs deleted file mode 100644 index c68490c..0000000 Binary files a/examples/FCFiles/sample011.fcs and /dev/null differ diff --git a/examples/FCFiles/sample012.fcs b/examples/FCFiles/sample012.fcs deleted file mode 100644 index f7ebb82..0000000 Binary files a/examples/FCFiles/sample012.fcs and /dev/null differ diff --git a/examples/FCFiles/sample013.fcs b/examples/FCFiles/sample013.fcs deleted file mode 100644 index d959490..0000000 Binary files a/examples/FCFiles/sample013.fcs and /dev/null differ diff --git a/examples/FCFiles/sample014.fcs b/examples/FCFiles/sample014.fcs deleted file mode 100644 index c49c348..0000000 Binary files a/examples/FCFiles/sample014.fcs and /dev/null differ diff --git a/examples/FCFiles/sample015.fcs b/examples/FCFiles/sample015.fcs deleted file mode 100644 index 7578c33..0000000 Binary files a/examples/FCFiles/sample015.fcs and /dev/null differ diff --git a/examples/FCFiles/sample029.fcs b/examples/FCFiles/sample029.fcs new file mode 100644 index 0000000..3f60515 Binary files /dev/null and b/examples/FCFiles/sample029.fcs differ diff --git a/examples/FCFiles/sample030.fcs b/examples/FCFiles/sample030.fcs new file mode 100644 index 0000000..cd3df6f Binary files /dev/null and b/examples/FCFiles/sample030.fcs differ diff --git a/examples/FCFiles/sample031.fcs b/examples/FCFiles/sample031.fcs new file mode 100644 index 0000000..156b551 Binary files /dev/null and b/examples/FCFiles/sample031.fcs differ diff --git a/examples/FCFiles/sample032.fcs b/examples/FCFiles/sample032.fcs new file mode 100644 index 0000000..2a85c5a Binary files /dev/null and b/examples/FCFiles/sample032.fcs differ diff --git a/examples/FCFiles/sample033.fcs b/examples/FCFiles/sample033.fcs new file mode 100644 index 0000000..56e1de9 Binary files /dev/null and b/examples/FCFiles/sample033.fcs differ diff --git a/examples/FCFiles/sample034.fcs b/examples/FCFiles/sample034.fcs new file mode 100644 index 0000000..aa0dab7 Binary files /dev/null and b/examples/FCFiles/sample034.fcs differ diff --git a/examples/FCFiles/sample035.fcs b/examples/FCFiles/sample035.fcs new file mode 100644 index 0000000..72ad08d Binary files /dev/null and b/examples/FCFiles/sample035.fcs differ diff --git a/examples/FCFiles/sample036.fcs b/examples/FCFiles/sample036.fcs new file mode 100644 index 0000000..4204af9 Binary files /dev/null and b/examples/FCFiles/sample036.fcs differ diff --git a/examples/FCFiles/sample037.fcs b/examples/FCFiles/sample037.fcs new file mode 100644 index 0000000..8b551fb Binary files /dev/null and b/examples/FCFiles/sample037.fcs differ diff --git a/examples/analyze_excel_ui.py b/examples/analyze_excel_ui.py index 29ca2e3..4aeb155 100644 --- a/examples/analyze_excel_ui.py +++ b/examples/analyze_excel_ui.py @@ -7,7 +7,9 @@ Excel file. The exact operations performed are identical to when normally using the Excel UI. However, instead of generating an output Excel file, an OrderedDict of objects representing gated and transformed flow cytometry -samples is obtained. +samples is obtained. In addition, multi-color compensation is performed on +all samples using data from no-fluorophore and single-fluorophore control +samples (NFC and SFCs). Part two exemplifies how to use the processed cell sample data with FlowCal's plotting and statistics modules to produce interesting plots. @@ -16,6 +18,7 @@ consult readme.txt. """ +import six import numpy as np import pandas as pd import matplotlib as mpl @@ -81,6 +84,23 @@ plot=True, plot_dir='plot_samples') + # Perform multi-color compensation + # ``FlowCal.compensate.get_transform_fxn()`` generates a transformation + # function that performs multi-color compensation on a specified set of + # channels in order to remove fluorophore bleedthrough. + # This function requires data from single-fluorophore controls (SFCs), one + # per channel to compensate, each from cells containing only one + # fluorophore. This function can optionally use data from a no-fluorophore + # control (NFC). + compensation_fxn = FlowCal.compensate.get_transform_fxn( + nfc_sample=samples['NFC'], + sfc_samples=[samples['SFC1'], samples['SFC2']], + comp_channels=['FL1', 'FL2'], + ) + # Compensate all samples + samples_compensated = {s_id: compensation_fxn(s, ['FL1', 'FL2']) + for s_id, s in six.iteritems(samples)} + ### # Part 2: Examples on how to use processed cell sample data ### @@ -88,44 +108,52 @@ # Each entry in the Excel table has a corresponding ID, which can be used # to reference information associated with or the sample loaded from a # row in the Excel file. Collect the IDs of the non-control samples (i.e., - # 'S0001', ..., 'S0010'). - sample_ids = ['S00{:02}'.format(n) for n in range(1,10+1)] + # 'S001', ..., 'S009'). + sample_ids = ['S0{:02}'.format(n) for n in range(1, 9+1)] - # We will read DAPG concentrations from the Excel file. ``samples_table`` + # We will read aTc concentrations from the Excel file. ``samples_table`` # contains all the data from sheet "Samples", including data not directly # used by ``FlowCal.excel_ui.process_samples_table()``. ``samples_table`` # is a pandas dataframe, with each column having the same name as the # corresponding header in the Excel file. - dapg = samples_table.loc[sample_ids,'DAPG (uM)'] + atc = samples_table.loc[sample_ids, 'aTc (ng/mL)'] + + # We now show how to generate plots using the processed flow cytometry + # data we just obtained. + print("\nGenerating plots...") # Plot 1: Histogram of all samples # - # Here, we plot the fluorescence histograms of all ten samples in the same + # Here, we plot the fluorescence histograms of all nine samples in the same # figure using ``FlowCal.plot.hist1d``. Note how this function can be used # in the context of accessory matplotlib functions to modify the axes # limits and labels and to add a legend, among other things. - # Color each histogram according to its DAPG concentration. Linearize the - # color transitions using a logarithmic normalization to match the - # logarithmic spacing of the DAPG concentrations. (Concentrations are also - # augmented slightly to move the 0.0 concentration into the log - # normalization range.) - cmap = mpl.cm.get_cmap('gray_r') - norm = mpl.colors.LogNorm(vmin=1e0, vmax=3500.) - colors = [cmap(norm(dapg_i+4.)) for dapg_i in dapg] - - plt.figure(figsize=(6,3.5)) + plt.figure(figsize=(6, 5.5)) + plt.subplot(2, 1, 1) FlowCal.plot.hist1d([samples[s_id] for s_id in sample_ids], channel='FL1', histtype='step', - bins=128, - edgecolor=colors) + bins=128) plt.ylim((0,2500)) plt.xlim((0,5e4)) plt.xlabel('FL1 (Molecules of Equivalent Fluorescein, MEFL)') - plt.legend(['{:.1f} $\mu M$ DAPG'.format(i) for i in dapg], + plt.legend(['{:.1f} ng/mL aTc'.format(i) for i in atc], loc='upper left', fontsize='small') + + plt.subplot(2, 1, 2) + FlowCal.plot.hist1d([samples[s_id] for s_id in sample_ids], + channel='FL2', + histtype='step', + bins=128) + plt.ylim((0,2500)) + plt.xlim((0,5e4)) + plt.xlabel('FL2 (Molecules of Equivalent PE, MEPE)') + plt.legend(['{:.1f} ng/mL aTc'.format(i) for i in atc], + loc='upper left', + fontsize='small') + plt.tight_layout() plt.savefig('histograms.png', dpi=200) plt.close() @@ -135,40 +163,67 @@ # Here, we illustrate how to obtain statistics from the fluorescence of # each sample and how to use them in a plot. The stats module contains # functions to calculate different statistics such as mean, median, and - # standard deviation. In this example, we calculate the mean from channel - # FL1 of each sample and plot them against the corresponding DAPG + # standard deviation. In this example, we calculate the mean from channels + # FL1 and FL2 of each sample and plot them against the corresponding aTc # concentrations. - samples_fluorescence = [FlowCal.stats.mean(samples[s_id], channels='FL1') - for s_id in sample_ids] - min_fluorescence = FlowCal.stats.mean(samples['min'], channels='FL1') - max_fluorescence = FlowCal.stats.mean(samples['max'], channels='FL1') - - dapg_color = '#ffc400' # common color used for DAPG-related plots + samples_fl1 = [FlowCal.stats.mean(samples[s_id], channels='FL1') + for s_id in sample_ids] + samples_fl2 = [FlowCal.stats.mean(samples[s_id], channels='FL2') + for s_id in sample_ids] + # No fluorescence control (NFC) will give the minimum fluorescence level in + # both channels. Single fluorescence controls (SFCs) containing sfGFP or + # mCherry only will give the maximum levels in channels FL1 and FL2. + min_fl1 = FlowCal.stats.mean(samples['NFC'], channels='FL1') + max_fl1 = FlowCal.stats.mean(samples['SFC1'], channels='FL1') + min_fl2 = FlowCal.stats.mean(samples['NFC'], channels='FL2') + max_fl2 = FlowCal.stats.mean(samples['SFC2'], channels='FL2') - plt.figure(figsize=(3,3)) - plt.plot(dapg, - samples_fluorescence, + plt.figure(figsize=(6, 3)) + + plt.subplot(1, 2, 1) + plt.plot(atc, + samples_fl1, marker='o', - color=dapg_color) - - # Illustrate min and max bounds - plt.axhline(min_fluorescence, + color='tab:green') + plt.axhline(min_fl1, color='gray', linestyle='--', zorder=-1) - plt.text(s='Min', x=2e2, y=1.6e2, ha='left', va='bottom', color='gray') - plt.axhline(max_fluorescence, + plt.text(s='Min', x=3e1, y=2e2, ha='left', va='bottom', color='gray') + plt.axhline(max_fl1, color='gray', linestyle='--', zorder=-1) - plt.text(s='Max', x=-0.7, y=5.2e3, ha='left', va='top', color='gray') - + plt.text(s='Max', x=-0.8, y=5.2e3, ha='left', va='top', color='gray') plt.yscale('log') plt.ylim((5e1,1e4)) plt.xscale('symlog') - plt.xlim((-1e0, 1e3)) - plt.xlabel('DAPG Concentration ($\mu M$)') + plt.xlim((-1e0, 1e2)) + plt.xlabel('aTc Concentration (ng/mL)') plt.ylabel('FL1 Fluorescence (MEFL)') + + plt.subplot(1, 2, 2) + plt.plot(atc, + samples_fl2, + marker='o', + color='tab:orange') + plt.axhline(min_fl2, + color='gray', + linestyle='--', + zorder=-1) + plt.text(s='Min', x=3e1, y=3.5e1, ha='left', va='bottom', color='gray') + plt.axhline(max_fl2, + color='gray', + linestyle='--', + zorder=-1) + plt.text(s='Max', x=-0.8, y=2.e3, ha='left', va='top', color='gray') + plt.yscale('log') + plt.ylim((1e1,5e3)) + plt.xscale('symlog') + plt.xlim((-1e0, 1e2)) + plt.xlabel('aTc Concentration (ng/mL)') + plt.ylabel('FL2 Fluorescence (MEPE)') + plt.tight_layout() plt.savefig('dose_response.png', dpi=200) plt.close() @@ -176,62 +231,114 @@ # Plot 3: Dose response violin plot # # Here, we use a violin plot to show the fluorescence of (almost) all - # cells as a function of DAPG. (The `upper_trim_fraction` and + # cells as a function of aTc. (The `upper_trim_fraction` and # `lower_trim_fraction` parameters eliminate the top and bottom 1% of # cells from each violin for aesthetic reasons. The summary statistic, # which is illustrated as a horizontal line atop each violin, is - # calculated before cells are removed, though.) We set `yscale` to 'log' + # calculated before cells are removed, though). We set `yscale` to 'log' # because the cytometer used to collect this data produces positive # integer data (as opposed to floating-point data, which can sometimes be # negative), so the added complexity of a logicle y-scale (which is the # default) is not necessary. + plt.figure(figsize=(8, 3.5)) - # FlowCal violin plots can also illustrate a mathematical model alongside - # the violins. To take advantage of this feature, we first recapitulate a - # model of the fluorescent protein (sfGFP) fluorescence produced by this - # DAPG sensor as a function of DAPG (from this study: - # https://doi.org/10.15252/msb.20209618). sfGFP fluorescence is cellular - # fluorescence minus autofluorescence. - def dapg_sensor_output(dapg_concentration): - mn = 86. - mx = 3147. - K = 20. - n = 3.57 - if dapg_concentration <= 0: - return mn - else: - return mn + ((mx-mn)/(1+((K/dapg_concentration)**n))) - - # To model cellular fluorescence, which we are plotting with this violin - # plot, we must add autofluorescence back to the sfGFP signal. For this - # model, autofluorescence is the mean fluorescence of an E. coli strain - # lacking sfGFP, which our min control is. - autofluorescence = FlowCal.stats.mean(samples['min'], channels='FL1') - def dapg_sensor_cellular_fluorescence(dapg_concentration): - return dapg_sensor_output(dapg_concentration) + autofluorescence - - plt.figure(figsize=(4,3.5)) + plt.subplot(1, 2, 1) FlowCal.plot.violin_dose_response( data=[samples[s_id] for s_id in sample_ids], channel='FL1', - positions=dapg, - min_data=samples['min'], - max_data=samples['max'], - model_fxn=dapg_sensor_cellular_fluorescence, - violin_kwargs={'facecolor':dapg_color, - 'edgecolor':'black'}, - violin_width_to_span_fraction=0.075, + positions=atc, + min_data=samples['NFC'], + max_data=samples['SFC1'], + xlabel='aTc Concentration (ng/mL)', xscale='log', yscale='log', - ylim=(1e1, 3e4), + ylim=(1e1,1e4), + violin_width=0.12, + violin_kwargs={'facecolor': 'tab:green', + 'edgecolor':'black'}, draw_model_kwargs={'color':'gray', 'linewidth':3, 'zorder':-1, - 'solid_capstyle':'butt'}) - plt.xlabel('DAPG Concentration ($\mu M$)') + 'solid_capstyle':'butt'}, + ) plt.ylabel('FL1 Fluorescence (MEFL)') + + plt.subplot(1, 2, 2) + FlowCal.plot.violin_dose_response( + data=[samples[s_id] for s_id in sample_ids], + channel='FL2', + positions=atc, + min_data=samples['NFC'], + max_data=samples['SFC2'], + xlabel='aTc Concentration (ng/mL)', + xscale='log', + yscale='log', + ylim=(1e0,1e4), + violin_width=0.12, + violin_kwargs={'facecolor': 'tab:orange', + 'edgecolor':'black'}, + draw_model_kwargs={'color':'gray', + 'linewidth':3, + 'zorder':-1, + 'solid_capstyle':'butt'}, + ) + plt.ylabel('FL2 Fluorescence (MEPE)') + plt.tight_layout() plt.savefig('dose_response_violin.png', dpi=200) plt.close() + # Plot 4: Dose response violin plot of compensated data + # + # Here, we redraw the previous violin plot but using compensated data. + # y axis will now be plotted in ``logicle`` scale since histograms will + # be centered around zero due to compensation. + plt.figure(figsize=(8, 3.5)) + + plt.subplot(1, 2, 1) + FlowCal.plot.violin_dose_response( + data=[samples_compensated[s_id] for s_id in sample_ids], + channel='FL1', + positions=atc, + min_data=samples_compensated['NFC'], + max_data=samples_compensated['SFC1'], + xlabel='aTc Concentration (ng/mL)', + xscale='log', + yscale='logicle', + ylim=(-3e2, 1e4), + violin_width=0.12, + violin_kwargs={'facecolor': 'tab:green', + 'edgecolor':'black'}, + draw_model_kwargs={'color':'gray', + 'linewidth':3, + 'zorder':-1, + 'solid_capstyle':'butt'}, + ) + plt.ylabel('FL1 Fluorescence (MEFL)') + + plt.subplot(1, 2, 2) + FlowCal.plot.violin_dose_response( + data=[samples_compensated[s_id] for s_id in sample_ids], + channel='FL2', + positions=atc, + min_data=samples_compensated['NFC'], + max_data=samples_compensated['SFC2'], + xlabel='aTc Concentration (ng/mL)', + xscale='log', + yscale='logicle', + ylim=(-1e2, 1e4), + violin_width=0.12, + violin_kwargs={'facecolor': 'tab:orange', + 'edgecolor':'black'}, + draw_model_kwargs={'color':'gray', + 'linewidth':3, + 'zorder':-1, + 'solid_capstyle':'butt'}, + ) + plt.ylabel('FL2 Fluorescence (MEPE)') + + plt.tight_layout() + plt.savefig('dose_response_violin_compensated.png', dpi=200) + plt.close() + print("\nDone.") diff --git a/examples/analyze_mef.py b/examples/analyze_mef.py index 17ae9be..70ee110 100644 --- a/examples/analyze_mef.py +++ b/examples/analyze_mef.py @@ -8,9 +8,11 @@ Molecules of Equivalent Fluorophore (MEF). At several points in this process, plots are generated that give insight into the relevant steps. -Part two processes data from twelve cell samples and uses the MEF +Part two processes data from nine cell samples and uses the MEF transformation functions from part one to convert fluorescence of these -samples to MEF. Plots are also generated in this stage. +samples to MEF. Plots are also generated in this stage. In addition, +multi-color compensation is performed on all samples using data from +no-fluorophore and single-fluorophore control samples (NFC and SFCs). Part three exemplifies how to use the processed cell sample data with FlowCal's plotting and statistics modules to produce interesting plots. @@ -30,39 +32,46 @@ # Definition of constants ### -# Name of the FCS files containing calibration beads data. The min and max -# controls were measured on separate days with their own beads samples, and -# one control was measured with a different cytometer gain setting, but we can -# still compare the min and max data to our samples because we are calibrating -# all fluorescence measurements to MEF units. +# Name of the FCS files containing calibration beads data. The no-fluorophore +# control (NFC) and single-fluorophore controls (SFC) were measured on separate +# days with their own beads samples, and one control was measured with a +# different cytometer gain setting. However, we can still compare these to our +# samples because we are calibrating all fluorescence measurements to MEF units. beads_filename = 'FCFiles/sample001.fcs' -min_beads_filename = 'FCFiles/min/sample001.fcs' -max_beads_filename = 'FCFiles/max/sample002.fcs' +nfc_beads_filename = 'FCFiles/controls/nfc/sample003.fcs' +sfc1_beads_filename = 'FCFiles/controls/sfc1/sample003.fcs' +sfc2_beads_filename = 'FCFiles/controls/sfc2/sample001.fcs' # Names of the FCS files containing data from cell samples -samples_filenames = ['FCFiles/sample006.fcs', - 'FCFiles/sample007.fcs', - 'FCFiles/sample008.fcs', - 'FCFiles/sample009.fcs', - 'FCFiles/sample010.fcs', - 'FCFiles/sample011.fcs', - 'FCFiles/sample012.fcs', - 'FCFiles/sample013.fcs', - 'FCFiles/sample014.fcs', - 'FCFiles/sample015.fcs'] -min_sample_filename = 'FCFiles/min/sample004.fcs' -max_sample_filename = 'FCFiles/max/sample008.fcs' +samples_filenames = ['FCFiles/sample029.fcs', + 'FCFiles/sample030.fcs', + 'FCFiles/sample031.fcs', + 'FCFiles/sample032.fcs', + 'FCFiles/sample033.fcs', + 'FCFiles/sample034.fcs', + 'FCFiles/sample035.fcs', + 'FCFiles/sample036.fcs', + 'FCFiles/sample037.fcs'] +nfc_sample_filename = 'FCFiles/controls/nfc/sample004.fcs' +sfc1_sample_filename = 'FCFiles/controls/sfc1/sample007.fcs' +sfc2_sample_filename = 'FCFiles/controls/sfc2/sample019.fcs' # Fluorescence values of each bead subpopulation, in MEF. # These values should be taken from the datasheet provided by the bead # manufacturer. We take Molecules of Equivalent Fluorescein (MEFL) to calibrate # the FL1 (GFP) channel. -mefl_values = np.array([0, 792, 2079, 6588, 16471, 47497, 137049, 271647]) -min_mefl_values = np.array([0, 771, 2106, 6262, 15183, 45292, 136258, 291042]) -max_mefl_values = np.array([0, 792, 2079, 6588, 16471, 47497, 137049, 271647]) +mefl_values = np.array([0, 792, 2079, 6588, 16471, 47497, 137049, None]) +nfc_mefl_values = np.array([0, 792, 2079, 6588, 16471, 47497, 137049, None]) +sfc1_mefl_values = np.array([0, 792, 2079, 6588, 16471, 47497, 137049, 271647]) +sfc2_mefl_values = np.array([0, 792, 2079, 6588, 16471, 47497, 137049, None]) -# DAPG concentration of each cell sample, in micromolar. -dapg = np.array([0, 2.33, 4.36, 8.16, 15.3, 28.6, 53.5, 100, 187, 350]) +mepe_values = np.array([0, 531, 1504, 4819, 12506, 36159, 109588, 250892]) +nfc_mepe_values = np.array([0, 531, 1504, 4819, 12506, 36159, 109588, 250892]) +sfc1_mepe_values = np.array([0, 531, 1504, 4819, 12506, 36159, 109588, 250892]) +sfc2_mepe_values = np.array([0, 531, 1504, 4819, 12506, 36159, 109588, 250892]) + +# aTc concentration of each cell sample, in ng/mL. +atc = np.array([0, 0.5, 1, 1.5, 2, 3, 4, 7.5, 20]) # Plots will be generated at various stages of analysis. The following are the # names of the folders in which we will store these plots. @@ -87,8 +96,12 @@ # cytometry data loaded from file ``filename``. print("Loading file \"{}\"...".format(beads_filename)) beads_sample = FlowCal.io.FCSData(beads_filename) - min_beads_sample = FlowCal.io.FCSData(min_beads_filename) - max_beads_sample = FlowCal.io.FCSData(max_beads_filename) + print("Loading file \"{}\"...".format(nfc_beads_filename)) + nfc_beads_sample = FlowCal.io.FCSData(nfc_beads_filename) + print("Loading file \"{}\"...".format(sfc1_beads_filename)) + sfc1_beads_sample = FlowCal.io.FCSData(sfc1_beads_filename) + print("Loading file \"{}\"...".format(sfc2_beads_filename)) + sfc2_beads_sample = FlowCal.io.FCSData(sfc2_beads_filename) # Data loaded from an FCS file is in "Channel Units", the raw numbers # reported from the instrument's detectors. The FCS file also contains @@ -97,8 +110,9 @@ # The function ``FlowCal.transform.to_rfi()`` performs this conversion. print("Performing data transformation...") beads_sample = FlowCal.transform.to_rfi(beads_sample) - min_beads_sample = FlowCal.transform.to_rfi(min_beads_sample) - max_beads_sample = FlowCal.transform.to_rfi(max_beads_sample) + nfc_beads_sample = FlowCal.transform.to_rfi(nfc_beads_sample) + sfc1_beads_sample = FlowCal.transform.to_rfi(sfc1_beads_sample) + sfc2_beads_sample = FlowCal.transform.to_rfi(sfc2_beads_sample) # Gating @@ -112,12 +126,15 @@ beads_sample_gated = FlowCal.gate.start_end(beads_sample, num_start=250, num_end=100) - min_beads_sample_gated = FlowCal.gate.start_end(min_beads_sample, - num_start=250, - num_end=100) - max_beads_sample_gated = FlowCal.gate.start_end(max_beads_sample, + nfc_beads_sample_gated = FlowCal.gate.start_end(nfc_beads_sample, num_start=250, num_end=100) + sfc1_beads_sample_gated = FlowCal.gate.start_end(sfc1_beads_sample, + num_start=250, + num_end=100) + sfc2_beads_sample_gated = FlowCal.gate.start_end(sfc2_beads_sample, + num_start=250, + num_end=100) # ``FlowCal.gate.high_low()`` removes events outside a range specified by # a ``low`` and a ``high`` value. If these are not specified (as shown @@ -137,10 +154,12 @@ # channels. beads_sample_gated = FlowCal.gate.high_low(beads_sample_gated, channels=['FSC','SSC']) - min_beads_sample_gated = FlowCal.gate.high_low(min_beads_sample_gated, - channels=['FSC','SSC']) - max_beads_sample_gated = FlowCal.gate.high_low(max_beads_sample_gated, + nfc_beads_sample_gated = FlowCal.gate.high_low(nfc_beads_sample_gated, channels=['FSC','SSC']) + sfc1_beads_sample_gated = FlowCal.gate.high_low(sfc1_beads_sample_gated, + channels=['FSC','SSC']) + sfc2_beads_sample_gated = FlowCal.gate.high_low(sfc2_beads_sample_gated, + channels=['FSC','SSC']) # ``FlowCal.gate.density2d()`` preserves only the densest population as # seen in a 2D density diagram of two channels. This helps remove particle @@ -162,23 +181,32 @@ beads_sample_gated = density_gate_output.gated_data gate_contour = density_gate_output.contour - min_density_gate_output = FlowCal.gate.density2d( - data=min_beads_sample_gated, + density_gate_output = FlowCal.gate.density2d( + data=nfc_beads_sample_gated, + channels=['FSC','SSC'], + gate_fraction=0.85, + sigma=5., + full_output=True) + nfc_beads_sample_gated = density_gate_output.gated_data + nfc_gate_contour = density_gate_output.contour + + density_gate_output = FlowCal.gate.density2d( + data=sfc1_beads_sample_gated, channels=['FSC','SSC'], gate_fraction=0.85, sigma=5., full_output=True) - min_beads_sample_gated = min_density_gate_output.gated_data - min_gate_contour = min_density_gate_output.contour + sfc1_beads_sample_gated = density_gate_output.gated_data + sfc1_gate_contour = density_gate_output.contour - max_density_gate_output = FlowCal.gate.density2d( - data=max_beads_sample_gated, + density_gate_output = FlowCal.gate.density2d( + data=sfc2_beads_sample_gated, channels=['FSC','SSC'], gate_fraction=0.85, sigma=5., full_output=True) - max_beads_sample_gated = max_density_gate_output.gated_data - max_gate_contour = max_density_gate_output.contour + sfc2_beads_sample_gated = density_gate_output.gated_data + sfc2_gate_contour = density_gate_output.contour # Plot forward/side scatter 2D density plot and 1D fluorescence histograms print("Plotting density plot and histogram...") @@ -200,16 +228,18 @@ # etc.) by just changing the extension. plot_filename = '{}/density_hist_{}.png'.format(beads_plot_dir, 'beads') - min_plot_filename = '{}/min_density_hist_{}.png'.format(beads_plot_dir, + nfc_plot_filename = '{}/nfc_density_hist_{}.png'.format(beads_plot_dir, + 'beads') + sfc1_plot_filename = '{}/sfc1_density_hist_{}.png'.format(beads_plot_dir, 'beads') - max_plot_filename = '{}/max_density_hist_{}.png'.format(beads_plot_dir, + sfc2_plot_filename = '{}/sfc2_density_hist_{}.png'.format(beads_plot_dir, 'beads') # Plot and save # The function ``FlowCal.plot.density_and_hist()`` plots a combined figure # with a 2D density plot at the top and an arbitrary number of 1D # histograms below. In this case, we will plot the forward/side scatter - # channels in the density plot and the fluorescence channels FL1 and FL3 + # channels in the density plot and the fluorescence channels FL1 and FL2 # below as two separate histograms. # Note that we are providing data both before (``beads_sample``) and after # (``beads_sample_gated``) gating. Each 1D histogram will display the @@ -222,67 +252,80 @@ beads_sample, beads_sample_gated, density_channels=['FSC', 'SSC'], - hist_channels=['FL1', 'FL3'], + hist_channels=['FL1', 'FL2'], gate_contour=gate_contour, density_params=density_params, savefig=plot_filename) FlowCal.plot.density_and_hist( - min_beads_sample, - min_beads_sample_gated, + nfc_beads_sample, + nfc_beads_sample_gated, + density_channels=['FSC', 'SSC'], + hist_channels=['FL1', 'FL2'], + gate_contour=nfc_gate_contour, + density_params=density_params, + savefig=nfc_plot_filename) + FlowCal.plot.density_and_hist( + sfc1_beads_sample, + sfc1_beads_sample_gated, density_channels=['FSC', 'SSC'], - hist_channels=['FL1', 'FL3'], - gate_contour=min_gate_contour, + hist_channels=['FL1', 'FL2'], + gate_contour=sfc1_gate_contour, density_params=density_params, - savefig=min_plot_filename) + savefig=sfc1_plot_filename) FlowCal.plot.density_and_hist( - max_beads_sample, - max_beads_sample_gated, + sfc2_beads_sample, + sfc2_beads_sample_gated, density_channels=['FSC', 'SSC'], - hist_channels=['FL1', 'FL3'], - gate_contour=max_gate_contour, + hist_channels=['FL1', 'FL2'], + gate_contour=sfc2_gate_contour, density_params=density_params, - savefig=max_plot_filename) + savefig=sfc2_plot_filename) # Use beads data to obtain a MEF transformation function - print("\nCalculating standard curve for channel FL1...") + print("\nCalculating standard curves for channels FL1 and FL2...") # ``FlowCal.mef.get_transform_fxn()`` generates a transformation function # that converts fluorescence from relative fluorescence units (RFI) to MEF. # This function uses bead data from ``beads_sample_gated``. We generate a - # MEF transformation function for channel FL1, with corresponding MEF - # fluorescence values specified by the array ``mefl_values``. In addition, - # we specify that clustering (subpopulation recognition) should be performed - # using information from both FL1 and FL3 channels. We also enable the + # MEF transformation function for channels FL1 and FL2, with corresponding + # MEF fluorescence values specified by the arrays ``mefl_values`` and + # ``mepe_values``. By default, clustering (subpopulation recognition) will + # be performed using information from both channels. We also enable the # ``verbose`` mode, which prints information of each step. Finally, we # instruct the function to generate plots of each step in the folder # specified in ``beads_plot_dir`` with the suffix "beads". mef_transform_fxn = FlowCal.mef.get_transform_fxn( beads_sample_gated, - mef_channels='FL1', - mef_values=mefl_values, - clustering_channels=['FL1', 'FL3'], + mef_channels=['FL1', 'FL2'], + mef_values=[mefl_values, mepe_values], verbose=True, plot=True, plot_dir=beads_plot_dir, plot_filename='beads') - min_mef_transform_fxn = FlowCal.mef.get_transform_fxn( - min_beads_sample_gated, - mef_channels='FL1', - mef_values=min_mefl_values, - clustering_channels=['FL1', 'FL3'], + nfc_mef_transform_fxn = FlowCal.mef.get_transform_fxn( + nfc_beads_sample_gated, + mef_channels=['FL1', 'FL2'], + mef_values=[nfc_mefl_values, nfc_mepe_values], verbose=True, plot=True, plot_dir=beads_plot_dir, - plot_filename='min_beads') - max_mef_transform_fxn = FlowCal.mef.get_transform_fxn( - max_beads_sample_gated, - mef_channels='FL1', - mef_values=max_mefl_values, - clustering_channels=['FL1', 'FL3'], + plot_filename='nfc_beads') + sfc1_mef_transform_fxn = FlowCal.mef.get_transform_fxn( + sfc1_beads_sample_gated, + mef_channels=['FL1', 'FL2'], + mef_values=[sfc1_mefl_values, sfc1_mepe_values], verbose=True, plot=True, plot_dir=beads_plot_dir, - plot_filename='max_beads') + plot_filename='sfc1_beads') + sfc2_mef_transform_fxn = FlowCal.mef.get_transform_fxn( + sfc2_beads_sample_gated, + mef_channels=['FL1', 'FL2'], + mef_values=[sfc2_mefl_values, sfc2_mepe_values], + verbose=True, + plot=True, + plot_dir=beads_plot_dir, + plot_filename='sfc2_beads') ### # Part 2: Processing cell sample data @@ -308,7 +351,7 @@ # We now use the transformation function obtained from calibration beads # to transform FL1 data from RFI to MEF. - sample = mef_transform_fxn(sample, channels=['FL1']) + sample = mef_transform_fxn(sample, channels=['FL1', 'FL2']) # Gating print("Performing gating...") @@ -322,7 +365,7 @@ # We only do this for the forward/side scatter channels and for # fluorescence channel FL1. sample_gated = FlowCal.gate.high_low(sample_gated, - channels=['FSC','SSC','FL1']) + channels=['FSC','SSC','FL1','FL2']) # Apply density gating on the forward/side scatter channels. Preserve # 85% of the events. Return also a contour around the gated region. @@ -343,9 +386,11 @@ density_params['mode'] = 'scatter' # Parameters for the fluorescence histograms - hist_params = {} - hist_params['xlabel'] = 'FL1 ' + \ + hist_params = [{}, {}] + hist_params[0]['xlabel'] = 'FL1 ' + \ '(Molecules of Equivalent Fluorescein, MEFL)' + hist_params[1]['xlabel'] = 'FL2 ' + \ + '(Molecules of Equivalent PE, MEPE)' # Plot filename # The figure can be saved in any format supported by matplotlib (svg, @@ -362,7 +407,7 @@ sample, sample_gated, density_channels=['FSC','SSC'], - hist_channels=['FL1'], + hist_channels=['FL1','FL2'], gate_contour=gate_contour, density_params=density_params, hist_params=hist_params, @@ -371,101 +416,160 @@ # Save cell sample object samples.append(sample_gated) - # Now, process the min and max control samples + # Now, process the nfc and sfc control samples print("\nProcessing control samples...") # Load, transform, and gate control samples - min_sample = FlowCal.io.FCSData(min_sample_filename) - max_sample = FlowCal.io.FCSData(max_sample_filename) + print("Loading file \"{}\"...".format(nfc_sample_filename)) + nfc_sample = FlowCal.io.FCSData(nfc_sample_filename) + print("Loading file \"{}\"...".format(sfc1_sample_filename)) + sfc1_sample = FlowCal.io.FCSData(sfc1_sample_filename) + print("Loading file \"{}\"...".format(sfc2_sample_filename)) + sfc2_sample = FlowCal.io.FCSData(sfc2_sample_filename) - min_sample = FlowCal.transform.to_rfi(min_sample) - max_sample = FlowCal.transform.to_rfi(max_sample) + print("Performing data transformation...") + nfc_sample = FlowCal.transform.to_rfi(nfc_sample) + sfc1_sample = FlowCal.transform.to_rfi(sfc1_sample) + sfc2_sample = FlowCal.transform.to_rfi(sfc2_sample) - min_sample = min_mef_transform_fxn(min_sample, channels=['FL1']) - max_sample = max_mef_transform_fxn(max_sample, channels=['FL1']) + nfc_sample = nfc_mef_transform_fxn(nfc_sample, channels=['FL1','FL2']) + sfc1_sample = sfc1_mef_transform_fxn(sfc1_sample, channels=['FL1','FL2']) + sfc2_sample = sfc2_mef_transform_fxn(sfc2_sample, channels=['FL1','FL2']) - min_sample_gated = FlowCal.gate.start_end(min_sample, - num_start=250, - num_end=100) - max_sample_gated = FlowCal.gate.start_end(max_sample, + print("Performing gating...") + nfc_sample_gated = FlowCal.gate.start_end(nfc_sample, num_start=250, num_end=100) + sfc1_sample_gated = FlowCal.gate.start_end(sfc1_sample, + num_start=250, + num_end=100) + sfc2_sample_gated = FlowCal.gate.start_end(sfc2_sample, + num_start=250, + num_end=100) + + nfc_sample_gated = FlowCal.gate.high_low(nfc_sample_gated, + channels=['FSC','SSC','FL1','FL2']) + sfc1_sample_gated = FlowCal.gate.high_low(sfc1_sample_gated, + channels=['FSC','SSC','FL1','FL2']) + sfc2_sample_gated = FlowCal.gate.high_low(sfc2_sample_gated, + channels=['FSC','SSC','FL1','FL2']) - min_sample_gated = FlowCal.gate.high_low(min_sample_gated, - channels=['FSC','SSC','FL1']) - max_sample_gated = FlowCal.gate.high_low(max_sample_gated, - channels=['FSC','SSC','FL1']) + print("Plotting density plot and histogram...") + density_gate_output = FlowCal.gate.density2d( + data=nfc_sample_gated, + channels=['FSC','SSC'], + gate_fraction=0.85, + full_output=True) + nfc_sample_gated = density_gate_output.gated_data + nfc_gate_contour = density_gate_output.contour - min_density_gate_output = FlowCal.gate.density2d( - data=min_sample_gated, + density_gate_output = FlowCal.gate.density2d( + data=sfc1_sample_gated, channels=['FSC','SSC'], gate_fraction=0.85, full_output=True) - min_sample_gated = min_density_gate_output.gated_data - min_gate_contour = min_density_gate_output.contour + sfc1_sample_gated = density_gate_output.gated_data + sfc1_gate_contour = density_gate_output.contour - max_density_gate_output = FlowCal.gate.density2d( - data=max_sample_gated, + density_gate_output = FlowCal.gate.density2d( + data=sfc2_sample_gated, channels=['FSC','SSC'], gate_fraction=0.85, full_output=True) - max_sample_gated = max_density_gate_output.gated_data - max_gate_contour = max_density_gate_output.contour + sfc2_sample_gated = density_gate_output.gated_data + sfc2_gate_contour = density_gate_output.contour # Plot and save - min_plot_filename = '{}/density_hist_min.png'.format(samples_plot_dir) - max_plot_filename = '{}/density_hist_max.png'.format(samples_plot_dir) + nfc_plot_filename = '{}/density_hist_nfc.png'.format(samples_plot_dir) + sfc1_plot_filename = '{}/density_hist_sfc1.png'.format(samples_plot_dir) + sfc2_plot_filename = '{}/density_hist_sfc2.png'.format(samples_plot_dir) FlowCal.plot.density_and_hist( - min_sample, - min_sample_gated, + nfc_sample, + nfc_sample_gated, + density_channels=['FSC','SSC'], + hist_channels=['FL1','FL2'], + gate_contour=nfc_gate_contour, + density_params=density_params, + hist_params=hist_params, + savefig=nfc_plot_filename) + FlowCal.plot.density_and_hist( + sfc1_sample, + sfc1_sample_gated, density_channels=['FSC','SSC'], - hist_channels=['FL1'], - gate_contour=min_gate_contour, + hist_channels=['FL1','FL2'], + gate_contour=sfc1_gate_contour, density_params=density_params, hist_params=hist_params, - savefig=min_plot_filename) + savefig=sfc1_plot_filename) FlowCal.plot.density_and_hist( - max_sample, - max_sample_gated, + sfc2_sample, + sfc2_sample_gated, density_channels=['FSC','SSC'], - hist_channels=['FL1'], - gate_contour=max_gate_contour, + hist_channels=['FL1','FL2'], + gate_contour=sfc2_gate_contour, density_params=density_params, hist_params=hist_params, - savefig=max_plot_filename) + savefig=sfc2_plot_filename) + + # Perform multi-color compensation + # ``FlowCal.compensate.get_transform_fxn()`` generates a transformation + # function that performs multi-color compensation on a specified set of + # channels in order to remove fluorophore bleedthrough. + # This function requires data from single-fluorophore controls (SFCs), one + # per channel to compensate, each from cells containing only one + # fluorophore. This function can optionally use data from a no-fluorophore + # control (NFC). + print("\nPerforming multi-color compensation...") + compensation_fxn = FlowCal.compensate.get_transform_fxn( + nfc_sample=nfc_sample_gated, + sfc_samples=[sfc1_sample_gated, sfc2_sample_gated], + comp_channels=['FL1', 'FL2'], + ) + # Compensate all samples + samples_compensated = [compensation_fxn(s, ['FL1', 'FL2']) for s in samples] + nfc_sample_compensated = compensation_fxn(nfc_sample_gated, ['FL1', 'FL2']) + sfc1_sample_compensated = compensation_fxn(sfc1_sample_gated, ['FL1', 'FL2']) + sfc2_sample_compensated = compensation_fxn(sfc2_sample_gated, ['FL1', 'FL2']) ### # Part 3: Examples on how to use processed cell sample data ### + # We now show how to generate plots using the processed flow cytometry + # data we just obtained. + print("\nGenerating plots...") # Plot 1: Histogram of all samples # - # Here, we plot the fluorescence histograms of all ten samples in the same - # figure, using ``FlowCal.plot.hist1d``. Note how this function can be used + # Here, we plot the fluorescence histograms of all nine samples in the same + # figure using ``FlowCal.plot.hist1d``. Note how this function can be used # in the context of accessory matplotlib functions to modify the axes # limits and labels and to add a legend, among other things. - # Color each histogram according to its DAPG concentration. Linearize the - # color transitions using a logarithmic normalization to match the - # logarithmic spacing of the DAPG concentrations. (Concentrations are also - # augmented slightly to move the 0.0 concentration into the log - # normalization range.) - cmap = mpl.cm.get_cmap('gray_r') - norm = mpl.colors.LogNorm(vmin=1e0, vmax=3500.) - colors = [cmap(norm(dapg_i+4.)) for dapg_i in dapg] - - plt.figure(figsize=(6,3.5)) + plt.figure(figsize=(6, 5.5)) + plt.subplot(2, 1, 1) FlowCal.plot.hist1d(samples, channel='FL1', histtype='step', - bins=128, - edgecolor=colors) + bins=128) plt.ylim((0,2500)) plt.xlim((0,5e4)) plt.xlabel('FL1 (Molecules of Equivalent Fluorescein, MEFL)') - plt.legend(['{} $\mu M$ DAPG'.format(i) for i in dapg], + plt.legend(['{:.1f} ng/mL aTc'.format(i) for i in atc], loc='upper left', fontsize='small') + + plt.subplot(2, 1, 2) + FlowCal.plot.hist1d(samples, + channel='FL2', + histtype='step', + bins=128) + plt.ylim((0,2500)) + plt.xlim((0,5e4)) + plt.xlabel('FL2 (Molecules of Equivalent PE, MEPE)') + plt.legend(['{:.1f} ng/mL aTc'.format(i) for i in atc], + loc='upper left', + fontsize='small') + plt.tight_layout() plt.savefig('histograms.png', dpi=200) plt.close() @@ -475,42 +579,65 @@ # Here, we illustrate how to obtain statistics from the fluorescence of # each sample and how to use them in a plot. The stats module contains # functions to calculate different statistics such as mean, median, and - # standard deviation. In this example, we calculate the mean from channel - # FL1 of each sample and plot them against the corresponding DAPG + # standard deviation. In this example, we calculate the mean from channels + # FL1 and FL2 of each sample and plot them against the corresponding aTc # concentrations. - samples_fluorescence = [FlowCal.stats.mean(s, channels='FL1') - for s in samples] - min_fluorescence = FlowCal.stats.mean(min_sample_gated, - channels='FL1') - max_fluorescence = FlowCal.stats.mean(max_sample_gated, - channels='FL1') - - dapg_color = '#ffc400' # common color used for DAPG-related plots - - plt.figure(figsize=(3,3)) - plt.plot(dapg, - samples_fluorescence, + samples_fl1 = [FlowCal.stats.mean(s, channels='FL1') for s in samples] + samples_fl2 = [FlowCal.stats.mean(s, channels='FL2') for s in samples] + # No fluorescence control (NFC) will give the minimum fluorescence level in + # both channels. Single fluorescence controls (SFCs) containing sfGFP or + # mCherry only will give the maximum levels in channels FL1 and FL2. + min_fl1 = FlowCal.stats.mean(nfc_sample_gated, channels='FL1') + max_fl1 = FlowCal.stats.mean(sfc1_sample_gated, channels='FL1') + min_fl2 = FlowCal.stats.mean(nfc_sample_gated, channels='FL2') + max_fl2 = FlowCal.stats.mean(sfc2_sample_gated, channels='FL2') + + plt.figure(figsize=(6, 3)) + + plt.subplot(1, 2, 1) + plt.plot(atc, + samples_fl1, marker='o', - color=dapg_color) - - # Illustrate min and max bounds - plt.axhline(min_fluorescence, + color='tab:green') + plt.axhline(min_fl1, color='gray', linestyle='--', zorder=-1) - plt.text(s='Min', x=2e2, y=1.6e2, ha='left', va='bottom', color='gray') - plt.axhline(max_fluorescence, + plt.text(s='Min', x=3e1, y=2e2, ha='left', va='bottom', color='gray') + plt.axhline(max_fl1, color='gray', linestyle='--', zorder=-1) - plt.text(s='Max', x=-0.7, y=5.2e3, ha='left', va='top', color='gray') - + plt.text(s='Max', x=-0.8, y=5.2e3, ha='left', va='top', color='gray') plt.yscale('log') plt.ylim((5e1,1e4)) plt.xscale('symlog') - plt.xlim((-1e0, 1e3)) - plt.xlabel('DAPG Concentration ($\mu M$)') + plt.xlim((-1e0, 1e2)) + plt.xlabel('aTc Concentration (ng/mL)') plt.ylabel('FL1 Fluorescence (MEFL)') + + plt.subplot(1, 2, 2) + plt.plot(atc, + samples_fl2, + marker='o', + color='tab:orange') + plt.axhline(min_fl2, + color='gray', + linestyle='--', + zorder=-1) + plt.text(s='Min', x=3e1, y=3.5e1, ha='left', va='bottom', color='gray') + plt.axhline(max_fl2, + color='gray', + linestyle='--', + zorder=-1) + plt.text(s='Max', x=-0.8, y=2.e3, ha='left', va='top', color='gray') + plt.yscale('log') + plt.ylim((1e1,5e3)) + plt.xscale('symlog') + plt.xlim((-1e0, 1e2)) + plt.xlabel('aTc Concentration (ng/mL)') + plt.ylabel('FL2 Fluorescence (MEPE)') + plt.tight_layout() plt.savefig('dose_response.png', dpi=200) plt.close() @@ -518,62 +645,114 @@ # Plot 3: Dose response violin plot # # Here, we use a violin plot to show the fluorescence of (almost) all - # cells as a function of DAPG. (The `upper_trim_fraction` and + # cells as a function of aTc. (The `upper_trim_fraction` and # `lower_trim_fraction` parameters eliminate the top and bottom 1% of # cells from each violin for aesthetic reasons. The summary statistic, # which is illustrated as a horizontal line atop each violin, is - # calculated before cells are removed, though.) We set `yscale` to 'log' + # calculated before cells are removed, though). We set `yscale` to 'log' # because the cytometer used to collect this data produces positive # integer data (as opposed to floating-point data, which can sometimes be # negative), so the added complexity of a logicle y-scale (which is the # default) is not necessary. + plt.figure(figsize=(8, 3.5)) - # FlowCal violin plots can also illustrate a mathematical model alongside - # the violins. To take advantage of this feature, we first recapitulate a - # model of the fluorescent protein (sfGFP) fluorescence produced by this - # DAPG sensor as a function of DAPG (from this study: - # https://doi.org/10.15252/msb.20209618). sfGFP fluorescence is cellular - # fluorescence minus autofluorescence. - def dapg_sensor_output(dapg_concentration): - mn = 86. - mx = 3147. - K = 20. - n = 3.57 - if dapg_concentration <= 0: - return mn - else: - return mn + ((mx-mn)/(1+((K/dapg_concentration)**n))) - - # To model cellular fluorescence, which we are plotting with this violin - # plot, we must add autofluorescence back to the sfGFP signal. For this - # model, autofluorescence is the mean fluorescence of an E. coli strain - # lacking sfGFP, which our min control is. - autofluorescence = FlowCal.stats.mean(min_sample_gated, channels='FL1') - def dapg_sensor_cellular_fluorescence(dapg_concentration): - return dapg_sensor_output(dapg_concentration) + autofluorescence - - plt.figure(figsize=(4,3.5)) + plt.subplot(1, 2, 1) FlowCal.plot.violin_dose_response( data=samples, channel='FL1', - positions=dapg, - min_data=min_sample_gated, - max_data=max_sample_gated, - model_fxn=dapg_sensor_cellular_fluorescence, - violin_kwargs={'facecolor':dapg_color, - 'edgecolor':'black'}, - violin_width_to_span_fraction=0.075, + positions=atc, + min_data=nfc_sample_gated, + max_data=sfc1_sample_gated, + xlabel='aTc Concentration (ng/mL)', xscale='log', yscale='log', - ylim=(1e1, 3e4), + ylim=(1e1,1e4), + violin_width=0.12, + violin_kwargs={'facecolor': 'tab:green', + 'edgecolor':'black'}, draw_model_kwargs={'color':'gray', 'linewidth':3, 'zorder':-1, - 'solid_capstyle':'butt'}) - plt.xlabel('DAPG Concentration ($\mu M$)') + 'solid_capstyle':'butt'}, + ) plt.ylabel('FL1 Fluorescence (MEFL)') + + plt.subplot(1, 2, 2) + FlowCal.plot.violin_dose_response( + data=samples, + channel='FL2', + positions=atc, + min_data=nfc_sample_gated, + max_data=sfc2_sample_gated, + xlabel='aTc Concentration (ng/mL)', + xscale='log', + yscale='log', + ylim=(1e0,1e4), + violin_width=0.12, + violin_kwargs={'facecolor': 'tab:orange', + 'edgecolor':'black'}, + draw_model_kwargs={'color':'gray', + 'linewidth':3, + 'zorder':-1, + 'solid_capstyle':'butt'}, + ) + plt.ylabel('FL2 Fluorescence (MEPE)') + plt.tight_layout() plt.savefig('dose_response_violin.png', dpi=200) plt.close() + # Plot 4: Dose response violin plot of compensated data + # + # Here, we repeat the previous violin plot but using compensated data. + # y axis will now be plotted in ``logicle`` scale since histograms will + # be centered around zero due to compensation. + plt.figure(figsize=(8, 3.5)) + + plt.subplot(1, 2, 1) + FlowCal.plot.violin_dose_response( + data=samples_compensated, + channel='FL1', + positions=atc, + min_data=nfc_sample_compensated, + max_data=sfc1_sample_compensated, + xlabel='aTc Concentration (ng/mL)', + xscale='log', + yscale='logicle', + ylim=(-3e2, 1e4), + violin_width=0.12, + violin_kwargs={'facecolor': 'tab:green', + 'edgecolor':'black'}, + draw_model_kwargs={'color':'gray', + 'linewidth':3, + 'zorder':-1, + 'solid_capstyle':'butt'}, + ) + plt.ylabel('FL1 Fluorescence (MEFL)') + + plt.subplot(1, 2, 2) + FlowCal.plot.violin_dose_response( + data=samples_compensated, + channel='FL2', + positions=atc, + min_data=nfc_sample_compensated, + max_data=sfc2_sample_compensated, + xlabel='aTc Concentration (ng/mL)', + xscale='log', + yscale='logicle', + ylim=(-1e2, 1e4), + violin_width=0.12, + violin_kwargs={'facecolor': 'tab:orange', + 'edgecolor':'black'}, + draw_model_kwargs={'color':'gray', + 'linewidth':3, + 'zorder':-1, + 'solid_capstyle':'butt'}, + ) + plt.ylabel('FL2 Fluorescence (MEPE)') + + plt.tight_layout() + plt.savefig('dose_response_violin_compensated.png', dpi=200) + plt.close() + print("\nDone.") diff --git a/examples/analyze_no_mef.py b/examples/analyze_no_mef.py index eb3e167..41b3bb6 100644 --- a/examples/analyze_no_mef.py +++ b/examples/analyze_no_mef.py @@ -2,8 +2,10 @@ """ FlowCal Python API example, without using calibration beads data. -This script is divided in two parts. Part one processes data from ten cell -samples and generates plots of each one. +This script is divided in two parts. Part one processes data from nine +cell samples and generates plots of each one. In addition, multi-color +compensation is performed on all samples using data from no-fluorophore and +single-fluorophore control samples (NFC and SFCs). Part two exemplifies how to use the processed cell sample data with FlowCal's plotting and statistics modules to produce interesting plots. @@ -24,19 +26,21 @@ ### # Names of the FCS files containing data from cell samples -samples_filenames = ['FCFiles/sample006.fcs', - 'FCFiles/sample007.fcs', - 'FCFiles/sample008.fcs', - 'FCFiles/sample009.fcs', - 'FCFiles/sample010.fcs', - 'FCFiles/sample011.fcs', - 'FCFiles/sample012.fcs', - 'FCFiles/sample013.fcs', - 'FCFiles/sample014.fcs', - 'FCFiles/sample015.fcs'] - -# DAPG concentration of each cell sample, in micromolar. -dapg = np.array([0, 2.33, 4.36, 8.16, 15.3, 28.6, 53.5, 100, 187, 350]) +samples_filenames = ['FCFiles/sample029.fcs', + 'FCFiles/sample030.fcs', + 'FCFiles/sample031.fcs', + 'FCFiles/sample032.fcs', + 'FCFiles/sample033.fcs', + 'FCFiles/sample034.fcs', + 'FCFiles/sample035.fcs', + 'FCFiles/sample036.fcs', + 'FCFiles/sample037.fcs'] +nfc_sample_filename = 'FCFiles/controls/nfc/sample004.fcs' +sfc1_sample_filename = 'FCFiles/controls/sfc1/sample007.fcs' +sfc2_sample_filename = 'FCFiles/controls/sfc2/sample019.fcs' + +# aTc concentration of each cell sample, in ng/mL. +atc = np.array([0, 0.5, 1, 1.5, 2, 3, 4, 7.5, 20]) # Plots will be generated after gating and transforming cell samples. These # will be stored in the following folder. @@ -105,7 +109,7 @@ # We will remove saturated events in the forward/side scatter channels, # and in the fluorescence channel FL1. sample_gated = FlowCal.gate.high_low(sample_gated, - channels=['FSC','SSC','FL1']) + channels=['FSC','SSC','FL1','FL2']) # ``FlowCal.gate.density2d()`` preserves only the densest population as # seen in a 2D density diagram of two channels. This helps remove @@ -136,8 +140,9 @@ density_params['mode'] = 'scatter' # Parameters for the fluorescence histograms - hist_params = {} - hist_params['xlabel'] = 'FL1 Fluorescence (a.u.)' + hist_params = [{}, {}] + hist_params[0]['xlabel'] = 'FL1 Fluorescence (a.u.)' + hist_params[1]['xlabel'] = 'FL2 Fluorescence (a.u.)' # Plot filename # The figure can be saved in any format supported by matplotlib (svg, @@ -163,7 +168,7 @@ sample, sample_gated, density_channels=['FSC','SSC'], - hist_channels=['FL1'], + hist_channels=['FL1','FL2'], gate_contour=gate_contour, density_params=density_params, hist_params=hist_params, @@ -172,9 +177,124 @@ # Save cell sample object samples.append(sample_gated) + # Now, process the nfc and sfc control samples + print("\nProcessing control samples...") + + # Load, transform, and gate control samples + print("Loading file \"{}\"...".format(nfc_sample_filename)) + nfc_sample = FlowCal.io.FCSData(nfc_sample_filename) + print("Loading file \"{}\"...".format(sfc1_sample_filename)) + sfc1_sample = FlowCal.io.FCSData(sfc1_sample_filename) + print("Loading file \"{}\"...".format(sfc2_sample_filename)) + sfc2_sample = FlowCal.io.FCSData(sfc2_sample_filename) + + print("Performing data transformation...") + nfc_sample = FlowCal.transform.to_rfi(nfc_sample) + sfc1_sample = FlowCal.transform.to_rfi(sfc1_sample) + sfc2_sample = FlowCal.transform.to_rfi(sfc2_sample) + + print("Performing gating...") + nfc_sample_gated = FlowCal.gate.start_end(nfc_sample, + num_start=250, + num_end=100) + sfc1_sample_gated = FlowCal.gate.start_end(sfc1_sample, + num_start=250, + num_end=100) + sfc2_sample_gated = FlowCal.gate.start_end(sfc2_sample, + num_start=250, + num_end=100) + + nfc_sample_gated = FlowCal.gate.high_low(nfc_sample_gated, + channels=['FSC','SSC','FL1','FL2']) + sfc1_sample_gated = FlowCal.gate.high_low(sfc1_sample_gated, + channels=['FSC','SSC','FL1','FL2']) + sfc2_sample_gated = FlowCal.gate.high_low(sfc2_sample_gated, + channels=['FSC','SSC','FL1','FL2']) + + print("Plotting density plot and histogram...") + density_gate_output = FlowCal.gate.density2d( + data=nfc_sample_gated, + channels=['FSC','SSC'], + gate_fraction=0.85, + full_output=True) + nfc_sample_gated = density_gate_output.gated_data + nfc_gate_contour = density_gate_output.contour + + density_gate_output = FlowCal.gate.density2d( + data=sfc1_sample_gated, + channels=['FSC','SSC'], + gate_fraction=0.85, + full_output=True) + sfc1_sample_gated = density_gate_output.gated_data + sfc1_gate_contour = density_gate_output.contour + + density_gate_output = FlowCal.gate.density2d( + data=sfc2_sample_gated, + channels=['FSC','SSC'], + gate_fraction=0.85, + full_output=True) + sfc2_sample_gated = density_gate_output.gated_data + sfc2_gate_contour = density_gate_output.contour + + # Plot and save + nfc_plot_filename = '{}/density_hist_nfc.png'.format(samples_plot_dir) + sfc1_plot_filename = '{}/density_hist_sfc1.png'.format(samples_plot_dir) + sfc2_plot_filename = '{}/density_hist_sfc2.png'.format(samples_plot_dir) + + FlowCal.plot.density_and_hist( + nfc_sample, + nfc_sample_gated, + density_channels=['FSC','SSC'], + hist_channels=['FL1','FL2'], + gate_contour=nfc_gate_contour, + density_params=density_params, + hist_params=hist_params, + savefig=nfc_plot_filename) + FlowCal.plot.density_and_hist( + sfc1_sample, + sfc1_sample_gated, + density_channels=['FSC','SSC'], + hist_channels=['FL1','FL2'], + gate_contour=sfc1_gate_contour, + density_params=density_params, + hist_params=hist_params, + savefig=sfc1_plot_filename) + FlowCal.plot.density_and_hist( + sfc2_sample, + sfc2_sample_gated, + density_channels=['FSC','SSC'], + hist_channels=['FL1','FL2'], + gate_contour=sfc2_gate_contour, + density_params=density_params, + hist_params=hist_params, + savefig=sfc2_plot_filename) + + # Perform multi-color compensation + # ``FlowCal.compensate.get_transform_fxn()`` generates a transformation + # function that performs multi-color compensation on a specified set of + # channels in order to remove fluorophore bleedthrough. + # This function requires data from single-fluorophore controls (SFCs), one + # per channel to compensate, each from cells containing only one + # fluorophore. This function can optionally use data from a no-fluorophore + # control (NFC). + print("\nPerforming multi-color compensation...") + compensation_fxn = FlowCal.compensate.get_transform_fxn( + nfc_sample=nfc_sample_gated, + sfc_samples=[sfc1_sample_gated, sfc2_sample_gated], + comp_channels=['FL1', 'FL2'], + ) + # Compensate all samples + samples_compensated = [compensation_fxn(s, ['FL1', 'FL2']) for s in samples] + nfc_sample_compensated = compensation_fxn(nfc_sample_gated, ['FL1', 'FL2']) + sfc1_sample_compensated = compensation_fxn(sfc1_sample_gated, ['FL1', 'FL2']) + sfc2_sample_compensated = compensation_fxn(sfc2_sample_gated, ['FL1', 'FL2']) + ### # Part 3: Examples on how to use processed cell sample data ### + # We now show how to generate plots using the processed flow cytometry + # data we just obtained. + print("\nGenerating plots...") # Plot 1: Histogram of all samples # @@ -183,27 +303,31 @@ # in the context of accessory matplotlib functions to modify the axes # limits and labels and to add a legend, among other things. - # Color each histogram according to its DAPG concentration. Linearize the - # color transitions using a logarithmic normalization to match the - # logarithmic spacing of the DAPG concentrations. (Concentrations are also - # augmented slightly to move the 0.0 concentration into the log - # normalization range.) - cmap = mpl.cm.get_cmap('gray_r') - norm = mpl.colors.LogNorm(vmin=1e0, vmax=3500.) - colors = [cmap(norm(dapg_i+4.)) for dapg_i in dapg] - - plt.figure(figsize=(6,3.5)) + plt.figure(figsize=(6, 5.5)) + plt.subplot(2, 1, 1) FlowCal.plot.hist1d(samples, channel='FL1', histtype='step', - bins=128, - edgecolor=colors) + bins=128) plt.ylim((0,2500)) plt.xlim((0,5e3)) plt.xlabel('FL1 Fluorescence (a.u.)') - plt.legend(['{} $\mu M$ DAPG'.format(i) for i in dapg], + plt.legend(['{:.1f} ng/mL aTc'.format(i) for i in atc], + loc='upper left', + fontsize='small') + + plt.subplot(2, 1, 2) + FlowCal.plot.hist1d(samples, + channel='FL2', + histtype='step', + bins=128) + plt.ylim((0,2500)) + plt.xlim((0,5e3)) + plt.xlabel('FL2 Fluorescence (a.u.)') + plt.legend(['{:.1f} ng/mL aTc'.format(i) for i in atc], loc='upper left', fontsize='small') + plt.tight_layout() plt.savefig('histograms.png', dpi=200) plt.close() @@ -214,40 +338,68 @@ # each sample and how to use them in a plot. The stats module contains # functions to calculate different statistics such as mean, median, and # standard deviation. In this example, we calculate the mean from channel - # FL1 of each sample and plot them against the corresponding DAPG + # FL1 of each sample and plot them against the corresponding aTc # concentrations. - samples_fluorescence = [FlowCal.stats.mean(s, channels='FL1') - for s in samples] - - dapg_color = '#ffc400' # common color used for DAPG-related plots - plt.figure(figsize=(3,3)) - plt.plot(dapg, - samples_fluorescence, + # Because some of our control samples were measured at a different cytometer + # gain setting and we aren't using MEF calibration here, we will use the 0 + # and 20 ng/mL aTc concentration samples instead. + samples_fl1 = [FlowCal.stats.mean(s, channels='FL1') for s in samples] + samples_fl2 = [FlowCal.stats.mean(s, channels='FL2') for s in samples] + # No fluorescence control (NFC) will give the minimum fluorescence level in + # both channels. Single fluorescence controls (SFCs) containing sfGFP or + # mCherry only will give the maximum levels in channels FL1 and FL2. + min_fl1 = FlowCal.stats.mean(nfc_sample_gated, channels='FL1') + max_fl1 = FlowCal.stats.mean(sfc1_sample_gated, channels='FL1') + min_fl2 = FlowCal.stats.mean(nfc_sample_gated, channels='FL2') + max_fl2 = FlowCal.stats.mean(sfc2_sample_gated, channels='FL2') + + plt.figure(figsize=(6,3)) + + plt.subplot(1, 2, 1) + plt.plot(atc, + samples_fl1, marker='o', - color=dapg_color) - - # Illustrate min and max bounds. Because some of our control samples were - # measured at a different cytometer gain setting and we aren't using MEF - # calibration here, we will use the 0uM and 350uM DAPG concentration - # samples instead. - plt.axhline(samples_fluorescence[0], + color='tab:green') + plt.axhline(min_fl1, color='gray', linestyle='--', zorder=-1) - plt.text(s='Min', x=2e2, y=2.0e1, ha='left', va='bottom', color='gray') - plt.axhline(samples_fluorescence[-1], + plt.text(s='Min', x=3e1, y=1.4e1, ha='left', va='bottom', color='gray') + plt.axhline(max_fl1, color='gray', linestyle='--', zorder=-1) - plt.text(s='Max', x=-0.7, y=2.1e2, ha='left', va='top', color='gray') - + plt.text(s='Max', x=-0.8, y=3.3e2, ha='left', va='top', color='gray') plt.yscale('log') - plt.ylim((5e0,5e2)) + plt.ylim((5e0, 5e2)) plt.xscale('symlog') - plt.xlim((-1e0, 1e3)) - plt.xlabel('DAPG Concentration ($\mu M$)') + plt.xlim((-1e0, 1e2)) + plt.xlabel('aTc Concentration (ng/mL)') plt.ylabel('FL1 Fluorescence (a.u.)') + + plt.subplot(1, 2, 2) + plt.plot(atc, + samples_fl2, + marker='o', + color='tab:orange') + plt.axhline(min_fl2, + color='gray', + linestyle='--', + zorder=-1) + plt.text(s='Min', x=3e1, y=0.9e1, ha='left', va='bottom', color='gray') + plt.axhline(max_fl2, + color='gray', + linestyle='--', + zorder=-1) + plt.text(s='Max', x=-0.8, y=5e2, ha='left', va='top', color='gray') + plt.yscale('log') + plt.ylim((4e0, 1.5e3)) + plt.xscale('symlog') + plt.xlim((-1e0, 1e2)) + plt.xlabel('aTc Concentration (ng/mL)') + plt.ylabel('FL2 Fluorescence (a.u.)') + plt.tight_layout() plt.savefig('dose_response.png', dpi=200) plt.close() @@ -255,33 +407,104 @@ # Plot 3: Dose response violin plot # # Here, we use a violin plot to show the fluorescence of (almost) all - # cells as a function of DAPG. (The `upper_trim_fraction` and + # cells as a function of aTc. (The `upper_trim_fraction` and # `lower_trim_fraction` parameters eliminate the top and bottom 1% of # cells from each violin for aesthetic reasons. The summary statistic, # which is illustrated as a horizontal line atop each violin, is - # calculated before cells are removed, though.) We again use the 0uM and - # 350uM DAPG concentration samples as the min and max data in lieu of - # controls. We also set `yscale` to 'log' because the cytometer used to - # collect this data produces positive integer data (as opposed to - # floating-point data, which can sometimes be negative), so the added - # complexity of a logicle y-scale (which is the default) is not necessary. - plt.figure(figsize=(4,3.5)) + # calculated before cells are removed, though.) We set `yscale` to 'log' + # because the cytometer used to collect this data produces positive + # integer data (as opposed to floating-point data, which can sometimes be + # negative), so the added complexity of a logicle y-scale (which is the + # default) is not necessary. + plt.figure(figsize=(8, 3.5)) + + plt.subplot(1, 2, 1) FlowCal.plot.violin_dose_response( data=samples, channel='FL1', - positions=dapg, - min_data=samples[0], - max_data=samples[-1], - violin_kwargs={'facecolor':dapg_color, + positions=atc, + min_data=nfc_sample_gated, + max_data=sfc1_sample_gated, + violin_kwargs={'facecolor':'tab:green', 'edgecolor':'black'}, violin_width_to_span_fraction=0.075, xscale='log', yscale='log', - ylim=(1e0,2e3)) - plt.xlabel('DAPG Concentration ($\mu M$)') + ylim=(1e0,1e3)) + plt.xlabel('aTc Concentration (ng/mL)') plt.ylabel('FL1 Fluorescence (a.u.)') + + plt.subplot(1, 2, 2) + FlowCal.plot.violin_dose_response( + data=samples, + channel='FL2', + positions=atc, + min_data=nfc_sample_gated, + max_data=sfc2_sample_gated, + violin_kwargs={'facecolor':'tab:orange', + 'edgecolor':'black'}, + violin_width_to_span_fraction=0.075, + xscale='log', + yscale='log', + ylim=(1e0,2e3)) + plt.xlabel('aTc Concentration (ng/mL)') + plt.ylabel('FL2 Fluorescence (a.u.)') + plt.tight_layout() plt.savefig('dose_response_violin.png', dpi=200) plt.close() + # Plot 4: Dose response violin plot of compensated data + # + # Here, we repeat the previous violin plot but using compensated data. + # y axis will now be plotted in ``logicle`` scale since histograms will + # be centered around zero due to compensation. + plt.figure(figsize=(8, 3.5)) + + plt.subplot(1, 2, 1) + FlowCal.plot.violin_dose_response( + data=samples_compensated, + channel='FL1', + positions=atc, + min_data=nfc_sample_compensated, + max_data=sfc1_sample_compensated, + xlabel='aTc Concentration (ng/mL)', + xscale='log', + yscale='logicle', + ylim=(-3e1, 1e3), + violin_width=0.12, + violin_kwargs={'facecolor': 'tab:green', + 'edgecolor':'black'}, + draw_model_kwargs={'color':'gray', + 'linewidth':3, + 'zorder':-1, + 'solid_capstyle':'butt'}, + ) + plt.ylabel('FL1 Fluorescence (a.u.)') + + plt.subplot(1, 2, 2) + FlowCal.plot.violin_dose_response( + data=samples_compensated, + channel='FL2', + positions=atc, + min_data=nfc_sample_compensated, + max_data=sfc2_sample_compensated, + xlabel='aTc Concentration (ng/mL)', + xscale='log', + yscale='logicle', + ylim=(-3e1, 2e3), + violin_width=0.12, + violin_kwargs={'facecolor': 'tab:orange', + 'edgecolor':'black'}, + draw_model_kwargs={'color':'gray', + 'linewidth':3, + 'zorder':-1, + 'solid_capstyle':'butt'}, + ) + plt.ylabel('FL2 Fluorescence (a.u.)') + + plt.tight_layout() + plt.savefig('dose_response_violin_compensated.png', dpi=200) + plt.close() + print("\nDone.") diff --git a/examples/experiment.xlsx b/examples/experiment.xlsx index 5efdf77..8c9bdf2 100644 Binary files a/examples/experiment.xlsx and b/examples/experiment.xlsx differ diff --git a/test/test_transform.py b/test/test_transform.py index f1467ba..7ee03db 100644 --- a/test/test_transform.py +++ b/test/test_transform.py @@ -560,5 +560,276 @@ def test_mef_bins_channels(self): ] np.testing.assert_array_equal(vit, vo) +class TestCompensateArray(unittest.TestCase): + def setUp(self): + self.d = np.array([ + [1, 7, 2], + [2, 8, 3], + [3, 9, 4], + [4, 10, 5], + [5, 1, 6], + [6, 2, 7], + [7, 3, 8], + [8, 4, 9], + [9, 5, 10], + [10, 6, 1], + ]) + self.a0 = np.array([1, 0.5, 2.1]) + self.A = np.array([[1, 0.3, 0.2], [0.6, 1, 0.5], [0.7, 0.8, 1]]) + self.dt_complete = np.linalg.solve(self.A, (self.d - self.a0).T).T + + def test_length_error_1(self): + self.assertRaises(ValueError, FlowCal.transform.to_compensated, + self.d, [0, 1, 2], self.a0, self.A, [0, 1]) + + def test_length_error_2(self): + self.assertRaises(ValueError, FlowCal.transform.to_compensated, + self.d, [0, 1, 2], self.a0[0:1], self.A, [0, 1, 2]) + + def test_length_error_3(self): + self.assertRaises(ValueError, FlowCal.transform.to_compensated, + self.d, [0, 1, 2], self.a0, self.A[0:2,0:2], [0, 1, 2]) + + def test_channel_error(self): + self.assertRaises(ValueError, FlowCal.transform.to_compensated, + self.d, [0, 2], self.a0[0:2], self.A[0:2,0:2], [0, 1]) + + def test_compensation_full(self): + dt = FlowCal.transform.to_compensated( + self.d, [0, 1, 2], self.a0, self.A, [0, 1, 2]) + np.testing.assert_array_almost_equal( + dt, + np.linalg.solve(self.A, (self.d - self.a0).T).T) + + def test_compensation_full_channels_none(self): + dt = FlowCal.transform.to_compensated( + self.d, None, self.a0, self.A, [0, 1, 2]) + np.testing.assert_array_almost_equal(dt, self.dt_complete) + + def test_compensation_full_comp_channels_none(self): + dt = FlowCal.transform.to_compensated( + self.d, None, self.a0, self.A, None) + np.testing.assert_array_almost_equal(dt, self.dt_complete) + + def test_compensation_subset_channels_1(self): + dt = FlowCal.transform.to_compensated( + self.d, [0], self.a0, self.A, [0, 1, 2]) + np.testing.assert_array_almost_equal(dt[:,0], self.dt_complete[:,0]) + np.testing.assert_array_almost_equal(dt[:,[1,2]], self.d[:,[1,2]]) + + def test_compensation_subset_channels_2(self): + dt = FlowCal.transform.to_compensated( + self.d, [0, 1], self.a0, self.A, [0, 1, 2]) + np.testing.assert_array_almost_equal( + dt[:,[0, 1]], self.dt_complete[:,[0, 1]]) + np.testing.assert_array_almost_equal(dt[:,2], self.d[:,2]) + + def test_compensation_subset_channels_3(self): + dt = FlowCal.transform.to_compensated( + self.d, [0, 2], self.a0, self.A, [0, 1, 2]) + np.testing.assert_array_almost_equal( + dt[:,[0, 2]], self.dt_complete[:,[0, 2]]) + np.testing.assert_array_almost_equal(dt[:,1], self.d[:,1]) + +class TestCompensateFCS(unittest.TestCase): + def setUp(self): + # Data + self.d = FlowCal.io.FCSData('test/Data001.fcs') + self.channel_names = ['FSC-H', 'SSC-H', 'FL1-H', + 'FL2-H', 'FL3-H', 'Time'] + self.n_samples = self.d.shape[0] + # Compensation coefficients + self.a0 = np.array([1, 0.5, 2.1]) + self.A = np.array([[1, 0.3, 0.2], [0.6, 1, 0.5], [0.7, 0.8, 1]]) + # Expected results + self.dt_complete = np.linalg.solve( + self.A, + (self.d[:, ['FL1-H', 'FL2-H', 'FL3-H']] - self.a0).T, + ).T + range_initial = self.d.range(['FL1-H', 'FL2-H', 'FL3-H']) + self.range_complete = np.linalg.solve( + self.A, + (np.array(range_initial).T - self.a0).T, + ) + + def test_length_error_1(self): + self.assertRaises( + ValueError, + FlowCal.transform.to_compensated, + self.d, + ['FL1-H', 'FL2-H', 'FL3-H'], + self.a0, + self.A, + ['FL1-H', 'FL2-H'], + ) + + def test_length_error_2(self): + self.assertRaises( + ValueError, + FlowCal.transform.to_compensated, + self.d, + ['FL1-H', 'FL2-H', 'FL3-H'], + self.a0[0:1], + self.A, + ['FL1-H', 'FL2-H', 'FL3-H'], + ) + + def test_length_error_3(self): + self.assertRaises( + ValueError, + FlowCal.transform.to_compensated, + self.d, + ['FL1-H', 'FL2-H', 'FL3-H'], + self.a0, + self.A[0:2,0:2], + ['FL1-H', 'FL2-H', 'FL3-H'], + ) + + def test_channel_error(self): + self.assertRaises( + ValueError, + FlowCal.transform.to_compensated, + self.d, + ['FL1-H', 'FL3-H'], + self.a0[0:2], + self.A[0:2,0:2], + ['FL1-H', 'FL2-H'], + ) + + def test_compensation_full(self): + dt = FlowCal.transform.to_compensated( + self.d, + ['FL1-H', 'FL2-H', 'FL3-H'], + self.a0, + self.A, + ['FL1-H', 'FL2-H', 'FL3-H'], + ) + # Check events + np.testing.assert_array_almost_equal( + dt[:, ['FL1-H', 'FL2-H', 'FL3-H']], + self.dt_complete, + ) + np.testing.assert_array_almost_equal( + dt[:, ['FSC-H', 'SSC-H', 'Time']], + self.d[:, ['FSC-H', 'SSC-H', 'Time']], + ) + # Check range + np.testing.assert_array_almost_equal( + dt.range(['FL1-H', 'FL2-H', 'FL3-H']), + self.range_complete, + ) + np.testing.assert_array_almost_equal( + dt.range(['FSC-H', 'SSC-H', 'Time']), + self.d.range(['FSC-H', 'SSC-H', 'Time']), + ) + + def test_compensation_full_channels_none(self): + dt = FlowCal.transform.to_compensated( + self.d, + None, + self.a0, + self.A, + ['FL1-H', 'FL2-H', 'FL3-H'], + ) + # Check events + np.testing.assert_array_almost_equal( + dt[:, ['FL1-H', 'FL2-H', 'FL3-H']], + self.dt_complete, + ) + np.testing.assert_array_almost_equal( + dt[:, ['FSC-H', 'SSC-H', 'Time']], + self.d[:, ['FSC-H', 'SSC-H', 'Time']], + ) + # Check range + np.testing.assert_array_almost_equal( + dt.range(['FL1-H', 'FL2-H', 'FL3-H']), + self.range_complete, + ) + np.testing.assert_array_almost_equal( + dt.range(['FSC-H', 'SSC-H', 'Time']), + self.d.range(['FSC-H', 'SSC-H', 'Time']), + ) + + def test_compensation_subset_channels_1(self): + dt = FlowCal.transform.to_compensated( + self.d, + ['FL1-H'], + self.a0, + self.A, + ['FL1-H', 'FL2-H', 'FL3-H'], + ) + # Check events + np.testing.assert_array_almost_equal( + dt[:, 'FL1-H'], + self.dt_complete[:, 0], + ) + np.testing.assert_array_almost_equal( + dt[:, ['FSC-H', 'SSC-H', 'FL2-H', 'FL3-H', 'Time']], + self.d[:, ['FSC-H', 'SSC-H', 'FL2-H', 'FL3-H', 'Time']], + ) + # Check range + np.testing.assert_array_almost_equal( + dt.range('FL1-H'), + self.range_complete[0,:], + ) + np.testing.assert_array_almost_equal( + dt.range(['FSC-H', 'SSC-H', 'FL2-H', 'FL3-H', 'Time']), + self.d.range(['FSC-H', 'SSC-H', 'FL2-H', 'FL3-H', 'Time']), + ) + + def test_compensation_subset_channels_2(self): + dt = FlowCal.transform.to_compensated( + self.d, + ['FL1-H', 'FL2-H'], + self.a0, + self.A, + ['FL1-H', 'FL2-H', 'FL3-H'], + ) + # Check events + np.testing.assert_array_almost_equal( + dt[:, ['FL1-H', 'FL2-H']], + self.dt_complete[:, [0, 1]], + ) + np.testing.assert_array_almost_equal( + dt[:, ['FSC-H', 'SSC-H', 'FL3-H', 'Time']], + self.d[:, ['FSC-H', 'SSC-H', 'FL3-H', 'Time']], + ) + # Check range + np.testing.assert_array_almost_equal( + dt.range(['FL1-H', 'FL2-H']), + self.range_complete[[0, 1], :], + ) + np.testing.assert_array_almost_equal( + dt.range(['FSC-H', 'SSC-H', 'FL3-H', 'Time']), + self.d.range(['FSC-H', 'SSC-H', 'FL3-H', 'Time']), + ) + + def test_compensation_subset_channels_3(self): + dt = FlowCal.transform.to_compensated( + self.d, + ['FL1-H', 'FL3-H'], + self.a0, + self.A, + ['FL1-H', 'FL2-H', 'FL3-H'], + ) + # Check events + np.testing.assert_array_almost_equal( + dt[:, ['FL1-H', 'FL3-H']], + self.dt_complete[:, [0, 2]], + ) + np.testing.assert_array_almost_equal( + dt[:, ['FSC-H', 'SSC-H', 'FL2-H', 'Time']], + self.d[:, ['FSC-H', 'SSC-H', 'FL2-H', 'Time']], + ) + # Check range + np.testing.assert_array_almost_equal( + dt.range(['FL1-H', 'FL3-H']), + self.range_complete[[0, 2], :], + ) + np.testing.assert_array_almost_equal( + dt.range(['FSC-H', 'SSC-H', 'FL2-H', 'Time']), + self.d.range(['FSC-H', 'SSC-H', 'FL2-H', 'Time']), + ) + if __name__ == '__main__': unittest.main()