Part1.Rmd

---
title: "Automated cell type annotation - Benchmarking different Machine Learning approaches"
author: "Dimitrios Kleftogiannis"
date: "2023-2-13"
output: html_document
---
  
### Utility
  
This code is part of the study titled "Automated cell type annotation and exploration of single-cell signalling dynamics using mass cytometry".

The utility of this code is to benchmark different computational approaches for automated cell type annotation.  

Please note that term "automated cell type annotation" describe machine learning-based approaches that go beyond the state-of-the-art unsupervised clustering, or gating schemes  The problem is more precisely formulated as classification task where labelled data are used as "landmarks/anchors" to predict the cell type of unknown cells. 

In the most extreme scenario the task is to use one annotated dataset as reference and predict the cell types in a completely different dataset that come from different experiments. 
  
### Contact
  
Comments and bug reports are welcome, please email: Dimitrios Kleftogiannis (dimitrios.kleftogiannis@uib.no)

We are also interested to know about how you have used our framework, including any improvements that you have implemented.

You are free to modify, extend or distribute our source codes, as long as our copyright notice remains unchanged and included in its entirety. 

### License

This code is licensed under the MIT License.

Copyright 2023, University of Bergen (UiB) and NeuroSysMed, Norway


### Load R libraries.
```{r load packages, echo=TRUE, eval=TRUE, error=TRUE, warning=FALSE,cache=TRUE}

library(HDCytoData)
library(ggplot2)
library(caret)
library(MASS)
library(dplyr)
library(ggridges)

```

### Initialise the working space and datasets included in the comparison analysis
We use several datasets for evaluating the effectiveness of different supervised and semi-supervised cell type annotation methods.

To build a reference of human hematopoiesis we utilise the REFERENCE DATASET from 
Velten, Lars; Triana, Sergio; Vonficht, Dominik; Jopp-Saile, Lea; Paulsen, Malte; Haas, Simon (2021). 

The dataset is named **Healthy.rds** and contains expression of 97 surface markers and 462 mRNAs in 49057 cells from healthy young and healthy old bone marrow. 

The dataset can be downloaded from https://doi.org/10.6084/m9.figshare.13397651.v4

We also utilize different benchmark data sets namely the **AML_benchmark** and **BMMC_benchmark** data sets included in the study:
Abdelaal et al. titled Predicting Cell Populations in Single Cell Mass Cytometry Data  https://doi.org/10.1002/cyto.a.23738

All data sets can be downloaded from Flow Repository (http://flowrepository.org/id/FR-FCM-ZYTT)

Lastly we use the **Samusik_all_SE** data set from the R package HDCytoData, available here https://rdrr.io/github/lmweber/HDCytoData/

```{r initialise workspace and load dataset,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}
setwd("/Users/kleftogi/Desktop/CyTOF_paper_corrections")
source('loadReferenceData.R')
source('extractRefUMAPCoordinates.R')
source('filterAndTabulate.R')

#first we read the Healthy.rds data set. It is generated using AB-seq technology, measured simultaneously surface markers and mRNAs.
#however, here we will use antibodies and we will neglect the information from mRNAs since it is not compatible with cyTOF technology

#we use a utility function called loadReferenceData that was sourced in the start of this chunk
filename <- 'Healthy.rds'
o <- loadReferenceData(filename)
healthy.AB.data <- o[[1]]
healthyCellTypes <- o[[2]]
rm(o)

#we can extract the original UMAP coordinates. We use function extractRefUMAPCoordinates that has been loaded too
umap <- extractRefUMAPCoordinates(filename)

#since the dataset contains many "rare" cell populations with not so many cells we set a cutoff and we discard those cell populations with less than 300 cells. For this purpose we use the function filterAndTabulate that has been loaded already. Please note that the function changes the names of the markers, by removing suffix -AB. In the next step we will see why we need to use the antibody names as column names in the dataset. 
myCellCutoff <- 300
o <- filterAndTabulate(healthy.AB.data,healthyCellTypes,umap,myCellCutoff)
healthy.AB.data.tab <- o[[1]]
healthy.AB.data.filtered <- o[[2]]
umap.filtered <- o[[3]]
rm(o)

#visualisation of the reference data set, cells are coloured using CD45 and CD34 markers expression as an example

#initialise colors to be used in the umap colorbar
color_grad_flow2 <- c("#331820", "#4d4f55", "#55626b", "#5a767e", "#628b8a", "#709b91", "#82aa96", "#98b89a", "#b0c6a2", "#c9d3ab", "#e4e0b6", "#feedc3")

tmp <- cbind(umap.filtered,healthy.AB.data.filtered$CD45)
colnames(tmp)[3] <- 'CD45'
umap_healthy_ref_1.1 <- ggplot(tmp, aes(x = UMAP1, y = UMAP2,color=CD45)) +
    geom_point(size = .1) +
    coord_fixed(ratio = 1) +
    ggtitle("Roadmap of human hematopoiesis")+
    scale_colour_gradientn(colours = color_grad_flow2, limits = c(0,1.5),name='CD45')+
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank(),
          legend.position = "bottom",
          legend.background = element_rect())

tmp <- cbind(umap.filtered,healthy.AB.data.filtered$CD34)
colnames(tmp)[3] <- 'CD34'
umap_healthy_ref_1.2 <- ggplot(tmp, aes(x = UMAP1, y = UMAP2,color=CD34)) +
    geom_point(size = .1) +
    coord_fixed(ratio = 1) +
    ggtitle("Roadmap of human hematopoiesis")+
    scale_colour_gradientn(colours = color_grad_flow2, limits = c(0,1.5),name='CD34')+
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank(),
          legend.position = "bottom",
          legend.background = element_rect())

# plot umap annotated cell types
db31 <- c('#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#911eb4', '#46f0f0', '#f032e6', 
          '#bcf60c', '#fabebe', '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000', 'slategray3',
          'khaki3','bisque3','coral1','mediumaquamarine','royalblue1','gold4','orange3',
          '#aaffc3', '#808000', '#ffd8b1', '#000075', '#808080','gray88' ,'#ffffff', '#000000')

tmp <- cbind(umap.filtered,healthy.AB.data.filtered$cellType)
colnames(tmp)[3] <- "cellType"
umap_healthy_ref_2 <- ggplot(tmp, aes(x = UMAP1, y = UMAP2,color=cellType)) +
    geom_point(size = .1) +
    coord_fixed(ratio = 1) +
    ggtitle("Roadmap of human hematopoiesis")+
    scale_color_manual(values=db31)+
    theme(axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank(),
          legend.position = "bottom",
          legend.background = element_rect())+
  guides(color = guide_legend(override.aes = list(size=1.4),ncol=3,title=""))

#produce some statistics about the cell populations and their relative frequency

 cellFreq <- healthy.AB.data.filtered %>% 
  group_by(cellType) %>% 
  summarise(samples = n())%>% 
  mutate(freq = samples/sum(samples))

 healthy_ref_cellFreq <- ggplot(cellFreq) + aes(x=reorder(cellType,-freq),y=freq,fill=cellType)+
  geom_bar(stat='identity',width = 0.58,alpha=1,color='black',size=0.2)+
  ylab('Relative abundance of cell types')+
  theme_bw()+
  theme(axis.text.y = element_text( size = 10 ),
        axis.text.x = element_text(angle = 0, vjust = 0.5, hjust = 0.5, size = 10),
        axis.title.x = element_text( size = 12,face = 'bold' ),
        axis.title.y = element_blank(),
        strip.text = element_text(size = 12,face='bold',lineheight=1),
        legend.position = "none",aspect.ratio = 1)+
  scale_fill_manual(values = db31)+
  guides(fill = guide_legend(override.aes = list(size=3),nrow=1,title=""))+
  coord_flip()

 print(umap_healthy_ref_1.1)
 print(umap_healthy_ref_1.2)
 print(umap_healthy_ref_2)
 print(healthy_ref_cellFreq)
 
 
```

### Self-consistency comparison using training and testing sets from the Healthy.rds dataset with random split
First we will perform a self consistency test focusing on **Healthy.rds** dataset. 

This means that we will generate random disjoint training and testing sets by sub-sampling randomly cells from healthy.AB.data.filtered data frame  without replacement. 

Since the original data set contains 97 antibodies that are not possible to incorporate in conventional CyTOF panels, we will select a subset of 30 surface markers that will serve as cell identity markers for our experimentation.

Will train models using different approaches namely: 
1) the Scaffold approach  originally described by Spitzer et al  DOI: 10.1126/science.1259425
2) two different versions of K-Nearest-Neighbors (KNN); 
3) Linear Discriminant Analysis (LDA) and
4) CyAnno which deploys a sophisticated ensemble of classifiers (XGboost, LDA and SVM) published by Kaushik et al. DOI: 10.1093/bioinformatics/btab409

```{r self-consistency analysis using Healthy.rds,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}

#load functions required for the scaffold approach
source('tabulateCellTypes.R')
source('unsupervised.R')

#define antibodies/features for the self consistency test
surfaceMarkersSubset <- c("CD11b","CD123","CD14","CD16","CD20",
                          "CD235a-b","CD25","CD3","CD33","CD34",
                          "CD38","CD4","CD45","CD45RA","CD56",
                          "CD61","CD8","IgG","CD194","CD195",
                          "CD196","CD133","CD10","CD336","CD9",
                          "CD5","CD155","CD7","CD279","HLAC")

#select specific columns for the rest of the analysis
data <- healthy.AB.data.filtered[,c(surfaceMarkersSubset,"cellType")]
#suffle the rows and get the indexes
IDX <- sample(nrow(data))
#cells for training is 80% and for testing is the remaining 20%
Ncell <- floor(0.8*nrow(data))
selectTrainIdx <- IDX[1:Ncell]
trainingSet <- data[selectTrainIdx,]
selectTestIdx <- IDX[(Ncell+1):length(IDX)]
testingSet <- data[selectTestIdx,]

#check the dimensions of the data
dim(trainingSet)
dim(testingSet)
nrow(testingSet)+nrow(trainingSet)==nrow(data)

#APPROACH #1: Cosine similarity to the scaffold, we assume that the testing data are already clustered
#We use the maximum distance from the landmarks (i.e., our training set) to annotate cell types. 

testingSetLabels <- testingSet$cellType
testingSetValues <- as.matrix(testingSet[,1:(ncol(testingSet)-1)])
testingSet.tab <- tabulateCellTypes(testingSetValues,testingSetLabels)
m <- as.matrix(testingSet.tab[, surfaceMarkersSubset])
rownames(m) <- testingSet.tab$cellType

trainingSetLabels <- trainingSet$cellType
trainingSetValues <- as.matrix(trainingSet[,1:(ncol(trainingSet)-1)])
trainingSet.tab <- tabulateCellTypes(trainingSetValues,trainingSetLabels)
#perform the same computation for the training set
att <- as.matrix(trainingSet.tab[,surfaceMarkersSubset])
row.names(att) <- trainingSet.tab$cellType
#compute the distance from the testing set to the training set which is the landmark using cosine similarity
dd_controls_to_landmarks <- t(apply(m, 1, function(x, att) {cosine_similarity_from_matrix(x, att)}, att))
#careful here, columns of the distance matrix are the rows of the landmark matrix
colnames(dd_controls_to_landmarks) <- rownames(att)

#apply hard filtering to remove connections where distance lower than the threshold --> adopted from the original implementation
hard_thres <- 0.75
for(i in 1:nrow(dd_controls_to_landmarks)){
  for(j in 1:ncol(dd_controls_to_landmarks)){
    if(dd_controls_to_landmarks[i,j]<hard_thres){
      dd_controls_to_landmarks[i,j] <- 0
    }
  }
}
#we parse dd_controls_to_landmarks and for each row/celltype we find the column with the maximum score. We use the column name as label for the predicted cell type, since it represents the closest landmark
predictionsScaffold <- data.frame()
for(idx in 1:nrow(dd_controls_to_landmarks)){
  
  realLabel <- rownames(dd_controls_to_landmarks)[idx]
  #find the predicted label
  a <- which.max(dd_controls_to_landmarks[idx,])
  predictedLabel <- colnames(dd_controls_to_landmarks)[a]
  distanceToLandmark <- dd_controls_to_landmarks[idx,a]
  dt <- data.frame(realLabel=realLabel,
                   predictedLabel=predictedLabel,
                   Dist=distanceToLandmark)
  predictionsScaffold <- rbind(predictionsScaffold,dt)
}
#generate a confusion matrix for this multi-class problem
#"TP of C1 is all  instances classified as C1 that are really C1.
#"TN of C1 is all instances not classified as C1 that are not really C1.
#"FP of C1 is all instances classified as C1 that are not really C1.
#"FN of C1 is all instances not classified as C1 that are really C1.
performanceScaffold <- data.frame()
for(idx in 1:nrow(predictionsScaffold)){
  
  currentType <- predictionsScaffold[idx,'realLabel']
  
  tp_row <- which(predictionsScaffold$realLabel==currentType & predictionsScaffold$predictedLabel==currentType)
  tn_row <- which(predictionsScaffold$realLabel!=currentType & predictionsScaffold$predictedLabel!=currentType)
  fp_row <- which(predictionsScaffold$realLabel!=currentType & predictionsScaffold$predictedLabel==currentType)
  fn_row <- which(predictionsScaffold$realLabel==currentType & predictionsScaffold$predictedLabel!=currentType)
  
  tp_row <- length(tp_row)
  tn_row <- length(tn_row)
  fp_row <- length(fp_row)
  fn_row <- length(fn_row)
  
  a <- data.frame(Sen=tp_row/(tp_row+fn_row),
                  Spe=tn_row/(tn_row+fp_row),
                  #PPV=tp_row/(tp_row+fp_row),
                  F1=2*tp_row/(2*tp_row+fp_row+fn_row),
                  Acc=(tp_row+tn_row)/(tp_row+tn_row+fp_row+fn_row))
  performanceScaffold <- rbind(performanceScaffold,a)
}
performanceScaffold$CellType <- predictionsScaffold$realLabel

#APPROACH #2: KNN with k=3
#We use the KNN classifier to predict the labels for the test set. 
#Caret R package implementation of KNN is used, and the package has been loaded from the beginning. 
myK <- 3
model_KNN3 <- train(trainingSetValues, trainingSetLabels,
                   method = 'knn',
                   tuneLength = myK)
predictedLabel <- predict(model_KNN3, newdata = testingSetValues)
predictionsKNN3 <- data.frame(realLabel=testingSetLabels,
                   predictedLabel=predictedLabel)
performanceKNN3 <- data.frame()
for(idx in unique(testingSetLabels)){
  
  currentType <- idx
  
  tp_row <- which(predictionsKNN3$realLabel==currentType & predictionsKNN3$predictedLabel==currentType)
  tn_row <- which(predictionsKNN3$realLabel!=currentType & predictionsKNN3$predictedLabel!=currentType)
  fp_row <- which(predictionsKNN3$realLabel!=currentType & predictionsKNN3$predictedLabel==currentType)
  fn_row <- which(predictionsKNN3$realLabel==currentType & predictionsKNN3$predictedLabel!=currentType)
  
  tp_row <- length(tp_row)
  tn_row <- length(tn_row)
  fp_row <- length(fp_row)
  fn_row <- length(fn_row)
  
  a <- data.frame(Sen=tp_row/(tp_row+fn_row),
                  Spe=tn_row/(tn_row+fp_row),
                  #PPV=tp_row/(tp_row+fp_row),
                  F1=2*tp_row/(2*tp_row+fp_row+fn_row),
                  Acc=(tp_row+tn_row)/(tp_row+tn_row+fp_row+fn_row))
  performanceKNN3 <- rbind(performanceKNN3,a)
}
performanceKNN3$CellType <- unique(testingSetLabels)

#APPROACH #3: KNN with k=10
#Same as approach #2, here we use the KNN classifier to predict the labels for the test set with different K
myK <- 10
model_KNN10 <- train(trainingSetValues, trainingSetLabels,
                   method = 'knn',
                   tuneLength = myK)
predictedLabel <- predict(model_KNN10, newdata = testingSetValues)
predictionsKNN10 <- data.frame(realLabel=testingSetLabels,
                   predictedLabel=predictedLabel)
performanceKNN10 <- data.frame()
for(idx in unique(testingSetLabels)){
  
   currentType <- idx
  
  tp_row <- which(predictionsKNN10$realLabel==currentType & predictionsKNN10$predictedLabel==currentType)
  tn_row <- which(predictionsKNN10$realLabel!=currentType & predictionsKNN10$predictedLabel!=currentType)
  fp_row <- which(predictionsKNN10$realLabel!=currentType & predictionsKNN10$predictedLabel==currentType)
  fn_row <- which(predictionsKNN10$realLabel==currentType & predictionsKNN10$predictedLabel!=currentType)
  
  tp_row <- length(tp_row)
  tn_row <- length(tn_row)
  fp_row <- length(fp_row)
  fn_row <- length(fn_row)
  
  a <- data.frame(Sen=tp_row/(tp_row+fn_row),
                  Spe=tn_row/(tn_row+fp_row),
                  #PPV=tp_row/(tp_row+fp_row),
                  F1=2*tp_row/(2*tp_row+fp_row+fn_row),
                  Acc=(tp_row+tn_row)/(tp_row+tn_row+fp_row+fn_row))
  performanceKNN10 <- rbind(performanceKNN10,a)
}
performanceKNN10$CellType <- unique(testingSetLabels)

#APPROACH #4: LDA 
tmp <- data.frame(trainingSetValues)
tmp$CellType <- trainingSetLabels
modelLDA <- lda(CellType~.,data=tmp)
predictedLabel <- predict(modelLDA, newdata = data.frame(testingSetValues))
#save the max posterior probability
prob_matrix<-matrix(0L, nrow = nrow(testingSetValues), ncol =1)
prob_matrix <- as.matrix(apply(predictedLabel$posterior,1,max))

predictionsLDA <- data.frame(realLabel=testingSetLabels,
                   predictedLabel=predictedLabel$class,
                   Posterior=prob_matrix[,1])
performanceLDA <- data.frame()
for(idx in unique(testingSetLabels)){
  
 currentType <- idx
  
  tp_row <- which(predictionsLDA$realLabel==currentType & predictionsLDA$predictedLabel==currentType)
  tn_row <- which(predictionsLDA$realLabel!=currentType & predictionsLDA$predictedLabel!=currentType)
  fp_row <- which(predictionsLDA$realLabel!=currentType & predictionsLDA$predictedLabel==currentType)
  fn_row <- which(predictionsLDA$realLabel==currentType & predictionsLDA$predictedLabel!=currentType)
  
  tp_row <- length(tp_row)
  tn_row <- length(tn_row)
  fp_row <- length(fp_row)
  fn_row <- length(fn_row)
  
  a <- data.frame(Sen=tp_row/(tp_row+fn_row),
                  Spe=tn_row/(tn_row+fp_row),
                  #PPV=tp_row/(tp_row+fp_row),
                  F1=2*tp_row/(2*tp_row+fp_row+fn_row),
                  Acc=(tp_row+tn_row)/(tp_row+tn_row+fp_row+fn_row))
  performanceLDA <- rbind(performanceLDA,a)
}
performanceLDA$CellType <- unique(testingSetLabels)

##APPROACH #5: CytoAnno
#We also try the CyAnno method published in DOI: 10.1093/bioinformatics/btab409
#Unfortunately this is a python package and its execution has not been incorporated in this R markdown.
#Instead, we have created a conda environment, and it can be executed outside of this R-based session
#To do so, we need to create the appropriate input files using the training and testing sets used in the above.
#The input file formats required for CyAnno are described here https://github.com/abbioinfo/CyAnno/wiki/CyAnno

#In short, for training/testing cyAnno we need the following files that have been created in the folder cytoAnno_training_example/:

# 1. training_handgated.csv
# 2. LivecellsTraining.csv
# 3. LivecellsTesting.csv

#once CyAnno is executed it produces the following file which contains the predicted labels
#Method_x__Posterior_Probability_Prediction_Result.csv
#we need to read this file and estimate the classification performance in a similar fashion as before

predictionsCyAnno <- read.csv('Method_x__Posterior_Probability_Prediction_Result.csv')

#convert the cell type names to their original format
predictionsCyAnno$CellTypeOrignal <- gsub("_", " ", predictionsCyAnno$CellTypeOrignal)
predictionsCyAnno$CellTypePredicted <- gsub("_", " ", predictionsCyAnno$CellTypePredicted)
#make the column names consistent same as before
colnames(predictionsCyAnno)[3] <- "realLabel"
colnames(predictionsCyAnno)[4] <- "predictedLabel"
performanceCyAnno <- data.frame()
for(idx in unique(testingSetLabels)){
  
  currentType <- idx
  
  tp_row <- which(predictionsCyAnno$realLabel==currentType & predictionsCyAnno$predictedLabel==currentType)
  tn_row <- which(predictionsCyAnno$realLabel!=currentType & predictionsCyAnno$predictedLabel!=currentType)
  fp_row <- which(predictionsCyAnno$realLabel!=currentType & predictionsCyAnno$predictedLabel==currentType)
  fn_row <- which(predictionsCyAnno$realLabel==currentType & predictionsCyAnno$predictedLabel!=currentType)
  
  tp_row <- length(tp_row)
  tn_row <- length(tn_row)
  fp_row <- length(fp_row)
  fn_row <- length(fn_row)
  
  a <- data.frame(Sen=tp_row/(tp_row+fn_row),
                  Spe=tn_row/(tn_row+fp_row),
                  #PPV=tp_row/(tp_row+fp_row),
                  F1=2*tp_row/(2*tp_row+fp_row+fn_row),
                  Acc=(tp_row+tn_row)/(tp_row+tn_row+fp_row+fn_row))
  performanceCyAnno <- rbind(performanceCyAnno,a)
}
performanceCyAnno$CellType <- unique(testingSetLabels)

```

### Summarising the performance of the self-consistency analysis using Healthy.rds dataset 
The results from the previous section are aggregated and visualized using boxplots to compare the relative performances.

It becomes apparent that the Scaffold and CyAnno are ranked best compared to all othe competitor methods. 

```{r self-consistency analysis visualisation,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}
  
  library(reshape2)

  performanceScaffold$Method <- 'Scaffold'
  performanceLDA$Method <- 'LDA'
  performanceKNN3$Method <- 'KNN3'
  performanceKNN10$Method <- 'KNN10'
  performanceCyAnno$Method <- 'CyAnno'
  
  plotData <- rbind(performanceScaffold,performanceLDA,performanceKNN3,performanceKNN10,performanceCyAnno)
  
  
  color_qual_flow2 <- c('dodgerblue','coral1',"#FFD258","#3F1B03", "#894F36",'#9a6324',"#F4AD31",'khaki3',
                      "#F4800C",'darkorange3' , 'brown2',"#CC242A",'firebrick4',
                      "#EF8ECC",'#e6beff','#fabebe','darkorchid4','lightpink3','mediumpurple2',
                      "#1C750C","#557F7A","#4DB23B","#066970",'#008080',"#b0c6a2",
                      '#808000',"#CBD49C",'olivedrab2','mediumaquamarine','paleturquoise1','steelblue3',
                      "#6471E2",'skyblue','#4363d8','#46f0f0','#000075','slategray3',
                      'mediumblue','cornflowerblue','#ffe119','#ffd8b1',"#FEEDC3",
                      'bisque3','#808080','black','gray88','gray','burlywood4')
  
  tmp <- melt(plotData,id.vars = c("CellType","Method" ))
  tmp$Method <- factor(tmp$Method,levels = c('Scaffold','CyAnno','LDA','KNN3','KNN10'))
  colnames(tmp)[3] <- 'Metric'
  tmp$Metric <- factor(tmp$Metric,levels = c('Acc','Sen','Spe','F1'))
  performance_sanityTest <- ggplot(tmp, aes(x = Method, y = value, fill = Method)) +
              geom_jitter(height=0,width=0.1, size=0.48,alpha=0.8)+
              geom_boxplot(width=0.28,
                           alpha=0.8,
                           size=0.3,
                    outlier.colour = "black",
                    outlier.shape = NA,
                    outlier.fill = "red",
                    outlier.size = 0.05,
                    notch = FALSE,color='black')+
              facet_wrap(~Metric,nrow = 1)+
              ggtitle('Self-consistency test performances')+
              theme_bw() +
              ylim(0,1)+
              scale_fill_manual(values = color_qual_flow2)+
              theme(axis.text.y = element_text( size = 12 ),
                      axis.text.x = element_text(angle = 90, vjust = 0.05, hjust = 0.95, size = 12),
                      axis.title.x = element_blank(),
                      axis.title.y = element_blank(),
                      strip.text = element_text(size = 12,face='bold',lineheight=1),
                      legend.position = "none",aspect.ratio = 0.8)
  print(performance_sanityTest)
  
   myfile <- paste('Sanity_check_perf.pdf',sep ='')
  pdf(myfile)
  print(performance_sanityTest)
  dev.off()
  
  
  
  tmp %>%
    group_by(Method) %>%
    dplyr::summarize(Mean = mean(value, na.rm=TRUE))
  
  
```

### Comparison analysis using independent cytof datasets
In this subsection we aim at benchmarking the Scaffold and CyAnno approaches using completely independent data sets, namely:

1) **AML_benchmark**, 
2) **BMMC_benchmark**, 
3) **Samusik_all_SE**,

This means that the training process will be perform using **Healthy_rds**, and the testing data will be completely unseen from the training.

To make the training/testing data sets compatible, we focus on a smaller number of cell types, that has been selected after manually inpsecting all datasets. It was also necessary to harmonise the names allowing automated analyses and comparisons using strings.

In this way our comparison analysis, although focuses on a smaller set of the available cell types in Healthy_rds, it provides a completely unbiased way to evaluate the effectiveness of Scaffold and cyAnno methods on annotating cells with completely unknown identify. 

```{r comparison analysis using independent datasets,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}
  
library(lsa)    
library(umap)
library(ggthemes)
library(ConsensusClusterPlus)
library(FlowSOM)
library(reshape2)

db13 <- c('#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#911eb4', '#46f0f0', '#f032e6', 
          '#bcf60c', '#fabebe', '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000', 'slategray3',
          'khaki3','bisque3','coral1','mediumaquamarine',
          '#aaffc3', '#808000', '#ffd8b1', '#000075', '#808080','gray88' ,'#ffffff', '#000000')

###############################################################################################################
##############################################################################################################
#work with AML_benchmark dataset
 AML_benchmark <- read.csv('AML_benchmark.csv')

  #1 SCAFFOLD approach
  #"Plasma cells"
  selected_cellTypes_Healthy <- c("Plasma cells","CD56dimCD16+ NK cells","HSCs & MPPs"  ,
  "CD4+ naive T cells","CD4+ cytotoxic T cells" ,"CD4+ memory T cells"  ,
  "CD8+ naive T cells","CD8+CD103+ tissue resident memory T cells","CD8+ effector memory T cells" ,"CD8+ central memory T cells",
  "Mature naive B cells","Classical Monocytes", "Non-classical monocytes" , 
  "Plasmacytoid dendritic cells","Plasmacytoid dendritic cell progenitors","Pre-B cells","Pro-B cells")  
  #"Plasma B cells"
  selected_cellTypes_AML <- c("Plasma B cells","CD16+ NK cells",
                          "CD34+CD38+CD123- HSPCs" ,"CD34+CD38+CD123+ HSPCs","CD34+CD38lo HSCs",
                          "CD4 T cells","CD8 T cells",
                          "Mature B cells",
                          "Monocytes","pDCs",
                          "Pre B cells","Pro B cells")
  
  idx <- which(AML_benchmark$cell_type %in% selected_cellTypes_AML)
  AML_benchmark <- AML_benchmark[idx,]
  idx <- which(healthy.AB.data.filtered$cellType %in% selected_cellTypes_Healthy)
  healthy.AB.data.filtered <- healthy.AB.data.filtered[idx,]
  #harmonise the cell type names
  healthy.AB.data.filtered$HarmonisedName <- NA
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="CD4+ cytotoxic T cells" | 
                            healthy.AB.data.filtered$cellType=="CD4+ memory T cells" |
                            healthy.AB.data.filtered$cellType=="CD4+ naive T cells" ,'HarmonisedName'] <- 'CD4_T_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="CD8+ central memory T cells" | 
                            healthy.AB.data.filtered$cellType=="CD8+ effector memory T cells" |
                            healthy.AB.data.filtered$cellType=="CD8+ naive T cells" | 
                            healthy.AB.data.filtered$cellType=="CD8+CD103+ tissue resident memory T cells"  ,'HarmonisedName'] <- 'CD8_T_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Classical Monocytes" | 
                            healthy.AB.data.filtered$cellType=="Non-classical monocytes",'HarmonisedName'] <- 'Monocytes'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="HSCs & MPPs",'HarmonisedName'] <- 'HSCs'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Mature naive B cells",'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Plasma cells",'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Plasmacytoid dendritic cell progenitors"|
                            healthy.AB.data.filtered$cellType=="Plasmacytoid dendritic cells"  ,'HarmonisedName'] <- 'pDCs'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Pre-B cells"  ,'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Pro-B cells"  ,'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="CD56dimCD16+ NK cells"  ,'HarmonisedName'] <- 'NK_CD16pos_cells'
  
  #do the same harmonisation for the other dataset
  AML_benchmark$HarmonisedName <- NA
  AML_benchmark[AML_benchmark$cell_type=="CD16+ NK cells"  ,'HarmonisedName'] <- 'NK_CD16pos_cells'
  AML_benchmark[ AML_benchmark$cell_type=="CD34+CD38+CD123- HSPCs" | 
                            AML_benchmark$cell_type=="CD34+CD38+CD123+ HSPCs" |
                            AML_benchmark$cell_type=="CD34+CD38lo HSCs"  ,'HarmonisedName'] <- 'HSCs'
  AML_benchmark[AML_benchmark$cell_type=="CD4 T cells"  ,'HarmonisedName'] <- 'CD4_T_cells'
  AML_benchmark[AML_benchmark$cell_type=="CD8 T cells"  ,'HarmonisedName'] <- 'CD8_T_cells'
  AML_benchmark[AML_benchmark$cell_type=="Mature B cells"  ,'HarmonisedName'] <- 'B_cells'
  AML_benchmark[AML_benchmark$cell_type=="Monocytes"  ,'HarmonisedName'] <- 'Monocytes'
  AML_benchmark[AML_benchmark$cell_type=="pDCs"  ,'HarmonisedName'] <- 'pDCs'
  AML_benchmark[AML_benchmark$cell_type=="Pre B cells"  ,'HarmonisedName'] <- 'B_cells'
  AML_benchmark[AML_benchmark$cell_type=="Pro B cells"  ,'HarmonisedName'] <- 'B_cells'
  AML_benchmark[AML_benchmark$cell_type=="Plasma B cells"  ,'HarmonisedName'] <- 'B_cells'
  
  #find common markers between the panels --> this is required for cyAnno to run, the scaffold since it computes cosine distance, it can, in principle run, with different markers
  colnames(AML_benchmark)
  colnames(healthy.AB.data.filtered)
  #after inspecting the markers we see that CD235ab is named CD235a-b in the healthy cohort, thus we need to change
  colnames(healthy.AB.data.filtered)[41] <- "CD235ab"
  #same for CD49d in the AML_benchmark is called CD49d2, we will harnonise
  colnames(AML_benchmark)[33] <- "CD49d"
  common.col.names <- intersect(colnames(AML_benchmark)[1:(ncol(AML_benchmark)-1)],colnames(healthy.AB.data.filtered)[1:(ncol(healthy.AB.data.filtered)-1)])
  
############################################################################
#generate the training data that will be used as landmarks for the scaffold
############################################################################
trainingSetLabels <- healthy.AB.data.filtered$HarmonisedName
trainingSetValues <- as.matrix(healthy.AB.data.filtered[,common.col.names])
trainingSet.tab <- tabulateCellTypes(trainingSetValues,trainingSetLabels)
  
#select randomly 15000 cells to perform UMAP
AML_benchmark$cell_id <- 1:nrow(AML_benchmark)

a <- sample(nrow(AML_benchmark))
Ncell <- 15000
asel <- a[1:Ncell]
data_sub <- AML_benchmark[asel,]
data_dimred <- data_sub[ , common.col.names]
# UMAP all together
start_time = Sys.time()
data_umap <-  umap(data_dimred, random_state=123, verbose =T)
end_time = Sys.time()
end_time-start_time
umap <- as.data.frame(data_umap$layout)
colnames(umap) <- c("UMAP1", "UMAP2")

#visualise the umap in black color to see the shape
umap_black_AML_benchmark <- ggplot(umap, aes(x = UMAP1, y = UMAP2)) +
  geom_point(size = 0.5) +
  coord_fixed(ratio = 1)+
  ggtitle('AML_benchmark - UMAP')

## plot UMAP with expression overlayed
dim(data_dimred)
data_dimred_exp <- cbind(data_dimred, umap)
data_dimred_melt <- melt(data_dimred_exp, id.vars = c("UMAP1", "UMAP2"))
data_dimred_comb  <- cbind(data_sub[,c("cell_id","HarmonisedName")], data_dimred_exp)

#visualise the UMAP and the selected cell types
umap_annot_AML_benchmark <- ggplot(data_dimred_comb, aes(x = UMAP1, y = UMAP2, color = HarmonisedName)) +
  geom_point(size = 0.5) +
  coord_fixed(ratio = 1) +
  scale_color_manual(values=db13)+
  #facet_wrap(~ Sample, ncol = 4) +
  theme_bw() +
  ggtitle('AML_benchmark - UMAP')+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        legend.position = "bottom")+
  guides(color = guide_legend(override.aes = list(size=3),nrow=2,title=""))

#estimate relative abundance of cell types in this dataset
cellFreq <- AML_benchmark %>% 
  group_by(HarmonisedName) %>% 
  summarise(samples = n())%>% 
  mutate(freq = samples/sum(samples))

 AML_benchmark_cellFreq <- ggplot(cellFreq) + aes(x=reorder(HarmonisedName,-freq),y=freq,fill=HarmonisedName)+
  geom_bar(stat='identity',width = 0.58,alpha=1,color='black',size=0.2)+
  ylab('AML_benchmark - Relative abundance')+
  theme_bw()+
  theme(axis.text.y = element_text( size = 10 ),
        axis.text.x = element_text(angle = 0, vjust = 0.5, hjust = 0.5, size = 10),
        axis.title.x = element_text( size = 12,face = 'bold' ),
        axis.title.y = element_blank(),
        strip.text = element_text(size = 12,face='bold',lineheight=1),
        legend.position = "none",aspect.ratio = 1)+
  scale_fill_manual(values = db31)+
  guides(fill = guide_legend(override.aes = list(size=3),nrow=1,title=""))+
  coord_flip()

#visualise the expression of markers in the dataset 
data <- AML_benchmark[,c('CD45','CD4','CD8','CD20','CD16','CD34','CD19','CD123','CD38')]
data <- as.matrix(data)
rng <- colQuantiles(data, probs = c(0.01, 0.99))
data01 <- t((t(data) - rng[, 1]) / (rng[, 2] - rng[, 1]))
data01[data01 < 0] <- 0
data01[data01 > 1] <- 1
test <- data.frame(data01)
test$Class <- AML_benchmark$HarmonisedName
test <- melt(test,id.vars = c('Class'))
AML_benchmark_ridgesPlot <- ggplot(test, aes(x = value, y = Class, fill = variable)) + 
  geom_density_ridges(scale=1)+
  facet_wrap(~ variable, scales = "free", nrow = 4)+
  theme_ridges() + 
  theme(legend.position = "none",
        axis.text.x = element_text(angle=45, hjust=1,size = 7),
        axis.text.y = element_text(hjust=1,size = 7),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),aspect.ratio = 1.58)+
  scale_x_continuous(breaks=c(0,0.2,0.4,0.6,0.8,1))
 

 myfile <- paste('AML_benchmark_cellFreq.pdf',sep ='')
 pdf(myfile)
 print(AML_benchmark_cellFreq)
 dev.off()

 myfile <- paste('AML_benchmark_ridgesPlot.pdf',sep ='')
 pdf(myfile)
 print(AML_benchmark_ridgesPlot)
 dev.off()
 

#the scaffold approach requires the testing data to be grouped together. In our framework this can be performed using
#unsupervised clustering with FlowSOM. The clustering itself is  automated and does not require from users to take any decision about the identify of the cells and/or how to merge clusters and form metaclusters.

data_flow <- as.matrix(AML_benchmark[,common.col.names])
# generate flowframe from the data
ff_new <- flowFrame(exprs = data_flow, desc = list(FIL = 1))

start_time = Sys.time()
# run FlowSOM (with set.seed for reproducibility)
set.seed(123)
out_fSOM <- FlowSOM::ReadInput(ff_new, transform = F, scale = F, compensate = F)
out_fSOM <- FlowSOM::BuildSOM(out_fSOM, colsToUse = common.col.names)
out_fSOM <- FlowSOM::BuildMST(out_fSOM)
labels <- out_fSOM$map$mapping[,1]
end_time = Sys.time()
end_time-start_time

# make heatmap of 100 FlowSOM clusters
heat_mat <- matrix(NA, nrow = 100, ncol = length(common.col.names))
for(i in 1:100) {
  temp_mat <- data_flow[labels == i, common.col.names]
  heat_mat[i,] <- as.matrix(apply(temp_mat, 2, function (x) mean(x, na.rm = T)))
}
rownames(heat_mat) <- paste("C",1:100,sep='')
colnames(heat_mat) <- common.col.names

#ready to produce the testing dataset for mapping to the scaffold
testingSetLabels <- labels
testingSetValues <- as.matrix(AML_benchmark[,common.col.names])
testingSet.tab <- tabulateCellTypes(testingSetValues,testingSetLabels)
m <- as.matrix(testingSet.tab[, common.col.names])
rownames(m) <- testingSet.tab$cellType

############################################################################
#prepare the training data that will be used as landmarks for the scaffold
############################################################################
trainingSetLabels <- healthy.AB.data.filtered$HarmonisedName
trainingSetValues <- as.matrix(healthy.AB.data.filtered[,common.col.names])
trainingSet.tab <- tabulateCellTypes(trainingSetValues,trainingSetLabels)

#perform the same computation for the training set
att <- as.matrix(trainingSet.tab[,common.col.names])
row.names(att) <- trainingSet.tab$cellType
#compute the distance from the testing set to the training set which is the landmark using cosine similarity
dd_controls_to_landmarks <- t(apply(m, 1, function(x, att) {cosine_similarity_from_matrix(x, att)}, att))
#careful here, columns of the distance matrix are the rows of the landmark matrix
colnames(dd_controls_to_landmarks) <- rownames(att)

predictionsMapping <- data.frame()
#parse the distance matrix and find the max value per row
for(idx in 1:nrow(dd_controls_to_landmarks)){
  
  a <- which.max(dd_controls_to_landmarks[idx,])
  predictedLabel <- colnames(dd_controls_to_landmarks)[a]
  dt <- data.frame(clusterID=rownames(dd_controls_to_landmarks)[idx],
                   predictedLabel=predictedLabel,
                   Dist=dd_controls_to_landmarks[idx,a])
  predictionsMapping <- rbind(predictionsMapping,dt)
}
predictionsMapping$clusterID <- paste("C",predictionsMapping$clusterID,sep='')

#combine the predictions and the real labels of cell types
predictionsScaffold <- data.frame(clusterID=paste("C",labels,sep=''),
                                          realLabel=AML_benchmark$HarmonisedName)
predictionsScaffold$predictedLabel <- 'Unknown'
for(idx in 1:nrow(predictionsMapping)){
  m <- predictionsMapping[idx,'clusterID']
  predictionsScaffold[predictionsScaffold$clusterID==m,'predictedLabel'] <- predictionsMapping[idx,'predictedLabel']
}
#estimate the performance
performanceScaffold <- data.frame()
cellTypesToCheck <- unique(trainingSetLabels)

for(idx in cellTypesToCheck){
  
  currentType <- idx
  
  tp_row <- which(predictionsScaffold$realLabel==currentType & predictionsScaffold$predictedLabel==currentType)
  tn_row <- which(predictionsScaffold$realLabel!=currentType & predictionsScaffold$predictedLabel!=currentType)
  fp_row <- which(predictionsScaffold$realLabel!=currentType & predictionsScaffold$predictedLabel==currentType)
  fn_row <- which(predictionsScaffold$realLabel==currentType & predictionsScaffold$predictedLabel!=currentType)
  
  tp_row <- length(tp_row)
  tn_row <- length(tn_row)
  fp_row <- length(fp_row)
  fn_row <- length(fn_row)
  
  a <- data.frame(Sen=tp_row/(tp_row+fn_row),
                  Spe=tn_row/(tn_row+fp_row),
                  #PPV=tp_row/(tp_row+fp_row),
                  F1=2*tp_row/(2*tp_row+fp_row+fn_row),
                  Acc=(tp_row+tn_row)/(tp_row+tn_row+fp_row+fn_row))
  performanceScaffold <- rbind(performanceScaffold,a)
}
performanceScaffold$CellType <- cellTypesToCheck

#APPROACH #2: CytoAnno
#For training we need the following files that have been created in the folder cytoAnno_training_example/:
# 1. training_handgated.csv
# 2. LivecellsTraining.csv
# 3. LivecellsTesting.csv
#below we use lines of code to convert the appropriate training/testing data into the CyAnno format:
#tmp <- data.frame(testingSetValues)
#tmp$labels <- testingSetLabels
#write.csv(tmp,file='cytoAnno_AML_benchmark/live/Testing.csv',row.names=F,col.names = T)
#tmp <- data.frame(trainingSetValues)
#tmp$labels <- trainingSetLabels
#write.csv(tmp,file='cytoAnno_AML_benchmark/live/Training.csv',row.names=F,col.names = T)
#for(idx in unique(testingSetLabels)){
#  a <- which(trainingSetLabels==idx)
#  tmp1 <- data.frame(trainingSetValues[a,])
#  tmp2 <- trainingSetLabels[a]
#  str <- paste('cytoAnno_AML_benchmark/manual/',idx,'.csv',sep='')
#  write.csv(tmp1,file=str,row.names=F,col.names = T)
#}      

#once CyAnno is executed it produces the following file which contains the predicted labels
#Method_x__Posterior_Probability_Prediction_Result.csv --> renamed to AML_benchmarkr__Posterior_Probability_Prediction_Result.csv
#we need to read this file and estimate the classification performance same as before

predictionsCyAnno <- read.csv('AML_benchmark__Posterior_Probability_Prediction_Result.csv')

testingSetLabels <- AML_benchmark$HarmonisedName
testingSetValues <- as.matrix(AML_benchmark[,common.col.names])

#make the column names consistent same as before
colnames(predictionsCyAnno)[3] <- "realLabel"
colnames(predictionsCyAnno)[4] <- "predictedLabel"
performanceCyAnno <- data.frame()
for(idx in unique(testingSetLabels)){
  
 currentType <- idx
  
  tp_row <- which(predictionsCyAnno$realLabel==currentType & predictionsCyAnno$predictedLabel==currentType)
  tn_row <- which(predictionsCyAnno$realLabel!=currentType & predictionsCyAnno$predictedLabel!=currentType)
  fp_row <- which(predictionsCyAnno$realLabel!=currentType & predictionsCyAnno$predictedLabel==currentType)
  fn_row <- which(predictionsCyAnno$realLabel==currentType & predictionsCyAnno$predictedLabel!=currentType)
  
  tp_row <- length(tp_row)
  tn_row <- length(tn_row)
  fp_row <- length(fp_row)
  fn_row <- length(fn_row)
  
  a <- data.frame(Sen=tp_row/(tp_row+fn_row),
                  Spe=tn_row/(tn_row+fp_row),
                  #PPV=tp_row/(tp_row+fp_row),
                  F1=2*tp_row/(2*tp_row+fp_row+fn_row),
                  Acc=(tp_row+tn_row)/(tp_row+tn_row+fp_row+fn_row))
  performanceCyAnno <- rbind(performanceCyAnno,a)
}
performanceCyAnno$CellType <- unique(testingSetLabels)

  #visualise the relative performance of both methods
  performanceScaffold$Method <- 'Scaffold'
  performanceCyAnno$Method <- 'CyAnno'
  
  plotData <- rbind(performanceScaffold,performanceCyAnno)
  
  color_qual_flow2 <- c('dodgerblue','coral1',"#FFD258","#3F1B03", "#894F36",'#9a6324',"#F4AD31",'khaki3',
                      "#F4800C",'darkorange3' , 'brown2',"#CC242A",'firebrick4',
                      "#EF8ECC",'#e6beff','#fabebe','darkorchid4','lightpink3','mediumpurple2',
                      "#1C750C","#557F7A","#4DB23B","#066970",'#008080',"#b0c6a2",
                      '#808000',"#CBD49C",'olivedrab2','mediumaquamarine','paleturquoise1','steelblue3',
                      "#6471E2",'skyblue','#4363d8','#46f0f0','#000075','slategray3',
                      'mediumblue','cornflowerblue','#ffe119','#ffd8b1',"#FEEDC3",
                      'bisque3','#808080','black','gray88','gray','burlywood4')
  
  tmp <- melt(plotData,id.vars = c("CellType","Method" ))
  tmp$Method <- factor(tmp$Method,levels = c('Scaffold','CyAnno'))
  colnames(tmp)[3] <- 'Metric'
  tmp$Metric <- factor(tmp$Metric,levels = c('Acc','Sen','Spe','F1'))
  comparison_AML_benchmark <- ggplot(tmp, aes(x = Method, y = value, fill = Method)) +
                   geom_jitter(height=0,width=0.1, size=0.48,alpha=0.8)+
                     geom_boxplot(width=0.28,
                           alpha=0.8,
                           size=0.3,
                          outlier.colour = "black",
                          outlier.shape = NA,
                          outlier.fill = "red",
                          outlier.size = 0.05,
                          notch = FALSE,color='black')+
                facet_wrap(~Metric,nrow = 1)+
                ggtitle('AML_benchmark performance')+
                theme_bw() +
                ylim(0,1)+
                scale_fill_manual(values = color_qual_flow2)+
                theme(axis.text.y = element_text( size = 12 ),
                      axis.text.x = element_text(angle = 90, vjust = 0.05, hjust = 0.95, size = 12),
                      axis.title.x = element_blank(),
                      axis.title.y = element_blank(),
                      strip.text = element_text(size = 12,face='bold',lineheight=1),
                      legend.position = "none",aspect.ratio = 0.8)

  
  print(comparison_AML_benchmark)

   myfile <- paste('comparison_AML_benchmark.pdf',sep ='')
   pdf(myfile)
   print(comparison_AML_benchmark)
   dev.off()
  
###############################################################################################################
###############################################################################################################

  #work with the BMMC benchmark dataset
  BMMC_benchmark <- read.csv('BMMC_benchmark.csv')
  
  #1 SCAFFOLD approach
  #from the previous analysis the Healthy_rds dataset has been modified and thus we need to load it and process it again
  filename <- 'Healthy.rds'
  o <- loadReferenceData(filename)
  healthy.AB.data <- o[[1]]
  healthyCellTypes <- o[[2]]
  rm(o)
  
  umap <- extractRefUMAPCoordinates(filename)
  myCellCutoff <- 300
  o <- filterAndTabulate(healthy.AB.data,healthyCellTypes,umap,myCellCutoff)
  healthy.AB.data.tab <- o[[1]]
  healthy.AB.data.filtered <- o[[2]]
  umap.filtered <- o[[3]]
  rm(o)

  selected_cellTypes_Healthy <- c("Plasma cells","CD56dimCD16+ NK cells","HSCs & MPPs"  ,
  "CD4+ naive T cells","CD4+ cytotoxic T cells" ,"CD4+ memory T cells"  ,
  "CD8+ naive T cells","CD8+CD103+ tissue resident memory T cells","CD8+ effector memory T cells" ,"CD8+ central memory T cells",
  "Mature naive B cells","Classical Monocytes", "Non-classical monocytes" , 
  "Plasmacytoid dendritic cells","Plasmacytoid dendritic cell progenitors","Pre-B cells","Pro-B cells")  
  
  selected_cellTypes_BMMC <- c("CD11b- Monocyte","CD11bhi Monocyte","CD11bmid Monocyte",
                               "CMP","GMP","HSC","MEP","MEP",
                               "Immature B","Mature CD38lo B","Mature CD38mid B","Plasma cell","Pre-B I","Pre-B II",
                               "Mature CD4+ T","Naive CD4+ T",
                               "Mature CD8+ T","Naive CD8+ T",
                               "NK","Plasmacytoid DC")
  
  idx <- which(BMMC_benchmark$cell_type %in% selected_cellTypes_BMMC)
  BMMC_benchmark <- BMMC_benchmark[idx,]
  
  idx <- which(healthy.AB.data.filtered$cellType %in% selected_cellTypes_Healthy)
  healthy.AB.data.filtered <- healthy.AB.data.filtered[idx,]
  
  #harmonise the cell type names
  healthy.AB.data.filtered$HarmonisedName <- NA
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="CD4+ cytotoxic T cells" | 
                            healthy.AB.data.filtered$cellType=="CD4+ memory T cells" |
                            healthy.AB.data.filtered$cellType=="CD4+ naive T cells" ,'HarmonisedName'] <- 'CD4_T_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="CD8+ central memory T cells" | 
                            healthy.AB.data.filtered$cellType=="CD8+ effector memory T cells" |
                            healthy.AB.data.filtered$cellType=="CD8+ naive T cells" | 
                            healthy.AB.data.filtered$cellType=="CD8+CD103+ tissue resident memory T cells"  ,'HarmonisedName'] <- 'CD8_T_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Classical Monocytes" | 
                            healthy.AB.data.filtered$cellType=="Non-classical monocytes",'HarmonisedName'] <- 'Monocytes'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="HSCs & MPPs",'HarmonisedName'] <- 'HSCs_MPPs'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Mature naive B cells",'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Plasma cells",'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Plasmacytoid dendritic cell progenitors"|
                            healthy.AB.data.filtered$cellType=="Plasmacytoid dendritic cells"  ,'HarmonisedName'] <- 'pDCs'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Pre-B cells"  ,'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Pro-B cells"  ,'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="CD56dimCD16+ NK cells"  ,'HarmonisedName'] <- 'NK_CD16pos_cells'
  
  #do the same harmonisation for the other dataset
  BMMC_benchmark$HarmonisedName <- NA
  BMMC_benchmark[BMMC_benchmark$cell_type=="NK"  ,'HarmonisedName'] <- 'NK_CD16pos_cells'
  BMMC_benchmark[BMMC_benchmark$cell_type=="CMP" | 
                            BMMC_benchmark$cell_type=="GMP" |
                            BMMC_benchmark$cell_type=="HSC" |
                            BMMC_benchmark$cell_type=="MEP",'HarmonisedName'] <- 'HSCs_MPPs'
  BMMC_benchmark[BMMC_benchmark$cell_type=="Mature CD4+ T" |
                BMMC_benchmark$cell_type=="Naive CD4+ T",'HarmonisedName'] <- 'CD4_T_cells'
  BMMC_benchmark[BMMC_benchmark$cell_type=="Mature CD8+ T" |
                BMMC_benchmark$cell_type=="Naive CD8+ T",'HarmonisedName'] <- 'CD8_T_cells'
  BMMC_benchmark[BMMC_benchmark$cell_type=="Mature CD38lo B" |
                BMMC_benchmark$cell_type=="Mature CD38mid B" |
                BMMC_benchmark$cell_type=="Plasma cell" |
                BMMC_benchmark$cell_type=="Immature B" |
                BMMC_benchmark$cell_type=="Pre-B I" |
                BMMC_benchmark$cell_type=="Pre-B II"  ,'HarmonisedName'] <- 'B_cells'
  BMMC_benchmark[BMMC_benchmark$cell_type=="CD11b- Monocyte" |
                BMMC_benchmark$cell_type=="CD11bhi Monocyte" |
                BMMC_benchmark$cell_type=="CD11bmid Monocyte"  ,'HarmonisedName'] <- 'Monocytes'
  BMMC_benchmark[BMMC_benchmark$cell_type=="Plasmacytoid DC"  ,'HarmonisedName'] <- 'pDCs'
   
  #find common markers between the panels --> this is required for cyAnno to run, the scaffold since it computes cosine distance can in principle run with different markers
  colnames(BMMC_benchmark)
  colnames(healthy.AB.data.filtered)

  common.col.names <- intersect(colnames(BMMC_benchmark)[1:(ncol(BMMC_benchmark)-1)],colnames(healthy.AB.data.filtered)[1:(ncol(healthy.AB.data.filtered)-1)])
  
############################################################################
#generate the training data that will be used as landmarks for the scaffold
############################################################################
trainingSetLabels <- healthy.AB.data.filtered$HarmonisedName
trainingSetValues <- as.matrix(healthy.AB.data.filtered[,common.col.names])
trainingSet.tab <- tabulateCellTypes(trainingSetValues,trainingSetLabels)
  
#select randomly 15000 cells to perform UMAP
BMMC_benchmark$cell_id <- 1:nrow(BMMC_benchmark)

a <- sample(nrow(BMMC_benchmark))
Ncell <- 15000
asel <- a[1:Ncell]
data_sub <- BMMC_benchmark[asel,]
data_dimred <- data_sub[ , common.col.names]
# UMAP all together
start_time = Sys.time()
data_umap <-  umap(data_dimred, random_state=123, verbose =T)
end_time = Sys.time()
end_time-start_time
umap <- as.data.frame(data_umap$layout)
colnames(umap) <- c("UMAP1", "UMAP2")

#visualise the umap in black color to see the shape
umap_black_BMCC <- ggplot(umap, aes(x = UMAP1, y = UMAP2)) +
  geom_point(size = 0.5) +
  coord_fixed(ratio = 1)

## plot UMAP with expression overlayed
dim(data_dimred)
data_dimred_exp <- cbind(data_dimred, umap)
data_dimred_melt <- melt(data_dimred_exp, id.vars = c("UMAP1", "UMAP2"))
data_dimred_comb  <- cbind(data_sub[,c("cell_id","HarmonisedName")], data_dimred_exp)

#visualise the UMAP and the selected cell types
umap_annot_BMMC_benchmark <- ggplot(data_dimred_comb, aes(x = UMAP1, y = UMAP2, color = HarmonisedName)) +
  geom_point(size = 0.5) +
  coord_fixed(ratio = 1) +
  scale_color_manual(values=db13)+
  #facet_wrap(~ Sample, ncol = 4) +
  theme_bw() +
  ggtitle('BMMC_benchmark - UMAP')+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        legend.position = "bottom")+
  guides(color = guide_legend(override.aes = list(size=3),nrow=2,title=""))

#estimate relative abundance of cell types in this dataset
cellFreq <- BMMC_benchmark %>% 
  group_by(HarmonisedName) %>% 
  summarise(samples = n())%>% 
  mutate(freq = samples/sum(samples))

BMMC_benchmark_cellFreq <- ggplot(cellFreq) + aes(x=reorder(HarmonisedName,-freq),y=freq,fill=HarmonisedName)+
  geom_bar(stat='identity',width = 0.58,alpha=1,color='black',size=0.2)+
  ylab('BMMC_benchmark - Relative abundance')+
  theme_bw()+
  theme(axis.text.y = element_text( size = 10 ),
        axis.text.x = element_text(angle = 0, vjust = 0.5, hjust = 0.5, size = 10),
        axis.title.x = element_text( size = 12,face = 'bold' ),
        axis.title.y = element_blank(),
        strip.text = element_text(size = 12,face='bold',lineheight=1),
        legend.position = "none",aspect.ratio = 1)+
  scale_fill_manual(values = db31)+
  guides(fill = guide_legend(override.aes = list(size=3),nrow=1,title=""))+
  coord_flip()

#visualise the expression of markers in the dataset 
data <- BMMC_benchmark[,c('CD45','CD4','CD8','CD20','CD33','CD34','CD19','CD123','CD38')]
data <- as.matrix(data)
rng <- colQuantiles(data, probs = c(0.01, 0.99))
data01 <- t((t(data) - rng[, 1]) / (rng[, 2] - rng[, 1]))
data01[data01 < 0] <- 0
data01[data01 > 1] <- 1
test <- data.frame(data01)
test$Class <- BMMC_benchmark$HarmonisedName
test <- melt(test,id.vars = c('Class'))
BMMC_benchmark_ridgesPlot <- ggplot(test, aes(x = value, y = Class, fill = variable)) + 
  geom_density_ridges(scale=1)+
  facet_wrap(~ variable, scales = "free", nrow = 4)+
  theme_ridges() + 
  theme(legend.position = "none",
        axis.text.x = element_text(angle=45, hjust=1,size = 7),
        axis.text.y = element_text(hjust=1,size = 7),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),aspect.ratio = 1.58)+
  scale_x_continuous(breaks=c(0,0.2,0.4,0.6,0.8,1))

 myfile <- paste('BMMC_benchmark_cellFreq.pdf',sep ='')
 pdf(myfile)
 print(BMMC_benchmark_cellFreq)
 dev.off()

 myfile <- paste('BMMC_benchmark_ridgesPlot.pdf',sep ='')
 pdf(myfile)
 print(BMMC_benchmark_ridgesPlot)
 dev.off()


#clustering with flowSOM
data_flow <- as.matrix(BMMC_benchmark[,common.col.names])
# generate flowframe from the data
ff_new <- flowFrame(exprs = data_flow, desc = list(FIL = 1))

start_time = Sys.time()
# run FlowSOM (with set.seed for reproducibility)
set.seed(123)
out_fSOM <- FlowSOM::ReadInput(ff_new, transform = F, scale = F, compensate = F)
out_fSOM <- FlowSOM::BuildSOM(out_fSOM, colsToUse = common.col.names)
out_fSOM <- FlowSOM::BuildMST(out_fSOM)
labels <- out_fSOM$map$mapping[,1]
end_time = Sys.time()
end_time-start_time

# make heatmap of 100 FlowSOM clusters
heat_mat <- matrix(NA, nrow = 100, ncol = length(common.col.names))
for(i in 1:100) {
  temp_mat <- data_flow[labels == i, common.col.names]
  if(is.vector(temp_mat)){
    heat_mat[i,] <- as.matrix(t(temp_mat))
  }else{
     heat_mat[i,] <- as.matrix(apply(temp_mat, 2, function (x) mean(x, na.rm = T)))
  }
  
}
rownames(heat_mat) <- paste("C",1:100,sep='')
colnames(heat_mat) <- common.col.names

#produce the testing dataset for scaffold mapping
testingSetLabels <- labels
testingSetValues <- as.matrix(BMMC_benchmark[,common.col.names])
testingSet.tab <- tabulateCellTypes(testingSetValues,testingSetLabels)
m <- as.matrix(testingSet.tab[, common.col.names])
rownames(m) <- testingSet.tab$cellType

############################################################################
#prepare the training data that will be used as landmarks for the scaffold
############################################################################
trainingSetLabels <- healthy.AB.data.filtered$HarmonisedName
trainingSetValues <- as.matrix(healthy.AB.data.filtered[,common.col.names])
trainingSet.tab <- tabulateCellTypes(trainingSetValues,trainingSetLabels)

#perform the same computation for the training set
att <- as.matrix(trainingSet.tab[,common.col.names])
row.names(att) <- trainingSet.tab$cellType
#compute the distance from the testing set to the training set which is the landmark using cosine similarity
dd_controls_to_landmarks <- t(apply(m, 1, function(x, att) {cosine_similarity_from_matrix(x, att)}, att))
#careful here, columns of the distance matrix are the rows of the landmark matrix
colnames(dd_controls_to_landmarks) <- rownames(att)

predictionsMapping <- data.frame()
#parse the distance matrix and find the max value per row
for(idx in 1:nrow(dd_controls_to_landmarks)){
  a <- which.max(dd_controls_to_landmarks[idx,])
  predictedLabel <- colnames(dd_controls_to_landmarks)[a]
  dt <- data.frame(clusterID=rownames(dd_controls_to_landmarks)[idx],
                   predictedLabel=predictedLabel,
                   Dist=dd_controls_to_landmarks[idx,a])
  predictionsMapping <- rbind(predictionsMapping,dt)
}
predictionsMapping$clusterID <- paste("C",predictionsMapping$clusterID,sep='')

#combine the predictions and the real information about cell types
predictionsScaffold <- data.frame(clusterID=paste("C",labels,sep=''),
                                          realLabel=BMMC_benchmark$HarmonisedName)
predictionsScaffold$predictedLabel <- 'Unknown'
for(idx in 1:nrow(predictionsMapping)){
  m <- predictionsMapping[idx,'clusterID']
  predictionsScaffold[predictionsScaffold$clusterID==m,'predictedLabel'] <- predictionsMapping[idx,'predictedLabel']
}

performanceScaffold <- data.frame()
cellTypesToCheck <- unique(trainingSetLabels)
for(idx in cellTypesToCheck){
  
 currentType <- idx
  
  tp_row <- which(predictionsScaffold$realLabel==currentType & predictionsScaffold$predictedLabel==currentType)
  tn_row <- which(predictionsScaffold$realLabel!=currentType & predictionsScaffold$predictedLabel!=currentType)
  fp_row <- which(predictionsScaffold$realLabel!=currentType & predictionsScaffold$predictedLabel==currentType)
  fn_row <- which(predictionsScaffold$realLabel==currentType & predictionsScaffold$predictedLabel!=currentType)
  
  tp_row <- length(tp_row)
  tn_row <- length(tn_row)
  fp_row <- length(fp_row)
  fn_row <- length(fn_row)
  
  a <- data.frame(Sen=tp_row/(tp_row+fn_row),
                  Spe=tn_row/(tn_row+fp_row),
                  #PPV=tp_row/(tp_row+fp_row),
                  F1=2*tp_row/(2*tp_row+fp_row+fn_row),
                  Acc=(tp_row+tn_row)/(tp_row+tn_row+fp_row+fn_row))
  performanceScaffold <- rbind(performanceScaffold,a)
}
performanceScaffold$CellType <- cellTypesToCheck

#APPROACH #2: CytoAnno
#below we use lines of code to convert the appropriate training/testing data into the CyAnno format:
#tmp <- data.frame(testingSetValues)
#tmp$labels <- testingSetLabels
#write.csv(tmp,file='cytoAnno_BMMC_benchmark/live/Testing.csv',row.names=F,col.names = T)

#tmp <- data.frame(trainingSetValues)
#tmp$labels <- trainingSetLabels
#write.csv(tmp,file='cytoAnno_BMMC_benchmark/live/Training.csv',row.names=F,col.names = T)

#for(idx in unique(testingSetLabels)){
#  a <- which(trainingSetLabels==idx)
#  tmp1 <- data.frame(trainingSetValues[a,])
#  tmp2 <- trainingSetLabels[a]
#  str <- paste('cytoAnno_BMMC_benchmark/manual/',idx,'.csv',sep='')
#  write.csv(tmp1,file=str,row.names=F,col.names = T)
#}      

#once CyAnno is executed it produces the following file which contains the predicted labels
#Method_x__Posterior_Probability_Prediction_Result.csv --> renamed to BMMC_benchmarkr__Posterior_Probability_Prediction_Result.csv
#we need to read this file and estimate the classification performance similar to the previous approaches

predictionsCyAnno <- read.csv('BMMC_benchmark__Posterior_Probability_Prediction_Result.csv')
testingSetLabels <- BMMC_benchmark$HarmonisedName
testingSetValues <- as.matrix(BMMC_benchmark[,common.col.names])

#make the column names consistent same as before
colnames(predictionsCyAnno)[3] <- "realLabel"
colnames(predictionsCyAnno)[4] <- "predictedLabel"
performanceCyAnno <- data.frame()
for(idx in unique(testingSetLabels)){
  
  currentType <- idx
  
  tp_row <- which(predictionsCyAnno$realLabel==currentType & predictionsCyAnno$predictedLabel==currentType)
  tn_row <- which(predictionsCyAnno$realLabel!=currentType & predictionsCyAnno$predictedLabel!=currentType)
  fp_row <- which(predictionsCyAnno$realLabel!=currentType & predictionsCyAnno$predictedLabel==currentType)
  fn_row <- which(predictionsCyAnno$realLabel==currentType & predictionsCyAnno$predictedLabel!=currentType)
  
  tp_row <- length(tp_row)
  tn_row <- length(tn_row)
  fp_row <- length(fp_row)
  fn_row <- length(fn_row)
  
  a <- data.frame(Sen=tp_row/(tp_row+fn_row),
                  Spe=tn_row/(tn_row+fp_row),
                  #PPV=tp_row/(tp_row+fp_row),
                  F1=2*tp_row/(2*tp_row+fp_row+fn_row),
                  Acc=(tp_row+tn_row)/(tp_row+tn_row+fp_row+fn_row))
  performanceCyAnno <- rbind(performanceCyAnno,a)
}
performanceCyAnno$CellType <- unique(testingSetLabels)
  
#visualise the relative performance of both methods
  performanceScaffold$Method <- 'Scaffold'
  performanceCyAnno$Method <- 'CyAnno'
  
  plotData <- rbind(performanceScaffold,performanceCyAnno)
  
  tmp <- melt(plotData,id.vars = c("CellType","Method" ))
  tmp$Method <- factor(tmp$Method,levels = c('Scaffold','CyAnno'))
  colnames(tmp)[3] <- 'Metric'
  tmp$Metric <- factor(tmp$Metric,levels = c('Acc','Sen','Spe','F1'))
  comparison_BMMC_benchmark <- ggplot(tmp, aes(x = Method, y = value, fill = Method)) +
                   geom_jitter(height=0,width=0.1, size=0.48,alpha=0.8)+
                     geom_boxplot(width=0.28,
                           alpha=0.8,
                           size=0.3,
                          outlier.colour = "black",
                          outlier.shape = NA,
                          outlier.fill = "red",
                          outlier.size = 0.05,
                          notch = FALSE,color='black')+
                facet_wrap(~Metric,nrow = 1)+
                ggtitle('BMMC_benchmark performance')+
                theme_bw() +
                ylim(0,1)+
                scale_fill_manual(values = color_qual_flow2)+
                theme(axis.text.y = element_text( size = 12 ),
                      axis.text.x = element_text(angle = 90, vjust = 0.05, hjust = 0.95, size = 12),
                      axis.title.x = element_blank(),
                      axis.title.y = element_blank(),
                      strip.text = element_text(size = 12,face='bold',lineheight=1),
                      legend.position = "none",aspect.ratio = 0.8)

   myfile <- paste('comparison_BMMC_benchmark.pdf',sep ='')
   pdf(myfile)
   print(comparison_BMMC_benchmark)
   dev.off()

  
###############################################################################################################v
###############################################################################################################

#work with the PANORAMA benchmark dataset
d_SE <-  Samusik_01_SE(metadata = FALSE)
#fetch the expression data  
PANORAMA_benchmark <- assay(d_SE)
population <- as.character(rowData(d_SE)$population_id)

# transform data using asinh with cofactor 5
cofactor <- 5
PANORAMA_benchmark <- asinh(PANORAMA_benchmark / cofactor)

#combine together with the labels
PANORAMA_benchmark <- as.data.frame(PANORAMA_benchmark)
PANORAMA_benchmark$cell_type <- population

 #1 SCAFFOLD approach
  #from the previous analysis the Healthy_rds dataset has been modified and thus we need to load it again
  filename <- 'Healthy.rds'
  o <- loadReferenceData(filename)
  healthy.AB.data <- o[[1]]
  healthyCellTypes <- o[[2]]
  rm(o)
  umap <- extractRefUMAPCoordinates(filename)
  myCellCutoff <- 300
  o <- filterAndTabulate(healthy.AB.data,healthyCellTypes,umap,myCellCutoff)
  healthy.AB.data.tab <- o[[1]]
  healthy.AB.data.filtered <- o[[2]]
  umap.filtered <- o[[3]]
  rm(o)

  selected_cellTypes_Healthy <- c("Plasma cells","CD56dimCD16+ NK cells","HSCs & MPPs"  ,
  "CD4+ naive T cells","CD4+ cytotoxic T cells" ,"CD4+ memory T cells"  ,
  "CD8+ naive T cells","CD8+CD103+ tissue resident memory T cells","CD8+ effector memory T cells" ,"CD8+ central memory T cells",
  "Mature naive B cells","Classical Monocytes", "Non-classical monocytes" , 
  "Plasmacytoid dendritic cells","Plasmacytoid dendritic cell progenitors","Pre-B cells","Pro-B cells")  
  
  selected_cellTypes_PANORAMA <- c("B-cell Frac A-C (pro-B cells)","CD4 T cells","CD8 T cells",
                                   "Classical Monocytes","CMP","GMP","HSC","IgD- IgMpos B cells",
                                   "IgDpos IgMpos B cells","IgM- IgD- B-cells","Intermediate Monocytes","MEP","MPP",
                                   "NK cells","Non-Classical Monocytes","pDCs","Plasma Cells")
  
  idx <- which(PANORAMA_benchmark$cell_type %in% selected_cellTypes_PANORAMA)
  PANORAMA_benchmark <- PANORAMA_benchmark[idx,]
  
  idx <- which(healthy.AB.data.filtered$cellType %in% selected_cellTypes_Healthy)
  healthy.AB.data.filtered <- healthy.AB.data.filtered[idx,]
  
  #harmonise the cell type names
  healthy.AB.data.filtered$HarmonisedName <- NA
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="CD4+ cytotoxic T cells" | 
                            healthy.AB.data.filtered$cellType=="CD4+ memory T cells" |
                            healthy.AB.data.filtered$cellType=="CD4+ naive T cells" ,'HarmonisedName'] <- 'CD4_T_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="CD8+ central memory T cells" | 
                            healthy.AB.data.filtered$cellType=="CD8+ effector memory T cells" |
                            healthy.AB.data.filtered$cellType=="CD8+ naive T cells" | 
                            healthy.AB.data.filtered$cellType=="CD8+CD103+ tissue resident memory T cells"  ,'HarmonisedName'] <- 'CD8_T_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Classical Monocytes" | 
                            healthy.AB.data.filtered$cellType=="Non-classical monocytes",'HarmonisedName'] <- 'Monocytes'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="HSCs & MPPs",'HarmonisedName'] <- 'HSCs_MPPs'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Mature naive B cells",'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Plasma cells",'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Plasmacytoid dendritic cell progenitors"|
                            healthy.AB.data.filtered$cellType=="Plasmacytoid dendritic cells"  ,'HarmonisedName'] <- 'pDCs'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Pre-B cells"  ,'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="Pro-B cells"  ,'HarmonisedName'] <- 'B_cells'
  
  healthy.AB.data.filtered[ healthy.AB.data.filtered$cellType=="CD56dimCD16+ NK cells"  ,'HarmonisedName'] <- 'NK_CD16pos_cells'
  
  #do the same harmonisation for the other dataset
  PANORAMA_benchmark$HarmonisedName <- NA
  PANORAMA_benchmark[PANORAMA_benchmark$cell_type=="NK cells"  ,'HarmonisedName'] <- 'NK_CD16pos_cells'
  PANORAMA_benchmark[PANORAMA_benchmark$cell_type=="CMP" | 
                            PANORAMA_benchmark$cell_type=="CMP" |
                            PANORAMA_benchmark$cell_type=="GMP" |
                            PANORAMA_benchmark$cell_type=="HSC" | 
                            PANORAMA_benchmark$cell_type=="MPP" |
                            PANORAMA_benchmark$cell_type=="MEP",'HarmonisedName'] <- 'HSCs_MPPs'
  PANORAMA_benchmark[PANORAMA_benchmark$cell_type=="CD4 T cells",'HarmonisedName'] <- 'CD4_T_cells'
  PANORAMA_benchmark[PANORAMA_benchmark$cell_type=="CD8 T cells",'HarmonisedName'] <- 'CD8_T_cells'
  PANORAMA_benchmark[PANORAMA_benchmark$cell_type=="B-cell Frac A-C (pro-B cells)" |
                PANORAMA_benchmark$cell_type=="IgD- IgMpos B cells" |
                PANORAMA_benchmark$cell_type=="IgDpos IgMpos B cells" |
                PANORAMA_benchmark$cell_type=="IgM- IgD- B-cells" |
                PANORAMA_benchmark$cell_type=="Plasma Cells",'HarmonisedName'] <- 'B_cells'
  PANORAMA_benchmark[PANORAMA_benchmark$cell_type=="Classical Monocytes" |
                PANORAMA_benchmark$cell_type=="Intermediate Monocytes" |
                PANORAMA_benchmark$cell_type=="Non-Classical Monocytes"  ,'HarmonisedName'] <- 'Monocytes'
  PANORAMA_benchmark[PANORAMA_benchmark$cell_type=="pDCs"  ,'HarmonisedName'] <- 'pDCs'
   
  #find common markers between the panels --> this is required for cyAnno to run, the scaffold since it computes cosine distance can in principle run with different markers
  colnames(PANORAMA_benchmark)
  colnames(PANORAMA_benchmark)[10] <- 'CD45'
  colnames(PANORAMA_benchmark)[27] <- 'CD16'
  colnames(healthy.AB.data.filtered)

  common.col.names <- intersect(colnames(PANORAMA_benchmark)[1:(ncol(PANORAMA_benchmark)-1)],colnames(healthy.AB.data.filtered)[1:(ncol(healthy.AB.data.filtered)-1)])
  
############################################################################
#generate the training data that will be used as landmarks for the scaffold
############################################################################
trainingSetLabels <- healthy.AB.data.filtered$HarmonisedName
trainingSetValues <- as.matrix(healthy.AB.data.filtered[,common.col.names])
trainingSet.tab <- tabulateCellTypes(trainingSetValues,trainingSetLabels)
  
#select randomly 15000 cells to perform UMAP
PANORAMA_benchmark$cell_id <- 1:nrow(PANORAMA_benchmark)

a <- sample(nrow(PANORAMA_benchmark))
Ncell <- 15000
asel <- a[1:Ncell]
data_sub <- PANORAMA_benchmark[asel,]
data_dimred <- data_sub[ , common.col.names]
# UMAP all together
start_time = Sys.time()
data_umap <-  umap(data_dimred, random_state=123, verbose =T)
end_time = Sys.time()
end_time-start_time
umap <- as.data.frame(data_umap$layout)
colnames(umap) <- c("UMAP1", "UMAP2")

#visualise the umap in black color to see the shape
umap_black_PANORAMA_benchmark <- ggplot(umap, aes(x = UMAP1, y = UMAP2)) +
  geom_point(size = 0.5) +
  coord_fixed(ratio = 1)+
  ggtitle('PANORAMA_benchmark - UMAP')

## plot UMAP with expression overlayed
dim(data_dimred)
data_dimred_exp <- cbind(data_dimred, umap)
data_dimred_melt <- melt(data_dimred_exp, id.vars = c("UMAP1", "UMAP2"))
data_dimred_comb  <- cbind(data_sub[,c("cell_id","HarmonisedName")], data_dimred_exp)

#visualise the UMAP and the selected cell types
umap_annot_PANORAMA_benchmark <- ggplot(data_dimred_comb, aes(x = UMAP1, y = UMAP2, color = HarmonisedName)) +
  geom_point(size = 0.5) +
  coord_fixed(ratio = 1) +
  scale_color_manual(values=db13)+
  #facet_wrap(~ Sample, ncol = 4) +
  theme_bw() +
  ggtitle('PANORAMA_benchmark - UMAP')+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        legend.position = "bottom")+
  guides(color = guide_legend(override.aes = list(size=3),nrow=2,title=""))

#estimate relative abundance of cell types in this dataset
cellFreq <- PANORAMA_benchmark %>% 
  group_by(HarmonisedName) %>% 
  summarise(samples = n())%>% 
  mutate(freq = samples/sum(samples))

PANORAMA_benchmark_cellFreq <- ggplot(cellFreq) + aes(x=reorder(HarmonisedName,-freq),y=freq,fill=HarmonisedName)+
  geom_bar(stat='identity',width = 0.58,alpha=1,color='black',size=0.2)+
  ylab('PANORAMA_benchmark - Relative abundance')+
  theme_bw()+
  theme(axis.text.y = element_text( size = 10 ),
        axis.text.x = element_text(angle = 0, vjust = 0.5, hjust = 0.5, size = 10),
        axis.title.x = element_text( size = 12,face = 'bold' ),
        axis.title.y = element_blank(),
        strip.text = element_text(size = 12,face='bold',lineheight=1),
        legend.position = "none",aspect.ratio = 1)+
  scale_fill_manual(values = db31)+
  guides(fill = guide_legend(override.aes = list(size=3),nrow=1,title=""))+
  coord_flip()

#visualise the expression of markers in the dataset 
data <- PANORAMA_benchmark[,c('CD45','CD4','CD8','CD11c','CD16','CD34','CD19','CD27','CD25')]
data <- as.matrix(data)
rng <- colQuantiles(data, probs = c(0.01, 0.99))
data01 <- t((t(data) - rng[, 1]) / (rng[, 2] - rng[, 1]))
data01[data01 < 0] <- 0
data01[data01 > 1] <- 1
test <- data.frame(data01)
test$Class <- PANORAMA_benchmark$HarmonisedName
test <- melt(test,id.vars = c('Class'))
PANORAMA_benchmark_ridgesPlot <- ggplot(test, aes(x = value, y = Class, fill = variable)) + 
  geom_density_ridges(scale=1)+
  facet_wrap(~ variable, scales = "free", nrow = 4)+
  theme_ridges() + 
  theme(legend.position = "none",
        axis.text.x = element_text(angle=45, hjust=1,size = 7),
        axis.text.y = element_text(hjust=1,size = 7),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),aspect.ratio = 1.58)+
  scale_x_continuous(breaks=c(0,0.2,0.4,0.6,0.8,1))


 myfile <- paste('PANORAMA_benchmark_cellFreqq.pdf',sep ='')
 pdf(myfile)
 print(PANORAMA_benchmark_cellFreq)
 dev.off()

 myfile <- paste('PANORAMA_benchmark_ridgesPlot.pdf',sep ='')
 pdf(myfile)
 print(PANORAMA_benchmark_ridgesPlot)
 dev.off()



data_flow <- as.matrix(PANORAMA_benchmark[,common.col.names])
# generate flowframe from the data
ff_new <- flowFrame(exprs = data_flow, desc = list(FIL = 1))

start_time = Sys.time()
# run FlowSOM (with set.seed for reproducibility)
set.seed(123)
out_fSOM <- FlowSOM::ReadInput(ff_new, transform = F, scale = F, compensate = F)
out_fSOM <- FlowSOM::BuildSOM(out_fSOM, colsToUse = common.col.names)
out_fSOM <- FlowSOM::BuildMST(out_fSOM)
labels <- out_fSOM$map$mapping[,1]
end_time = Sys.time()
end_time-start_time

# make heatmap of 100 FlowSOM clusters
heat_mat <- matrix(NA, nrow = 100, ncol = length(common.col.names))
for(i in 1:100) {
  temp_mat <- data_flow[labels == i, common.col.names]
  if(is.vector(temp_mat)){
    heat_mat[i,] <- as.matrix(t(temp_mat))
  }else{
     heat_mat[i,] <- as.matrix(apply(temp_mat, 2, function (x) mean(x, na.rm = T)))
  }
}
rownames(heat_mat) <- paste("C",1:100,sep='')
colnames(heat_mat) <- common.col.names

#ready to prepare the testing dataset for scaffold mapping
testingSetLabels <- labels
testingSetValues <- as.matrix(PANORAMA_benchmark[,common.col.names])
testingSet.tab <- tabulateCellTypes(testingSetValues,testingSetLabels)
m <- as.matrix(testingSet.tab[, common.col.names])
rownames(m) <- testingSet.tab$cellType

############################################################################
#generate the training data that will be used as landmarks for the scaffold
############################################################################
trainingSetLabels <- healthy.AB.data.filtered$HarmonisedName
trainingSetValues <- as.matrix(healthy.AB.data.filtered[,common.col.names])
trainingSet.tab <- tabulateCellTypes(trainingSetValues,trainingSetLabels)

#perform the same computation for the training set
att <- as.matrix(trainingSet.tab[,common.col.names])
row.names(att) <- trainingSet.tab$cellType
#compute the distance from the testing set to the training set which is the landmark using cosine similarity
dd_controls_to_landmarks <- t(apply(m, 1, function(x, att) {cosine_similarity_from_matrix(x, att)}, att))
#careful here, columns of the distance matrix are the rows of the landmark matrix
colnames(dd_controls_to_landmarks) <- rownames(att)

predictionsMapping <- data.frame()
#parse the distance matrix and find the max value per row
for(idx in 1:nrow(dd_controls_to_landmarks)){
  a <- which.max(dd_controls_to_landmarks[idx,])
  predictedLabel <- colnames(dd_controls_to_landmarks)[a]
  dt <- data.frame(clusterID=rownames(dd_controls_to_landmarks)[idx],
                   predictedLabel=predictedLabel,
                   Dist=dd_controls_to_landmarks[idx,a])
  predictionsMapping <- rbind(predictionsMapping,dt)
}
predictionsMapping$clusterID <- paste("C",predictionsMapping$clusterID,sep='')

#combine the predictions and the real information about cell types
predictionsScaffold <- data.frame(clusterID=paste("C",labels,sep=''),
                                          realLabel=PANORAMA_benchmark$HarmonisedName)
predictionsScaffold$predictedLabel <- 'Unknown'
for(idx in 1:nrow(predictionsMapping)){
  m <- predictionsMapping[idx,'clusterID']
  predictionsScaffold[predictionsScaffold$clusterID==m,'predictedLabel'] <- predictionsMapping[idx,'predictedLabel']
}

performanceScaffold <- data.frame()
cellTypesToCheck <- unique(trainingSetLabels)
for(idx in cellTypesToCheck){
  
  currentType <- idx
  
  tp_row <- which(predictionsScaffold$realLabel==currentType & predictionsScaffold$predictedLabel==currentType)
  tn_row <- which(predictionsScaffold$realLabel!=currentType & predictionsScaffold$predictedLabel!=currentType)
  fp_row <- which(predictionsScaffold$realLabel!=currentType & predictionsScaffold$predictedLabel==currentType)
  fn_row <- which(predictionsScaffold$realLabel==currentType & predictionsScaffold$predictedLabel!=currentType)
  
  tp_row <- length(tp_row)
  tn_row <- length(tn_row)
  fp_row <- length(fp_row)
  fn_row <- length(fn_row)
  
  a <- data.frame(Sen=tp_row/(tp_row+fn_row),
                  Spe=tn_row/(tn_row+fp_row),
                  #PPV=tp_row/(tp_row+fp_row),
                  F1=2*tp_row/(2*tp_row+fp_row+fn_row),
                  Acc=(tp_row+tn_row)/(tp_row+tn_row+fp_row+fn_row))
  performanceScaffold <- rbind(performanceScaffold,a)
}
performanceScaffold$CellType <- cellTypesToCheck


#APPROACH #2: CytoAnno
#below we use lines of code to convert the appropriate training/testing data into the CyAnno format:
#tmp <- data.frame(testingSetValues)
#tmp$labels <- testingSetLabels
#write.csv(tmp,file='cytoAnno_PANORAMA_benchmark/live/Testing.csv',row.names=F,col.names = T)

#tmp <- data.frame(trainingSetValues)
#tmp$labels <- trainingSetLabels
#write.csv(tmp,file='cytoAnno_PANORAMA_benchmark/live/Training.csv',row.names=F,col.names = T)

#for(idx in unique(testingSetLabels)){
#  a <- which(trainingSetLabels==idx)
#  tmp1 <- data.frame(trainingSetValues[a,])
#  tmp2 <- trainingSetLabels[a]
#  str <- paste('cytoAnno_PANORAMA_benchmark/manual/',idx,'.csv',sep='')
#  write.csv(tmp1,file=str,row.names=F,col.names = T)
#}      


#once CyAnno is executed it produces the following file which contains the predicted labels
#Method_x__Posterior_Probability_Prediction_Result.csv --> renamed to PANORAMAL_benchmarkr__Posterior_Probability_Prediction_Result.csv
#we need to read this file and estimate the classification performance similar to the previous approaches

predictionsCyAnno <- read.csv('PANORAMA__Posterior_Probability_Prediction_Result.csv')
testingSetLabels <- PANORAMA_benchmark$HarmonisedName
testingSetValues <- as.matrix(PANORAMA_benchmark[,common.col.names])

#make the column names consistent same as before
colnames(predictionsCyAnno)[3] <- "realLabel"
colnames(predictionsCyAnno)[4] <- "predictedLabel"
performanceCyAnno <- data.frame()
for(idx in unique(testingSetLabels)){
  
  currentType <- idx
  
  tp_row <- which(predictionsCyAnno$realLabel==currentType & predictionsCyAnno$predictedLabel==currentType)
  tn_row <- which(predictionsCyAnno$realLabel!=currentType & predictionsCyAnno$predictedLabel!=currentType)
  fp_row <- which(predictionsCyAnno$realLabel!=currentType & predictionsCyAnno$predictedLabel==currentType)
  fn_row <- which(predictionsCyAnno$realLabel==currentType & predictionsCyAnno$predictedLabel!=currentType)
  
  tp_row <- length(tp_row)
  tn_row <- length(tn_row)
  fp_row <- length(fp_row)
  fn_row <- length(fn_row)
  
  
  a <- data.frame(Sen=tp_row/(tp_row+fn_row),
                  Spe=tn_row/(tn_row+fp_row),
                  #PPV=tp_row/(tp_row+fp_row),
                  F1=2*tp_row/(2*tp_row+fp_row+fn_row),
                  Acc=(tp_row+tn_row)/(tp_row+tn_row+fp_row+fn_row))
  performanceCyAnno <- rbind(performanceCyAnno,a)
}
performanceCyAnno$CellType <- unique(testingSetLabels)
  
#visualise the relative performance of both methods
  performanceScaffold$Method <- 'Scaffold'
  performanceCyAnno$Method <- 'CyAnno'
  
  plotData <- rbind(performanceScaffold,performanceCyAnno)
  
  tmp <- melt(plotData,id.vars = c("CellType","Method" ))
  tmp$Method <- factor(tmp$Method,levels = c('Scaffold','CyAnno'))
  colnames(tmp)[3] <- 'Metric'
  tmp$Metric <- factor(tmp$Metric,levels = c('Acc','Sen','Spe','F1'))
   comparison_PANORAMA_benchmark <- ggplot(tmp, aes(x = Method, y = value, fill = Method)) +
                   geom_jitter(height=0,width=0.1, size=0.48,alpha=0.8)+
                     geom_boxplot(width=0.28,
                           alpha=0.8,
                           size=0.3,
                          outlier.colour = "black",
                          outlier.shape = NA,
                          outlier.fill = "red",
                          outlier.size = 0.05,
                          notch = FALSE,color='black')+
                facet_wrap(~Metric,nrow = 1)+
                ggtitle('PANORAMA_benchmark performance')+
                theme_bw() +
                ylim(0,1)+
                scale_fill_manual(values = color_qual_flow2)+
                theme(axis.text.y = element_text( size = 12 ),
                      axis.text.x = element_text(angle = 90, vjust = 0.05, hjust = 0.95, size = 12),
                      axis.title.x = element_blank(),
                      axis.title.y = element_blank(),
                      strip.text = element_text(size = 12,face='bold',lineheight=1),
                      legend.position = "none",aspect.ratio = 0.8)

   
   myfile <- paste('comparison_PANORAMA_benchmark.pdf',sep ='')
   pdf(myfile)
   print(comparison_PANORAMA_benchmark)
   dev.off()

   #show all plots --> TODO save the plots and load them using knitr::include_graphics(myfile) for better quality
   print(umap_annot_AML_benchmark)
   print(AML_benchmark_cellFreq) 
   print(AML_benchmark_ridgesPlot)
   print(comparison_AML_benchmark)
   
   print(umap_annot_BMMC_benchmark)
   print(BMMC_benchmark_cellFreq) 
   print(BMMC_benchmark_ridgesPlot)
   print(comparison_BMMC_benchmark)
   
   print(umap_annot_PANORAMA_benchmark)
   print(PANORAMA_benchmark_cellFreq) 
   print(PANORAMA_benchmark_ridgesPlot)
   print(comparison_PANORAMA_benchmark)
  
```