Part2.Rmd

---
title: "Processing Benedicte's cohort - Patient data 0h (pre-treatment) and healthy controls of the study"
author: "Dimitrios Kleftogiannis"
date: "2024-02-22"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

### 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 process and phenotype patient data as well as the healthy controls from the study "Early response evaluation by single cell signaling profiling in acute myeloid leukemia", doi: 10.1038/s41467-022-35624-4
  
### 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(ncdfFlow) 
library(ggplot2)
library(readxl)
library(dplyr)
library(ggridges)
library(MASS)
library(matrixStats)
library(reshape2)
library(ggthemes)
library(RColorBrewer)
library(umap)
library(ggbeeswarm)

```

## Setting workspace and loading in-house functions

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


## Load the patient data of Benedicte's cohort
To identify the patient data files, we need to load an external file that contains the metadata and filter the ones names with 0h in the column called condition. This corresponds to pre-treatment data. 

```{r load patient dataset,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}

filename <- "AML_Benedicte_cohort/barcode_metadata_leukemia_cohort.xlsx"
#read the samples
md <- read_excel(filename)
#find the Ref or Hv 
patient_samples <- md[md$condition=='0h',]
file_names_patients <- patient_samples$file_name
file_names_patients <- paste('AML_Benedicte_cohort/raw_data/',file_names_patients,sep='')

file_barcode_patients <- patient_samples$barcode

#read them
fs_patients  <- read.ncdfFlowSet(files = file_names_patients,
                              transformation = F,
                              phenoData = ,
                              truncate_max_range = F)

#generate the expression matrix --> here the names are the metal codes and we need to change that
data_fs_patient <- fsApply(fs_patients, exprs)
#fetch the panel info
panel_info <- pData(parameters(fs_patients[[1]]))
#rename the matrix
colnames(data_fs_patient) <- as.character(panel_info$desc)
#add the patient label name 
patient_labels <- data.frame()
for(idx in 1:length(file_names_patients)){
  a <- dim(fs_patients[[idx]])[1]
  a <-rep(paste('P',idx,sep=''),a)
  a <- data.frame(Sample=a)
  patient_labels <- rbind(patient_labels,a)
}

#retrieve the markers of interest - including the Iridium
col.names <- c("CD11b","pAxl","CD8a","CD33","CyclinB1","pNFkB",
               "pErk","CD34","CD3","CD123","pSTAT1","CD56",
               "pP38","pSTAT3","CD14","pCREB",
               "pHistone3","Caspase3_cleaved","CD117","CD38","CD4","CD64","CD16",
               "pSTAT5","CD66b","CD20","p4E_BP1","pAkt","pRB",
               "CD45","pS6","HLA-DR","CD7","191Ir_DNA1","193Ir_DNA2")

#fetch the channels of interest
data_fs_patient <- data_fs_patient[,col.names]
#rename long channel names and possibly problematic
col.names[17] <- 'pHist3'
col.names[18] <- 'Casp3'
col.names[27] <- 'p4EBP1'
col.names[32] <- 'HLA.DR'
col.names[34] <- 'X191Ir'
col.names[35] <- 'X193Ir'

colnames(data_fs_patient)[17] <- 'pHist3'
colnames(data_fs_patient)[18] <- 'Casp3'
colnames(data_fs_patient)[27] <- 'p4EBP1'
colnames(data_fs_patient)[32] <- 'HLA.DR'
colnames(data_fs_patient)[34] <- 'X191Ir'
colnames(data_fs_patient)[35] <- 'X193Ir'

```

## Load the heathly control data of Benedicte's cohort
The published study contains 7 healthy individuals that were analysed together with the patient data for comparisons. 
To identify the healthy control data files, we need to load an external file that contains the metadata and filter the ones names with Hv in the column called condition. This corresponds to healthy individuals included in the publication. 

```{r load healthy cohort dataset,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}

filename <- "AML_Benedicte_cohort/barcode_metadata_leukemia_cohort.xlsx"
#read the samples
md <- read_excel(filename)
#find the Ref or Hv 
control_samples <- md[md$condition=='Hv',]
file_names_ctrl <- control_samples$file_name
file_names_ctrl <- paste('AML_Benedicte_cohort/raw_data/',file_names_ctrl,sep='')


file_barcode_ctrl <- control_samples$barcode

#read them
fs_ctrl  <- read.ncdfFlowSet(files = file_names_ctrl,
                              transformation = F,
                              phenoData = ,
                              truncate_max_range = F)

#generate the expression matrix --> here the names are the metal codes and we need to change that
data_fs_ctrl <- fsApply(fs_ctrl, exprs)
#fetch the panel info
panel_info <- pData(parameters(fs_ctrl[[1]]))
#rename the matrix
colnames(data_fs_ctrl) <- as.character(panel_info$desc)
#add the control label name 
ctrl_labels <- data.frame()
for(idx in 1:length(file_names_ctrl)){
  a <- dim(fs_ctrl[[idx]])[1]
  a <-rep(paste('S',idx,sep=''),a)
  a <- data.frame(Sample=a)
  ctrl_labels <- rbind(ctrl_labels,a)
}

#retrieve the markers of interest - including the Iridium
col.names <- c("CD11b","pAxl","CD8a","CD33","CyclinB1","pNFkB",
               "pErk","CD34","CD3","CD123","pSTAT1","CD56",
               "pP38","pSTAT3","CD14","pCREB",
               "pHistone3","Caspase3_cleaved","CD117","CD38","CD4","CD64","CD16",
               "pSTAT5","CD66b","CD20","p4E_BP1","pAkt","pRB",
               "CD45","pS6","HLA-DR","CD7","191Ir_DNA1","193Ir_DNA2")

#fetch the channels of interest
data_fs_ctrl <- data_fs_ctrl[,col.names]

#rename long channel names and possibly problematic
col.names[17] <- 'pHist3'
col.names[18] <- 'Casp3'
col.names[27] <- 'p4EBP1'
col.names[32] <- 'HLA.DR'
col.names[34] <- 'X191Ir'
col.names[35] <- 'X193Ir'

colnames(data_fs_ctrl)[17] <- 'pHist3'
colnames(data_fs_ctrl)[18] <- 'Casp3'
colnames(data_fs_ctrl)[27] <- 'p4EBP1'
colnames(data_fs_ctrl)[32] <- 'HLA.DR'
colnames(data_fs_ctrl)[34] <- 'X191Ir'
colnames(data_fs_ctrl)[35] <- 'X193Ir'

```

## Load the reference data to build a scaffold and annotate the cells from Benedicte's cohort
To build a scaffold or equally denoted a roadmap 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.

It is the same dataset that was used in the self consistency test, and in the comparison analysis with other methods. 

```{r load scaffold data,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}

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

# 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')


#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)

#need to use the antibody names as column names in the dataset. 
myCellCutoff <- 299
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)
#we rename CD8 to become compatible with our AML data
colnames(healthy.AB.data.filtered)[82] <- 'CD8a'
#find the common markers between the scaffold dataset and the AML cohort
common.col.names <- intersect(col.names,colnames(healthy.AB.data.filtered)[1:(ncol(healthy.AB.data.filtered)-1)])

#since the scaffold dataset is of high resolution, that cannot be achieved using a conventional cytof panel we will merge cell types of the same lineage and we will generate fewer cell populations that reflect the major ones.

selected_cellTypes_Healthy <- c("CD4+ cytotoxic T cells", "CD4+ memory T cells","CD4+ naive T cells",
                                "CD56brightCD16- NK cells","CD56dimCD16+ NK cells",
                                "CD8+ central memory T cells","CD8+ effector memory T cells","CD8+ naive T cells","CD8+CD103+ tissue resident memory T cells",
                                "Class switched memory B cells","Immature B cells","Mature naive B cells","Nonswitched memory B cells","Pre-B cells","Pro-B cells",
                                "Classical Monocytes","Non-classical monocytes",
                                "HSCs & MPPs",
                                "Plasmacytoid dendritic cell progenitors","Plasmacytoid dendritic cells")  


idx <- which(healthy.AB.data.filtered$cellType %in% selected_cellTypes_Healthy)
healthy.AB.data.filtered <- healthy.AB.data.filtered[idx,]
umap.filtered <- umap.filtered[idx,]


#rename the cell types
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'


healthy.AB.data.filtered[healthy.AB.data.filtered$cellType=="CD56dimCD16+ NK cells"  | 
                            healthy.AB.data.filtered$cellType=="CD56brightCD16- NK cells",'HarmonisedName'] <- 'NK'

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'

healthy.AB.data.filtered[healthy.AB.data.filtered$cellType=="Class switched memory B cells" |
                         healthy.AB.data.filtered$cellType=="Immature B cells"  |
                         healthy.AB.data.filtered$cellType=="Mature naive B cells" |
                         healthy.AB.data.filtered$cellType=="Nonswitched memory B cells" |
                         healthy.AB.data.filtered$cellType=="Pre-B cells" |
                         healthy.AB.data.filtered$cellType=="Pro-B cells"  ,'HarmonisedName'] <- 'B'

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=="Plasmacytoid dendritic cell progenitors"|
                            healthy.AB.data.filtered$cellType=="Plasmacytoid dendritic cells"  ,'HarmonisedName'] <- 'pDCs'


tmp <- cbind(umap.filtered,healthy.AB.data.filtered$HarmonisedName)
colnames(tmp)[3] <- "cellType"
umap_scaffold_original <- ggplot(tmp, aes(x = UMAP1, y = UMAP2,color=cellType)) +
    geom_point(size = .1) +
    coord_fixed(ratio = 1) +
    ggtitle("Roadmap of human hematopoiesis - Original UMAP")+
    scale_color_manual(values=db31)+
  theme_few()+
    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(),aspect.ratio = 1)+
  guides(color = guide_legend(override.aes = list(size=1.4),ncol=4,title=""))

#compute statistics about the scaffold data
cellFreq <- healthy.AB.data.filtered %>% 
  group_by(HarmonisedName) %>% 
  summarise(samples = n())%>% 
  mutate(freq = samples/sum(samples))

 scaffold_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('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()

 #visualise the expression of markers in the scaffold dataset 
data <- healthy.AB.data.filtered[,c('CD45','CD4','CD8a','CD20','CD16','CD34','CD11b','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 <- healthy.AB.data.filtered$HarmonisedName
test <- melt(test,id.vars = c('Class'))
scaffold_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))


# plot UMAP with expression overlayed -> use the original coordinates of the cohort
data_dimred_exp <- cbind(healthy.AB.data.filtered[,common.col.names], umap.filtered)
data_dimred_melt <- melt(data_dimred_exp, id.vars = c("UMAP1", "UMAP2"))

#overlay colors
umap_scaffold_original_overlay <- ggplot(data_dimred_melt, aes(x = UMAP1, y = UMAP2, color = value)) +
  geom_point(size = 0.05) +
  coord_fixed(ratio = 1) +
  scale_colour_gradientn(colours = color_grad_flow2, limits = c(0,1.5)) +
  facet_wrap(~ variable, ncol = 6) +
  theme_few() +
  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(),strip.text = element_text(size = 6,face='bold',lineheight=1))

library(umap)
#generate new umap using the filtered data
data_dimred <- healthy.AB.data.filtered[ , common.col.names]
# UMAP all together the scaffold data
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")
# plot UMAP with expression overlayed
data_dimred_exp <- cbind(healthy.AB.data.filtered[,common.col.names], umap)
colnames(data_dimred_exp)[ncol(data_dimred_exp)-1] <- 'UMAP1'
colnames(data_dimred_exp)[ncol(data_dimred_exp)] <- 'UMAP2'
data_dimred_melt <- melt(data_dimred_exp, id.vars = c("UMAP1", "UMAP2"))

save(data_umap,file='data_umap.RDa')

#overlay colors
umap_scaffold_new <- ggplot(data_dimred_melt, aes(x = UMAP1, y = UMAP2, color = value)) +
  geom_point(size = 0.05) +
  coord_fixed(ratio = 1) +
  scale_colour_gradientn(colours = color_grad_flow2, limits = c(0,1.5)) +
  facet_wrap(~ variable, ncol = 6) +
  theme_few() +
  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(),strip.text = element_text(size = 6,face='bold',lineheight=1))

#visualise again UMAP with the annotated cell types
tmp <- cbind(umap,healthy.AB.data.filtered$HarmonisedName)
colnames(tmp)[3] <- "cellType"
umap_scaffold_new_annot <- ggplot(tmp, aes(x = UMAP1, y = UMAP2,color=cellType)) +
    geom_point(size = .1) +
    coord_fixed(ratio = 1) +
    ggtitle("Roadmap of human hematopoiesis - new UMAP")+
    scale_color_manual(values=db31)+
    theme_few()+
    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(),
          aspect.ratio = 1)+
  guides(color = guide_legend(override.aes = list(size=1.4),ncol=4,title=""))

```

## Visualise the scaffold data
This is the reference we will use to annotate the AML cohort from doi: 10.1038/s41467-022-35624-4
```{r show figures,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}
print(umap_scaffold_original)
print(scaffold_cellFreq)
print(scaffold_ridgesPlot)
print(umap_scaffold_original_overlay)
print(umap_scaffold_new)
print(umap_scaffold_new_annot)

   myfile <- paste('Supplementary_Figure_5a.pdf',sep ='')
   pdf(myfile)
   print(umap_scaffold_new_annot)
   dev.off()
  
   myfile <- paste('Supplementary_Figure_6a.pdf',sep ='')
   pdf(myfile)
   print(umap_scaffold_new)
   dev.off()
   
   myfile <- paste('Supplementary_Figure_6b.pdf',sep ='')
   pdf(myfile)
   print(scaffold_ridgesPlot)
   dev.off()
   
```

## Transform the healthy data using barcode specific and channel specific cofactors.
Harmonize the values in the control samples of the cohort.

```{r scale healthy data with cofactors,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}

  #first transform the data with global cofactor of 5 
  asinh_scale <- 5
  data_fs_ctrl <- asinh(data_fs_ctrl/ asinh_scale)
  data_fs_ctrl <- as.data.frame(data_fs_ctrl)
  data_fs_ctrl$Sample <- ctrl_labels
  
  #load the cofactor data that derive from separate analysis of control samples
  load('cofactors_v2.RDa')
  summariseData <- data.frame()
  for(myCtrl in 1:length(file_names_ctrl)){
    promptKey <- paste('S',myCtrl,sep='')
    idx <- which(ctrl_labels==promptKey)
    myCtrlData <- data_fs_ctrl[idx,]
    barcodeKey <- file_barcode_ctrl[myCtrl]
    #scale all markers with their corresponding factor
    for(myMarker in col.names){
        currentCof <- coFactorMatrix[barcodeKey,myMarker]
        myCtrlData[,myMarker] <- asinh(5*sinh(myCtrlData[,myMarker])/currentCof)
    }
    summariseData <- rbind(summariseData,myCtrlData)
    gc()
  }
  # normalize based on the 99.99th percentile to have everything between 0 and 1
  data_fs_ctrl <- summariseData[,1:(ncol(summariseData)-1)]
  rm(summariseData)
  per.tmp <- apply(data_fs_ctrl, 2, function(x) quantile(x, 0.9999, names = F))
  data_fs_ctrl <- t(t(data_fs_ctrl) / as.numeric(per.tmp))
  #clean memory
  gc()
  
#gate using CD45 and CD66b to identify cutoffs for filtering
  tmp <- data.frame(CD45=data_fs_ctrl[,'CD45'],
                    CD66b=data_fs_ctrl[,'CD66b'])
  
 o <-  ggplot(tmp, aes(x = CD66b, y = CD45)) +
      geom_hex(bins = 100) +
      ylab('CD45')+
      xlab('CD66b')+
      ggtitle('Healthy controls - Inspecting granulocytes')+
      geom_hline(yintercept = 0.25,color='red')+
       geom_vline(xintercept = 0.12,color='red')+
      scale_x_continuous(limits = c(-0.1, 1.5))+
      scale_y_continuous(limits = c(-0.1, 1.5))+
      coord_fixed(ratio = 1) +
      scale_fill_gradientn(colours = db31, trans = "sqrt") +
      theme_bw()+ theme(legend.position="none", axis.text = element_text(size = 10), axis.title = element_text(size = 10),aspect.ratio = 1)
  
 print(o)
 #filter out cells with higher expression of CD66b that are negative for CD45 --> these are the granulocytes and do not match with the scaffold data
 idx <- which(data_fs_ctrl[,'CD66b']>=0.12 & data_fs_ctrl[,'CD45']<0.2508)
 data_fs_ctrl <- data_fs_ctrl[-idx,]
 ctrl_labels <- ctrl_labels[-idx,]
 
```

## Visualise and cluster the control data of the cohort
```{r visualise and cluster the control data,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}

data <- data_fs_ctrl
data <- as.data.frame(data)
#add cell ID
data$Sample <- ctrl_labels
data$cell_id <- 1:nrow(data)
#subsample for umap
# subsample cohort for dimensionality reduction
n_samples <- unique(data[,"Sample"])
#save the subsampled data
data_sub <- data.frame()
#original n_sub_i <- 1500
n_sub_i <- 3500
for(i in n_samples){
  data_i <- data[data[,"Sample"]==i,]
  n_i <- nrow(data_i)
  set.seed(123)
  if(n_sub_i > n_i){
    ix_i <- sample(1:n_i, n_i)
  }
  else{
    ix_i <- sample(1:n_i, n_sub_i)
  }
  data_sub_i <- data_i[ix_i,]
  data_sub <- rbind(data_sub, data_sub_i)
}


clustering_cols <- common.col.names
data_dimred <- data_sub[ , clustering_cols]

## 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")

## plot UMAP black
#umap_ctrl_black <- ggplot(umap, aes(x = UMAP1, y = UMAP2)) +
#  geom_point(size = 0.5) +
#  coord_fixed(ratio = 1)
#umap_black

## 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("Sample","cell_id")], data_dimred_exp)

##overlay colors
#umap_ctrl_overlay <- ggplot(data_dimred_melt, aes(x = UMAP1, y = UMAP2, color = value)) +
#  geom_point(size = 0.05) +
#  coord_fixed(ratio = 1) +
#  scale_colour_gradientn(colours = color_grad_flow2, limits = c(0,1.5)) +
#  facet_wrap(~ variable, ncol = 6) +
#  theme_few() +
#  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(),strip.text = element_text(size = 6,face='bold',lineheight=1))
#umap_overlay_ctrl

library(umap)
load('data_umap.RDa')
#project into the scaffold UMAP --> we use data_umap which is the trained umap for the scaffold
start_time = Sys.time()
ctrl.prjct.umap <- predict(data_umap, data_dimred)
end_time = Sys.time()
end_time-start_time

dim(data_dimred)
data_dimred_exp <- cbind(data_dimred, ctrl.prjct.umap)
colnames(data_dimred_exp)[ncol(data_dimred_exp)] <- 'UMAP2'
colnames(data_dimred_exp)[ncol(data_dimred_exp)-1] <- 'UMAP1'
data_dimred_melt <- melt(data_dimred_exp, id.vars = c("UMAP1", "UMAP2"))

data_dimred_comb  <- cbind(data_sub[,c("Sample","cell_id")], data_dimred_exp)


#overlay colors
umap_projected_ctrl_overlay <- ggplot(data_dimred_melt, aes(x = UMAP1, y = UMAP2, color = value)) +
  geom_point(size = 0.05) +
  coord_fixed(ratio = 1) +
  scale_colour_gradientn(colours = color_grad_flow2, limits = c(0,1.5)) +
  facet_wrap(~ variable, ncol = 6) +
  ggtitle('Healthy control data')+
  theme_few() +
  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(),strip.text = element_text(size = 6,face='bold',lineheight=1))

#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(data_fs_ctrl[,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(data_fs_ctrl[,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='')

#compute the quantiles and give "unknown" predictions below 75th
myQT <- quantile(dd_controls_to_landmarks)

predictionsMapping[predictionsMapping$Dist < myQT[4],'predictedLabel'] <- 'Unknown'

#combine the cluster predictions and the estimated labels of cell types
labels <- paste("C",labels,sep='')
for(idx in 1:nrow(predictionsMapping)){
  m <- predictionsMapping[idx,'clusterID']
  labels[labels==m] <- predictionsMapping[idx,'predictedLabel']
}

# plot resulting clusters on umap
data_clust <- cbind(data, labels)
data_umap_clust <- merge(data_dimred_comb, data_clust[,c("labels", "cell_id")], by = "cell_id")
data_umap_clust$labels <- as.factor(data_umap_clust$labels)

umap_projected_ctrl_annot <- ggplot(data_umap_clust, aes(x = UMAP1, y = UMAP2, color = labels)) +
  geom_point(size = .1) +
  coord_fixed(ratio = 1) +
  scale_color_manual(values = db31) +
   theme_few() +
  ggtitle('Healthy control data')+
  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(),aspect.ratio = 1)+
  guides(colour = guide_legend(override.aes = list(size=6), title="Cell type"))

#save the healthy controls with annotated labels
#save(data_clust,file='ctrl_annotated.RDa')

#visualise the expression of markers in the controls dataset 
data <- data_clust[,c('CD45','CD4','CD8a','CD20','CD16','CD34','CD56','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 <- data_clust$labels
test <- melt(test,id.vars = c('Class'))
ctrl_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))


print(umap_projected_ctrl_overlay)
print(umap_projected_ctrl_annot)
print(ctrl_ridgesPlot)


   myfile <- paste('Supplementary_Figure_5b.pdf',sep ='')
   pdf(myfile)
   print(umap_projected_ctrl_annot)
   dev.off()
  
   myfile <- paste('Supplementary_Figure_6c.pdf',sep ='')
   pdf(myfile)
   print(umap_projected_ctrl_overlay)
   dev.off()
   
   myfile <- paste('Supplementary_Figure_6d.pdf',sep ='')
   pdf(myfile)
   print(ctrl_ridgesPlot)
   dev.off()

```



## Transform the patient data using barcode specific and channel specific cofactors.
We use the scaling factors found from healthy controls analysis to harmonise the values in patient samples.

```{r scale patient data with cofactors,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}

  #first transform the data with global cofactor of 5 
  asinh_scale <- 5
  data_fs_patient <- asinh(data_fs_patient/ asinh_scale)
  data_fs_patient <- as.data.frame(data_fs_patient)
  data_fs_patient$Sample <- patient_labels
  
  #load the cofactor data that derive from separate analysis of samples
  load('cofactors_v2.RDa')
  summariseData <- data.frame()
  for(myPatient in 1:length(file_names_patients)){
    promptKey <- paste('P',myPatient,sep='')
    idx <- which(patient_labels==promptKey)
    myPatientData <- data_fs_patient[idx,]
    barcodeKey <- file_barcode_patients[myPatient]
    #scale all markers with their corresponding factor
    for(myMarker in col.names){
        currentCof <- coFactorMatrix[barcodeKey,myMarker]
        myPatientData[,myMarker] <- asinh(5*sinh(myPatientData[,myMarker])/currentCof)
    }
    summariseData <- rbind(summariseData,myPatientData)
    gc()
  }
  # normalize based on the 99.99th percentile to have everything between 0 and 1
  data_fs_patient <- summariseData[,1:(ncol(summariseData)-1)]
  rm(summariseData)
  per.tmp <- apply(data_fs_patient, 2, function(x) quantile(x, 0.9999, names = F))
  data_fs_patient <- t(t(data_fs_patient) / as.numeric(per.tmp))
  #clean memory
  gc()
  
  #filter out the granulocytes --> cells with high CD66b and low CD45 are excluded from the analysis
  tmp <- data.frame(CD45=data_fs_patient[,'CD45'],
                    CD66b=data_fs_patient[,'CD66b'])
  
 o <-  ggplot(tmp, aes(x = CD66b, y = CD45)) +
      geom_hex(bins = 100) +
      ylab('CD45')+
      xlab('CD66b')+
    ggtitle('Patient data - Inspecting granulocytes')+
      geom_hline(yintercept = 0.2,color='red')+
       geom_vline(xintercept = 0.08,color='red')+
      scale_x_continuous(limits = c(-0.1, 1.5))+
      scale_y_continuous(limits = c(-0.1, 1.5))+
      coord_fixed(ratio = 1) +
      scale_fill_gradientn(colours = db31, trans = "sqrt") +
      theme_bw()+ theme(legend.position="none", axis.text = element_text(size = 10), axis.title = element_text(size = 10),aspect.ratio = 1)
 
 print(o) 
 
 #filter out cells with higher expression of CD66b that are negative for CD45 --> these are the granulocytes
 idx <- which(data_fs_patient[,'CD66b']>=0.08 & data_fs_patient[,'CD45']<0.2)
 data_fs_patient <- data_fs_patient[-idx,]
 patient_labels <- patient_labels[-idx,]
 
```


## Visualise and cluster the patient data of the cohort
```{r visualise and cluster the patient data,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}

data <- data_fs_patient
data <- as.data.frame(data)
#add cell ID
data$Sample <- patient_labels
data$cell_id <- 1:nrow(data)
#subsample for umap
# subsample cohort for dimensionality reduction
n_samples <- unique(data[,"Sample"])
#save the subsampled data
data_sub <- data.frame()
#original n_sub_i <- 1500
n_sub_i <- 1500
for(i in n_samples){
  data_i <- data[data[,"Sample"]==i,]
  n_i <- nrow(data_i)
  set.seed(123)
  if(n_sub_i > n_i){
    ix_i <- sample(1:n_i, n_i)
  }
  else{
    ix_i <- sample(1:n_i, n_sub_i)
  }
  data_sub_i <- data_i[ix_i,]
  data_sub <- rbind(data_sub, data_sub_i)
}

clustering_cols <- common.col.names
data_dimred <- data_sub[ , clustering_cols]

library(umap)
load('data_umap.RDa')

#project into the scaffold UMAP --> we use data_umap which is the trained umap for the scaffold
start_time = Sys.time()
patient.prjct.umap <- predict(data_umap, data_dimred)
end_time = Sys.time()
end_time-start_time

dim(data_dimred)
data_dimred_exp <- cbind(data_dimred, patient.prjct.umap)
colnames(data_dimred_exp)[ncol(data_dimred_exp)] <- 'UMAP2'
colnames(data_dimred_exp)[ncol(data_dimred_exp)-1] <- 'UMAP1'
data_dimred_melt <- melt(data_dimred_exp, id.vars = c("UMAP1", "UMAP2"))

data_dimred_comb  <- cbind(data_sub[,c("Sample","cell_id")], data_dimred_exp)

#overlay colors
umap_projected_patient_overlay <- ggplot(data_dimred_melt, aes(x = UMAP1, y = UMAP2, color = value)) +
  geom_point(size = 0.05) +
  coord_fixed(ratio = 1) +
  scale_colour_gradientn(colours = color_grad_flow2, limits = c(0,1.5)) +
  facet_wrap(~ variable, ncol = 6) +
  ggtitle('Patient data')+
  theme_few() +
  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(),strip.text = element_text(size = 6,face='bold',lineheight=1))

#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(data_fs_patient[,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(data_fs_patient[,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='')

#compute the quantiles and give "unknown" predictions below 75th
myQT <- quantile(dd_controls_to_landmarks)

predictionsMapping[predictionsMapping$Dist < myQT[4],'predictedLabel'] <- 'Unknown'

#combine the cluster predictions and the estimated labels of cell types
labels <- paste("C",labels,sep='')
for(idx in 1:nrow(predictionsMapping)){
  m <- predictionsMapping[idx,'clusterID']
  labels[labels==m] <- predictionsMapping[idx,'predictedLabel']
}

# plot resulting clusters on umap
data_clust <- cbind(data, labels)
data_umap_clust <- merge(data_dimred_comb, data_clust[,c("labels", "cell_id")], by = "cell_id")
data_umap_clust$labels <- as.factor(data_umap_clust$labels)

data_umap_clust <- data_umap_clust[data_umap_clust$labels!='Unknown',]
data_umap_clust$labels <- as.factor(data_umap_clust$labels)

umap_projected_patient_annot <- ggplot(data_umap_clust, aes(x = UMAP1, y = UMAP2, color = labels)) +
  geom_point(size = .1) +
  coord_fixed(ratio = 1) +
  scale_color_manual(values = db31) +
   theme_few() +
  ggtitle('Patient data')+
  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(),aspect.ratio = 1)+
  guides(colour = guide_legend(override.aes = list(size=6), title="Cell type"))

#filter out the "Unknown" cells to improve data quality
idx <- which(labels=='Unknown')
data_clust <- data_clust[-idx,]

#save the healthy controls with annotated labels
#save(data_clust,file='patient_annotated.RDa')

#visualise the expression of markers in the controls dataset 
data <- data_clust[,c('CD45','CD4','CD8a','CD20','CD16','CD34','CD56','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 <- data_clust$labels
test <- melt(test,id.vars = c('Class'))
patient_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))


print(umap_projected_patient_overlay)
print(umap_projected_patient_annot)
print(patient_ridgesPlot)

   myfile <- paste('Supplementary_Figure_5c.pdf',sep ='')
   pdf(myfile)
   print(umap_projected_patient_annot)
   dev.off()
  
   myfile <- paste('Supplementary_Figure_6e.pdf',sep ='')
   pdf(myfile)
   print(umap_projected_patient_overlay)
   dev.off()
   
   myfile <- paste('Supplementary_Figure_6f.pdf',sep ='')
   pdf(myfile)
   print(patient_ridgesPlot)
   dev.off()

```

## Summarise the cell types and perform comparison between AML vs. Control data
In this subsection we perform exploratory analysis considering the relative cellular abundances.
First we compare Leukemia patients versus healthy controls.
Then we perform stratification of patients given their survival time and we explore their predictive value
```{r compare cellular abundances,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}
library(ggbeeswarm)

#erase everything before loading again
rm(data_clust)
load('ctrl_annotated.RDa')
#healthy control data annotated
ctrl_data <- data_clust

rm(data_clust)
#patient data annotated before treatment
load('patient_annotated.RDa')
patient_data <- data_clust


cellFreq_ctrl <- ctrl_data %>% 
  group_by(Sample,labels) %>% 
  summarise(samples = n())%>% 
  mutate(freq = samples/sum(samples))

cellFreq_ctrl$Condition <- 'Controls'

cellFreq_patient <- patient_data %>% 
  group_by(Sample,labels) %>% 
  summarise(samples = n())%>% 
  mutate(freq = samples/sum(samples))

cellFreq_patient$Condition <- 'Leukemia'

plot_data <- rbind(cellFreq_patient,cellFreq_ctrl)
plot_data <- plot_data[,c('Condition','labels','freq')]
pt <- melt(plot_data,id.vars = c('Condition','labels'))

pt$Condition <- factor(pt$Condition,levels = c('Leukemia','Controls'))

cellularAbundance_patient_ctrl <- ggplot(pt,aes(x=Condition,y=(value),fill=Condition),color=Condition)+
  geom_beeswarm(size=0.78,alpha=0.8)+
  geom_boxplot(width=0.38 ,size=0.58,alpha=0.8,fatten=TRUE,
               position=position_dodge(width=0.78),
               outlier.colour = "black",
               outlier.shape = 21,
               outlier.fill = "gray",
               outlier.size = 0.4)+
  facet_wrap(~labels,scales = "free_x",ncol = 2)+
  theme_bw()+
  theme(axis.text.y = element_text( size = 12 ),
        axis.text.x = element_text(angle = 0, vjust = 1, hjust = 1, size = 12),
        axis.title.y = element_blank(),
        axis.title.x = element_text( size = 12 ),
        strip.text = element_text(size = 14,face='bold',lineheight=1),
        legend.position = "none",
        legend.text = element_text(size=12),
        strip.background = element_rect(colour = "black", fill = "white"),
        plot.title = element_text(size = 14, face = "bold",hjust = 0.5),aspect.ratio = 0.6)+
  ylab('Cells % in population')+
  scale_fill_manual(values=c("#D73027","#74ADD1"),name='')+
  theme(panel.spacing = unit(1.8, "lines"))+coord_flip()

#estimate p values
cellTypes <- c('B','CD4_T','CD8_T','HSCs_MPPs', 'Monocytes','NK','pDCs') 
for(idx in cellTypes){
  
  str <- paste('Test for ',idx,sep='')
  cat(str)
  tmp <- plot_data[plot_data$labels==idx ,]
  res <- wilcox.test(freq ~ Condition, data = tmp, exact = FALSE)
  print(res)
}

#save the abundance plot
myfile <- paste('cellularAbundance_patient_ctrl.pdf',sep ='')
pdf(myfile)
print(cellularAbundance_patient_ctrl)
dev.off()

####################################################################################################
######## load clinical data to genereate plots that link cellular abundance ratio with clinical info 
#the clinical info file cannot be shared in its current form; please contact the authors

#load the clinical info
metadata_filename <- "AML_Benedicte_cohort/leukemia_cohort_info.xlsx"
md_0h <- read_excel(metadata_filename)
summary_patient_data <- data.frame()
for(idx in 1:nrow(patient_samples)){
  
  myPatientID <- patient_samples[idx,'patient_id']
  a <- which(md_0h$Patient_nr_cytobank==myPatientID$patient_id)
  if(length(a)>0){
    survivalTime <- md_0h[a,"5-year survival (days)"]
    chromAbber <- md_0h[a,'Karyotype']
    chromAbber <- chromAbber$Karyotype
    if(chromAbber=='NA'){
      chromAbber <- NA
    }
    mutFLT3 <- md_0h[a,'FLT3-ITD']
    mutFLT3 <- mutFLT3$`FLT3-ITD`
    if(mutFLT3=='Present'){
      mutFLT3 <- 'Mutated'
    }else{
      mutFLT3 <- 'Wt'
    }
    #fetch the patient cells
    promptKey <- paste('P',idx,sep='')
    data <- patient_data[patient_data$Sample==promptKey,]
    #compute the relative cellular abundance
    a1 <- which(data$labels=='B')
    a2 <- which(data$labels=='CD4_T')
    a3 <- which(data$labels=='CD8_T')
    a4 <- which(data$labels=='HSCs_MPPs')
    a5 <- which(data$labels=='Monocytes')
    a6 <- which(data$labels=='NK')
    a7 <- which(data$labels=='pDCs')
  
    if(survivalTime$`5-year survival (days)`=='Alive'){
      survivalTime <- 5*365
    }else{
      survivalTime <- survivalTime$`5-year survival (days)`
    }
  
    a <- data.frame(PatientNr=promptKey,
                    PatientID=myPatientID$patient_id,
                    Ncells=nrow(data),
                    B=length(a1)/nrow(data),
                    CD4_T=length(a2)/nrow(data),
                    CD8_T=length(a3)/nrow(data),
                    HSCs_MPPs=length(a4)/nrow(data),
                    Monocytes=length(a5)/nrow(data),
                    NK=length(a6)/nrow(data),
                    pDCs=length(a7)/nrow(data),
                    SurvivalTime=survivalTime,
                    karyotype=chromAbber,
                    FLT3=mutFLT3)
  summary_patient_data <- rbind(summary_patient_data,a)
  }
}

summary_patient_data$SurvivalTime <- as.numeric(summary_patient_data$SurvivalTime)
summary_patient_data$SurvivalStatus <- ifelse(summary_patient_data$SurvivalTime<5*365,1,0)

#first stratify patients based on survival time and compare their relative cellular abundances
summary_patient_data$SurvivalGroup <- ifelse(summary_patient_data$SurvivalTime<5*365,'STS','LTS')
#plot 1 based on survival
plot_data <- summary_patient_data[,c(4:10,ncol(summary_patient_data))]
plot_data <- melt(plot_data,id.vars = 'SurvivalGroup')

plot_data$SurvivalGroup <- factor(plot_data$SurvivalGroup,levels = c('STS','LTS'))

cellularAbundance_STS_LTS <- ggplot(plot_data,aes(x=SurvivalGroup,y=(value),fill=SurvivalGroup),color=SurvivalGroup)+
  geom_beeswarm(size=0.78,alpha=0.8)+
  geom_boxplot(width=0.38 ,size=0.58,alpha=0.8,fatten=TRUE,
               position=position_dodge(width=0.78),
               outlier.colour = "black",
               outlier.shape = 21,
               outlier.fill = "gray",
               outlier.size = 0.4)+
  facet_wrap(~variable,scales = "free_x",ncol = 2)+
  theme_bw()+
  theme(axis.text.y = element_text( size = 12 ),
        axis.text.x = element_text(angle = 0, vjust = 1, hjust = 1, size = 12),
        axis.title.y = element_blank(),
        axis.title.x = element_text( size = 12 ),
        strip.text = element_text(size = 14,face='bold',lineheight=1),
        legend.position = "none",
        legend.text = element_text(size=12),
        strip.background = element_rect(colour = "black", fill = "white"),
        plot.title = element_text(size = 14, face = "bold",hjust = 0.5),aspect.ratio = 0.6)+
  ylab('Cells % in population')+
  scale_fill_manual(values=c('#e6194b', '#3cb44b'),name='')+
  theme(panel.spacing = unit(1.8, "lines"))+coord_flip()

#save the abundance plot
myfile <- paste('cellularAbundance_STS_LTS.pdf',sep ='')
pdf(myfile)
print(cellularAbundance_STS_LTS)
dev.off()

#estimate p values
cellTypes <- c('B','CD4_T','CD8_T','HSCs_MPPs', 'Monocytes','NK','pDCs') 
for(idx in cellTypes){
  
  str <- paste('Test for ',idx,sep='')
  cat(str)
  tmp <- plot_data[plot_data$variable==idx ,]
  res <- wilcox.test(value ~ SurvivalGroup, data = tmp, exact = FALSE)
  print(res)
}


#plot 2 based on karyotype
plot_data <- summary_patient_data[,c(4:10,12)]
plot_data <- melt(plot_data,id.vars = 'karyotype')

plot_data$karyotype <- factor(plot_data$karyotype,levels = c('Normal','Complex','Chromosomal alter.'))

plot_data <- na.omit(plot_data)

cellularAbundance_Karyotype <- ggplot(plot_data,aes(x=karyotype,y=(value),fill=karyotype),color=karyotype)+
  geom_beeswarm(size=0.78,alpha=0.8)+
  geom_boxplot(width=0.38 ,size=0.58,alpha=0.8,fatten=TRUE,
               position=position_dodge(width=0.78),
               outlier.colour = "black",
               outlier.shape = 21,
               outlier.fill = "gray",
               outlier.size = 0.4)+
  facet_wrap(~variable,scales = "free_x",ncol = 2)+
  theme_bw()+
  theme(axis.text.y = element_text( size = 12 ),
        axis.text.x = element_text(angle = 0, vjust = 1, hjust = 1, size = 12),
        axis.title.y = element_blank(),
        axis.title.x = element_text( size = 12 ),
        strip.text = element_text(size = 14,face='bold',lineheight=1),
        legend.position = "none",
        legend.text = element_text(size=12),
        strip.background = element_rect(colour = "black", fill = "white"),
        plot.title = element_text(size = 14, face = "bold",hjust = 0.5),aspect.ratio = 0.6)+
  ylab('Cells % in population')+
  scale_fill_manual(values=c('coral1','bisque3','gray88'),name='')+
  theme(panel.spacing = unit(1.8, "lines"))+coord_flip()

#save the abundance plot
myfile <- paste('cellularAbundance_Karyotype.pdf',sep ='')
pdf(myfile)
print(cellularAbundance_Karyotype)
dev.off()

#estimate p values with anova
cellTypes <- c('B','CD4_T','CD8_T','HSCs_MPPs', 'Monocytes','NK','pDCs') 
for(idx in cellTypes){
  
  str <- paste('Test for ',idx,sep='')
  cat(str)
  cat('\n')
  tmp <- plot_data[plot_data$variable==idx ,]
  res <- aov(value ~ karyotype, data = tmp, exact = FALSE)
  posthoc <- TukeyHSD(x=res, 'karyotype', conf.level=0.95)
  print(posthoc)
  
}

#plot 3 based on FLT3
plot_data <- summary_patient_data[,c(4:10,13)]
plot_data <- melt(plot_data,id.vars = 'FLT3')
plot_data$FLT3 <- factor(plot_data$FLT3,levels = c('Wt','Mutated'))
plot_data <- na.omit(plot_data)

cellularAbundance_FLT3 <- ggplot(plot_data,aes(x=FLT3,y=(value),fill=FLT3),color=FLT3)+
  geom_beeswarm(size=0.78,alpha=0.8)+
  geom_boxplot(width=0.38 ,size=0.58,alpha=0.8,fatten=TRUE,
               position=position_dodge(width=0.78),
               outlier.colour = "black",
               outlier.shape = 21,
               outlier.fill = "gray",
               outlier.size = 0.4)+
  facet_wrap(~variable,scales = "free_x",ncol = 2)+
  theme_bw()+
  theme(axis.text.y = element_text( size = 12 ),
        axis.text.x = element_text(angle = 0, vjust = 1, hjust = 1, size = 12),
        axis.title.y = element_blank(),
        axis.title.x = element_text( size = 12 ),
        strip.text = element_text(size = 14,face='bold',lineheight=1),
        legend.position = "none",
        legend.text = element_text(size=12),
        strip.background = element_rect(colour = "black", fill = "white"),
        plot.title = element_text(size = 14, face = "bold",hjust = 0.5),aspect.ratio = 0.6)+
  ylab('Cells % in population')+
  scale_fill_manual(values=c('#f032e6','gray88'),name='')+
  theme(panel.spacing = unit(1.8, "lines"))+coord_flip()

#save the abundance plot
myfile <- paste('cellularAbundance_FLT3.pdf',sep ='')
pdf(myfile)
print(cellularAbundance_FLT3)
dev.off()

#estimate p values
cellTypes <- c('B','CD4_T','CD8_T','HSCs_MPPs', 'Monocytes','NK','pDCs') 
for(idx in cellTypes){
  
  str <- paste('Test for ',idx,sep='')
  cat(str)
  tmp <- plot_data[plot_data$variable==idx ,]
  res <- wilcox.test(value ~ FLT3, data = tmp, exact = FALSE)
  print(res)
}

####################################################################################################################################
#we can also perform Cox proportional Hazard's modeling 
#the cellular abundances are binarised to high and low based on median value
test <- summary_patient_data
a1 <- quantile(test$B)
test$B <- ifelse(test$B>a1[3],'High','Low')
test$B <- factor(test$B,levels = c('Low','High'))

a1 <- quantile(test$CD4_T)
test$CD4_T <- ifelse(test$CD4_T>a1[3],'High','Low')
test$CD4_T <- factor(test$CD4_T,levels = c('Low','High'))

a1 <- quantile(test$CD8_T)
test$CD8_T <- ifelse(test$CD8_T>a1[3],'High','Low')
test$CD8_T <- factor(test$CD8_T,levels = c('Low','High'))

a1 <- quantile(test$HSCs_MPPs)
test$HSCs_MPPs <- ifelse(test$HSCs_MPPs>a1[3],'High','Low')
test$HSCs_MPPs <- factor(test$HSCs_MPPs,levels = c('Low','High'))

a1 <- quantile(test$Monocytes)
test$Monocytes <- ifelse(test$Monocytes>a1[3],'High','Low')
test$Monocytes <- factor(test$Monocytes,levels = c('Low','High'))

a1 <- quantile(test$NK)
test$NK <- ifelse(test$NK>a1[3],'High','Low')
test$NK <- factor(test$NK,levels = c('Low','High'))

a1 <- quantile(test$pDCs)
test$pDCs <- ifelse(test$pDCs>a1[3],'High','Low')
test$pDCs <- factor(test$pDCs,levels = c('Low','High'))

library(survival)
library(survminer)

model_comb<-coxph(Surv(test$SurvivalTime,test$SurvivalStatus)~(B+CD4_T+CD8_T+HSCs_MPPs+Monocytes+NK+pDCs), 
                  data = test)

plot_comb <- ggforest(model_comb,data=test,noDigits = 2)
test.ph.comb <- cox.zph(model_comb)
s_comb <- ggcoxzph(test.ph.comb,font.main = 4,font.x=4,font.y=4)

myfile<-paste('Cox_model_cellular_abundance.pdf',sep='')
pdf(myfile,onefile = TRUE,width = 6,height = 5)
print(plot_comb)
dev.off()

myfile<-paste('Residuals_cellular_abundance.pdf',sep='')
pdf(myfile,onefile = TRUE)
print(s_comb)
dev.off()


myfile <- paste('cellularAbundance_patient_ctrl.pdf',sep ='')
knitr::include_graphics(myfile)

myfile <- paste('cellularAbundance_STS_LTS.pdf',sep ='')
knitr::include_graphics(myfile)

myfile <- paste('cellularAbundance_Karyotype.pdf',sep ='')
knitr::include_graphics(myfile)

myfile <- paste('cellularAbundance_FLT3.pdf',sep ='')
knitr::include_graphics(myfile)

#print(cellularAbundance_patient_ctrl)
#print(cellularAbundance_STS_LTS)
print(s_comb)
print(plot_comb)

```