umap_plot_generation_public.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cc09bfca",
   "metadata": {},
   "source": [
    "# Readme\n",
    "\n",
    "This script contains codes for generating the UMAP plots in the manuscript titled 'Benchmarking algorithms for joint integration of unpaired and paired single-cell RNA-seq and ATAC-seq data'. Most of these figures are in the Additional File 1. \n",
    "\n",
    "Please run steps described in evaluate_vary_situations_public.ipynb before this script. \n",
    "\n",
    "The same procedure is applied to all scenarios to generate the UMAP plot. One is demonstrated here as an example\n",
    "\n",
    "There are two separate parts. First part is run using a python kernal in which the umap coordinates and cell labels are saved into csv files. Second part is run using a R kernel that reads in the csv files and generate the UMAP plots. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd403e8b",
   "metadata": {},
   "source": [
    "# Part 1 - Python"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44d99303",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baec6d06",
   "metadata": {},
   "outputs": [],
   "source": [
    "# use scib2 python kernel \n",
    "import sys \n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pickle\n",
    "import utils_eval\n",
    "import matplotlib.pyplot as plt\n",
    "from anndata import AnnData\n",
    "import scanpy as sc\n",
    "import pickle\n",
    "import anndata as ad \n",
    "os.environ['R_HOME'] = '/home/myylee/anaconda3/envs/scib2/lib/R/'\n",
    "import scib\n",
    "from copy import deepcopy \n",
    "import tempfile\n",
    "import shutil\n",
    "\n",
    "import random\n",
    "\n",
    "# Function\n",
    "# Generate UMAP plot: one with single modality cells, the other one with all cells. Will be different depending on what cells should be evaluated \n",
    "# inputs: \n",
    "# @folder_dir: folder_dir\n",
    "# @dir_path: dir_path\n",
    "# @output_folder: output_folder \n",
    "# @cond_key: cond_key\n",
    "# @var_list: var_list\n",
    "# @rep_list: rep_list \n",
    "# @method_list: method_list\n",
    "# @res_list: res_list\n",
    "# @ct_ref: ct_ref\n",
    "# @nclust: nclust=None\n",
    "# @plot_out_dir: plot_out_dir \n",
    "\n",
    "# single modality cells only \n",
    "def umap_plot(dir_path,ct_ref,nclust_target,plot_out_dir,folder_dir,output_folder,cond_key,var_list,rep_list,res_list,\n",
    "              with_multi_folder=\"results_single_same_cell_number\",legend_loc=None):\n",
    "    pd_comb_list = [None]*len(res_list)\n",
    "    os.makedirs(os.path.join(\"figures\",\"umap\",plot_out_dir),exist_ok=True)\n",
    "    for i in range(len(var_list)):\n",
    "        print(i)\n",
    "        nclust = nclust_target\n",
    "        if \"seurat4\" in res_list[i] or \"liger\" in res_list[i]:\n",
    "            print(\"not clustering\")\n",
    "            nclust = None\n",
    "\n",
    "        random.seed(1234)\n",
    "\n",
    "        path = os.path.join(folder_dir,dir_path,\n",
    "                            \"{}{}_{}\".format(cond_key,var_list[i],rep_list[i]),\n",
    "                            output_folder[i],\n",
    "                             res_list[i])\n",
    "        res_df = pd.read_csv(path,index_col=0)\n",
    "\n",
    "        col_sel = [col for col in res_df.columns if 'latent' in col]\n",
    "        sub_df = res_df.loc[: , col_sel]\n",
    "\n",
    "        cell_type = pd.read_csv(ct_ref,index_col=0)\n",
    "        cell_type.columns = ['cell_type']\n",
    "        # Assuming everything before '_' is not useful \n",
    "        cell_bc = [sub.split('_')[len(sub.split('_'))-1] for sub in list(res_df.index)]\n",
    "        res_df['cell_type'] = list(cell_type.loc[cell_bc,'cell_type'])\n",
    "\n",
    "        adata_plot = AnnData(res_df.loc[: , col_sel].to_numpy(),\n",
    "                             obs=res_df.loc[:,['dataset','cell_type']],\n",
    "                             dtype=np.float32)\n",
    "        sc.pp.neighbors(adata_plot, n_neighbors=15, use_rep=\"X\")\n",
    "        umap_col_sel = [col for col in res_df.columns if 'umap' in col]\n",
    "        if len(umap_col_sel) > 0: \n",
    "            print(\"using saved UMAP loading\")\n",
    "            adata_plot.obsm['umap'] = res_df.loc[: , umap_col_sel].to_numpy()\n",
    "        else:\n",
    "            print(\"calculating UMAP loading\")\n",
    "            sc.tl.umap(adata_plot, min_dist=0.2)\n",
    "\n",
    "        if nclust is not None:\n",
    "            print(\"clustering with Louvain\")\n",
    "            resolution = utils_eval.find_resolution_louvain(adata_plot,nclust)\n",
    "            sc.tl.louvain(adata_plot, resolution = resolution, random_state = 0)\n",
    "            adata_plot.obs['predicted_ct'] = adata_plot.obs['louvain']\n",
    "        else:\n",
    "            print(\"loading saved clustering result\")\n",
    "            adata_plot.obs['predicted_ct'] = res_df.loc[:,['predicted_ct']]\n",
    "            adata_plot.obs['predicted_ct'] = adata_plot.obs['predicted_ct'].astype('category')\n",
    "\n",
    "        adata_plot.obsm['embed'] = adata_plot.X\n",
    "        adata_plot.obsm['X_emb'] = adata_plot.X\n",
    "\n",
    "        add_str = \"\"\n",
    "        if output_folder[i] == with_multi_folder:\n",
    "            add_str = \"2\"\n",
    "        plot_name = os.path.join(plot_out_dir,\"{}{}_{}_{}\".format(cond_key,var_list[i],rep_list[i],\n",
    "                                                                    res_list[i].replace(\"/\", \"-\").replace(\".csv\", add_str+\".svg\")))\n",
    "        plot_name2 = os.path.join(plot_out_dir,\"{}{}_{}_{}\".format(cond_key,var_list[i],rep_list[i],\n",
    "                                                                     res_list[i].replace(\"/\", \"-\").replace(\".csv\", \"_with_multiome\"+add_str+\".svg\")))\n",
    "        idx = list(adata_plot.obs['dataset'].isin([\"scRNA\",\"snATAC\"]))\n",
    "        adata_plot_sel = deepcopy(adata_plot)[idx,:]\n",
    "        adata_plot_sel.obs = deepcopy(adata_plot.obs).loc[idx,:]\n",
    "        print(adata_plot_sel)\n",
    "        if len(umap_col_sel) > 0: \n",
    "            sc.pl.umap(adata_plot_sel, \n",
    "                       color=['cell_type','dataset','predicted_ct'],\n",
    "                       use_raw =False,\n",
    "                       layer='umap',\n",
    "                       wspace=0.2,\n",
    "                       save=\"/\"+plot_name,\n",
    "                       legend_fontsize=\"x-small\",\n",
    "                       legend_loc=legend_loc,\n",
    "                      size = 10)\n",
    "            pd_umap = pd.DataFrame(adata_plot_sel.obsm['umap'])\n",
    "        else:\n",
    "            sc.pl.umap(adata_plot_sel, \n",
    "                       color=['cell_type','dataset','predicted_ct'],\n",
    "                       wspace=0.2,\n",
    "                       save=\"/\"+plot_name,\n",
    "                       legend_fontsize=\"x-small\",\n",
    "                       legend_loc=legend_loc,\n",
    "                      size = 10)\n",
    "            pd_umap = pd.DataFrame(adata_plot_sel.obsm['X_umap'])\n",
    "\n",
    "        # if multiome dataset is used, plot all cells as well in a different UMAP plot\n",
    "        if adata_plot_sel.shape[0] != adata_plot.shape[0]:\n",
    "            if len(umap_col_sel) > 0: \n",
    "                sc.pl.umap(adata_plot, \n",
    "                           color=['cell_type','dataset','predicted_ct'],\n",
    "                           use_raw =False,\n",
    "                           layer='umap',\n",
    "                           wspace=0.2,\n",
    "                           save=\"/\"+plot_name2,\n",
    "                       legend_fontsize=\"x-small\",\n",
    "                       legend_loc=legend_loc,\n",
    "                      size = 10)\n",
    "            else:\n",
    "                sc.pl.umap(adata_plot, \n",
    "                           color=['cell_type','dataset','predicted_ct'],\n",
    "                           wspace=0.2,\n",
    "                           save=\"/\"+plot_name2,\n",
    "                       legend_fontsize=\"x-small\",\n",
    "                       legend_loc=legend_loc,\n",
    "                      size = 10)\n",
    "        \n",
    "        pd_umap.index = adata_plot_sel.obs.index.tolist()\n",
    "        method_name = res_list[i].split('/', 1) [0]\n",
    "        if output_folder[i] in with_multi_folder:\n",
    "            method_name =method_name+\"_2\"\n",
    "        pd_umap['key'] = \"{}{}_{}_{}\".format(cond_key,var_list[i],rep_list[i],method_name)\n",
    "        pd_comb_list[i] = pd.concat([pd_umap,adata_plot_sel.obs], axis=1)\n",
    "\n",
    "    pd_comb_all = pd.concat(pd_comb_list, axis=0)\n",
    "    return(pd_comb_all)\n",
    "\n",
    "\n",
    "# Function\n",
    "# not just the single-modality cells, plot the multiome cells as well, and return the embedding.\n",
    "# half multiomeRNA + half multiomeATAC for the unpaired multiome-splitted methods\n",
    "def umap_plot_all_cells(dir_path,ct_ref,nclust_target,plot_out_dir,folder_dir,output_folder,cond_key,var_list,rep_list,res_list,\n",
    "              with_multi_folder=\"results_single_same_cell_number\",legend_loc=None):\n",
    "    from numpy.random import choice\n",
    "    subsample_size = lambda vector, size: choice(vector, size = size, replace = False) \n",
    "\n",
    "    pd_comb_list = [None]*len(res_list)\n",
    "    os.makedirs(os.path.join(\"figures\",\"umap\",plot_out_dir),exist_ok=True)\n",
    "    for i in range(len(var_list)):\n",
    "        print(i)\n",
    "        nclust = nclust_target\n",
    "        if \"seurat4\" in res_list[i] or \"liger\" in res_list[i]:\n",
    "            print(\"not clustering\")\n",
    "            nclust = None\n",
    "\n",
    "        random.seed(1234)\n",
    "\n",
    "        path = os.path.join(folder_dir,dir_path,\n",
    "                            \"{}{}_{}\".format(cond_key,var_list[i],rep_list[i]),\n",
    "                            output_folder[i],\n",
    "                             res_list[i])\n",
    "        res_df = pd.read_csv(path,index_col=0)\n",
    "\n",
    "        col_sel = [col for col in res_df.columns if 'latent' in col]\n",
    "        sub_df = res_df.loc[: , col_sel]\n",
    "\n",
    "        cell_type = pd.read_csv(ct_ref,index_col=0)\n",
    "        cell_type.columns = ['cell_type']\n",
    "        # Assuming everything before '_' is not useful \n",
    "        cell_bc = [sub.split('_')[len(sub.split('_'))-1] for sub in list(res_df.index)]\n",
    "        res_df['cell_type'] = list(cell_type.loc[cell_bc,'cell_type'])\n",
    "\n",
    "        adata_plot = AnnData(res_df.loc[: , col_sel].to_numpy(),\n",
    "                             obs=res_df.loc[:,['dataset','cell_type']],\n",
    "                             dtype=np.float32)\n",
    "        sc.pp.neighbors(adata_plot, n_neighbors=15, use_rep=\"X\")\n",
    "        umap_col_sel = [col for col in res_df.columns if 'umap' in col]\n",
    "        if len(umap_col_sel) > 0: \n",
    "            print(\"using saved UMAP loading\")\n",
    "            adata_plot.obsm['umap'] = res_df.loc[: , umap_col_sel].to_numpy()\n",
    "        else:\n",
    "            print(\"calculating UMAP loading\")\n",
    "            sc.tl.umap(adata_plot, min_dist=0.2)\n",
    "\n",
    "        if nclust is not None:\n",
    "            print(\"clustering with Louvain\")\n",
    "            resolution = utils_eval.find_resolution_louvain(adata_plot,nclust)\n",
    "            sc.tl.louvain(adata_plot, resolution = resolution, random_state = 0)\n",
    "            adata_plot.obs['predicted_ct'] = adata_plot.obs['louvain']\n",
    "        else:\n",
    "            print(\"loading saved clustering result\")\n",
    "            adata_plot.obs['predicted_ct'] = res_df.loc[:,['predicted_ct']]\n",
    "            adata_plot.obs['predicted_ct'] = adata_plot.obs['predicted_ct'].astype('category')\n",
    "\n",
    "        adata_plot.obsm['embed'] = adata_plot.X\n",
    "        adata_plot.obsm['X_emb'] = adata_plot.X\n",
    "\n",
    "        add_str = \"\"\n",
    "        if output_folder[i] == with_multi_folder:\n",
    "            add_str = \"2\"\n",
    "        plot_name2 = os.path.join(plot_out_dir,\"{}{}_{}_{}\".format(cond_key,var_list[i],rep_list[i],\n",
    "                                                                     res_list[i].replace(\"/\", \"-\").replace(\".csv\", \"_with_multiome\"+add_str+\".svg\")))\n",
    "        \n",
    "\n",
    "        idx = list(adata_plot.obs['dataset'].isin([\"Multiome-RNA\",\"Multiome-ATAC\",\"multiomeRNA\",\"multiomeATAC\"]))\n",
    "        adata_plot.obs['dataset2'] = adata_plot.obs['dataset'].tolist()\n",
    "        adata_plot.obs['dataset2'][idx] = 'Multiome'\n",
    "        if len(idx) > 0:\n",
    "            # select half of the rna cells and half of the atac cells, instead of randomly selecting from the giant list, which does not gaureentee always get half-half\n",
    "            # get the barcode of every cell \n",
    "            adata_plot.obs['bc2'] = [s.split(\"_\") [len(s.split(\"_\"))-1] for s in adata_plot.obs_names]\n",
    "            # get the id for RNA \n",
    "            idx_rna = list(adata_plot.obs['dataset'].isin([\"Multiome-RNA\",\"multiomeRNA\"]))\n",
    "            # get id for ATAC \n",
    "            idx_atac = list(adata_plot.obs['dataset'].isin([\"Multiome-ATAC\",\"multiomeATAC\"]))\n",
    "            # create barcode column with the right prefix\n",
    "            adata_plot.obs['bc3'] = adata_plot.obs['bc2']\n",
    "            adata_plot.obs['bc3'][idx_rna] = \"prna_\"+adata_plot.obs['bc3'][idx_rna] \n",
    "            adata_plot.obs['bc3'][idx_atac] = \"patac_\"+adata_plot.obs['bc3'][idx_atac]\n",
    "            # select bc, get \n",
    "            npaired = round(sum(idx)/2)\n",
    "            bc_select = (['prna_']*round(npaired/2) + ['patac_']*(npaired-round(npaired/2))) + adata_plot.obs['bc2'][idx].sample(frac=1).unique()\n",
    "            idx_multi_kept = adata_plot.obs['bc3'].isin(bc_select)\n",
    "            idx_multi_kept = [i for i, x in enumerate(idx_multi_kept) if x]\n",
    "            print(len(idx_multi_kept))\n",
    "            idx_kept = list(~adata_plot.obs['dataset'].isin([\"Multiome-RNA\",\"Multiome-ATAC\",\"multiomeRNA\",\"multiomeATAC\"]))\n",
    "            idx_kept = [i for i, x in enumerate(idx_kept) if x]\n",
    "            idx_f = idx_kept + idx_multi_kept\n",
    "            adata_plot_sel = deepcopy(adata_plot)[idx_f,:]\n",
    "            adata_plot_sel.obs = deepcopy(adata_plot.obs).iloc[idx_f,:]\n",
    "        else:\n",
    "            adata_plot_sel = adata_plot\n",
    "\n",
    "        if len(umap_col_sel) > 0: \n",
    "            sc.pl.umap(adata_plot_sel, \n",
    "                       color=['cell_type','dataset','predicted_ct'],\n",
    "                       use_raw =False,\n",
    "                       layer='umap',\n",
    "                       wspace=0.4,\n",
    "                       save=\"/\"+plot_name2,\n",
    "                       legend_loc=legend_loc,\n",
    "                       legend_fontsize=\"xx-small\")\n",
    "            pd_umap = pd.DataFrame(adata_plot_sel.obsm['umap'])\n",
    "        else:\n",
    "            sc.pl.umap(adata_plot_sel, \n",
    "                       color=['cell_type','dataset','predicted_ct'],\n",
    "                       wspace=0.4,\n",
    "                       save=\"/\"+plot_name2,\n",
    "                       legend_loc=legend_loc,\n",
    "                       legend_fontsize=\"xx-small\")\n",
    "            pd_umap = pd.DataFrame(adata_plot_sel.obsm['X_umap'])\n",
    "\n",
    "        pd_umap.index = adata_plot_sel.obs.index.tolist()\n",
    "        method_name = res_list[i].split('/', 1) [0]\n",
    "        if output_folder[i] in with_multi_folder:\n",
    "            method_name =method_name+\"_2\"\n",
    "        pd_umap['key'] = \"{}{}_{}_{}\".format(cond_key,var_list[i],rep_list[i],method_name)\n",
    "        pd_comb_list[i] = pd.concat([pd_umap,adata_plot_sel.obs], axis=1)\n",
    "\n",
    "    pd_comb_all = pd.concat(pd_comb_list, axis=0)\n",
    "    return(pd_comb_all)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d371d5c",
   "metadata": {},
   "source": [
    "## Applying on SHARE-seq mouse skin dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9863cb6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# SHARE-seq mouse skin\n",
    "# working_dir \n",
    "folder_dir=\"/home/myylee/scmint/methods_eval/\"\n",
    "output_folder = ([\"results_single_mod\"]*5+[\"results_single_same_cell_number\"]*5+[\"results_single_mod\"]*4)*2\n",
    "cond_key = \"nmulti\"\n",
    "var_list = [15000]*14+[5000]*14\n",
    "rep_list = [2]*14+[1]*14\n",
    "res_list = ([\"seurat3/seurat3_result.csv\",\"rliger/rliger_result.csv\",\"rfigr/rfigr_result.csv\",\"rbindsc/rbindsc_result.csv\",\"glue/glue_result.csv\", # unpaired \n",
    "            \"seurat3/seurat3_result.csv\",\"rliger/rliger_result.csv\",\"rfigr/rfigr_result.csv\",\"rbindsc/rbindsc_result.csv\",\"glue/glue_result.csv\", # unpaired with multiome\n",
    "           \"seurat4/seurat4_result.csv\",\"multivi/multivi_result.csv\",\"cobolt/cobolt_result.csv\" ,\"scmomat/scmomat_result.csv\"#multiome-guided\n",
    "           ])*2\n",
    "\n",
    "print(res_list)\n",
    "\n",
    "#----- shareseq mouse skin 22 ct -----#\n",
    "dir_path = \"/project/mingyaolpc/myylee/scmint/methods_eval/dataset/mouse_skin/multiome_ncells_pmat/\"\n",
    "ct_ref = \"/project/mingyaolpc/myylee/scmint/methods_eval/dataset/mouse_skin/mouse_skin_shareseq_bc_ct3.csv\"\n",
    "nclust_target = 22\n",
    "plot_out_dir = \"mouse_skin_22ct_all/\"\n",
    "\n",
    "\n",
    "\n",
    "df_to_plot = umap_plot(dir_path,\n",
    "          ct_ref,\n",
    "          nclust_target,\n",
    "          plot_out_dir,\n",
    "          folder_dir,\n",
    "          output_folder,\n",
    "          cond_key,\n",
    "          var_list,\n",
    "          rep_list,\n",
    "          res_list,\n",
    "          with_multi_folder=\"results_single_same_cell_number\",\n",
    "          legend_loc=None)\n",
    "df_to_plot.to_csv(os.path.join(folder_dir,\"figures/umap/\",plot_out_dir,\"df_to_plot.csv\"))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e73b77e4",
   "metadata": {},
   "source": [
    "# Part 2 - R "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f0304ce",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fb4620f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ----- Function ----- # \n",
    "# load saved UMAP loading and metadata csv from umap_plot function in evaluate_vary_situations_all.ipynb and organize umap plota into one plot in R  \n",
    "library(gridExtra)        \n",
    "generate_umap_plots_raster <- function(df_plot,\n",
    "                                       method_names,\n",
    "                                       tops,\n",
    "                                       column_names=c(\"cell_type\",\"predicted_ct\",\"dataset\"),\n",
    "                                       columns_to_plot= c(\"cell_type\",\"predicted_ct\",\"dataset\"),\n",
    "                                       sample_cp = NA,\n",
    "                                       cell_type_cp = NA,\n",
    "                                       pixels = c(128,128),\n",
    "                                       pointsize = 0.1){\n",
    "    require('scattermore')\n",
    "    if((length(column_names)!=length(columns_to_plot))){print(\"columns names ane column to plot list are not the same size. Return without running\"); return()}\n",
    "    p_outer = list()\n",
    "    counter_outer = 1\n",
    "    for(x in columns_to_plot){\n",
    "        df_plot[,x] <- as.factor(df_plot[,x])\n",
    "    }\n",
    "    add_left_title_n = length(unique(df_plot$key))/2\n",
    "    add_top_title_2 = length(unique(df_plot$key))/2+1\n",
    "    for (i in unique(df_plot$key)){\n",
    "        p = list()\n",
    "        counter = 1\n",
    "        for (c_i in 1:length(columns_to_plot)){\n",
    "            c = columns_to_plot[c_i]\n",
    "            p[[counter]]<- df_plot %>% filter(key == i) %>%\n",
    "                ggplot(aes_string(x=\"X0\", y = \"X1\",color=c)) +\n",
    "                geom_scattermore(\n",
    "                    pointsize = pointsize,\n",
    "                    pixels = pixels\n",
    "                ) + \n",
    "                theme(\n",
    "                    panel.grid.major = element_blank(),\n",
    "                    panel.grid.minor = element_blank(),\n",
    "                    strip.background = element_blank(),\n",
    "                    panel.border = element_rect(colour = \"black\", fill = NA),\n",
    "                    panel.background = element_rect(fill='transparent'),\n",
    "                    legend.position=\"None\",\n",
    "                    axis.title.x=element_blank(),\n",
    "                    axis.title.y=element_blank(),\n",
    "                    axis.text.x=element_blank(), #remove x axis labels\n",
    "                    axis.ticks.x=element_blank(), #remove x axis ticks\n",
    "                    axis.text.y=element_blank(),  #remove y axis labels\n",
    "                    axis.ticks.y=element_blank(),  #remove y axis ticks\n",
    "                    plot.title = element_text(hjust = 0.5,size=6)\n",
    "                ) \n",
    "            # use custom color palette for the 'Sample' plot\n",
    "            if(!is.na(sample_cp)&& c ==\"sample\"){\n",
    "                p[[counter]]<- p[[counter]] + scale_color_manual(values=sample_cp)\n",
    "            }\n",
    "            if(!is.na(cell_type_cp)&& c ==\"cell_type\"){\n",
    "                p[[counter]]<- p[[counter]] + scale_color_manual(values=cell_type_cp)\n",
    "            }\n",
    "            \n",
    "            counter = counter +1 \n",
    "        }\n",
    "\n",
    "        p_outer[[counter_outer]]<- grid.arrange(grobs=p,nrow=1)\n",
    "        counter_outer = counter_outer + 1\n",
    "    }\n",
    "    half_n =  length(unique(df_plot$key))/2\n",
    "    pos = list()\n",
    "    for(i in 1:half_n){\n",
    "        pos[[i]] = c(i, i+half_n)\n",
    "    }\n",
    "    lay <- do.call(rbind,pos)\n",
    "    left_text_grobs <- lapply(method_names,function(t){\n",
    "        textGrob(t, gp=gpar(fontsize=6))\n",
    "    })\n",
    "    top_text_grobs <-  lapply(tops,function(t){\n",
    "        textGrob(t, gp=gpar(fontsize=9))\n",
    "    })\n",
    "    p_outer_all<- c(p_outer,top_text_grobs,left_text_grobs)\n",
    "        \n",
    "    # add blank space to top and bottom\n",
    "    # suppose there is two columns \n",
    "    lay <- rbind(c(NA,(length(p_outer)+1),NA,(length(p_outer)+2)),#,ncol(lay)+2),# top text \n",
    "                 cbind(c((length(p_outer_all)-length(method_names)+1):length(p_outer_all)),#rep(NA,nrow(lay)),#left text \n",
    "                       lay[,1], \n",
    "                       rep(NA,nrow(lay)),\n",
    "                       lay[,2]))\n",
    "    p_return <-arrangeGrob(grobs=p_outer_all, layout_matrix = lay,\n",
    "                           widths = c(2,9,1,9),\n",
    "                           height = c(1,rep(8,nrow(lay)-1)))\n",
    "    return(p_return)\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cda1a702",
   "metadata": {},
   "source": [
    "## Applying on SHARE-seq mouse skin dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7109a130",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plot <- read.csv(\"figures/umap/mouse_skin_22ct_all/df_to_plot.csv\")\n",
    "unique(df_plot$key)\n",
    "\n",
    "sceanrio = c(\"nmulti5000_1\",\"nmulti15000_2\")\n",
    "methods = c(\"seurat3\",\"rliger\",\"rfigr\",\"rbindsc\",\"glue\",\n",
    "            \"seurat3_2\",\"rliger_2\",\"rfigr_2\",\"rbindsc_2\",\"glue_2\",\n",
    "            \"seurat4\",\"multivi\",\"cobolt\",\"scmomat\")\n",
    "levels_order = paste0(rep(sceanrio,each=length(methods)),\n",
    "                      \"_\",\n",
    "                      rep(methods,length(sceanrio)))\n",
    "levels_order\n",
    "\n",
    "#create a factor, level it by d2\n",
    "df_plot$key <- factor(df_plot$key, levels=levels_order)\n",
    "\n",
    "#order\n",
    "df_plot <- df_plot[order(df_plot$key),] \n",
    "unique(df_plot$key)\n",
    "\n",
    "# load csv \n",
    "method_names = c(\"Seurat3\\nUnpaired\",\"LIGER\\nUnpaired\",\"FigR\\nUnpaired\",\"BindSC\\nUnpaired\",\"Glue\\nUnpaired\",\n",
    "                \"Seurat3\\nUnpaired\\n(Multiome-split)\",\"LIGER\\nUnpaired\\n(Multiome-split)\",\n",
    "                 \"FigR\\nUnpaired\\n(Multiome-split)\",\"BindSC\\nUnpaired\\n(Multiome-split)\",\n",
    "                 \"Glue\\nUnpaired\\n(Multiome-split)\",\n",
    "                \"Seurat4\\nMultiome-guided\",\"MultiVI\\nMultiome-guided\",\"Cobolt\\nMultiome-guided\",\n",
    "                \"scMoMaT\\nMultiome-guided\")\n",
    "\n",
    "#method_names = c(\"scMoMaT\\nMultiome-guided\",\"Glue\\nUnpaired\\n(Multiome-splitted)\")\n",
    "tops = c(\"n_multiome = 5000 cells\",\"n_multiome = 15000 cells\")\n",
    "column_names=c(\"Cell type\",\"Predicted\",\"Modality\")\n",
    "df_plot$dataset <- factor(df_plot$dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cbc67c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig <- generate_umap_plots_raster(df_plot,method_names,tops,column_names)\n",
    "\n",
    "ggsave(file=\"figures/umap/mouse_skin_22ct_all/df_to_plot_dot2_shape19_rev1_raster.pdf\", fig, width=6, height=11) #saves g\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c0f7933",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": "341.594px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}