---
title: "Cytofast - vignette for the Spitzer study"
shorttitle: "Cytofast"
author: 
 - name: "K.A. Stam"
   affiliation: "Department of Parasitology, Leiden University Medical Center, Leiden, The Netherlands"
 - name: "G. Beyrend"
   affiliation: "Department of Immunohematology and Blood Transfusion, Leiden University Medical Center, Leiden The Netherlands"
package: cytofast
date: '`r Sys.Date()`'
output:
  BiocStyle::html_document
bibliography: REFERENCES.bib
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Spitzer} 
  %\VignetteEncoding{UTF-8}
---

```{r style, echo=FALSE, results = "asis"}
BiocStyle::markdown()
```


# Introduction
 
*Cytofast* is a R-package aimed at quickly visualizing and analyzing the output of [Cytosplore](https://www.cytosplore.org/); software used to characterizes and cluster immune subsets from flow or mass cytometry (CyTOF) data. Cytosplore is mainly used for its fine implementation of t-sne (and HSNE) algorithm to identify and cluster distinct populations of cells. However, further downstream analysis is unsupported, such as relating these population of cells to a phenotype and/or ascribing any experimental conditions. Therefore, we introduce the *cytofast* package, which provides a quick visualization and quantification of clusters (or combination of clusters) playing a key role in your study. 

In this vignette will guide the users through the basic workflow designed around the package. As an example we will use the mass cytometry dataset generated by Spitzer [-@Spitzer2017]. Here we work with a small example to explain the functionality of the package, for the complete workfow and a more thorough analysis we refer to our peer reviewed paper [@Beyrend2018]. 

The workflow is designed around several steps, which are addressed per section. In each section we showcase some functions, which may or may not be of your interest. A summary of the most used functions is as follows: 

1. `readCytosploreFCS`: reading *.fcs* files from Cytosplore into R
2. `cytoHeatmaps`: generating two aligned heatmaps 
+ the upper panel represents the phenotype (based on the markers) of each cluster
+ the lower panel represents the relative abundance of the clusters in respect to the samples 
3. `msiPlot`: generating histograms of the signal intensity for each marker, per group or cluster
4. `cytoBoxPlots`: creates boxplots (percentage of total given cells) and representing each cluster as a percentage of the total analyzed immune cells  
5. `cytottest` and `cytogt`, functions to quantify clusters in respect to experimental conditions. The latter function
is still under development.

# Installation

To install *cytofast* run the following code in the console:
```{r, eval=FALSE}
if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("cytofast")
```


# Creating clusters

Before *cytofast* can be applied, the data should be clustered. We recommend using Cytosplore for this task, however it is certainly possible to use any other clustering algorithm. As an example, we will also show how data clustered by [FlowSOM](https://bioconductor.org/packages/release/bioc/html/FlowSOM.html) can be parsed for usage with cytofast.

## Creating clusters in *Cytosplore*

After downloading *Cytosplore*, hosted at www.cytosplore.org, upload the FCS files by clicking on *File* then *Open FCS File(s)*. Be sure to add unique sample tag as channel when prompted and select a cofactor for hyperbolic arcsinh transformation (default is 5).

Click on *Run H-SNE* and wait for the map being generated. This step may require some time depending on the number of cells analyzed.

## Exporting the data
 
Once the tSNE map is generated, we can save the clusters defined by *Cytosplore* by right-clicking on the tSNE map and save the clusters as **"Save Clusters ..."**". We can choose the directory of the output files as prompted by *Cytosplore*; make sure to note this location, because we will use this directory to load the .fcs files into R.

We advise using a simple name (characters only) when renaming the output files, this makes identification and further handling easier. After saving, we will have a folder with the created .fcs files, with each file corresponding to your identified clusters in *Cytosplore*. The next step will be loading the files into R with the help of `cytofast`. 

# Reading the data 

Let's start with loading the package.

```{r}
library(cytofast)
```

To follow this vignette we added some pre-clustered (by Cytosplore) .fcs files to the package. These are ten downsampled clusters from the Spitzer study and should be treated as a toy-example to discover the package. A complete analysis of the Spitzer data is reported and published in *Beyrend et al., CSBJ, 2018*.

To import the *.fcs* files created by Cytosplore, we use the function `readCytosploreFCS`. Simply give the directory where you saved your .fcs files and the function will store the data in a `cfList` (cytofast list). When not using Cytosplore, check the function `cfList` to create such a list yourself. It is an S4 class object with several slots, a data frame `expr` containing the expression of the markers for all measured cells, `counts` a data frame with cell counts per cluster per sample and `samples` a data frame with all the meta information.

```{r echo=FALSE}
cfList(samples = data.frame(sampleID = as.factor(1:10)),
       expr = data.frame(sampleID = as.factor(1:10), clusterID=letters[1:10]), 
       counts = data.frame())
```

Note that the `expr` slot needs both a column named "sampleID" and "clusterID", which are used to identify the origin of each individual cell. When using `readCytosploreFCS` this is all done automatically. 

From here you can choose to use your own data, or it's possible to follow along with the added data from Spitzer. So, let's load the data into R:

```{r, eval=TRUE}
dirFCS <- system.file("extdata", package="cytofast")
cfData <- readCytosploreFCS(dir = dirFCS, colNames = "description")
```

In total there were 10 .fcs files (also meaning 10 clusters) and are now stored into a `cfList`. We can acces the list and have a quick look at the marker expression of the cells.

```{r, eval=T}
cfData@expr[1:5, 1:5]
```

As shown, both clusterID and sampleID are added to the expression data. We also see that there are some 'markers' that we don't want in our further analysis. Therefore, in the next section we will clean-up the data and add some meta information to the `cfList`. 

## Cleaning up the data

Some channels are not needed (e.g. DNA1, DNA2, barcoding, etc...). We can remove those and only keep the markers that we are actually intrested in. 

```{r, eval=T}
cfData@expr <- cfData@expr[,-c(3:10, 13:16, 55:59, 61:63)]
```

Next we will add some meta data to the `cfList`. Make sure that the sampleID corresponds with the sampleID of the meta data! When using Cytosplore, the sampleID is based on the sample tag (CSPLR_ST). 

```{r, eval=T}
data(spitzer) 
meta <- spitzer[match(row.names(cfData@samples), spitzer[,"CSPLR_ST"]),] # match sampleID
cfData@samples <- cbind(cfData@samples, meta[,2:3])
```

Lastly, we will make some small changes to the labeling of the data. This helps with the visualization of the heatmaps and discovering relations of clusters.

```{r, eval=T}
levels(cfData@expr[,"clusterID"]) <- gsub("[^0-9]", "", levels(cfData@expr[,"clusterID"]))  
```

# 3. Visualization

## Heatmaps

Before we can create the heatmaps, we first address the last slot `counts`, which is still empty. To fill this slot we can call the function `cellCounts`. The default function will simply add the raw cell counts per cluster per sample into the cfList. 

```{r, eval=F}
cfData <- cellCounts(cfData)
head(cfData@counts)
```

In many cases the abundance of the cells could be biased in respect to the sample size of the donors, therefore it is better to look at the frequency of the clusters in respect to the total amount of cells per sample. We will also standardize the data (mean zero, unit variance) for a better comparison between clusters.

```{r, eval=T}
cfData <- cellCounts(cfData, frequency = TRUE, scale = TRUE)
head(cfData@counts)
```

Now everything is ready to visualize the data. One of the main functions of this package is `cytoHeatmaps`, which is used to visualize both the phenotype of the created clusters and their heterogeneity in respect to the samples.

```{r fig.width=10, fig.height=12, fig.cap="\\label{fig:fig2}Cytofast heatmap"}
cytoHeatmaps(cfData, group="group", legend=TRUE)
```

The function will plot two heatmaps:
+ on the upper panel the phenotype (based on the markers) of each cluster
+ on the lower panel the relative abundance of the clusters in respect to the samples  
A hierarchical clustering was performed on subset frequencies using the Spearmans rank correlation, resulting in a dendogram displayed on the side of the sample abundancy heatmap. This will group samples sharing similarities.  
Groups can pre-defined as an input to 'cytoHeatmaps'. 

\clearpage

## boxplots  

We can also represents the data obtained in a quantitative way by calling the function `cytoBoxplots`. The output of this function represents the proportion of each sample for every cluster.  

```{r, eval=T}
cytoBoxplots(cfData, group = "group")
```

## density plots (median signal intensity)

According to the *Spitzer et al., Cell, 2017* study, an effective therapy lead to an upregulation of MHC-II. This actually can be seen by drawing the median intensity plots for this marker, as shown below together with CD45 and CD4. 

```{r, eval=T}
msiPlot(cfData, markers = c("MHC.II", "CD45", "CD4"), byGroup='group')
```

On this MSI plot, it indeed appears that the MHC-II signal is higher in the Effective therapy day 3 and day 8 compared to the Ineffective therapy day 3 and day 8.  

# 4. Statistical testing

## t.test

```{r, eval=F}
cfData@samples$effect <- gsub("_D\\d", "", spitzer$group) # difference between effective and ineffective
cfData <- cytottest(cfData, group="effect", adjustMethod = "bonferroni")
```


## global test

# 5. Extra

## clustering with flowSOM

It is also possible to use other clustering algorithms 

```{r, eval=F}
library(FlowSOM)
fSOM <- FlowSOM(input = dirFCS, 
                transform = FALSE,
                scale = FALSE,
                colsToUse = c(9:11, 15:52),
                nClus = 10, # We simply choose 10 clusters here
                seed = 123)
```

From the flowSOM object we can extract the mapping and relevel each individual cell to its meta cluster. 

```{r, eval=F}
clusterID_FS <- as.factor(fSOM$FlowSOM$map$mapping[,1])
levels(clusterID_FS) <- fSOM$metaclustering
```

From here it is of course possible to construct a new cfList and go through the whole workflow, but as simplification we just swap `clusterID` with `clusterID_FS` and only show the heatmap.

```{r fig.width=10, fig.height=12, fig.cap="\\label{fig:fig3}heatmap based on flowSOM", eval=F}
cfData@expr$clusterID <- clusterID_FS
cfData <- cellCounts(cfData) # Update cellCounts with new clusters
cytoHeatmaps(cfData, group='group', legend=TRUE)
```

\clearpage

## References
