In JonathanIrish/MEMv3: Marker Enrichment Modeling (MEM)

# Load all libraries
# If you get an error message, you will need to try re-installing packages
library(FlowSOM)
library(flowCore)
library(Biobase)
library(gplots)
library(ggplot2)
library(hexbin)
library(MEM)
library(tidyverse)
library(Rtsne)
library(uwot)
library(viridis)
library(ggExtra)
library(RColorBrewer)

# read FCS files into R
setwd(paste(getwd(), "/datafiles/MCBead", sep = ""))
file <-  dir(pattern = "*.fcs")

# Run MEM
MEM.values.orig = MEM(
  file,
  transform = TRUE,
  cofactor = 5,
  choose.markers = FALSE,
  markers = "2,13,14,21,23:55",
  choose.ref = FALSE,
  zero.ref = FALSE,
  rename.markers = FALSE,
  file.is.clust = FALSE,
  add.fileID = FALSE,
  IQR.thresh = NULL
)

# build MEM heatmap and output enrichment scores
build.heatmaps(
  MEM.values.orig,
  cluster.MEM = "none",
  cluster.medians = "none",
  display.thresh = 1,
  newWindow.heatmaps = FALSE,
  output.files = FALSE,
  labels = FALSE,
  only.MEMheatmap = TRUE
)

# prepare data for use in UMAP
FCS = read.FCS(file)
data.df = as.data.frame(FCS@exprs)
chosen.markers = data.df[, c(2, 13, 14, 21, 23:55)]
transformed.chosen.markers <- chosen.markers %>%
  mutate_all(function(x)
    asinh(x / 5))
overall_seed = 43

# Run UMAP on chosen markers
set.seed(overall_seed)
myumap <- umap(transformed.chosen.markers, 
               ret_model = TRUE, 
               n_threads = 1,
               verbose = TRUE)
umap.data = as.data.frame(myumap$embedding)

range <- apply(apply(umap.data, 2, range), 2, diff)
graphical.ratio <- (range[1]/range[2])

# UMAP flat dot plot and density dot plot
UMAP.plot <- data.frame(x = umap.data[,1], y = umap.data[,2])

ggplot(UMAP.plot) + coord_fixed(ratio=graphical.ratio) + geom_point(aes(x=x, y=y), cex = 1) + labs( x = "UMAP 1", y = "UMAP 2", title = "Mass Cytometry Beads") + theme_bw()

ggplot(UMAP.plot, aes(x=x, y=y)) + coord_fixed(ratio = graphical.ratio)  + geom_bin2d(bins = 128) + 
  scale_fill_viridis_c(option = "A", trans = "sqrt") + scale_x_continuous(expand = c(0.1,0)) + 
  scale_y_continuous(expand = c(0.1,0)) + labs(x = "UMAP 1", y = "UMAP 2",title = "Mass Cytometry Beads") + theme_bw()

# Run FlowSOM on the UMAP axes
umap.mat <- as.matrix(umap.data)

# create flowFrame
UMAP.metadata <-
  data.frame(name = dimnames(umap.mat)[[2]],
             desc = paste('UMAP', dimnames(umap.mat)[[2]]))
UMAP.metadata$range <- apply(apply(umap.mat, 2, range), 2, diff)
UMAP.metadata$minRange <- apply(umap.mat, 2, min)
UMAP.metadata$maxRange <- apply(umap.mat, 2, max)
umap.flowframe <- new("flowFrame",
                      exprs = umap.mat,
                      parameters = AnnotatedDataFrame(UMAP.metadata))

# implement the FlowSOM on the data
fsom <-
  FlowSOM(
    umap.flowframe,
    compensate = FALSE,
    transform = FALSE,
    toTransform = c(1:2),
    scale = TRUE,
    colsToUse = c(1:2),
    nClus = 38,
    seed = overall_seed
  )
FlowSOM.clusters <-
  as.matrix(fsom[[2]][fsom[[1]]$map$mapping[, 1]])

qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, 
                           rownames(qual_col_pals)))
col_vector = col_vector[-c(4,17,19,27,29:45)]
values = sample(col_vector)

# plot FlowSOM clusters on UMAP axes
ggplot(UMAP.plot) + coord_fixed(ratio=graphical.ratio) + 
  geom_point(aes(x=x, y=y, color=FlowSOM.clusters),cex = 1.5) + 
  guides(colour = guide_legend(override.aes = list(size=5), nrow = 13)) +
  labs(x = "UMAP 1", y = "UMAP 2",title = "FlowSOM Clustering on UMAP Axes", 
       color = "FlowSOM Cluster") + theme_bw() + 
  scale_color_manual(values = values)

# Run MEM on the FlowSOM clusters from UMAP
cluster = as.numeric(as.vector((FlowSOM.clusters)))
MEM.data = cbind(transformed.chosen.markers, cluster)

MEM.values.uf = MEM(
  MEM.data,
  transform = FALSE,
  cofactor = 5,
  choose.markers = FALSE,
  markers = "all",
  choose.ref = FALSE,
  zero.ref = FALSE,
  rename.markers = FALSE,
  file.is.clust = FALSE,
  add.fileID = FALSE,
  IQR.thresh = NULL
)

# build MEM heatmap and output enrichment scores
build.heatmaps(
  MEM.values.uf,
  cluster.MEM = "none",
  cluster.medians = "none",
  display.thresh = 1,
  newWindow.heatmaps = FALSE,
  output.files = FALSE,
  labels = FALSE,
  only.MEMheatmap = TRUE
)

# RMSD to compare labels from all populations
orig.MEM.scores = as.data.frame(MEM.values.orig[[5]])
rownames(orig.MEM.scores) = paste0(rownames(orig.MEM.scores), " Manual")
uf.MEM.scores = as.data.frame(MEM.values.uf[[5]])
rownames(uf.MEM.scores) = paste0(rownames(uf.MEM.scores), ' UMAP')
all.MEM.values = as.matrix(rbind(uf.MEM.scores, orig.MEM.scores))

RMSD_vals <-
  MEM_RMSD(
    all.MEM.values,
    format = NULL,
    newWindow.heatmaps = FALSE,
    output.matrix = FALSE
  )