Source code for cytoflow.operations.autofluorescence

#!/usr/bin/env python3.8
# coding: latin-1

# (c) Massachusetts Institute of Technology 2015-2018
# (c) Brian Teague 2018-2022
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

"""
cytoflow.operations.autofluorescence
------------------------------------

The `cytoflow.operations.autofluorescence` module contains two classes:

`AutofluorescenceOp` -- corrects an `Experiment` for
autofluorescence
  
`AutofluorescenceDiagnosticView` -- a diagnostic to make sure that
`AutofluorescenceOp` estimated its parameters correctly.
"""

from warnings import warn
from traits.api import (HasStrictTraits, Str, Float, File, Dict,
                        Instance, List, Constant, Tuple, Array, provides)
                       
import numpy as np

import cytoflow.views
import cytoflow.utility as util

from .i_operation import IOperation
from .import_op import Tube, ImportOp, check_tube

[docs]@provides(IOperation) class AutofluorescenceOp(HasStrictTraits): """ Apply autofluorescence correction to a set of fluorescence channels. The `estimate` function loads a separate FCS file (not part of the input `Experiment`) and computes the untransformed median and standard deviation of the blank cells. Then, `apply` subtracts the median from the experiment data. To use, set the `blank_file` property to point to an FCS file with unstained or nonfluorescing cells in it; set the `channels` property to a list of channels to correct. `apply` also adds the ``af_median`` and ``af_stdev`` metadata to the corrected channels, representing the median and standard deviation of the measured blank distributions. Attributes ---------- channels : List(Str) The channels to correct. blank_file : File The filename of a file with "blank" cells (not fluorescent). Used to `estimate` the autofluorescence. blank_file_conditions : Dict Occasionally, you'll need to specify the experimental conditions that the blank tube was collected under (to apply the operations in the history.) Specify them here. Examples -------- Create a small experiment: .. plot:: :context: close-figs >>> import cytoflow as flow >>> import_op = flow.ImportOp() >>> import_op.tubes = [flow.Tube(file = "tasbe/rby.fcs")] >>> ex = import_op.apply() Create and parameterize the operation .. plot:: :context: close-figs >>> af_op = flow.AutofluorescenceOp() >>> af_op.channels = ["Pacific Blue-A", "FITC-A", "PE-Tx-Red-YG-A"] >>> af_op.blank_file = "tasbe/blank.fcs" Estimate the model parameters .. plot:: :context: close-figs >>> af_op.estimate(ex) Plot the diagnostic plot .. plot:: :context: close-figs >>> af_op.default_view().plot(ex) Apply the operation to the experiment .. plot:: :context: close-figs >>> ex2 = af_op.apply(ex) """ # traits id = Constant('edu.mit.synbio.cytoflow.operations.autofluorescence') friendly_id = Constant("Autofluorescence correction") name = Constant("Autofluorescence") channels = List(Str) blank_file = File(exists = True) blank_file_conditions = Dict({}) _af_median = Dict(Str, Float, transient = True) _af_stdev = Dict(Str, Float, transient = True) _af_histogram = Dict(Str, Tuple(Array, Array), transient = True)
[docs] def estimate(self, experiment, subset = None): """ Estimate the autofluorescence from `blank_file` in channels specified in `channels`. Parameters ---------- experiment : `Experiment` The experiment to which this operation is applied subset : Str (default = "") An expression that specifies the events used to compute the autofluorescence """ if experiment is None: raise util.CytoflowOpError('experiment', "No experiment specified") if not self.blank_file: raise util.CytoflowOpError('blank_file', "No blank tube specified") if not self.channels: raise util.CytoflowOpError('channels', "No channels specified") if not set(self.channels) <= set(experiment.channels): raise util.CytoflowOpError('channels', "Specified channels that weren't found " "in the experiment.") # don't have to validate that blank_file exists; should crap out on # trying to set a bad value # clear the current values self._af_median.clear() self._af_stdev.clear() self._af_histogram.clear() # make a little Experiment check_tube(self.blank_file, experiment) exp_conditions = {k: experiment.data[k].dtype.name for k in self.blank_file_conditions.keys()} blank_exp = ImportOp(tubes = [Tube(file = self.blank_file, conditions = self.blank_file_conditions)], conditions = exp_conditions, channels = {experiment.metadata[c]["fcs_name"] : c for c in experiment.channels}, name_metadata = experiment.metadata['name_metadata']).apply() # apply previous operations for op in experiment.history: if hasattr(op, 'by'): for by in op.by: if 'experiment' in experiment.metadata[by]: warn("Found experimental metadata {} in experiment history; " "make sure it's specified in blank_file_conditions." .format(by), util.CytoflowOpWarning) blank_exp = op.apply(blank_exp) # subset it if subset: try: blank_exp = blank_exp.query(subset) except Exception as exc: raise util.CytoflowOpError('subset', "Subset string '{0}' isn't valid" .format(subset)) from exc if len(blank_exp.data) == 0: raise util.CytoflowOpError('subset', "Subset string '{0}' returned no events" .format(subset)) af_median = {} af_stdev = {} for channel in self.channels: self._af_histogram[channel] = np.histogram(blank_exp[channel], bins = 250) channel_min = blank_exp[channel].quantile(0.025) channel_max = blank_exp[channel].quantile(0.975) blank_exp[channel] = blank_exp[channel].clip(channel_min, channel_max) af_median[channel] = np.median(blank_exp[channel]) af_stdev[channel] = np.std(blank_exp[channel]) # set atomically to support GUI self._af_stdev = af_stdev self._af_median = af_median
[docs] def apply(self, experiment): """ Applies the autofluorescence correction to channels in an experiment. Parameters ---------- experiment : `Experiment` the experiment to which this op is applied Returns ------- Experiment a new experiment with the autofluorescence median subtracted. The corrected channels have the following metadata added to them: - **af_median** : Float The median of the non-fluorescent distribution - **af_stdev** : Float The standard deviation of the non-fluorescent distribution """ if experiment is None: raise util.CytoflowOpError('experiment', "No experiment specified") if not self.channels: raise util.CytoflowOpError('channels', "No channels specified") if not self._af_median: raise util.CytoflowOpError(None, "Autofluorescence values aren't set. Did " "you forget to run estimate()?") if not set(self._af_median.keys()) <= set(experiment.channels) or \ not set(self._af_stdev.keys()) <= set(experiment.channels): raise util.CytoflowOpError(None, "Autofluorescence estimates aren't set, or are " "different than those in the experiment " "parameter. Did you forget to run estimate()?") if not set(self._af_median.keys()) == set(self._af_stdev.keys()): raise util.CytoflowOpError(None, "Median and stdev keys are different! " "What the hell happened?!") if not set(self.channels) == set(self._af_median.keys()): raise util.CytoflowOpError('channels', "Estimated channels differ from the channels " "parameter. Did you forget to (re)run estimate()?") new_experiment = experiment.clone(deep = True) for channel in self.channels: new_experiment[channel] = \ experiment[channel] - self._af_median[channel] new_experiment.metadata[channel]['af_median'] = self._af_median[channel] new_experiment.metadata[channel]['af_stdev'] = self._af_stdev[channel] new_experiment.history.append(self.clone_traits(transient = lambda t: True)) return new_experiment
[docs] def default_view(self, **kwargs): """ Returns a diagnostic plot to see if the autofluorescence estimation is working. Returns ------- IView An diagnostic view, call `AutofluorescenceDiagnosticView.plot` to see the diagnostic plots """ v = AutofluorescenceDiagnosticView(op = self) v.trait_set(**kwargs) return v
[docs]@provides(cytoflow.views.IView) class AutofluorescenceDiagnosticView(HasStrictTraits): """ Plots a histogram of each channel, and its median in red. Serves as a diagnostic for the autofluorescence correction. Attributes ---------- op : Instance(`AutofluorescenceOp`) The `AutofluorescenceOp` whose parameters we're viewing. Set automatically if you created the instance using `AutofluorescenceOp.default_view`. """ # traits id = Constant('edu.mit.synbio.cytoflow.view.autofluorescencediagnosticview') friendly_id = Constant("Autofluorescence Diagnostic") op = Instance(AutofluorescenceOp)
[docs] def plot(self, experiment, **kwargs): """ Plot a faceted histogram view of a channel """ if experiment is None: raise util.CytoflowViewError('experiment', "No experiment specified") if not self.op.channels: raise util.CytoflowViewError('op', "No channels specified") if not self.op._af_median: raise util.CytoflowViewError('op', "Autofluorescence values aren't set. Did " "you forget to run estimate()?") if not set(self.op._af_median.keys()) <= set(experiment.channels) or \ not set(self.op._af_stdev.keys()) <= set(experiment.channels) or \ not set(self.op._af_histogram.keys()) <= set(experiment.channels): raise util.CytoflowViewError('op', "Autofluorescence estimates aren't set, or are " "different than those in the experiment " "parameter. Did you forget to run estimate()?") if not set(self.op._af_median.keys()) == set(self.op._af_stdev.keys()): raise util.CytoflowOpError('op', "Median and stdev keys are different! " "What the heck happened?!") if not set(self.op.channels) == set(self.op._af_median.keys()): raise util.CytoflowOpError('channels', "Estimated channels differ from the channels " "parameter. Did you forget to (re)run estimate()?") import matplotlib.pyplot as plt plt.figure() for idx, channel in enumerate(self.op.channels): hist, bin_edges = self.op._af_histogram[channel] hist = hist[1:-1] bin_edges = bin_edges[1:-1] plt.subplot(len(self.op.channels), 1, idx+1) plt.title(channel) plt.bar(bin_edges[:-1], hist, width = bin_edges[2] - bin_edges[1], linewidth = 0) plt.axvline(self.op._af_median[channel], color = 'r') plt.tight_layout(pad = 0.8)