FlowCore.Rmd

---
title: "FlowCore analysis"
output:
  html_document:
    code_folding: hide
---


```{r, include=FALSE}
knitr::opts_chunk$set(cache=TRUE, warning = FALSE)
library(flowCore)
library(CATALYST)
library(umap)
library(R.utils)
library(yaml)
library(tidyverse)
library(SingleCellExperiment)
library(scater)
library(diffcyt)
library(BiocParallel)
config <- read_yaml("config.yml")
register(MulticoreParam())
bpparam()
```

# Input data

Input the data from design file into a SingleCellExperiment class object using the wrapper function prepData from CATALYST. 

```{r input}
# design_<dir>_<sample>.csv
design_files <- list.files(pattern = "design_")

.unlink <- function(x, recursive = FALSE, force = FALSE) {
  if (unlink(x, recursive, force) == 0)
    return(invisible(TRUE))
  stop(sprintf("Failed to remove [%s]", x))
}

unlink("input/", recursive = TRUE)

for (i in design_files){
  x <- strsplit(i, "_")
  condition <- x[[1]][2]
  value <- x[[1]][3]
  value <- gsub(".csv","", value)
  dir_files <- paste0("input/", value)
  dir.create(file.path("input"), showWarnings = FALSE)
  dir.create(file.path(dir_files), showWarnings = FALSE)
  
  md <- read.csv(i)

  file_names <- as.vector(md$file_name)
  linked_files <- gsub("./","", file_names)
  linked_files <- paste0(getwd(), "/" ,linked_files)
  
  
  output_dir <- paste0(getwd(), "/", dir_files, "/")
  
  
  for (y in linked_files){
    file.symlink(y, output_dir)
  } 
  
  
  fcs_files <- list.files(path = output_dir, pattern = ".fcs$", full.names = T)
  md$file_name <- fcs_files
  fs <- read.flowSet(fcs_files, transformation = FALSE, truncate_max_range = FALSE)
  
  
  panel <- read.csv("metals_meta_expressed.csv", stringsAsFactors = T)

  md$condition <- as.factor(md$condition)
  md$sample_id <- as.factor(md$sample_id)
  
  sce <- prepData(fs, panel, md)
  
  # Define a downsampling ceiling
  sampling.ceiling <- as.integer(config$downsample)
  # Being reproducible is a plus
  set.seed(0)

  ## BUILD A DOWNSAMPLED FLOWSET

  fs.ds <- fsApply(fs, function(ff) {
    idx <- sample.int(nrow(ff), min(sampling.ceiling, nrow(ff)))
    ff[idx,]  # alt. ff[order(idx),]
  })
  fs.ds
  
  sce.downsample <- prepData(fs.ds, panel, md)
  
  assign(paste("sce", i, sep = "."), sce)
  assign(paste("sce.downsample", i, sep = "."), sce.downsample)
}
```

## Plot counts {.tabset .tabset-fade}

```{r expression, results='asis', echo = FALSE}
for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.SAMPLE_FILE"))
  cat("### ",i,"\n")
  p <- plotCounts(sce, color_by = "condition")

  print(p)
  cat('\n\n')
}
```


## Heatmap of (scaled) median marker expressions {.tabset .tabset-fade}

plotExprHeatmap will show a heatmap on median marker intensities with hierarchically clustered rows (samples) and columns (markers). This plot should give an idea of which markers will drive sample clustering, and how similiar samples are in their expression profile. We specify bin_anno = TRUE to display expression values inside each bin, and row_anno = TRUE to include row annotations for each factor in metadata(daf).

```{r marker expression, results='asis', echo = FALSE}

for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.SAMPLE_FILE"))
  cat("### ",i,"\n")
  p <- CATALYST::plotExprHeatmap(sce, bin_anno = TRUE, row_anno = TRUE)
  print(p)
  cat('\n\n')
}
```


## Cluster cells {.tabset .tabset-fade}

Cluster cells based on the use_for_clustering in the panel data table.

```{r counts, results='asis', echo = FALSE}
# specify markers to use for clustering

for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.SAMPLE_FILE"))
  lineage_markers <- panel %>% 
    filter(use_for_clustering == TRUE) %>% 
    select(antigen) 

  sce <- CATALYST::cluster(sce, features = lineage_markers$antigen, 
      xdim = 10, ydim = 10, maxK = 20, verbose = FALSE, seed = 1)  

  p <- CATALYST::plotMedExprs(sce, k = "meta8", facet = "cluster_id", 
    group_by = "condition", shape_by = "patient_id")
  p$facet$params$ncol <- 4
  
  cat("### ",i,"\n")
  print(p)
  cat('\n\n')

  assign(paste("sce", i, sep = "."), sce)
}
```

## Delta area plot {.tabset .tabset-fade}

The delta area represents the amount of extra cluster stability gained when clustering into k groups as compared to k-1 groups. It can be expected that high stability of clusters can be reached when clustering into the number of groups that best fits the data. The “natural” number of clusters present in the data should thus corresponds to the value of k where there is no longer a considerable increase in stability (pleateau onset). For more details, the user can refer to the original description of the consensus clustering method (Monti et al. 2003).

```{r, delta_plot, results='asis', echo = FALSE}

for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.SAMPLE_FILE"))
  
  cat("### ",i,"\n")
  print(metadata(sce)$delta_area)
  cat('\n\n')
}
```

## Median marker-expressions by cluster {.tabset .tabset-fade}

We can facet the above plot by antigen in order to compare marker expressions calculated over all cells across conditions

```{r, results='asis', echo = FALSE}
for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.SAMPLE_FILE"))
  
  cat("### ",i,"\n")
  p <- CATALYST::plotMedExprs(sce, facet = "antigen", 
    group_by = "condition", shape_by = "patient_id")
p$facet$params$ncol <- 4
print(p)
  cat('\n\n')
}
```

## Marker-densities by cluster {.tabset .tabset-fade}

Distributions of marker intensities (arcsinh-transformed) across cell populations of interest can be plotted with plotClusterExprs. We specify features = "type" (equivalent to type_markers(sce)), to include type-markers only. Here, blue densities are calculated over all cells and serve as a reference.

```{r, results='asis', echo = FALSE}
for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.SAMPLE_FILE"))
  
  cat("### ",i,"\n")
  p <- plotClusterExprs(sce, k = "meta2", features = "type")

  print(p)
  cat('\n\n')
}
```

<!--Remember to add overriding clustering names here -->

# Heatmap {.tabset .tabset-fade}

Clusterings and metaclusters maybe be viewing with the plotClusterHeatmap. In its 1st panel, the function will display median (arcsinh-transformed and optionally scaled) cell-type marker expressions (across all samples). Depending on argument hm2, the 2nd panel will vary as follows:

"abundances": cluster frequencies by sample;
"state_markers": median cell state marker expressions across clusters (analogous to the left-hand side heatmap);
a character string/vector corresponding to one/multiple marker(s): median marker expressions by sample.
Argument scale (default TRUE) specifies whether scaled values should be plotted. These correspond to arcsinh-transformed expression values scaled between 0 and 1 using low (1%) and high (99%) percentiles as boundaries. Note that, in any case, hierarchical clustering is performed on the unscaled data.
While it is only used here for visualization, this additional transformation of the arcsinh-transformed data can sometimes give a better representation of relative differences in marker expression between cell populations.

```{r, results='asis', echo = FALSE}
for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.SAMPLE_FILE"))
  
  cat("### ",i,"\n")
  p <- CATALYST::plotClusterHeatmap(sce, hm2 = "abundances", draw_freqs = TRUE, cluster_anno = FALSE)
  print(p)
  cat('\n\n')
}
```


## Relative population abundances {.tabset .tabset-fade}

Relative population abundances for any clustering of interest can be plotted with plotAbundances. Argument by will specify whether to plot proportions for each sample or cluster.
If by = "sample_id", the function displays each sample’s cell type composition, and the size of a given stripe reflects the proportion of the corresponding cell type the given sample. Argument group then specifies the facetting. If by = "cluster_id", argument group then specifies the grouping and color coding.

```{r, results='asis', echo = FALSE}
for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.SAMPLE_FILE"))
  
  cat("### ",i,"\n")
  p <- CATALYST::plotAbundances(sce, k = "meta12", by = "sample_id", group_by = "condition")
  print(p)
  cat('\n\n')
}
```


```{r, results='asis', echo = FALSE}
for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.SAMPLE_FILE"))
  
  cat("### ",i,"\n")
  p <- plotAbundances(sce, k = "meta2", by = "cluster_id", group_by = "condition", shape_by = "patient_id")
  print(p)
  cat('\n\n')
}
```



# Dimensionality reduction 

## Multi-dimensional reduction 

```{r}
for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.SAMPLE_FILE"))
  set.seed(1601)
  
  cat("### ",i,"\n")
  p <- CATALYST::plotMDS(sce, color_by = "condition")
  print(p)
  cat('\n\n')
}
```


## UMAP {.tabset .tabset-fade}

```{r, results='asis', echo = FALSE}

# Need to downsample all of the datasets because CyTOF data is too large
# Will be performed using DropletUtils downsampleMatrix

for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.downsample.SAMPLE_FILE"))
  set.seed(1601)

  sce <- scater::runUMAP(sce, exprs_values="exprs", BPPARAM=bpparam())
  cat("### ",i,"\n")
  p <- plotDR(sce, "UMAP") + facet_wrap("condition")
  print(p)
  cat('\n\n')
}
```


## TSNE {.tabset .tabset-fade}

```{r, results='asis', echo = FALSE, include=FALSE}
for (i in design_files){
  sce <- get(gsub("SAMPLE_FILE",i , "sce.downsample.SAMPLE_FILE"))
  set.seed(1601)
  sce <- scater::runTSNE(sce, exprs_values = "exprs", BPPARAM=bpparam())
  cat("### ",i,"\n")
  p <- plotDR(sce, "TSNE") + facet_wrap("condition")
  print(p)
  cat('\n\n')
}
```

# Differential analysis {.tabset .tabset-fade}

Modify (and uncomment) the following document below to include bespoke analysis code for your  project.

```{r}
#design <- createDesignMatrix(ei(sce), cols_design = "condition")
#contrast <- createContrast(c(0, 1))
```