Part3_DREMI.Rmd

---
title: "Feature generation using DREMI scores - Application to Benedicte's cohort"
author: "Dimitrios Kleftogiannis"
date: "2024-02-27"
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 generate features using DREMI. The data used here are described at "Early response evaluation by single cell signaling profiling in acute myeloid leukemia", doi: 10.1038/s41467-022-35624-4

The aim is to generate features that quantify the "strenght" of signaling relationships between proteins in the andibody panel and use them later on to predict long-term survivors from short-term survivors at time of diagnosis using machine learning based .

Please note that the DREMI code has been adopted from its original publication in Science "Conditional density-based analysis of T cell signaling in single-cell data" doi:10.1126/science.1250689, and it is translated to R using in-house functions. 

### 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.
We also define functions that are required
```{r load packages, echo=TRUE, eval=TRUE, error=TRUE, warning=FALSE,cache=TRUE}
library(ggplot2)
library(readxl)
library(dplyr)
library(ggridges)
library(MASS)
library(matrixStats)
library(reshape2)
library(ggthemes)
library(RColorBrewer)
library(ggbeeswarm)
library(DMwR)
library(xgboost)
library(pROC)

library(pracma)
library(matlab)
library(ggExtra)
library(minpack.lm)

#here select the channel name in our example:
remove_outliers <- function(x, na.rm = TRUE, ...) {
  qnt <- quantile(x, probs=c(0.02, 0.98), na.rm = na.rm, ...)
  H <- 1.5 * IQR(x, na.rm = na.rm)
  y <- x
  y[x < (qnt[1] - H)] <- NA
  y[x > (qnt[2] + H)] <- NA
  y
}

source('/Users/kleftogi/Desktop/Benedicte_data/find_data_cutoffs.R')
source('/Users/kleftogi/Desktop/Benedicte_data/compute_dremi_v2.R')
source('/Users/kleftogi/Desktop/Benedicte_data/pairwise_visualise_v2.R')
source('/Users/kleftogi/Desktop/Benedicte_data/delta_entropyweight_rawdata.R')

```

## 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")
set.seed(1234)
```

## Load the annotated patient data from Benedicte's cohort

We use DREMI score to assess all possible pairs of signaling markers. And this  

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

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

#retrieve the markers of interest - for this example we focus on signaling markers only
signal.col.names <- c("pAxl","CyclinB1","pNFkB","pErk","pSTAT1",
                      "pP38","pSTAT3","pCREB","pHist3","Casp3",
                      "pSTAT5","p4EBP1","pAkt","pRB","pS6")

patient_data <- patient_data[,c(signal.col.names,"Sample","cell_id","labels")]

#we generate all possible pairs of signaling markers in the panel
x = signal.col.names
testRelation <- with(subset(expand.grid(x,x),Var1!=Var2),paste0(Var1,'-',Var2))
#generate 210 features

#set limits for the number of cells allowed for DREMI score computation
Ncells_min <- 200
Ncells_max <- 10000
  

####################################################################################################

#load the clinical info from external files; we are not allowed to share these files please contact the authors
filename <- "AML_Benedicte_cohort/barcode_metadata_leukemia_cohort.xlsx"
#read the samples
md <- read_excel(filename)
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
metadata_filename <- "AML_Benedicte_cohort/leukemia_cohort_info.xlsx"
md_0h <- read_excel(metadata_filename)

start_time = Sys.time()
#save all DREMI features
DREMI_feature_vector <- data.frame()
for(currentPatient in 1:nrow(patient_samples)){
  
  myPatientID <- patient_samples[currentPatient,'patient_id']
  
  msg<-paste('Processing patient: ',myPatientID,'\n',sep='')
  cat(msg)
  
  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',currentPatient,sep='')
    data <- patient_data[patient_data$Sample==promptKey,]
    #parse the cell types one by one
    #estimate p values
    cellTypes <- c('B','CD4_T','CD8_T','HSCs_MPPs', 'Monocytes','NK','pDCs') 
    for(currentCellType in cellTypes){
      
      currentDataValues <- as.matrix(data[data$labels==currentCellType,signal.col.names])
      #only Ncells_max max are allowed
      if(nrow(currentDataValues)>=Ncells_max){
        a <- sample(nrow(currentDataValues))
        asel <- a[1:Ncells_max]
        currentDataValues <- currentDataValues[asel,]
      }
      
      #check for the lower limit
      if(nrow(currentDataValues)>=Ncells_min){
        
            featVector  <- matrix(-1,nrow = 1,ncol=length(testRelation))
            myK <- 1
            #parse all pairs of markers
            for(myRel in testRelation){
                
                if(myK %% 10 ==0){
                #str <- paste('    Relationship ',myRel,' ',myK,'/',length(testRelation),sep='')
                #cat(str)
                #cat('\n')
                }
                a1 <- strsplit(myRel, "-")
                a1 <- a1[[1]]
                markerX <- a1[1]
                markerY <- a1[2]
                Y <- currentDataValues[,markerY]
                X <- currentDataValues[,markerX]
                X <- remove_outliers(X)
                Y <- remove_outliers(Y)
                mB<- as.matrix(cbind(X,Y))
                mB <- na.omit(mB)
                test1 <- max(mB[,1], na.rm = TRUE)
                test2 <- max(mB[,2],na.rm=TRUE)
                  if(round(test1)==0 | round(test2)==0){
                  #if one of the markers has zero max value we return dremi equal to 0  
                  dremi <- runif(1,0,0.0001)
                  }else{
                    myLabel <- promptKey
                    maxy <- 0;
                    #mB[,1],mB[,2
                    myCutoffs <- find_data_cutoffs(mB, 50, 255)
                    minx1 <- myCutoffs[[1]]
                    miny1 <- myCutoffs[[2]]
                    maxx1 <- myCutoffs[[3]]
                    maxy1 <- myCutoffs[[4]]
                    if(maxy1>maxy){
                      maxy <- maxy1;
                    }
                    noise_threshold <- 0.8;
                    set_maxy <- 0
                    makePlot <- 0
                    dremi <- compute_dremi_v2(mB,markerX,markerY,noise_threshold,set_maxy,maxy,makePlot,myLabel)
                    dremi <- dremi[[1]]
                  }
                  featVector[1,myK]<- dremi
                  myK <- myK + 1
            }
        }else{
        #here write the mock features that are everywhere -1
        featVector  <- matrix(-1,nrow = 1,ncol=length(testRelation))
        }
          if(survivalTime$`5-year survival (days)`=='Alive'){
              survivalTimeValue <- 5*365
          }else{
              survivalTimeValue <- survivalTime$`5-year survival (days)`
          }
    
          DT <- data.frame(featVector)
          colnames(DT) <- testRelation
          DT$PatientNr <- promptKey
          DT$PatientID <- myPatientID$patient_id
          DT$SurvivalTime <- survivalTimeValue
          DT$karyotype <- chromAbber
          DT$FLT3 <- mutFLT3
          DT$CellType <- currentCellType
          DREMI_feature_vector <- rbind(DREMI_feature_vector,DT)
      }#closing for-loop for cell type
  }
  if(currentPatient %% 10==0){
    gc()
  }
}
end_time = Sys.time()
end_time-start_time

#add survival data and status STS vs. LTS
DREMI_feature_vector$SurvivalTime <- as.numeric(DREMI_feature_vector$SurvivalTime)
DREMI_feature_vector$SurvivalStatus <- ifelse(DREMI_feature_vector$SurvivalTime<5*365,1,0)
#first stratify patients based on survival time and compare their relative cellular abundances
DREMI_feature_vector$SurvivalGroup <- ifelse(DREMI_feature_vector$SurvivalTime<5*365,'STS','LTS')

#save the data
#save(DREMI_feature_vector,file='DREMI_feature_vector_v2.RDa')
```

## Predict survival using XGBoost and DREMI scores 
To predict short from long term survival using as input pairwise DREMI scores.

Same as in the previous comparison analysis, since the data are imbalanced with STS class having almost double the size of LTS we will also use the Synthetic Minority Oversampling Technique, often abbreviated SMOTE to generate synthetic datasets with different ratios between positive and negative classes.

```{r predicting survival using XGBoost and SMOTE on DREMI,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE}

####################################################################################################
#xgboost
#remove CyclinB1,pP38,pS6 because they have 0 median everywhere

load('DREMI_feature_vector_v2.RDa')

#perform cell type specific modeling
cellTypes <- c('B','CD4_T','CD8_T','HSCs_MPPs', 'Monocytes','NK','pDCs') 
totalPerformanceXGBoostDREMI <- data.frame()

#repeat the same multiple times with random sets, this is to assess reproducibility using independent training,validation and testing sets
N <- 100
for(iter in 1:N){
    
         for(currentCellType in cellTypes){
          #fetch the data
          currentDataValues <- DREMI_feature_vector[DREMI_feature_vector$CellType==currentCellType,]
          #first we check the the cases where there are -1 in the dataset.this means that DREMI was not computed so we filter out those samples
          a <- which(currentDataValues[,1]==-1)
          if(length(a)<=0){
            a <- 1
          }else{
            currentDataValues <- currentDataValues[-a,]
          }
          
          
          #generate random disjoint sets for training, validation and testing
          cut_point <- round(nrow(currentDataValues)*0.8)
          a <- sample(nrow(currentDataValues))
          asel <- a[1:cut_point]
          train_valid_Set <- currentDataValues[asel,]
          asel <- a[(cut_point+1):length(a)]
          testSet <- currentDataValues[asel,]
          
          #if the test set has only CellType from one class AUC cannot be built, thus we ship
          check1 <- which(testSet$SurvivalStatus==0)
          check2 <- which(testSet$SurvivalStatus==1)
          if(length(check1)>0 & length(check2)>0){
                     
                  #split the train_valid_Set to validation and actual train with ratios again 60-40
                  cut_point <- round(nrow(train_valid_Set)*0.6)
                  a <- sample(nrow(train_valid_Set))
                  asel <- a[1:cut_point]
                  trainSet <- train_valid_Set[asel,]
                  asel <- a[(cut_point+1):length(a)]
                  validationSet <- train_valid_Set[asel,]
                  
                  depthTunning <- data.frame()
                  for(depthVal in c(2,4,6,8,10,12,14,16,18,20,30,40,50,60)){
                    
                     Model <- xgboost(data = as.matrix(trainSet[,testRelation]),
                                    label = trainSet$SurvivalStatus,
                                    max_depth = depthVal,
                                    eta = 0.3,
                                    nthread = 2,
                                    nrounds = 4,
                                    objective = "binary:logistic",
                                    verbose = 0)
                    
                     pred <- predict(Model, as.matrix(validationSet[,testRelation]))
                     prediction <- as.numeric(pred > 0.5)
                     
                     #if you want to estimate error
                     #err <- mean(as.numeric(pred > 0.5) != validationSet$SurvivalStatus)
                     #print(paste("test-error=", err))
                     #assess performance
                     prediction <- factor(prediction,levels = c(0,1))
                     validationSet$SurvivalStatus <- factor(validationSet$SurvivalStatus,levels = c(0,1))
                     perf <- confusionMatrix(prediction, validationSet$SurvivalStatus, positive = c('1'))
                     perf <- as.data.frame(t(perf$byClass))
                     perf$depthVal <- depthVal
                     depthTunning <- rbind(depthTunning,perf)
                  }
                  #find the depthVal that maximises F1 in the validation set
                 a <- which.max(depthTunning$F1)
                 if(length(a)>0){
                   bestdepthVal <- depthTunning[a,'depthVal']
                 }else{
                   bestdepthVal <- 4
                 }
                 
                 #final model after tuning, tested on the test set
                 bstModel <- xgboost(data = as.matrix(train_valid_Set[,testRelation]),
                                    label = train_valid_Set$SurvivalStatus,
                                    max_depth = bestdepthVal,
                                    eta = 0.3,
                                    nthread = 2,
                                    nrounds = 4,
                                    objective = "binary:logistic",
                                    verbose = 0)
                    
                 pred <- predict(bstModel, as.matrix(testSet[,testRelation]))
                 prediction <- as.numeric(pred > 0.5)
                 prediction <- factor(prediction,levels = c(0,1))
                 testSet$SurvivalStatus <- factor(testSet$SurvivalStatus,levels = c(0,1))
                 perf <- confusionMatrix(prediction, testSet$SurvivalStatus, positive = c('1'))
                 perf <- as.data.frame(t(perf$byClass)) 
                 
                 #ROC using pROC package
                 tS <- as.numeric(testSet$SurvivalStatus)-1
                 pR <- as.numeric(prediction)-1
                 a <- roc(tS, pR,levels=c("1","0"))
                
                 perf$AUC <- a$auc     
                 perf$depthVal <- bestdepthVal
                 perf$CellType <- currentCellType
                 perf$Iter <- iter
                 totalPerformanceXGBoostDREMI <- rbind(totalPerformanceXGBoostDREMI,perf)
          }
      }
}
  

####################################################################################################
#SMOTE XGBoost SMOTE 1:1 ratio
totalPerformanceXGBoostDREMI_1_1 <- data.frame()
#repeat the same multiple times with random sets, this is to assess reproducibility using independent training,validation and testing sets
N <- 100
for(iter in 1:N){
    
         for(currentCellType in cellTypes){
          #fetch the data
          currentDataValues <- DREMI_feature_vector[DREMI_feature_vector$CellType==currentCellType,]
          #first we check the the cases where there are -1 in the dataset.this means that DREMI was not computed so we filter out those samples
          a <- which(currentDataValues[,1]==-1)
          if(length(a)<=0){
            a <- 1
          }else{
            currentDataValues <- currentDataValues[-a,]
          }
          
          #fetch the data and filter to perform SMOTE sampling
          currentDataValues_filt <- currentDataValues[,c(1:210,218)]
          currentDataValues_filt$SurvivalGroup <- as.factor(currentDataValues_filt$SurvivalGroup)
          #1000 vs 110 gives balanced data 
          currentDataValues_new <- SMOTE(SurvivalGroup ~ ., currentDataValues_filt, perc.over = 1000, perc.under = 110)
          currentDataValues_new$SurvivalStatus <- ifelse(currentDataValues_new$SurvivalGroup=='STS',1,0)

          #generate random disjoint sets for training, validation and testing
          cut_point <- round(nrow(currentDataValues_new)*0.8)
          a <- sample(nrow(currentDataValues_new))
          asel <- a[1:cut_point]
          train_valid_Set <- currentDataValues_new[asel,]
          asel <- a[(cut_point+1):length(a)]
          testSet <- currentDataValues_new[asel,]
          
          #if the test set has only CellType from one class AUC cannot be built, thus we ship
          check1 <- which(testSet$SurvivalStatus==0)
          check2 <- which(testSet$SurvivalStatus==1)
          if(length(check1)>0 & length(check2)>0){
                     
                  #split the train_valid_Set to validation and actual train with ratios again 60-40
                  cut_point <- round(nrow(train_valid_Set)*0.6)
                  a <- sample(nrow(train_valid_Set))
                  asel <- a[1:cut_point]
                  trainSet <- train_valid_Set[asel,]
                  asel <- a[(cut_point+1):length(a)]
                  validationSet <- train_valid_Set[asel,]
                  
                  depthTunning <- data.frame()
                  for(depthVal in c(2,4,6,8,10,12,14,16,18,20,30,40,50,60)){
                    
                     Model <- xgboost(data = as.matrix(trainSet[,testRelation]),
                                    label = trainSet$SurvivalStatus,
                                    max_depth = depthVal,
                                    eta = 0.3,
                                    nthread = 2,
                                    nrounds = 4,
                                    objective = "binary:logistic",
                                    verbose = 0)
                    
                     pred <- predict(Model, as.matrix(validationSet[,testRelation]))
                     prediction <- as.numeric(pred > 0.5)
                     
                     #if you want to estimate error
                     #err <- mean(as.numeric(pred > 0.5) != validationSet$SurvivalStatus)
                     #print(paste("test-error=", err))
                     #assess performance
                     prediction <- factor(prediction,levels = c(0,1))
                     validationSet$SurvivalStatus <- factor(validationSet$SurvivalStatus,levels = c(0,1))
                     perf <- confusionMatrix(prediction, validationSet$SurvivalStatus, positive = c('1'))
                     perf <- as.data.frame(t(perf$byClass))
                     perf$depthVal <- depthVal
                     depthTunning <- rbind(depthTunning,perf)
                  }
                  #find the depthVal that maximises F1 in the validation set
                 a <- which.max(depthTunning$F1)
                 if(length(a)>0){
                   bestdepthVal <- depthTunning[a,'depthVal']
                 }else{
                   bestdepthVal <- 4
                 }
                 
                 #final model after tuning, tested on the test set
                 bstModel <- xgboost(data = as.matrix(train_valid_Set[,testRelation]),
                                    label = train_valid_Set$SurvivalStatus,
                                    max_depth = bestdepthVal,
                                    eta = 0.3,
                                    nthread = 2,
                                    nrounds = 4,
                                    objective = "binary:logistic",
                                    verbose = 0)
                    
                 pred <- predict(bstModel, as.matrix(testSet[,testRelation]))
                 prediction <- as.numeric(pred > 0.5)
                 prediction <- factor(prediction,levels = c(0,1))
                 testSet$SurvivalStatus <- factor(testSet$SurvivalStatus,levels = c(0,1))
                 perf <- confusionMatrix(prediction, testSet$SurvivalStatus, positive = c('1'))
                 perf <- as.data.frame(t(perf$byClass)) 
                 
                 #ROC using pROC package
                 tS <- as.numeric(testSet$SurvivalStatus)-1
                 pR <- as.numeric(prediction)-1
                 a <- roc(tS, pR,levels=c("1","0"))
                
                 perf$AUC <- a$auc     
                 perf$depthVal <- bestdepthVal
                 perf$CellType <- currentCellType
                 perf$Iter <- iter
                 totalPerformanceXGBoostDREMI_1_1 <- rbind(totalPerformanceXGBoostDREMI_1_1,perf)
          }
      }
}

gc()
####################################################################################################
#SMOTE XGBoost 1:2 ratio
totalPerformanceXGBoostDREMI_1_2 <- data.frame()
N <- 100
for(iter in 1:N){
    
         for(currentCellType in cellTypes){
          #fetch the data
          currentDataValues <- DREMI_feature_vector[DREMI_feature_vector$CellType==currentCellType,]
          #first we check the the cases where there are -1 in the dataset.this means that DREMI was not computed so we filter out those samples
          a <- which(currentDataValues[,1]==-1)
          if(length(a)<=0){
            a <- 1
          }else{
            currentDataValues <- currentDataValues[-a,]
          }
          
          #fetch the data and filter to perform SMOTE sampling
          currentDataValues_filt <- currentDataValues[,c(1:210,218)]
          currentDataValues_filt$SurvivalGroup <- as.factor(currentDataValues_filt$SurvivalGroup)
          #1000 vs 110 gives balanced data 
          currentDataValues_new <- SMOTE(SurvivalGroup ~ ., currentDataValues_filt, perc.over = 1000, perc.under = 220)
          currentDataValues_new$SurvivalStatus <- ifelse(currentDataValues_new$SurvivalGroup=='STS',1,0)

          #generate random disjoint sets for training, validation and testing
          cut_point <- round(nrow(currentDataValues_new)*0.8)
          a <- sample(nrow(currentDataValues_new))
          asel <- a[1:cut_point]
          train_valid_Set <- currentDataValues_new[asel,]
          asel <- a[(cut_point+1):length(a)]
          testSet <- currentDataValues_new[asel,]
          
          #if the test set has only CellType from one class AUC cannot be built, thus we ship
          check1 <- which(testSet$SurvivalStatus==0)
          check2 <- which(testSet$SurvivalStatus==1)
          if(length(check1)>0 & length(check2)>0){
                     
                  #split the train_valid_Set to validation and actual train with ratios again 60-40
                  cut_point <- round(nrow(train_valid_Set)*0.6)
                  a <- sample(nrow(train_valid_Set))
                  asel <- a[1:cut_point]
                  trainSet <- train_valid_Set[asel,]
                  asel <- a[(cut_point+1):length(a)]
                  validationSet <- train_valid_Set[asel,]
                  
                  depthTunning <- data.frame()
                  for(depthVal in c(2,4,6,8,10,12,14,16,18,20,30,40,50,60)){
                    
                     Model <- xgboost(data = as.matrix(trainSet[,testRelation]),
                                    label = trainSet$SurvivalStatus,
                                    max_depth = depthVal,
                                    eta = 0.3,
                                    nthread = 2,
                                    nrounds = 4,
                                    objective = "binary:logistic",
                                    verbose = 0)
                    
                     pred <- predict(Model, as.matrix(validationSet[,testRelation]))
                     prediction <- as.numeric(pred > 0.5)
                     
                     #if you want to estimate error
                     #err <- mean(as.numeric(pred > 0.5) != validationSet$SurvivalStatus)
                     #print(paste("test-error=", err))
                     #assess performance
                     prediction <- factor(prediction,levels = c(0,1))
                     validationSet$SurvivalStatus <- factor(validationSet$SurvivalStatus,levels = c(0,1))
                     perf <- confusionMatrix(prediction, validationSet$SurvivalStatus, positive = c('1'))
                     perf <- as.data.frame(t(perf$byClass))
                     perf$depthVal <- depthVal
                     depthTunning <- rbind(depthTunning,perf)
                  }
                  #find the depthVal that maximises F1 in the validation set
                 a <- which.max(depthTunning$F1)
                 if(length(a)>0){
                   bestdepthVal <- depthTunning[a,'depthVal']
                 }else{
                   bestdepthVal <- 4
                 }
                 
                 #final model after tuning, tested on the test set
                 bstModel <- xgboost(data = as.matrix(train_valid_Set[,testRelation]),
                                    label = train_valid_Set$SurvivalStatus,
                                    max_depth = bestdepthVal,
                                    eta = 0.3,
                                    nthread = 2,
                                    nrounds = 4,
                                    objective = "binary:logistic",
                                    verbose = 0)
                    
                 pred <- predict(bstModel, as.matrix(testSet[,testRelation]))
                 prediction <- as.numeric(pred > 0.5)
                 prediction <- factor(prediction,levels = c(0,1))
                 testSet$SurvivalStatus <- factor(testSet$SurvivalStatus,levels = c(0,1))
                 perf <- confusionMatrix(prediction, testSet$SurvivalStatus, positive = c('1'))
                 perf <- as.data.frame(t(perf$byClass)) 
                 
                 #ROC using pROC package
                 tS <- as.numeric(testSet$SurvivalStatus)-1
                 pR <- as.numeric(prediction)-1
                 a <- roc(tS, pR,levels=c("1","0"))
                
                 perf$AUC <- a$auc     
                 perf$depthVal <- bestdepthVal
                 perf$CellType <- currentCellType
                 perf$Iter <- iter
                 totalPerformanceXGBoostDREMI_1_2 <- rbind(totalPerformanceXGBoostDREMI_1_2,perf)
          }
      }
}

gc()


####################################################################################################
#SMOTE XGBoost 1:3 ratio
totalPerformanceXGBoostDREMI_1_3 <- data.frame()
N <- 100
for(iter in 1:N){
    
         for(currentCellType in cellTypes){
          #fetch the data
          currentDataValues <- DREMI_feature_vector[DREMI_feature_vector$CellType==currentCellType,]
          #first we check the the cases where there are -1 in the dataset.this means that DREMI was not computed so we filter out those samples
          a <- which(currentDataValues[,1]==-1)
          if(length(a)<=0){
            a <- 1
          }else{
            currentDataValues <- currentDataValues[-a,]
          }
          
          #fetch the data and filter to perform SMOTE sampling
          currentDataValues_filt <- currentDataValues[,c(1:210,218)]
          currentDataValues_filt$SurvivalGroup <- as.factor(currentDataValues_filt$SurvivalGroup)
          #1000 vs 110 gives balanced data 
          currentDataValues_new <- SMOTE(SurvivalGroup ~ ., currentDataValues_filt, perc.over = 1000, perc.under = 330)
          currentDataValues_new$SurvivalStatus <- ifelse(currentDataValues_new$SurvivalGroup=='STS',1,0)

          #generate random disjoint sets for training, validation and testing
          cut_point <- round(nrow(currentDataValues_new)*0.8)
          a <- sample(nrow(currentDataValues_new))
          asel <- a[1:cut_point]
          train_valid_Set <- currentDataValues_new[asel,]
          asel <- a[(cut_point+1):length(a)]
          testSet <- currentDataValues_new[asel,]
          
          #if the test set has only CellType from one class AUC cannot be built, thus we ship
          check1 <- which(testSet$SurvivalStatus==0)
          check2 <- which(testSet$SurvivalStatus==1)
          if(length(check1)>0 & length(check2)>0){
                     
                  #split the train_valid_Set to validation and actual train with ratios again 60-40
                  cut_point <- round(nrow(train_valid_Set)*0.6)
                  a <- sample(nrow(train_valid_Set))
                  asel <- a[1:cut_point]
                  trainSet <- train_valid_Set[asel,]
                  asel <- a[(cut_point+1):length(a)]
                  validationSet <- train_valid_Set[asel,]
                  
                  depthTunning <- data.frame()
                  for(depthVal in c(2,4,6,8,10,12,14,16,18,20,30,40,50,60)){
                    
                     Model <- xgboost(data = as.matrix(trainSet[,testRelation]),
                                    label = trainSet$SurvivalStatus,
                                    max_depth = depthVal,
                                    eta = 0.3,
                                    nthread = 2,
                                    nrounds = 4,
                                    objective = "binary:logistic",
                                    verbose = 0)
                    
                     pred <- predict(Model, as.matrix(validationSet[,testRelation]))
                     prediction <- as.numeric(pred > 0.5)
                     
                     #if you want to estimate error
                     #err <- mean(as.numeric(pred > 0.5) != validationSet$SurvivalStatus)
                     #print(paste("test-error=", err))
                     #assess performance
                     prediction <- factor(prediction,levels = c(0,1))
                     validationSet$SurvivalStatus <- factor(validationSet$SurvivalStatus,levels = c(0,1))
                     perf <- confusionMatrix(prediction, validationSet$SurvivalStatus, positive = c('1'))
                     perf <- as.data.frame(t(perf$byClass))
                     perf$depthVal <- depthVal
                     depthTunning <- rbind(depthTunning,perf)
                  }
                  #find the depthVal that maximises F1 in the validation set
                 a <- which.max(depthTunning$F1)
                 if(length(a)>0){
                   bestdepthVal <- depthTunning[a,'depthVal']
                 }else{
                   bestdepthVal <- 4
                 }
                 
                 #final model after tuning, tested on the test set
                 bstModel <- xgboost(data = as.matrix(train_valid_Set[,testRelation]),
                                    label = train_valid_Set$SurvivalStatus,
                                    max_depth = bestdepthVal,
                                    eta = 0.3,
                                    nthread = 2,
                                    nrounds = 4,
                                    objective = "binary:logistic",
                                    verbose = 0)
                    
                 pred <- predict(bstModel, as.matrix(testSet[,testRelation]))
                 prediction <- as.numeric(pred > 0.5)
                 prediction <- factor(prediction,levels = c(0,1))
                 testSet$SurvivalStatus <- factor(testSet$SurvivalStatus,levels = c(0,1))
                 perf <- confusionMatrix(prediction, testSet$SurvivalStatus, positive = c('1'))
                 perf <- as.data.frame(t(perf$byClass)) 
                 
                 #ROC using pROC package
                 tS <- as.numeric(testSet$SurvivalStatus)-1
                 pR <- as.numeric(prediction)-1
                 a <- roc(tS, pR,levels=c("1","0"))
                
                 perf$AUC <- a$auc     
                 perf$depthVal <- bestdepthVal
                 perf$CellType <- currentCellType
                 perf$Iter <- iter
                 totalPerformanceXGBoostDREMI_1_3 <- rbind(totalPerformanceXGBoostDREMI_1_3,perf)
          }
      }
}

gc()



####################################################################################################
####################################################################################################
#aggregate the performances and visualise

method_colors <- c('olivedrab2',"#4DB23B","#066970","#b0c6a2")
metric_colors <-  c('dodgerblue','khaki3',"#FFD258",'brown2')


  totalPerformanceXGBoostDREMI$Method <- 'XGBoost'
  totalPerformanceXGBoostDREMI_1_1$Method <- 'XGBoost_SMOTE_1:1'
  totalPerformanceXGBoostDREMI_1_2$Method <- 'XGBoost_SMOTE_1:2'
  totalPerformanceXGBoostDREMI_1_3$Method <- 'XGBoost_SMOTE_1:3'
  
  #plot #1
  tmp1 <- rbind(totalPerformanceXGBoostDREMI,totalPerformanceXGBoostDREMI_1_1,
                    totalPerformanceXGBoostDREMI_1_2,totalPerformanceXGBoostDREMI_1_3)
  plotData <- tmp1[,c("Sensitivity","Specificity","F1","AUC","CellType","Method")] 
  colnames(plotData)[1:4] <- c('SEN','SPE','F1','AUC')
  tmp <- melt(plotData,id.vars = c("CellType","Method" ))
  tmp$Method <- factor(tmp$Method,levels = c('XGBoost','XGBoost_SMOTE_1:1','XGBoost_SMOTE_1:2','XGBoost_SMOTE_1:3'))
  colnames(tmp)[3] <- 'Metric'
  tmp$Metric <- factor(tmp$Metric,levels = c('SEN','SPE','F1','AUC'))
  
  performance_Medians_DREMI <- ggplot(tmp, aes(x = Metric, y = value, fill = Method)) +
              #geom_jitter(height=0,width=0.1, size=0.48,alpha=0.8)+
              geom_boxplot(width=0.88,
                           alpha=0.6,
                           size=0.2,
                    outlier.colour = "black",
                    outlier.shape = NA,
                    outlier.fill = "red",
                    outlier.size = 0.05,
                    notch = FALSE,color='black')+
              facet_wrap(~CellType,ncol = 2)+
              ggtitle('DREMI features')+
              theme_few() +
              ylim(0,1)+
              scale_fill_manual(values = method_colors,name='')+
              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 = 10,face='bold',lineheight=1),
                      legend.position = "bottom",aspect.ratio = 0.4,
                      legend.text=element_text(size=12))+
              guides(fill = guide_legend(override.aes = list(size=1.4),ncol=2,title=""))
    
   
  #save the plot
  myfile <- paste('performance_Medians_DREMI.pdf',sep ='')
  pdf(myfile)
  print(performance_Medians_DREMI)
  dev.off()

  meanPerformanceDREMI <- data.frame()
  #compute the median per cell type and per method
  for(currentCellType in cellTypes){
   
    tmp <-  plotData[plotData$CellType==currentCellType,]
    #compute the medians
    tmp <- tmp %>%
    group_by(Method) %>%
    dplyr::summarize(meanSEN = mean(SEN, na.rm=TRUE),
                     meanSPE = mean(SPE, na.rm=TRUE),
                     meanF1 = mean(F1, na.rm=TRUE),
                     meanAUC = mean(AUC, na.rm=TRUE))
    tmp$CellType <- currentCellType
    meanPerformanceDREMI <- rbind(meanPerformanceDREMI,tmp)
  }
  
  tmp <- melt(meanPerformanceDREMI,id.vars = c("CellType","Method" ))
  tmp$Method <- factor(tmp$Method,levels = c('XGBoost','XGBoost_SMOTE_1:1','XGBoost_SMOTE_1:2','XGBoost_SMOTE_1:3'))
  colnames(tmp)[3] <- 'Metric'
  tmp$Metric <- factor(tmp$Metric,levels = c('meanSEN','meanSPE','meanF1','meanAUC'))
  
  performance_summary_DREMI <- ggplot(tmp, aes(x = CellType, y = value, fill = Metric)) +
              geom_bar(stat="identity",size=0.3,width = 0.5,color='black')+
              facet_wrap(~Method+Metric,ncol = 4)+
              geom_hline(yintercept = 0.75,color='red')+
              ggtitle('DREMI features')+
              theme_few() +
              ylim(0,1)+
              scale_fill_manual(values = metric_colors,name='')+
              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 = 10,face='bold',lineheight=1),
                      legend.position = "none",aspect.ratio = 0.5,
                      legend.text=element_text(size=12))

  #save the  plot
  myfile <- paste('performance_summary_DREMI.pdf',sep ='')
  pdf(myfile)
  print(performance_summary_DREMI)
  dev.off()
  
  #save(meanPerformanceDREMI,file='meanPerformanceDREMI.RDa')
  
```

## Visualise the results from the comparison analysis using DREMI
```{r show comparison figures,echo=TRUE,eval=TRUE,error=TRUE, warning=FALSE,cache=TRUE,fig.align="center",out.width = "600px",out.height = "600px"}

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

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



#find the best ranked method based on average performance across all cell types
  tmp <- meanPerformanceDREMI
  tmp <- tmp %>%
    group_by(Method) %>%
    dplyr::summarize(meanSEN = mean(meanSEN, na.rm=TRUE),
                     meanSPE = mean(meanSPE, na.rm=TRUE),
                     meanF1 = mean(meanF1, na.rm=TRUE),
                     meanAUC = mean(meanAUC, na.rm=TRUE))
  tmp


```