summary_metric_plot_public.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d842ac76",
   "metadata": {},
   "source": [
    "# Readme\n",
    "\n",
    "This script contains codes for figure generation described in the manuscript titled 'Benchmarking algorithms for joint integration of unpaired and paired single-cell RNA-seq and ATAC-seq data'. \n",
    "\n",
    "Please run steps described in evaluate_vary_situations_public.ipynb before this script. \n",
    "\n",
    "For scenario 1-3, one summary plot is generated for each challenge. For senario 4 and 5, multiple plots could be generated based on the number of missing cell types. \n",
    "\n",
    "Additional plots are generated using umap_plot_generation_public.ipynb and hpap_result_plot_public.ipynb\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3c3f429",
   "metadata": {},
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc74e06b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run with an R kernel\n",
    "library(dplyr)\n",
    "library(ggplot2)\n",
    "library(tidyverse)\n",
    "library(ComplexHeatmap)\n",
    "library(gridExtra)\n",
    "\n",
    "plot_metric <- function(res_df, folder_path,file_name,ct_metric=c('ARI','NMI','ct_aws','clisi'),\n",
    "                        batch_metric=c('b_saw','kbet','ilisi','gconn'),var='var',x_label=\"\"){\n",
    "    pdf(paste0(folder_path,file_name))\n",
    "    if(x_label==\"\"){\n",
    "        x_label=var\n",
    "    }\n",
    "    # ct metrics\n",
    "    print(\"===== cell type annotation metrics =====\")\n",
    "    ct_glist = list()\n",
    "    \n",
    "    \n",
    "    for (col in ct_metric){\n",
    "        ct_glist[[col]] = ggplot(res_df,aes_string(var,col,color=\"method\")) + \n",
    "            geom_point() +\n",
    "            stat_summary(aes(group = method),fun=mean, geom=\"line\") +\n",
    "            xlab(x_label)\n",
    "    }\n",
    "\n",
    "    gridExtra::grid.arrange(grobs=ct_glist)\n",
    "\n",
    "    batch_glist = list()\n",
    "    for (col in batch_metric){\n",
    "        batch_glist[[col]] = ggplot(res_df,aes_string(var,col,color=\"method\")) + \n",
    "            geom_point() +\n",
    "            stat_summary(aes(group = method),fun=mean, geom=\"line\") +\n",
    "            xlab(x_label)\n",
    "    }\n",
    "    print(\"===== batch removal metrics =====\")\n",
    "    gridExtra::grid.arrange(grobs=batch_glist)\n",
    "    dev.off()\n",
    "    return(list(ct_glist,batch_glist))\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "plot_metric_boxplot <- function(res_df, folder_path,file_name,ct_metric=c('ARI','NMI','ct_aws','clisi'),\n",
    "                        batch_metric=c('b_saw','kbet','ilisi','gconn'),var='var',title=\"\"){\n",
    "    pdf(paste0(folder_path,file_name))\n",
    "\n",
    "    # ct metrics\n",
    "    print(\"===== cell type annotation metrics =====\")\n",
    "    ct_glist = list()\n",
    "    \n",
    "    \n",
    "    for (col in ct_metric){\n",
    "        ct_glist[[col]] = ggplot(res_df,aes_string(x=\"method\",y=col,fill=\"method\")) + \n",
    "            geom_boxplot() +\n",
    "            ggtitle(title) +\n",
    "            theme(plot.title = element_text(hjust = 0.5),\n",
    "                 axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))\n",
    "    }\n",
    "\n",
    "    gridExtra::grid.arrange(grobs=ct_glist)\n",
    "\n",
    "    batch_glist = list()\n",
    "    for (col in batch_metric){\n",
    "        batch_glist[[col]] = ggplot(res_df,aes_string(x=\"method\",y=col,fill=\"method\")) + \n",
    "            geom_boxplot() +\n",
    "            ggtitle(title) +\n",
    "            theme(plot.title = element_text(hjust = 0.5),\n",
    "                 axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))\n",
    "    }\n",
    "    print(\"===== batch removal metrics =====\")\n",
    "    gridExtra::grid.arrange(grobs=batch_glist)\n",
    "    dev.off()\n",
    "    return(list(ct_glist,batch_glist))\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "plot_result <- function(dir_path,res_folder_filter=NULL,runtime_exists=T,prediction_exists=T,metric_file=\"metric.csv\"){\n",
    "    file_paths = list.files(dir_path,metric_file,recursive = T)\n",
    "    res_df <- c()\n",
    "    for(i in file_paths){\n",
    "        #print(i)\n",
    "        m_i = read.csv(file.path(dir_path,i))\n",
    "        elements = stringr::str_split(i,\"/\")[[1]]\n",
    "        result_out = elements[2]\n",
    "        method = elements[3]\n",
    "        ct_ref = elements[4]\n",
    "        orig = gsub(\"[0-9]*\",\"\",gsub(\"_.*\",\"\",elements[1]))\n",
    "        var = gsub(\"[a-z]*[A-Z]*\",\"\",gsub(\"_.*\",\"\",elements[1]))\n",
    "        rep = gsub(\".*_\",\"\",elements[1])\n",
    "        res_i = as.data.frame(c(m_i,c(\"orig\"=orig,\"var\"=var,\"rep\"=rep,\"method\"=method,\"out\"=result_out,\"ct_ref\"=ct_ref)))\n",
    "        res_df <- rbind(res_df,res_i)\n",
    "    }\n",
    "    res_df$var <- as.integer(res_df$var)\n",
    "    \n",
    "    # save intemediate result to output variable \n",
    "    res_all <- res_df\n",
    "    if(runtime_exists){\n",
    "        runtimes = list.files(dir_path,\".*runtime.txt\",recursive = T)\n",
    "        res_df2 <- c()\n",
    "        for(i in runtimes){\n",
    "            m_i = read.table(file.path(dir_path,i))\n",
    "            elements = stringr::str_split(i,\"/\")[[1]]\n",
    "            result_out = elements[2]\n",
    "            method = gsub(\"_.*\",\"\",elements[4])\n",
    "            orig = gsub(\"[0-9]*\",\"\",gsub(\"_.*\",\"\",elements[1]))\n",
    "            var = gsub(\"[a-z]*[A-Z]*\",\"\",gsub(\"_.*\",\"\",elements[1]))\n",
    "            rep = gsub(\".*_\",\"\",elements[1])\n",
    "            res_i = as.data.frame(c(m_i,c(\"orig\"=orig,\"var\"=var,\"rep\"=rep,\"method\"=method,\"out\"=result_out)))\n",
    "            res_df2 <- rbind(res_df2,res_i)\n",
    "        }\n",
    "\n",
    "        rownames(res_df2) <- NULL\n",
    "        res_df2<-data.frame(res_df2)\n",
    "        res_df2$var<-as.integer(res_df2$var)\n",
    "        colnames(res_df2)[1]=\"runtime\"\n",
    "\n",
    "        res_df2$runtime<-as.numeric(res_df2$runtime)\n",
    "        \n",
    "        res_all <- left_join(res_df,res_df2)\n",
    "    }\n",
    "    \n",
    "    if(prediction_exists){\n",
    "        prediciton_eval = list.files(dir_path,\".*prediction_eval.txt\",recursive = T)\n",
    "        res_df3 <- c()\n",
    "        for(i in prediciton_eval){\n",
    "            m_i = read.csv(file.path(dir_path,i),header = F,sep='\\t')\n",
    "            vec = m_i[,2]\n",
    "            names(vec) = unlist(m_i[1])\n",
    "            elements = stringr::str_split(i,\"/\")[[1]]\n",
    "            result_out = elements[2]\n",
    "            method = elements[3]\n",
    "            orig = gsub(\"[0-9]*\",\"\",gsub(\"_.*\",\"\",elements[1]))\n",
    "            var = gsub(\"[a-z]*[A-Z]*\",\"\",gsub(\"_.*\",\"\",elements[1]))\n",
    "            rep = gsub(\".*_\",\"\",elements[1])\n",
    "            res_i = t(data.frame(c(vec,\"orig\"=orig,\"var\"=var,\"rep\"=rep,\"method\"=method,\"out\"=result_out)))\n",
    "            res_df3 <- rbind(res_df3,res_i)\n",
    "        }\n",
    "\n",
    "        rownames(res_df3) <- NULL\n",
    "        res_df3<-data.frame(res_df3)\n",
    "        res_df3$var<-as.integer(res_df3$var)\n",
    "        res_df3$percent_recovered_50kb<-as.numeric(res_df3$percent_recovered_50kb)\n",
    "        res_df3$f1<-as.numeric(res_df3$f1)\n",
    "        res_all <- left_join(res_all,res_df3)\n",
    "    }\n",
    "    if(!is.null(res_folder_filter)){\n",
    "        res_all = res_all %>% filter(out %in% res_folder_filter)\n",
    "    }\n",
    "   return(res_all)\n",
    "}\n",
    "\n",
    "\n",
    "# format unpaired integration method result, print out one average value despite changes in multiome cells, because they did not use the multiome cells \n",
    "\n",
    "format_unpaired_res <- function(res, n_times){\n",
    "    output <- do.call(cbind,apply(unpaired_method_perform,2,function(x){\n",
    "        y<-as.data.frame(x)\n",
    "        do.call(cbind,rep(list(y),n_times))})\n",
    "    )\n",
    "    return(output)\n",
    "}\n",
    "\n",
    "# get observed (truth) peak-gene association accuracy result (percentage and F1)\n",
    "## Assuming there is \"truth_prediction_eval.txt\" file under the dir_path folder. Read in results and format the table so it shows which repeat it is. \n",
    "\n",
    "get_pred_truth <- function(dir_path){\n",
    "    prediciton_eval = list.files(dir_path,\"truth_prediction_eval.txt\",recursive = T)\n",
    "    res_df3<-c()\n",
    "    for(i in prediciton_eval){\n",
    "        m_i = read.csv(file.path(dir_path,i),header = F,sep='\\t')\n",
    "        vec = m_i[,2]\n",
    "        names(vec) = unlist(m_i[1])\n",
    "        elements = stringr::str_split(i,\"/\")[[1]]\n",
    "        result_out = elements[2]\n",
    "        method = elements[3]\n",
    "        orig = gsub(\"[0-9]*\",\"\",gsub(\"_.*\",\"\",elements[1]))\n",
    "        var = gsub(\"[a-z]*[A-Z]*\",\"\",gsub(\"_.*\",\"\",elements[1]))\n",
    "        rep = gsub(\".*_\",\"\",elements[1])\n",
    "        res_i = t(data.frame(c(vec,\"orig\"=orig,\"var\"=var,\"rep\"=rep,\"method\"=method,\"out\"=result_out)))\n",
    "        res_df3 <- rbind(res_df3,res_i)\n",
    "    }\n",
    "\n",
    "    rownames(res_df3) <- NULL\n",
    "    res_df3<-data.frame(res_df3)\n",
    "    res_df3$var<-as.integer(res_df3$var)\n",
    "    res_df3$percent_recovered_50kb<-as.numeric(res_df3$percent_recovered_50kb)\n",
    "    res_df3$f1<-as.numeric(res_df3$f1)\n",
    "    return(res_df3)\n",
    "}\n",
    "\n",
    "\n",
    "plot_rare_cell_result <- function(dir_path,res_folder_filter=NULL,runtime_exists=T,prediction_exists=T,metric_file=\"metric.csv\"){\n",
    "    file_paths = list.files(dir_path,metric_file,recursive = T)\n",
    "    res_df <- c()\n",
    "    for(i in file_paths){\n",
    "        #print(i)\n",
    "        m_i = read.table(file.path(dir_path,i),sep=',',header = T)[1,]\n",
    "        elements = stringr::str_split(i,\"/\")[[1]]\n",
    "        result_out = elements[2]\n",
    "        method = elements[3]\n",
    "        ct_ref = elements[4]\n",
    "        elements_sep = stringr::str_split(elements[1],\"_\")[[1]]\n",
    "        orig = elements_sep[1]\n",
    "        var = gsub(\"[a-z]*[A-Z]*\",\"\",gsub(\"_.*\",\"\",elements_sep[2]))\n",
    "        rep = elements_sep[3]\n",
    "        res_i = as.data.frame(c(m_i,c(\"orig\"=orig,\"var\"=var,\"rep\"=rep,\"method\"=method,\"out\"=result_out,\"ct_ref\"=ct_ref,out_path = elements[1])))\n",
    "        res_df <- rbind(res_df,res_i)\n",
    "    }\n",
    "    res_df$var <- as.integer(res_df$var)\n",
    "    \n",
    "    # save intemediate result to output variable \n",
    "    res_all <- res_df\n",
    "    if(runtime_exists){\n",
    "        runtimes = list.files(dir_path,\".*runtime.txt\",recursive = T)\n",
    "        res_df2 <- c()\n",
    "        for(i in runtimes){\n",
    "            m_i = read.table(file.path(dir_path,i))\n",
    "            elements = stringr::str_split(i,\"/\")[[1]]\n",
    "            method = method = gsub(\"_.*\",\"\",elements[4])\n",
    "            result_out = elements[2]\n",
    "            res_i = as.data.frame(c(m_i,c(\"method\"=method,\"out\"=result_out,out_path = elements[1])))\n",
    "            res_df2 <- rbind(res_df2,res_i)\n",
    "        }\n",
    "\n",
    "        rownames(res_df2) <- NULL\n",
    "        res_df2<-data.frame(res_df2)\n",
    "        colnames(res_df2)[1]=\"runtime\"\n",
    "\n",
    "        res_df2$runtime<-as.numeric(res_df2$runtime)\n",
    "        \n",
    "        res_all <- left_join(res_df,res_df2,by = c(\"method\",\"out\",\"out_path\"))\n",
    "    }\n",
    "    \n",
    "    if(prediction_exists){\n",
    "        prediciton_eval = list.files(dir_path,\".*prediction_eval.txt\",recursive = T)\n",
    "        res_df3 <- c()\n",
    "        for(i in prediciton_eval){\n",
    "            m_i = read.csv(file.path(dir_path,i),header = F,sep='\\t')\n",
    "            vec = m_i[,2]\n",
    "            names(vec) = unlist(m_i[1])\n",
    "            elements = stringr::str_split(i,\"/\")[[1]]\n",
    "            result_out = elements[2]\n",
    "            method = elements[3]\n",
    "            orig = gsub(\"[0-9]*\",\"\",gsub(\"_.*\",\"\",elements[1]))\n",
    "            var = gsub(\"[a-z]*[A-Z]*\",\"\",gsub(\"_.*\",\"\",elements[1]))\n",
    "            rep = gsub(\".*_\",\"\",elements[1])\n",
    "            res_i = t(data.frame(c(vec,\"orig\"=orig,\"var\"=var,\"rep\"=rep,\"method\"=method,\"out\"=result_out)))\n",
    "            res_df3 <- rbind(res_df3,res_i)\n",
    "        }\n",
    "\n",
    "        rownames(res_df3) <- NULL\n",
    "        res_df3<-data.frame(res_df3)\n",
    "        res_df3$var<-as.integer(res_df3$var)\n",
    "        res_df3$percent_recovered_50kb<-as.numeric(res_df3$percent_recovered_50kb)\n",
    "        res_df3$f1<-as.numeric(res_df3$f1)\n",
    "        res_all <- left_join(res_all,res_df3)\n",
    "    }\n",
    "    if(!is.null(res_folder_filter)){\n",
    "        res_all = res_all %>% filter(out %in% res_folder_filter)\n",
    "    }\n",
    "   return(res_all)\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94d900c1",
   "metadata": {},
   "source": [
    "# Color pallete"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f80518ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\",\"scMoMaT\",\"Seurat v4 integrate\")\n",
    "hex <- c(\"#cc6677\",\"#999933\",\"#44aa99\",\"#117733\",\"#88ccee\",\"#882255\",\"#332288\",\"#ddcc77\",\"#aa4499\",\"#555555\")\n",
    "names(hex)<-method_names\n",
    "hex"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a485ea30",
   "metadata": {},
   "source": [
    "# Scenario 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db8fed0e",
   "metadata": {},
   "source": [
    "## PBMC"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f374db6d",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e463e9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = \"dataset/multiome_pbmc_10k/pbmc_vary_cell_test/\"\n",
    "\n",
    "res_all<- plot_result(dir_path,res_folder_filter = c(\"results_single_mod\",\"results_single_same_cell_number\"),runtime_exists=T,prediction_exists=T)\n",
    "\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "unpaired_method_list<-c(\"seurat3\",\"rfigr\",\"rbindsc\",\"rliger\",\"glue\")\n",
    "res_plot <- res_all %>% \n",
    "    filter(var %in% c(1000,3000,8000)) %>% \n",
    "    mutate(method_name = ifelse(out==\"results_single_same_cell_number\",paste0(method,\"_2\"),method)) %>%\n",
    "    # assign method type \n",
    "    mutate(type = ifelse(out!=\"results_single_same_cell_number\",\n",
    "                         ifelse(method%in%unpaired_method_list,\"Unpaired\",\"Multiome-guided\"),\n",
    "                         \"Unpaired_multiome_splitted\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"percent_recovered_50kb\",\"f1\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "res_plot = res_plot %>% filter(ct_ref ==\"ct3_metric.csv\")\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Peak-gene \\n pair % recovered\",\"F1\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n",
    "#res_plot_2$metric_name <- dplyr::recode(res_plot_2$metric_name,\"Gene-peak \\n pair % recovered\" =\"Peak-gene \\n pair % recovered\")\n",
    "\n",
    "# get observed accuracy\n",
    "plot_ref_accuracy = TRUE\n",
    "#dir_path = \"dataset/multiome_pbmc_10k/multiome_ncells_pmat/\"\n",
    "truth_tbl <- get_pred_truth(dir_path)\n",
    "truth_df <- data.frame(metric_name = factor(c(\"Peak-gene \\n pair % recovered\",\"F1\"),\n",
    "                                            levels=c(\"Peak-gene \\n pair % recovered\",\"F1\")),\n",
    "                       Z = as.numeric(round(colMeans(truth_tbl[,1:2]),digits=2)))\n",
    "truth_df\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32efb319",
   "metadata": {},
   "source": [
    "### Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af125258",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"pbmc_single_mod_n_splitted_metric_summary19_7ct_w_observed_rev1.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),c(\"Cell type ASW\",\"cLisi\"),c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Peak-gene \\n pair % recovered\",\"F1\"),c('log2(Runtime)',\"Graph \\n connectivity\"))\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "x_interval = c(1000,3000,8000) \n",
    "\n",
    "ylim_list = list(c(0.2,1),c(0.4,0.9),c(0.5,0.8),c(0.8,1),c(0.6,1),c(0,1),c(0.25,0.5),c(0.25,0.45),c(0,1),c(6,14))\n",
    "nbreaks = list(4,5,6,4,4,5,5,4,5,8)\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=35,height=15)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13ce2514",
   "metadata": {},
   "source": [
    "## BMMC"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d6d2971",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a2740e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path =\"dataset/bmmc/bmmc_vary_cell_test/\"\n",
    "res_all<- plot_result(dir_path,res_folder_filter = c(\"results_single_mod\",\"results_single_same_cell_number\"),runtime_exists=T,prediction_exists=T)\n",
    "\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "unpaired_method_list<-c(\"seurat3\",\"rfigr\",\"rbindsc\",\"rliger\",\"glue\")\n",
    "res_plot <- res_all %>% \n",
    "    mutate(method_name = ifelse(out==\"results_single_same_cell_number\",paste0(method,\"_2\"),method)) %>%\n",
    "    # assign method type \n",
    "    mutate(type = ifelse(out!=\"results_single_same_cell_number\",\n",
    "                         ifelse(method%in%unpaired_method_list,\"Unpaired\",\"Multiome-guided\"),\n",
    "                         \"Unpaired_multiome_splitted\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"percent_recovered_50kb\",\"f1\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "res_plot = res_plot %>% filter(ct_ref ==\"ct3_metric.csv\")\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Peak-gene \\n pair % recovered\",\"F1\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n",
    "res_plot_2$metric_name <- dplyr::recode(res_plot_2$metric_name,\"Gene-peak \\n pair % recovered\" =\"Peak-gene \\n pair % recovered\")\n",
    "\n",
    "\n",
    "# get observed accuracy\n",
    "plot_ref_accuracy = TRUE\n",
    "#dir_path = \"dataset/multiome_pbmc_10k/multiome_ncells_pmat/\"\n",
    "truth_tbl <- get_pred_truth(dir_path)\n",
    "truth_df <- data.frame(metric_name = factor(c(\"Peak-gene \\n pair % recovered\",\"F1\"),\n",
    "                                            levels=c(\"Peak-gene \\n pair % recovered\",\"F1\")),\n",
    "                       Z = as.numeric(round(colMeans(truth_tbl[,1:2]),digits=2)))\n",
    "truth_df\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc8c1745",
   "metadata": {},
   "source": [
    "### Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ce2e9d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"bmmc_single_mod_n_splitted_metric_summary19_21ct_w_observed_rev1.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),c(\"Cell type ASW\",\"cLisi\"),c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Peak-gene \\n pair % recovered\",\"F1\"),c('log2(Runtime)',\"Graph \\n connectivity\"))\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "x_interval = c(1000,2000,4000) # unique(res_plot_3$var)\n",
    "ylim_list = list(c(0,0.5),c(0.2,0.7),c(0.45,0.7),c(0.85,1),c(0.4,1),c(0,1),c(0.05,0.3),c(0.05,0.3),c(0,1),c(6,14))\n",
    "nbreaks = list(5,5,5,6,6,5,5,5,5,8)\n",
    "counter=1\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                           \n",
    "pdf(file.path(folder_path,file_name),width=35,height=15)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d85574d",
   "metadata": {},
   "source": [
    "## SHARE-seq "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01463b1e",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da23b8cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = \"dataset/mouse_skin/multiome_ncells_pmat/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "res_all<- plot_result(dir_path,res_folder_filter = c(\"results_single_mod\",\"results_single_same_cell_number\"),runtime_exists=T,prediction_exists=T)\n",
    "table(res_all$method)\n",
    "table(res_all$out)\n",
    "\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "unpaired_method_list<-c(\"seurat3\",\"rfigr\",\"rbindsc\",\"rliger\",\"glue\")\n",
    "res_plot <- res_all %>% \n",
    "    mutate(method_name = ifelse(out==\"results_single_same_cell_number\",paste0(method,\"_2\"),method)) %>%\n",
    "    # assign method type \n",
    "    mutate(type = ifelse(out!=\"results_single_same_cell_number\",\n",
    "                         ifelse(method%in%unpaired_method_list,\"Unpaired\",\"Multiome-guided\"),\n",
    "                         \"Unpaired_multiome_splitted\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"percent_recovered_50kb\",\"f1\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "res_plot = res_plot %>% filter(ct_ref ==\"ct3_metric.csv\")\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Peak-gene \\n pair % recovered\",\"F1\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n",
    "res_plot_2$metric_name <- dplyr::recode(res_plot_2$metric_name,\"Gene-peak \\n pair % recovered\" =\"Peak-gene \\n pair % recovered\")\n",
    "\n",
    "\n",
    "# get observed accuracy\n",
    "plot_ref_accuracy = TRUE\n",
    "\n",
    "truth_tbl <- get_pred_truth(dir_path)\n",
    "truth_df <- data.frame(metric_name = factor(c(\"Peak-gene \\n pair % recovered\",\"F1\"),\n",
    "                                            levels=c(\"Peak-gene \\n pair % recovered\",\"F1\")),\n",
    "                       Z = as.numeric(round(colMeans(truth_tbl[,1:2]),digits=2)))\n",
    "truth_df\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23e52f06",
   "metadata": {},
   "source": [
    "### Summary Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd012c62",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"mouse_skin_single_mod_n_splitted_metric_summary19_21ct_w_observed_rev1.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),c(\"Cell type ASW\",\"cLisi\"),c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Peak-gene \\n pair % recovered\",\"F1\"),c('log2(Runtime)',\"Graph \\n connectivity\"))\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "x_interval = c(5000,10000,15000) # unique(res_plot_3$var)\n",
    "ylim_list = list(c(0,0.6),c(0.2,0.7),c(0.45,0.7),c(0.85,1),c(0.4,1),c(0,1),c(0.1,0.7),c(0.05,0.5),c(0,1),c(6,15))\n",
    "nbreaks = list(5,5,5,6,6,5,8,5,5,8)\n",
    "counter=1\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "        #     stat_summary(fun.y=mean, geom=\"point\",aes(col=method_name)) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                           \n",
    "pdf(file.path(folder_path,file_name),width=35,height=15)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab63b44e",
   "metadata": {},
   "source": [
    "# Scenario 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8890704",
   "metadata": {},
   "source": [
    "## PBMC 2000 cells"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d87d691",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbbe719a",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = \"dataset/multiome_pbmc_10k/nmulti2000_7ct_vdepth_test/\"\n",
    "\n",
    "res_all <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3_metric.csv\")\n",
    "\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "# add _2 to method name and add type column\n",
    "multiome_guided_method_list<-c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\")\n",
    "unpaired_method_list<-c(\"seurat3_2\",\"rfigr_2\",\"rbindsc_2\",\"rliger_2\",\"glue_2\")\n",
    "res_plot <- res_all %>%  \n",
    "    mutate(method_name = ifelse(method %in% multiome_guided_method_list,method,paste0(method,\"_2\")))%>% \n",
    "    # assign method type\n",
    "    mutate(type = ifelse(method_name %in% unpaired_method_list,\"Multiome-guided\",\"Unpaired_multiome_splitted\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"percent_recovered_50kb\",\"f1\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Peak-gene \\n pair % recovered\",\"F1\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7fa88c8",
   "metadata": {},
   "source": [
    "### Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5db4e8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"pbmc_vdepth_summary19_7ct_fair_rev1.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),c(\"Cell type ASW\",\"cLisi\"),c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Peak-gene \\n pair % recovered\",\"F1\"),c('log2(Runtime)',\"Graph \\n connectivity\"))\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "x_interval = unique(res_plot_3$var)\n",
    "\n",
    "ylim_list = list(c(0.2,1),c(0.4,0.9),c(0.5,0.8),c(0.8,1),c(0.6,1),c(0,1),c(0.25,0.5),c(0.25,0.45),c(0,1),c(6,12))\n",
    "nbreaks = list(4,5,6,4,4,5,5,4,5,8)\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=35,height=15)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ad34261",
   "metadata": {},
   "source": [
    "## BMMC 2000 cells"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60fbc936",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c82f278d",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = \"dataset/bmmc/nmulti2000_all_ct_vdepth/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "res_all <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3wbatch_metric.csv\")\n",
    "table(res_all$method)\n",
    "table(res_all$out)\n",
    "\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "multiome_guided_method_list<-c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\")\n",
    "unpaired_method_list<-c(\"seurat3_2\",\"rfigr_2\",\"rbindsc_2\",\"rliger_2\",\"glue_2\")\n",
    "res_plot <- res_all %>% \n",
    "    mutate(method_name = ifelse(method %in% multiome_guided_method_list,method,paste0(method,\"_2\")))%>% \n",
    "    # assign method type\n",
    "    mutate(type = ifelse(method_name %in% unpaired_method_list,\"Multiome-guided\",\"Unpaired_multiome_splitted\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"percent_recovered_50kb\",\"f1\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Peak-gene \\n pair % recovered\",\"F1\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7eba469",
   "metadata": {},
   "source": [
    "### Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f848bea",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"bmmc_nmulti2000_vdepth_metric_summary19_21ct_fair_rev1.pdf\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),c(\"Cell type ASW\",\"cLisi\"),c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Peak-gene \\n pair % recovered\",\"F1\"),c('log2(Runtime)',\"Graph \\n connectivity\"))\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "x_interval = unique(res_plot_3$var)\n",
    "ylim_list = list(c(0,0.5),c(0.2,0.7),c(0.45,0.7),c(0.85,1),c(0.4,1),c(0,1),c(0.05,0.3),c(0.05,0.3),c(0,1),c(6,12))\n",
    "nbreaks = list(5,5,5,6,6,5,5,5,5,8)\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "        #     stat_summary(fun.y=mean, geom=\"point\",aes(col=method_name)) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=35,height=15)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e37ac01",
   "metadata": {},
   "source": [
    "## BMMC 4000 cells"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07a3c070",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acf1217b",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = \"/project/mingyaolpc/myylee/scmint/methods_eval/dataset/bmmc/nmulti4000_all_ct_vdepth/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "res_all <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3_metric.csv\")\n",
    "table(res_all$method)\n",
    "table(res_all$out)\n",
    "\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "multiome_guided_method_list<-c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\")\n",
    "unpaired_method_list<-c(\"seurat3_2\",\"rfigr_2\",\"rbindsc_2\",\"rliger_2\",\"glue_2\")\n",
    "res_plot <- res_all %>% \n",
    "    #filter(var %in% c(1000,3000,5000)) %>% \n",
    "    mutate(method_name = ifelse(method %in% multiome_guided_method_list,method,paste0(method,\"_2\")))%>% \n",
    "    # assign method type\n",
    "    mutate(type = ifelse(method_name %in% unpaired_method_list,\"Multiome-guided\",\"Unpaired_multiome_splitted\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"percent_recovered_50kb\",\"f1\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Peak-gene \\n pair % recovered\",\"F1\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b99fb50",
   "metadata": {},
   "source": [
    "### Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73ef2c02",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"bmmc_nmulti4000_vdepth_metric_summary19_21ct_fair_rev1.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),c(\"Cell type ASW\",\"cLisi\"),c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Peak-gene \\n pair % recovered\",\"F1\"),c('log2(Runtime)',\"Graph \\n connectivity\"))\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "x_interval = unique(res_plot_3$var)\n",
    "ylim_list = list(c(0,0.5),c(0.2,0.7),c(0.45,0.7),c(0.85,1),c(0.4,1),c(0,1),c(0.05,0.3),c(0.05,0.3),c(0,1),c(6,12))\n",
    "\n",
    "nbreaks = list(5,5,5,6,6,5,5,5,5,8)\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=35,height=15) \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e05ef5d3",
   "metadata": {},
   "source": [
    "## BMMC Increasing number of cells; 100% depth"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2996debc",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "979e004f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in metric in ct3_metric.csv\n",
    "dir_path = \"dataset/bmmc/bmmc_vcells_intervals/\"\n",
    "\n",
    "res_all <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3_metric.csv\")\n",
    "\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "multiome_guided_method_list<-c(\"multivi\",\"seurat4\",\"cobolt\")\n",
    "unpaired_method_list<-c(\"seurat3_2\",\"rfigr_2\",\"rbindsc_2\",\"rliger_2\")\n",
    "res_plot <- res_all %>% \n",
    "    mutate(method_name = ifelse(method %in% multiome_guided_method_list,method,paste0(method,\"_2\")))%>% \n",
    "    # assign method type\n",
    "    mutate(type = ifelse(method_name %in% unpaired_method_list,\"Multiome-guided\",\"Unpaired_multiome_split\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"percent_recovered_50kb\",\"f1\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"Seurat v3\",\"Seurat v4\")\n",
    "row_order = c(\"seurat3_2\",\"seurat4\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Unpaired with \\n Multiome-split\",1),rep(\"Multiome-guided\",1))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Peak-gene \\n pair % recovered\",\"F1\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-split\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f246cb8",
   "metadata": {},
   "source": [
    "### Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b80556fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"bmmc_nmulti4000_vdepth_10interv_summary19_21ct_fair.pdf\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),c(\"Cell type ASW\",\"cLisi\"),c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Peak-gene \\n pair % recovered\",\"F1\"),c('log2(Runtime)',\"Graph \\n connectivity\"))\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "x_interval = unique(res_plot_3$var)\n",
    "\n",
    "ylim_list = list(c(0,0.5),c(0.2,0.7),c(0.45,0.7),c(0.85,1),c(0.4,1),c(0,1),c(0.05,0.3),c(0.05,0.3),c(0,1),c(6,14))\n",
    "nbreaks = list(5,5,5,6,6,5,5,5,5,8)\n",
    "counter=1\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            scale_color_manual(values= hex,labels=names(hex)) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x){x$name}) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=35,height=15)              \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71ee924c",
   "metadata": {},
   "source": [
    "## BMMC Increasing depth; 4000 cells"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34d5e5fc",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0152d8f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = \"dataset/bmmc/bmmc_vdepth_intervals/\"\n",
    "\n",
    "res_all <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3_metric.csv\")\n",
    "\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "multiome_guided_method_list<-c(\"multivi\",\"seurat4\",\"cobolt\")\n",
    "unpaired_method_list<-c(\"seurat3_2\",\"rfigr_2\",\"rbindsc_2\",\"rliger_2\")\n",
    "res_plot <- res_all %>% \n",
    "    mutate(method_name = ifelse(method %in% multiome_guided_method_list,method,paste0(method,\"_2\")))%>% \n",
    "    # assign method type\n",
    "    mutate(type = ifelse(method_name %in% unpaired_method_list,\"Multiome-guided\",\"Unpaired_multiome_split\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"percent_recovered_50kb\",\"f1\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\n",
    "                 \"Seurat v3\",\"bindSC\",\"FigR\", \"Liger\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",3),rep(\"Unpaired with \\n Multiome-split\",4))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Peak-gene \\n pair % recovered\",\"F1\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-split\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37b585c9",
   "metadata": {},
   "source": [
    "### Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "829765f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"bmmc_nmulti4000_vdepth_10interv_summary19_21ct_fair.pdf\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),c(\"Cell type ASW\",\"cLisi\"),c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Peak-gene \\n pair % recovered\",\"F1\"),c('log2(Runtime)',\"Graph \\n connectivity\"))\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "x_interval = unique(res_plot_3$var)\n",
    "\n",
    "ylim_list = list(c(0,0.5),c(0.2,0.7),c(0.45,0.7),c(0.85,1),c(0.4,1),c(0,1),c(0.05,0.3),c(0.05,0.3),c(0,1),c(6,14))\n",
    "nbreaks = list(5,5,5,6,6,5,5,5,5,8)\n",
    "counter=1\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            scale_color_manual(values= hex,labels=names(hex)) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x){x$name}) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "#pdf(file.path(folder_path,file_name),width=35,height=15) \n",
    "pdf(\"tmp_fig4.pdf\",width=35,height=15) \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab996031",
   "metadata": {},
   "source": [
    "# Scenario 3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1591de2b",
   "metadata": {},
   "source": [
    "## Technical batch effect"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd355f60",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eadaccdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in metric in ct3wbatch_metric.csv, then ct3wbatch_metric_sample_batch.csv, and then combine\n",
    "dir_path = \"dataset/bmmc/multiome_ncells_pmat_batch/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "res_all <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3wbatch_metric.csv\")\n",
    "table(res_all$method)\n",
    "table(res_all$out)\n",
    "\n",
    "\n",
    "res_all_batch <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3wbatch_metric_sample_batch.csv\")\n",
    "table(res_all_batch$method)\n",
    "table(res_all_batch$out)\n",
    "\n",
    "res_all_comb <- left_join(res_all,res_all_batch[2:10])\n",
    "res_all_comb$log2_runtime <- log2(res_all_comb$runtime)\n",
    "\n",
    "table(is.na(res_all_comb$batch_b_saw))\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "multiome_guided_method_list<-c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\")\n",
    "unpaired_method_list<-c(\"seurat3_2\",\"rfigr_2\",\"rbindsc_2\",\"rliger_2\",\"glue_2\")\n",
    "res_plot <- res_all_comb %>% \n",
    "    filter(var %in% c(1000,3000,5000)) %>% \n",
    "    mutate(method_name = ifelse(method %in% multiome_guided_method_list,method,paste0(method,\"_2\")))%>% \n",
    "    # assign method type\n",
    "    mutate(type = ifelse(method_name %in% unpaired_method_list,\"Multiome-guided\",\"Unpaired_multiome_split\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"batch_b_saw\",\"batch_kbet\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "res_plot = res_plot %>% filter(ct_ref ==\"ct3wbatch_metric.csv\")\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-split\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Site ASW\",\"Site kBET\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-split\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0765e0d9",
   "metadata": {},
   "source": [
    "### Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd25300d",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"bmmc_tbatch_metric_summary19_21ct_fair_rev1.pdf\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),\n",
    "                   c(\"Cell type ASW\",\"cLisi\"),\n",
    "                   c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Site ASW\",\"Site kBET\"),\n",
    "                   c(\"Graph \\n connectivity\",'log2(Runtime)'))\n",
    "\n",
    "ylim_list = list(c(0.1,0.7),c(0.3,0.8),\n",
    "                 c(0.4,0.7), c(0.8,1.0),\n",
    "                 c(0.5,1),c(0.1,0.8),\n",
    "                 c(0.2 ,1),c(0.1,1),\n",
    "                 c(0,1),c(6,12))\n",
    "nbreaks = list(6,5,\n",
    "               6,5,\n",
    "               5,6,\n",
    "               7,8,\n",
    "               5,6)\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "x_interval = unique(res_plot_3$var)\n",
    "subplots = list()\n",
    "counter = 1\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "            theme(#strip.text.x = element_text(size = 35,angle=90,colour = \"black\"),\n",
    "                  strip.text.x = element_blank(),\n",
    "                  strip.text.y = element_text(size = 35,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 30,angle = 90),axis.text.y = element_text(size = 30),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(2, \"lines\"),panel.spacing.x = unit(1, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=30,height=10)             \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6cbf167",
   "metadata": {},
   "source": [
    "## Biological batch effect"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d948679",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d80b71b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in metric in ct3wbatch_metric.csv, then ct3wbatch_metric_sample_batch.csv, and then combine\n",
    "dir_path = \"dataset/bmmc/bmmc_biological_batch_test/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "res_all <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3wbatch_metric.csv\")\n",
    "\n",
    "res_all_batch <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3wbatch_metric_sample_batch.csv\")\n",
    "\n",
    "res_all_comb <- left_join(res_all,res_all_batch[2:10])\n",
    "res_all_comb$log2_runtime <- log2(res_all_comb$runtime)\n",
    "\n",
    "table(is.na(res_all_comb$batch_b_saw))\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "multiome_guided_method_list<-c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\")\n",
    "unpaired_method_list<-c(\"seurat3_2\",\"rfigr_2\",\"rbindsc_2\",\"rliger_2\",\"glue_2\")\n",
    "res_plot <- res_all_comb %>% \n",
    "    filter(var %in% c(1000,3000,5000)) %>% \n",
    "    mutate(method_name = ifelse(method %in% multiome_guided_method_list,method,paste0(method,\"_2\")))%>% \n",
    "    # assign method type\n",
    "    mutate(type = ifelse(method_name %in% unpaired_method_list,\"Multiome-guided\",\"Unpaired_multiome_split\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"batch_b_saw\",\"batch_kbet\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "res_plot = res_plot %>% filter(ct_ref ==\"ct3wbatch_metric.csv\")\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-split\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Donor ASW\",\"Donor kBET\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-split\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48b67588",
   "metadata": {},
   "source": [
    "### Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65cb3c33",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"summary_bmmc_biological_batch_test_rev1.pdf\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),\n",
    "                   c(\"Cell type ASW\",\"cLisi\"),\n",
    "                   c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Donor ASW\",\"Donor kBET\"),\n",
    "                   c(\"Graph \\n connectivity\",'log2(Runtime)'))\n",
    "\n",
    "ylim_list = list(c(0.1,0.7),c(0.3,0.8),\n",
    "                 c(0.4,0.7), c(0.8,1.0),\n",
    "                 c(0.5,1),c(0.1,0.8),\n",
    "                 c(0.2,1),c(0.1,1),\n",
    "                 c(0,1),c(6,12))\n",
    "nbreaks = list(6,5,\n",
    "               6,5,\n",
    "               5,6,\n",
    "               8,9,\n",
    "               5,6)\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "x_interval = unique(res_plot_3$var)\n",
    "subplots = list()\n",
    "counter = 1\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "            theme(#strip.text.x = element_text(size = 35,angle=90,colour = \"black\"),\n",
    "                  strip.text.x = element_blank(),\n",
    "                  strip.text.y = element_text(size = 35,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 30,angle = 90),axis.text.y = element_text(size = 30),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(2, \"lines\"),panel.spacing.x = unit(1, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=30,height=10)             \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85fbad34",
   "metadata": {},
   "source": [
    "## Complex test #1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba743bbc",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ae0d1ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in metric in ct3wbatch_metric.csv, then ct3wbatch_metric_sample_batch.csv, and then combine\n",
    "dir_path = \"dataset/bmmc/bmmc_complex1_test/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "res_all <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3_metric.csv\")\n",
    "\n",
    "res_all_batch <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3_metric_sample_batch.csv\")\n",
    "\n",
    "res_all_comb <- left_join(res_all,res_all_batch[2:10])\n",
    "res_all_comb$log2_runtime <- log2(res_all_comb$runtime)\n",
    "\n",
    "table(is.na(res_all_comb$batch_b_saw))\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "multiome_guided_method_list<-c(\"multivi\",\"seurat4\",\"cobolt\",\"seurat4int\",\"scmomat\")\n",
    "unpaired_method_list<-c(\"seurat3_2\",\"rfigr_2\",\"rbindsc_2\",\"rliger_2\",\"glue_2\")\n",
    "res_plot <- res_all_comb %>% \n",
    "    mutate(method_name = ifelse(method %in% multiome_guided_method_list,method,paste0(method,\"_2\")))%>% \n",
    "    # assign method type\n",
    "    mutate(type = ifelse(method_name %in% unpaired_method_list,\"Multiome-guided\",\"Unpaired_multiome_splitted\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"batch_b_saw\",\"batch_kbet\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "res_plot = res_plot %>% filter(ct_ref ==\"ct3_metric.csv\")\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v4 integrate\",\"Seurat v3\",\"bindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"bindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat4int\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",5),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Donor ASW\",\"Donor kBET\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9a161e8",
   "metadata": {},
   "source": [
    "### Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cab2307",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"summary_bmmc_complex1_test_rev1.pdf\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),\n",
    "                   c(\"Cell type ASW\",\"cLisi\"),\n",
    "                   c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Donor ASW\",\"Donor kBET\"),\n",
    "                   c(\"Graph \\n connectivity\",'log2(Runtime)'))\n",
    "\n",
    "#ylim_list = list(c(0,0.7),c(0,0.8),\n",
    "ylim_list = list(c(0.2,0.65),c(0.4,0.8),\n",
    "                 c(0.4,0.7), c(0.7,1.0),\n",
    "                 c(0.5,1),c(0,0.8),\n",
    "                 c(0.4,1),c(0,0.9),\n",
    "                 c(0.5,1),c(6,12))\n",
    "#nbreaks = list(6,6,\n",
    "nbreaks = list(6,6,\n",
    "               6,6,\n",
    "               5,8,\n",
    "               8,9,\n",
    "               5,6)\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "res_plot_3$method <- factor(res_plot_3$method,\n",
    "                           levels=c(\"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\",\n",
    "                                   \"cobolt\",\"multivi\",\"seurat4\",\"seurat4int\",\"scmomat\"))\n",
    "\n",
    "x_interval = unique(res_plot_3$var)\n",
    "subplots = list()\n",
    "counter = 1\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=value,group=method,col=method_name,shape=type)) +\n",
    "            #geom_point(size=3) +\n",
    "            geom_boxplot() +\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            #scale_color_manual(values= hex,labels=names(hex)) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "            theme(#strip.text.x = element_text(size = 35,angle=90,colour = \"black\"),\n",
    "                  strip.text.x = element_blank(),\n",
    "                  strip.text.y = element_text(size = 35,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 30,angle = 90),axis.text.y = element_text(size = 30),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(2, \"lines\"),panel.spacing.x = unit(1, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=30,height=10)             \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27529135",
   "metadata": {},
   "source": [
    "## Complex test #2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdb2f1a2",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6c07abe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in metric in ct3wbatch_metric.csv, then ct3wbatch_metric_sample_batch.csv, and then combine\n",
    "dir_path = \"dataset/bmmc/bmmc_complex2_test/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "res_all <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3_metric.csv\")\n",
    "\n",
    "res_all_batch <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3_metric_sample_batch.csv\")\n",
    "\n",
    "res_all_comb <- left_join(res_all,res_all_batch[2:10])\n",
    "res_all_comb$log2_runtime <- log2(res_all_comb$runtime)\n",
    "\n",
    "table(is.na(res_all_comb$batch_b_saw))\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "multiome_guided_method_list<-c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat4int\")\n",
    "unpaired_method_list<-c(\"seurat3_2\",\"rfigr_2\",\"rbindsc_2\",\"rliger_2\",\"glue_2\")\n",
    "res_plot <- res_all_comb %>% \n",
    "    mutate(method_name = ifelse(method %in% multiome_guided_method_list,method,paste0(method,\"_2\")))%>% \n",
    "    # assign method type\n",
    "    mutate(type = ifelse(method_name %in% unpaired_method_list,\"Multiome-guided\",\"Unpaired_multiome_splitted\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"batch_b_saw\",\"batch_kbet\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "res_plot = res_plot %>% filter(ct_ref ==\"ct3_metric.csv\")\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v4 integrate\",\"Seurat v3\",\"bindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"bindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat4int\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",5),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Donor ASW\",\"Donor kBET\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "827e4a73",
   "metadata": {},
   "source": [
    "### Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a776b72",
   "metadata": {},
   "outputs": [],
   "source": [
    "# color palette for methods \n",
    "file_name <- \"summary_bmmc_complex2_test.pdf\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "subplot_idx = list(c(\"ARI\",\"NMI\"),\n",
    "                   c(\"Cell type ASW\",\"cLisi\"),\n",
    "                   c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Donor ASW\",\"Donor kBET\"),\n",
    "                   c(\"Graph \\n connectivity\",'log2(Runtime)'))\n",
    "\n",
    "ylim_list = list(c(0.2,0.7),c(0.45,0.8), \n",
    "                 c(0.4,0.7), c(0.85,1.0), \n",
    "                 c(0.5,1),c(0,0.8),\n",
    "                 c(0.4,1),c(0,0.9),\n",
    "                 c(0.5,1),c(6,14))\n",
    "\n",
    "nbreaks = list(6,6, \n",
    "               6,8,\n",
    "               5,8,\n",
    "               8,9,\n",
    "               5,8)\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "res_plot_3$method <- factor(res_plot_3$method,\n",
    "                           levels=c(\"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\",\n",
    "                                   \"cobolt\",\"multivi\",\"seurat4\",\"seurat4int\",\"scmomat\"))\n",
    "\n",
    "x_interval = unique(res_plot_3$var)\n",
    "subplots = list()\n",
    "counter = 1\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=value,group=method,col=method_name,shape=type)) +\n",
    "            #geom_point(size=3) +\n",
    "            geom_boxplot() +\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "            theme(#strip.text.x = element_text(size = 35,angle=90,colour = \"black\"),\n",
    "                  strip.text.x = element_blank(),\n",
    "                  strip.text.y = element_text(size = 35,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 30,angle = 90),axis.text.y = element_text(size = 30),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(2, \"lines\"),panel.spacing.x = unit(1, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]]),\n",
    "            scale_y_continuous(limits = ylim_list[[counter+1]],n.breaks = nbreaks[[counter+1]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 2\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x){x$name}) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=30,height=10)             \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()\n",
    "                              "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81bd5c6d",
   "metadata": {},
   "source": [
    "# Scenario 4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccfe875b",
   "metadata": {},
   "source": [
    "## PBMC"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9347039b",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f4b6c4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = \"dataset/pbmc/single_modality_fixed_missing_ct/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "\n",
    "res_all<- plot_rare_cell_result(dir_path,\n",
    "                      runtime_exists=T,prediction_exists=F)\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "res_all \n",
    "\n",
    "plot_ref_accuracy = FALSE\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "unpaired_method_list<-c(\"seurat3\",\"rfigr\",\"rbindsc\",\"rliger\",\"glue\")\n",
    "res_plot <- res_all %>% \n",
    "    filter(var %in% c(1000,3000,6000)) %>% \n",
    "    mutate(method_name = ifelse(out==\"results_single_same_cell_number\",paste0(method,\"_2\"),method)) %>%\n",
    "    # assign method type \n",
    "    mutate(type = ifelse(out!=\"results_single_same_cell_number\",\n",
    "                         ifelse(method%in%unpaired_method_list,\"Unpaired\",\"Multiome-guided\"),\n",
    "                         \"Unpaired_multiome_splitted\"),\n",
    "           ct = gsub(\"_metric.csv\",\"\",ct_ref),\n",
    "          metric=paste0(orig,\"_\",ct))\n",
    "\n",
    "metric_cond = c(\"noMiss\",\"rnaMissNK\",\"atacMissNK\")\n",
    "cts = c(\"ct3_NK\")\n",
    "cts_name = c(\"NK\")\n",
    "\n",
    "metric_sel = paste0(rep(metric_cond,each=length(cts)),\"_\",rep(cts_name,length(metric_cond)))\n",
    "\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"bindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"bindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = paste0(rep(metric_cond,each=length(cts_name)),\"_\",rep(cts_name,length(metric_cond)))\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    separate(metric, c('cond', 'fill',\"ct\",\"dataset\"),sep = \"_\") %>%\n",
    "    mutate(nTarget = TP+FN) %>%\n",
    "    mutate(perAcc = TP/nTarget) %>%\n",
    "    select(F1,nTarget,perAcc,log2_runtime,method_name,var,cond,ct,dataset) %>%\n",
    "    gather(key = \"key\",\n",
    "               value = \"value\", -method_name,-var,-cond,-ct,-dataset)%>%\n",
    "    mutate(method = method_name,metric_type = paste0(cond,\"_\",ct)) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "684e4e5a",
   "metadata": {},
   "source": [
    "### Plots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90074707",
   "metadata": {},
   "source": [
    "#### RNA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe0977d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"pbmc_rare_cell_1_f1_try2_rna.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"noMiss_NK\"),\n",
    "                   c(\"atacMissNK_NK\")\n",
    "                  )\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"scRNA\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "x_interval = c(1000,3000,6000)\n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "        #     stat_summary(fun.y=mean, geom=\"point\",aes(col=method_name)) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "    \n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=15,height=10)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdc3eb2a",
   "metadata": {},
   "source": [
    "#### ATAC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7affb0af",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"pbmc_rare_cell_1_f1_try2_atac.pdf\"\n",
    "library(gridExtra)\n",
    "\n",
    "subplot_idx = list(c(\"noMiss_NK\"),\n",
    "                   c(\"rnaMissNK_NK\")\n",
    "                  )\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"snATAC\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "x_interval = c(1000,3000,6000) \n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "    \n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=15,height=10)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2235e6fa",
   "metadata": {},
   "source": [
    "## SHARE-seq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee5d30fe",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad892351",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = \"/dataset/mouse_skin/single_modality_fixed_missing_ct/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "res_all<- plot_rare_cell_result(dir_path,\n",
    "                      runtime_exists=T,prediction_exists=F)\n",
    "plot_ref_accuracy = FALSE\n",
    "\n",
    "res_all <- res_all %>% filter(out %in% c(\"results_single_mod\",\"results_single_same_cell_number\"))\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "res_all \n",
    "\n",
    "# add _2 to method name and add type column\n",
    "unpaired_method_list<-c(\"seurat3\",\"rfigr\",\"rbindsc\",\"rliger\",\"glue\")\n",
    "res_plot <- res_all %>% \n",
    "    filter(var %in% c(1000,3000,6000)) %>% \n",
    "    mutate(method_name = ifelse(out==\"results_single_same_cell_number\",paste0(method,\"_2\"),method)) %>%\n",
    "    # assign method type \n",
    "    mutate(type = ifelse(out!=\"results_single_same_cell_number\",\n",
    "                         ifelse(method%in%unpaired_method_list,\"Unpaired\",\"Multiome-guided\"),\n",
    "                         \"Unpaired_multiome_splitted\"),\n",
    "           ct = gsub(\"_metric.csv\",\"\",ct_ref),\n",
    "          metric=paste0(orig,\"_\",ct))\n",
    "\n",
    "metric_cond = c(\"noMiss\",\"rnaMissHS\",\"atacMissHS\",\"rnaMissEndo\",\"atacMissEndo\",\n",
    "                \"rnaMissTwo\",\"atacMissTwo\",\"eachMissOne\",\"eachMissOneAlt\")\n",
    "cts = c(\"10k_HS\",\"10k_Endo\")\n",
    "cts_name = c(\"HS\",\"Endo\")\n",
    "\n",
    "metric_sel = paste0(rep(metric_cond,each=length(cts)),\"_\",rep(cts_name,length(metric_cond)))\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"bindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"bindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = paste0(rep(metric_cond,each=length(cts_name)),\"_\",rep(cts_name,length(metric_cond)))\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    separate(metric, c('cond', 'fill',\"ct\",\"dataset\"),sep = \"_\") %>%\n",
    "    mutate(nTarget = TP+FN) %>%\n",
    "    mutate(perAcc = TP/nTarget) %>%\n",
    "    select(F1,nTarget,perAcc,log2_runtime,method_name,var,cond,ct,dataset) %>%\n",
    "    gather(key = \"key\",\n",
    "               value = \"value\", -method_name,-var,-cond,-ct,-dataset)%>%\n",
    "    mutate(method = method_name,metric_type = paste0(cond,\"_\",ct)) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n",
    "\n",
    "res_plot_2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ba2f9de",
   "metadata": {},
   "source": [
    "### Plots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a183c30",
   "metadata": {},
   "source": [
    "#### scRNA F1 HS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3faddae9",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"shareseq_mouse_skin_rare_cell_1_f1_try2_rna_HS.pdf\"\n",
    "subplot_idx = list(c(\"noMiss_HS\"),\n",
    "                   c(\"atacMissHS_HS\"),\n",
    "                   c(\"atacMissTwo_HS\")\n",
    "                  )\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"scRNA\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "x_interval = c(1000,3000,6000) \n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "    \n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=20,height=10)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()\n",
    "                              \n",
    "   \n",
    "                    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b40592d",
   "metadata": {},
   "source": [
    "#### scRNA F1 Endo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c43f9c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"shareseq_mouse_skin_rare_cell_1_f1_try2_rna_Endo.pdf\"\n",
    "subplot_idx = list(c(\"noMiss_Endo\"),\n",
    "                   c(\"atacMissEndo_Endo\"),\n",
    "                   c(\"atacMissTwo_Endo\")\n",
    "                  )\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"scRNA\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "x_interval = c(1000,3000,6000) \n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "    \n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=20,height=10)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()\n",
    "                              \n",
    "   \n",
    "                    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "806f88b0",
   "metadata": {},
   "source": [
    "#### snATAC F1 HS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40e168dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "file_name <- \"shareseq_mouse_skin_rare_cell_1_f1_try2_atac_HS.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"noMiss_HS\"),\n",
    "                   c(\"rnaMissHS_HS\"),\n",
    "                   c(\"rnaMissTwo_HS\")\n",
    "                  )\n",
    "\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"snATAC\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "x_interval = c(1000,3000,6000) \n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "        #     stat_summary(fun.y=mean, geom=\"point\",aes(col=method_name)) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "    \n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=20,height=10)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f13a5daf",
   "metadata": {},
   "source": [
    "#### snATAC F1 Endo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15bf7a43",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "file_name <- \"shareseq_mouse_skin_rare_cell_1_f1_try2_atac_Endo.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"noMiss_Endo\"),\n",
    "                   c(\"multiMissEndo_Endo\"),\n",
    "                   c(\"atacOnlyEndo_Endo\")\n",
    "                  )\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"snATAC\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "x_interval = c(1000,3000,6000) \n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=mean_f,group=method,col=method_name,shape=type)) + \n",
    "            geom_point(size=3) +\n",
    "            geom_line() +\n",
    "            geom_errorbar(aes(ymin=mean_f-sd_f, ymax=mean_f+sd_f),\n",
    "                         position=position_dodge(0.05),alpha=0.2)+\n",
    "            scale_shape_manual(values = c(15,16,17))+#c(0,1,2))+\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "    \n",
    "            theme(strip.text.x = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  strip.text.y = element_text(size = 30,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 25,angle = 90),axis.text.y = element_text(size = 25),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(5, \"lines\"),panel.spacing.x = unit(2, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=20,height=10)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca257fdf",
   "metadata": {},
   "source": [
    "# Scenario 5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "630fb295",
   "metadata": {},
   "source": [
    "## PBMC"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57db0437",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bad2169",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = \"dataset/pbmc/multiome_fixed_missing_ct/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "res_all<- plot_rare_cell_result(dir_path,runtime_exists=T,prediction_exists=F)\n",
    "plot_ref_accuracy = FALSE\n",
    "\n",
    "res_all <- res_all %>% filter(out %in% c(\"results_single_mod\",\"results_single_same_cell_number\"))\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "res_all \n",
    "\n",
    "# add _2 to method name and add type column\n",
    "unpaired_method_list<-c(\"seurat3\",\"rfigr\",\"rbindsc\",\"rliger\",\"glue\")\n",
    "res_plot <- res_all %>% \n",
    "    filter(var %in% c(1000,3000,6000)) %>% \n",
    "    mutate(method_name = ifelse(out==\"results_single_same_cell_number\",paste0(method,\"_2\"),method)) %>%\n",
    "    # assign method type \n",
    "    mutate(type = ifelse(out!=\"results_single_same_cell_number\",\n",
    "                         ifelse(method%in%unpaired_method_list,\"Unpaired\",\"Multiome-guided\"),\n",
    "                         \"Unpaired_multiome_splitted\"),\n",
    "           ct = gsub(\"_metric.csv\",\"\",ct_ref),\n",
    "          metric=paste0(orig,\"_\",ct))\n",
    "\n",
    "metric_cond = c(\"noMiss\",\"multiMissNK\",\"rnaOnlyNK\",\"atacOnlyNK\")\n",
    "\n",
    "cts_name = c(\"NK\")\n",
    "\n",
    "metric_sel = paste0(rep(metric_cond,each=length(cts_name)),\"_\",rep(cts_name,length(metric_cond)))\n",
    "\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = paste0(rep(metric_cond,each=length(cts_name)),\"_\",rep(cts_name,length(metric_cond)))\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    separate(metric, c('cond', 'fill',\"ct\",\"dataset\"),sep = \"_\") %>%\n",
    "    mutate(nTarget = TP+FN) %>%\n",
    "    mutate(perAcc = TP/nTarget) %>%\n",
    "    select(F1,nTarget,perAcc,log2_runtime,method_name,var,cond,ct,dataset) %>%\n",
    "    gather(key = \"key\",\n",
    "               value = \"value\", -method_name,-var,-cond,-ct,-dataset)%>%\n",
    "    mutate(method = method_name,metric_type = paste0(cond,\"_\",ct)) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n",
    "\n",
    "res_plot_2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d63745d3",
   "metadata": {},
   "source": [
    "### Plots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ad4d234",
   "metadata": {},
   "source": [
    "#### boxplot RNA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "502b2a8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"pbmc_rare_cell_2_f1_try2_rna_boxplot.pdf\"\n",
    "library(gridExtra)\n",
    "\n",
    "subplot_idx = list(c(\"noMiss_NK\"),\n",
    "                   c(\"multiMissNK_NK\"),\n",
    "                   c(\"rnaOnlyNK_NK\")\n",
    "                  )\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"scRNA\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "# plot only the unpaired and multiome-guided categories\n",
    "res_plot_3 <- res_plot_3 %>% filter(type %in% c(\"Unpaired\",\"Multiome-guided\"))\n",
    "\n",
    "x_interval = c(3000) \n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=value,group=method,col=method_name,shape=type)) +\n",
    "            geom_boxplot(lwd=1) +\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            scale_color_manual(values=hex[method_names]) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_blank(),\n",
    "                  strip.text.y = element_text(size = 50,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 30,angle = 90),\n",
    "                  axis.text.y = element_text(size = 50),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(2, \"lines\"),panel.spacing.x = unit(1, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=25,height=8)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "847f802f",
   "metadata": {},
   "source": [
    "#### boxplot ATAC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e43a648",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"pbmc_rare_cell_2_f1_try2_atac_boxplot.pdf\"\n",
    "library(gridExtra)\n",
    "\n",
    "subplot_idx = list(c(\"noMiss_NK\"),\n",
    "                   c(\"multiMissNK_NK\"),\n",
    "                   c(\"atacOnlyNK_NK\")\n",
    "                  )\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"snATAC\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "# plot only the unpaired and multiome-guided categories\n",
    "res_plot_3 <- res_plot_3 %>% filter(type %in% c(\"Unpaired\",\"Multiome-guided\"))\n",
    "\n",
    "x_interval = c(3000) \n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=value,group=method,col=method_name,shape=type)) +\n",
    "            geom_boxplot(lwd=1) +\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            scale_color_manual(values=hex[method_names]) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_blank(),\n",
    "                  strip.text.y = element_text(size = 50,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 30,angle = 90),\n",
    "                  axis.text.y = element_text(size = 50),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(2, \"lines\"),panel.spacing.x = unit(1, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=25,height=8)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8288f041",
   "metadata": {},
   "source": [
    "## SHARE-seq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5da3ac7a",
   "metadata": {},
   "source": [
    "### Load result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f68d7888",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = \"/project/mingyaolpc/myylee/scmint/methods_eval/dataset/mouse_skin/multiome_fixed_missing_ct/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    " \n",
    "res_all<- plot_rare_cell_result(dir_path,runtime_exists=T,prediction_exists=F)\n",
    "plot_ref_accuracy = FALSE\n",
    "\n",
    "res_all <- res_all %>% filter(out %in% c(\"results_single_mod\",\"results_single_same_cell_number\"))\n",
    "res_all$log2_runtime <- log2(res_all$runtime)\n",
    "res_all \n",
    "\n",
    "# add _2 to method name and add type column\n",
    "unpaired_method_list<-c(\"seurat3\",\"rfigr\",\"rbindsc\",\"rliger\",\"glue\")\n",
    "res_plot <- res_all %>% \n",
    "    filter(var %in% c(1000,3000,6000)) %>% \n",
    "    mutate(method_name = ifelse(out==\"results_single_same_cell_number\",paste0(method,\"_2\"),method)) %>%\n",
    "    # assign method type \n",
    "    mutate(type = ifelse(out!=\"results_single_same_cell_number\",\n",
    "                         ifelse(method%in%unpaired_method_list,\"Unpaired\",\"Multiome-guided\"),\n",
    "                         \"Unpaired_multiome_splitted\"),\n",
    "           ct = gsub(\"_metric.csv\",\"\",ct_ref),\n",
    "          metric=paste0(orig,\"_\",ct))\n",
    "\n",
    "metric_cond = c(\"noMiss\",\"multiMissHS\",\"rnaOnlyHS\",\"atacOnlyHS\",\"multiOnlyHS\",\n",
    "                \"multiMissEndo\",\"rnaOnlyEndo\",\"atacOnlyEndo\",\"multiOnlyEndo\")\n",
    "cts = c(\"10k_HS\",\"10k_Endo\")\n",
    "cts_name = c(\"HS\",\"Endo\")\n",
    "\n",
    "metric_sel = paste0(rep(metric_cond,each=length(cts)),\"_\",rep(cts_name,length(metric_cond)))\n",
    "\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v3\",\"bindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"bindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",4),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = paste0(rep(metric_cond,each=length(cts_name)),\"_\",rep(cts_name,length(metric_cond)))\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    separate(metric, c('cond', 'fill',\"ct\",\"dataset\"),sep = \"_\") %>%\n",
    "    mutate(nTarget = TP+FN) %>%\n",
    "    mutate(perAcc = TP/nTarget) %>%\n",
    "    select(F1,nTarget,perAcc,log2_runtime,method_name,var,cond,ct,dataset) %>%\n",
    "    gather(key = \"key\",\n",
    "               value = \"value\", -method_name,-var,-cond,-ct,-dataset)%>%\n",
    "    mutate(method = method_name,metric_type = paste0(cond,\"_\",ct)) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n",
    "\n",
    "\n",
    "res_plot_2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95e4f442",
   "metadata": {},
   "source": [
    "### Plots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "127f002a",
   "metadata": {},
   "source": [
    "#### box plot - RNA HS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b90cbd03",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"shareseq_mouse_skin_rare_cell_2_f1_try2_boxplot_rna_HS.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"noMiss_HS\"),\n",
    "                   c(\"multiMissHS_HS\"),\n",
    "                   c(\"rnaOnlyHS_HS\")\n",
    "                  )\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"scRNA\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "\n",
    "# plot only the unpaired and multiome-guided categories\n",
    "res_plot_3 <- res_plot_3 %>% filter(type %in% c(\"Unpaired\",\"Multiome-guided\"))\n",
    "\n",
    "x_interval = c(3000) \n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=value,group=method,col=method_name,shape=type)) +\n",
    "            geom_boxplot(lwd=1) +\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            scale_color_manual(values=hex[method_names]) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_blank(),\n",
    "                  strip.text.y = element_text(size = 50,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 30,angle = 90),\n",
    "                  axis.text.y = element_text(size = 50),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(2, \"lines\"),panel.spacing.x = unit(1, \"lines\"),\n",
    "                  panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=25,height=8)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b33c64e3",
   "metadata": {},
   "source": [
    "#### box plot - RNA Endo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2083cee",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"shareseq_mouse_skin_rare_cell_2_f1_try2_boxplot_rna_Endo.pdf\"\n",
    "subplot_idx = list(c(\"noMiss_Endo\"),\n",
    "                   c(\"multiMissEndo_Endo\"),\n",
    "                   c(\"rnaOnlyEndo_Endo\")\n",
    "                  )\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"scRNA\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    # mean after grouping by metric and method\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "\n",
    "# plot only the unpaired and multiome-guided categories\n",
    "res_plot_3 <- res_plot_3 %>% filter(type %in% c(\"Unpaired\",\"Multiome-guided\"))\n",
    "\n",
    "x_interval = c(3000) \n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=value,group=method,col=method_name,shape=type)) +\n",
    "            geom_boxplot(lwd=1) +\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            scale_color_manual(values=hex[method_names]) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_blank(),\n",
    "                  strip.text.y = element_text(size = 50,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 30,angle = 90),\n",
    "                  axis.text.y = element_text(size = 50),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(2, \"lines\"),panel.spacing.x = unit(1, \"lines\"),\n",
    "                  panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=25,height=8)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "591ae0ae",
   "metadata": {},
   "source": [
    "#### box plot - ATAC HS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b928049",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"shareseq_mouse_skin_rare_cell_2_f1_try2_boxplot_atac_HS.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"noMiss_HS\"),\n",
    "                   c(\"multiMissHS_HS\"),\n",
    "                   c(\"atacOnlyHS_HS\")\n",
    "                  )\n",
    "\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"snATAC\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "# plot only the unpaired and multiome-guided categories\n",
    "res_plot_3 <- res_plot_3 %>% filter(type %in% c(\"Unpaired\",\"Multiome-guided\"))\n",
    "\n",
    "x_interval = c(3000) \n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=value,group=method,col=method_name,shape=type)) +\n",
    "            geom_boxplot(lwd=1) +\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            scale_color_manual(values=hex[method_names]) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_blank(),\n",
    "                  strip.text.y = element_text(size = 50,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 30,angle = 90),\n",
    "                  axis.text.y = element_text(size = 50),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(2, \"lines\"),panel.spacing.x = unit(1, \"lines\"),\n",
    "                  panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=25,height=8)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb55d918",
   "metadata": {},
   "source": [
    "#### box plot - ATAC Endo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82cfc503",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"shareseq_mouse_skin_rare_cell_2_f1_try2_boxplot_atac_Endo.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"noMiss_Endo\"),\n",
    "                   c(\"multiMissEndo_Endo\"),\n",
    "                   c(\"atacOnlyEndo_Endo\")\n",
    "                  )\n",
    "# add summarized values: unpaired method - summarize across all var (despite 1000/3000/8000 multiome cells); for others, summarize across repeats\n",
    "\n",
    "res_plot_3 <- res_plot_2 %>% \n",
    "    filter(dataset==\"snATAC\") %>%\n",
    "    filter(key==\"F1\") %>%\n",
    "    group_by(metric_name,var,method) %>% \n",
    "    # mean after grouping by metric, method, and var\n",
    "    dplyr::mutate(mean_var = mean(value),\n",
    "                  sd_var = sd(value)) %>%\n",
    "    group_by(metric_name,method) %>%\n",
    "    dplyr::mutate(mean_method = mean(value),\n",
    "                  sd_method = sd(value)) %>%\n",
    "    mutate(mean_f = ifelse(type==\"Unpaired\",mean_method,mean_var),\n",
    "           sd_f = ifelse(type==\"Unpaired\",0,sd_var))\n",
    "\n",
    "# plot only the unpaired and multiome-guided categories\n",
    "res_plot_3 <- res_plot_3 %>% filter(type %in% c(\"Unpaired\",\"Multiome-guided\"))\n",
    "\n",
    "x_interval = c(3000) \n",
    "ylim_list = rep(list(c(0,1)), length(unlist(subplot_idx)))\n",
    "nbreaks = rep(list(8),length(unlist(subplot_idx)))\n",
    "counter=1\n",
    "\n",
    "subplots = list()\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_3 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=value,group=method,col=method_name,shape=type)) +\n",
    "            geom_boxplot(lwd=1) +\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            scale_color_manual(values=hex[method_names]) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_blank(),\n",
    "                  strip.text.y = element_text(size = 50,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 30,angle = 90),\n",
    "                  axis.text.y = element_text(size = 50),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(2, \"lines\"),panel.spacing.x = unit(1, \"lines\"),\n",
    "                  panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "    if(plot_ref_accuracy){\n",
    "        if(\"F1\" %in% subplot_idx[[i]]){\n",
    "            subplots[[i]] = subplots[[i]] +\n",
    "                geom_hline(data = truth_df, aes(yintercept = Z),linetype =2,color=\"red\")\n",
    "        }\n",
    "    }\n",
    "\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=25,height=8)                 \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f4d82bc",
   "metadata": {},
   "source": [
    "# HPAP integration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68d8bdd2",
   "metadata": {},
   "source": [
    "## Load results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe6b7ca0",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path =  \"dataset/hpap/real_data/\"\n",
    "folder_path <- \"figures/metric_plots/\"\n",
    "dir.create(folder_path)\n",
    "\n",
    "res_all <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3_metric.csv\")\n",
    "\n",
    "res_all_batch <- plot_result(dir_path,res_folder_filter = c(\"results\"),runtime_exists=T,prediction_exists=T,metric_file = \"ct3_metric_sample_batch.csv\")\n",
    "\n",
    "res_all_comb <- left_join(res_all,res_all_batch[2:10])\n",
    "res_all_comb$log2_runtime <- log2(res_all_comb$runtime)\n",
    "\n",
    "table(is.na(res_all_comb$batch_b_saw))\n",
    "\n",
    "\n",
    "# add _2 to method name and add type column\n",
    "multiome_guided_method_list<-c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat4int\")\n",
    "unpaired_method_list<-c(\"seurat3_2\",\"rfigr_2\",\"rbindsc_2\",\"rliger_2\",\"glue_2\")\n",
    "res_plot <- res_all_comb %>% \n",
    "    mutate(method_name = ifelse(method %in% multiome_guided_method_list,method,paste0(method,\"_2\")))%>% \n",
    "    # assign method type\n",
    "    mutate(type = ifelse(method_name %in% unpaired_method_list,\"Multiome-guided\",\"Unpaired_multiome_splitted\"))\n",
    "\n",
    "metric_sel = c(\"ARI\",\"NMI\",\"ct_aws\",\"clisi\",\"b_saw\",\"kbet\",\"gconn\",\"batch_b_saw\",\"batch_kbet\",\n",
    "               \"log2_runtime\")\n",
    "\n",
    "res_plot = res_plot %>% filter(ct_ref ==\"ct3_metric.csv\")\n",
    "table(res_plot$method_name)\n",
    "\n",
    "method_names = c(\"MultiVI\",\"Seurat v4\", \"Cobolt\",\"scMoMaT\",\"Seurat v4 integrate\",\"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\",\n",
    "                 \"Seurat v3\",\"BindSC\",\"FigR\", \"LIGER\",\"GLUE\")\n",
    "row_order = c(\"multivi\",\"seurat4\",\"cobolt\",\"scmomat\",\"seurat4int\",\"seurat3\",\"rbindsc\",\"rfigr\",\"rliger\",\"glue\",\n",
    "              \"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\")\n",
    "\n",
    "names(method_names)  = row_order\n",
    "\n",
    "method_type = c(rep(\"Multiome-guided\",5),rep(\"Unpaired\",5),rep(\"Unpaired with \\n Multiome-splitted\",5))\n",
    "\n",
    "names(method_type) <- row_order\n",
    "\n",
    "title_list = c(\"ARI\",\"NMI\",\"Cell type ASW\",\"cLisi\", \"Batch ASW\",\"kBET\",\"Graph \\n connectivity\", \n",
    "               \"Donor ASW\",\"Donor kBET\",\"log2(Runtime)\")\n",
    "names(title_list)=metric_sel\n",
    "\n",
    "\n",
    "res_plot_2 <- res_plot %>% \n",
    "    select(ARI:gconn,percent_recovered_50kb:log2_runtime,method_name,var) %>%\n",
    "    gather(key = \"metric_type\",\n",
    "               value = \"value\", -method_name,-var,)%>%\n",
    "    mutate(method = method_name) %>%\n",
    "    mutate(type=factor(method_type[method],levels = c(\"Unpaired\",\"Unpaired with \\n Multiome-splitted\",\"Multiome-guided\"))) %>% \n",
    "    mutate(method_name= factor(method_names[method],levels=unique(method_names))) %>% \n",
    "    mutate(metric_name = factor(title_list[metric_type],levels = title_list))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aae0a305",
   "metadata": {},
   "source": [
    "## Summary plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb826b95",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name <- \"hpap_fair_long_scmomatCorrected.pdf\"\n",
    "\n",
    "subplot_idx = list(c(\"Batch ASW\",\"kBET\"),\n",
    "                   c(\"Donor ASW\",\"Donor kBET\"),\n",
    "                   c('log2(Runtime)'))\n",
    "subplot_idx = as.list(unlist(subplot_idx))\n",
    "\n",
    "\n",
    "ylim_list = list(c(0.5,1),c(0,0.5),c(0.5,1),c(0,0.7),c(8,16))\n",
    "\n",
    "nbreaks = list(5,8,8,9,8)\n",
    "\n",
    "\n",
    "res_plot_2$method <- factor(res_plot_2$method,\n",
    "                           levels=c(\"seurat3_2\",\"rbindsc_2\",\"rfigr_2\",\"rliger_2\",\"glue_2\",\n",
    "                                   \"cobolt\",\"multivi\",\"seurat4\",\"seurat4int\",\"scmomat\"))\n",
    "\n",
    "x_interval = unique(res_plot_2$var)\n",
    "subplots = list()\n",
    "counter = 1\n",
    "for(i in 1:length(subplot_idx)){\n",
    "    print(i)\n",
    "    subplots[[i]] = res_plot_2 %>%  \n",
    "            filter(metric_name %in% subplot_idx[[i]]) %>%  \n",
    "            ggplot(.,aes(x=var,y=value,group=method,col=method_name,shape=type)) +\n",
    "            geom_point(size=3) +\n",
    "            scale_x_continuous(breaks=x_interval) +\n",
    "            facet_grid(cols = vars(type),rows=vars(metric_name),scales = \"free_y\", switch=\"y\") +\n",
    "            scale_color_manual(values=hex[method_names])+\n",
    "            theme_light()+\n",
    "            theme(strip.text.x = element_blank(),\n",
    "                  strip.text.y = element_text(size = 35,angle=90,colour = \"black\"),\n",
    "                  axis.text.x  = element_text(size = 30,angle = 90),axis.text.y = element_text(size = 30),\n",
    "                  axis.title.x = element_blank(),axis.title.y = element_blank(),\n",
    "                  strip.background = element_blank(), strip.placement = \"outside\",\n",
    "                  panel.spacing.y = unit(2, \"lines\"),panel.spacing.x = unit(1, \"lines\"),\n",
    "                 panel.border=element_blank(), axis.line=element_line())+\n",
    "    ggh4x::facetted_pos_scales(\n",
    "        y = list(\n",
    "            scale_y_continuous(limits = ylim_list[[counter]],n.breaks = nbreaks[[counter]])\n",
    "        )\n",
    "    )\n",
    "    counter = counter + 1\n",
    "}\n",
    "subplots_no_bdg = lapply(subplots,function(x){x+theme(legend.position = \"none\")})\n",
    "pcomb <- grid.arrange(grobs=subplots_no_bdg,nrow=1)\n",
    "\n",
    "\n",
    "# function to extract legend from plot\n",
    "get_only_legend <- function(plot) {\n",
    "  plot_table <- ggplot_gtable(ggplot_build(plot))\n",
    "  legend_plot <- which(sapply(plot_table$grobs, function(x){x$name}) == \"guide-box\")\n",
    "  legend <- plot_table$grobs[[legend_plot]]\n",
    "  return(legend)\n",
    "}\n",
    "                            \n",
    "# extract legend from plot1 using above function\n",
    "legend <- get_only_legend(subplots[[1]])   \n",
    "                                                          \n",
    "pdf(file.path(folder_path,file_name),width=30,height=8)             \n",
    "grid.arrange(pcomb, legend, ncol = 2, widths = c(10, 1))\n",
    "dev.off()\n",
    "\n",
    "                              "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "r_py-r",
   "language": "R",
   "name": "r_py-r"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "322.188px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}