plotEmbedding_Starmap_newdata.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "compliant-photograph",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('/home/xinyiz/pamrats')\n",
    "\n",
    "import time\n",
    "import os\n",
    "\n",
    "import scanpy\n",
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "\n",
    "import torch\n",
    "from torch import optim\n",
    "\n",
    "\n",
    "import gae.gae.optimizer as optimizer\n",
    "import gae.gae.model\n",
    "import gae.gae.preprocessing as preprocessing\n",
    "\n",
    "import pickle\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import umap\n",
    "import pandas as pd\n",
    "from sklearn.preprocessing import scale\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "from sklearn.cluster import DBSCAN,MiniBatchKMeans,AgglomerativeClustering\n",
    "from sklearn import metrics\n",
    "\n",
    "import anndata as ad\n",
    "import gc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a48a04ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "olddata=scanpy.read_h5ad('/home/xinyiz/2021-01-13-mAD-test-dataset/2020-12-27-starmap-mAD-raw.h5ad')\n",
    "oldgenes=olddata.var.index.to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f5423685",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "# Settings\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"2\" \n",
    "use_cuda=True\n",
    "seed=3\n",
    "useSavedMaskedEdges=True\n",
    "maskedgeName='knn20_connectivity'\n",
    "\n",
    "ifplot=True\n",
    "ifcluster=False\n",
    "\n",
    "inverseAct='leakyRelu'\n",
    "# inverseAct=None\n",
    "plottype='pca'\n",
    "pca=PCA()\n",
    "minCells=15 #min number of cells for analysis\n",
    "# clustermethod=['kmeanbatch']\n",
    "# clustermethod=['leiden','agglomerative','kmeanbatch']\n",
    "clustermethod=['leiden']\n",
    "#umap/leiden clustering parameters\n",
    "n_neighbors=10\n",
    "min_dist=0.25\n",
    "n_pcs=40 #for clustering\n",
    "# resolution=[0.1]\n",
    "resolution=[0.081,0.082,0.083,0.085,0.087]\n",
    "plotepoch=9990\n",
    "savenameAdd=''\n",
    "#DBscan\n",
    "epslist= [6,8,10]\n",
    "min_sampleslist=[15,30,45] \n",
    "#agglomerative\n",
    "nclusterlist=[2,3,4,5,8,10,15]\n",
    "aggMetric=['euclidean']\n",
    "\n",
    "hidden1=6000 #Number of units in hidden layer 1\n",
    "hidden2=6000 #Number of units in hidden layer 2\n",
    "# hidden3=2048\n",
    "# hidden4=2048\n",
    "# hidden5=128\n",
    "fc_dim1=6000\n",
    "# fc_dim2=128\n",
    "# fc_dim3=128\n",
    "# fc_dim4=128\n",
    "# gcn_dim1=2600\n",
    "# clf_hidden=256\n",
    "adv_hidden=128\n",
    "\n",
    "dropout=0.01\n",
    "testNodes=0.1 #fraction of total nodes for testing\n",
    "valNodes=0.05 #fraction of total nodes for validation\n",
    "XreconWeight=20\n",
    "# clfweight=20\n",
    "advWeight=2\n",
    "# randFeatureSubset=None\n",
    "model_str='gcn_vae_xa_e2_d1_dca'\n",
    "clf=None\n",
    "adv=None  #'clf_fc1_eq'  #'clf_fc1_control_eq' #'clf_fc1_control'  #'clf_fc1'\n",
    "protein=None #'nearest' #None #'scaled_binary'\n",
    "plaqueMaskName='PlaqueMask'\n",
    "proteinWeights=0.001\n",
    "adj_decodeName=None #gala or None\n",
    "ridgeL=0.01\n",
    "shareGenePi=True\n",
    "\n",
    "pretrainedAE=None #{'name':'controlphy5XAbin_01_dca','epoch':9990}\n",
    "targetBatch=None\n",
    "plot_sample_X=['logminmax']\n",
    "switchFreq=10\n",
    "standardizeX=False\n",
    "# name='allk20XA_01_dca_noD8' ###rename and retrain D8 C8\n",
    "name='allk20XA_02_dca_over'\n",
    "logsavepath='/mnt/external_ssd/xinyi/log/train_gae_starmap/'+name\n",
    "modelsavepath='/mnt/external_ssd/xinyi/models/train_gae_starmap/'+name\n",
    "# plotsavepath='/mnt/external_ssd/xinyi/plots/train_gae_starmap/'+name\n",
    "plotsavepath='/mnt/external_ssd/xinyi/plots/train_gae_starmap/'+name+'/newold'\n",
    "\n",
    "datadir='/mnt/external_ssd/xinyi/starmap_new/2022-06-14-Xinyi-mAD'\n",
    "\n",
    "#Load data\n",
    "# plot_samples={'disease13':'ADmouse_11346','control13':'ADmouse_11351','disease8':'ADmouse_9723_2','control8':'ADmouse_9707'}\n",
    "plot_samples={'disease13old':'ADmouse_9494','control13old':'ADmouse_9498','disease8old':'ADmouse_9723','control8old':'ADmouse_9735','disease13':'ADmouse_11346','control13':'ADmouse_11351','disease8':'ADmouse_9723_2','control8':'ADmouse_9707'}\n",
    "plot_samples_old={'disease13old':'AD_mouse9494','control13old':'AD_mouse9498','disease8old':'AD_mouse9723','control8old':'AD_mouse9735'}\n",
    "adj_dir=os.path.join(datadir,'a')\n",
    "\n",
    "featureslist={}\n",
    "if 'dca' in model_str:\n",
    "    features_raw_list={}\n",
    "if plot_sample_X[0] in ['corrected','scaled']:\n",
    "    scaleddata=scanpy.read_h5ad(datadir+'/2022-04-06-Hu-AD-stardist-scaled.h5ad')\n",
    "    \n",
    "    for s in plot_samples.keys():\n",
    "        featureslist[s+'X_'+'corrected']=torch.tensor(scaleddata.layers['corrected'][scaleddata.obs['sample']==plot_samples[s]])\n",
    "        featureslist[s+'X_'+'scaled']=torch.tensor(scaleddata.layers['scaled'][scaleddata.obs['sample']==plot_samples[s]])\n",
    "\n",
    "else:\n",
    "    scaleddata=scanpy.read_h5ad(datadir+'/2022-04-06-Hu-AD-stardist-scaled.h5ad') #new model\n",
    "    #old model\n",
    "    scaleddata_all=scanpy.read_h5ad(datadir+'/2022-04-06-Hu-AD-stardist-scaled.h5ad')\n",
    "    _,idx,_=np.intersect1d(np.char.upper(scaleddata_all.var.index.to_numpy().astype(str)),oldgenes,return_indices=True)\n",
    "    scaleddata=scaleddata_all[:,idx]\n",
    "    \n",
    "    for s in plot_samples.keys():\n",
    "        if 'old' in s:\n",
    "            scaleddata_train=olddata.X[olddata.obs['sample']==plot_samples_old[s]]\n",
    "            features_raw_list[s+'X_'+'raw']=torch.tensor(scaleddata_train)\n",
    "        else:\n",
    "            scaleddata_train=scaleddata.layers['raw'][scaleddata.obs['sample']==plot_samples[s]]\n",
    "            features_raw_list[s+'X_'+'raw']=torch.tensor(scaleddata_train)\n",
    "        \n",
    "        if plot_sample_X[0]=='logminmax':\n",
    "            featurelog_train=np.log2(scaleddata_train+1/2)\n",
    "            scaler = MinMaxScaler()\n",
    "            featurelog_train_minmax=np.transpose(scaler.fit_transform(np.transpose(featurelog_train)))\n",
    "            featureslist[s+'X_'+plot_sample_X[0]]=torch.tensor(featurelog_train_minmax)\n",
    "        elif plot_sample_X[0]=='logminmax10':\n",
    "            featurelog_train=np.log2(scaleddata_train+1/2)\n",
    "            scaler = MinMaxScaler(feature_range=(0,10))\n",
    "            featurelog_train_minmax=np.transpose(scaler.fit_transform(np.transpose(featurelog_train)))\n",
    "            featureslist[s+'X_'+plot_sample_X[0]]=torch.tensor(featurelog_train_minmax)\n",
    "\n",
    "    \n",
    "adj_list={}\n",
    "for s in plot_samples.keys():\n",
    "    if not 'old' in s:\n",
    "        adj_list[s]=sp.load_npz(os.path.join(adj_dir,maskedgeName+'_'+plot_samples[s]+'.npz'))\n",
    "savedir=os.path.join('/home/xinyiz/starmap')\n",
    "adj_dir=os.path.join(savedir,'a')\n",
    "adj_list['disease13old']=sp.load_npz(os.path.join(adj_dir,maskedgeName+'_AD_mouse9494.npz'))\n",
    "adj_list['control13old']=sp.load_npz(os.path.join(adj_dir,maskedgeName+'_AD_mouse9498.npz'))\n",
    "adj_list['disease8old']=sp.load_npz(os.path.join(adj_dir,maskedgeName+'_AD_mouse9723.npz'))\n",
    "adj_list['control8old']=sp.load_npz(os.path.join(adj_dir,maskedgeName+'_AD_mouse9735.npz'))\n",
    "\n",
    "adjnormlist={}\n",
    "pos_weightlist={}\n",
    "normlist={}\n",
    "for ai in adj_list.keys():\n",
    "    adjnormlist[ai]=preprocessing.preprocess_graph(adj_list[ai])\n",
    "    \n",
    "    pos_weightlist[ai] = torch.tensor(float(adj_list[ai].shape[0] * adj_list[ai].shape[0] - adj_list[ai].sum()) / adj_list[ai].sum()) #using full unmasked adj\n",
    "    normlist[ai] = adj_list[ai].shape[0] * adj_list[ai].shape[0] / float((adj_list[ai].shape[0] * adj_list[ai].shape[0] - adj_list[ai].sum()) * 2)\n",
    "    \n",
    "    adj_label=adj_list[ai] + sp.eye(adj_list[ai].shape[0])\n",
    "    adj_list[ai]=torch.tensor(adj_label.todense())\n",
    "    \n",
    "\n",
    "\n",
    "if clf:\n",
    "    ct_train=ctlist[training_samples].astype(int)\n",
    "if standardizeX:\n",
    "    features=torch.tensor(scale(features,axis=0, with_mean=True, with_std=True, copy=True))\n",
    "\n",
    "# Set cuda and seed\n",
    "np.random.seed(seed)\n",
    "if use_cuda and (not torch.cuda.is_available()):\n",
    "    print('cuda not available')\n",
    "    use_cuda=False\n",
    "torch.manual_seed(seed)\n",
    "if use_cuda:\n",
    "    torch.cuda.manual_seed(seed)\n",
    "    torch.backends.cudnn.enabled = True\n",
    "\n",
    "\n",
    "plotRecon=''\n",
    "num_features=scaleddata.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "planned-graduation",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set cuda and seed\n",
    "np.random.seed(seed)\n",
    "if use_cuda and (not torch.cuda.is_available()):\n",
    "    print('cuda not available')\n",
    "    use_cuda=False\n",
    "torch.manual_seed(seed)\n",
    "if use_cuda:\n",
    "    torch.cuda.manual_seed(seed)\n",
    "    torch.backends.cudnn.enabled = True\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "grateful-farming",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<All keys matched successfully>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# load model\n",
    "num_nodes,num_features = list(featureslist.values())[0].shape\n",
    "if model_str=='gcn_vae_xa':\n",
    "    model  = gae.gae.model.GCNModelVAE_XA(num_features, hidden1, hidden2,fc_dim1,fc_dim2,fc_dim3,fc_dim4, dropout)\n",
    "elif model_str=='fc1':\n",
    "    model  = gae.gae.model.FCVAE1(num_features, hidden1,dropout)\n",
    "elif model_str == 'gcn_vae_xa_e2_d1':\n",
    "    model  = gae.gae.model.GCNModelVAE_XA_e2_d1(num_features, hidden1,hidden2, dropout)\n",
    "elif model_str == 'gcn_vae_gcnX_inprA':\n",
    "    model = gae.gae.model.GCNModelVAE_gcnX_inprA(num_features, hidden1, hidden2,gcn_dim1, dropout)\n",
    "elif model_str=='fc1_dca':\n",
    "    model = gae.gae.model.FCVAE1_DCA(num_features, hidden1,fc_dim1, dropout)\n",
    "elif model_str=='gcn_vae_xa_e2_d1_dca':\n",
    "    model = gae.gae.model.GCNModelVAE_XA_e2_d1_DCA(num_features, hidden1,hidden2,fc_dim1, dropout)\n",
    "elif model_str=='gcn_vae_xa_e2_d1_dcaFork':\n",
    "    model = gae.gae.model.GCNModelVAE_XA_e2_d1_DCAfork(num_features, hidden1,hidden2,fc_dim1, dropout)\n",
    "elif model_str=='gcn_vae_xa_e2_d1_dcaElemPi':\n",
    "    model = gae.gae.model.GCNModelVAE_XA_e2_d1_DCAelemPi(num_features, hidden1,hidden2,fc_dim1, dropout,shareGenePi)\n",
    "elif model_str=='gcn_vae_xa_e2_d1_dcaConstantDisp':\n",
    "    model = gae.gae.model.GCNModelVAE_XA_e2_d1_DCA_constantDisp(num_features, hidden1,hidden2,fc_dim1, dropout,shareGenePi)\n",
    "else:\n",
    "    print('model not found')\n",
    "model.load_state_dict(torch.load(os.path.join(modelsavepath,str(plotepoch)+'.pt')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "alive-michigan",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed)\n",
    "def plotembeddingbyCT_str(ctlist,savename,excludelist,embedding,savepath,plotname,plotdimx=0,plotdimy=1,savenameAdd=''):\n",
    "    celltypes=np.unique(ctlist)\n",
    "    celltypes_dict={}\n",
    "    idx=0\n",
    "    for ct in celltypes:\n",
    "        celltypes_dict[ct]=idx\n",
    "        idx+=1\n",
    "        \n",
    "    colortest=sns.color_palette(\"husl\", celltypes.size)\n",
    "#     colortest=sns.color_palette(\"husl\", 4)\n",
    "#     np.random.shuffle(colortest)\n",
    "    fig, ax = plt.subplots(dpi=400)\n",
    "    for ct in celltypes:\n",
    "        if ct in excludelist:\n",
    "            continue\n",
    "        idx=(ctlist==ct)\n",
    "        ax.scatter(\n",
    "            embedding[idx, plotdimx],\n",
    "            embedding[idx, plotdimy],\n",
    "            color=colortest[celltypes_dict[ct]],label=ct,s=1.5,alpha=0.5\n",
    "#             color=colortest[int(ct)],label=ct,s=1.5,alpha=0.5\n",
    "            )\n",
    "\n",
    "    plt.gca().set_aspect('equal', 'datalim')\n",
    "    fig.set_figheight(5)\n",
    "    fig.set_figwidth(5)\n",
    "    box = ax.get_position()\n",
    "    ax.set_position([box.x0, box.y0 + box.height * 0.1,\n",
    "                     box.width, box.height * 0.9])\n",
    "    # Put a legend below current axis\n",
    "    ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.05),\n",
    "              fancybox=True, shadow=True, ncol=5)\n",
    "#     ax.legend(ncol=3)\n",
    "    plt.title(plotname+' embedding', fontsize=24)\n",
    "    plt.savefig(os.path.join(savepath,savename+savenameAdd+'.jpg'))\n",
    "#     plt.show()\n",
    "    \n",
    "#     fig.clf()\n",
    "    plt.close('all')\n",
    "    \n",
    "    gc.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "02af9abf",
   "metadata": {},
   "outputs": [],
   "source": [
    "seed=3\n",
    "\n",
    "np.random.seed(seed)\n",
    "def plotembeddingbyCT(ctlist,savename,excludelist,embedding,savepath,plotname,plotdimx=0,plotdimy=1,savenameAdd='',img=None,ncolors=None):\n",
    "    \n",
    "    celltypes=np.unique(ctlist)\n",
    "# #     print(celltypes)\n",
    "#     celltypes_dict={}\n",
    "#     idx=0\n",
    "#     for ct in celltypes:\n",
    "#         celltypes_dict[ct]=idx\n",
    "#         idx+=1\n",
    "    if ncolors is None:\n",
    "        colortest=sns.color_palette(\"husl\", celltypes.size)\n",
    "    else:\n",
    "        colortest=sns.color_palette(\"husl\", ncolors)\n",
    "#     colortest=sns.color_palette(\"husl\", 9)\n",
    "#     np.random.shuffle(colortest)\n",
    "    fig, ax = plt.subplots(dpi=400)\n",
    "    if not img is None:\n",
    "        plt.imshow(img)\n",
    "    for ct in celltypes:\n",
    "        if ct in excludelist:\n",
    "            continue\n",
    "        idx=(ctlist==ct)\n",
    "        if not img is None:\n",
    "            ax.scatter(\n",
    "                embedding[idx, plotdimy],\n",
    "                embedding[idx, plotdimx],\n",
    "#                 color=colortest[celltypes_dict[ct]],label=ct,s=1.5,alpha=0.5\n",
    "                color=colortest[int(ct)],label=ct,s=1.5,alpha=0.5\n",
    "                )\n",
    "        else:\n",
    "            ax.scatter(\n",
    "                embedding[idx, plotdimx],\n",
    "                embedding[idx, plotdimy],\n",
    "#                 color=colortest[celltypes_dict[ct]],label=ct,s=1.5,alpha=0.5\n",
    "                color=colortest[int(ct)],label=ct,s=1.5,alpha=0.5\n",
    "                )\n",
    "\n",
    "    plt.gca().set_aspect('equal', 'datalim')\n",
    "    fig.set_figheight(5)\n",
    "    fig.set_figwidth(5)\n",
    "    box = ax.get_position()\n",
    "    ax.set_position([box.x0, box.y0 + box.height * 0.1,\n",
    "                     box.width, box.height * 0.9])\n",
    "    # Put a legend below current axis\n",
    "#     ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.05),\n",
    "#               fancybox=True, shadow=True, ncol=2,prop={'size': 6})\n",
    "    ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.05),\n",
    "              fancybox=True,ncol=5, shadow=True,prop={'size': 6})\n",
    "#     ax.legend(ncol=3)\n",
    "    plt.title(plotname+' embedding', fontsize=12)\n",
    "    plt.savefig(os.path.join(savepath,savename+savenameAdd+'.jpg'))\n",
    "#     plt.show()\n",
    "    \n",
    "#     fig.clf()\n",
    "    plt.close('all')\n",
    "    \n",
    "    gc.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "daily-converter",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed)\n",
    "def plotembeddingbyCT_contrast(ctlist,savename,excludelist,embedding,savepath,plotname,plotdimx=0,plotdimy=1,savenameAdd='',maxplot=None): \n",
    "    celltypes=np.unique(ctlist)\n",
    "    celltypes_dict={}\n",
    "    idx=0\n",
    "    for ct in celltypes:\n",
    "        celltypes_dict[ct]=idx\n",
    "        idx+=1\n",
    "\n",
    "    colortest=sns.color_palette(\"tab10\")\n",
    "    if not os.path.exists(os.path.join(savepath)):\n",
    "        os.makedirs(savepath)\n",
    "\n",
    "    for ct in celltypes:\n",
    "        if maxplot and int(ct)>maxplot:\n",
    "            continue\n",
    "        fig, ax = plt.subplots()\n",
    "        if ct == 'Unassigned':\n",
    "            continue\n",
    "\n",
    "        idx=(ctlist!=ct)\n",
    "        ax.scatter(\n",
    "            embedding[idx, plotdimx],\n",
    "            embedding[idx, plotdimy],\n",
    "            color=colortest[1],label='others',s=1,alpha=0.5\n",
    "            )\n",
    "\n",
    "        idx=(ctlist==ct)\n",
    "        ax.scatter(\n",
    "            embedding[idx, plotdimx],\n",
    "            embedding[idx, plotdimy],\n",
    "            color=colortest[0],label=ct,s=3,alpha=0.5\n",
    "            )\n",
    "\n",
    "        plt.gca().set_aspect('equal', 'datalim')\n",
    "        fig.set_figheight(10)\n",
    "        fig.set_figwidth(10)\n",
    "        ax.legend()\n",
    "        plt.title(plotname+' embedding', fontsize=24)\n",
    "        plt.gcf().savefig(os.path.join(savepath,savename+'_'+str(ct)+savenameAdd+'.jpg'))\n",
    "#         plt.show()\n",
    "#         nplot+=1\n",
    "        \n",
    "    \n",
    "#         fig.clf()\n",
    "        plt.close('all')\n",
    "        gc.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "essential-treasure",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed)\n",
    "def inverseLeakyRelu(v,slope=0.01):\n",
    "    vnegidx=(v<0)\n",
    "    v[vnegidx]=1/slope*v[vnegidx]\n",
    "    return v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "floral-consortium",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed)\n",
    "def clusterLeiden_single(inArray,n_neighbors,n_pcs,min_dist,resolution,randseed=seed):\n",
    "    n_pcs=np.min([inArray.shape[0]-1,inArray.shape[1]-1,n_pcs])\n",
    "    adata=ad.AnnData(inArray)\n",
    "    scanpy.tl.pca(adata, svd_solver='arpack')\n",
    "    scanpy.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=n_pcs)\n",
    "    scanpy.tl.umap(adata,min_dist=min_dist,random_state=randseed)\n",
    "    scanpy.tl.leiden(adata,resolution=resolution,random_state=randseed)\n",
    "    return adata.obs['leiden'].to_numpy()\n",
    "\n",
    "def clusterLeiden(inArray,n_neighbors,n_pcs,min_dist,resolution,sobj_coord_np,randseed=seed):\n",
    "    for r in resolution:\n",
    "        clusterRes=clusterLeiden_single(inArray,n_neighbors,n_pcs,min_dist,r,randseed=seed)\n",
    "#         print(clusterRes.shape)\n",
    "        savenamecluster='leiden_nn'+str(n_neighbors)+'mdist0'+str(int(min_dist*100))+'n_pcs'+str(n_pcs)+'res'+str(r)+'epoch'+str(plotepoch)\n",
    "        with open(os.path.join(clustersavedir,savenamecluster), 'wb') as output:\n",
    "            pickle.dump(clusterRes, output, pickle.HIGHEST_PROTOCOL)\n",
    "        plotembeddingbyCT(clusterRes,'leiden',[],embedding,savedir,plottype+' of '+s,savenameAdd=savenamecluster)\n",
    "        plotembeddingbyCT_contrast(clusterRes,'leiden',[],embedding,os.path.join(savedir,'contrast'),plottype+' of '+s,savenameAdd=savenamecluster)\n",
    "\n",
    "        plotembeddingbyCT(clusterRes,'leiden_location',[],sobj_coord_np,savedir,'location'+' of '+s,savenameAdd=savenamecluster)\n",
    "        plotembeddingbyCT_contrast(clusterRes,'leiden_location',[],sobj_coord_np,os.path.join(savedir,'contrast'),'location'+' of '+s,savenameAdd=savenamecluster)\n",
    "\n",
    "def clusterLeiden_allsample(embedding,savedir,clustersavedir,inArray,n_neighbors,n_pcs,min_dist,resolution,sobj_coord_np,samplenameList,randseed=seed):\n",
    "    for r in resolution:\n",
    "        clusterRes=clusterLeiden_single(inArray,n_neighbors,n_pcs,min_dist,r,randseed=seed)\n",
    "        savenamecluster='leiden_nn'+str(n_neighbors)+'mdist0'+str(int(min_dist*100))+'n_pcs'+str(n_pcs)+'res'+str(r)+'epoch'+str(plotepoch)\n",
    "        with open(os.path.join(clustersavedir,savenamecluster), 'wb') as output:\n",
    "            pickle.dump(clusterRes, output, pickle.HIGHEST_PROTOCOL)\n",
    "        plotembeddingbyCT(clusterRes,'leiden',[],embedding,savedir,plottype+' of all samples',savenameAdd=savenamecluster)\n",
    "        plotembeddingbyCT_contrast(clusterRes,'leiden',[],embedding,os.path.join(savedir,'contrast'),plottype+' of all samples',savenameAdd=savenamecluster,maxplot=50)\n",
    "\n",
    "        for s in plot_samples.keys():\n",
    "            sidx=(samplenameList==s)\n",
    "            plotembeddingbyCT(clusterRes[sidx],'leiden_location'+s,[],sobj_coord_np[sidx],savedir,'location'+' of '+s,savenameAdd=savenamecluster,ncolors=np.unique(clusterRes).size)\n",
    "            plotembeddingbyCT_contrast(clusterRes[sidx],'leiden_location'+s,[],sobj_coord_np[sidx],os.path.join(savedir,'contrast'),'location'+' of '+s,savenameAdd=savenamecluster,maxplot=50)\n",
    "\n",
    "         "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "active-junior",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed)\n",
    "def clusterDBscan_single(inArray,eps,min_samples,n_pcs):\n",
    "    n_pcs=np.min([inArray.shape[0]-1,inArray.shape[1]-1,n_pcs])\n",
    "    inArray=pca.fit_transform(inArray)\n",
    "    labels = DBSCAN(eps=eps, min_samples=min_samples).fit_predict(inArray[:,:n_pcs])\n",
    "#     db = DBSCAN(eps=eps, min_samples=min_samples).fit(inArray[:,:n_pcs])\n",
    "#     core_samples_mask = np.zeros_like(db.labels_, dtype=bool)\n",
    "#     core_samples_mask[db.core_sample_indices_] = True\n",
    "#     labels = db.labels_\n",
    "    return labels\n",
    "\n",
    "def clusterDBscan(inArray,epsL,min_samplesL,n_pcs,sobj_coord_np):\n",
    "    for eps in epsL:\n",
    "        for min_samples in min_samplesL:\n",
    "            clusterRes=clusterDBscan_single(inArray,eps,min_samples,n_pcs)\n",
    "    #         print(clusterRes.shape)\n",
    "            savenamecluster='dbscan_eps'+str(eps)+'msamples'+str(min_samples)+'n_pcs'+str(n_pcs)+'epoch'+str(plotepoch)\n",
    "            with open(os.path.join(clustersavedir,savenamecluster), 'wb') as output:\n",
    "                pickle.dump(clusterRes, output, pickle.HIGHEST_PROTOCOL)\n",
    "            plotembeddingbyCT(clusterRes,'dbscan',[],embedding,savedir,plottype+' of '+s,savenameAdd=savenamecluster)\n",
    "            plotembeddingbyCT_contrast(clusterRes,'dbscan',[],embedding,os.path.join(savedir,'contrast'),plottype+' of '+s,savenameAdd=savenamecluster)\n",
    "\n",
    "            plotembeddingbyCT(clusterRes,'dbscan_location',[],sobj_coord_np,savedir,'location'+' of '+s,savenameAdd=savenamecluster)\n",
    "            plotembeddingbyCT_contrast(clusterRes,'dbscan_location',[],sobj_coord_np,os.path.join(savedir,'contrast'),'location'+' of '+s,savenameAdd=savenamecluster)\n",
    "\n",
    "def clusterDBscan_allsample(embedding,savedir,clustersavedir,inArray,epsL,min_samplesL,n_pcs,sobj_coord_np,samplenameList):\n",
    "    for eps in epsL:\n",
    "        for min_samples in min_samplesL:\n",
    "            clusterRes=clusterDBscan_single(inArray,eps,min_samples,n_pcs)\n",
    "            savenamecluster='dbscan_eps'+str(eps)+'msample'+str(min_samples)+'n_pcs'+str(n_pcs)+'epoch'+str(plotepoch)\n",
    "            with open(os.path.join(clustersavedir,savenamecluster), 'wb') as output:\n",
    "                pickle.dump(clusterRes, output, pickle.HIGHEST_PROTOCOL)\n",
    "            plotembeddingbyCT(clusterRes,'dbscan',[],embedding,savedir,plottype+' of all samples',savenameAdd=savenamecluster)\n",
    "            plotembeddingbyCT_contrast(clusterRes,'dbscan',[],embedding,os.path.join(savedir,'contrast'),plottype+' of all samples',savenameAdd=savenamecluster,maxplot=50)\n",
    "\n",
    "            for s in plot_samples.keys():\n",
    "                sidx=(samplenameList==s)\n",
    "                plotembeddingbyCT(clusterRes[sidx],'dbscan_location'+s,[],sobj_coord_np[sidx],savedir,'location'+' of '+s,savenameAdd=savenamecluster)\n",
    "                plotembeddingbyCT_contrast(clusterRes[sidx],'dbscan_location'+s,[],sobj_coord_np[sidx],os.path.join(savedir,'contrast'),'location'+' of '+s,savenameAdd=savenamecluster,maxplot=50)\n",
    "\n",
    "         "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "entire-values",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed)\n",
    "def clusterAgg_single(inArray,ncluster,aggmetric,n_pcs):\n",
    "    n_pcs=np.min([inArray.shape[0]-1,inArray.shape[1]-1,n_pcs])\n",
    "    inArray=pca.fit_transform(inArray)\n",
    "    labels = AgglomerativeClustering(n_clusters=ncluster,affinity=aggmetric).fit_predict(inArray[:,:n_pcs])\n",
    "#     labels = agg.labels_\n",
    "    return labels\n",
    "\n",
    "def clusterAgg(inArray,nclusterL,aggmetricL,n_pcs,sobj_coord_np):\n",
    "    for ncluster in nclusterL:\n",
    "        for aggmetric in aggmetricL:\n",
    "            clusterRes=clusterAgg_single(inArray,ncluster,aggmetric,n_pcs)\n",
    "    #         print(clusterRes.shape)\n",
    "            savenamecluster='agg_ncluster'+str(ncluster)+aggmetric+'n_pcs'+str(n_pcs)+'epoch'+str(plotepoch)\n",
    "            with open(os.path.join(clustersavedir,savenamecluster), 'wb') as output:\n",
    "                pickle.dump(clusterRes, output, pickle.HIGHEST_PROTOCOL)\n",
    "            plotembeddingbyCT(clusterRes,'agg',[],embedding,savedir,plottype+' of '+s,savenameAdd=savenamecluster)\n",
    "            plotembeddingbyCT_contrast(clusterRes,'agg',[],embedding,os.path.join(savedir,'contrast'),plottype+' of '+s,savenameAdd=savenamecluster)\n",
    "\n",
    "            plotembeddingbyCT(clusterRes,'agg_location',[],sobj_coord_np,savedir,'location'+' of '+s,savenameAdd=savenamecluster)\n",
    "            plotembeddingbyCT_contrast(clusterRes,'agg_location',[],sobj_coord_np,os.path.join(savedir,'contrast'),'location'+' of '+s,savenameAdd=savenamecluster)\n",
    "\n",
    "def clusterAgg_allsample(embedding,savedir,clustersavedir,inArray,nclusterL,aggmetricL,n_pcs,sobj_coord_np,samplenameList):\n",
    "    for ncluster in nclusterL:\n",
    "        for aggmetric in aggmetricL:\n",
    "            clusterRes=clusterAgg_single(inArray,ncluster,aggmetric,n_pcs)\n",
    "            savenamecluster='agg_ncluster'+str(ncluster)+aggmetric+'n_pcs'+str(n_pcs)+'epoch'+str(plotepoch)\n",
    "            with open(os.path.join(clustersavedir,savenamecluster), 'wb') as output:\n",
    "                pickle.dump(clusterRes, output, pickle.HIGHEST_PROTOCOL)\n",
    "            plotembeddingbyCT(clusterRes,'agg',[],embedding,savedir,plottype+' of all samples',savenameAdd=savenamecluster)\n",
    "            plotembeddingbyCT_contrast(clusterRes,'agg',[],embedding,os.path.join(savedir,'contrast'),plottype+' of all samples',savenameAdd=savenamecluster,maxplot=50)\n",
    "\n",
    "            for s in plot_samples.keys():\n",
    "                sidx=(samplenameList==s)\n",
    "                plotembeddingbyCT(clusterRes[sidx],'agg_location'+s,[],sobj_coord_np[sidx],savedir,'location'+' of '+s,savenameAdd=savenamecluster)\n",
    "                plotembeddingbyCT_contrast(clusterRes[sidx],'agg_location'+s,[],sobj_coord_np[sidx],os.path.join(savedir,'contrast'),'location'+' of '+s,savenameAdd=savenamecluster,maxplot=50)\n",
    "\n",
    "         "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "liberal-clock",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed)\n",
    "def clusterMinibatchKmean_single(inArray,ncluster,n_pcs,batchsize=100):\n",
    "    n_pcs=np.min([inArray.shape[0]-1,inArray.shape[1]-1,n_pcs])\n",
    "    batchsize=int(np.min([(inArray.shape[0]-1)/3,(inArray.shape[1]-1)/3,batchsize]))\n",
    "    inArray=pca.fit_transform(inArray)\n",
    "    labels = MiniBatchKMeans(n_clusters=ncluster,random_state=seed,batch_size=batchsize).fit_predict(inArray[:,:n_pcs])\n",
    "    return labels\n",
    "\n",
    "def clusterMinibatchKmean(inArray,nclusterL,n_pcs,sobj_coord_np):\n",
    "    for ncluster in nclusterL:\n",
    "        clusterRes=clusterMinibatchKmean_single(inArray,ncluster,n_pcs)\n",
    "#         print(clusterRes.shape)\n",
    "        savenamecluster='minibatchkmean_ncluster'+str(ncluster)+'n_pcs'+str(n_pcs)+'epoch'+str(plotepoch)\n",
    "        with open(os.path.join(clustersavedir,savenamecluster), 'wb') as output:\n",
    "            pickle.dump(clusterRes, output, pickle.HIGHEST_PROTOCOL)\n",
    "        plotembeddingbyCT(clusterRes,'minibatchkmean',[],embedding,savedir,plottype+' of '+s,savenameAdd=savenamecluster)\n",
    "        plotembeddingbyCT_contrast(clusterRes,'minibatchkmean',[],embedding,os.path.join(savedir,'contrast'),plottype+' of '+s,savenameAdd=savenamecluster)\n",
    "\n",
    "        plotembeddingbyCT(clusterRes,'minibatchkmean_location',[],sobj_coord_np,savedir,'location'+' of '+s,savenameAdd=savenamecluster)\n",
    "        plotembeddingbyCT_contrast(clusterRes,'minibatchkmean_location',[],sobj_coord_np,os.path.join(savedir,'contrast'),'location'+' of '+s,savenameAdd=savenamecluster)\n",
    "\n",
    "def clusterMinibatchKmean_allsample(embedding,savedir,clustersavedir,inArray,nclusterL,n_pcs,sobj_coord_np,samplenameList):\n",
    "    for ncluster in nclusterL:\n",
    "        clusterRes=clusterMinibatchKmean_single(inArray,ncluster,n_pcs)\n",
    "        savenamecluster='minibatchkmean_ncluster'+str(ncluster)+'n_pcs'+str(n_pcs)+'epoch'+str(plotepoch)\n",
    "        with open(os.path.join(clustersavedir,savenamecluster), 'wb') as output:\n",
    "            pickle.dump(clusterRes, output, pickle.HIGHEST_PROTOCOL)\n",
    "        plotembeddingbyCT(clusterRes,'minibatchkmean',[],embedding,savedir,plottype+' of all samples',savenameAdd=savenamecluster)\n",
    "        plotembeddingbyCT_contrast(clusterRes,'minibatchkmean',[],embedding,os.path.join(savedir,'contrast'),plottype+' of all samples',savenameAdd=savenamecluster,maxplot=50)\n",
    "\n",
    "        for s in plot_samples.keys():\n",
    "            sidx=(samplenameList==s)\n",
    "            plotembeddingbyCT(clusterRes[sidx],'minibatchkmean_location'+s,[],sobj_coord_np[sidx],savedir,'location'+' of '+s,savenameAdd=savenamecluster)\n",
    "            plotembeddingbyCT_contrast(clusterRes[sidx],'minibatchkmean_location'+s,[],sobj_coord_np[sidx],os.path.join(savedir,'contrast'),'location'+' of '+s,savenameAdd=savenamecluster,maxplot=50)\n",
    "\n",
    "         "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "exact-discovery",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "0\n",
      "0\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "#compute embeddings\n",
    "mulist={}\n",
    "for s in plot_samples.keys():\n",
    "    adj=adj_list[s]\n",
    "    adj_norm = preprocessing.preprocess_graph(adj)\n",
    "    adj_decode=None\n",
    "    if adj_decodeName == 'gala':\n",
    "        adj_decode=preprocessing.preprocess_graph_sharp(adj)\n",
    "    for xcorr in plot_sample_X:\n",
    "        samplename=s+'X_'+xcorr\n",
    "        features=featureslist[samplename]\n",
    "        if standardizeX:\n",
    "            features=torch.tensor(scale(features,axis=0, with_mean=True, with_std=True, copy=True))\n",
    "        if use_cuda:\n",
    "            model.cuda()\n",
    "            features = features.cuda().float()\n",
    "            adj_norm=adj_norm.cuda()\n",
    "            if adj_decodeName:\n",
    "                adj_decode=adj_decode.cuda()\n",
    "        \n",
    "        model.eval()\n",
    "        if adj_decodeName==None:\n",
    "            adj_recon,mu,logvar,z, features_recon = model(features, adj_norm)\n",
    "        else:\n",
    "            adj_recon,mu,logvar,z, features_recon = model(features, adj_norm,adj_decode)\n",
    "        if inverseAct=='leakyRelu':\n",
    "            muplot=inverseLeakyRelu(mu.cpu().detach().numpy())\n",
    "        else:\n",
    "            muplot=mu.cpu().detach().numpy()\n",
    "        if plotRecon:\n",
    "            if plotRecon=='meanRecon':\n",
    "                mulist[samplename]=features_recon[3].cpu().detach().numpy()\n",
    "        else:\n",
    "            mulist[samplename]=muplot\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "sitting-nelson",
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine all latents to one plot \n",
    "np.random.seed(seed)\n",
    "for xcorr in plot_sample_X:\n",
    "    latents=None\n",
    "    celltype_broad=None\n",
    "    celltype_sub=None\n",
    "    region=None\n",
    "    samplenameList=None\n",
    "    sobj_coord_np=None\n",
    "    \n",
    "    for s in plot_samples.keys():\n",
    "        if 'old' in s:\n",
    "            sampleidx=plot_samples_old[s]\n",
    "        else:\n",
    "            sampleidx=plot_samples[s]        \n",
    "        samplename=s+'X_'+xcorr\n",
    "        muplot=np.copy(mulist[samplename])\n",
    "            \n",
    "        if latents is None:\n",
    "            latents=muplot\n",
    "            if 'old' in s:\n",
    "                celltype_broad=olddata.obs.loc[olddata.obs['sample']==sampleidx,'top_level']\n",
    "                celltype_sub=olddata.obs.loc[olddata.obs['sample']==sampleidx,'cell_type_label']\n",
    "                region=olddata.obs.loc[olddata.obs['sample']==sampleidx,'region']\n",
    "                sobj_coord_np=olddata.obs.loc[olddata.obs['sample']==sampleidx,['x','y']].to_numpy()\n",
    "                samplenameList=np.repeat(s,muplot.shape[0])\n",
    "            else:\n",
    "                celltype_broad=scaleddata.obs.loc[scaleddata.obs['sample']==sampleidx,'top_level']\n",
    "                celltype_sub=scaleddata.obs.loc[scaleddata.obs['sample']==sampleidx,'cell_type']\n",
    "                region=scaleddata.obs.loc[scaleddata.obs['sample']==sampleidx,'region']\n",
    "                sobj_coord_np=scaleddata.obs.loc[scaleddata.obs['sample']==sampleidx,['x','y']].to_numpy()\n",
    "                samplenameList=np.repeat(s,muplot.shape[0])\n",
    "        else:\n",
    "            if 'old' in s:\n",
    "                latents=np.vstack((latents,muplot))\n",
    "                celltype_broad=np.concatenate((celltype_broad,olddata.obs.loc[olddata.obs['sample']==sampleidx,'top_level']),axis=None)\n",
    "                celltype_sub=np.concatenate((celltype_sub,olddata.obs.loc[olddata.obs['sample']==sampleidx,'cell_type_label']),axis=None)\n",
    "                region=np.concatenate((region,olddata.obs.loc[olddata.obs['sample']==sampleidx,'region']),axis=None)\n",
    "                sobj_coord_np=np.concatenate((sobj_coord_np,olddata.obs.loc[olddata.obs['sample']==sampleidx,['x','y']].to_numpy()),axis=0)\n",
    "                samplenameList=np.concatenate((samplenameList,np.repeat(s,muplot.shape[0])),axis=None)\n",
    "            else:\n",
    "                latents=np.vstack((latents,muplot))\n",
    "                celltype_broad=np.concatenate((celltype_broad,scaleddata.obs.loc[scaleddata.obs['sample']==sampleidx,'top_level']),axis=None)\n",
    "                celltype_sub=np.concatenate((celltype_sub,scaleddata.obs.loc[scaleddata.obs['sample']==sampleidx,'cell_type']),axis=None)\n",
    "                region=np.concatenate((region,scaleddata.obs.loc[scaleddata.obs['sample']==sampleidx,'region']),axis=None)\n",
    "                sobj_coord_np=np.concatenate((sobj_coord_np,scaleddata.obs.loc[scaleddata.obs['sample']==sampleidx,['x','y']].to_numpy()),axis=0)\n",
    "                samplenameList=np.concatenate((samplenameList,np.repeat(s,muplot.shape[0])),axis=None)\n",
    "       \n",
    "    origCT=np.unique(celltype_broad)\n",
    "#     celltypeplot=np.concatenate((origCT,list(combineCelltype.keys())),axis=None)\n",
    "    sampledir=os.path.join(plotsavepath,'combined'+xcorr)\n",
    "    if inverseAct:\n",
    "        sampledir+='_beforeAct'\n",
    "    savedir=os.path.join(sampledir,'embedding_'+plottype)\n",
    "    clustersavedir=os.path.join(sampledir,'cluster')\n",
    "    if not os.path.exists(sampledir):\n",
    "        os.mkdir(sampledir)\n",
    "    if not os.path.exists(savedir):\n",
    "        os.mkdir(savedir)\n",
    "    if not os.path.exists(clustersavedir):\n",
    "        os.mkdir(clustersavedir)\n",
    "    \n",
    "    if plottype=='umap':\n",
    "        reducer = umap.UMAP(n_neighbors=n_neighbors,min_dist=min_dist,random_state=seed)\n",
    "        embedding = reducer.fit_transform(latents)\n",
    "        savenameAdd='_nn'+str(n_neighbors)+'mdist0'+str(int(min_dist*100))+'epoch'+str(plotepoch)\n",
    "    elif plottype=='pca':\n",
    "        embedding=pca.fit_transform(latents)\n",
    "        savenameAdd='_epoch'+str(plotepoch)\n",
    "    if ifplot:\n",
    "        plotembeddingbyCT_str(samplenameList,'sample',[],embedding,savedir,plottype+'of all samples',savenameAdd=savenameAdd)\n",
    "        plotembeddingbyCT_str(celltype_broad,'celltype_broad',[],embedding,savedir,plottype+'all samples',savenameAdd=savenameAdd)\n",
    "        plotembeddingbyCT_str(celltype_sub,'celltype_sub',[],embedding,savedir,plottype+'all samples',savenameAdd=savenameAdd)\n",
    "        plotembeddingbyCT_str(region,'region',[],embedding,savedir,plottype+'all samples',savenameAdd=savenameAdd)\n",
    "\n",
    "        plotembeddingbyCT_contrast(celltype_sub,'celltype_sub',[],embedding,os.path.join(savedir,'contrast'),plottype+'all samples',savenameAdd=savenameAdd)    \n",
    "        plotembeddingbyCT_contrast(region,'region',[],embedding,os.path.join(savedir,'contrast'),plottype+'all samples',savenameAdd=savenameAdd)    \n",
    "        plotembeddingbyCT_contrast(samplenameList,'sample',[],embedding,os.path.join(savedir,'contrast'),plottype+'all samples',savenameAdd=savenameAdd)    \n",
    "\n",
    "\n",
    "    \n",
    "    if embedding.shape[0]<minCells:\n",
    "        continue\n",
    "    if ifcluster:\n",
    "        if 'leiden' in clustermethod:\n",
    "#             clusterLeiden_allsample(embedding,latents,n_neighbors,n_pcs,min_dist,resolution,sobj_coord_np,samplenameList,randseed=seed)\n",
    "            clusterLeiden_allsample(embedding,savedir,clustersavedir,latents,n_neighbors,n_pcs,min_dist,resolution,sobj_coord_np,samplenameList,randseed=seed)\n",
    "            assert np.sum(muplot-np.copy(mulist[s+'X_'+xcorr]))==0\n",
    "        if 'dbscan' in clustermethod:\n",
    "            clusterDBscan_allsample(embedding,latents,epslist,min_sampleslist,n_pcs,sobj_coord_np,samplenameList)\n",
    "            assert np.sum(muplot-np.copy(mulist[s+'X_'+xcorr]))==0\n",
    "        if 'agglomerative' in clustermethod:\n",
    "            clusterAgg_allsample(embedding,latents,nclusterlist,aggMetric,n_pcs,sobj_coord_np,samplenameList)\n",
    "            assert np.sum(muplot-np.copy(mulist[s+'X_'+xcorr]))==0\n",
    "        if 'kmeanbatch' in clustermethod:\n",
    "            clusterMinibatchKmean_allsample(embedding,latents,nclusterlist,n_pcs,sobj_coord_np,samplenameList)\n",
    "            assert np.sum(muplot-np.copy(mulist[s+'X_'+xcorr]))==0\n",
    "    #by region\n",
    "    for reg in np.unique(region):\n",
    "        savedir=os.path.join(sampledir,'embedding_'+plottype+'_'+reg)\n",
    "        clustersavedir=os.path.join(sampledir,'cluster'+'_'+reg)\n",
    "        if not os.path.exists(savedir):\n",
    "            os.mkdir(savedir)\n",
    "        if not os.path.exists(clustersavedir):\n",
    "            os.mkdir(clustersavedir)\n",
    "\n",
    "        reg_idx=region==reg\n",
    "\n",
    "        if plottype=='umap':\n",
    "            reducer = umap.UMAP(n_neighbors=n_neighbors,min_dist=min_dist,random_state=seed)\n",
    "            embedding = reducer.fit_transform(latents[reg_idx])\n",
    "            savenameAdd='_nn'+str(n_neighbors)+'mdist0'+str(int(min_dist*100))+'epoch'+str(plotepoch)\n",
    "        elif plottype=='pca':\n",
    "            embedding=pca.fit_transform(latents[reg_idx])\n",
    "            savenameAdd='_epoch'+str(plotepoch)\n",
    "        if ifplot:\n",
    "#             plotembeddingbyCT_str(samplenameList[reg_idx],'sample',[],embedding,savedir,plottype+' of '+'all samples'+' '+reg,savenameAdd=savenameAdd)\n",
    "#             plotembeddingbyCT_str(celltype_broad[reg_idx],'celltype_broad',[],embedding,savedir,plottype+' of '+'all samples'+' '+reg,savenameAdd=savenameAdd)\n",
    "#             plotembeddingbyCT_str(celltype_sub[reg_idx],'celltype_sub',[],embedding,savedir,plottype+' of '+'all samples'+' '+reg,savenameAdd=savenameAdd)\n",
    "#     #             plotembeddingbyCT(region,'region',[],embedding[reg_idx],savedir,'UMAP of '+s)\n",
    "\n",
    "#             plotembeddingbyCT_contrast(celltype_sub[reg_idx],'celltype_sub',[],embedding,os.path.join(savedir,'contrast'),plottype+' of '+'all samples'+' '+reg,savenameAdd=savenameAdd)\n",
    "            plotembeddingbyCT_contrast(samplenameList[reg_idx],'sample',[],embedding,os.path.join(savedir,'contrast'),plottype+' of '+'all samples'+' '+reg,savenameAdd=savenameAdd)\n",
    "\n",
    "        \n",
    "        if embedding.shape[0]<minCells:\n",
    "            continue\n",
    "        if ifcluster:\n",
    "            if 'leiden' in clustermethod:\n",
    "                clusterLeiden_allsample(embedding,latents[reg_idx],n_neighbors,n_pcs,min_dist,resolution,sobj_coord_np[reg_idx],samplenameList[reg_idx],randseed=seed)\n",
    "                assert np.sum(muplot-np.copy(mulist[s+'X_'+xcorr]))==0\n",
    "            if 'dbscan' in clustermethod:\n",
    "                clusterDBscan_allsample(embedding,latents[reg_idx],epslist,min_sampleslist,n_pcs,sobj_coord_np[reg_idx],samplenameList[reg_idx])\n",
    "                assert np.sum(muplot-np.copy(mulist[s+'X_'+xcorr]))==0\n",
    "            if 'agglomerative' in clustermethod:\n",
    "                clusterAgg_allsample(embedding,latents[reg_idx],nclusterlist,aggMetric,n_pcs,sobj_coord_np[reg_idx],samplenameList[reg_idx])\n",
    "                assert np.sum(muplot-np.copy(mulist[s+'X_'+xcorr]))==0\n",
    "            if 'kmeanbatch' in clustermethod:\n",
    "                clusterMinibatchKmean_allsample(embedding,latents[reg_idx],nclusterlist,n_pcs,sobj_coord_np[reg_idx],samplenameList[reg_idx])\n",
    "                assert np.sum(muplot-np.copy(mulist[s+'X_'+xcorr]))==0\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}