% NOTE -- ONLY EDIT THE .Rnw FILE!!! The .tex file is % likely to be overwritten. % %\VignetteIndexEntry{Basic Functions for Flow Cytometry Data} %\VignetteDepends{flowCore, flowViz} %\VignetteKeywords{} %\VignettePackage{flowCore} \documentclass[11pt]{article} \usepackage{times} \usepackage{hyperref} \usepackage[authoryear,round]{natbib} \usepackage{times} \usepackage{comment} \usepackage{graphicx} \usepackage{subfigure} \textwidth=6.2in \textheight=8.5in \oddsidemargin=.1in \evensidemargin=.1in \headheight=-.3in \newcommand{\scscst}{\scriptscriptstyle} \newcommand{\scst}{\scriptstyle} \newcommand{\Rfunction}[1]{{\texttt{#1}}} \newcommand{\Rcode}[1]{{\texttt{#1}}} \newcommand{\Robject}[1]{{\texttt{#1}}} \newcommand{\Rpackage}[1]{{\textsf{#1}}} \newcommand{\Rclass}[1]{{\textit{#1}}} \title{flowCore: data structures package for flow cytometry data} \author{N. Le Meur F. Hahne B. Ellis P. Haaland} \begin{document} \maketitle \begin{abstract} \noindent \textbf{Background} The recent application of modern automation technologies to staining and collecting flow cytometry (FCM) samples has led to many new challenges in data management and analysis. We limit our attention here to the associated problems in the analysis of the massive amounts of FCM data now being collected. From our viewpoint, see two related but substantially different problems arising. On the one hand, there is the problem of adapting existing software to apply standard methods to the increased volume of data. The second problem, which we intend to address here, is the absence of any research platform which bioinformaticians, computer scientists, and statisticians can use to develop novel methods that address both the volume and multidimensionality of the mounting tide of data. In our opinion, such a platform should be Open Source, be focused on visualization, support rapid prototyping, have a large existing base of users, and have demonstrated suitability for development of new methods. We believe that the Open Source statistical software R in conjunction with the Bioconductor Project fills all of these requirements. Consequently we have developed a Bioconductor package that we call \Rpackage{flowCore}. The \Rpackage{flowCore} package is not intended to be a complete analysis package for FCM data, rather, we see it as providing a clear object model and a collection of standard tools that enable R as an informatics research platform for flow cytometry. One of the important issues that we have addressed in the \Rpackage{flowCore} package is that of using a standardized representation that will insure compatibility with existing technologies for data analysis and will support collaboration and interoperability of new methods as they are developed. In order to do this, we have followed the current standardized descriptions of FCM data analysis as being developed under NIH Grant xxxx [n]. We believe that researchers will find flowCore to be a solid foundation for future development of new methods to attack the many interesting open research questions in FCM data analysis.\\ \noindent \textbf{Methods} We propose a variety different data structures. We have implemented the classes and methods in the Bioconductor package flowCore. We illustrate their use with X case studies.\\ \noindent \textbf{Results} We hope that those proposed data structures will be the base for the development of many tools for the analysis of high throughput flow cytometry.\\ \noindent \textbf{keywords} Flow cytometry, high throughput, software, standard \end{abstract} <>= library(flowCore) @ \section{Introduction} Traditionally, flow cytometry has been a tube-based technique limited to small-scale laboratory and clinical studies. High throughput methods for flow cytometry have recently been developed for drug discovery and advanced research methods \citep{Gasparetto2004}. As an example, the flow cytometry high content screening (FC-HCS) can process up to a thousand samples daily at a single workstation, and the results have been equivalent or superior to traditional manual multi-parameter staining and analysis techniques. The amount of information generated by high throughput technologies such as FC-HCS need to be transformed into executive summaries (which are brief enough) for creative studies by a human researcher \citep{Brazma2001}. Standardization is critical when developing new high throughput technologies and their associated information services \citep{Brazma2001, Chicurel2002, Boguski2003}. Standardization efforts have been made in clinical cell analysis by flow cytometry \citep{Keeney2004}, however data interpretation has not been standardized for even low throughput FCM. It is one of the most difficult and time consuming aspects of the entire analytical process as well as a primary source of variation in clinical tests, and investigators have traditionally relied on intuition rather than standardized statistical inference \citep{Bagwell2004,Braylan2004,Parks1997,Suni2003}. In the development of standards in high throughput FCM, few progress has been done in term of Open Source software. In this article we propose R data structures to handle flow cytometry data through the main steps of preprocessing: compensation, transformation, filtering. The aim is to merge both prada and rflowcyt \citep{LeMeur2006} into one core package which is compliant with the data exchange standards that are currently developed in the community \citep{Spidlen2006}. Visualization as well as quality control will than be part of the utility packages that depend on the data structures defined in the flowCore package. \section{Representing Flow Cytometry Data} \Rpackage{flowCore}'s primary task is the representation and basic manipulation of flow cytometry (or similar) data. This is accomplished through a data model very similar to that adopted by other Bioconductor packages using the \Robject{expressionSet} and \Robject{AnnotatedDataFrame} structures familiar to most Bioconductor users. \subsection{The \Rclass{flowFrame} Class} The basic unit of manipulation in \Rpackage{flowCore} is the \Rclass{flowFrame}, which corresponds roughly with a single ``FCS'' file exported from the flow cytometer's acquisition software. At the moment we support FCS file versions 2.0 through 3.0, and we expect to support FCS4/ACS1 as soon as the specification has been ratified. \subsubsection{Data elements} The primary elements of the \Robject{flowFrame} are the \Robject{exprs} and \Robject{parameters} slots, which contain the event-level information and column metadata respectively. The event information, stored as a single matrix, is accessed and manipulated via the \Rfunction{exprs()} and \Rfunction{exprs<-} methods, allowing \Robject{flowFrame}s to be stitched together if necessary (for example, if the same tube has been collected in two acquisition files for memory reasons). The \Robject{parameters} slot is an \Rclass{AnnotatedDataFrame} that contains information derived from an FCS file's ``\$P'' keywords, which describe the detector and stain information. The entire list is available via the \Rfunction{parameter()} method, but more commonly this information is accessed through the \Rfunction{names}, \Rfunction{featureNames} and \Rfunction{colnames} methods. The \Rfunction{names} function returns a concatenated version of \Rfunction{names} and \Rfunction{featureNames} using a format similar to the one employed by most flow cytometry analysis software. The \Rfunction{colnames} method returns the detector names, often named for the fluorochrome detected, while the \Rfunction{featureNames} methods returns the description field of the parameters, which will typically be an identifier for the antibody. The \Rfunction{keyword} method allows access to the raw FCS keywords, which are a mix of standard entries such as ``SAMPLE ID,'' vendor specific keywords and user-defined keywords that add more information about an experiment. In the case of plate-based experiments, there are also one or more keywords that identify the specific well on the plate. Most vendor software also include some sort of unique identifier for the file itself. The specialized methods \Rfunction{identifier} attempts to locate an appropriate globally unique identifier that can be used to uniquely identify a frame. Failing that, this method will return the original file name offering some assurance that this frame is at least unique to a particular session. \subsubsection{Reading a \Robject{flowFrame}} FCS files are read into the R environment via the \Rfunction{read.FCS} function using the standard connection interface---allowing for the possibility of accessing FCS files hosted on a remote resource as well as those that have been compressed or even retrieved as a blob from a database interface. FCS files (version 2.0 and 3.0) and LMD (List Mode Data) extensions are currently supported. There are also several immediate processing options available in this function, the most important of which is the \Rfunction{transformation} parameter, which can either ``linearize'' (the default) or ``linearize-with-PnG-scaling'' or ``scale'' our data. To see how this works, first we will examine an FCS file without any transformation at all: <>= file.name <- system.file("extdata","0877408774.B08", package="flowCore") x <- read.FCS(file.name, transformation=FALSE) summary(x) @ As we can see, in this case the values from each parameter seem to run from $0$ to $1023$ ($2^{10}-1$). However, inspection of the ``exponentiation'' keyword (\$PE) reveals that some of the parameters (3 and 4) have been stored in the format of the form $a\times{10^{x/R}}$ where $a$ is given by the first element of the string. <>= keyword(x,c("$P1E", "$P2E", "$P3E", "$P4E")) @ The default ``linearize'' transformation option will convert these to, effectively, have a ``\$PE'' of ``0,0'': <>= summary(read.FCS(file.name)) @ The ``linearize-with-PnG-scaling'' option will perform the previous transformation and it will also apply a ``division by gain'' to pamameters stored on linear scale with specified gain. The gain is specified in the \$PnG keywords. This option has been introduced as part of Gating-ML 2.0 compliance. Finally, the ``scale'' option will both linearize values as well as ensure that output values are contained in $[0,1]$, which is the proposed method of data storage for the ACS1.0/FCS4.0 specification: <>= summary(read.FCS(file.name,transformation="scale")) @ Another parameter of interest is the \Rfunction{alter.names} parameter, which will convert the parameter names into more ``R friendly'' equivalents, usually by replacing ``-'' with ``.'': <>= read.FCS(file.name,alter.names=TRUE) @ When only a particular subset of parameters is desired the \Rfunction{column.pattern} parameter allows for the specification of a regular expression and only parameters that match the regular expression will be included in the frame. For example, to include on the Height parameters: <>= x <- read.FCS(file.name, column.pattern="-H") x @ \textit{Note that \Rfunction{column.pattern} is applied after \Rfunction{alter.names} if it is used.} Finally, only a sample of lines can be read in case you need a quick overview of a large series of files. <>= lines <- sample(100:500, 50) y <- read.FCS(file.name, which.lines = lines) y @ \subsubsection{Visualizing a \Rclass{flowFrame}} Much of the more sophisticated visualization of \Rclass{flowFrame} and \Rclass{flowSet} objects, including an interface to the \Rpackage{lattice} graphics system is implemented by the \Rpackage{flowViz} package, also included as part of Bioconductor. Here, we will only introduce the standard \Rfunction{plot} function. The basic plot provides a simple pairs plot of all parameters: <>= library(flowViz) plot(x) @ To control the parameters being plotted we can supply a \Robject{y} value in the form of a character vector. If we choose exactly two parameters this will create a bivariate density plot. <>= plot(x,c("FL1-H", "FL2-H")) @ However, if we only supply a single parameter we instead get a univariate histogram, which also accepts the usual histogram arguments: <>= plot(x, "FL1-H", breaks=256) @ \subsection{The \Rclass{flowSet} Class} Most experiments consist of several \Rclass{flowFrame} objects, which are organized using a \Rclass{flowSet} object. This class provides a mechanism for efficiently hosting the \Rclass{flowFrame} objects with minimal copying, reducing memory requirements, as well as ensuring that experimental metadata stays properly to the appropriate \Rclass{flowFrame} objects. \subsubsection{Creating a \Rclass{flowSet}} To facilitate the creation of \Rclass{flowSet} objects from a variety of sources, we provide a means to coerce \Rclass{list} and \Rclass{environment} objects to a \Rclass{flowSet} object using the usual coercion mechanisms. For example, if we have a directory containing FCS files we can read in a list of those files and create a \Rclass{flowSet} out of them: <>= frames <- lapply(dir(system.file("extdata", "compdata", "data", package="flowCore"), full.names=TRUE), read.FCS) as(frames, "flowSet") @ Note that the original list is unnamed and that the resulting sample names are not particularly meaningful. If the list is named, the list constructed is much more meaningful. One such approach is to employ the \Rfunction{keyword} method for \Rclass{flowFrame} objects to extract the ``SAMPLE ID'' keyword from each frame: <>= names(frames) <- sapply(frames, keyword, "SAMPLE ID") fs <- as(frames, "flowSet") fs @ \subsubsection{Working with experimental metadata} Like most Bioconductor organizational classes, the \Rclass{flowSet} has an associated \Rclass{AnnotatedDataFrame} that provides metadata not contained within the \Rclass{flowFrame} objects themselves. This data frame is accessed and modified via the usual \Rfunction{phenoData} and \Rfunction{phenoData<-} methods. You can also generally treat the phenotypic data as a normal data frame to add new descriptive columns. For example, we might want to track the original filename of the frames from above in the phenotypic data for easier access: <>= phenoData(fs)$Filename <- fsApply(fs,keyword, "$FIL") pData(phenoData(fs)) @ Note that we have used the \Rclass{flowSet}-specific iterator, \Rfunction{fsApply}, which acts much like \Rfunction{sapply} or \Rfunction{lapply}. Additionally, we should also note that the \Rfunction{phenoData} data frame \textbf{must} have row names that correspond to the original names used to create the \Rclass{flowSet}. \subsubsection{Bringing it all together: \Rfunction{read.flowSet}} Much of the functionality described above has been packaged into the \Rfunction{read.flowSet} convenience function. In it's simplest incarnation, this function takes a \Rfunction{path}, that defaults to the current working directory, and an optional \Rfunction{pattern} argument that allows only a subset of files contained within the working directory to be selected. For example, to read a \Rclass{flowSet} of the files read in by \Robject{frame} above: <>= read.flowSet(path = system.file("extdata", "compdata", "data", package="flowCore")) @ \Rfunction{read.flowSet} will pass on additional arguments meant for the underlying \Rfunction{read.FCS} function, such as \Rfunction{alter.names} and \Rfunction{column.pattern}, but also supports several other interesting arguments for conducting initial processing: \begin{description} \item[files]{ An alternative to the \Rfunction{pattern} argument, you may also supply a vector of filenames to read.} \item[name.keyword]{ Like the example in the previous section, you may specify a particular keyword to use in place of the filename when creating the \Rclass{flowSet}.} \item[phenoData]{ If this is an \Rclass{AnnotatedDataFrame}, then this will be used in place of the data frame that is ordinarily created. Additionally, the row names of this object will be taken to be the filenames of the FCS files in the directory specified by \Robject{path}. This argument may also be a named list made up of a combination of \Robject{character} and \Robject{function} objects that specify a keyword to extract from the FCS file or a function to apply to each frame that will return a result. } \end{description} To recreate the \Rclass{flowSet} that we created by hand from the last section we can use \Rfunction{read.flowSet}s advanced functionality: <>= fs <- read.flowSet(path=system.file("extdata", "compdata", "data", package="flowCore"), name.keyword="SAMPLE ID", phenoData=list(name="SAMPLE ID", Filename="$FIL")) fs pData(phenoData(fs)) @ \subsubsection{Manipulating a \Rclass{flowSet}} You can extract a \Rclass{flowFrame} from a \Rclass{flowSet} object in the usual way using the \Rfunction{[[} or \Rfunction{\$} extraction operators. On the other hand using the \Rfunction{[} extraction operator returns a new \Rclass{flowSet} by \textbf{copying} the environment. However, simply assigning the \Rclass{flowFrame} to a new variable will \textbf{not} copy the contained frames. The primary iterator method for a \Rclass{flowSet} is the \Rfunction{fsApply} method, which works more-or-less like \Rfunction{sapply} or \Rfunction{lapply} with two extra options. The first argument, \Robject{simplify}, which defaults to \Robject{TRUE}, instructs \Rfunction{fsApply} to attempt to simplify it's results much in the same way as \Rfunction{sapply}. The primary difference is that if all of the return values of the iterator are \Rclass{flowFrame} objects, \Rfunction{fsApply} will create a new \Rclass{flowSet} object to hold them. The second argument, \Robject{use.exprs}, which defaults to \Robject{FALSE} instructs \Rfunction{fsApply} to pass the expression matrix of each frame rather than the \Rclass{flowFrame} object itself. This allows functions to operate directly on the intensity information without first having to extract it. As an aid to this sort of operation we also introduce the \Rfunction{each\_row} and \Rfunction{each\_col} convenience functions that take the place of \Rfunction{apply} in the \Rfunction{fsApply} call. For example, if we wanted the median value of each parameter of each \Rclass{flowFrame} we might write: <>= fsApply(fs, each_col, median) @ % which is equivalent to the less readable <>= fsApply(fs,function(x) apply(x, 2, median), use.exprs=TRUE) @ % In this case, the \Robject{use.exprs} argument is not required in the first case because \Rfunction{each\_col} and \Rfunction{each\_row} are methods and have been defined to work on \Rclass{flowFrame} objects by first extracting the intensity data. \section{Transformation} \Rpackage{flowCore} features two methods of transforming parameters within a \Rclass{flowFrame}: inline and out-of-line. The inline method, discussed in the next section has been developed primarily to support filtering features and is strictly more limited than the out-of-line transformation method, which uses R's \Rfunction{transform} function to accomplish the filtering. Like the normal \Rfunction{transform} function, the \Rclass{flowFrame}is considered to be a data frame with columns named for parameters of the FCS file. For example, if we wished to plot our first \Rclass{flowFrame}'s first two fluorescence parameters on the log scale we might write: <>= plot(transform(fs[[1]], `FL1-H`=log(`FL1-H`), `FL2-H`=log(`FL2-H`)), c("FL1-H","FL2-H")) @ Like the usual \Rfunction{transform} function, we can also create new parameters based on the old parameters, without destroying the old <>= plot(transform(fs[[1]], log.FL1.H=log(`FL1-H`), log.FL2.H=log(`FL2-H`)), c("log.FL1.H", "log.FL2.H")) @ \subsection{Standard Transforms} Though any function can be used as a transform in both the out-of-line and inline transformation techniques, \Rpackage{flowCore} provides a number of parameterized transform generators that correspond to the transforms commonly found in flow cytometry and defined in the Transformation Markup Language (Transformation-ML, see \url{http://www.ficcs.org/} and \cite{Spidlen2006} for more details). Briefly, the predefined transforms are: \begin{description} \item[truncateTransform]{ $y = \left\{\begin{array}{ll} a & x < a \\ x & x \geq a \\ \end{array}\right.$ } \item[scaleTransform]{$ f(x) = \frac{x-a}{b-a}$} \item[linearTransform]{ $ f(x) = a+bx$ } \item[quadraticTransform]{ $ f(x) = ax^2 + bx + c$ } \item[lnTransform]{$ f(x) = \log\left(x\right)\frac{r}{d} $} \item[logTransform]{$ f(x) = \log_b\left(x\right)\frac{r}{d} $} \item[biexponentialTransform]{ $ f^{-1}(x) = ae^{bx}-ce^{dx}+f$} \item[logicleTransform]{ A special form of the biexponential transform with parameters selected by the data.} \item[arcsinhTransform]{ $ f(x) = asinh\left(a+bx\right)+c$} \end{description} To use a standard transform, first we create a transform function via the constructors supplied by \Rpackage{flowCore}: <>= aTrans <- truncateTransform("truncate at 1", a=1) aTrans @ which we can then use on the parameter of interest in the usual way <>= transform(fs,`FL1-H`=aTrans(`FL1-H`)) @ \section{Gating} The most common task in the analysis of flow cytometry data is some form of filtering operation, also known as gating, either to obtain summary statistics about the number of events that meet a certain criteria or to perform further analysis on a subset of the data. Most filtering operations are a composition of one or more common filtering operations. The definition of gates in \Rpackage{flowCore} follows the Gating Markup Language Candidate Recommendation \cite{Spidlen2008}, thus any \Rpackage{flowCore} gating strategy can be reproduced by any other software that also adheres to the standard and \textit{vice versa}. \subsection{Standard gates and filters} Like transformations, \Rpackage{flowCore} includes a number of built-in common flow cytometry gates. The simplest of these gates are the geometric gates, which correspond to those typically found in interactive flow cytometry software: \begin{description} \item[rectangleGate]{Describes a cubic shape in one or more dimensions--a rectangle in one dimension is simply an interval gate.} \item[polygonGate]{Describes an arbitrary two dimensional polygonal gate.} \item[polytopeGate]{Describes a region that is the convex hull of the given points. This gate can exist in dimensions higher than $2$, unlike the \Robject{polygonGate}.} \item[ellipsoidGate]{Describes an ellipsoidal region in two or more dimensions} \end{description} These gates are all described in more or less the same manner (see man pages for more details): <>= rectGate <- rectangleGate(filterId="Fluorescence Region", "FL1-H"=c(0, 12), "FL2-H"=c(0, 12)) @ In addition, we introduce the notion of data-driven gates, or filters, not usually found in flow cytometry software. In these approaches, the necessary parameters are computed based on the properties of the underlying data, for instance by modelling data distribution or by density estimation : \begin{description} \item[norm2Filter]{ A robust method for finding a region that most resembles a bivariate Normal distribution. } \item[kmeansFilter]{ Identifies populations based on a one dimensional k-means clustering operation. Allows the specification of \textbf{multiple} populations. } \end{description} \subsection{Count Statistics} When we have constructed a filter, we can apply it in two basic ways. The first is to collect simple summary statistics on the number and proportion of events considered to be contained within the gate or filter. This is done using the \Rfunction{filter} method. The first step is to apply our filter to some data <>= result = filter(fs[[1]],rectGate) result @ As we can see, we have returned a \Rclass{filterResult} object, which is in turn a filter allowing for reuse in, for example, subsetting operations. To obtain count and proportion statistics, we take the summary of this \Rclass{filterResult}, which returns a list of summary values: <>= summary(result) summary(result)$n summary(result)$true summary(result)$p @ A filter which contains multiple populations, such as the \Rclass{kmeansFilter}, can return a list of summary lists: <>= summary(filter(fs[[1]], kmeansFilter("FSC-H"=c("Low", "Medium", "High"), filterId="myKMeans"))) @ A filter may also be applied to an entire \Rclass{flowSet}, in which case it returns a list of \Rclass{filterResult} objects: <>= filter(fs,rectGate) @ \subsection{Subsetting} To subset or split a \Rclass{flowFrame} or \Rclass{flowSet}, we use the \Rfunction{Subset} and \Rfunction{split} methods respectively. The first, \Rfunction{Subset}, behaves similarly to the standard R \Rfunction{subset} function, which unfortunately could not used. For example, recall from our initial plots of this data that the morphology parameters, Forward Scatter and Side Scatter contain a large more-or-less ellipse shaped population. If we wished to deal only with that population, we might use \Rfunction{Subset} along with a \Rclass{norm2Filter} object as follows: <>= morphGate <- norm2Filter("FSC-H", "SSC-H", filterId="MorphologyGate", scale=2) smaller <- Subset(fs, morphGate) fs[[1]] smaller[[1]] @ Notice how the \Robject{smaller} \Rclass{flowFrame} objects contain fewer events. Now imagine we wanted to use a \Rclass{kmeansFilter} as before to split our first fluorescence parameter into three populations. To do this we employ the \Rfunction{split} function: <>= split(smaller[[1]], kmeansFilter("FSC-H"=c("Low","Medium","High"), filterId="myKMeans")) @ or for an entire \Rclass{flowSet} <>= split(smaller, kmeansFilter("FSC-H"=c("Low", "Medium", "High"), filterId="myKMeans")) @ \subsection{Combining Filters} Of course, most filtering operations consist of more than one gate. To combine gates and filters we use the standard R Boolean operators: \&, | and ! to construct an intersection, union and complement respectively: <>= rectGate & morphGate rectGate | morphGate !morphGate @ we also introduce the notion of the subset operation, denoted by either \%subset\% or \%\&\%. This combination of two gates first performs a subsetting operation on the input \Rclass{flowFrame} using the right-hand filter and then applies the left-hand filter. For example, <>= summary(filter(smaller[[1]],rectGate %&% morphGate)) @ first calculates a subset based on the \Robject{morphGate} filter and then applies the \Robject{rectGate}. \subsection{Transformation Filters} Finally, it is sometimes desirable to construct a filter with respect to transformed parameters. To allow for this in our filtering constructs we introduce a special form of the \Rfunction{transform} method along with another filter combination operator \%on\%, which can be applied to both filters and \Rclass{flowFrame} or \Rclass{flowSet} objects. To specify our transform filter we must first construct a transform list using a simplified version of the \Rfunction{transform} function: <>= tFilter <- transform("FL1-H"=log,"FL2-H"=log) tFilter @ Note that this version of the transform filter does not take parameters on the right-hand side--the functions can only take a single vector that is specified by the parameter on the left-hand side. In this case those parameters are ``FL1-H'' and ``FL2-H.'' The function also does not take a specific \Rclass{flowFrame} or \Rclass{flowSet} allowing us to use this with any appropriate data. We can then construct a filter with respect to the transform as follows: <>= rect2 <- rectangleGate(filterId="Another Rect", "FL1-H"=c(1,2), "FL2-H"=c(2,3)) %on% tFilter rect2 @ Additionally, we can use this construct directly on a \Rclass{flowFrame} or \Rclass{flowSet} by moving the transform to the left-hand side and placing the data on the right-hand side: <>= plot(tFilter %on% smaller[[1]],c("FL1-H","FL2-H")) @ which has the same effect as the $\log$ transform used earlier. \section{filterSet: Organizing Filtering Strategies} \subsection{Building Filter Sets} In addition to allowing for sequential filtering, \Rpackage{flowCore} also provides a \Rclass{filterSet} object that serves as an analogue to \Rclass{flowSets} for filters. The primary use of this object is to organize and manipulate complex gating strategies. For example, recall from the filtering section we had two gates, a ``morphologyGate'' and a ``Fluorescence Region'' gate that we used to demonstrate the various logical operators available for gates. To recreate this example as a \Rclass{filterSet}, we would do the following: <>= fset1 = filterSet( rectangleGate(filterId="Fluorescence Region","FL1-H"=c(50,100), "FL2-H"=c(50,100)), norm2Filter("FSC-H","SSC-H",filterId="Morphology Gate", scale=2), ~ `Fluorescence Region` & `Morphology Gate`, ~ `Fluorescence Region` | `Morphology Gate`, Debris ~ ! `Morphology Gate`, ~ `Fluorescence Region` %&% `Morphology Gate` ) fset1 @ There are two features of note in \Rfunction{filterSet}, which can also take a single list argument for programmatic creation of \Rclass{filterSet} objects. The first is that there is a formula interface for the creation of \Rclass{filter} objects from operators. The formula interface can be used with or without a left-hand side, which specifies an optional filter identifier. If the filter identifier is not specified, as is the case with all but the ``Debris'' gate in the example above the usual default filter identifier is constructed. Non-formula gates can also have an optional name specified which overrides the filter identifier specified in the gate constructor. This is discouraged at this time as it leads to confusion when creating new filters. \subsection{Manipulating Filter Sets} Manipulating a \Rclass{filterSet} is done using the normal list element replacement methods, though we do provide a special case where replacing the ``""'' element or the ``NULL'' element causes the \Rclass{filterSet} to use the ``filterId'' slot of the incoming filter to choose a name. Similarly, specifying an identifier will override any ``filterId'' slot in the original filter. Additionally, there are several convenience functions and methods available for manipulating a \Rclass{filterSet}. The \Rfunction{names} function provides a list of filter identifiers in an arbitrary order. The \Rfunction{sort} lists the filter identifiers according to a topological sort of the filter dependencies such that parent filters are always evaluated before their children. Filters without dependencies are provided in an arbitrary relative ordering. To obtain the adjacency matrix for this sorting, the option ``dependencies=TRUE'' supplied to the \Rfunction{sort} function will attach an ``AdjM'' attribute to the resulting character vector. \subsection{Using Filter Sets} Though they are not explicitly filters, meaning that they cannot play a role in the construction of sub filters via operators, \Rclass{filterSet} objects do define \Rfunction{filter} and \Rfunction{split} methods that allow for the two most common use cases. Additionally, a filter can be selected from the \Rclass{filterSet} using ``[['' to perform a \Rfunction{Subset} operation if desired. For standard filtering operations, the \Rclass{filterSet} can yield better performance than a manual filtering operation because it ensures that each parent filter is only evaluated once per call to \Rfunction{filter} so that if many sub filters rely on a particular, likely expensive, filtering operation they use the original result rather than recalculating each time. The individual \Rclass{filterResult} objects may be later extracted for further analysis of desire, making the \Rclass{filterSet} the logical choice for most filtering operations. Additionally, though not presently implemented, it can allow visualization methods to make more intelligent choices about display when rendering gating strategies because the disposition of the entire filtering graph is known and the filter results are already available. Like any other filtering operation, \Rfunction{summary} methods are also available to collect the usual sorts of count statistics for a filtering operation: <>= f = filter(fs,fset1) f @ When splitting a \Rclass{flowFrame} using \Rfunction{split}, we have two options. The first is to include all of the resultant subsets in the usual way: <>= split(fs[[1]], fset1, flowSet=TRUE) @ While this provides access to all subsets, for automated analysis we will often only wish to work with the ``leaf'' gates (i.e. those without children of their own). In this case, we can specify ``drop=TRUE'' to obtain only the leaf gates: <>= split(fs[[1]],fset1,drop=TRUE,flowSet=TRUE) @ Note that in both cases, we use the ``flowSet=TRUE'' option that has been added to the \Rfunction{split} function to automatically create a \Rclass{flowSet} object from the resulting \Rclass{flowFrame} objects. For \Rfunction{filterSet} objects, the \Rfunction{split} function can take advantage of the added structure in this case to provide some information about the filtering operation itself in the \Rclass{flowSet}'s metadata as well. \section{Work flows} \Rclass{filterSets} are very limited in their use for complex analysis work flows. They are result-centric and it is hard to access intermediate results. \Rpackage{flowCore} offers much more versatile tools for such tasks though the \Rclass{workFlow} class. The general idea is to let the software handle the organization of intermediate results, naming schemes and operations and to provide a unified API to access and summarize these operations. \subsection{Abstraction of work flows} There are three classes in \Rpackage{flowCore} that are used to abstract work flows: \Rclass{workFlow} objects are the basic container holding all the necessary bits and pieces and they are the main structure for user interaction. \Rclass{actionItem} objects are abstractions of data analysis operations like gating, transformation or normalization. There are sub-classes for different types of operations. Finally, \Rclass{view} objects are created when applying \Rclass{actionItems} to a \Rclass{workFlow}. One can think of \Rclass{views} as separate data sets, for instance a subset that is created by a gating operation or the modified data following a transformation. Again, different types of \Rclass{views} are available through various sub-classes. This structure allows for a unified API in which the user interacts with either \Rclass{views} (the most common application) or \Rclass{actionItems}. It is important to know that \Rclass{workFlows} use a reference semantic instead of the pass-by-value semantic that is usually found in the R language. This design was chosen to minimize memory consumption and to keep the framework as flexible as possible. The only main consequence on the user-level is the fact that direct assignments to a \Rclass{workFlow} object are usually not necessary; i.e., functions that operate on the \Rclass{workFlow} have the potential side-effect of modifying the object. \subsection{Naming schemes} One major obstacle of using scripts as a work flow representation is the need to keep track of object symbols or names of list items. Choosing arbitrary names makes it hard to access objects and there is no way to easily query the structure of the work flow other than reading through the complete script. It is also very hard to appreciate the hierarchical structure of the work flow from a linear listing of R commands. Most objects in \Rpackage{flowCore} provide useful names or identifiers, and the work flow framework tries to utilize these names as much as possible. For example, a \Rclass{rectangleGate} called ``Lymphocytes'' will create two subpopulations, namely the cells within the gate and their complement. Intuitive names for these subpopulation would therefore be ``Lymphocyte+'' and ``Lymphocyte-''. These names don't necessarily need to be unique, and in order to be able to unambiguously address each item in a work flow we need additional unique identifiers. Our software keeps an alias table that maps these internal unique identifiers to more human-readable and user friendly names, and unless there are naming collisions, those names will also be the primary identifiers that are exposed in the API. Only if and alias is not unique, objects need to be addressed by their internal unique identifier. \subsection{Creating \Rclass{workFlow} objects} Objects of class \Rclass{workFlow} have to be created using the constructor \Rfunction{workFlow}. The single mandatory argument is a flow data object, either a \Rclass{flowFrame} or a \Rclass{flowSet}. <>= data(GvHD) wf <- workFlow(GvHD[1:5], name="myWorkflow") wf @ % Printing the objects shows its structure: the work flow contains a single \Rclass{view} ``base view'' which essentially is the \Rclass{flowSet} \Robject{fs}. Since the data was added by the constructor and is not the result of a particular data analysis operation, no \Rclass{actionItem} is associated to the \Rclass{view}. We can query the available views in the \Rclass{workFlow} using the \Rfunction{views} method: <>= views(wf) @ \subsection{Adding \Rclass{actionItems}} Currently there are four different types of \Rclass{actionItems} that can be added to a \Rclass{workFlow}: \begin{itemize} \item{\Rclass{compensateActionItem}: A \Rclass{compensation} object} \item{\Rclass{transformActionItem}: A \Rclass{transformList} object} \item{\Rclass{normalizeActionItem}: A \Rclass{normalization} object} \item{\Rclass{gateActionItem}: A \Rclass{filter} object} \end{itemize} The user doesn't have to bother with the details of these classes, they will be created internally when adding one of the above source objects to the \Rclass{workFlow} using the \Rfunction{add} method. <>= tf <- transformList(colnames(GvHD[[1]])[3:6], asinh, transformationId="asinh") add(wf, tf) wf views(wf) @ % There are several things to note here. First, we didn't assign the return value of \Rfunction{add} back to \Robject{wf}. This is the aforementioned feature of the reference semantic. Internally, \Rclass{workFlows} are to the most part \Rclass{environments}, and all modifications change the content of these \Rclass{environments} rather than the object itself. We also used the \Rfunction{add} methods with two arguments only: the \Rclass{workFlow} and the \Rclass{transformList}. In a hierarchical work flow we usually want to control to which subset (or rather to which \Rclass{view}, to stick to the work flow lingo) we want to add the \Rclass{actionItem} to. The default is to use the base view, and alternative parent \Rclass{views} have to be specified using the \Robject{parent} argument; \Robject{add(wf,tf,parent="base view")} would have been equivalent. Adding the transformation created a new \Rclass{view} called ``asinh'' and we can see that this name was automatically taken from the \Rclass{transformationList} object. If we want to control the name of the new \Rclass{view}, we can do so using the \Robject{name} argument: <>= add(wf, tf, name="another asinh transform") wf @ % This operation created a warning since we didn't use a leaf node to add the transformation. The software can't be sure if any of the downstream leafs in the work flow hierarchy are affected by the transformation (as it would certainly be true for most gates). One solution would be to update all downstream leafs in the work flow, but this feature is not supported yet. Also note that we used the same \Rclass{transformList} object twice. While we controlled the name that was used for the second \Rclass{view}, the alias for the \Rclass{actionItem} is the same in both cases, and we would have to use the internal unique identifier to unambiguously access the respective object. Transformation operations (and also compensations or normalizations) only create single \Rclass{views}. This is not true for gating operations. Even the most simple gate creates two populations: those cells in the gate and the complement. Accordingly, two separate \Rclass{views} have to be created as well. <>= rg <- rectangleGate("FSC-H"=c(200,400), "SSC-H"=c(250, 400), filterId="rectangle") add(wf, rg, parent="asinh") wf @ % As we see, there are now two new \Rclass{views} under the ``asinh'' node called ``rectangle+'' and ``rectangle-''. Other filter types may create more than two populations, for instance a \Rclass{quadrantGate} always results in four sub-populations. There is no restriction on the number of populations that are allowed to be created by a single gating operation, however, when gating a \Rclass{flowSet}, the number of populations have to be the same for each frame in the set. <>= qg <- quadGate("FL1-H"=2, "FL2-H"=4) add(wf,qg,parent="rectangle+") wf @ % For complex work flows it becomes increasingly hard to convey all information by printing on the screen. There is a plotting function for \Rclass{workFlow} objects which plots the underlying tree. The function depends on the \Rpackage{Rgraphviz} package. <>= plot(wf) @ <>= if(suppressWarnings(require(Rgraphviz))){ plot(wf) }else{ plot(1,1, type="n", axes=FALSE, ann=FALSE) text(1,1,"Need to install Rgraphviz") } @ \subsection{Accessing items in the \Rclass{workFlow} object} The main reason for having the \Rclass{workFlow} class is to be able to easily access all elements of a potentially complex analysis work flow. The predominant use case here is to access individual \Rclass{views}, either for plotting, to create summary statistics or to use the underlying data for subsequent analysis steps. You can do that using the familiar list subsetting syntax: <>= wf[["rectangle+"]] wf$asinh @ % This also works for \Rclass{actionItems}: <>= wf[["action_rectangle"]] @ % In order to retreive the underlying data for a \Rclass{view} you can use the \Rfunction{Data} method. <>= Data(wf[["rectangle-"]]) @ % There are \Rfunction{summary} methods defined for the different \Rclass{view} subclasses, which basically create the respective summaries of the associated \Rpackage{flowCore} objects, e.g. a \Rclass{filterResult} for a filtering operation. <>= summary(wf[["action_rectangle"]]) summary(wf[["CD15 FITC+CD45 PE+"]]) @ The \Rpackage{flowViz} package also defines \Rfunction{xyplot} and \Rfunction{densityplot} methods for \Rclass{view} objects. Sensible defaults are usually chosen automatically, for instance to add the gate boundaries to the plot of a \Rclass{gateView}. <>= densityplot(wf[["base view"]]) xyplot(`FL1-H` ~ `FL2-H`, wf[["CD15 FITC+CD45 PE+"]]) @ \subsection{Removing items from \Rclass{workFlow} object} The hierarchical structur of the \Rclass{workFlow} object introduces dependencies between \Rclass{views} or between \Rclass{views} and \Rclass{actionItems}. Thus, removing a particular view means also removing all of its associated child \Rclass{views} and \Rclass{actionItems}. The easiest way to remove objects from a \Rclass{workFlow} is through the \Robject{undo} mechanism. The function takes a \Rclass{workFlow} object as first argument and an optional second integer argument giving the number of operartions that are supposed to be roled back. <>= undo(wf) wf @ % Sometimes it is more convenient to specify a particular \Rclass{view} or \Rclass{action} item that is to be removed. This can be archived using the \Rfunction{Rm} methods for \Rclass{views} or \Rclass{actionItems}. As a side-effect these methods will also remove all dependent items and the user should be well aware of this potentially dangerous behaviour. <>= Rm('rectangle-', wf) @ <>= Rm('asinh', wf) @ \clearpage \bibliographystyle{plainnat} \bibliography{cytoref} \end{document}