---
title: "Cross-platform  evaluation"
author: "Won-Min Song"
date: '2025-07-17'
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 Ding et al. 2020 data set. These data sets have PBMC cell types and subtypes as the silver-standard ground-truth clusters. It requires to download the data from "Ding_et_al_2020_crossplatform" folder under the Synapse project (SynID: syn52966803), with DOI: [https://doi.org/10.7303/syn52966803](https://doi.org/10.7303/syn52966803). Once all data are downloaded,;\
- Make sure to download 'cross_platform_functions.R' to the working directory from the synapse folder to run the codes.;\
- Make sure to unzip 'results_data.zip' under "individual_data" subfolder in the working directory.;\

[Section 1. UMAP visualization](#S1): Seurat object to house the integrated data from Ding et al. 2020 is loaded (cross_platform.seurat.RDS), and respective UMAP plots are plotted.;\
[Section 2. Loading and evaluating clustering results](#S2): Clustering results for each platform is stored under "individual_data" subfolder. These data will be loaded for clustering results per platform. In this section, we only make use of pbmc1 and discard pbmc2 data sets among the two biological replicates studied in Ding et al. 2020, due to scarce cells and smaller pool of platforms used for pbmc2. Also, we load cross_platform.seurat_spit.RDS, which is a list of seurat object and annotation table for individual data.;\

## [Section 1. UMAP Visualization]{#S1}

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

library(MASS)
library(Matrix)
library(ggbeeswarm)
library(MSC)
library(doParallel)
library(igraph)
library(GGally)
library(SC3)
library(aKNNO)
library(scCAN)
library(RaceID)
library(Seurat)
library(SingleCellExperiment)
library(cidr)
library(R.utils)
library(mclust)
library(cowplot)
library(SC3)
library(SingleCellExperiment)

#### Set up folders
root.dir <- "/media/won-min/My Passport1/SingleCell/MSC/PBMC_crossplatform/Results";setwd(root.dir)
out.dir <- getwd()
source("cross_platform_functions.R")


#### load data
seu = readRDS(file = "cross_platform.seurat.RDS") # integrated single-cell data across platforms and replicates per Ding et al. 2020.

# update cell type labels
seu@meta.data$CellType[seu@meta.data$CellType == "Unassigned"] = NA
seu@meta.data$CellType = gsub("monocyte","Monocyte",gsub("T cell","T-cell",gsub("B cell","B-cell",seu@meta.data$CellType)))
seu@meta.data$CellType[seu@meta.data$CellType == "Plasmacytoid dendritic cell"] = "pDC"
seu@meta.data$CellType[seu@meta.data$CellType == "Dendritic cell"] = "mDC"
seu@meta.data$CellType[seu@meta.data$CellType == "Natural killer cell"] = "NK cell"
major.ct = list("B-cell" = c("B-cell"),"T-cell" = c("CD4+ T-cell","Cytotoxic T-cell"),Monocyte = c("CD14+ Monocyte","CD16+ Monocyte"),"Megakaryocyte" = "Megakaryocyte","pDC" = "pDC","mDC" = "mDC","NK cell" = "NK cell")

# update major cell type
vec = rep(NA,ncol(seu))
for (n in 1:length(major.ct)) vec[seu@meta.data$CellType %in% major.ct[[n]]] = names(major.ct)[n]
seu@meta.data$Major.CellType = vec

# create dotplot for numbers of cells used per setting
pdat = reshape::melt(table(seu@meta.data[,c("Experiment","Method")]))
dobj = ggplot(data = pdat,aes(x = Experiment,y = Method,size = value)) + geom_point() + 
  guides(size = guide_legend(title = "#. cells",ncol = 2)) + scale_size_continuous(range = c(0,4)) + 
  theme_classic() + 
  theme(axis.title = element_blank(),axis.text.x = element_text(angle = 45,vjust = 1,hjust = 1,size = 14),axis.text.y = element_text(size = 14),
        legend.position = "bottom",legend.title = element_text(size = 13),legend.text = element_text(size = 12))

# Create UMAP plots 
colmap = c("Megakaryocyte" = "aquamarine","B-cell" = "bisque1","pDC" = "brown1","mDC" = "darkolivegreen1","Monocyte" = "deeppink2","NK cell" = "cyan1","T-cell" = "slateblue",
           "CD4+ T-cell" = "cornflowerblue","Cytotoxic T-cell" = "darkorchid1","CD14+ Monocyte" = "darkorange","CD16+ Monocyte" = "darkred")
umap.mct = DimPlot(seu,group.by = "Major.CellType",cols = colmap,label = T,label.box = T,repel= T,label.color= "black") + guides(colour = "none") + 
  theme(plot.title = element_blank(),axis.line = element_blank(),axis.ticks = element_blank(),axis.text = element_blank(),axis.title = element_blank())
umap.ct = DimPlot(seu,group.by = "CellType",cols = colmap,label = T,label.box = T,repel= T,label.color= "black") + guides(colour = "none") + 
  theme(plot.title = element_blank(),axis.line = element_blank(),axis.ticks = element_blank(),axis.text = element_blank(),axis.title = element_blank())

umap.obj = umap.mct + umap.ct

print(umap.obj)

```

## [Section 2. Load and evaluate individual data]{#S2}

```{r cluster analysis, message = FALSE}
##### Create functions for generating the benchmark clusters
setwd(root.dir)
resfiles = list.files(path = "data_for_synapse",pattern = "clustering_results.v3.RDS",full.names = T,recursive = T)
names(resfiles) = sapply(strsplit(resfiles,"/"),function(x) x[2])

# load annotated individual data
output = readRDS("cross_platform.seurat_split.RDS")
output$data.annotation$datanames = gsub(":","__",gsub("_$","",gsub("__","_",gsub(" |\\(|\\)|-","_",output$data.annotation[[1]]))))
output$data.annotation$cell.num = sapply(output$data.list,ncol)
  
#### Evaluate discrete partitions
if (TRUE)
{
  ii = match(names(resfiles),output$data.annotation$datanames)
 
  rm.row= c(paste0("aKNN_O_res.",c(0.6,1,1.4,1.6,1.8,2)),paste0("seurat_snn_res.",c(0.6,1,1.4,1.6,1.8,2)))
  for (i in 1:length(resfiles))
  {
    
    results = readRDS(resfiles[i])
    results$benchmark.result = results$benchmark.result[!(names(results$benchmark.result) %in% rm.row)]
    seu.obj = output$data.list[[ii[i]]]
    
    seu.obj@meta.data$CellType[seu.obj@meta.data$CellType == "Unassigned"] = NA
    seu.obj@meta.data$CellType = gsub("monocyte","Monocyte",gsub("T cell","T-cell",gsub("B cell","B-cell",seu.obj@meta.data$CellType)))
    seu.obj@meta.data$CellType[seu.obj@meta.data$CellType == "Plasmacytoid dendritic cell"] = "pDC"
    seu.obj@meta.data$CellType[seu.obj@meta.data$CellType == "Dendritic cell"] = "mDC"
    seu.obj@meta.data$CellType[seu.obj@meta.data$CellType == "Natural killer cell"] = "NK cell"
    
    # update major cell type
    vec = rep(NA,ncol(seu.obj))
    for (n in 1:length(major.ct)) vec[seu.obj@meta.data$CellType %in% major.ct[[n]]] = names(major.ct)[n]
    seu.obj@meta.data$Major.CellType = vec

    ###### Coerce discrete partitions
    if (TRUE)
    {
      # now, get discrete partitions from benchmark results
      benchmark.mem = do.call('cbind',lapply(results$benchmark.result,function(x,bg) mods_to_vec(x,bg),bg = colnames(seu.obj)))
      
      ## now, get first splits from MSC
      # Euclidean
      vec.euc = rep(1,ncol(seu.obj))
      split.euc = extract_first_split(msc.obj = results$msc.result$euclidean,use.valid.gate = F,bg = colnames(seu.obj))
      if (!is.null(split.euc)) vec.euc = split.euc[[2]]
      
      # correlation
      vec.cor = rep(1,ncol(seu.obj))
      split.cor = extract_first_split(msc.obj = results$msc.result$correlation,use.valid.gate = F,bg = colnames(seu.obj))
      if (!is.null(split.cor)) vec.cor = split.cor[[2]]
      
      msc.mem = cbind(vec.euc,vec.cor)
      colnames(msc.mem) = c("MSC.EUC.1L","MSC.COR.1L")
      
      # total 
      all.mem = cbind(benchmark.mem,msc.mem)
    }
    
    if (i == 1) 
    {
      adjr.matrix.ct <- adjr.matrix.mct <- 
        nmi.matrix.ct <- nmi.matrix.mct <- 
        prt.matrix.ct <- prt.matrix.mct  <- 
        acc.matrix.ct <- acc.matrix.mct <- 
        matrix(0,nrow = ncol(all.mem),ncol = length(resfiles))
      
      incl.matrix.ct <- incl.matrix.mct <- 
        covr.matrix.ct <- covr.matrix.mct <- 
        det.matrix.ct <- det.matrix.mct <- 
        matrix(0,nrow = ncol(all.mem) + 2,ncol = length(resfiles))
      
      colnames(adjr.matrix.ct) <- colnames(adjr.matrix.mct) <- 
        colnames(nmi.matrix.ct) <- colnames(nmi.matrix.mct) <- 
        colnames(prt.matrix.ct) <-  colnames(prt.matrix.mct) <- 
        colnames(acc.matrix.ct) <-  colnames(acc.matrix.mct) <-
        colnames(incl.matrix.ct) <- colnames(incl.matrix.mct) <- 
        colnames(covr.matrix.ct) <- colnames(covr.matrix.mct) <- 
        colnames(det.matrix.ct) <- colnames(det.matrix.mct) <- 
        output$data.annotation$data.name[match(names(resfiles),output$data.annotation$datanames)]
      
      rownames(adjr.matrix.ct) <- rownames(adjr.matrix.mct) <- 
        rownames(nmi.matrix.ct) <- rownames(nmi.matrix.mct) <- 
        rownames(prt.matrix.ct) <- rownames(prt.matrix.mct) <- 
        rownames(acc.matrix.ct) <- rownames(acc.matrix.mct) <- 
        colnames(all.mem)
      
      rownames(incl.matrix.ct) <- rownames(incl.matrix.mct) <- 
        rownames(covr.matrix.ct) <- rownames(covr.matrix.mct) <- 
        rownames(det.matrix.ct) <- rownames(det.matrix.mct) <- 
        c(colnames(all.mem),c("MSC.EUC.ML","MSC.COR.ML"))
    }
    
    ###### Discrete partition measures
    if (TRUE)
    {
      ### Calculate ARI
      adjr.matrix.ct[,i] = as.vector(apply(all.mem,2,function(x,cls.true) adjustedRandIndex(x = cls.true,y = x),cls.true = seu.obj@meta.data$CellType))
      adjr.matrix.mct[,i] = as.vector(apply(all.mem,2,function(x,cls.true) adjustedRandIndex(x = cls.true,y = x),cls.true = seu.obj@meta.data$Major.CellType))
      
      ### Caculate NMI
      nmi.matrix.ct[,i] = as.vector(apply(all.mem,2,function(x,cls.true) NMI(c = cls.true,t = x),cls.true = seu.obj@meta.data$CellType))
      nmi.matrix.mct[,i] = as.vector(apply(all.mem,2,function(x,cls.true) NMI(c = cls.true,t = x),cls.true = seu.obj@meta.data$Major.CellType))
      
      ### Calculate entropy measures
      prt.matrix.ct[,i] = apply(all.mem,2,function(x,y) get_cluster_purity(x,y)$overall.purity,y = seu.obj@meta.data$CellType)
      acc.matrix.ct[,i] = apply(all.mem,2,function(x,y) get_cluster_accuracy(x,y)$overall.accuracy,y = seu.obj@meta.data$CellType)
      
      prt.matrix.mct[,i] = apply(all.mem,2,function(x,y) get_cluster_purity(x,y)$overall.purity,y = seu.obj@meta.data$Major.CellType)
      acc.matrix.mct[,i] = apply(all.mem,2,function(x,y) get_cluster_accuracy(x,y)$overall.accuracy,y = seu.obj@meta.data$Major.CellType)
    }
    
    ###### Coerce non-discrete measures
    if (TRUE)
    {
      mod.lst = apply(all.mem,2,function(x) split(names(x),factor(x)))
      mod.calc = c(mod.lst,list(MSC.EUC.ML = results$msc.result$euclidean$modules,MSC.COR.ML = results$msc.result$correlation$modules))
      mod.true.ct = split(colnames(seu.obj),seu.obj@meta.data$CellType)
      mod.true.mct = split(colnames(seu.obj),seu.obj@meta.data$Major.CellType)
      
      incl.matrix.ct[,i] = sapply(mod.calc,function(x,cls.true) inclusion_rate(mod.o = cls.true,mod.f = x),cls.true = mod.true.ct)
      covr.matrix.ct[,i] = sapply(mod.calc,function(x,cls.true) coverage_rate(mod.o = cls.true,mod.f = x),cls.true = mod.true.ct)
      incl.matrix.mct[,i] = sapply(mod.calc,function(x,cls.true) inclusion_rate(mod.o = cls.true,mod.f = x),cls.true = mod.true.mct)
      covr.matrix.mct[,i] = sapply(mod.calc,function(x,cls.true) coverage_rate(mod.o = cls.true,mod.f = x),cls.true = mod.true.mct)
      
      detect.m = do.call('cbind',lapply(mod.calc,function(x,cls.true) apply(calculate_jaccard(cls.true,mod.f = x),1,max),cls.true = mod.true.ct))
      len = sapply(mod.true.ct,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(resfiles)[i],nrow(pdat))
      sum.data = pdat %>% group_by(method,data.id) %>% summarize(overall.detection.rate = sum(detection.rate*cluster.size)/sum(cluster.size))
      det.matrix.ct[,i] = sum.data$overall.detection.rate[match(rownames(det.matrix.ct),sum.data$method)]
      
      detect.m = do.call('cbind',lapply(mod.calc,function(x,cls.true) apply(calculate_jaccard(cls.true,mod.f = x),1,max),cls.true = mod.true.mct))
      len = sapply(mod.true.mct,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(resfiles)[i],nrow(pdat))
      sum.data = pdat %>% group_by(method,data.id) %>% summarize(overall.detection.rate = sum(detection.rate*cluster.size)/sum(cluster.size))
      det.matrix.mct[,i] = sum.data$overall.detection.rate[match(rownames(det.matrix.ct),sum.data$method)]
    }
    
  }
  
  # set name labels 
  nm = rownames(det.matrix.ct)
  nm = gsub("aKNN_O_res\\.","aKNNO, \u03B3=",gsub("seurat_snn_res\\.","SNN, \u03B3=",nm))
  nm = as.list(nm);names(nm) = rownames(adjr.matrix.ct)
  nm['MSC.EUC.1L'][[1]] = bquote(MSC * scriptstyle(atop("EUC", "1L")))
  nm['MSC.COR.1L'][[1]] = bquote(MSC * scriptstyle(atop("COR", "1L")))
  nm['MSC.EUC.ML'][[1]] = bquote(MSC * scriptstyle(atop("EUC", "ML")))
  nm['MSC.COR.ML'][[1]] = bquote(MSC * scriptstyle(atop("COR", "ML")))
  
  ###### create plots for discrete partitions
  # ARI
  if (TRUE)
  {
    pdat1 = reshape::melt(adjr.matrix.ct);pdat1$cell.type.level = "Cell subtypes"
    pdat2 = reshape::melt(adjr.matrix.mct);pdat2$cell.type.level = "Major cell types"
    pdata = rbind.data.frame(pdat1,pdat2)
    colnames(pdata)[1:3] = c("method","data.id","mi")
    ii = match(pdata$data.id,output$data.annotation$data.name)
    pdata$sequencing.method = output$data.annotation$method[ii]
    pdata$experiment = output$data.annotation$experiment[ii]
    
    # subset for pbmc1 
    pdata = subset(pdata,experiment == "pbmc1")
    
    pobj.ARI = ggplot() + geom_quasirandom(data = pdata,aes(x = method,y = mi,colour = sequencing.method)) +
      labs(x = "Methods",y = "ARI") + 
      guides(colour = guide_legend(title = "Data ID",ncol = 4)) + 
      scale_x_discrete(labels = nm) + 
      facet_grid(. ~ cell.type.level,scale = "free",space = "free") + 
      geom_boxplot(data = pdata,aes(x = method,y = mi),width = 0.6,alpha = 0.2,outlier.colour = NA) + theme_bw() + 
      theme(axis.text.x = element_text(size = 16,angle = 45,vjust = 1,hjust = 1),axis.text.y = element_text(size = 15),
            axis.title = element_text(size = 20),legend.title = element_text(size = 19),legend.text = element_text(size = 15),
            plot.subtitle = element_text(size = 27,face = "bold",hjust = -0.05),
            strip.text = element_text(size = 17),
            legend.position = "bottom")
  }
  
  # NMI
  if (TRUE)
  {
    pdat1 = reshape::melt(nmi.matrix.ct);pdat1$cell.type.level = "Cell subtypes"
    pdat2 = reshape::melt(nmi.matrix.mct);pdat2$cell.type.level = "Major cell types"
    pdata = rbind.data.frame(pdat1,pdat2)
    colnames(pdata)[1:3] = c("method","data.id","mi")
    ii = match(pdata$data.id,output$data.annotation$data.name)
    pdata$sequencing.method = output$data.annotation$method[ii]
    pdata$experiment = output$data.annotation$experiment[ii]
    
    # subset for pbmc1 
    pdata = subset(pdata,experiment == "pbmc1")
    
    pobj.NMI = ggplot() + geom_quasirandom(data = pdata,aes(x = method,y = mi,colour = sequencing.method)) +
      labs(x = "Methods",y = "NMI") + 
      guides(colour = guide_legend(title = "Data ID",ncol = 4)) + 
      scale_x_discrete(labels = nm) + 
      facet_grid(. ~ cell.type.level) + 
      geom_boxplot(data = pdata,aes(x = method,y = mi),width = 0.6,alpha = 0.2,outlier.colour = NA) + theme_bw() + 
      theme(axis.text.x = element_text(size = 16,angle = 45,vjust = 1,hjust = 1),axis.text.y = element_text(size = 15),
            axis.title = element_text(size = 20),legend.title = element_text(size = 19),legend.text = element_text(size = 15),
            plot.subtitle = element_text(size = 27,face = "bold",hjust = -0.05),
            strip.text = element_text(size = 17),
            legend.position = "bottom")
  }
  
  ### Cluster Purity
  if (TRUE)
  {
    pdat1 = reshape::melt(prt.matrix.ct);pdat1$cell.type.level = "Cell subtypes"
    pdat2 = reshape::melt(prt.matrix.mct);pdat2$cell.type.level = "Major cell types"
    pdata = rbind.data.frame(pdat1,pdat2)
    colnames(pdata)[1:3] = c("method","data.id","mi")
    ii = match(pdata$data.id,output$data.annotation$data.name)
    pdata$sequencing.method = output$data.annotation$method[ii]
    pdata$experiment = output$data.annotation$experiment[ii]
    # subset for pbmc1 
    pdata = subset(pdata,experiment == "pbmc1")
    
    pobj.prt = ggplot() + geom_quasirandom(data = pdata,aes(x = method,y = mi,colour = sequencing.method)) +
      labs(x = "Methods",y = "Cluster Purity") + 
      guides(colour = guide_legend(title = "Data ID",ncol = 4)) + 
      scale_x_discrete(labels = nm) + 
      facet_grid(. ~ cell.type.level) + 
      geom_boxplot(data = pdata,aes(x = method,y = mi),width = 0.6,alpha = 0.2,outlier.colour = NA) + theme_bw() + 
      theme(axis.text.x = element_text(size = 16,angle = 45,vjust = 1,hjust = 1),axis.text.y = element_text(size = 15),
            axis.title = element_text(size = 20),legend.title = element_text(size = 19),legend.text = element_text(size = 15),
            plot.subtitle = element_text(size = 27,face = "bold",hjust = -0.05),
            strip.text = element_text(size = 17),
            legend.position = "bottom")
  }
  
  ### Cluster Accuracy
  if (TRUE)
  {
    pdat1 = reshape::melt(acc.matrix.ct);pdat1$cell.type.level = "Cell subtypes"
    pdat2 = reshape::melt(acc.matrix.mct);pdat2$cell.type.level = "Major cell types"
    pdata = rbind.data.frame(pdat1,pdat2)
    colnames(pdata)[1:3] = c("method","data.id","mi")
    ii = match(pdata$data.id,output$data.annotation$data.name)
    pdata$sequencing.method = output$data.annotation$method[ii]
    pdata$experiment = output$data.annotation$experiment[ii]
    # subset for pbmc1 
    pdata = subset(pdata,experiment == "pbmc1")
    
    pobj.acc = ggplot() + geom_quasirandom(data = pdata,aes(x = method,y = mi,colour = sequencing.method)) +
      labs(x = "Methods",y = "Cluster Accuracy") + 
      guides(colour = guide_legend(title = "Data ID",ncol = 4)) + 
      scale_x_discrete(labels = nm) + 
      facet_grid(. ~ cell.type.level) + 
      geom_boxplot(data = pdata,aes(x = method,y = mi),width = 0.6,alpha = 0.2,outlier.colour = NA) + theme_bw() + 
      theme(axis.text.x = element_text(size = 16,angle = 45,vjust = 1,hjust = 1),axis.text.y = element_text(size = 15),
            axis.title = element_text(size = 20),legend.title = element_text(size = 19),legend.text = element_text(size = 15),
            plot.subtitle = element_text(size = 27,face = "bold",hjust = -0.05),
            strip.text = element_text(size = 17),
            legend.position = "bottom")
  }
  
  ###### Create plots for non-discrete partitions
  ### Detection Accuracy
  if (TRUE)
  {
    pdat1 = reshape::melt(det.matrix.ct);pdat1$cell.type.level = "Cell subtypes"
    pdat2 = reshape::melt(det.matrix.mct);pdat2$cell.type.level = "Major cell types"
    pdata = rbind.data.frame(pdat1,pdat2)
    colnames(pdata)[1:3] = c("method","data.id","mi")
    ii = match(pdata$data.id,output$data.annotation$data.name)
    pdata$sequencing.method = output$data.annotation$method[ii]
    pdata$experiment = output$data.annotation$experiment[ii]
    
    # subset for pbmc1 
    pdata = subset(pdata,experiment == "pbmc1")
    
    pobj.det = ggplot() + geom_quasirandom(data = pdata,aes(x = method,y = mi,colour = sequencing.method)) +
      labs(x = "Methods",y = "Detection Accuracy") + 
      guides(colour = guide_legend(title = "Data ID",ncol = 4)) + 
      scale_x_discrete(labels = nm) + 
      facet_grid(. ~ cell.type.level) + 
      geom_boxplot(data = pdata,aes(x = method,y = mi),width = 0.6,alpha = 0.2,outlier.colour = NA) + theme_bw() + 
      theme(axis.text.x = element_text(size = 16,angle = 45,vjust = 1,hjust = 1),axis.text.y = element_text(size = 15),
            axis.title = element_text(size = 20),legend.title = element_text(size = 19),legend.text = element_text(size = 15),
            plot.subtitle = element_text(size = 27,face = "bold",hjust = -0.05),
            strip.text = element_text(size = 17),
            legend.position = "bottom")
  }
  
  ### Inclusion
  if (TRUE)
  {
    pdat1 = reshape::melt(incl.matrix.ct);pdat1$cell.type.level = "Cell subtypes"
    pdat2 = reshape::melt(incl.matrix.mct);pdat2$cell.type.level = "Major cell types"
    pdata = rbind.data.frame(pdat1,pdat2)
    colnames(pdata)[1:3] = c("method","data.id","mi")
    ii = match(pdata$data.id,output$data.annotation$data.name)
    pdata$sequencing.method = output$data.annotation$method[ii]
    pdata$experiment = output$data.annotation$experiment[ii]
    # subset for pbmc1 
    pdata = subset(pdata,experiment == "pbmc1")
    
    pobj.incl = ggplot() + geom_quasirandom(data = pdata,aes(x = method,y = mi,colour = sequencing.method)) +
      labs(x = "Methods",y = "Inclusion") + 
      guides(colour = guide_legend(title = "Data ID",ncol = 4)) + 
      scale_x_discrete(labels = nm) + 
      facet_grid(. ~ cell.type.level) + 
      geom_boxplot(data = pdata,aes(x = method,y = mi),width = 0.6,alpha = 0.2,outlier.colour = NA) + theme_bw() + 
      theme(axis.text.x = element_text(size = 16,angle = 45,vjust = 1,hjust = 1),axis.text.y = element_text(size = 15),
            axis.title = element_text(size = 20),legend.title = element_text(size = 19),legend.text = element_text(size = 15),
            plot.subtitle = element_text(size = 27,face = "bold",hjust = -0.05),
            strip.text = element_text(size = 17),
            legend.position = "bottom")
  }
  
  ### Coverage
  if (TRUE)
  {
    pdat1 = reshape::melt(covr.matrix.ct);pdat1$cell.type.level = "Cell subtypes"
    pdat2 = reshape::melt(covr.matrix.mct);pdat2$cell.type.level = "Major cell types"
    pdata = rbind.data.frame(pdat1,pdat2)
    colnames(pdata)[1:3] = c("method","data.id","mi")
    ii = match(pdata$data.id,output$data.annotation$data.name)
    pdata$sequencing.method = output$data.annotation$method[ii]
    pdata$experiment = output$data.annotation$experiment[ii]
    
    # subset for pbmc1 
    pdata = subset(pdata,experiment == "pbmc1")
    
    pobj.covr = ggplot() + geom_quasirandom(data = pdata,aes(x = method,y = mi,colour = sequencing.method)) +
      labs(x = "Methods",y = "Coverage") + 
      guides(colour = guide_legend(title = "Data ID",ncol = 4)) + 
      scale_x_discrete(labels = nm) + 
      facet_grid(. ~ cell.type.level) + 
      geom_boxplot(data = pdata,aes(x = method,y = mi),width = 0.6,alpha = 0.2,outlier.colour = NA) + theme_bw() + 
      theme(axis.text.x = element_text(size = 16,angle = 45,vjust = 1,hjust = 1),axis.text.y = element_text(size = 15),
            axis.title = element_text(size = 20),legend.title = element_text(size = 19),legend.text = element_text(size = 15),
            plot.subtitle = element_text(size = 27,face = "bold",hjust = -0.05),
            strip.text = element_text(size = 17),
            legend.position = "bottom")
  }
}

##### Now, combine discrete partition results into one plot
if (TRUE)
{
  library(cowplot)
  library(ggpubr)
  lst = list(pobj.ARI,pobj.NMI,pobj.prt,pobj.acc)
  lst[1:(length(lst)-1)] = lapply(lst[1:(length(lst)-1)],function(x) x + guides(color = "none",shape = "none",fill = "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)] = 2.3;len = len/sum(len)
  p.metrics.disc = plot_grid(plotlist = c(lst),ncol = 1,align = "v",rel_heights = len)
  print(p.metrics.disc)
  
}


##### Now, combine overlap results into one plot
if (TRUE)
{
  library(cowplot)
  library(ggpubr)
  lst = list(pobj.incl,pobj.covr,pobj.det)
  lst[1:(length(lst)-1)] = lapply(lst[1:(length(lst))],function(x) x + guides(color = "none",shape = "none",fill = "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)] = 2.1;len = len/sum(len)
  p.metrics.ndisc = plot_grid(plotlist = c(lst),ncol = 1,align = "v",rel_heights = len)
  

  print(p.metrics.ndisc)


}

```

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