data_prep/spring_notebook_10X_CITEseq.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cellranger import\n",
    "#import modules, define some functions for loading, saving and processing a gene-barcode matrix\n",
    "%matplotlib inline\n",
    "import pkg_resources\n",
    "import collections\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.sparse as sp_sparse\n",
    "import tables\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import h5py\n",
    "\n",
    "np.random.seed(0)\n",
    "\n",
    "GeneBCMatrix = collections.namedtuple('GeneBCMatrix', ['gene_ids', 'gene_names', 'barcodes', 'matrix'])\n",
    "\n",
    "def get_matrix_from_h5(filename, genome):\n",
    "    with tables.open_file(filename, 'r') as f:\n",
    "        try:\n",
    "            dsets = {}\n",
    "            for node in f.walk_nodes('/' + genome, 'Array'):\n",
    "                dsets[node.name] = node.read()\n",
    "            matrix = sp_sparse.csc_matrix((dsets['data'], dsets['indices'], dsets['indptr']), shape=dsets['shape'])\n",
    "            if 'id' in dsets.keys():\n",
    "                return GeneBCMatrix(dsets['id'], dsets['name'], dsets['barcodes'], matrix)\n",
    "            else:\n",
    "                return GeneBCMatrix(dsets['genes'], dsets['gene_names'], dsets['barcodes'], matrix)\n",
    "        except tables.NoSuchNodeError:\n",
    "            raise Exception(\"Genome %s does not exist in this file.\" % genome)\n",
    "        except KeyError:\n",
    "            raise Exception(\"File is missing one or more required datasets.\")\n",
    "\n",
    "def save_matrix_to_h5(gbm, filename, genome):\n",
    "    flt = tables.Filters(complevel=1)\n",
    "    with tables.open_file(filename, 'w', filters=flt) as f:\n",
    "        try:\n",
    "            group = f.create_group(f.root, genome)\n",
    "            f.create_carray(group, 'genes', obj=gbm.gene_ids)\n",
    "            f.create_carray(group, 'gene_names', obj=gbm.gene_names)\n",
    "            f.create_carray(group, 'barcodes', obj=gbm.barcodes)\n",
    "            f.create_carray(group, 'data', obj=gbm.matrix.data)\n",
    "            f.create_carray(group, 'indices', obj=gbm.matrix.indices)\n",
    "            f.create_carray(group, 'indptr', obj=gbm.matrix.indptr)\n",
    "            f.create_carray(group, 'shape', obj=gbm.matrix.shape)\n",
    "        except:\n",
    "            raise Exception(\"Failed to write H5 file.\")\n",
    "        \n",
    "def subsample_matrix(gbm, barcode_indices):\n",
    "    return GeneBCMatrix(gbm.gene_ids, gbm.gene_names, gbm.barcodes[barcode_indices], gbm.matrix[:, barcode_indices])\n",
    "\n",
    "def get_expression(gbm, gene_name):\n",
    "    gene_indices = np.where(gbm.gene_names == gene_name)[0]\n",
    "    if len(gene_indices) == 0:\n",
    "        raise Exception(\"%s was not found in list of gene names.\" % gene_name)\n",
    "    return gbm.matrix[gene_indices[0], :].toarray().squeeze()\n",
    "\n",
    "def save_count_matrix_to_h5(E, gene_list, filename):\n",
    "    with h5py.File(filename, 'w') as hf:\n",
    "        hf.create_dataset(\"indptr\" ,  data= E.indptr)\n",
    "        hf.create_dataset(\"indices\",  data= E.indices)\n",
    "        hf.create_dataset(\"data\"   ,  data= E.data)\n",
    "        hf.create_dataset(\"shape\"  ,  data= E.shape)\n",
    "        hf.create_dataset(\"genes\"  ,  data= gene_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from helper_functions import *\n",
    "#from doublet_detector import *\n",
    "from collections import defaultdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.rcParams['font.family'] = 'sans-serif'\n",
    "plt.rcParams['font.family'] = 'sans-serif'\n",
    "plt.rcParams['font.sans-serif'] = 'DejaVu Sans'\n",
    "plt.rc('font', size=14)\n",
    "plt.rcParams['pdf.fonttype'] = 42"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2021-02-10 18:50:02--  http://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_10k_protein_v3/pbmc_10k_protein_v3_raw_feature_bc_matrix.h5\n",
      "Resolving cp-tinyproxy.default.svc.cluster.local (cp-tinyproxy.default.svc.cluster.local)... 10.105.184.90\n",
      "Connecting to cp-tinyproxy.default.svc.cluster.local (cp-tinyproxy.default.svc.cluster.local)|10.105.184.90|:3128... connected.\n",
      "Proxy request sent, awaiting response... 301 Moved Permanently\n",
      "Location: https://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_10k_protein_v3/pbmc_10k_protein_v3_raw_feature_bc_matrix.h5 [following]\n",
      "--2021-02-10 18:50:04--  https://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_10k_protein_v3/pbmc_10k_protein_v3_raw_feature_bc_matrix.h5\n",
      "Connecting to cp-tinyproxy.default.svc.cluster.local (cp-tinyproxy.default.svc.cluster.local)|10.105.184.90|:3128... connected.\n",
      "Proxy request sent, awaiting response... 200 OK\n",
      "Length: 154515158 (147M) [binary/octet-stream]\n",
      "Saving to: 'pbmc_10k_protein_v3_raw_feature_bc_matrix.h5.2'\n",
      "\n",
      "pbmc_10k_protein_v3 100%[===================>] 147.36M  59.7MB/s    in 2.5s    \n",
      "\n",
      "2021-02-10 18:50:08 (59.7 MB/s) - 'pbmc_10k_protein_v3_raw_feature_bc_matrix.h5.2' saved [154515158/154515158]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# download raw data from 10X Genomics\n",
    "!wget http://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_10k_protein_v3/pbmc_10k_protein_v3_raw_feature_bc_matrix.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pbmc_10k_protein_v3_raw_feature_bc_matrix.h5\n",
      "(6794880, 33555)\n",
      "6794880\n"
     ]
    }
   ],
   "source": [
    "# D stores all the data; one entry per library\n",
    "sample_list = [\"Sample1\"]\n",
    "D = {}\n",
    "\n",
    "for j, s in enumerate(sample_list):\n",
    "    filename = 'pbmc_10k_protein_v3_raw_feature_bc_matrix.h5'\n",
    "    raw_matrix_h5 = filename\n",
    "    print raw_matrix_h5\n",
    "    genome = \"matrix\"\n",
    "    D[s] = {}\n",
    "    D[s]['meta'] = {}\n",
    "    gbm = get_matrix_from_h5(raw_matrix_h5, genome)\n",
    "    D[s]['E'] = transpose(gbm.matrix)\n",
    "    D[s]['meta']['gene_list']=gbm.gene_names\n",
    "    D[s]['meta']['gene_id']=gbm.gene_ids\n",
    "    print D[s]['E'].shape\n",
    "    D[s]['cell_index']=gbm.barcodes\n",
    "    print (len(D[s]['cell_index']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_list = D[s]['meta']['gene_list']\n",
    "gene_id = D[s]['meta']['gene_id']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['CD3_TotalSeqB', 'CD4_TotalSeqB', 'CD8a_TotalSeqB', 'CD14_TotalSeqB', 'CD15_TotalSeqB', 'CD16_TotalSeqB', 'CD56_TotalSeqB', 'CD19_TotalSeqB', 'CD25_TotalSeqB', 'CD45RA_TotalSeqB', 'CD45RO_TotalSeqB', 'PD-1_TotalSeqB', 'TIGIT_TotalSeqB', 'CD127_TotalSeqB', 'IgG2a_control_TotalSeqB', 'IgG1_control_TotalSeqB', 'IgG2b_control_TotalSeqB']\n"
     ]
    }
   ],
   "source": [
    "# get protein names and gene names\n",
    "\n",
    "adt_ix = [i for i,g in enumerate(gene_list) if g.endswith('TotalSeqB')]\n",
    "gn_ix = [i for i,g in enumerate(gene_list) if not g.endswith('TotalSeqB')]\n",
    "\n",
    "print [gene_list[i] for i in adt_ix]\n",
    "proteins = gene_list[adt_ix]\n",
    "gns = gene_list[gn_ix]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter cells by total counts and number of genes detected\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEiCAYAAADNgWQ8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmcXFWZ//HP000DHbYOEoS0QBAhYAgQyQgYkO2HKKgTYdBBRJxxcEVliwaJEEaUYGSbEUdgdEAWBQF72GQzoMimwQ6ELSqyTSWBsHSApANN5/n9cW51bm7fqr7VXXt/369XvbrrbnVuL/XUOec555i7IyIiUmkttS6AiIiMDgo4IiJSFQo4IiJSFQo4IiJSFQo4IiJSFQo4IiJSFQo4Ik3IzGabmcY8SF1RwBEpwswmm9k1ZvaMma0ys8Vm9jszm13rslWKmW1oZmeY2S1mtszM3Mxm1rpc0vgUcEQKMLNpwHzg/cClwHHAT4AVwKzalaziNgNOAyYD3TUuizSRdWpdAJE69m1CcJnq7i/Fd5jZ+NoUqSqWAJ3uvtjMJgBP17Y40ixUwxEpbDvgsWSwAXD3xfnvzezjZnaTmeXM7E0ze9bMzjaz9eLnmNmlUbNcp5l1mdnrZrbEzE6M9m9vZrea2RvRtb6YOH+/qHnrqKjJa7GZrTSzO8xsYpYbMrODzOyu6LVXRN9/IHFvb8bvT6RcFHBECnsG2M3Mdh3iuH8F3gb+E/g6cDdwMvA/Kce2AL8BXgS+CSwCzjGzfwPuAJ6Iti8BfmJmU1Ku8S3gn4AfAucAewB3mdk7ihXSzI4EbgX6CU2CswjNZ/PMbI8h7lFk5NxdDz30SHkABxDenPuBBwlv8IcC6yeOG5Ny7ixgNfCu2LZLAQe+E9u2IfBadOy/xraPJwpisW37Ree/AHQkyunA92PbZod/74HnGwCvAJcmy05oMvttgZ/BhOjaM2v9+9Cj8R+q4YgU4O7zgH2AG4H3AicBNwEvmNm/xI5bCWBmLWa2iZltBvweMOB9KZf+79i5bwCPEYLa5bHti4HnCc16ST93955EOR8DDilyOwcBY4ErzWyz/IMQcO4E9jGztiLni4yYkgZEinD3+4DpZtYK7AJ8lBB4fmZmz7r7PDObBMwl1EDaE5foSDzvc/cliW3LgaXu3peyfWxKsf6asu0vwP5FbmWH6OvtRY7pAJYV2S8yIgo4Ihm4ez8hRbjbzO4Ffgt8xsweAu4CVhKy2p4CeoFOQhNashVhdYGX6C+w3dKKk/G4uHw5PgfkChyzfIhriIyIAo5I6f4YfR1PqFWMA/Zz99/lDzCzgyr4+jukbNseeLbIOU9FX5e5+53lL5LI0NSHI1KAmR1gZmn/I/m+kidZUzMZqGFE55xYwaJ91swGmurM7ABgEiH7rZBbgR5gVjJdO7rGuLKXUiRBNRyRwv4D2NDMfk1IV24hJAEcDbwMnE/IMHsZuMzM/hPoI6Qsb1jBcr0A3GtmPyX0uxwfbTun0Anu/rqZfQH4BfCwmV0JLAXeRailrQA+kj/ezI6Lrp0PbPubWf794j/dXc1vUjIFHJHCTgYOBw4GPg+sBywGrgS+5+7PAJjZoYQ3+zOAN4DrgP8CHqlQuc4GJkblGwvcC3zNUwaoxrn7r8wsR+hrOoGQobaEkPL934nDTwa2iT3/UPQAuAL198gwmLsmlBVpBGa2HyFB4Uh3/2WNiyNSMvXhiIhIVSjgiIhIVSjgiIhIVagPR0REqkJZajGbbbaZT5gwodbFEBFpKA899NBL7j7kWC4FnJgJEyYwf/78WhdDRKShmFmxWS4GqA9HRESqQgFHRESqQgFHRESqQgFHRESqQgFHRESqQllqIiINrqs7x9zbFrG4p5fxHe3MOHgi06d01rpYgyjgiIg0sK7uHKdcv5DevrA0U66nl1OuXwhQd0FHTWoiIg1s7m2LBoJNXm9fP3NvW1SjEhWmgCMi0sAW9/SWtL2WFHBERBrY+I72krbXkgKOiEgDm3HwRNrbWtfa1t7WyoyDJ9aoRIUpaUBEpIHlEwOUpSYiIhWXDDr5hIF6CzoKOCIiDa5RUqPVhyMi0uAaJTVaNRwRkQYQn02gY0wb7rC8t4/xHe3kGiQ1WgFHRKTOJZvMXl3ZN7Av19OLAZ5yXr2lRqtJTUSkzqU1mcU5YIlt9ZgarYAjIlLnCjWZxTnQ2dGORV/POmxyXSUMgJrURETqWld3LtNxrWbcO/OACpdmZFTDERGpY1kzzfo9rRenvijgiIjUsayZZp11liCQRgFHRKSOZck0q8cEgTQKOCIidSxtcs64ek0QSKOkARGROpYPJGfc+Nha429gTc2mEYINqIYjItIQXut9e9C23r5+TrrmYbadeTPT5szLnNFWK6rhiIjUsfwsA4Wy0PLb63XCzjjVcERE6thQswzE1eOEnXGq4YiI1In4BJ35hdRKnYCz3ibsjFPAERGpA4XWtOkY0zYoWaCYepuwM66qTWpmtqWZXWZmy8xslZk9bmb7xvabmc02s8Vm1mtmd5vZpMQ1xprZ5Wa2PHpcbmYdiWMmm9nvomvkzOw0M0vObSciUjcKrWnz6sq+1Ik5P7Pn1oPSpet9PE7VAk4UFO4lTGp6KLAT8DXgxdhh3wROirb/Q7TvDjPbKHbMVcD7gI8AH46+vzz2OhsDdwAvRNf4OjADOLES9yUiUg7FmsLis0Hnx92cOX0yZx02ue4n7IyrZpPaN4El7v7Z2Lan899ENZDjgTnufl207RhC0Pk0cJGZ7UQIMnu7+33RMV8E7jGzie6+CDgKGAMc4+69wKPReSea2bnuDTDhkIiMOsUWUoMQdFrNWNzTO5AYMH1KZ10HmKRqNqlNBx40s6vN7EUzW2Bmx8WaurYFtgBuz58QBYzfAx+INu0FvAHcF7vuvcCKxDH3ROfm3QaMByYkC2VmXzCz+WY2f9myZSO9RxGRYRlqRgEIKdDOmv6deh93k1TNgPNu4CvA34GDgQuAOcBXo/1bRF9fSJz3QmzfFsCyeC0l+v7FxDFp14i/xgB3v9jdp7r71HHjxpV6TyIiZTF9SieH7569tlLvKdBpqtmk1gLMd/dToufdZrY9IeD8KHZcsskruXpqWpPYUMdYge0iInXjpoeXlHR8PadAp6lmDWcJ8Hhi2xPA1tH3S6OvyVrI5qypoSwFNo9nnEXfj0sck3YNGFzzERGpC13dOXp6s6c/Q32nQKepZsC5F0jm6+0APBt9/zQhWByU32lm6wP7sKbP5n5gQ0I/Td5ewAaJY/aJzs07CFgMPDPSmxARKaeu7hzT5szj+KsXFD2uJZEbXe8p0GmqGXDOA/Y0s1PN7D1mdgQhZflCGOiLOR+YaWaHmdnOwKWEJIGromOeAG4lZKztaWZ7ARcBN0UZakTHrgQuNbOdzewwYCagDDURqSv5wZ7FstMAtt98AzZpbxt43tHeVvcp0Gmq1ofj7n8ys+nA94HvAM9FX38cO+wHQDshCI0FHgQ+5O6vx445CvgP1mSz3QAcF3ud5WZ2UHSN+cCrwDnAuRW4LRGRYcsyT5oBf31xxVrb3nx7dQVLVTmmD/1rTJ061efPn1/rYojIKLHtzJuHncnU2dHOvTMPKGt5hsvMHnL3qUMdp7nURERqZKjBnsXkenqZNmfeWhN91nsTmwKOiEgNdHXnWPHm4EXVSpEPVo2wFg5oPRwRkarLJwuUmgZdTCMMBFXAERGpslIWVStFvQ8EVcAREamyYoEhPvvzZ/bcms4SBnfW+0BQ9eGIiFRZoWSBQplnycXZ0jTCQFDVcEREqixtZuhiAWP6lM5Ba9/kaz+NshYOqIYjIlIxXd055t62aFDqcj4wpO0rpNHWvkmjgCMiUgHJZrBk6nKxAFIoUDU6NamJiFRAWiZaltTlWV0LOeHqBeR6eht6sbU0CjgiIhVQKBOtWIZaV3eOKx94btB0N40wxiYLBRwRkQoolKJcLHX5jBsfKzi3Wr2PsclCAUdEpAJKzUTr6s7x6srCMw/U+xibLJQ0ICJSAaVmohVrMjOo+zE2WRQMOGb2NGSbOdvd3122EomINIlSUpmLNZk5awJSI2erFavh/Cj2/YbAicAfCUs4Q1ja+f2Exc1ERKSALGnOQy1V0CgzQhdTMOC4+0AgMbNLgbPd/fvxY8zsFGBSxUonItLghhqPkzfj4IlDTl+Tz1Zr1ICTNWngMOCalO2/Aj5evuKIiDSXrONxktPXFNLI2WpZkwZWAPsBf0ts3w9YWcbyiIg0lVLG48T7fKbNmZfaxNbI2WpZA855wIVmNhV4INq2J3AMMLsC5RIRaVjxPpsWM/p9cP5VMnAk+3n233Ec1z2UW6t21AgzQhdjnvKDSD3Q7JPAN4Cdok1PABe4e1pTW0OaOnWqz58/v9bFEJEG1dWd44wbHys6ngZC4DjrsMlAaHLL9fRirJ0W3N7WyuG7d3LXk8vqfk41M3vI3acOdVzmcThRYGma4CIiUk5DrVnTasZq94HAAax1fNp0Nnc9uSx1fZxGlTngmNn6wEeB7YCL3L3HzLYDXnX3VypVQBGRRjDUstGr3Xl6zqEDz6fNmTfkMtONnCCQJlPAMbP3AHcSxuN0ELLTeoAvR8//rVIFFBFpBEMFh2SfTZZg0sgJAmmypkWfD9wOvBOI/5RuAPYvd6FERBpNseDQ1mqDOvuzBJP9dxw34nLVk6wB5wPAD909Wf97Dhhf3iKJiDSeGQdPLDh+Zp0WG9TZnza5Z9JdTy4rU+nqQymzRbelbNsaWF6msoiINKzpUzoLTj7Z27d60AJq06d0cvjuxTPOmq0PJ2vAuZ0wl1qem9nGwBnAzWUvlYhIA+nqzjFtzryixyRnFujqznHdQ8VX8RytfTgnAnub2SJgfeBq4BlgC2BmZYomIlL/urpzzPjVw0Un3oQwh9p2p9zCrK4wj9pQWW3NsiRBXKYsNXdfbGa7AUcC7yMEqouBK929uep8IiIlmH3DY/StzjaAvt+dKx54DijeXGbAUXtuXZeDPEcia1r0B4H73P1nwM9i29cxsw+6++8rVUARkXrW01t8VoE0v3jw+YLLEbSacc4nd226YAPZm9TuAjZN2b5JtE9ERDLqdy+4BHWzBhvIPtNAcpqfvHcQZpIWERmVxo5pG3LutKRWs5KXoG4GRQOOmd0QfevAFWb2Zmx3K7AzcF+FyiYiUvdO/9gkZlz7MH392fpxAI7cYyugtCWom8FQNZyXo68GvMraswy8BfwBuKQC5RIRqXv52aGHCjat0RIFrWYcucdWnDl9cpVKWF+KBhx3/xcAM3sGmOvuWmxNRIQoHTpDzabVjKfOOqRKpapvWZMGrgXek9xoZruY2XvLWyQRkfo397ZFmZrR0hZfG62yBpyLCf01Se+N9omIjCqlTDuTH+w52mUNOLsAf0zZ/idgdDZGisio1dWdo8UKTdU52C8efH6tc6fNmce2M29m2px5g+ZYa2ZZA04/YcxN0lgoOEGqiEjTya/sWUpTWf7Y/Lm5nl6cMN3N8VcvYMIoCT5Zx+H8DjjVzI7IL1FgZusApwKaZUBERo2h5kBL0xrVhoqdm+vpZca1DzP/2Ve468llTTk2J2vA+SYhBfpvZvaHaNvehBVAP1iJgomI1KPhLBlw5B5b0dWdG3KCz77+NXOtQQhCp1wf+n+aIehkalJz90WEfpyrCFPcvAO4EtjV3Z+oXPFEROpHqX03ANtvvgFTt9l0IHCUqrevf9DSBo0qaw0Hd19CaEITERl1htN3M227Tbny2L2YNmdeyc1wcc2yEFvmgGNmk4EvAu8GPu/uS8xsOvCsu3dXqoAiIvWglL6bsWPaOP1jkwaawUYaMJplIbasyxN8CLgB+A1wIJC/++2AzwHTK1E4EZF6MVT/S975n9ptUH9LoaUI0rQAq2PP29tam2Yhtqxp0d8FTnT3TxDmUMu7G3j/cF7YzL5tZm5mP4ptMzObbWaLzazXzO42s0mJ88aa2eVmtjx6XG5mHYljJpvZ76Jr5MzsNLMSG15FRGKyvIO0FjgobSmCQjYZ00ZnRzsGdHa0c9Zhk5siYQCyN6lNAm5J2f4K6evkFGVmewLHAo8kdn0TOIlQa1oEnAbcYWYT3f316JirgK2BjxBmsf5v4HLgY9G1NwbuIKRr/wMwEbiUsIzCOaWWVUQEIEvXTb97wayy9dtaMjXJ9azso/u0Dw2rjPUuaw3nVSAtxL4P+L9SXtDMNiFkuH0+um5+uwHHA3Pc/Tp3fxQ4BtgI+HR0zE7Ah4EvuPt97n4/oV/po2aWr3MeBYwBjnH3R939OuBs4ETVckSkVF3dOab8++2Zj09mleWTDbKumdMs/TVpsgacq4C5ZvYuQq1iHTPbF/gh8PMSX/Ni4Fp3n5fYvi2wBTDwm3X3XkJN5QPRpr2AN1h7DZ57CbWX+DH3ROfm3QaMByaUWFYRGcXyM0KXusBarqd3YPqak655OHOyQTP116TJGnBmAU8DzxIGez4OzCMMBv1e1hczs2MJs05/J2X3FtHXFxLbX4jt2wJY5r6mcht9/2LimLRrxF8jXqYvmNl8M5u/bNmyrLciIqNA1hmhkwwGpq/JmkZtRlP116TJ1Ifj7n3AUWb2HUIzWgvQ7e5/zfpCUZPX94F93P2tIocmfzvJ5a3TfntDHWMFtuPuFxPNeD116lTNIy4iA4abzjysNxJvjtkEislUw4myvMa4+9/d/Vp3v8bd/2pm7WZ2WsbX2gvYDHjUzN42s7eBfYGvRN/nVxdN1kI2Z00NZSmwebwvJvp+XOKYtGvA4JqPiEhB1exPaea+m7ysTWqnE5rSksZE+7LoIixlsFvsMR/4ZfT9XwjB4qD8CWa2PrAPa/ps7o/KsVfsunsBGySO2Sc6N+8gYDHwTMayioiU3J/S3tZKR3tbya/T1mpN3XeTlzUtOtlklTeFkBo9JHfvAXrWuqjZCuCVKCMNMzufMCv1k4QANIuQJHBVdI0nzOxW4KKoP8iAi4CbovneiI49HbjUzM4EdgBmAmfE+35ERIrp6s5xxo2PZT6+JeqDAZjxq4fpW13C8gUlHNvIitZwzOx1M3uNEGz+bmavxR4rCNlf15SxPD8AzgUuJNR+tgQ+FBuDAyHt+WFCNttt0fdH53e6+3JCjWZ8dI0LCeNvzi1jOUWkiZWaytzaYpz7ydgMAyUOwFjtlBTcGtVQNZzjCD+6nxEm7lwe2/cW8Ew0FmZY3H2/xHMHZkePQue8AnxmiOsuRMsmiMgwlbrmTf9qZ/6zrzB9SuewM9tKTb1uREUDjrtfBmBmTwP3RdlqIiJNbTjZaVc+8BxTt9l0RBN1dnXnmjpTLet6OL8DNjWzk83sv8xsMwAzm2Zm21a0hCIiVdYxpvSOfyc0i5W6Xk7cSdc83NTLTGdNi96dMLfZUYQpaTaOdh1ECQM/RUTqXVd3jp5hNm+9urKvpPVykvJzsTVr0MmaFv1D4AJ3nwK8Gdt+GzCt7KUSEamB/FQ2tcwZa6YVPpOyBpzdgctSti8B3lm+4oiI1M4ZNz42rA7/cmuWFT6Tso7D6QXGpmzfkTCPmYhIw+rqzjH3tkXDzhQbO6YNd+jpLU9eVbPOOpC1hvO/wOlmtl703M1sAmHa/+sqUC4RkarIj7nJuiJnmjHrrsPsj0/KvMhaXEsix6CZZ4zOGnBOJiy0towwnc0fgL8RZg6YVZmiiYhUXqljbtIs7ull+pROzjps8sBqnYVW/4zraG/j03ts3bQrfCZlnS36NWBvMzuANbNF/9nd76xk4UREKq0c/SX5JrDpUzoHgsW2M28e8rye3j6ueyjX1EEmLmsfDgDRomnJhdNERBrW+I72ETWnxZvA8n1Bi3t6aTHLlCKdz0ob1QGnhGUHcPd/L09xRESqa8I7RhZwDt891GryfUH55rlSxuM0a1ZaUrEazhGJ59sQ+m8WR8/HAysJU/4r4IhIw+nqznHfU5kmvC/opoeXcOb0ySPqC2rWrLSkggHH3SfnvzezfwE+Cxzj7s9F27YG/ge4stKFFBGphNk3PDbiQZ75VOjh1lKaOSstKWuW2mnA8flgAxB9fxLZF2ATEakbXd25so2bgcK1lGLZamPHtLHeOi2ccPUCps2Z17RT2uRlDTjvBNJ+musTlo0WEWko5Zo+Zmw00eeMgycOGofT3tbKkXtslbr9M3tuzaq+1fT09uFArqe3qedRg+wB5w7gEjPb08xao8eehNU276hc8UREKqMcHfVtrcbpH5sEMGgcTn5MzZnTJ6duv+vJZYP6fJp5HjXInhb9b4S51O4D8j+hFsLkncdWoFwiIhXVMaZtRIuedXa0s/+O45h72yJOuHoB4zvamXHwRO6deUDBcxxYunwVx1+9oOAxzZyxlnXg5zLgEDPbHtiJsAroE+7+l0oWTkSk3EY6b1re/juO47qHcgO1lHyTGLDWmJpS06WbOWOt1IGffwX+WqGyiIhUVPLNfyR+8eDzg4JHvEms1AGg0PwZayUFHBGRRlaOedPyCgWRfE2n1AGgxppBpM1KAUdEml6+GW0kMwoktRaoubSaDSuoOXDXk8vKULL6lTVLTUSkIZVj+YGkYunOI1liOtfTq7RoEZFGVc5mtLx8uvPhu3cODOxsNePw3TvpHMYA0LhmHouTKeCY2XvNbGLs+UFmdoWZnWJmpa84JCJSJZVIM5572yJmdS3kuodyAzWafneueyjH/juOS635nPPJXckScpp5LE7WGs5PgSkAZvYuwgqgmwJfBc6sTNFEREauEmnGuZ5ernjgudSBm3c9uSy15jN9Sicd0awEQ2nWsThZA85OwJ+j748AHnT3Q4CjgSMrUTARkXKYcfDETDWLcsn19KbWfGZ1LeSNVW9nukazjsXJGnBagbei7w8Ebom+f4owz5qISF2aPqWTD2y3adVeLy1Lrbevn188+Dx9q4dOKMiPxenqzjFtzjy2nXlz00zsmTXgPAp82cz2IQScW6PtncBLlSiYiEi5PPNydZqoimWpFcteS86zBgxk1jXTxJ5Zx+F8C+gCTgYuc/eF0faPA3+sRMFERMqhqzs34pTodVuN1paWotlu+b6aKx94LnWNnULjdjo72gfNvzZtzryCE3s28sDQTDUcd/89MA7YzN3/NbbrIuDLlSiYiMhIzepayAlFJsrMor2tlR/8066cddjkgaUI0o4555O7cteTy1KDjUHBcTtpU9kUShpo9GSCzONw3L0faDWzPcxsvWjbM+7+YsVKJyIyTF3duYK1jVKcddjkgVrFqr7Vg/Z3tLcNHFMoIDgUXKYgrcZSKGmg0ZMJMjWpmdlGwM+Awwk/u+2Bv5vZT4Cl7j67YiUUERmGciwfDWtmfi40gHSD9dYZOKbQkgf5waDTp2SbK23GwRMHTTLaDBN7Zq3hnA2MB94HxEP4TcAnyl0oEZGRKPfy0VC4OSvfP9TVnUtNe25rtZIDRaHF3Bq5/wayJw18HPiEuy8ws/iHhieAd5e/WCIiw1eukfobrLumz2V8R3tq8oGxZnLQtLTnDdZdZ1iBImttqJFkreGMBV5O2b4Ra1YAFRGpC+WaqLOttWUgFbnQAFJnzdo3aZaXuabVyLIGnD8Rajl5+TD+RcKy0yIiTaent29g/Mv0KZ0F+4QW9/QW7NB3aJqBmyOVNeB8G/iumV1CaIY70czmEaa2mVWpwomIlKKrO8cOp94y9IEliE+mWWgm6PEd7cw4eOKgtOe8Zhm4OVJZx+HcB3wAWJcwnc2BwGJgL3f/c7FzRUSqoas7x/FXL+Ct/nLkpq0t30SXFlTy2WPxjv40zTwLdFaZV/yMZhc4poJlEREZtjNufKxi184nBsRTpPPNaPlgA2s6+redeXNq81ujD9wcqZKWmDaz8cDmJGpGquWISK2ljX8pl3xiQD6gDJU9ViijrdEHbo5U1oGfU4ArgB1hUKKGE2aTFhFpWqXUTrIM3MynUqfVlJpV1hrOxcDzwLGEvpvyN5KKiAzTrK6FQx+U0dgCswWUUjsZqumtqzu3VkDKJxXEz21GWQPOe4Ep7v6XShZGRKRUXd05rnjguRFfp72thbMO2wWgLNPKFGt6S5smpxlmgx5K1oCzENgCUMARkbpSrmSBJ777kbWeV7K5q1lngx5K1oDzbeAHZjaLEHzWqm+6+yvlLpiIyFBmdS0sS7JAq63dNV3paWVGa1JB1oGfdwLvB24n9OEsix4vRV9FRKpqVtfCsjSlQVirppqKjedpZllrOPtXtBQiIiUqV7Bpa4Gp22xalmtlNVRSQbPKFHDc/XcjfSEzOwU4DJgIvAk8AJzi7o/GjjHgdOALhAlDHwS+6u6PxY4ZC/wHa+Z2uwH4mrv3xI6ZDPyIUCt7hbAy6XfdiywqLiINo5xTxPStpiYZYs04G/RQMq/4aWaTzexHZvYbM9sy2jY9GqOTxX7AjwlT5BwAvA3caWbxjxbfBE4Cvgb8A/AicEe0AFzeVYR1eT4CfDj6/vJYOTcG7gBeiK7xdWAGcGLWexWR+lbuWQU07Ux1ZB34+SFCTeI3hGCR79naDvgcMH2oa7j7wYlrHg0sB6YBN0a1m+OBOe5+XXTMMYSg82ngIjPbiRBk9o7md8PMvgjcY2YT3X0RcBQwBjjG3XuBR6PzTjSzc1XLEWlsXd25iswq0OwZYvUgaw3nu8CJ7v4J4K3Y9rsJzVbDsVH0+q9Gz7clpF7fnj8gChi/J9SKAPYC3mDtJRHuBVYkjrknOjfvNsKKpROShTCzL5jZfDObv2yZ8h9E6tlRl9zP8VcvqMi1mz1DrB5kDTiTgLQ5v18BhtvbdgGwALg/er5F9PWFxHEvxPZtASyL11Ki719MHJN2jfhrDHD3i919qrtPHTdu3HDuQ0SqYFbXQu59qjIjMIww2n+odWu6unNMmzOPbWferDVuhiFrltqrQCfwTGL7+4D/K/VFzexcYG9C01hyxdBkk5cltqU1iQ11jBXYLiINolxZaUnxN49iU8yM1uloyilrDecqYK6ZvYvwu1nHzPYFfgj8vJQXNLPzgCOBA9z977FdS6OvyVrI5qypoSwFNo/6e/LXM2Bc4pi0a8Dgmo8vPZDxAAAXjUlEQVSINIBK1CTO/9RudLS3DfoUWiiBoNh0NJJN1oAzC3gaeBbYEHgcmAf8Afhe1hczswsICQAHuPuTid1PE4LFQbHj1wf2YU2fzf3R6+8VO28vYIPEMftE5+YdRBiw+kzWsopI/Zh9Q3mz0j6z59ZAWEI6TVoCwWidjqacsq742efuRwE7AJ8kBI0d3f3olCaxVGZ2IfAvhNrNq2a2RfTYMHoNB84HZprZYWa2M3ApIUngquiYJ4BbCRlre5rZXoQxNjdFGWpEx64ELjWznc3sMGAmoAw1kQbU1Z0rGBiGo72thTOnTy5aM0lLICiUVKBkg+xKWoDN3Z8iLDE9HF+Jvv42sf0MYHb0/Q8IKdcXsmbg54fc/fXY8UcRBn7ms9luAI6LlXG5mR0UXWM+of/pHODcYZZbRGqgqzvHjF8toG91ea+bnxG6WM0kbYqZLGvcSHFZx+H8rMAuB1YBfwOudvfFha7h7smF29KOcULwmV3kmFeAzwxxnYXAB4d6PRGpT13dOU68egFljjW02poO/kITaI4d05aaBDBap6Mpp6w1nHGEvpTVQH4qmp0JCR4PEaas+Xcz28fdK5MkLyKjxtzbFpU92AD0xxrVC9VYTv/YpILnj8bpaMopa8C5l9CX8nl3XwlgZmOAS4CHgUMI2WrnAAdWoJwiMkp0dedSax7lphpL9WUNON8gZJatzG9w95Vm9j3gt+7+AzM7m7CMgYjIsHR15yo2k0Aa1ViqK2ta9IbAlinbt4j2AbxGiUkIIiJxlQ42LUP2JEslZQ04vwZ+amZHmNkEM9vGzI4AfgpcHx3zfrQEtYgM0y6n31rx11itgRE1lbVG8iVCWvEVsXPeBn4GnBw9fwI4tqylE5FRYcdTb2FVf+WjQafGzNRU1gXYVgJfMrOTCEsSGPA3d18RO0bZaSJSsqMuub8qwSY/ZqarO6dEgRopdeDnCuCRCpVFREaZru5cxWaAjuuMAgugCThrSJ38IlIT1chI+8yeW3Pm9MkDz6fNmVdwAk4FnMrLvMS0iEg5VTrYrLdOy1rBBjQBZ62phiMiVZPvP6nGwM6zD99lUH9Nx5i21OWpNQFndSjgiEhVhMk4H6avCrnJY9pC402yv6atxWhrNfpiSQqagLN61KQmIlVxyvWPVCXYtLYY3z9sl9QF0/pWOxusuw6dHe0YIZngrMMmq/+mSlTDEZGK6+rO0VvudQYKOOeIXZk+pZMTCvQRLe/tY8HpH6pKWWRtquGISMWVe8XOYuLLD6RRf03tKOCISEWVe8XOYuIzCcw4eCLtba1r7Vd/TW2pSU1EKuqMG6tTu2lrtbWCiZYfqD8KOCJSEfmU5LQ05EqY+0+7DgomWn6gvijgiEjZVTMFGkJTmgJL/VPAEZGyO/HqBRVZIjpNW4upX6ZBKGlARMrqqEvur1qwMWDuEYOb0qQ+KeCISNnM6lpYldmfIWScnfep3RRsGoia1ERkxI665P6qBRpYs9yAgk1jUcARkRGpZrBJLjcgjUUBR0RGpBrBpqO9jdkfn6QaTYNTwBGRYQmpz5VfWb6tBc191iQUcESkJF3dOU66ZgH9VRhi0wLMPWK3yr+QVIUCjogMqas7x+wbHqvanGgAY8e0cfrH1IzWTBRwRKSoru5cxZeDjlMGWvNSwBGRVNVcDnr7zTfgjhP3q/jrSG0p4IjIIF3duapNT6NU59FDAUdE1nLQuXfz1xdXVPx1Nli3le99Qss7jyYKOCJS1eYz1WhGLwUckVGuqzvHjGsfpq/Cec7KOhMFHJFRalbXQq588Dm8CuNpnplzaOVfROqeAo7IKDOrayFXPPBc1V6vs6O9aq8l9U0BR2QUqfaszu1trVocTQYo4Ig0mXwCwOKeXsZ3tLP/juO4+ZElvLqyerME5PX29TP3tkUA6rsRBRyRZtLVneOU6xfS29cPQK6nt6rNZ2lyPb2ccv1CQEFntFPAEWlQ1RovU0xHexsf3XXLtWpQBiTzEPI1nawBJ1lL01Q3zUEBR6QBzOpayC8efJ7+aqSUFdDe1sKqvtWpASA+rmbbmTennr844xiftFqaakjNQQFHpI50def49vWPsLKvGpPKZFfKYM3xHe2pA0jHZ8xWm3vbooFgk1dqDUnqkwKOSJV1def41nWP8Obb9RVU0hhwVIkzA8w4eOJaNRQoLVutUE0oaw1J6pd5Davo9Wbq1Kk+f/78YZ17xo2P8fji18pcImlkDz5dvfTjcmoxWO2wbmsLW23azmYbrlfyNV56402ef6WXt/pXF7xOoWO6n+vhrf7BwXjd1hambN0x7PuS4t47fmNO/9ikYZ1rZg+5+9ShjlMNR2QYGjWYFLPHtpuW7Vqbbbhe0UD10htv8vRLK1gdfd59q381T78UEiC22rR9rX0QguBWm2oAaaNTwCmT4X4ykNqaUKCDe7Tp7Gjn6i/uVbXXmzZn3loBBUKtalXfau6deYCy1JpU0wYcM/sKMAPYEngMON7d76ltqaQYvflXVlq6cl6up5ftTrmFfvdBK25W4s1/qH6a6VM6FWCaUFMGHDP7FHAB8BXgD9HX35jZe929tqPgYtLeYEud5DDrNYodpzf60WGo3tp8ynU8DRmoSIpyx5i21JkPOsa0DfuaUv9aal2ACjkRuNTdL3H3J9z9a8AS4Ms1LteAQm/ypbz5Z71GseMUbCRNPg25WIrySBTKVVIOU3NruhqOma0L7A78MLHrduAD1S+RSGMqloY80hTl5b3p87oV2i7NoRlrOJsBrcALie0vAFskDzazL5jZfDObv2zZsmqUT6QhjO9oLzhYM+sgzmLXrsR1pb41Y8DJS1bOU/tM3f1id5/q7lPHjRtXnZKJ1Ln8QM0ZB0+kva01dd9IVOq6Ut+arkkNeAnoZ3BtZnMG13pERo34J6784M6xY9pwh57ePlrNUrPUgLJnqeXPV+rz6NKUMw2Y2YPAw+7+hdi2vwDXufsphc4byUwDw6Estfo1bbtNOWLq1sy+4TF6on6FsWPaOHSXLbnryWVDvkkWSiUeaYqxxqdIPco600CzBpxPAZcT0qHvBb4EfB6Y5O7PFjqv2gFHRKQZjOqpbdz9ajN7BzCLMPDzUeCQYsFGREQqqykDDoC7/xj4ca3LISIiQTNnqYmISB1RwBERkapQwBERkapoyiy14TKzZUAtEgs2I4wfamS6h/rRDPehe6gPWe9hG3cfcuS8Ak4dMLP5WVIK65nuoX40w33oHupDue9BTWoiIlIVCjgiIlIVCjj14eJaF6AMdA/1oxnuQ/dQH8p6D+rDERGRqlANR0REqkIBR0REqkIBR0REqkIBp4bM7Ctm9rSZrTKzh8xsn1qXKSszO8XM/mRmr5nZMjO70cx2rnW5RsLMvm1mbmY/qnVZSmVmW5rZZdHvYpWZPW5m+9a6XFmZWauZfTf2//C0mZ1pZnU9wbCZfdDMbjCzXPS387nEfjOz2Wa22Mx6zexuM5tUo+KmKnYPZtZmZmeb2SNmtsLMlpjZVWa29XBeSwGnRqI1ey4Avg9MAe4DfjPcX2QN7EeYjfsDwAHA28CdZrZpLQs1XGa2J3As8Eity1IqM+sgrPtkwKHATsDXgBdrWa4SfQv4KvB1YEfgG9Hzggsm1okNCcuffAPoTdn/TeAkwu/jHwi/kzvMbKOqlXBoxe5hDPA+4HvR138EtgJuHc6HAWWp1Ui0Kukj7n5sbNtfgWuLrUpar8xsQ2A5MN3db6x1eUphZpsAfyYEnNOAR939uNqWKjsz+z6wr7tPq3VZhsvMbgJedvdjYtsuA97h7h+tXcmyM7M3gOPc/dLouQGLgR+5+/eibe2EoHOyu19Uq7IWkryHAse8F3gM2MXdF5ZyfdVwasDM1gV2B25P7LqdUGNoRBsR/p5erXVBhuFiQqCfV+uCDNN04EEzu9rMXjSzBWZ2XPSG1yj+AOxvZjvCwJvaAcAtNS3VyGwLbEHs/9zde4Hf07j/5wAbR19L/l+v6/bRJrYZ0Aq8kNj+AvD/ql+csrgAWADcX+uClMLMjgXeAxxd67KMwLsJy6mfB8wBdgP+M9rXKP1RZxM+tDxuZv2E96bvRQspNqotoq9p/+edVS5LWUQfls8BbnT3/yv1fAWc2kq2Z1rKtrpnZucCewN7u3t/rcuTlZlNJPSh7ePub9W6PCPQAsyPNcV2m9n2hD6QRgk4nwI+C3ya0FyzG3CBmT3t7j+taclGrln+z9cBrgA6gI8P5xpqUquNl4B+1nwCytucwZ+G6pqZnQccCRzg7n+vdXlKtBehtvmomb1tZm8D+wJfiZ6vV9viZbYEeDyx7QmgURJQAOYCP3T3X7r7Qne/HDiX+k8aKGZp9LUZ/s/XAX4B7AIc6O4vD+c6Cjg1EH2afgg4KLHrIEK2WkMwswsIn0gPcPcna12eYegCJhM+Tecf84FfRt83Sq3nXmBiYtsO1GZtp+EaQ/gQFtdPY79HPU0IOgP/52a2PrAPjfV/3gZcTQg2+7v70iFOKUhNarVzLnC5mf2R8IbxJWA88JOaliojM7uQ0O8xHXjVzPKf4t5w9zdqV7Ls3L0H6IlvM7MVwCvu/mhtSjUs5wH3mdmphDeGKYT04m/XtFSluRGYaWZPE5rUpgAnAj+vaamGEGVnvid62gJsbWa7Ef6GnjOz84FTzexJ4C/ALOAN4KqaFDhFsXsgZNn9ipDS/THAY//ry6MkiOzcXY8aPQgdvc8AbxJqPB+sdZlKKLsXeMyuddlGeF93E9JYa16WEst9KPAwsIrwxvZ1omEPjfAgJAycT6iV9QJ/J/SvrV/rsg1R7v0K/B9cGu03YDah2XMV8Dtg51qXO+s9ABOK/K9/rtTX0jgcERGpikZuHxURkQaigCMiIlWhgCMiIlWhgCMiIlWhgCMiIlWhgCMiIlWhgCMNz8zWjxaOaohp7IfLzHaM7rNhFrozszvNbMSDmc1sXDQT9lbR87r9WZjZVDN7NppVQGIUcEaZ6J+02OPSWpexVO6+CtgSuCPrOWb24eh+N6xcyaSMTgOuc/fnK/kiZvZJM7vdzF6K/j72TDlmjJn9l5m9bGZvmNn1ZrZlfr+7zycsaPa1Spa1ESngjD5bxh7Hpmz7RtpJ0XxKdSeaLh13X+rub9a6POVgZi1m1lrF11u3Wq81HGa2MfA5oBqzRm9IWJtnRpFjLiTM7HAEYbLXLYCuxPpD/wN81cz0HhtX62kV9KjdA/in8CcwaPuOhKkrjiBMxbEK+DfgnYS5unLASsKnuKMS5z5AmNtrLmEupqWEKUosdsynonN7gZeBuwgrO+b3/yNhEs1VhJm1u4B1on1LCXOE/ZywwugVwPpReT+aKP+nCOvzrCLMz7V/Yn/88ZNo34HAH4EVhHnWHgAmFvkZLgVOJUz4uYIw99TXE8dsSnizXAa8BswDdovt/1J0n/9ImPX5beA9RX4vqfcVHbMe4c3umejn+xfghMTP/5fAtcB3ovI+H21fH/gB8DxhuqW/AV+KnXcg8Kdo35Lo2LbY/o2i38eKaP/JwJ35n23sNc6J/oZWAA8SJn8t9nf6GWBxgZ/FziWUb2PCHGb58p2ULF/s2HdF198zsX2z6PdzeGzbe6Jj941tGwP0EZbsqPn/er08al4APWr4yx864DxFmJxzW8LEohOiN6/dCIt+fTX659s7du4DhEAwC9geOIow6+8nov3bROccF11vMvBFooATvd7bhPmndgJ2Jax33xbtX0oIBCcA20X/7IUCzrPAYdF1LiJMmrg5YfG7I6Nj3k34hLox4c36dUKAfHd03tHA9kV+hkuj+/0mYYbm46I3mkOi/a2EANYFTI1+Jj8gBONx0TFfit4k/0BYMmEisEGR30vqfUXHbACcHr3WttF9vk7sgwEh4LxOCEyTiN60gV9H154e3f+B+fOi39UqwsJuO0bHLCMskpa/7s+i8w+Mfq//Swiw8YBzHXAPYf2k7QgTdL4J7FTkZ3wRcEOBn8XOJZTvUsIcbQcAOwPXR7+7UgLOIdH2jRLbnwJOSWxbAJxa6//zenrUvAB61PCXP3TA+WqGa3QRm+ySEHDuShxzT/4YwtK6q4EtClzvIaKJDwvsXwr8KrGtUMA5KXZMK+FT/6zo+YejYzaMHTM+2rZHCT/DpYTVD+PbrgDujL4/hBBc1k0c8yRRTYgQcByYNMRrDXlfBc47H7gp9vyXhBpG/NP/5Oja+xW4xjmE2lS8pvQlQk23jVCLS37y7yAEw3zt8b2EDx/vTFz7VuDcIuW/FbiowM9i5xLLNz22fxNC4C0l4PwrsDLl+PuACxLbbgEuyfq3NBoeWp5AipkffxItwnQqIVB1AusSagW/SZz3SOL5YkLNAkKTxz3AIjO7ndDRf527vxy1ge9KeIPMXK4iBpa7dvd+M/sT4U0vlbsvNrNfAneb2Z3AbwnBLZf1dWLPT42+353wxvbK2k38rE/4hJ/X6+6PDfE6g14v7b7M7GuEPo+tgXbCG+6ixDUecfe+2PMphJrZPQVecyfgPo/eSSN/iK6/LTCWEPziZesxsydix+9O6Dd+KvGzWI9QyymknVB7KSZr+f4YK9/yaNmAUnnKtrRVPHuj15eIAo4UsyLx/FRCM9rxhE+TKwifLJMrY/Ylnjvhnx137zOz/QlNRx8CvgzMMbNphE/9wylX2bj7kWb2PuBg4HDg+2Z2qLvfNcxLtgD/R2hmSloe+37lMK+/FjM7Bjib0FT1IOET/Akpr5/8GRrFFVsW2TOcD+Fn0UcIbslrFfudvkQIGOUoX6FjsloKjDGzjdz99dj2tFU8NyU0q0lEGRRSir2BX7v7Ve7+MKE9fIdSL+Luq939Xnc/nfCp91XgiOjT6cOkvzkPx0BKa5T1NZWw9DKsWc1zUDaYu//Z3c9y930Ib9pHZ32d2PP86/yZ0FT3prv/LfFYVtrtDH69lPvaG7jH3X/i7t3u/jdCv9FQ/kyoCe1TYP/jwLREJtbehE/xzxCSE/oTZduEUPNIvsZmKT+LJUXK1k2RmmmJ5Xt/rHwbE5rmSvHH6DrxVTy3I/R5JVfxnES4Z4mohiOl+AtwqJntRei4P5HwZpp5KWMz24fwRnAH8CJhJcEtCW8YAN8DrolWfryGEBA+TGgfT9achvJ1M/s7oeb0DcKn0Iujfc9EXz8aNe2tJLTbfxa4idAMuD3hje7mIV5nXzM7mdBJfhDwz8Anon23EN50bjCzbxF+huMJfTs3uvuDJd7TUPf1F+AIMzuI8Hs5GtiDkJVVkLsvNLMbgMvM7BuEwL810OnuVxE6478KXGBmPyYkNnwXOC/6vbxiZpcD55hZD+F3+11Cf138Na4DrjSzkwif/jcjdOI/7u43FijebcBpZraxu79W4Jgs5bsiKt9yQkLBGVH5Bmo9ZvYOYCvWNAFvb2arCFlyL7r7S9F9nmdmrxKSIi4gNBX/PnadHQk1nDsLlHd0qnUnkh61ezB00sDOie2bATcQOoJfIASHnwK3xo55APhh4rxfAtdG308mvIG8yJrVKU9IHH844VPtm4Q3hl+zdlr0cYnjCyUN/DPhE+mbhIB2YOK8M6P7WE1Y2ruTkASxODrnmegeW4v8DPNp0dcSmoWWptzPJoSxG4sJNavngCuBbaL9XwJeyvD7GvK+op/FZYQPBK8SMrzOBJ5M+30krt9OSGnP3/9fgS/G9h9I6D97M7rPYmnHSwnZhcm06PWi8jwd/SyWRL/fXYe4927g88X+RjOW7xeEDxdLCB+Y7iUEpfwx+QSO5GNm7Jgx0d/LK9G9/hrYMlHeM4CuWv+P19tDK35K04k+XT4BTHb3Ryv8WkuBM939R5V8ndHOzKYTaiy7eJnetMysndC/dpq7X1iOa0bXHUMYw/RxD7MOSERNaiJS99y9y8zy48GGyhpMZWbvJ2SszSfUOk8l9CldW65yRiYA31GwGUwBR0QagrufN8JLGKGZbwdCc143sI+7J7PLRsTdH2dNn6TEqElNRESqQmnRIiJSFQo4IiJSFQo4IiJSFQo4IiJSFQo4IiJSFf8fxYeM2Pp/KVIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot total counts histograms - don't actually filter out any barcodes yet\n",
    "\n",
    "\n",
    "# adjust total counts thresholds\n",
    "\n",
    "for j,s in enumerate(sample_list):\n",
    "    D[s]['total_counts'] = np.sum(D[s]['E'][:,gn_ix], axis=1).A[:,0]\n",
    "    D[s]['genes_detected'] = np.sum(D[s]['E'][:,gn_ix]>0, axis=1).A[:,0]\n",
    "    D[s]['meta']['min_tot']= np.mean(D[s]['total_counts'])+np.std(D[s]['total_counts'])\n",
    "    D[s]['meta']['min_genes_detected']=200;\n",
    "\n",
    "    fig3 = plt.figure()\n",
    "    ax3 = fig3.add_subplot(111)\n",
    "    ax3.set_xlabel('Transcripts per barcode (log10)')\n",
    "    ax3.set_ylabel('genes detected')\n",
    "\n",
    "    ax3.scatter(log(1 + D[s]['total_counts']),D[s]['genes_detected']);\n",
    "    ax3.plot(ax3.get_xlim(),[D[s]['meta']['min_genes_detected'],D[s]['meta']['min_genes_detected']]);\n",
    "    title(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6794880\n",
      "Sample1 7965 / 6794880 3628.0 4664.67507847\n"
     ]
    }
   ],
   "source": [
    "D_orig=D\n",
    "for j,s in enumerate(sample_list):\n",
    "    ix = (D[s]['genes_detected'] > 200) & (D[s]['total_counts'] > 500)\n",
    "    print D_orig[s]['E'].shape[0]\n",
    "    if np.sum(ix) > 0:\n",
    "        print s, np.sum(ix), '/', D[s]['E'].shape[0], np.median(D[s]['total_counts'][ix]), np.mean(D[s]['total_counts'][ix])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---  Sample1 ---\n",
      "Pre-filter: 6794880 barcodes\n",
      "6794372\n",
      "(6794880, 33555)\n",
      "Post-filter: 7965 barcodes\n",
      "(7965, 33555)\n"
     ]
    }
   ],
   "source": [
    "# Actually filter out low-count barcodes\n",
    "for j,s in enumerate(sample_list):\n",
    "    print '---  %s ---' %s\n",
    "    print 'Pre-filter: %i barcodes' %D[s]['E'].shape[0]\n",
    "    tmpfilt = np.nonzero((D[s]['genes_detected'] > 200) & (D[s]['total_counts'] > 500))[0]\n",
    "    print max(tmpfilt)\n",
    "    print D[s]['E'].shape\n",
    "    D[s] = filter_dict(D[s], tmpfilt)\n",
    "    print 'Post-filter: %i barcodes' %D[s]['E'].shape[0]\n",
    "    print D[s]['E'].shape\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter cells by mito fraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['MT-ND1', 'MT-ND2', 'MT-CO1', 'MT-CO2', 'MT-ATP8', 'MT-ATP6', 'MT-CO3', 'MT-ND3', 'MT-ND4L', 'MT-ND4', 'MT-ND5', 'MT-ND6', 'MT-CYB']\n"
     ]
    }
   ],
   "source": [
    "# get mitochondrial genes\n",
    "mt_ix = [i for i,g in enumerate(gene_list) if g.startswith('MT-')]\n",
    "print [gene_list[i] for i in mt_ix]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7965 7902\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAGPCAYAAAAwSojkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X28HVV97/HPF1GRB4WW0IgWQ1ulCPh0QY0aRdsoFK9t8bYqqGCrWK1QikpptRhtaxGFi3q5V2ir0SjaKq1XUXkQ5eECYmNVhCKigg+EhCgYCaAC/d0/Zo7Zbk6ydk7OU875vF+vee2916xZe828Dtlf1qyZSVUhSZK0KdvMdAckSdLsZ2CQJElNBgZJktRkYJAkSU0GBkmS1GRgkCRJTQYGSVudJMuSeE24NI0MDNI8lmS/JP+S5MYkP0myKsnFSZbNdN+mSpIdk7w5yaeTrE1SSU6Y6X5Js52BQZqnkjwVWAk8EVgOvAZ4D3AH8MaZ69mU2xU4EdgP+PIM90Xaamw70x2QNGP+ii4c7F9VPxhckWT3menStLgZeFhVrUqyCLhhZrsjbR0cYZDmr18HrhkOCwBVtWrsfZLnJTknyU1JfprkO0neluSBg9skWd6f1nhYko8nuT3JzUmO69c/Msm5Sdb3bb1yaPsD+9MDh/enDFYluTPJBUn2GmWHkixN8vn+u+/o3z9laN9+Orh/kkZjYJDmrxuBxyV5bKPeHwH3AO8GjgEuAl4HvG+cutsAnwFuAY4HrgNOSfJy4ALg2r78ZuA9SR4/Tht/AfwP4B3AKcCTgM8n+eVNdTLJi4BzgXvpTqm8ke70w+eSPKmxj5IaPCUhzV8n0/2I/0eSlcClwOeBC6vqJwP1DquqOwc+n5HkeuAtSY6vqu8PrLs/8NGq+huAJCuAVcCZwMur6r19+ceB79KFkaOH+vUrwF5V9aO+7ueBC4HX0p1GuY8kOwCnAyuq6siB8jOAa4C3Ar810lGRNC5HGKR5qqo+BywBPgk8mu4H+RxgTZKXDdS7EyDJNkkekmRX4BIgwBPGafofB7ZdT/eDfS+wYqB8FfA9utMiwz4wFhYG+nkN8Dub2J2lwC7Ah5LsOrYA2wOfBZYkuf8mtpfU4AiDNI9V1eXA7yW5H/AY4Ll0weG9Sb5TVZ9Lsg/wduBA4EFDTew89Pnuqrp5qGwdsLqq7h6nfJdxunX9OGXfAJ65iV15VP96/ibq7Ays3cR6SZtgYJBEVd1Ld4nhl5NcRncK4MVJvkR3muJOutMB3wLuAh5Gdynm8Cjlf23kK+7dSHnG686I9QaN9eNI4KaN1FnXaEPSJhgYJA37Yv+6O93/1S8ADqyqi8cqJFk6hd//qHHKHgl8ZxPbfKt/XVtVn538LklyDoM0TyV5VpLx/g0YmyvwdTaMDPz8//D7bY6bwq69NMnPT3UkeRawD93VFxtzLvAj4I3Dl3v2bSyY9F5K84wjDNL89S5gxyT/Rne54zZ0kxhfAvwQOA34cf/+/UneDdxNd8njjlPYrzXAZUn+iW7ewbF92Skb26Cqbk9yFPBh4KtJPgSsBh5ON0pyB3DwWP0kr+nbHgsmz0wy9u/hu6vK0xfSEAODNH+9Dng+8Bzgj4EH0l0C+SHg76rqRoAkh9D9WL8ZWA+cDfwf4Kop6tfbgL36/u0CXAYcPd4NpgZV1UeT3EQ31+LP6a6QuBm4koErN3qvAx4x8PnZ/QLwQZzvIN1Hqnzgm6SZl+RAugmWL6qqj8xwdyQNcQ6DJElqMjBIkqQmA4MkSWpyDoMkSWpyhEGSJDV5WeWAXXfdtRYtWjTT3ZAkadp86Utf+kFVNW9uZmAYsGjRIlauXDnT3ZAkadok2dRt13/OUxKSJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAya3T5zQrdIkmaUT6vU7Lb6azPdA0kSjjBIkqQRGBgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTj7eepRad8KkJbXfjSYdMck8kSXKEQZIkjcDAIEmSmgwMkiSpycAgSZKaDAySJKlp2gJDkmVJamhZPbA+fZ1VSe5KclGSfYba2CXJiiTr+mVFkp2H6uyX5OK+jZuSnJgk07WfkiTNRdM9wnAd8NCBZb+BdccDrwWOBg4AbgEuSLLTQJ2zgCcABwMH9e9XjK1M8mDgAmBN38YxwOuB46ZmdyRJmh+m+z4M91TV6uHCfgTgWOCkqjq7LzuCLjQcBpyRZG+6kPC0qrq8r/NK4NIke1XVdcDhwPbAEVV1F3B1v91xSU6tqpqGfZQkac6Z7hGGX+tPE9yQ5CNJfq0v3xNYCJw/VrH/wb8EeEpftBhYD1w+0N5lwB1DdS7ttx1zHrA7sGiS90WSpHljOgPDlcCRdKcTXkEXEC5P8sv9e+hOJQxaM7BuIbB2cJSgf3/LUJ3x2mCgzi9IclSSlUlWrl27dnP3SZKkeWHaTklU1WcGPyf5AvBt4AjgC2PVhjbLUNl4pxRadbKR8rF+nQmcCbD//vt7ykKSpHHM2GWVVbUeuAZ4JDA2r2F4FGA3NowQrAZ2G7zioX+/YKjOeG3AfUceJEnSiGYsMCTZDvhN4GbgBrof+6VD65ewYc7CFcCOdPMUxiwGdhiqs6TfdsxSYBVw46TvhCRJ88R03ofhHUmekWTPJE8CPkb3Y//+fi7CacAJSQ5Nsi+wnG6S41kAVXUtcC7dFRNPTrIYOAM4p79Cgr7uncDyJPsmORQ4AfAKCUmStsB0Xlb5cODDwK7AWrp5C0+uqu/0608GHgScDuxCN0ny2VV1+0AbhwPvYsPVFJ8AXjO2sqrWJVnat7ESuA04BTh1ivZJkqR5YTonPb6wsb6AZf2ysTq3Ai9utPM14Omb30NJkrQxPktCkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUtOMBYYkf5WkkvyvgbIkWZZkVZK7klyUZJ+h7XZJsiLJun5ZkWTnoTr7Jbm4b+OmJCcmyXTtmyRJc82MBIYkTwZeAVw1tOp44LXA0cABwC3ABUl2GqhzFvAE4GDgoP79ioG2HwxcAKzp2zgGeD1w3FTsiyRJ88G0B4YkDwE+BPwxcNtAeYBjgZOq6uyquho4AtgJOKyvszddSDiqqi6vqiuAVwLPTbJX39ThwPbAEVV1dVWdDbwNOM5RBkmSJmYmRhjOBD5WVZ8bKt8TWAicP1ZQVXcBlwBP6YsWA+uBywe2uwy4Y6jOpf22Y84DdgcWTc4uSJI0v0xrYEjyCuA3gL8eZ/XC/nXNUPmagXULgbVVVWMr+/e3DNUZr43B7xjs01FJViZZuXbt2lF3RZKkeWXaAkN/yuCtwOFV9bNNVK2hzxkqG14/Sp1spJyqOrOq9q+q/RcsWLCJbkmSNH9N5wjDYmBX4Ook9yS5B3gG8Or+/Q/7esOjALuxYYRgNbDb4FyE/v2CoTrjtQH3HXmQJEkjmM7A8HFgP+BxA8tK4CP9+2/Q/dgvHdsgyXbAEjbMWbgC2JEufIxZDOwwVGdJv+2YpcAq4MbJ3CFJkuaLbafri6rqR8CPBsuS3AHc2l8RQZLTgDck+TpdgHgj3STHs/o2rk1yLnBGPx8iwBnAOVV1Xd/sWcCbgOVJ/hZ4FHAC8ObBuQ+SJGl00xYYRnQy8CDgdGAX4Erg2VV1+0Cdw4F3seFqik8ArxlbWVXrkizt21hJd+nmKcCpU957SZLmqBkNDFV14NDnApb1y8a2uRV4caPdrwFP3+IOSpIkwGdJSJKkERgYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNU0oMCR5UJLfTvKIye6QJEmafUYKDEmWJ3l1//4BwBeB84Hrkhw8hf2TJEmzwKgjDM8BvtC/fx6wE7AQWNYvkiRpDhs1MOwC3NK/Pwg4u6puAT4CPHoqOiZJkmaPUQPDamDfJPejG234bF++I3D3VHRMkiTNHtuOWO+9wD8Dq4B7gQv78icBX5+CfkmSpFlkpMBQVW9Jcg2wB/DRqvpZv+oe4G1T1TlJkjQ7jDrCQFWdPU7Z+ye3O5IkaTbaaGBIcuiojVTVv05OdyRJ0my0qRGGj43YRgH3m4S+SJKkWWqjgaGqvG20JEkCfJaEJEkagXMYJElSk3MYJElSk3MYJElSk6FAkiQ1jRwYkhyc5FNJrk3yq33Zy5P81ojb/2mSq5L8uF+uSHLIwPokWZZkVZK7klyUZJ+hNnZJsiLJun5ZkWTnoTr7Jbm4b+OmJCcmyaj7KUmS7mukwJDkcOBfgG8Ai4D796vuBxw/4nd9H/gL4AnA/sDngI8neUy//njgtcDRwAF0T8e8IMlOA22c1W9/MN1TM58ArBjo54OBC4A1fRvHAK8Hjhuxj5IkaRyjjjAcD7yiqv6c7vkRY74APG6UBqrq/1bVZ6rqm1X1jap6A3A7sLgfATgWOKmqzq6qq4EjgJ2AwwCS7E0XEo6qqsur6grglcBzk+zVf83hwPbAEVV1dX8767cBxznKIEnSxI0aGB4JXDFO+XrgwZv7pUnul+SFdI/HvhzYE1gInD9Wp6ruAi4BntIXLe6/7/KBpi4D7hiqc2m/7ZjzgN3pRkbG68tRSVYmWbl27drN3RVJkuaFUQPDKuBR45Q/HfjWqF/Wzy9YD/wUeA/w+1X1NbqwAN2phEFrBtYtBNZWVY2t7N/fMlRnvDYYqPMLqurMqtq/qvZfsGDBqLsiSdK8MurTKs8E3pXk5f3nX02yBDgZWLYZ33cd3SmMnYHnA+9PcuDA+hqqn6Gy4fWj1MlGyiVJ0ohGCgxVdXKSh9BNKNwO+DzdKME7qur0Ub+sqn4GfLP/uDLJAcCfA3/Xly0EvjewyW5sGCFYDeyWJGOjDP28hAVDdYZHEnbrX4dHHiRJ0ohGvqyyn6S4K/BE4MnAgqr660n4/gcCN9D92C8dW5FkO2AJG+YsXEE352HxwPaLgR2G6izptx2zlO6Uyo1b2FdJkuatkUYYkiwEtq2q7wMrB8ofDtxdVc3/e09yEvApuhGEsasfDgQOqapKchrwhiRfp7t88410kxzPAqiqa5OcC5yR5BV0pxrOAM6pquv6rzkLeBOwPMnf0s27OAF48+DcB0mStHlGHWFYQXfvg2HPYeA+CA0LgQ/SzWO4kO4+CQdX1Wf69ScDpwKn04WShwLPrqrbB9o4HPgq3dUU5/XvXzK2sqrW0Y0o7N63cTpwSt+uJEmaoFEnPR4AvGac8kuBt4/SQFUd2VhfdBMol22izq3AixvtfI3u6g1JkjRJRh1h2JZursGw7TZSLkmS5pBRA8OVwKvGKf9T4N8nrzuSJGk2GvWUxBuAzyV5LN38A4BnAY8HfnsqOiZJkmaPkUYYquoLdJcwfhs4lO6mSzcAi6vq8k1tK0mStn6jjjBQVV+lMeFQkiTNTSPfuEmSJM1fBgZJktRkYJAkSU0GBkmS1GRgkCRJTVsUGJL8bpKXTlZnJEnS7LSlIwwnAe+bjI5IkqTZa+T7MGzEbwH3n4yOSJKk2WuLAkNVrZqsjkiSpNlr5MCQ5IHA4cCjgQKuAT5cVT+dor5JkqRZYqQ5DEkeDXwDOBV4EvBk4DTgG0n2nrruSZKk2WDUSY/vBL4C7FFVS6pqCbAH8FW64CBJkuawUU9JPBU4oKp+PFZQVT9O8gbgC1PSM0mSNGuMOsLwE2Dnccof0q+TJElz2KiB4ZPAPyR5apL79cvTgDOAT0xd9yRJ0mwwamD4M+B64FK6EYWfABfTTYQ8dmq6JkmSZouR5jBU1Y+A303ySOA3gQD/WVXfnMrOSZKk2WGzbtxUVdfTjTRIkqR5ZHNu3PQCultB78bQqYyqet4k90uSJM0iIwWGJG+nm6vweWAV3Z0eJUnSPDHqCMNLgRdV1cemsjOSJGl2GvUqiW3o7vQoSZLmoVEDw5nAi6eyI5IkafYa9ZTEzsBhSZYCVwF3D66sqmMmu2OSJGn2GDUwPJoNpyR+c2idEyAlSZrjRr1x0zOnuiOSJGn2GnUOgyRJmscMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmqYtMCT5yyT/nuTHSdYm+WSSfYfqJMmyJKuS3JXkoiT7DNXZJcmKJOv6ZUWSnYfq7Jfk4r6Nm5KcmCTTsZ+SJM1F0znCcCDwv4GnAM8C7gE+m+SXBuocD7wWOBo4ALgFuCDJTgN1zgKeABwMHNS/XzG2MsmDgQuANX0bxwCvB46bip2SJGk+2Ha6vqiqnjP4OclLgHXAU4FP9iMAxwInVdXZfZ0j6ELDYcAZSfamCwlPq6rL+zqvBC5NsldVXQccDmwPHFFVdwFX99sdl+TUqqrp2F9JkuaSmZzDsFP//bf1n/cEFgLnj1Xof/AvoRuVAFgMrAcuH2jnMuCOoTqX9tuOOQ/YHVg0qXsgSdI8MW0jDON4J/AV4Ir+88L+dc1QvTXAwwbqrB0cJaiqSnLLwPYLge+P08bYuhsGVyQ5CjgKYI899pjQjswmi0741IS2u/GkQya5J5KkuWRGRhiSnAo8DXh+Vd07tHr4lEGGysY7pdCqk42UU1VnVtX+VbX/ggULmn2XJGk+mvbAkOR/Ai8CnlVV3x5Ytbp/XTi0yW5sGCFYDew2eMVD/37BUJ3x2oD7jl5IkqQRTGtgSPJOugmMz6qqrw+tvoHux37pQP3tgCVsmLNwBbAj3TyFMYuBHYbqLOm3HbMUWAXcOCk7IknSPDOd92E4HXgZ3ejCbUkW9suO0M1FAE4DTkhyaH+PhuV0kxzP6utcC5xLd8XEk5MsBs4AzumvkKCveyewPMm+SQ4FTgC8QkKSpAmazkmPr+5fLxwqfzOwrH9/MvAg4HRgF+BK4NlVdftA/cOBd7HhaopPAK8ZW1lV65Is7dtYSXcVxinAqZO1I5IkzTfTeR+G5p0W+xGAZWwIEOPVuRV4caOdrwFP37weSpKkjfFZEpIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKapjUwJHl6kk8kuSlJJTlyaH2SLEuyKsldSS5Kss9QnV2SrEiyrl9WJNl5qM5+SS7u27gpyYlJMg27KEnSnDTdIww7AlcDfwbcNc7644HXAkcDBwC3ABck2WmgzlnAE4CDgYP69yvGViZ5MHABsKZv4xjg9cBxk7wvkiTNG9tO55dV1aeBTwMkWT64rh8BOBY4qarO7suOoAsNhwFnJNmbLiQ8raou7+u8Erg0yV5VdR1wOLA9cERV3QVc3W93XJJTq6qmYVclSZpTZtMchj2BhcD5YwX9D/4lwFP6osXAeuDyge0uA+4YqnNpv+2Y84DdgUVT0XFJkua62RQYFvava4bK1wysWwisHRwl6N/fMlRnvDYGv+PnkhyVZGWSlWvXrt2C7kuSNHfNpsAwZviUQYbKxjul0KqTjZRTVWdW1f5Vtf+CBQs2t6+SJM0LsykwrO5fh0cBdmPDCMFqYLfBKx769wuG6ozXBtx35EGSJI1gNgWGG+h+7JeOFSTZDljChjkLV9BdabF4YLvFwA5DdZb0245ZCqwCbpyKjkuSNNdN930YdkzyuCSP6797j/7zHv1chNOAE5IcmmRfYDndJMezAKrqWuBcuismnpxkMXAGcE5/hQR93TuB5Un2TXIocALgFRKSJE3QdI8w7A98uV8eBLy5f/+Wfv3JwKnA6cBK4KHAs6vq9oE2Dge+Snc1xXn9+5eMrayqdXQjCrv3bZwOnNK3K0mSJmC678NwERsmII63voBl/bKxOrcCL258z9eAp0+kj5Ik6b5m0xwGSZI0SxkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVLTtjPdgblu0QmfmukuSJK0xRxhkCRJTQYGSZLUZGCQJElNBgZJktRkYJAkSU0GBkmS1GRgkCRJTd6HQcDE7xdx40mHTHJPJEmzkSMMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKlpzt7pMcmrgdcDDwWuAY6tqktntldzj3eIlKT5YU6OMCR5AfBO4K3A44HLgc8k2WNGOyZJ0lZqro4wHAcsr6p/6D8fneQg4FXAX85ctzRm1JGJjzzghwC8cIIjGWMc0ZCkLTPnAkOSBwD/DXjH0KrzgadMf480G3jqRJK2zJwLDMCuwP2ANUPla4DfHq6c5CjgqP7j+iTXTUF/fjDJbc4bi7uXXeG5M3IM87aZ+NYp4d/hlvMYbjmP4ZabimP4iFEqzcXAMKaGPmecMqrqTODMqepEkpVVtf9UtT8feAy3nMdwy3kMt5zHcMvN5DGci5MefwDcCywcKt+N+446SJKkEcy5wFBVPwO+BCwdWrWU7moJSZK0mebqKYlTgRVJvghcBvwJsDvwnhnoy5Sd7phHPIZbzmO45TyGW85juOVm7Bim6j6n9eeE/sZNx9PduOlq4M+r6pKZ7ZUkSVunORsYJEnS5JlzcxgkSdLkMzBsgSSvTnJDkp8k+VKSJY36z+jr/STJt5P8yXT1dTbbnOOY5NAk5ydZm+T2JFcmed509nc22ty/xYHtnpbkniRXT3UfZ7sJ/Pf8gCRv6bf5aZLvJjlmuvo7G03gGB6W5CtJ7kyyOskHkwxf4TZvJHl6kk8kuSlJJTlyhG32S3Jxkrv67U5Mkqnon4Fhgjb3eRVJ9gQ+3dd7PPD3wLuTPH96ejw7TeC5H88APgcc0tf/NPBvo/5AzkUTfXZKkl2ADwAXTnknZ7kJHsMPAwfR3fhtL+APgKumuKuz1gT+TXwqsAJ4P7AP8HvAo4EPTUuHZ6cd6ebc/RlwV6tykgcDF9DdMuAA4Bi6hy4eNyW9qyqXCSzAlcA/DJVdD/z9Ruq/Dbh+qOwfgStmel+2puO4kTa+CJwy0/uytR1D4F+BNwHLgKtnej+2pmMIPBtYB+w6032fLcsEjuHrgO8Mlb0MWD/T+zIbFmA9cGSjzquAHwMPGih7I3AT/RzFyVwcYZiAgedVnD+0alPPq1g8Tv3zgP2T3H9ye7h1mOBxHM9OwG2T1a+tyUSPYX8V0ULgb6eud1uHCR7D3wP+HTguyfeTXJ/kXUl2nMKuzloTPIaXAQ9N8t/T2RV4Id2ooUazGLi0qgZHI86ju43Aosn+MgPDxGzqeRUbO/+2cCP1t+3bm48mchx/QZI/BR5ON7Q5H232MUyyH93IwuFVde/Udm+rMJG/w18DngY8Fng+8Bq60xPLp6aLs95mH8OqugJ4Ed0piJ8Ba+lu4X/E1HVzztnY78rYukllYNgyIz2volF/vPL5ZnOPY1epm//xdrofvu9MRce2IiMdwyQPBD4CvK6qbpiOjm1FNufvcJt+3WFVdWVVnUcXGp6f5FemsI+z3cjHMMmjgXcBf0M3OnEQ3Y/cGVPZwTlo2n5X5uqdHqfaRJ5XsXoj9e8Bfjipvdt6TPi5H31YWAG8tKo+MTXd2yps7jF8KN3EsvcleV9ftg2QJPcAv1NVw8PKc91E/g5vBm6qqnUDZdf2r3tsYru5aiLH8C+BL1bV2/vPVyW5A7g0yRuq6ntT09U5ZWO/KzAFf4OOMExATex5FVdw38drLwVWVtXdk9vDrcMEjyNJ/hD4IN2EoI9NXQ9nvwkcw5uA/YDHDSzvAb7Zv593z1uZ4N/hZcDuQ3MWHtW/zrvRrgkew+3pQsagsc9TclngHHQFsCTJdgNlS4FVwI2T/m0zPRN0a12AF9Cdd3s5sDfd5UTrgUf06z8AfGCg/p7AHcBpff2X99s/f6b3ZSs7ji8E7qa77GjhwPJLM70vW8sxHGf7ZXiVxOb+He4IfA/4KN0lgU+luxzuozO9L1vRMTyy/2/5VXRzQp5KN5H0SzO9LzN4DHdkQ5C/Ezixf79Hv/7vgQsH6j+EbpThI8C+wKF0V028dkr6N9MHaGtegFfTpbif0qXrpw+suwi4aKj+M4D/6OvfAPzJTO/DbFg25zj2n2uc5aLp7vdsWjb3b3Fo23kfGCZyDOkUou9dAAADFUlEQVTuvXB+/w/7TcDpwE4zvR9b2TE8GrimP4Y3A2cBD5/p/ZjB43fgRv59W96vXw7cOLTNfsAlwE/6Y/gmpuCSyqryWRKSJKnNOQySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySZlSSo5J8N8l/JVk20/2RND4Dg6RxJVmepJL84zjrTu7XndN/vrH/vLHloo18xy50Nzx6O/Aw4B1TuEuStoCBQdKmfA94QZIdxgqSbAu8BPjuQL0D6B5s9VC6pw4CPHGg7NCNtP8IuofgnVNVN1fV+uEKSbZN4rMFpBlmYJC0KVcB1wN/OFB2CN1taC8aK6iqtVW1uqpWs+Hpqz8vq6pbhxtOciTw5f7jt/uRiEVJliW5OsmRSb5Fd5vhHZIclOTSJLcluTXJeUn2Hmpz9yQfSvLDJHcm+UqSZ07SsZDmNQODpJZ/Av5o4PMfAe+ju8f9lvhn7jsaMfZI4z2Bw4A/AB5LF1B2oHt42xPp7rm/DvhkkgcA9KMgFwOLgN+nu8f+W7awj5J62850ByTNemcB70jySOB2uh/5o9nCH+OquivJL4xGAPRnHx4AvKSq1gxscvbg9kleRvdkvicC/48uYCwEFlfVD/pq39qSPkrawMAgaZOq6rYk/0Y3svAjuicOfneKpxV8fygskOTXgb8BngQsoBsh3QbYo6/yeOCqgbAgaRIZGCSN4r3A+4H1wInT8H13jFP2SbrHSL+yf70H+E+60QgAJ0ZKU8g5DJJGcSHwM2BX4OPT/eVJfhnYG3hrVX22qq4FduIX/6fnP4DHJNl1uvsnzQcGBklNVVXAY4A9q+qnM9CF24AfAK9I8htJngG8h26UYcxZwC3Ax5MsSbJnkueNXSWR5IlJvp7kidPee2kOMDBIGklV3V5VP56h7/4v4AV0oeVqups9/TXdJZdjde4AnkF3uuKTwDXAm9lwNcf2wF79q6TNlO5/HCRJkjbOEQZJktRkYJAkSU0GBkmS1GRgkCRJTQYGSZLUZGCQJElNBgZJktRkYJAkSU0GBkmS1PT/ATpGvhmy9i8sAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot mito-gene frac histograms - don't actually filter out any cells yet\n",
    "# set mito-gene frac threshold\n",
    "for j,s in enumerate(sample_list):\n",
    "    D[s]['meta']['max_mt'] = 0.2\n",
    "\n",
    "for j,s in enumerate(sample_list):\n",
    "    fig = plt.figure(figsize=(8, 6))\n",
    "    ax = fig.add_subplot(111, xscale='linear', yscale='linear',\n",
    "        xlabel='MT frac.', ylabel='no. cells')\n",
    "\n",
    "    D[s]['mito_frac'] = np.sum(D[s]['E'][:,mt_ix], axis=1).A[:,0] / np.sum(D[s]['E'][:,gn_ix], axis=1,dtype=float).A[:,0]\n",
    "\n",
    "    ax.hist(D[s]['mito_frac'], cumulative=False, \n",
    "            bins=np.linspace(0, 1, 30))\n",
    "\n",
    "    ax.plot([D[s]['meta']['max_mt'],D[s]['meta']['max_mt']],ax.get_ylim());\n",
    "    title(s)\n",
    "\n",
    "    print D[s]['E'].shape[0], np.sum(D[s]['mito_frac'] <= D[s]['meta']['max_mt'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---  Sample1 ---\n",
      "Pre-filter: 7965 barcodes\n",
      "Post-filter: 7902 barcodes\n"
     ]
    }
   ],
   "source": [
    "# Actually filter out mito-high cells \n",
    "D_filt=D\n",
    "for s in sample_list:\n",
    "    print '---  %s ---' %s\n",
    "    print 'Pre-filter: %i barcodes' %D[s]['E'].shape[0]\n",
    "    tmpfilt = np.nonzero(D[s]['mito_frac'] <= D[s]['meta']['max_mt'])[0]\n",
    "    D[s] = filter_dict(D[s], tmpfilt)\n",
    "    print 'Post-filter: %i barcodes' %D[s]['E'].shape[0]\n",
    "    \n",
    "#del tmpfilt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merge data, normalize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample1\n"
     ]
    }
   ],
   "source": [
    "# create master dataset (all SPRING subsets will refer back to this)\n",
    "\n",
    "samp_lookup = {}\n",
    "samp_id_flat = np.array([],dtype=str)\n",
    "\n",
    "for s in D.keys():\n",
    "    samp_id_flat = np.append(samp_id_flat, [s] * D[s]['E'].shape[0])\n",
    "\n",
    "        \n",
    "E = scipy.sparse.lil_matrix((len(samp_id_flat), len(gene_list)), dtype=int)\n",
    "total_counts = np.zeros(len(samp_id_flat), dtype=int)\n",
    "mito_frac = np.zeros(len(samp_id_flat), dtype=float)\n",
    "genes_detected = np.zeros(len(samp_id_flat), dtype=int)\n",
    "\n",
    "for s in D.keys():\n",
    "    print s\n",
    "    E[samp_id_flat == s, :] = D[s]['E']\n",
    "    total_counts[samp_id_flat == s] = D[s]['total_counts']\n",
    "    mito_frac[samp_id_flat == s] = D[s]['mito_frac']\n",
    "    genes_detected[samp_id_flat == s] = D[s]['genes_detected']\n",
    "\n",
    "E = E.tocsc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "E_full=E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['MT-ND1', 'MT-ND2', 'MT-CO1', 'MT-CO2', 'MT-ATP8', 'MT-ATP6', 'MT-CO3', 'MT-ND3', 'MT-ND4L', 'MT-ND4', 'MT-ND5', 'MT-ND6', 'MT-CYB']\n",
      "20922 / 33555\n"
     ]
    }
   ],
   "source": [
    "# remove genes that are not expressed by any cells\n",
    "# optionally remove mito and rps genes\n",
    "gene_list = D[sample_list[0]]['meta']['gene_list']\n",
    "gene_id = D[sample_list[0]]['meta']['gene_id']\n",
    "remove_crapgenes = 1\n",
    "\n",
    "if remove_crapgenes:\n",
    "    import re\n",
    "    mt_ix = [i for i,g in enumerate(gene_list) if g.startswith('MT-')]\n",
    "    print [gene_list[i] for i in mt_ix]\n",
    "    #rp_ix = [i for i,g in enumerate(gene_list) if g.startswith('RPS') or g.startswith('RPL')]\n",
    "    #rp_ix = [i for i,g in enumerate(gene_list) if g.startswith('RPS','RPL')]\n",
    "    rp_ix = [i for i,g in enumerate(gene_list) if g.startswith('RP1') \n",
    "             or g.startswith('RP2') or g.startswith('RP3')\n",
    "             or g.startswith('RP4') or g.startswith('RP5')\n",
    "             or g.startswith('RP6') or g.startswith('RP7')\n",
    "             or g.startswith('RP8') or g.startswith('RP9')\n",
    "             or g.startswith('RPL') or g.startswith('RPS')\n",
    "            ]\n",
    "    #or g.startswith('RP4','RP5','RP6') or g.startswith('RP7','RP8','RP9')]\n",
    "    keep_genes = (E.sum(0) > 0).A.squeeze() #* rp_ix * mt_ix\n",
    "    keep_genes[rp_ix] = 0\n",
    "    keep_genes[mt_ix] = 0\n",
    "    print sum(keep_genes), '/', len(keep_genes)\n",
    "else:\n",
    "    keep_genes = (E.sum(0) > 0).A.squeeze()\n",
    "    print sum(keep_genes), '/', len(keep_genes)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(7902, 20922)\n"
     ]
    }
   ],
   "source": [
    "# filter out genes\n",
    "E = E[:,keep_genes]\n",
    "print shape(E)\n",
    "gene_list = gene_list[keep_genes]\n",
    "gene_id = gene_id[keep_genes]\n",
    "merged_list=[x+\"_\"+gene_id[i] for i,x in enumerate(gene_list)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['CD3_TotalSeqB_CD3', 'CD4_TotalSeqB_CD4', 'CD8a_TotalSeqB_CD8a', 'CD14_TotalSeqB_CD14', 'CD15_TotalSeqB_CD15', 'CD16_TotalSeqB_CD16', 'CD56_TotalSeqB_CD56', 'CD19_TotalSeqB_CD19', 'CD25_TotalSeqB_CD25', 'CD45RA_TotalSeqB_CD45RA', 'CD45RO_TotalSeqB_CD45RO', 'PD-1_TotalSeqB_PD-1', 'TIGIT_TotalSeqB_TIGIT', 'CD127_TotalSeqB_CD127', 'IgG2a_control_TotalSeqB_IgG2a', 'IgG1_control_TotalSeqB_IgG1', 'IgG2b_control_TotalSeqB_IgG2b']\n"
     ]
    }
   ],
   "source": [
    "# get protein names and gene names\n",
    "adt_ix = [i for i,g in enumerate(gene_list) if g.endswith('TotalSeqB')]\n",
    "gn_ix = [i for i,g in enumerate(gene_list) if not g.endswith('TotalSeqB')]\n",
    "\n",
    "print [merged_list[i] for i in adt_ix]\n",
    "proteins = gene_list[adt_ix]\n",
    "gns = gene_list[gn_ix]\n",
    "merged_proteins = [merged_list[i] for i in adt_ix]\n",
    "\n",
    "# normalize by total counts seperately\n",
    "E[:,gn_ix] = tot_counts_norm_sparse(E[:,gn_ix])[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save base directory files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "main_spring_dir=os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if BCR data\n",
    "BCR = 1\n",
    "\n",
    "if BCR:\n",
    "    barcode_dir = main_spring_dir + \"/barcodes/\"\n",
    "    if not os.path.exists(barcode_dir):\n",
    "        os.makedirs(barcode_dir)\n",
    "        for s in sample_list:\n",
    "            print '---  %s ---' %s\n",
    "            np.savetxt(barcode_dir + s, D[s]['cell_index'], fmt='%s')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "MT = 1\n",
    "\n",
    "if MT:\n",
    "    MT_dir = main_spring_dir + \"/mito/\"\n",
    "    if not os.path.exists(MT_dir):\n",
    "        os.makedirs(MT_dir)\n",
    "        for s in sample_list:\n",
    "            print '---  %s ---' %s\n",
    "            np.savetxt(MT_dir + s, D[s]['mito_frac'], fmt='%s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_list_new=np.array(merged_list)\n",
    "np.savetxt(main_spring_dir + '/genes.txt', merged_list, fmt='%s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving hdf5 file for fast gene loading...\n",
      "AL627309.1_ENSG00000238009 0 / 20922\n",
      "TANC1_ENSG00000115183 3000 / 20922\n",
      "PCDHGA1_ENSG00000204956 6000 / 20922\n",
      "PPP1R3B_ENSG00000173281 9000 / 20922\n",
      "KCNMA1-AS3_ENSG00000225652 12000 / 20922\n",
      "USP3-AS1_ENSG00000259248 15000 / 20922\n",
      "C18orf65_ENSG00000275763 18000 / 20922\n",
      "Saving hdf5 file for fast cell loading...\n",
      "0 / 7902\n",
      "3000 / 7902\n",
      "6000 / 7902\n"
     ]
    }
   ],
   "source": [
    "# save master expression matrix in hdf5 format\n",
    "import h5py\n",
    "\n",
    "print 'Saving hdf5 file for fast gene loading...'\n",
    "E = E.tocsc()\n",
    "hf = h5py.File(main_spring_dir + 'counts_norm_sparse_genes.hdf5', 'w')\n",
    "counts_group = hf.create_group('counts')\n",
    "cix_group = hf.create_group('cell_ix')\n",
    "\n",
    "hf.attrs['ncells'] = E.shape[0]\n",
    "hf.attrs['ngenes'] = E.shape[1]\n",
    "\n",
    "for iG, g in enumerate(merged_list):\n",
    "    if iG % 3000 == 0:\n",
    "        print g, iG, '/', len(merged_list)\n",
    "    counts = E[:,iG].A.squeeze()\n",
    "    cell_ix = np.nonzero(counts)[0]\n",
    "    counts = counts[cell_ix]\n",
    "    counts_group.create_dataset(g, data = counts)\n",
    "    cix_group.create_dataset(g, data = cell_ix)\n",
    "\n",
    "hf.close()\n",
    "\n",
    "##############\n",
    "\n",
    "print 'Saving hdf5 file for fast cell loading...'\n",
    "E = E.tocsr()\n",
    "hf = h5py.File(main_spring_dir + 'counts_norm_sparse_cells.hdf5', 'w')\n",
    "counts_group = hf.create_group('counts')\n",
    "gix_group = hf.create_group('gene_ix')\n",
    "\n",
    "hf.attrs['ncells'] = E.shape[0]\n",
    "hf.attrs['ngenes'] = E.shape[1]\n",
    "\n",
    "for iC in range(E.shape[0]):\n",
    "    if iC % 3000 == 0:\n",
    "        print iC, '/', E.shape[0]\n",
    "    counts = E[iC,:].A.squeeze()\n",
    "    gene_ix = np.nonzero(counts)[0]\n",
    "    counts = counts[gene_ix]\n",
    "    \n",
    "    counts_group.create_dataset(str(iC), data = counts)\n",
    "    gix_group.create_dataset(str(iC), data = gene_ix)\n",
    "\n",
    "hf.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "scipy.io.mmwrite(main_spring_dir + 'matrix.mtx', E)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Save merged samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Sample1']\n",
      "1612984906.71\n",
      "________________ /FullDataset_v1_rna\n",
      "Filtering genes\n",
      "Using 986 genes\n",
      "Excluded 17 genes\n",
      "Running PCA\n",
      "Building kNN graph\n",
      "kNN built in 5.93095 sec\n",
      "Saving SPRING files to /home/I0369218/precision-immunology/CHAMBERLAIN.Mat/CITESEQ_EXPLORATORY_CITESEQ_5K_PBMCS/FullDataset_v1_rna\n",
      "Running ForceAtlas2\n",
      "190.111049891\n",
      "190.200725079\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAGPCAYAAACnLVLPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xt8U+X9B/DPN23TtNA2UJCmKgXxLigT0HmbU7b9dHPTTdztVzZ3syi6sTkF5wU75pTpNhTR1t1/sM055hXFqXjFObVVoXgXpIptgRbaUto0bfP8/khS0+QkOSc9yclJPu/XK6/Cybk8Sdp8z/M83+d5RCkFIiIisieH1QUgIiKi5DGQExER2RgDORERkY0xkBMREdkYAzkREZGNMZATERHZGAM50SiIyPUioiK2PS0iT1tUJCLKMQzklPVEpFlEaoP/flpEVIzHJ1NYhvDrDIrIHhFpFJFbReToUZ77GhE5z6yyjlY6yiMiJ4jIKhF5WUT6g+9rRSqvSZSp8q0uAFEqicgMAJMBPBy2uRXAlRq7v5fi4jwJ4E8ABEAZgOMAfBvAQhFZrJT6dZLnvQbA3QDuN6WUo5eO8nwewEUAtgB4B8D0FF6LKKMxkFO2+wKAXQBeDtvWrZRaY0FZ3o28rogsAfAQgFtE5E2l1CMWlMuO7gSwXCnVJyLXw+RALiLFSqleM89JlCpsWqdsdw6A9crAXMQicmGwqXZKxPYpwe0XmlU4pVQHgK8DGARwddi1nCLycxFpEJG9ItInIi+KyJciyqQAFAL4dljT/dPB58aLyC0isllE9olIj4g8JSKnaLzmrwabqbtFpEtEmkTkuoh9SoPn2x5szv5ARG4SkUI95UlERNYFu0EkxnMfhJ5TSu1USvXpOa+O6z4tIm+JyLEi8qSI7AdwR8Rzh4nIoyKyX0R2isiNIuKIOM9PReR5EWkXEW/wPfyeGWUkioc1cspaIjIewCcBrIh4Kk9EJkRs8ymlutNTspGUUh+IyDMAzhCR0mA5SgEsAPAPAH8EUATgfwHcLyJnK6X+HTx8fvD5/wK4K7htZ/DnIQDmAViLQLeBG8D3AWwQkdlKqS0AICKfQaAp/EkAVwEYAnAEgE+FyigiRQCeBjA1eJ2tAGYCuBzAkQBCfeLxypPI3Qi0oJwM4Pmwa48D8DkAtxq5ITPIDeAxAPcGy9EZ9lwZgCcArEOgu+AsAEsAvI+PXyMA/BjAegD3AFAAvgzg9yKSp5QK34/IXEopPvjIygeAbwLwASgN2/Y0Al+ykY91YftcGNw2JeJ8U4LbLwzbdn3gz2jEfk8DeDpimwJQF6esK4L7HBv8fx6Awoh9nADeAPB4xHYvgD9rnLMQQF7EtvEIdDX8LmzbbwF0Re4bcdzPAPQCOCpi+yXBcp+WqDw6Pq+xwWusjNj+/eA1ZsU47vrg8xVJ/p6Efid+GOe570dsfw3AyxHbijWOfxyBLhXL/x74yN4Hm9Ypm50DYKOKrml/COCzEY+rYa2e4M8SAFBKDSml+oHhZvbxCNTSnwEwS88JlVL9Sqmh4DlcIlKOQHfaSxHn6AYwBoFabyxfRaCWvFtEJoQeCNRUAeBMPWVKUN4eBJIS50U0W38NwHtKqcbRXiOOQYysXYfzIpCkGO4ZBFo8hqlgn7qIFAS7NSYg0MpxqIiUmVxeomFsWqesJCJ5AP4HwA0aT/cqpZ7Q2G6lscGf+0IbROT7CDTXHoVApnuIrublYDC8EoHs7qkRT78f9u87EQjUj4hICwLB+V4ADyqlQtc6HIEs+90xLneAnjLpcDcC3QGfBvCkiBwA4AwAN5p0/lhalFLeGM/tCN0QhdmLQOvGMBE5F8C1CHQ55EXsX4ZAqweR6RjIKVudhMAX7cOJdtQQK1BGfjmbaToCfdPvA4CIfAPA7wA8CGA5As3hgwC+g0CXgR6LAfwSwF8QGBLWEbzGVQCmhXZSSrWJyHEItEycFXx8C8B6EflCMJg7EKhdxgqoO/S+0AQeRqCF4OvB612AwPt+t0nnjyVe4lxkEI8STCC8D8BGADUIDHH0ITBM7sdgYjGlEAM5ZatzAGxVSr2dxLF7gz/dEdunjKpEMYjIZACnA/ivUipUI/8agG0AzgurFUNEvqNxilg3Hl9DoK/+wojr1UadQCkfAkH04WBm+I0I3AiEEs+2AijR2ZKRdEKaUsorIg8A+IqILEQgoG9RSr2e7DnT5AIEmuA/F16zF5EzrCsS5QreJVK2+gICWcbJCE0ME/klvDD54mgL9n3/HYGb6l+EPRWqBUrYvocgkAkdaT+AcRrbhzCySR4icjICrRXh28rD/x+8cXg1+N/Qee8GMCdy+Fvw+CIRGRu2KVZ59Po7gHIEkg5PQepr42YYQuAGZvg7NZht/13LSkQ5gzVyyjrBGu50AD9J5nil1BsishHADcFAuxPAlxDRJ5qEw0SkGoHgWopAn/MFCPSP/1gp9WjYvg8C+AqAB0XkQQAHIpAh/k7wuHANAD4jIj9FoIl7l1LqyeA5rheR/wPwHIDDEOgvfwMf98kDgSFSEwBsQCAR8EAAlyLQPPxMcJ9bELg5uldE1iCQMFeIwDC1ryLQhPzfBOXR63EEugF+HXyvogK5iFQhMNQN+HiY3CIR6QHQrJRabeB6ZngIgd+3x0VkNQK/Kz8A0AaAU8dSalmdNs8HH2Y/EAh4+wA4NZ57GsBbOs5RhcCY4F4A7QBWATgGoxt+FnoMIdB8/woCw86OjlGGKxBoXvciMBVpdYzrHYFAf3JP8PxPB7c7Eehf34FAH/DLCPR//xnA9rDjzwfwKAJBpx/AB8F9pkZcpzh4/beC+3UEz3kdAHei8hj8DOuCx74c4/lPR7yn4Q9D14v3OxHruRifw7cAvBn8vN5DILB/BxpDGfngw8yHKJWq+RWIrCEiDwPoV0p9xeqyEBGlGpvWKRs9g7CZwYiIshlr5ESUcpJ4iVGfUmqPydeciPhDBoeUUrHGxRPZBgM5EaVccDGVeJ5RSn3a5GtuRyDXIZZmpdQUM69JZAU2rRNROnw2wfN7EzyfjP9FYLGZWExZPY3IaqyRExER2ZgtauQTJkxQU6ZMsboYREREadHY2NiulJqoZ19bBPIpU6agoaHB6mIQERGlhYg0692XU7QSERHZGAM5ERGRjTGQExER2RgDORERkY0xkBMREdkYAzkREZGNMZATERHZGAM5ERGRjTGQExER2RgDORERkY0xkBMREdkYAzkREZGNMZATEVFWqK+vR2VlJRwOByorK1FfX291kdLCFqufERERxVNfX49FixbB6/UCAFpbW7Fo0SIAQE1NjZVFSzlRSlldhoRmz56tuIwpERHFUllZidbW1qjtHo8HLS0tFpRodESkUSk1W8++bFonIiLba2trM7Q9mzCQExGR7VVUVBjank0YyImIyPaWLl0Kl8s1YpvL5cLSpUstKlH6MJATEZHt1dTUYMWKFfB4PBAReDwerFixIusT3QAmuxEREWUcJrsRERHlCAZyIiKyXK5O5mIGTghDRESWyuXJXMzAPnIiIrJUtk3mYgb2kRMRkW3k8mQuZmAgJyIiS+XyZC5mYCAnIiJL5fJkLmZgICciIkvl8mQuZmCyGxERUYZhshsREVGOSFsgF5E8EVkmIu+LiDf48xciwrHsRERESUpnEF0MYCGAbwNoAnAsgL8A6AewLI3lICIiyhrpDOQnA3hIKfVQ8P/bReRBACemsQxERERZJZ195BsBnCEiRwKAiBwN4EwAj6SxDERERFklnTXy5QBKALwhIkPBa9+glLpDa2cRuQjARQAwefLktBWSiIjITtJZI/8agG8B+CaA44P/vkREvqe1s1LqLqXUbKXU7IkTJ6axmERERPaRzhr5zQBuUUrdHfx/k4hUAbgKwB/SWA4iIqKskc4aeTGAoYhtQ2kuAxER0bBsWAc9nTXyhwAsEZH3AbwO4BMAfgLg/9JYBiIiIgDZsw562qZoFZESBMaLfxnAAQBaAdwN4OdKKW+8YzlFKxERmS2T10E3MkVr2mrkSql9ABYFH0RERJbKlnXQ2T9NREQ5KVvWQWcgJyKinJQt66AzkBMRUU7KlnXQuR45ERFRhuF65ERERDmCgZyIiLJWNkz4kkg6J4QhIiJKm2yZ8CUR9pETEVFWyuQJXxJhHzkREeW8bJnwJREGciIiykqxJnYRkazqK2cgJyKirKQ14QsA+P1+LFq0KGuCOQM5ERFlpdCELw5HdKjzer2ora21oFTmYyAnIqKsVVNTg1hJ3dnSV85ATkREWS1bFkeJhYGciIiyWrYsjhILAzkREWW1bFkcJRZOCENERJRhOCEMERFRjmAgJyIisjEGciIiIhtjICciIrIxBnIiIiIbYyAnIqK0q6+vR2VlJRwOByorK7Nm3nMr5FtdACIiyi319fVYtGgRvF4vAKC1tRWLFi0CgKwZ251OHEdORERpVVlZidbW1qjtHo8HLS0tFpQo83AcORERZaxYi5VkyyIm6cZATkREaZXti5ikGwM5ERGlVbYvYpJuDORERJRW2b6ISbox2Y2IiCjDMNmNiIgoRzCQExER2RgDORERkY0xkBMRkW6cWjXzcIpWIiLShVOrZiZmrRMRkS6cWjV9mLVORJRCdmleNrucnFo1M7FpnYjIALs0L6einBUVFZo1ck6tai02rRMRGWCX5uVUlDPy5gAITK2a67OyDQ754VeAM9+8Rm42rRMRpYhdmpdTUU5Oraqt7pmtOHfV89jnHbDk+gzkRGRbVvRV22XlrlSVs6amBi0tLfD7/WhpaUk6iNslzyCRt9q6ceuGd3HIxDEocRVYUgYGciKypVAzb2trK5RSw33AqQ4Idlm5K5PLadVnZ7aBIT8uv2cTyooKsOzc6dYVRCmV8Y9Zs2YpIqJwHo9HAYh6eDyelF+7rq5OeTweJSLK4/Gourq6lF8zGZlaTis/OzPd+sQ7qmrxOrW+qcX0cwNoUDpjJJPdiMiWHA4HtL6/RAR+v9+CEpFe2fDZvdnajS/dvhFnTfdg5Tc+Yfr5mexGRFnPLn3VFM3un93AkB8//WegSb32S8dYXRwGciKyp0zuA6b47P7Z3fHUVrze0o1fnDcD48c4rS4OAzkR2ROHQtmXnT+711u6sPLJd/Gl4ypx1vTMaEFgHzkRZaT6+nrU1taira0NFRUVWLp0qS2+6Cl7+Qb9OHfV89i9rx+P//hTGJfC2riRPnJO0UpEGccu06BSbln11Ht4s7Ubd82fldIgbhSb1oko49TW1o6YBhQAvF4vamtrLSoR5brXW7qw6qn3cN7MSnzumMxoUg9hICeijGOXaVApN/gGAxO/jBvjxPUZkKUeiYGciDKO3YcnUXa5/an38FbbPtz45RlwF2dOk3oIAzkRZRy7D0+i7LHloy7c8dR7+MonDsRnjp5kdXE0MZATUcYZzfCkbFmMg6znGwxM/DJ+jBNLv5h5TeohzFonooxUU1NjOEOd2e5kppVPvou32vbhD9+ejbJia1Y204M1ciLKGsx2J7Ns3tGJO57eivOPPwhzj8rMJvUQBnIiyhrMdicz9A8O4af/3IQJY5247otHW12chBjIiShrMNudzHDrE+/inZ09uOkrx6KsKHOb1EPSGshFxCMifxGR3SLiFZE3ROT0dJaBiLIXs92NYWJgtE0fdqLuma2YN+sgnHHkAVYXR5e0JbuJiBvA8wA2AvgCgN0ADgGwK11lIKLsFkpo4xztiTExMJp3INCkfkCJC9eek/lN6iFpWzRFRH4J4HSl1ClGj+WiKUTpw8VKckNlZSVaW1ujtns8HrS0tFhQIustf/Qt3Pn0VvzpO3NwxhHW1saNLJqSzqb18wC8KCL/EJFdIvKaiFwqIpLGMhBRHKFaWmtrK5RSw7U0NrlmHyYGjvTqB3tR/8xWfHX2QZYHcaPSWSMPjQn5LYB7AMwEsBLAEqXU7Rr7XwTgIgCYPHnyrObm5rSUkyiXsZaWO/hZf8w7MIQv3PYcen1D+PePP4VSl/UJbplaI3cAeEUpdZVS6lWl1J8A3AZgodbOSqm7lFKzlVKzJ06cmMZiEuUu1tLSz6qEMyYGfuy3T7yDrbv346bzj82IIG5UOgN5K4A3Ira9CWByGstARHFw+FZ6WdmVMZppcLPJKx/sxe+e3YavzzkYpx9uz0pjOpvW/wbgYKXUaWHblgE4XykVNz2QyW5E6RGZyQwEamm5+AWfDmzetpZ3YAifv+059A/48eii01CSQbXxTG1a/y2AT4rI1SJyqIhcAOCHAFalsQxEFAdraenFrgxr/ebxd7Bt934sP//YjAriRqWtRg4AIvIFAL8EcASADwDcDmClSlAI1siJKBuxRm6dxuY9mFf3Ar5xwmT88sszrC5OlEytkUMp9bBS6jillEspdbhS6rZEQZyIKFsx4cwa3oEhXPHPzagsK8LPPn+U1cUZNc61TkRkEXZlWOOWf7+Nbe37cfO8YzG20P6reTOQExFZqKamBi0tLfD7/WhpadEM4pwT3TwN2/fgD8+/j+pPTsbJh06wujimsP+tCBFRFuOc6Obp8wXmUj9oXBGuOtv+TeohrJETEWWw2traEcMBAcDr9aK2ttaiEtnXr/79FrZ39OJX5x+HMVnQpB7CQE5ElME4RM0cL27rwJ+e345vn1SFk6aVW10cUzGQE+U49r9mNs62N3q9vkFcsXYzJo8vxuKzj7S6OKZjICfKYVztLPNxiNroLV//Fj7Y04ub5x2LYmf2NKmHMJAT5TD2v2Y+DlEbnRe2duAvLzTjO6dMwYmHZFeTekhaZ3ZLFmd2I0oNh8MBre8AEYHf77egRETm2d8/iP9Z8SzyHYL1P/oUipx5VhdJt4yd2Y2IMgv7Xymb3bj+TXzU2YebLzjOVkHcKAZyohzG/lfKVs+/1441//0A3ztlKuZMGW91cVKKgZwoh7H/lbLRPu8Arly7GYdMGIOf/s8RVhcn5bIvfY+IDKmpqWHgpqzyy0feQmtXH/654GS4CrK3ST2ENXIiIsoaz727G39/6QN8/7RDMKtqnNXFSQsGciIiygrd3gEsXrsZ0yaOwU8+e7jVxUkbNq0TEVFWuGHdm2jr9uJfF+dGk3oIa+RERGR7T729C/9o+BA1p0/DJybnRpN6CAM5ERHZWlffAK76VxMOnzQWiz5zmNXFSTs2rRMRka0tW/cGdvf0465vzUJhfu40qYewRk5ERLa14c2dWNu4A5d8ehqOPchtdXEswUBORES21Nnrw1X3NuHIihJcdmbuNamHsGmdiIhsqfahN7Bnvw9/vHAOnPm5Wy/N3VdORES29e/X23Dfqx/h0jMPxfQDy6wujqUYyImIckx9fT0qKyvhcDhQWVmJ+vp6q4tkyJ79Plx9XxOO9pRi4RmHWl0cy7FpnYgoh9TX12PRokXwer0AgNbWVixatAgAbDPn/tIHX0dX3wBWf+9EFOSxPsp3gIgoh9TW1g4H8RCv14va2lqLSmTMI02teGhTC3409zAc5Sm1ujgZgYGciCiHtLW1GdqeSdp7+nHN/Vsw48AyLDh9mtXFyRgM5EREOaSiosLQ9kyhlMK1929Bj3cQv/7qcchnk/owvhNERDlk6dKlcLlcI7a5XC4sXbrUohLp89DmVqzf0oYff/ZwHD6pxOriZBQGcqIsZ/cMZTJXTU0NVqxYAY/HAxGBx+PBihUrMjrRbdc+L657YAtmHuzGD06banVxMo4opYwfJOIA4FJK9ZpfpGizZ89WDQ0N6bgUUVaJzFAOqa6uxurVqy0qFZF+SilctLoRz7yzG4/88DQcesBYq4uUFiLSqJSarWdf3TVyETlbRP5PRJoB+ADsE5H9IvKciFwtIpXJFpiIUkMrQxkA1qxZw5o52cJ9r36Ex9/YiSs+d0TOBHGjEgZyETlPRN4B8EcAAwB+CeDLAP4HwHcAPAXgMwC2iUidiExMYXmJyIB4mch2GW5Euauty4vrH3wds6rG4bunskk9Fj018qsA/ATAgUqp7yml6pVSDymlnlBK3aOUuk4pdQaAaQC6AXwrlQUmIv3iZSLbYbhRJPb35w6lFK66dzN8Q37ccsFxyHOI1UXKWAkDuVLqRKXUOqWUP8F+HymlrlRK/dq84hFRSDJBLF4mcqYPN4oU6u9vbW2FUmp4RjIG8+z0z4YdeOrt3Vh81pGYOmGM1cXJaLr6yEWE0+cQWSjZIFZTU4Pq6uqo7XYYbhTJ7jOSmSUXWiU+6uzDsnVv4MSp4/Htk6ZYXZzMp5RK+AAwBOCAsP/fAWCCnmPNeMyaNUsR5TKPx6MARD08Ho+u4+vq6pTH41Eiojwej6qrq0txic0nIprvgYhYXbSkGf1c6urqlMvlGvH6XS6XLT/PWPx+v6r+/X/VUdeuV83t+60ujmUANCidMVLX8DMR8QOoUErtCv6/G8BMpdQ2E+8pYuLwM8p1DocDWn+rIgK/P26vV9aorKxEa2tr1HaPx4OWlhYLSjQ6WkMDXS5X3DHd2fYeaPnri824+r4tWHbedMz/ZJXVxbFMSoafRV4jyeOIKAl2nVbTTHadkSyWZLoK7DxPuh4f7unFDQ+/iVMOLUf1iZOtLo5tGAnkxmeOISJTZEsQG03/rh1nJIsnmaCczTd0fr/CFWs3wSGCX807DiKsL+qmp/0dgB/A7wHcFnz0A/hL2P9vA3Cb3vZ8ow/2kRPZv587F/p3jUgm7yGb38M/bdymqhavU3e/1Gx1UTICUtBH/jQS18iVUupMQ3cROrGPnMj+cqF/14hk+shDx9XW1qKtrQ0VFRVYunSpbVslQt5v34+zb30WJx1Sjj9eOIe1cRjrI09qrvV0YyAnsj8m7EXLxqBs1JBf4av1L+Ddnfvw2I9PR0WZK/FBOcBIIM/XecJSpVT36IpFRLmsoqJCs0aeDf27yaqpqcm5wB3pjxvfR2PzXvzmq8cxiCdJb7LbXhE5IPQfEblDRCakqExElIWyJWGPzPPern24+bG38dmjJ+HLnzjQ6uLYlt5AHtlhUQ2As70RkW7ZlnVOozM45Mfl92zCGGcefvnlGewXHwWOIyeitKmpqUFLSwv8fj9aWlqyOojnwlSqo1H/7DZs2tGFZedNx8SSQquLY2u6+siDMj8rjogoA0RmpIfmxgeQ1Tcver3Z2o0VT7yDLxzrwTnHVlpdHNszMkXrHwH0BjfVALgbQFf4fkqpH5pdQIBZ60RkLxxqF5tv0I9zVz2P3fu8eOzHp2P8GKfVRcpIpmetA3gWgfXGQ/4DIHL+PNbYiYiQ/VOpjsbtT76LN1u7cdf8WQziJtEVyJVSn05xOYiIsgaH2mnbvKMTq57eiq984kB87pjcfi/MlGyyGxERxcChdtG8A0O4/J5NmDi2EEu/eIzVxckqCQO5iFwjImP0nExEThGRL46+WERE9sWhdtF++/g7eHdXD246fwbKigusLk5W0dO0fiiAD0TkXwAeRGAi9zYAEBEXgKMBnIrA2PJyAN9OUVmJiGyDs7Z9rLF5D+56bhu+ccLB+PQRByQ+gAxJGMiVUheKyAwAlwJYDaBURBSAAQBOBMaUvwLgLgB/Vkr5UlheIiKykV7fIC6/ZxMOdBfh6i8cbXVxspLeZLcmADUicjGAYwFUASgC0A7gNaVUe+qKSEREdvWrR9/G9o5e/O0HJ2JsoZGpS0gvQ8luSim/Uuo1pdQDSqm7lVJPMIgTWYsziFGm+s977fjzf7bjwpOn4ORpXJ4jVQzdHolI5NjxEAXAq5TaPfoiEZFenEGMMtU+7wCuWLsZUyeMweKzjrS6OFnN6PCz7QDe13hsB9AmIntF5DcikvAGQUR+JiJKRG43WAYiCqqtrR0O4iFerxe1tbUWlYgo4IaH30RrVx9uueA4FDnzrC5OVjPaYfENAL8CUAfgxeC2EwFcBOB6AG4A1wDYByDmgEkR+SSAHwDYbPD6RBSGM4hRJnrqrV24++UPcfGnp2FW1Tiri5P1jAbyiwH8WCl1b9i2J0XkbQA/UkqdLiK7ANQiRiAXkTIAfwXwPQDXJVFmIgriDGKUaTp7fVj8r804YlIJFn3mMKuLkxOMNq2fCKBJY/sWAHOC/34BwEFxznEXgLVKqSfjXUhELhKRBhFp2L2bXe9EWjiDGGWapQ++jj37ffj1V49DYT6b1NPBaCBvRqAZPdIPAHwQ/PdEAHu0DhaRHyAwwcy1iS6klLpLKTVbKTV74sSJBotJZI5MzwjnDGKZLdN/f8z2SFMrHnitBZedeRimH1hmdXFyhtGm9csB/EtEPg/gZQSy1ecgsDLa+cF95gC4J/JAETkCwC8BnMZJY8gO7JIRzhnEMpNdfn/MsmufF1ff14RjDyrDJWdMS3wAmUbXeuQjDhA5GMAlAI5AYFa3NwHUKaU+SHDchQD+BGAobHMeAjcDfgBjlFL9WsdyPXKyAteUptHIpd8fpRR+8H8NePbddjzyw1Nx6AElVhfJ9lKxHvkwpdSHAK4yXCrgfgCR0fhPAN5FoKbOWjplFGaE02jk0u/P2sYdeOLNXbjmC0cxiFvA8DKmIjJJRH4uImtF5J8icr2ITEp0nFKqUym1JfwBYD+APcH/G2saIEqxWJnfzAjPXmb2aWf7709HRweamprw6HMvYekDTTj+oBJ895SpVhcrJxkK5CJyCoD3AHwTQB8ALwKrnr0rIieZXzwi6zAjPLeE+rRbW1uhlBru0042mGfz709HRweam5vh7e/Hqpe74FfAD451Ye9ezTxnSjFDfeQi8gICw88WKKX8wW0OBCaIma6UOjkVhWQfOVmlvr4etbW1aGtrQ0VFBZYuXZqViUqUmj7tbP39aWpqgs/nwyPv7scfXtuHBbNK8dlDiuF0OjFjxgyri5cVjPSRGw3kfQBmKqXejth+JIBXlVJFhkqqEwM5EaWaw+GA1vehiMDv91tQoszV2NiIln2DuPzxdhwz0YmrTx0HEQEAzJo1y+LSZQcjgdxoH3kXAK1OkKkAOg2ei4goY2R7n7aZ8vILcNtLXXA6BJfMLhsO4k6n0+KS5SajgfxuAH8Qkf8VkakiMkVEqgH8DsDfzS8eEVF6ZHOfttke+xB4d8+VRg76AAAgAElEQVQAfnB8KcYXBWZvExFUVlZaXLLcZDSQXwlgLYA/IpD0tg3A7wH8E8ASc4tGRJQ+emfJy7XZ2iK93tKF372wA589YjzOPDQwe5vT6URVVRXKy8stLl1uMjwhDACISDECs7kJgPeUUr1mFywc+8iJKBNEztYGBGrtuTItbv/gEL608nns6fXhsUWfwrgxbEpPFVP7yEXkwcgHAk3sNwD4BYC7w7YTZSWtWliu18xyUa6v//6bx9/B2zv3Yfn5MxjEM4iemd06Ul4KogymNWf2ZZddBqUUBgcHh7dl8zzaFJBLs7VFenn7Htz17DZ844SDceaRCecAozRKqmk93di0TlaKNb5YSzbOo00fy6X508P19A/i7FufBQCs/9GnMLbQ8OzeZFAqh58R5Rwjta1cqJnlsniZ7dnc1XLDw29ix94+/PqCmQziGYiBnCgBI+OIOeY4u8XKbAdg6vSumeSpt3bh7y99gIs+dQhOmDre6uKQBjatEyWglalcUFAwoo8cyK3sZRopW5vc9+734XMrnsX4YicevOwUFObnWV2knMGmdSITadXCVq5cidtvvz3hmGPKDdmYBKeUwtX3N6Gz14fffO04BvEMxkBOpENNTQ1aWlrg9/vR0tKCmpoazW12lc39u+mQjdO73v/aR3ikqQ0//uzhOKayzOriUBwM5GR7DEKjY/bynbkomeldM/n3tqWzD9c98DpmV41DzaemWV0cSkQplfGPWbNmKSItdXV1yuVyKQDDD5fLperq6pI6l9vtHj5PWVlZUuexG4/HM+L9Cz08Ho/VRUururo65fF4lIgoj8dj+LM3cryZv7dmGxryq2/c9YI66tr1ant7j9XFyVkAGpTOGGl5kNbzYCCnWMwKQnV1daqgoCDqPPn5+Rnx5ZpKIqL5HoqI1UVLm3QH1ky+efrDc9tU1eJ16m8vNltdlJxmJJAza51szaw1pONN+mL3zONEsjXj2oh0vweZuvb5e7v24Qu3bcSph07A7789e3h5Uko/Zq1TzjAryShedrGdM4/14PKd6c86z8TkON+gH4v+8RrGFObjpvOPZRC3EQZysjWzglC8L1Arv1zTkRCld/nObJbuwJqJN0+3bXgXWz7qxo1fmYGJJYWWlYOSoLcN3soH+8gpntEmKYXOkZeXp9lvWV1dnYJS6ytTpiZEZRsr3mszfm8TaW9vV5s3b1YNDQ1q8+bNqr29XXO/hu0dauqSdeqKf75mehkoOWAfOZFxbrcbXV1dUdut6itm33V61dfXo7a2Fm1tbaioqMDSpUtt3SrR0dGB5ubmEX3xIoKqqiqUl5cPb+vpH8Tnb30OCooLomQQ9pETJaG7u1tzu1V95Nk4W1gmS/UEP+keN97S0hKVUKeUiroJXPbQG9ixtxe/+SoXRLErBnKioExLQMq08mQLKyZisWLSHZ/Pl3D7Y6+34R8NH2LB6dMwZwoXRLErBnKioExLQMq08mQDq2axq62tHbHoDgB4vV7U1tam7JpOpzPu9l37vFhybxOO9pRi0WcOT1k5KPUYyImCMi17O9PKkw2sCKiANd0klZWVUUPIRASVlZVQSmHx2s3Y3z+IW78+E858hgI7Y7IbEeUMqyZisSpxsaOjAy0tLfD5fHA6naisrER5eTlW/7cZ196/Bdd/8WhceMrUlF2fkmck2Y2ZDUSUMyoqKjQDaqrzDpYuXRq1pn06uknKy8tHZKgDwHu7enDDw2/gU4dPxLdOmpLS61N6sD2FiHKGVXkHmdJNEpi97VUUFeTh5nnHwuHg7G3ZgDVyIsoZocBpxXjx0Br2Vrp1wzvY8lE36qqPx6RSV+IDyBZYIyeinJLsePFMXj9cjxe3deCOp7figlkH4azpHquLQyZiICfKYnYPPpnCqmFrZunqG8BP7tmEyeOLcf2XjrG6OGQyBnKiLGX34JNJNyFWDVszy7X3b0FbtxcrvjYTYzh7W9ZhIKeckknBIdXsHHwy7SbEztPl3v/qR3hwUwsWzT0Mn5g8zuriUAowkFPO0BscsiXY2zn4mHETYubnaNfpcj/c04tr79+C2VXjcMkZh1pdHEoVvcukWfngMqZkBo/Ho7lMqcfjGd4nm5YO1fN6M5WIaJZdRHQdb/bnaMffi4HBIXX+Hc+r6dc9qj7o2G91ccggGFjG1PIgrefBQE5m0BMc7Bz8Itkx+ISM9nOIdbzD4RhVME/1+uFm+u3jb6uqxevUfa/ssLoolAQGciINeoLDaGuCmcZuwSdktDchsT5HO93MjMbL73eoqUvWqUV3v2p1UShJRgI5+8gpayTqE9Uzq5dd+0JjSfUa26ky2pnQ4n1eevra7Zwn0e0dwI/ufg0HjivCz8/lULOcoDfiW/lgjZwS0VuDS1RDtXNzNH1M63OEzhYWO/8O+P1+ddnfXlGHXPWwamzeY3VxaBTApnXKNWb2bdu1OZoCQp9frCCe6PfCznkS/2r8UFUtXqdue+Idq4tCo2QkkLNpnVKuo6MDTU1NaGxsRFNTEzo6Oky/hllDrerr6y2ZhzsT2bF5OXyIYSyJFkmx67C97e37ce39W3DClPFRQ83S8TdI1mEgp5Tq6OhAc3MzfD4fAMDn86G5udn0LxIz+rYzbRISK2XCe5HMjYTW+HMgsA653r52O+ZJ+Ab9+OHdryI/z4EVX5+JvLBVzdL1N0jWYSCnlGppaQn04YRRSqGlpcXU65ixPKWdZ0Izm9XvRbI3ErFqzUop3Ql/Vi11Ohq3PPY2Nu/owvLzj0Wlu2jEc+n6GyTrMJBTSoVqAXq3J8uM9Z7t2qSaCla/F8neSJhRm86UtcP1euad3bjr2W2o/uRknDU9+nWm62+QrMNATinldDoNbQeS75sd7VArOzaphpjdn230vTD7+sneSJhVm7bLsL1d+7y4/J7XcMSkElzzhaM190nmb5DshYGcUqqyshIiMmKbiKCyslJzfyNNqmYHDzs2qQKp6c828l6k4vrJ3lRlQm06lUmC4UlrmzZvxg//2oB93kGs/OYn4CrI0zzG6N8g2ZDe9HYrHxx+Zm/t7e1q8+bNqqGhQW3evFm1t7fH3Ffv0J9UjfW149CzVA2X0vtepOL6Vo3lHu3nn8pyt7e3q8bGRtXQ0KAaGhrU1WueVlWL16n6J17Xdazev0HKDOA4crIrvVOk2nmsr9mSmVbWzBuWVE1rW1dXp8rKyobP53a7UxrIzQjCqfy9DAXihoYG9dd//0dNXbJOfX3lE2rTpk2jPjdlHiOBnE3rlFA6x6DqbVKN1Vfa2tpqmzHPZkmmP9vMpvBU5hb09/cP/7uzszOlQ+DMyNRPZZJgKDmtx+fHb/7biQnFebh4dikGBgZGfW6yNwZyiivdY1Dj9c2G9z1G9vmFy7Xx30b79kcTsLT6f1OVW5DuIXDxgrDefu9U3tQ4nU4opXBHQxf29vnxk0+WYUyBg0lrxEBO8SUag2p2bT1WshKAEbVIv98f9zxGvvDtOINZOKMJXsnWGmPV5AGkJMFsNOVM5vOMFWzLysp0t2CkMmGytLQUj27txYsf9aN6RgkOG+9k0hoF6G2Dt/LBPnLrhPrktB6RyTcNDQ2qsbExJYk0iebO1nro6aO18wIZyUq2HzeZ4yL74qurq3X3zSd7vWQ/z1jHhvfT6ylHKhIm29vb1T1P/FdNu2qd+sqKx9VLL7+sGhoa1Pbt20d9bspMYB85mSXeGNRYtfVU1HCT6WPU05yZzubbeO9LOlsFkq01Gq0h19fX49JLLx1Rk12zZo3uvvlkyrl48eKkP89YLRvd3d2a+8d63akYg/7u9h245T97UVbowKVzyuAIdi3FKhvlGL0R38oHa+TWiVfr1qqlL1u2LCU1XKM1cr3XTFXGdaTq6uqYZbSiVSCZWqPRGnKsmqze442Ws66ublStM2a9brP5/X71tdseV1OXrFNrHv1P1N8cZSdw+BmZpb29Xb366qvDXxqvvvrqcNN5+HCY0KOioiIlX3qJ1pcuKytLqjkzHV/S8QKMx+OxPFDoZfSGQ+9NV2SQTbZp2u12J3WzYPbrNtufNm5TVYvXqWv++kzU39vmzZvTUgZKPyOBnE3rOcRoYlooY31oaGh4W3iSmdaMUTt37tQ812iH34SaPd1ud9RzLpcLy5cvj2rO1NNcnY7Z3OI167a1tVk+r7leqZo1LbwLRCuhbsGCBXC73Qm7Gzo7O2M+19fXl3S3hZWzxW3e0YkbHnkTp01z49wjxox4joluNExvxLfywRr56CWTmKZV446sBUTOGGVkZrZkE4L0HGukFmW0+dZouWM134feF7vUyI2KV0OO9ZnE60JJVAtOdC2t82TyTH6d+33q1OUb1Mk3blB79/dzdrYcg0xsWgdwFYCXAXQD2A3gIQDT9RzLQD56eoJypHgZ67HoCaDpaKrMpGlD4wUnq/rI06Gurk4VFBSMeF0Oh0O53e6YgTPeTU+iz09vn3zoPJn8vg8N+dX3/vySOvRnD6vG5j1WF4cskKmB/N8AvgNgOoAZAO4D0AZgfKJjGchHL5mgnEzwVypxLScdNdBUJLElW+5Y/fvV1dUj9snUmqERoxluppS+pMZ4187Pz9cdzB0OR8a2hKx66l1VtXid+tPGbVYXhSySkYE86sLAWABDAL6YaF8G8tFLJiinapx4oqZmM4JYKm4WRnNzkK2BOlFLS0FBgSorK9P9uhMlNYZaMfSWT0/zvpk3e2Z4/r3dauqSderSv72i/H6/pWUh69glkHuCfzinJtqXgTx5TU1NMWvieoJyKvrlEtW6jEzgYSSojLbZNFv7svVK9J7qqU3r+Qzq6upMy0DXc2OQSZ9nW1efmrXsMXXmLU+pHu+AZeUg69klkN8D4FUAeTGevwhAA4CGyZMnp+SNynbxgriVyTJ6vlyTaa7W6os3sxY82psDu9fKE93IJOrfjvfZar03ZtWYQ+fWUzYr+8h9g0Nq3p3Pq6OuXa/eaeu2pAyUOTI+kAP4DYAWAIfo2Z818uQk0y+eLom+XBN9WVtVO042GGdyYpUeegKr3mCpNW7cSN+23huBSLHK53A4MuLm6voHt6iqxevU/a/usKwMlDkyOpAD+C2AVgBH6j2GgTw5mRzIQ5INyOmakc0sVjbLj7YlIFELSug16G3GDt/fjBn79N4kZfLN1H2v7FBVi9ep6x/cYnVRKENkbCAHcCuAnQCOMnIcA3ly4gXyhoYG9dZbb1ldxLhTl8ZjVmBMV3O3VTce8YLraKdm1fqswt/PoqKimMeVlZVFDU3TcwOglQWf6HchvExut1tX8l06u0Fe/6hLHXHNI+qCO/+jfINDKbsO2UtGBnIAqxAYQ34mgIqwx9hExzKQ65coeCcK5umcdELPsCwjxxqtXaWzhmZVjdyMxMJ4fd9ut1uz5mu0pp3oISIxP69kjkmU+Z6u34u9+/vVqcs3qBNueFzt7O4z/fxkX0YCuQT2T73gl4GWWqXU9fGOnT17tmpoaDC/UFmmsbExqeNmzZoF4OMpWcN/J0QE5eXl6O7uhs/ng9PpRGVlJcrLy0dd3srKSrS2tkZt93g8w+udx1NfX4/a2lq0tbWhoqICS5cuNTRt5mivb0Ro6tHwlblcLlfcqT5H+/oAwOFwINHfeKLXG+t9Cgl/HVqv0wwejwcANMvhcDg016ePd0y815yu34shv8J3/vwyXtjajrsvOgmzqsaZdm6yPxFpVErN1rNv2uZaV0pJjMf16SoDaQvNuR5rWdL29nb4fD4AgM/nQ3Nzc8J52vWINZd4+JdovPnSR7tcZDrnODc6X7fWnOPxlvyMRc9SrlqvN/x937t3b9zjw5cJ1VoWdrRCc9/H+ly0gniiY+J9xun6vfj1Y2/j2Xd24/ovHWM4iBtdN4GyGxdNoeHAHArWiSilTKmZxAsy9fX1pgUzo9dPFPwiby7mz5+vay1xIzceZq2TrrUoTKTI1xv5vusJzKEgZ0awKyoqgtvtjrrhKS0t1X2OefPmoaamRvdnXF9fj3HjxkFEYrZg6Lkp0uuhTS244+mt+Pqcg/HNEyYbOjbUcpaKm2uyJwZyglIK27dvN3SM3qAfT7wVxpYsWWJaMIt3faMrn2ndXKxZs8b0mw2zaoXhLQFatF5vMrXqUJArKyszdJyWvr4+dHZ2RnUnRK60F8+GDRsA6PuM6+vrcdlll8VdPc3MFfG2fNSFK9Zuwuyqcfj5udMNvS4gdsuZ2d1BZB8M5JQUp9M56nPEq5F2dnbG7JeN11+rZ+nS8OsbXZ5ST5Az42Yj2dYCLaGWAKUU6urqhl9vWVkZXC4XLr744hHvVTK16rlz56K+vh49PT1RzxUUFKCurs7wOSNviuIFWq1jQzkGXq8XDkfgq07rM66trcXAwEDMc5m5bGl7Tz9qVjdiXLETd1bPgjPf+FdwrJtoM26uyZ4YyG3utddeQ2NjY9KJbskKrYNcX18/3AwqIhg3bpyh2mismmI8oS/lSMk0xYc3dy9duhS1tbVxbwL0Bjm9+8W68YjVJD537lxd540l9HrvvPNO9Pf3o7Ozc8R7NX/+fMM1RCBQA66trcXg4GDUc8XFxaipqUnqsx7NTdFll102fNPn9/uHa9WRa9XHuzEEkFT+hRbfoB+XrHkF7T39uGv+bEwsKYzaR0/fd6ybaDNursme0pa1PhrMWtf22muvYWhoKGp7Xl4eZs6cOargHq+vEAhkus+fPx9r1qyJeq6goAArV67U9eUX6xyJaJVNb7axVjY4AF1Z5YkyuGNdU+u6U6ZMwQsvvDBin/Brar03iTLd9dL7OvQKBX+tz0VE4Pf7k85oDx2fzA1GJI/Hg6VLl+ouh8PhwK5du9DS0jKqURtKKfzsvib8/aUPcevXZ+LcmQdG7RNr1EhVVdWI6+ndj+zNSNY6A7mNxQvUs2bNSmktfePGjVi0aFHM5/UO1UkmoIS+jCOD8cUXXxw3kACxh4G5XC7Nplutm4BEQUAr2BoJYqFrpnIYlJ5haeGKi4vR29sb83m9Q71CNzNGPvPQ8WYEchFBRUWFoes3NjaOOmj+/rlt+MXDb+KST0/DlWcdqblPU1OTZvO40+nEjBkzRmzr6OgY9c0FZTYjgTw/1YUh6xQWFqK/vz8l577pppviPt/W1qbry8Zof6zL5cLcuXNHBMVQs3BpaSm6urqijgnvV46VQBcrwEaWLxScw28i5s6diw0bNsQc711fX49LLrlEc5hUvGvqGZ6XLCPBLNSfHiuQhyeCad0khSeJ1dTUGEpeCz/e7XYb6ifXUlFRYeh3rqKiImZimd7A+cQbO3HDI2/i7OkV+Onnjoi5n5G+7/LycgZuGsY+8iw2ffp0FBZG98OZYefOnXGfLy0t1TVExkjyVijpaMOGDZrBWEQ0+5VbW1uH+5+N3jholS9yGNnq1atj9rPPnz8fixYt0h3Ew6+ZaHhe5P+NDInTMywtvDzx3rdQ64PR5MFYfeYOh0Pz+JtuugkFBQW6ypyfn6+5b29vr+7MepfLhYULF2o+pzex7M3Wbvzo7lcxvbIMv/nqTDgcsW9e2PdNyWIgt7G8vLyYz4US4AYHB4dnbjPTpEmT4j4fmjowcltkk7CRgBJKOooVVLq6umIOtQrV2mN9iYcyuMPFG3KklaQWa2ia0T7h0DXjDXdavHjxiLIYHRKnFXSrq6tjvgd6b7iMjJWPNTTsjjvu0Dy+pqYGK1eu1JU0d/vtt2PlypVwu90jtnd1daGnpycqyLtcLlRXV0fdhJx77rma59cTXHft8+L7f2lAiasAv//2bBQ5Y/+9AoFupshWChEZTiwlioV95DYXK+EtXF5eXsJ9jBARPPfcc1iyZEncIFVRUYGdO3di0qRJWLhwIc4+++yo8uTl5eGZZ57B8uXL0dbWNqI/O1JdXR1qamp09R3H2ideEl9ZWRlEBF1dXXGnRDXaz25EdXU1Vq9ePaK8sYT6e3t7ezW7FCKVlZUlLF+saWHr6+uxYMECzWP09NnHO28y09DG6+N3u93Yu3dv3ETKsrIyFBcXJ7xusollvb5BfP2u/+LdnT3454KTMP1Afa0A7PumECa75ah0DEEL/xK75ZZbsHjxYl3NxgUFBbjuuuuGg3nkOZ977jksX748br9tKGDombfcaDJXrPNoMSvj+6STTsL27dvjBhMzErzCRd4ohIuVzR/aFuv9jHfzFTqv0XnmE4n3GZSVleGLX/xi3NEQicocLjy4hlrBhoaGYgbawSE/alY34qm3d+Gu+bPxmaMnaZ6LgZriYSDPUlqBOrzZPF1jySdMmICqqio0NjZi/fr1uOGGG3Q1H5eVlQ3PuBXOyDlCtfLwoFNaWhpVkzaaGR0uUQ0z2ZsEo9cBgHHjxo26lh+prKwM3d3dUTXjyGBbUFAApZTm2PBwiV6H2+3WbDGINfpAT3CP10IQmtMgXqD2eDxoamoyFFT11M6VUrj2gS1Y898PsOy86Zj/ySpDxxOFZOSiKTQ6sYJ0uieCAYDVq1fD4/Fgzpw5uPnmm0dM0BJvPuxYzb+rVq3S3Y8c6utNNLHJ3Llzdfe9R0qUEGfWnNt6Eu+MJHjp1dXVFdV3rpXNPzAwkDCIJ5q6dP78+TE/99D1k5netqamJma+Q0VFRdwg7nK5sHjxYsPzleuZGrXumW1Y898PsOD0aSOCuN7jiZLBQJ5F4iW/mSVUew41tXZ3d48YlpTMNJGJMuDDRc70FWs42YYNG0Ykc8WaDU6LVqAOT27r7e01FFxjNY+XlZUlnE5Wb4KX2+2OSuzSI/R+Gs3m17uCW7zmbYfDofnZLViwAPPnz09YhuXLlyM/P3oEbaKWmHnz5uHMM880HFQTDQ974LWPsPzRt/DF4ypx5f9EDzPj1KqUKgzkWWTmzJkpD+aJas/xnotVW0+UAR8pPOgkGmsdyqC+4447omroBQUFUYFAq4YZmRUeqtEWFRXpKq9WM3xeXh56enpG1Eb1BjAt55xzju7yRAo1a+vl8XhiZqWH3/Bccsklcc8Tr9a8Zs0aXe9FMjkEa9euxQMPPKD5nM/ni1krj/W35XQ68fTbu3D5PZtw4tTxuOWCY7F3756oqVY5vIxShYE8y8ycOROzZs3SPeTMaODXW3vWqil1d3fjnHPOwfr164e3rV+/XnOiEZfLFbOGGV6TjfdFHt5MqzXcauXKlbj99ttHbJs3b17UfOtatf7BwUH09fXpei+0DA0NaTZbr1mzZkTNPPwmIp4NGzYknRMQ6pvW2xURPt97eOB2u93D85srpQyNndfy17/+dcT/I6918cUXx13sJBav14tVq1bFfF6rib2jo0Nz5IeIoAMluHjNKzh8Ugl+9+3Z6Onq1Gy2D+VyRB7P4WU0Wkx2s4lE07EaPSZZ55xzjq4+5IULF2LVqlVx9y0uLobP54sKaG63e3jmuGQTsEKMTGcaK7va6Djw0QoNnwLMnxM9Un5+Pm6//fYRCYSJrhdv9IDZPB4P2traUFZWhp6eHt2feyIigpdffjnhfqEkuFBSXKSWHj+ufWYvyooKsHbByZhYUhh3qtXwczFrneJh1nqW0grMiRY3MVuiDHOXy4Wrr74aABIG8nhCGc3AyOlQ9Y6ZDtHKXo41djmZsedaHA7HqGujoex8vRny8a5ZUFAQs+bqdDqjpvHVc/OgJzM8k1VUVGDdunW69o31+e/aP4Srn+pAXn4B/rXgZEwuLwaQ3E03USRmrWepUJN56JHuIB4S3qdXVFQ0om+2sLAQmzZtGk6IS1YogxnAiJnCuru7DZ0nsu833lKnscprNIib8ZmEEvr09F27XK64AXXlypUxn9OqOer53MxoOrdKvKlXtWh9nnv7hvDzZ/fAO6jwf989YTiIE1mBgdzG0h3EQ7Xx8GDa398/oq+4q6sLa9euNaW51ev14sorrxyxLnOswOZ2uzX7ePv6+kb0OcfKcl+8eHHcoXORtHIAgEAClxmfS1tbG+rr6zX74QsKCoZnoQtljsfLaje6nrdZw+usENkHHZp6NjwP4pprrtGcmEivLu8Qrn9mD/b2+XHNaeNxlEf/7w1RKjCQk25aGeuprpV1d3fjgQceGE5AijU/90033YQVK1ZEZW53dnbi0ksvxbhx4+BwOGI2GXd1dRkaBmRWX20spaWlWLRoUdRkMGVlZVi5ciU6OztHZI7HS1Yz2sduJPEtUqIExFRyu9248847h+f5Dw0nC1/UpqWlBQsWLEg6U3yfz4/aZ/diV+8QfnbaOBxRXhCVGMfsdEo3BvIME1rsJPyRKYyM9zbTqlWrhr+UE62wpVWDHRwcHJ4wJp7RZKGbraurS7NVo7g40IQbOf48/H0xYtq0aSP+HytLXw+3221a10IiIhI1lr+zs1PXhDLl5eVR63vrsX/Aj2XP7kHLvkEsOWUcjpkYCMyRyZRc/ITSjcluGcRokswrr7yS1uZ1PRnrqRCeYRz+PkQmrfX19Zk+namZzEiCA7Qz6UP5Eslew+PxYO7cuaZ1i1gp3kiF8LnOjdg/4McNz+3Fe3sGcOXJbsyuHNliEfn3yTnVabSMJLtpd/SRLRx//PFpDeannnoq1q5dm5ZrhQtNGBPeNBk59MmMIVputxterzdlgcyMIK41Gxrwcb5EstcILX2aDWLdbGrNda5Hj8+PZc/uxfudA7j8pOggrtVkXl5ezsBNacOmdZs7/vjjDU0AMxobN27Utd+8efNMu2ZeXh4WLlwY1TSZbPNvPF6v19BUrumWKDudAiZNmjQiQTJEa67zRLr7/bj+mT3Y3jWAZWdNxicPGpmDwSZzygSZ+61Fllq/fj3OOecczJkzZ3g2tnh95KG1sZctW4YlS5YYznzW6vMMbXc4HFErRKWiid/r9WrOMpcJ3G53Un3guSY0tExrERSjzeld3iEsfWYPPto3hN9/ew6+efoMVFVVDdfAnU4nV1CGcfgAABnbSURBVC6jjMA+8gySzEQSzc3NaG9vN7UcWpO+uFwuFBYWak7GojW5xvr161FbW2tKdrdWn2eqZzzLRNXV1Tj11FNTPpuaXbhcLpxzzjnYuHEjdu7ciUmTJmHhwoVRQ8ucTidKS0sN/Z209wbGie/u9eOPF56AUw+bYHbxieJiH7lNzZo1K+6a4+nKYNcaZub1euF0OqMSreJNrmHWMCSt2vfcuXOzpk9XrzVr1mDNmjXDc9DncjAPTQOsZzy4z+czFMR3dA/i58/uQd+AwnWnl6OouxlNTa1MWKOMxUCeYcyaNz3WTYEesZrQu7u7UVpaisLCQnR3d2vWgNavXz+qqVm1KKVw6KGHYuvWraad0846OzvhcrlQXV2Nhx56yNCUtdlARHRPr2rUOx0+3LBxL/IdgmVnlGOKO/AVGWqqB8BgThmHgTxLjab2PmnSpJiBuLu7Gy6XCz//+c+jakOJ5mEfDQbxkbxeb861SIQYXfZWr1fb+nHzfzoxYawTS08fj/JC7fXKGcgp0zDZjaIsXLgw7sxesZaBTLRWeaZJdv1uSr3S0lLNEQQFBQWG5knX64ltvbhx414cPK4Q9116alQQDzGaMEeUDgzkFOXss8/G1VdfHTfzXKv53aqZ35LldDo1M+XJOgUFBVi2bBmefPJJ1NbWjpj/vqysDNddd92o5kmPNKQU/rKpG3c2dmPGAU784owJOKAkuelpiazCpnXSdPbZZ+Pss8+OOZtbqHnzpptuwr333mvL8c2Rfcv5+fkpn0OdYisrK8NPf/rT4UAd+h00w9ixY9HX14ehoaHhbd5BP1a82IWXW/px1rRifHdmCfLAz5/sh4HcJkaTvDYaCxcujOr3zs/PR19fH2bP1jUywjYYxK1jZH3wZEQG8V37h/Cr/+xFc+cgvjezBJ8/bMzwc83NzXA6nZrN6Fz4hDIRA7mNRGa0pyOwh2pEq1atws6dO1FSUoK+vr6cy5Sm1DG6PngywoP4a239+O2LnfD7gSWnjsMsT+GIfdvb2zFhwgR0dHSMmAkufBY3zqVOmYQTwqRRvDHiyR6fblYtnELZych48NHyK4V739qPu7f04ODSfFxxshuVJdp1mVBw1grWWnO2iwhneSNTcUKYDBQrCDc2NuoK5pkQxAH7JbRRZnA4HCgpKYk5/0Cq7fP5cftLXWho7cdpk11YMKsUrvzYub4+ny/mwidac7ZzaBpZiYGcDIk3xpxISzpr3VqadvXjtpe60O3143ufKMHZ04oTzjoYry881hA0Dk0jqzCQkyFayW8hDocDTqfTVmPJKTUiM9CtMOBXuHtLDx54ez88Y/OwZG45po1LPNww0YpmTISjTMNAToZEJr9FNpOuX78e1157rZVFJItYXfMOt6N7ELe91ImtewfxmalF+M7MkhFN6Xl5eTj44IOH+7yNJK5VVlZq9pFzOVOyCpPd0iSZlc30Hp9p5s6dy6z2LDRv3jwsWbJkeD79eCuOWWXAr3D/W/ux9s0eFOULLp5dhhMPHDnBS2jJXJ/Ph7y8PACBrHYj2efMWqdUM5LsxkCeRnqy1jNh9bPRYq3cngoKCjA4ODiipplJtexE3tszgFUNXfigaxCnHOzC92aWoMyVF7WfiEQlq4U/x+xzygQM5DZltNaeyYH9zDPPRHd3t9XFoDicTicGBgYyrlZt1D6fH/94vQf/fq8XbpcDF80qxZzK6GlWnU4nhoaGRowp1+J0OjFjxoxUFZdIFw4/I8tdccUVmklxDofDltO5ppqIYOrUqXj//fej+l6VUigtLYWIoKura/g9DNWWgdg5C9lsSCk8vq0Pd2/Zh/0+hc9OK8L/zijBmIJAX3iyN7/MPie7YSC3MaumbdUjVlJc5LZTTz0VGzduRFtb23CA0gpUN954I3p7e5MqS3FxMXp7e6POr+WEE07AHXfcgfXr1+Pmm2+O2aoQCrChMj700EN46aWXos7V1NSEvr6+4WPOP/98LFmyJKnXEU8uBO5wm3f24y+b9mF71yCOmViA784sxRR34oz0WBnnkfsQ2Qmb1jNIoqb1RH3smRrUiczydocPf2vqwZbdPkwoduBbx5bg5INcUePCYzWPa83KFo595JQp2LSehUY7MxyRnW3dO4B/vN6DxtZ+lBU68N2ZJfjsIcVw5kVP7BJvKFgoQIcyzpPNWifKJKyRZ5hYtW7WtinXKKXw2k4fHnh7P5p2+TCmQHDeEWPw+cOKY06vymBM2YI1chszq3bN4E921T+k8J8PvXjonf1o7hrEeFegCf0zhxQNJ7KFY/CmXMdAnqUYxMluPto3iMe29uLp7X3oGVA4qCQPC2eX4rSqIhQ4ovvAGbyJAhjIicgyPT4/XtjhxbPNfXijfQB5Apx4oAufm1aE6ROdcDgcXC6UKAEGcptgUzlli74BP15p68fGD7x4pa0fg36gcmwevjl9LM6cWoRxwdnY4q0JTkQfYyC3ET3TuRJloo6+ITS09OPlFi+advkw6AfcLgfOmlaM0yYX4dDxI8eAhzLPY60JTkQfYyC3MdbSKVP1DfjxRrsPm3YGHju6BwEAFWPycPahxZhT6cKREwqQJzK8Ehlr3kTJYSC3OU4IQ5mgo28Ib7X78Fb7AN7u8OH9zkH4FeB0AEdNdOKMKSWYe9QkjB3aFzV5S2g5UQZuouQwkGcR1tAp1fxKYdf+IXzQNYhtnQPYtncQW/cOoNMbmPLWmQccNr4A5x0xBjMOcOLICU6MLSocrmFz+U8i8zGQZxn2o5MZ+gb9aOsZQuu+IbT2DKJl3yA+7B7Eju4h9A8FssgdAA4szcfMSU5MHVeAI8oLMNVdgHyHxAzSrHkTmY+BPMux6Z1CQkO3+gf9eGPbDuzs9qKzH+joHUR77yB27R9Ce58f7b1+dHpHLvU53uXAZLcT58+cgGOrJuLwihIcVVGK3n2drGETWSztgVxELgFwBQAPgNcBLFJKPZfuciQj0aIlZhwfb59kA3Ho+NGeh1IrFAgBxF3YAwhMXzrgFxSWlGH3nm509vrgU3lwji2Dyndhb+8A2vZ0oaWjG119g9g3APQMCDr7dmNf/2DU+QrzHThwXBEOLC/CnEOLcNC4YkwpH4MpEwI/xxRqf1UUsYZNZLm0zrUuIl8DsAbAJQA2Bn9+B8DRSqkPYh2XCXOtJ1qZzIzjGWCtpZSCH4BfASr48EPBrxD2UCP+PeT/+N+OvAKMnzABQ36FHS1tGPT7MeQPrJs9pATu8RPgKi7GwKBC/5AfvkE/BoI/+weH0D/gR3/w3137+9C5rxfegSEM+AVDkgfvgB/7vAPwDvrhHQyUIx5nvgPjigswrtiJ8WOcKB9biAljnZgQ/Dmp1IWKMhcmlbjgLi6ISkIjIusYmWs93YH8RQCblVI/CNv2LoC1SqmrYh1nZiB/dEsbfvv4O4aPC60praWoqGjUxyso9PV5DZcrloSfqo7PPXKPxOccuV/U/kr7vypig4IK/DPifMP/VoGfUf8GABU4VkXs4x/+qSL+H/bTQnkOQWG+A4X5DjjzHSgqyINr+OHAGGc+ipx5wz/HFuZjrCsfYwrzMbYwD6WuApQWFaDUVYASVz7cxQUoKshjcCayqYxcNEVEnABmAbgl4qnHAJyssf9FAC4CgMmTJ5tWjhJXPqZOGGP4uL2dvpjPjXMnPp+e4zs7BwyXK57E3+GJv+Qj99AbF0K7Re4f6/DQfjJim2juLzLy/JHXEsiI7Q4JnUshLy8PrsJC9Pd7h29mAs8DDhG43W6UjCmGwyFwiMAR3J7nCDxEJDD22QHkORzDP/OD++c7BHl5ggKHA/l5goI8QX7w3848B/LzHCjIEzjzHXDmOYZ/5udpr+ZFRJRIOvvIJwDIA7AzYvtOAJ+J3FkpdReAu4BAjdysQpxy6ASccugEw8exaT0zhPqRP/zwQwwNBRKyQhOKhPpqm5ub0d7ePnzMhAkTUFVVZUl5iYhSzYqs9cigLBrbiKLonbazqqqKgZuIckY62/PaAQwBqIjYfgCia+kZJ1atW2/Wup7j4+1j1jrldiAimDJlCqZMmQKn0wkgUBPnqldERNGsSHbbpJS6KGzbOwD+la5kt1xjt+Z6jkUmIsrQZLeg3wBYLSIvAXgewAIAlQDq0lyOnJFLNXkiolyU1kCulPqHiJQDuAaBCWG2APi8Uqo5neUgIiLKFmlPdlNK3QHgjnRfl4iIKBtx8CoREZGNMZATERHZGAM5ERGRjTGQExER2RgDORERkY0xkBMREdkYAzkREZGNMZATERHZWFrnWk+WiOwGYHT2twkILNSSzfgas0MuvEYgN14nX2N2yITXWKWUmqhnR1sE8mSISIPeCeftiq8xO+TCawRy43XyNWYHu71GNq0TERHZGAM5ERGRjWVzIL/L6gKkAV9jdsiF1wjkxuvka8wOtnqNWdtHTkRElAuyuUZORESU9RjIiYiIbMyWgVxELhKRp0SkU0SUiEzR2GeciKwWka7gY7WIuBOcV0TkehFpEZE+EXlaRI5J1evQS0SmBF+n1uOKOMd9OsYxR6az/EYE3/PI8t6t47jzReQNEekP/vxyOsprlIiMF5GVIvJW8HfsQxG5U0TKExx3YYzP0pWusscjIpeIyPsi4hWRRhE5LcH+pwf384rINhFZkK6yJkNErhKRl0WkW0R2i8hDIjI9wTGx/m7PSle5jQh+90WWtS3BMTNE5Jng7/JHInKdiEi6ymyUiGyP8Zk8HOcYrf0z6vfVloEcQDGAxwBcH2efvwE4HsDZAM4K/nt1gvNeCeByAJcBmANgF4DHRaRklOUdrQ8BeCIelwBQANbqOP6YiGPfTU0xTfMnjCxvTbydReQkAP8A8FcAM4M//ykiJ6a4nMmoBHAgAr9rMwBUA/jU/7d37sFXVVUc/3wTHyiQkSJO5XMclcB8gZJPTIoos8ARC0fNNDUztUzHLHxOZJlCppk5isoI+JrxiS8MNUUEzWcCMghqIA9L8afydPXH2jcPl3vv73fv73c59+L6zJz5/c4+a++z1tnn7rX3PuucDYxrQ94PKboPzGxZnfRsM5KGAaOB3wJ7AE8BEyVtU0Z+e+D+JLcHMBK4UtLQdaNxTRwMXA18FTgEWAU8Iql7G/IOYs16e7ROOnYEM1lT1z7lBCV1Ax4GFuLt5c+AXwI/r7+aNdOXNe3bE29Hb20l34lF+W6so47VY2ZNuwF745WwXVH6ril9v0za/ilt5zJlCVgAnJdJ6wy8D5yUt60l9H0YeKgVmYOTzVvkrW8Vdk0G/lxlngnAw0VpjwDj8ranjfoPBj4GulWQOQ5oyVvXMrpNBf5WlPYaMLKM/KXAa0Vp1wFT8ralCpu7AKuBwyrIbJd+f3vnrW8bbboAeLkK+VOApUDnTNqvgX+TAqkbfQPOA94FNq0gY8AReetaaWvWEXlr9Ada8B5/gSeBD/AedSm2B3riI30AzOwj4PEKeXIhjWi+RttfkZguaYGkSZIG1FG1juIoSUskvSLpsjbMiPQnU2+JB2mweqtAN2A5PuKuRGdJ8yS9JeleSXusA90qImkjYC/Wvv4PUf76l6uvvSVt2LEa1o2u+Izmf9sge6ekRZKelHREnfVqLzukKfLXJY2XtEMF2f7AE6mdLPAgPuu0XT2V7AjSI4AfAWPNrLXf3ujUJk2TdLKkhvKdDaVMB9ITWGypOwWQ/l+UjpXLAz5NlGVhhTx5cSL+HeC7WpFbgPeahwJD8GmzSZIOrK967eIWYDgwALgY1/3OVvL0pDnqbS3kcRsX4yPaVRVEZwLHA4cD3weWAU9K2qn+WlZkC2ADqrv+5eqrUyqvGRgNPA9MqSDTApwFHInPukwCJkg6uv7q1cRUfObnm3gb0xN4qkL8Rrl6LBxrdAbiA7jrWpEbAQwDDgXGA38EflVf1aqjU94KFJB0CT7NUYkBZja5jUWWekFeZdIr5WtLnpqoxWZJnfAf2xgzW1kpo5nNxB1AgSnywMCz8JmGdUI1dppZdpbhJUlzgKmS9jSz5yrkX2f1Vooa63Iz4B58KvLsShnNbAoZpyHpKdyRnIY/m8ybaq9/KflS6Q2HpMvxR3X7m9nqcnJmtgRv9AtMl7QFXtdj66tl9ZjZxOy+pKeBOcCxwOXlshXtN0094p2VaWb2fCUhM7s4s/u8pA3w3/ol9VSuGhrGkQOjaP3mfqONZb0N9JCkwqg8TaNsydo9yGwe8J7km5n0HhXytJdabD4MD7ZorRdZjqnAUTXmrZX21O10/FnkTkA5R/42a48A6llvpajKRkld8IAvgG9blUFrZrZa0nT8uuTJErx+qrn+5eprFfBOh2rXwUi6Av/9DDCzOTUUMRX4YcdqVR/MrEXSK5S/x8rVI6zb317VSOqBz26dWkP2qUA3SVuZWUPY2TCOPPVeO2rZuCl4MEp/PnlO3h/YjDWfm2d5Hb8xBwLTANKrPQfgkZgdTo02nwg8Zmazajzt7viU+zqjnXXbB5+6raTzFLze/pBJG0j5uu5wqrExPfOfiI9eBplZS7XnSx3T3YAXqs3bkZjZCknP4tf7tsyhgcAdZbJNAb5blDYQmN7aLFOeSBqNO/GDzWxGjcWs899fraT2bxfg72VEpgCXStok0xEdCMwH5tZfw3ZxHB6X0uqrrSXYHX+09W5HKtQu8o62q2XDe4G7Az/Ap3AGp/3uGZmJwEvAvrgTfwm4J3P8C8AM4HuZtHPwKMwhQG+8kucDXfO2Oem3DT76GV7m+E3ATZn9M/AGcyf8FbSR6XoNyduWMvrviD+P2hsPlhkMvIqPxDfIyE0iExGNB1WtAs7FG55zgZXAPnnbVMLGrngDWBjp9MxsG1Ww8XzgG8AO6V6/PtnYrwFsGgasAE7A3xgZjT8f3rbMfbk9Hng6KsmfkPIPzduWCjZeldqGQ4rqrEtGZiQwKbN/bGqjdgV2xh9prQDOzNueMjZeBhyU6mcf4N5kc6Eei+37LD74GZ/ayyFJ/hd529KKnQJmUfSmRTr2U2BGZv8wfPDUO7VPJwDvAaPztmMNvfNWoMaKuCA5pOLtuIxMd3yqc2naxgKbZ45vVyKPUtkL8B7XY0DvvO3N6Hch8B9gkzLHJwOTM/tnA7OBj1K+J4DBedtRwb4vpWv+Dt5bnp2cQvciubl4jEA27Qi8Y7YCd/6N2lk5uMy9a/hIr6SNwBXAvHRdFuHRwf3ztiej30+SzsuBZ4EDy92XKe0gvIO2HJ8NOzlvG1qxr1ydXZCRGQPMzewfC/wL77QsxR8THZ23LRVsLAxcVuBxG3cAvcrZl9L64PE2y1K7eT4N/uoZHkhrlOgEp/bfMvuDgH/iryF/gA8ITwc65W1HdotFU4IgCIKgiVlfXz8LgiAIgk8F4ciDIAiCoIkJRx4EQRAETUw48iAIgiBoYsKRB0EQBEETE448CIIgCJqYcORB0IBIGiPp3jqVO6Kjy60HknpIWizpi3nrEgSNTDjyIPiUIKkP/n3pUXnr0hbMbBH+VbgL89YlCBqZcORB8OnhNOAOM1uatyJVcAMwXFL3vBUJgkYlHHkQNDiSNpY0StJCScskPS1p/yKZb0mamY4/LukoSZaWrSUtvXgkcHdRvrmSRqQp9/clvSlpmKTNJY2X1CLpNUlfL8rXS9J9Kc8iSeMk9cwc7yvpIUlLJC2V9A9J/YvKMEk/lnSbpA8kzSleq9vMXsY/Gzqk/VcyCNZPwpEHQePze3xhkuOBPfDvPT8gaWsASdsAdwL3AV8B/pTyZNkNX+RieonyzwCeAfYEbgVuBG7Bl1ndHf+W9ti0GhbpvI8DLwP9gEPx1QbvllRoU7oCN+OrB/bD106/P63HnWUEcFfSewJwvaRti2Sewb/NHgRBCcKRB0EDI2kz4BTgHDO7z8xeBU7G13surKV8CjAHX3VqppndDlxTVNS2+EIRb5c4zYNmdrWZvYYverExMNvMbjKz2cDFwJb4ClCF871gZueY2atm9iJwDNAXX7kOM3vUzG5Ox2fg0/rL8EUostxsZmPTeX6Dr2J3QJHMfHyRoyAIShCOPAgamx2BDYEnCwlmthpfCrVXStoFmGZrroA0taiczsBKM/u4xDlezJTdAnyIj/oLLEx/e6S/ewEHpmn3FkktwJsZfQsR53+VNEvSe/jqUT3wpXjLnXsVsDhzngIfJf2DIChBp7wVCIKgIkp/Sy1TaBmZ1pYxXAJsJGlTM/uw6NjKEuWuLNqHTzr+n8Gn8c8qcZ6C078R2Ao4k0+WN50EbNSGcxcPMLrjDj4IghLEiDwIGpvZ+PrQ/w9uS4Fr/fG1rsHXX+9blK9f0f7z6W8v2s9zwJeBeWY2u2h7P8nsD1yZHge8go/It67xfL3TOYMgKEE48iBoYMzsA+AvwO8kDZa0a9rfCrg6iV0D7CjpMkk7SxoCnFQoIpWzGHeGa0S718hVeODcBEn7SNpB0qGSrpXUNcnMAo5O0e19gfF4h6QqJG2KT+U/0AF6B8F6STjyIGh8zsGjyW/AR9a7AYPMbAGAmc0DhgLfAV7Ap7MLH1FZlinnWmB4e5Uxs/nAfsDHuIN9BXfuy9MGHmHfBXgWd+LX41Ps1XI48IaZPdE+rYNg/UVrxscEQbA+IOl04CLgc4UAN0kbAzOAY5rFMUp6BhhlZrfkrUsQNCoR7BYE6wGSTgWm4UFh++Kvco3JRqmb2XJJx+LBYw2PpB7A7cC4vHUJgkYmRuRBsB4g6Qr8y22fB97Cp7MvMrOqn0sHQdBchCMPgiAIgiYmgt2CIAiCoIkJRx4EQRAETUw48iAIgiBoYsKRB0EQBEETE448CIIgCJqYcORBEARB0MT8DwySmVLCqywQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "merge_setup = {'/FullDataset_v1_rna':sample_list}\n",
    "\n",
    "for s, smerge in merge_setup.items():\n",
    "    print smerge\n",
    "    t0 = time.time()\n",
    "    print t0\n",
    "    print '________________', s\n",
    "    \n",
    "    cell_ix = np.in1d(samp_id_flat, smerge)\n",
    "    \n",
    "    run_all_spring_1_6(E[cell_ix,:],\n",
    "                     list(merged_list), s, main_spring_dir, normalize=False, tot_counts_final = total_counts[cell_ix],\n",
    "                     min_counts = 3, min_cells = 3, min_vscore_pctl = 90, exclude_gene_names = list(merged_proteins),\n",
    "                     show_vscore_plot = True, num_pc = 60, pca_method = '', k_neigh=4, use_approxnn = False,\n",
    "                     output_spring = True, num_force_iter = 2000,\n",
    "                     cell_groupings = {'Sample': list(samp_id_flat[cell_ix])})\n",
    "    \n",
    "    print time.time() - t0\n",
    "    \n",
    "    np.save(main_spring_dir + s + '/cell_filter.npy', np.nonzero(cell_ix)[0])\n",
    "    np.savetxt(main_spring_dir + s + '/cell_filter.txt', np.nonzero(cell_ix)[0], fmt='%i')\n",
    "    \n",
    "    print time.time() - t0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Sample1']\n",
      "1612985097.21\n",
      "________________ /FullDataset_v1_protein\n",
      "Normalizing\n",
      "Using user-supplied gene filter\n",
      "Running PCA\n",
      "Building kNN graph\n",
      "kNN built in 0.83048 sec\n",
      "Saving SPRING files to /home/I0369218/precision-immunology/CHAMBERLAIN.Mat/CITESEQ_EXPLORATORY_CITESEQ_5K_PBMCS/FullDataset_v1_protein\n",
      "Running ForceAtlas2\n",
      "172.316283941\n",
      "172.402792931\n"
     ]
    }
   ],
   "source": [
    "merge_setup = {'/FullDataset_v1_protein':sample_list}\n",
    "\n",
    "for s, smerge in merge_setup.items():\n",
    "    print smerge\n",
    "    t0 = time.time()\n",
    "    print t0\n",
    "    print '________________', s\n",
    "    \n",
    "    cell_ix = np.in1d(samp_id_flat, smerge)\n",
    "    \n",
    "    run_all_spring_1_6(E[cell_ix,:],\n",
    "                     list(merged_list), s, main_spring_dir, normalize=True, tot_counts_final = total_counts[cell_ix],\n",
    "                     min_counts = 3, min_cells = 3, min_vscore_pctl = 90,\n",
    "                     show_vscore_plot = True, num_pc = 17, pca_method = '', k_neigh=4, use_approxnn = False,\n",
    "                     output_spring = True, num_force_iter = 2000, gene_filter = adt_ix,\n",
    "                     cell_groupings = {'Sample': list(samp_id_flat[cell_ix])})\n",
    "    \n",
    "    print time.time() - t0\n",
    "    \n",
    "    np.save(main_spring_dir + s + '/cell_filter.npy', np.nonzero(cell_ix)[0])\n",
    "    np.savetxt(main_spring_dir + s + '/cell_filter.txt', np.nonzero(cell_ix)[0], fmt='%i')\n",
    "    \n",
    "    print time.time() - t0\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}