---
title: "An Introduction to the openCyto package"
output: 
  html_document: 
    number_sections: false
    toc: true
    toc_float: true
    collapsed: false
vignette: >
  %\VignetteEngine{knitr}
  %\VignetteIndexEntry{An Introduction to the openCyto package}
---

```{r requirements, echo=FALSE}
if (!require(flowWorkspaceData)) {
  stop("Cannot build the vignettes without 'flowWorkspaceData'")
}
```


```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(out.extra='style="display:block; margin: auto"', fig.align="center", message = FALSE, warning = FALSE)
```


1. Introduction
------------


The  **openCyto** package is designed to facilitate the application of automated gating methods in a sequential way to mimic the construction of a manual gating scheme.


### 1.1. Manual gating 

Traditionally, scientists have to draw the gates for each individual sample on each 2-D projection (2 channels) within `flowJo`. Alternatively, they can draw template gates on one sample and replicate them to other samples, then manually inspect the gate on each sample
to do the correction if necessary. Either way is time consuming and subjective, thus not suitable for the large data sets
generated by high-throughput flow cytometry, CyTOF, or "cross-lab" data analysis.

Here is one `xml` workspace (manual gating scheme) exported from `flowJo`.
```{r load-flowWorkspace, echo=F}
library(flowWorkspace)
```
```{r load-xml, eval=TRUE}
flowDataPath <- system.file("extdata", package = "flowWorkspaceData")
wsfile <- list.files(flowDataPath, pattern="manual.xml",full = TRUE)
wsfile
```
By using the `CytoML` package, We can load it into R,
```{r open_flowjo_xml, eval=F}
library(CytoML)
ws <- open_flowjo_xml(wsfile)
```
apply the`manual gates`defined in`xml`to the raw`FSC`files,
```{r flowjo_to_gatingset, eval=F}
gs <- flowjo_to_gatingset(ws, name= "T-cell", subset =1)
```
```{r load_gs_manual, echo = FALSE}
gs <- load_gs(file.path(flowDataPath,"gs_manual"))
```

and then visualize the`Gating Hierarchy` 
```{r plot-manual-GatingHierarchy}
gh <- gs[[1]]
plot(gh)
```
and the`gates`: 
```{r plot-manual-gates, fig.width = 9}
library(ggcyto)
autoplot(gh)
```   
This is a gating scheme for a `T cell` panel, which tries to identify `T cell` sub-populations.
We can achieve the same results by using the automated gating pipeline provided by this package.

### 1.2. Automated Gating
`flowCore`,`flowStats`,`flowClust` and other packages provide many different gating methods to 
detect cell populations and draw gates automatically. 

The `flowWorkspace` package provides the `GatingSet` as an efficient data structure to store, query and visualize the hierarchical gated data.

By taking advantage of these tools, the `openCyto` package can create the automated gating pipeline by a `gatingTemplate`, which is essentially the same kind of hierarchical gating scheme 
used by scientists.
  

2. Create gating templates
-----------------------------


### 2.1. Template format
First of all, we need to describe the gating hierarchy in a spread sheet (a plain text format).
This spread sheet must have the following columns:

* `alias`: a name used to label the cell population, with the path composed of the alias and its precedent nodes (e.g. /root/A/B/alias) being uniquely identifiable.
* `pop`: population patterns of `+/-` or `+/-+/-`, which tell the algorithm which side (postive or negative) of a 1-D gate or which quadrant of a 2-D gate are to be kept.
* `parent`: the parent population alias, whose path also has to be uniquely identifiable.
* `dims`: characters seperated by commas specifying the dimensions (1-D or 2-D) used for gating. These can be either channel names or stained marker names.
* `gating_method`: the name of the gating function (e.g. `flowClust`). It is invoked by a wrapper function that has the identical function name prefixed with a dot.(e.g. `.flowClust`)
* `gating_args`: the named arguments passed to the gating function
* `collapseDataForGating`: When TRUE, data is collapsed (within groups if `groupBy` is specified) before gating and the gate is replicated across collapsed samples.
  When set FALSE (or blank), the `groupBy` argument is only used by `preprocessing` and ignored by gating.
* `groupBy`: If provided, samples are split into groups by the unique combinations of the named study variable (i.e. column names of pData, e.g."PTID:VISITNO").
  When this is numeric (N), samples are grouped by every N samples 
* `preprocessing_method`: the name of the preprocessing function (e.g. `prior_flowClust`). It is invoked by a wrapper function that has the identical function name prefixed with a dot (e.g. `.prior_flowClust`).
  The preprocessing results are then passed to the appropriate gating wrapper function through its `pps_res` argument.
* `preprocessing_args`: the named arguments passed to the preprocessing function.

### 2.2. Example template
Here is an example of a gating template.
```{r gatingTemplate, eval = T}
library(openCyto)
library(data.table)
gtFile <- system.file("extdata/gating_template/tcell.csv", package = "openCyto")
dtTemplate <- fread(gtFile)
dtTemplate
```
Each row is usually corresponding to one cell population and the gating method that is used to get that population.
We will try to explain how to create this gating template based on the manual gating scheme row by row.

#### 2.2.1. "nonDebris"
```{r gatingTemplate-nonDebris, eval = T}
dtTemplate[1,]
```  
* The population name is `"nonDebris"` (specified in the `alias` field).
* The `parent` node is `root` (which is always the first node of a `GatingHierarchy` by default). 
* We use `mindensity` (one of the `gating` functions provided by `openCyto` package) as the `gating_method` to gate on dimension (`dim`) of `FSC-A`.
* As a result, it will generate a 1-D gate on `FSC-A`. The `+` in the `pop` field indicates the
`positive` side of the 1-D gate is kept as the population of interest. 
* There is no `grouping` or `preprocessing` involved in this gate, so the other columns are left blank.

#### 2.2.2. "singlets" 
```{r gatingTemplate-singlets, eval = T}
dtTemplate[2,]
```
* The population name is `"singlets"` (the `alias` field).
* The `parent` node is `nonDebris`.
* The `gating_method` is `singletGate` (a function from the `flowStats` package)
* As a result, a `polygonGate` will be generated on `FSC-A` and `FSC-H` (specified by `dims`) for each sample.
* Again, the `+` in the `pop` field stands for `"singlets+"`. But here it is 2-D gate, which means we want to keep the area 
inside of the polygon. 

#### 2.2.3. "lymphocytes" 
```{r gatingTemplate-lympth, eval = T}
dtTemplate[3,]
```
* Similarly, `alias` specifies the name of population.
* `parent` points to `singlets`
* Since we are going to use `flowClust` as `gating_method` to do the 2-dimensional gating,
	`dims` is a comma-separated string: `x` axis (`FSC-A`) goes first, `y` (`SSC-A`) the second. 
	This order doesn't affect the gating process but will determine how the gates are displayed.     
* All the parameters that `flowClust` algorithm accepts can be put in `gating_args` as if they are typed in the `R console`.
	see `help(flowClust)` for more details of these arguments
* The `flowClust` algorithm accepts the extra argument `prior` that is calculated during the `preprocessing` stage (before the actual gating). Thus, we supply the `preprocessing_method` with `prior_flowClust`.

#### 2.2.4. "cd3+" (Tcells) 

```{r gatingTemplate-cd3, eval = T}
dtTemplate[4,]
```
This is similar to the `nonDebris` gate except that we specify `collapseDataForGating` as `TRUE`, 
which tells the pipeline to `collapse` all samples into one and apply `mindensity` to the collapsed data on `CD3` dimension.
Once the gate is generated, it is replicated across all samples. This is only useful when each individual sample does not have
enough events to deduce the gate. Here we do this just for the purpose of proof of concept. 

#### 2.2.5. CD4 and CD8
The fifth row specifies `pop` as `cd4+/-cd8+/-`, which will be expanded into 6 rows. 
```{r gatingTemplate-cd4cd8, eval = T}
dtTemplate[5,]
```
```{r gatingTemplate-expand, echo = F, results = F}
expanded <- openCyto:::.preprocess_csv(dtTemplate)
rownames(expanded) <- NULL
```

The first two rows are two 1-D gates that will be generated by `gating_method` on each 
dimension (`cd4` and `cd8`) independently:
```{r gatingTemplate-expand1, echo = F}
expanded[5:6,]
```

Then another 4 rows are 4 `rectangleGate`s that corresponds to the 4 `quadrants` in the 2-D projection (`cd4 vs cd8`).
```{r gatingTemplate-expand2, echo = F}
expanded[7:10,]
```
As we see here, `"refGate"` in `gating_method` indicates that they are constructed based on the 
`gate coordinates` of the previous two 1-D gates.
Those 1-D gates are thus considered as "reference gates" that are referred to by a colon-separated `alias` string in `gating_args`: `"cd4+:cd8+"`.

Alternatively, we can expand it into these 6 rows explicitly in the spreadsheet. 
But this convenient representation is recommended unless the user wants to have finer control on how the gating is done. 
For instance, sometimes we need to use different `gating_method`s to generate 1-D gates on `cd4` and `cd8`.
Or it could be the case that `cd8` gating needs to depend on `cd4` gating, i.e. the `parent` of `cd8+` is `cd4+`(or `cd4-`) instead of `cd3`.
Sometimes we want to have a customized `alias` other than the quadrant-like name (`x+y+`) that gets generated automatically. 
(e.g. 5th row of the gating template)

3. Load gating template
----------------------------- 
After the gating template is defined in the spreadsheet, it can be loaded into R:  	 	
```{r load-gt, eval = T}
gt_tcell <- gatingTemplate(gtFile)
gt_tcell
```
Besides looking at the spreadsheet, we can examine the gating scheme by visualizing it:
```{r plot-gt, eval = T}
plot(gt_tcell)
```
As we can see, the gating scheme has been expanded as we described above.
All the **colored** arrows source from a `parent` population and the **grey** arrows source from a `reference` population(/gate).  

4. Run the gating pipeline
-----------------------------
Once we are satisfied with the gating template, we can apply it to the actual flow data.

### 4.1. Load the raw data
First of all, we load the raw FCS files into R by `ncdfFlow::read.ncdfFlowSet` (it uses less memory than `flowCore::read.flowSet`) and create an empty `GatingSet` object.
```{r load-fcs}
fcsFiles <- list.files(pattern = "CytoTrol", flowDataPath, full = TRUE)
cs  <- load_cytoset_from_fcs(fcsFiles)
cf <- realize_view(cs[[1]])
gs <- GatingSet(cs)
gs
```
### 4.2. Compensation
Then, we compensate the data. If we have compensation controls (i.e. singly stained samples), we can calculate the
compensation matrix by using the `flowStats::spillover` function.
Here we simply use the compensation matrix defined in the `flowJo` workspace.
```{r compensate}
compMat <- gh_get_compensations(gh)
compensate(gs, compMat)
```

Here is one example showing the compensation outcome:
```{r compensate_plot, echo = F, fig.width = 4, fig.height = 4}
sub_chnl <- c("V545-A","V450-A")
cf <- cf[,sub_chnl]
cf_comp <- realize_view(gh_pop_get_data(gs[[1]])[,sub_chnl])
cs <- cytoset(list(cf = cf, cf_comp = cf_comp))
#transform data to better visualize the compensation effect
transform(cs, estimateLogicle(cf,sub_chnl))
gridExtra::grid.arrange(as.ggplot(autoplot(cs[[1]], "V545", "V450")), as.ggplot(autoplot(cs[[2]], "V545", "V450")), nrow = 1)
# Need to grab the post-compensation, pre-transformation cytoframe before next steps
cf_comp <- realize_view(gh_pop_get_data(gs[[1]])[,sub_chnl])
```

### 4.3. Transformation
All of the stained channels need to be transformed properly before the gating.
Here we use the `flowCore::estimateLogicle` method to determine the `logicle` transformation.
```{r transformation, eval = T}
chnls <- parameters(compMat)
trans <- estimateLogicle(gs[[1]], channels = chnls)
gs <- transform(gs, trans)
```
Here is one example showing the transformation outcome:
```{r transformation_plot, echo = F, fig.width = 5, fig.height = 5}
cf_trans <- gh_pop_get_data(gs[[1]])[,sub_chnl[1]]
cf_comp <- cf_comp[,sub_chnl[1]]
p1 <- as.ggplot(autoplot(cf_comp, "V545"))
p2 <- as.ggplot(autoplot(cf_trans, "V545"))
plot(gridExtra::arrangeGrob(p1,p2))
```

### 4.5. Gating
Now we can apply the gating template to the data:
```{r gating, eval = TRUE}
gt_gating(gt_tcell, gs)
```
Optionally, we can run the pipeline in parallel to speed up gating. e.g.
```{r gating_par, eval = FALSE}
gt_gating(gt_tcell, gs, mc.cores=2, parallel_type = "multicore")
```

### 4.6. Hide nodes
After gating, there are some extra populations generated automatically by the pipeline (e.g. `refGate`).
```{r plot_afterGating}
plot(gs[[1]])
```
We can hide these populations if we are not interested in them: 
```{r hideGate, results = "hide"}
nodesToHide <- c("cd8+", "cd4+"
				, "cd4-cd8-", "cd4+cd8+"
				, "cd4+cd8-/HLA+", "cd4+cd8-/CD38+"
				, "cd4-cd8+/HLA+", "cd4-cd8+/CD38+"
				, "CD45_neg/CCR7_gate", "cd4+cd8-/CD45_neg"
				, "cd4-cd8+/CCR7+", "cd4-cd8+/CD45RA+"
				)
lapply(nodesToHide, function(thisNode) gs_pop_set_visibility(gs, thisNode, FALSE))
```
### 4.7. Rename nodes
And rename the populations:
```{r rename, results = "hide"}
gs_pop_set_name(gs, "cd4+cd8-", "cd4")
gs_pop_set_name(gs, "cd4-cd8+", "cd8")
```

```{r plot_afterHiding}
plot(gs[[1]])
```


### 4.8. Visualize the gates

```{r plotGate_autoGate, fig.width = 9}
autoplot(gs[[1]])
```

### 4.9. Apply a gating method without csv template
Sometimes it will be helpful (especially when working with data that is already gated) to be able to interact with the `GatingSet` directly without the need to write the complete csv gating template. We can apply each automated gating method using the same fields as in the `gatingTemplate`, but provided as arguments to the `gs_add_gating_method` function. The populations added by each of these calls to `gs_add_gating_method` can be removed sequentially by `gs_remove_gating_method`, which will remove *all* populations added by the prior call to `gs_add_gating_method`. These two functions allow for interactive stagewise prototyping of a `gatingTemplate`.

For example, suppose we wanted to add a `CD38-/HLA-` sub-population to the `cd4+cd8-` population. We could do this as follows:


```{r gt_add_gating_method}
gs_add_gating_method(gs, alias = "non-activated cd4",
                         pop = "--",
                         parent = "cd4",
                         dims = "CD38,HLA",
                         gating_method = "mindensity")
plot(gs[[1]])
```

The addition of this population can then easily be undone by a call to `gs_remove_gating_method`:

```{r gs_remove_gating_method}
gs_remove_gating_method(gs)
plot(gs[[1]])
```


5. Conclusion
------------
The `openCyto` package allows users to specify their gating schemes and gate the data 
in a data-driven fashion. It frees the scientists from the labor-intensitive manual gating routines 
and increases the speed as well as the reproducibilty and objectivity of the data analysis work.    







