---
title: "pipeComp evaluation"
author: "Won-Min Song"
date: '2024-12-13'
output: html_document
---

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

#### Comparative analysis of MSC and other benchmark clustering methods
This is the R script used to perform comparative analysis of MSC with other benchmark methods with pipeComp data sets. These data sets have gold-standard ground-truth clusters, and the raw count matrix and pre-processed data by Seurat workflow can be downloaded from Synpase with DOI, [https://doi.org/10.7303/syn52966803](https://doi.org/10.7303/syn52966803), under 'pipeComp' folder. Make sure to download 'pipeComp_functions.R' to the working directory from the synapse folder to run the codes. 

[Section 1. pipeComp data pre-processing](#S1): The 'Data' subfolder contains the raw count matrices in .RDS format, and the pre-processed data is in 'quality_controlled_seurat.results.RDS' in the main folder in .RDS format, a gzipped format of R objects. Placing  'quality_controlled_seurat.results.RDS' in the working folder will bypass [Section 1](#S1).;\
[Section 2. Calculate clusters using benchmark methods](#S2): The main folder contains *results..RDS files to hold results for individual benchmark methods and MSCs. Place these results.RDS files under the working directory.;\
[Section 3. Calculate clusters using MSC](#S3): The main folder contains .RDS files to hold MSC results from Pearson's correlation (iterative_MSC.correlation.results.RDS) and Euclidean distance (iterative_MSC.euclidean.results.RDS). Also, this section will generate results from the single-layer split by AdaptSplit without performing the iterative top-down approach. ;\
[Section 4. Comparative evaluation of discrete partitions using adjusted Rand Index](#S4): This section reproduces the figures to calculate the adjusted Rand Index on single-layer clustering results including CIDR, SC3, SNN-based Louvain clustering and AdaptSplit;\
[Section 5. Comparative evaluation of iterative top-down clustering results](#S5): This section reproduces the figures to calculate inclusion, coverage rates and detection accuracy for performance metrics in the gold standard clusters.;

## [Section 1. pipeComp data pre-processing]{#S1}

```{r pre-processing, message = FALSE}
rm(list = ls())

library(Seurat)
library(scater)
library(scran)
library(cowplot)
library(CIDER)
library(cidr)
library(mclust)
library(MSC)
library(Matrix)
library(doParallel)
library(igraph)
library(SC3)
library(infotheo)
library(ggtext)
library(ggplot2)
library(Cairo)



#### Set up folders
root.dir <- "/home/won-min/Documents/R_packages/Single_Cell_Clustering/pipeComp_v4";setwd(root.dir)
out.dir <- getwd()
source("pipeComp_functions.R")

min.size = 5 # minimum cell cluster size
n.cores = 6 # number of cores to use

run.preprocess = FALSE # set this FALSE if 'quality_controlled_seurat.results.RDS' is already downloaded in the working folder. 
run.clustering = FALSE # set FALSE if want to bypass the clustering process. 

#### load pipeComp data collection and use Seurat workflow for pre-processing. 
if (run.preprocess)
{
  #############
  simfiles <- list.files(path = "pipeComp/Data",full.names = TRUE)
  names(simfiles) = gsub("^(.*)/|\\.SCE\\.rds","",simfiles)
  seu.list = vector("list",length(simfiles));names(seu.list) = names(simfiles)
  for (i in 1:length(simfiles))
  {
    sce= readRDS(file = simfiles[i])
    seu = CreateSeuratObject(counts = counts(sce),meta.data = as.data.frame(colData(sce)))
    if (any(colnames(seu@meta.data) == "pldf.class")) 
    {
      seu = subset(seu,pldf.class == "singlet" & pct_counts_Mt <= 20 & total_features >= 100)
    }else{
      seu = subset(seu,pct_counts_Mt <= 20 & total_features >= 100)
    }
    seu <- NormalizeData(seu)
    #seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = 3000)
    seu <- ScaleData(seu)
    VariableFeatures(seu) = rownames(subset(qc.results[[i]]$gene.data,bio > tech & p.value < 0.05))
    seu = RunPCA(seu,npcs = 20)
    seu.list[[i]] = seu;
    rm(seu)
  }
  
  saveRDS(seu.list,file = "quality_controlled_seurat.results.RDS")
  
}else{
  seu.list = readRDS("quality_controlled_seurat.results.RDS")
}

```

## [Section 2. Calculate clusters using benchmark methods]{#S2}

```{r cluster analysis, message = FALSE}
##### Create functions for generating the benchmark clusters
setwd(root.dir)
run.benchmark.clustering = TRUE
run.methods = c("SNN" = FALSE,"CIDR" = FALSE,"SC3" = FALSE,"aKNNO" = FALSE,"scCAN" = FALSE,"RaceID3" = FALSE) # indicate if each needs to run
if (run.benchmark.clustering)
{
  # update SNN results
  if (run.methods["SNN"])
  {
    get_snn_clusters <- function(seu,resol.vec = seq(0.4,1.2,0.1))
    {
      seu <- FindNeighbors(seu)
      seu <- FindClusters(seu, resolution = resol.vec)
      seu.cls.col = colnames(seu@meta.data)[grep("snn_res",colnames(seu@meta.data))]
      
      snn.data = seu@meta.data[,seu.cls.col]
      rownames(snn.data) = colnames(seu)
      colnames(snn.data) = gsub("^(.*)_snn_res","seurat_snn_res",colnames(snn.data))
      snn.mods = lapply(snn.data,function(x) split(rownames(snn.data),factor(x)))
      
      for (i in 1:length(seu.cls.col))
      {
        names(snn.mods[[i]]) = paste0(names(snn.mods)[i],"_C",names(snn.mods[[i]]))
      }
      names(snn.mods) = NULL
      snn.mods = do.call('c',snn.mods)
      
      output = list(modules = snn.mods,cluster.matrix = snn.data)
      
      return(output)
    }
    
    #### run analysis: SNN
    snn.results = vector("list",length(simfiles));
    names(snn.results) = names(simfiles)
    for (i in 1:length(simfiles))
    {
      snn.results[[i]] = get_snn_clusters(seu.list[[i]],resol.vec = seq(0.4,1.2,0.4))
      
    }
    
    saveRDS(snn.results,file = "cross_scale_SNN.results.RDS")
  }else{
    snn.results = readRDS(file = "cross_scale_SNN.results.RDS")
  }
  
  # update CIDR and SC3
  if (run.methods["CIDR"] | run.methods["SC3"])
  {
  
    sc3.lst = cidr.lst = vector("list",length(seu.list));
    names(sc3.lst) = names(cidr.lst) = names(seu.list)
    for (i in 1:length(simfiles))
    {
      ###### generate benchmark clustering results
      seu = seu.list[[i]]
      nc = min(c(20,ncol(Embeddings(seu,"pca"))))
      bench.out = run_benchmark_clusters(seu,n.cores = n.cores,dimred = "pca",dims = 1:nc,sc3.seed = 1234)
      vec = bench.out@meta.data$SC3_clusters;names(vec) = colnames(seu)
      sc3.lst[[i]] = vec;
      rm(vec)
      
      #### add cidr results
      cidr.obj = run_cidr(tags = as.matrix(GetAssayData(seu,layer = "counts")))
      cidr.res = factor(cidr.obj@clusters);names(cidr.res) = colnames(seu)
      cidr.lst[[i]] = cidr.res
      rm(bench.out,seu,cidr.obj,cidr,res)
    }
    
    saveRDS(sc3.lst,file = "SC3.results.RDS")
    saveRDS(cidr.lst,file = "CIDR.results.RDS")
    
  }else{
    sc3.lst = readRDS(file = "SC3.results.RDS")
    cidr.lst = readRDS(file = "CIDR.results.RDS")
  }
  
  # run aKNNO
  if (run.methods["aKNNO"])
  {
    library(aKNNO)
    resol.vec = c(0.4,0.8,1.2)
    aKNNO.results = vector("list",length(seu.list))
    names(aKNNO.results) = names(simfiles)
    for (i in 1:length(seu.list))
    {
      seu <- FindNeighbors_aKNNO(seu.list[[i]])
      seu <- FindClusters(seu,graph.name="aKNN_O", resolution = resol.vec)
      seu.cls.col = colnames(seu@meta.data)[grep("aKNN_O_res",colnames(seu@meta.data))]
      
      snn.data = seu@meta.data[,seu.cls.col]
      rownames(snn.data) = colnames(seu)
      #colnames(snn.data) = gsub("^(.*)_snn_res","seurat_snn_res",colnames(snn.data))
      
      aKNNO.results[[i]] = snn.data
      
    }
    saveRDS(aKNNO.results,file = "cross_scale_aKNNO.results.RDS")
  }else{
    aKNNO.results = readRDS(file = "cross_scale_aKNNO.results.RDS")
  }
  
  # run scCAN
  if (run.methods["scCAN"])
  {
    library(scCAN)
    scCAN.results = vector("list",length(seu.list))
    names(scCAN.results) = names(seu.list)
    for (i in 1:length(seu.list))
    {
      cat(paste0("processing:",names(seu.list)[i],"\n"))
      cnt = GetAssayData(seu.list[[i]],layer = "counts")
      cnt = t(as.matrix(cnt))
      result <- scCAN(cnt)
      names(result$cluster) = rownames(cnt)
      rownames(result$latent) = rownames(cnt)
      scCAN.results[[i]] = result
      rm(cnt,result)
      gc()
    }
    saveRDS(scCAN.results,file = "scCAN.results.RDS")
  }else{
    scCAN.results = readRDS(file ="scCAN.results.RDS")
  }
  
  # run RaceID3
  if (run.methods["RaceID3"])
  {
    library(RaceID)
    RaceID3.results = vector("list",length(seu.list));
    names(RaceID3.results) = names(seu.list)
    for (i in 1:length(seu.list))
    {
      cat(paste0("processing RaceID3:",names(seu.list)[i],"\n"))
      nk = length(unique(seu.list[[i]]@meta.data$phenoid))
      sc = SCseq(GetAssayData(seu.list[[i]],layer = "counts"))
      sc = filterdata(sc,mintotal=2000)
      fdata <- getfdata(sc)
      sc <- compdist(sc,metric="pearson")
      sc <- clustexp(sc)
      RaceID3.results[[i]] = sc@cluster$kpart
      rm(sc,fdata,nk)
      gc()
    }
    saveRDS(RaceID3.results,file = "RaceID3.results.RDS")
  }else{
    RaceID3.results = readRDS(file = "RaceID3.results.RDS")
  }
}

```

## [Section 3. Calculate clusters using MSC]{#S3}
```{r MSC, message = FALSE}
setwd(root.dir)
run.MSC.clustering = TRUE
run.each = c("msc.euc" = F,"msc.cor" = F) # indicate if each needs to run. 
if (run.MSC.clustering)
{
  # load pre-processed file if section 1 wasn't run
  if (!any(ls() == "seu.list")) seu.list = readRDS("quality_controlled_seurat.results.RDS")
  
  if (run.each["msc.euc"])
  {
    lena.euclid.results = split.euclid.results = vector("list",length(seu.list))
    names(lena.euclid.results) = names(split.euclid.results) = names(seu.list)
    
    for (i in 1:length(seu.list))
    {
      cat(paste0("processing:",names(seu.list)[i],"\n"))
      mat = as.matrix(t(Embeddings(seu.list[[i]],"pca")))
      if (T)
      {
        # get LEN
        len.out = generate_cell_network.wt.loess(
          mat,
          bcnt = NULL,
          dist.func = function(x,y) sqrt(sum((x-y)^2)),
          is.decreasing = FALSE,
          n.cores = 4
        );
        
        
        # run MSC
        d.func = function(x) x
        valid.gate = "both"
        
        cout = components(len.out)
        ci = which(cout$csize > log(vcount(len.out)))
        glst = lapply(ci,function(n) induced_subgraph(graph = len.out,vids = names(cout$membership)[cout$membership == n]))
        
        alpha.res = sapply(glst,function(gi) identify_scale_parameter(g = gi,nv = 20,nr.max = 3,seed = 1234,d.func = d.func,mode = "diameter"))[1]
        gc()
        
        iter.res = iterative_clustering(g.in = len.out,min.size = 5,d.func = d.func,alpha = alpha.res)
        
        lena.euclid.results[[i]] = iter.res
        
      }else{
        iter.res = lena.euclid.results[[i]]
      }
      
      
      ### get adapt split results
      # find root
      gh = graph.data.frame(iter.res$module.table[,c("parent.cluster","cluster.name")],directed = T)
      deg = igraph::degree(gh,mode = "in");rid = names(deg)[deg == 0];
      pid = sapply(iter.res$modules,length)[rid];pid = names(pid)[which.max(pid)]
      vec = get_rooted_membership(mtbl = iter.res$module.table,modules = iter.res$modules,pid = pid,bg = colnames(mat))
      split.euclid.results[[i]] = vec
      rm(cout,cu,glst,alpha.res,iter.res,len.out)
    }
    
    saveRDS(lena.euclid.results,file = "iterative_MSC.euclidean.results.RDS")
    
  }else{
    lena.euclid.results = readRDS(file = "iterative_MSC.euclidean.results.RDS")
    split.euclid.results = lapply(lena.euclid.results,function(x) extract_first_split(x)[[2]])
    
  }
  
  # get results: MSC with correlation (** do a fresh run)
  if (run.each["msc.cor"])
  {
    lena.cor.results = split.cor.results = vector("list",length(seu.list))
    names(lena.cor.results) = names(split.cor.results) = names(seu.list)
    for (i in 1:length(seu.list))
    {
      cat(paste0("processing:",names(seu.list)[i],"\n"))
      
      ### run MSC clustering
      # for imputed data for correlation computation 
      feats = rownames(subset(qc.results[[i]]$gene.data.alra,bio > tech & p.value < 0.05))
      cat(paste0("- #. features:",length(feats),"\n"))
      dat = GetAssayData(seu.list[[i]],layer = "data")
      dat = dat[feats,]
      bcnt = GetAssayData(seu.list[[i]],layer = "counts")
      bcnt = bcnt[rownames(dat),]
      
      # get LEN
      len.out = generate_cell_network.wt.loess(
        mat = dat,
        bcnt = bcnt,
        dist.func = function(a, b) cor(a,b,method ="pearson"),
        is.decreasing = TRUE,
        n.cores = 4
      );
      
      # run MSC
      d.func = function(x) sqrt(2*(1-x))
      valid.gate = "both"
      
      cout = components(len.out)
      ci = which(cout$csize > log(vcount(len.out)))
      glst = lapply(ci,function(n) induced_subgraph(graph = len.out,vids = names(cout$membership)[cout$membership == n]))
      
      alpha.res = sapply(glst,function(gi) identify_scale_parameter(g = gi,nv = 20,nr.max = 3,seed = 1234,d.func = d.func,mode = "diameter"))[1]
      gc()
      
      iter.res = iterative_clustering(g.in = len.out,min.size = 5,d.func = d.func,alpha = alpha.res)
      
      lena.cor.results[[i]] = iter.res
      
      ### get adapt split results
      # find root
      gh = graph.data.frame(iter.res$module.table[,c("parent.cluster","cluster.name")],directed = T)
      deg = igraph::degree(gh,mode = "in");rid = names(deg)[deg == 0];
      pid = sapply(iter.res$modules,length)[rid];pid = names(pid)[which.max(pid)]
      vec = get_rooted_membership(mtbl = iter.res$module.table,modules = iter.res$modules,pid = pid,bg = colnames(seu.list[[i]]))
      split.cor.results[[i]] = vec
      rm(cout,cu,glst,alpha.res,iter.res,len.out)
    }
    
    saveRDS(lena.cor.results,file = "iterative_MSC.correlation.results.RDS")
    
  }else{
    lena.cor.results = readRDS(file = "iterative_MSC.correlation.results.RDS")
    split.cor.results = lapply(lena.cor.results,function(x) extract_first_split(x)[[2]])
    
  }
}

```

## [Section 4. Comparative evaluation of discrete partitions using adjusted Rand Index]{#S4}
```{r adjusted rand index, message = FALSE}
setwd(root.dir)
# load inputs
if (!any(ls() == "seu.list")) seu.list = readRDS("quality_controlled_seurat.results.RDS")

## collect cluster results
clus.data.list = vector("list",length(seu.list));
names(clus.data.list) = names(seu.list)
for (i in 1:length(seu.list))
{
  clus.data = data.frame(cell.id = colnames(seu.list[[i]]),phenoid = seu.list[[i]]@meta.data$phenoid)
  
  # add SNN results
  snn.res = snn.results[[i]]$cluster.matrix;colnames(snn.res) = gsub("^(.*)snn_res.","SNN \u03B3=",colnames(snn.res));
  snn.res = snn.res[match(clus.data[[1]],rownames(snn.res)),]
  clus.data = cbind.data.frame(clus.data,snn.res)
  
  # add aKNNO results
  aknn.res = aKNNO.results[[i]];colnames(aknn.res) = gsub("aKNN_O_res.","aKNNO \u03B3=",colnames(aknn.res));
  aknn.res = aknn.res[match(clus.data[[1]],rownames(aknn.res)),]
  clus.data = cbind.data.frame(clus.data,aknn.res)
  
  # add MSC results
  clus.data$MSC.Euc.1L = split.euclid.results[[i]][colnames(seu.list[[i]])]
  clus.data$MSC.Cor.1L = split.cor.results[[i]][colnames(seu.list[[i]])]
  
  # add SC3
  clus.data$SC3 = sc3.lst[[i]][match(clus.data[[1]],names(sc3.lst[[i]]))]
  
  # add CIDR
  clus.data$CIDR = cidr.lst[[i]][match(clus.data[[1]],names(cidr.lst[[i]]))]
  
  # add scCAN
  clus.data$scCAN = scCAN.results[[i]]$cluster[match(clus.data[[1]],names(scCAN.results[[i]]$cluster))]
  
  # add RaceID3
  clus.data$RaceID3 = RaceID3.results[[i]][match(clus.data[[1]],names(RaceID3.results[[i]]))]
  
  clus.data.list[[i]] = clus.data
  rm(clus.data,snn.res,aknn.res)
  gc()
}

# set name labels for visualization
nm = as.list(colnames(clus.data.list[[1]])[-c(1,2)])
nm[[7]] = bquote(MSC * scriptstyle(atop("EUC", "1L")))
nm[[8]] = bquote(MSC * scriptstyle(atop("COR", "1L")))
names(nm) = colnames(clus.data.list[[1]])[-c(1,2)]

nm = c(nm,list("MSC.Euc.ML" = bquote(MSC * scriptstyle(atop("EUC", "ML")))),list("MSC.Cor.ML" = bquote(MSC * scriptstyle(atop("COR", "ML")))))

##### evaluate adjusted rand index for first level splits
if (TRUE)
{
  library(mclust)
  adj.matrix = matrix(0,nrow = length(clus.data.list),ncol = ncol(clus.data.list[[1]])-2)
  rownames(adj.matrix) = names(clus.data.list)
  for (i in 1:length(clus.data.list))
  {
    clus.res.f = clus.data.list[[i]]
    adj.vec = sapply(clus.res.f[,3:ncol(clus.res.f)],function(x) adjustedRandIndex(x = clus.res.f$phenoid,y = x))
    adj.matrix[i,] = adj.vec
    if (i == 1) colnames(adj.matrix) = colnames(clus.res.f)[3:ncol(clus.res.f)]
  }
  
  # plot out results
  adjr.matrix = adj.matrix
  
  # rename methods for simpler names
  #colnames(adjr.matrix) = c("SNN r=0.4","SNN r=0.8","SNN r=1.2","SC3","CIDR","AdaptSplit:Euclidean","AdaptSplit:Correlation")
  #adjr.matrix = cbind(apply(adjr.matrix[,1:3],1,max),adjr.matrix)
  #colnames(adjr.matrix)[1] = "SNN Best"
  adjr.p = make_summary_plots(dat = adjr.matrix,nm = nm,ylab = "ARI")

}

##### evaluate NMI
if (TRUE)
{
  library(infotheo)
  mi.matrix = matrix(0,nrow = length(clus.data.list),ncol = ncol(clus.data.list[[1]])-2)
  rownames(mi.matrix) = names(clus.data.list)
  for (i in 1:length(clus.data.list))
  {
    clus.res.f = clus.data.list[[i]]
    t = clus.res.f$phenoid;t = match(t,unique(t))
    #mi.vec = sapply(clus.res.f[,3:ncol(clus.res.f)],function(x) mutinformation(X = clus.res.f$phenoid,Y = x))
    mi.vec = sapply(clus.res.f[,3:ncol(clus.res.f)],function(x) NMI(c = clus.res.f$phenoid,t = x))
    mi.matrix[i,] = mi.vec
    if (i == 1) colnames(mi.matrix) = colnames(clus.res.f)[3:ncol(clus.res.f)]
  }
  
  mi.matrix[is.nan(mi.matrix)] = 0
  
  mi.p = make_summary_plots(dat = mi.matrix,nm = nm,ylab = "NMI")
}

##### evaluate cluster entropy
if (TRUE)
{
  prt.matrix = acc.matrix = matrix(0,nrow = length(clus.data.list),ncol = ncol(clus.data.list[[1]])-2)
  rownames(prt.matrix) = rownames(acc.matrix) = names(clus.data.list)
  
  for (i in 1:length(clus.data.list))
  {
    cls.t = clus.data.list[[i]][[2]];names(cls.t) = clus.data.list[[i]][[1]]
    prtvec = sapply(clus.data.list[[i]][,-c(1,2)],function(x,y) get_cluster_purity(x,y)$overall.purity,y = cls.t)
    accvec = sapply(clus.data.list[[i]][,-c(1,2)],function(x,y) get_cluster_accuracy(x,y)$overall.accuracy,y = cls.t)
    
    if (i == 1) colnames(prt.matrix) = colnames(acc.matrix) = colnames(clus.res.f)[3:ncol(clus.res.f)]
    
    prt.matrix[i,] = prtvec
    acc.matrix[i,] = accvec
    rm(cls.t,prtvec,accvec)
  }
  
  prt.p = make_summary_plots(dat = prt.matrix,nm = nm,ylab = "Cluster\nPurity")
  acc.p = make_summary_plots(dat = acc.matrix,nm = nm,ylab = "Cluster\nAccuracy")
}

##### Now, combine discrete partition results into one plot
if (TRUE)
{
  library(ggpubr)
  lst = list(adjr.p[[2]],mi.p[[2]],prt.p[[2]],acc.p[[2]])
  lgd = as_ggplot(get_legend(lst[[1]]))
  lst = lapply(lst,function(x) x + guides(color = "none"))
  for (i in 1:(length(lst)-1)) lst[[i]] = lst[[i]] + theme(axis.text.x = element_blank(),axis.title.x = element_blank())
  len = rep(1,length(lst));len[length(lst)] = 1.6;len = len/sum(len)
  p.metrics = plot_grid(plotlist = c(lst,list(lgd)),ncol = 1,align = "v",rel_heights = c(len,len[1]*0.6))
  
  print(p.metrics)
  
  ggsave(filename="Discrete_partition.results.pdf", plot=p.metrics, device=cairo_pdf,
         width = 8, height = 10, units = "in", dpi = 600)
}
  
```

## [Section 5. Comparative evaluation of iterative top-down clustering results]{#S5}
```{r evaluation, echo=FALSE}
setwd(root.dir)
##### evaluate coverage and inclusion rate
if (TRUE)
{
  library(dplyr)
  mod.bench.lst = lapply(clus.data.list,function(x) lapply(x[,3:ncol(x)],function(y) split(x[[1]],factor(y))))
  
  # add multiscale results
  for (i in 1:length(mod.bench.lst))
  {
    mod.bench.lst[[i]] = c(mod.bench.lst[[i]],list(MSC.Euc.ML = lena.euclid.results[[i]]$modules),list(MSC.Cor.ML = lena.cor.results[[i]]$modules))
  }

  incl.vec = covr.vec = matrix(NA,ncol = length(clus.data.list),nrow = length(mod.bench.lst[[i]]));
  colnames(incl.vec) = colnames(covr.vec) = names(clus.data.list)
  detection.lst = vector("list",length(clus.data.list)-1);names(detection.lst) = names(clus.data.list)[-2]
  detect.pdat = data.frame()
  for (i in 1:length(mod.bench.lst))
  {
    #mod.true = split(clus.data.list[[i]][[1]],clus.data.list[[i]][[2]]);
    mod.true = split(colnames(seu.list[[i]]),seu.list[[i]]@meta.data$phenoid)
    mod.calc = mod.bench.lst[[i]]
    
    if (i == 1) rownames(incl.vec) = rownames(covr.vec) = names(mod.calc)
    incl.vec[,i] = sapply(mod.calc,function(x) inclusion_rate(mod.o = mod.true,mod.f = x))
    covr.vec[,i] = sapply(mod.calc,function(x) coverage_rate(mod.o = mod.true,mod.f = x))
    
    detect.m = do.call('cbind',lapply(mod.calc,function(x) apply(calculate_jaccard(mod.true,mod.f = x),1,max)))
    detection.lst[[i]] = detect.m
    len = sapply(mod.true,length)
    
    pdat = do.call('rbind.data.frame',lapply(1:ncol(detect.m),function(x) data.frame(method = rep(colnames(detect.m)[x],nrow(detect.m)),detection.rate = detect.m[,x],nrep = rep(i,nrow(detect.m)),cluster.id = rownames(detect.m))))
    pdat$cluster.size = len[match(pdat$cluster.id,names(len))]
    pdat$data.id = rep(names(detection.lst)[i],nrow(pdat))
    detect.pdat = rbind.data.frame(detect.pdat,pdat)
    rm(mod.true,mod.lena,mod.bench,mod.calc,pdat)
  }
  
  library(reshape)
  pdat.incl = melt(incl.vec);colnames(pdat.incl) = c("method","data.id","inclusion.rate")
  pdat.incl$method = gsub("CIDER","CIDR",gsub("LENA","MSC",gsub("SC3_clusters","SC3",gsub("seurat_snn\\.","SNN, \u03B3=",pdat.incl$method))))
  pdat.covr = melt(covr.vec);colnames(pdat.covr) = c("method","data.id","coverage.rate")
  pdat.covr$method = gsub("CIDER","CIDR",gsub("LENA","MSC",gsub("SC3_clusters","SC3",gsub("seurat_snn\\.","SNN, \u03B3=",pdat.covr$method))))
  
  incl.obj = ggplot(data = pdat.incl) +
    geom_quasirandom(aes(colour = factor(data.id),x = method,y = inclusion.rate),dodge.width = 0.75,alpha = 0.6) + 
    geom_boxplot(aes(x = method,y = inclusion.rate),alpha = 0.2,outlier.colour = NA,width = 0.25) + 
    #stat_summary(aes(x = factor(data.id),group = method,colour = method,y = detection.rate),geom = "line",fun = function(x) mean(x)) + 
    #stat_summary(aes(x = factor(data.id),colour = method,y = detection.rate),geom = "errorbar",fun.data = mean_se,width = 0.2,dodge.width = 0.75) + 
    labs(y = "Inclusion Rate") + 
    scale_x_discrete(labels = nm[unique(pdat.incl$method)]) + 
    guides(colour = "none",fill = guide_legend(title = "Data ID",ncol = 3)) +
    theme_bw() + 
    theme(axis.text.x = element_text(size = 14,angle=45,vjust = 1,hjust = 1),
          axis.text.y = element_text(size = 14),
          axis.title.y = element_text(size = 17),
          axis.title.x = element_blank(),
          plot.title = element_text(size = 20,hjust = 0.5),
          plot.subtitle = element_text(size = 27,face = "bold",hjust = -0.05),
          legend.title = element_text(size = 17),legend.text = element_text(size = 13),
          legend.position = "bottom")
  covr.obj = ggplot(data = pdat.covr) +
    geom_quasirandom(aes(colour = factor(data.id),x = method,y = coverage.rate),dodge.width = 0.75,alpha = 0.6) + 
    geom_boxplot(aes(x = method,y = coverage.rate),alpha = 0.2,outlier.colour = NA,width = 0.25) + 
    #stat_summary(aes(x = factor(data.id),group = method,colour = method,y = detection.rate),geom = "line",fun = function(x) mean(x)) + 
    #stat_summary(aes(x = factor(data.id),colour = method,y = detection.rate),geom = "errorbar",fun.data = mean_se,width = 0.2,dodge.width = 0.75) + 
    labs(y = "Coverage Rate") + 
    scale_x_discrete(labels = nm) + 
    guides(colour = "none",fill = guide_legend(title = "Data ID",ncol = 3)) +
    theme_bw() + 
    theme(axis.text.x = element_text(size = 14,angle=45,vjust = 1,hjust = 1),
          axis.text.y = element_text(size = 14),
          axis.title.y = element_text(size = 17),
          axis.title.x = element_blank(),
          plot.title = element_text(size = 20,hjust = 0.5),
          plot.subtitle = element_text(size = 27,face = "bold",hjust = -0.05),
          legend.title = element_text(size = 17),legend.text = element_text(size = 13),
          legend.position = "bottom")
  
  sum.data = detect.pdat %>% group_by(method,data.id) %>% summarize(overall.detection.rate = sum(detection.rate*cluster.size)/sum(cluster.size))
  acr.obj = ggplot(data = sum.data) + geom_quasirandom(aes(x = method,colour= data.id,y = overall.detection.rate),dodge.width = 0.35,alpha = 0.6) + 
    geom_boxplot(aes(x = method,y = overall.detection.rate),alpha = 0.2,outlier.colour = NA,width =0.2) + 
    labs(y = "Detection Accuracy") + 
    scale_x_discrete(labels = nm) + 
    guides(colour = guide_legend(title = "Data ID",ncol = 3)) +
    theme_bw() + 
    theme(axis.text.x = element_text(size = 14,angle=45,vjust = 1,hjust = 1),
          axis.text.y = element_text(size = 14),
          axis.title.y = element_text(size = 17),
          axis.title.x = element_blank(),
          plot.title = element_text(size = 20,hjust = 0.5),
          plot.subtitle = element_text(size = 27,face = "bold",hjust = -0.05),
          legend.title = element_text(size = 17),legend.text = element_text(size = 13),
          legend.position = "bottom")
}
## combine plots
library(patchwork)
incl.obj = incl.obj + theme(axis.text.x = element_blank())
covr.obj = covr.obj + theme(axis.text.x = element_blank())

metric.obj = incl.obj/covr.obj/acr.obj

print(metric.obj)


```
Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
