data_prep/spring_notebook_10X.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#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 collections import defaultdict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "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-03 15:48:52--  http://cf.10xgenomics.com/samples/cell-exp/2.1.0/pbmc8k/pbmc8k_raw_gene_bc_matrices_h5.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/2.1.0/pbmc8k/pbmc8k_raw_gene_bc_matrices_h5.h5 [following]\n",
      "--2021-02-03 15:48:53--  https://cf.10xgenomics.com/samples/cell-exp/2.1.0/pbmc8k/pbmc8k_raw_gene_bc_matrices_h5.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: 32549555 (31M) [binary/octet-stream]\n",
      "Saving to: 'pbmc8k_raw_gene_bc_matrices_h5.h5.5'\n",
      "\n",
      "pbmc8k_raw_gene_bc_ 100%[===================>]  31.04M  88.8MB/s    in 0.3s    \n",
      "\n",
      "2021-02-03 15:48:54 (88.8 MB/s) - 'pbmc8k_raw_gene_bc_matrices_h5.h5.5' saved [32549555/32549555]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# download raw data from 10X Genomics\n",
    "!wget http://cf.10xgenomics.com/samples/cell-exp/2.1.0/pbmc8k/pbmc8k_raw_gene_bc_matrices_h5.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_list = [\"Sample1\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pbmc8k_raw_gene_bc_matrices_h5.h5\n",
      "(737280, 33694)\n",
      "737280\n"
     ]
    }
   ],
   "source": [
    "# D stores all the data; one entry per library\n",
    "D = {}\n",
    "\n",
    "for j, s in enumerate(sample_list):\n",
    "    filename = 'pbmc8k_raw_gene_bc_matrices_h5.h5'\n",
    "    raw_matrix_h5 = filename\n",
    "    print raw_matrix_h5\n",
    "    if filename == 'raw_feature_bc_matrix.h5':\n",
    "        genome = \"matrix\"\n",
    "    else:\n",
    "        genome = 'GRCh38'\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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_list = D[s]['meta']['gene_list']\n",
    "gene_id = D[s]['meta']['gene_id']"
   ]
  },
  {
   "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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3X2clXWd//HXm3GUUVQ0MYTWNDM0pDSpJMXMFq2sfqZrbZlre6O1ZWUapukqblkaabqbu2lba6GWpTarVqLmbd6GjYp3aAVagyCmoMKoCJ/fH9/rwMXFmZkzw7mZc+b9fDzOY865ru91Xd9rBq7P+d4rIjAzM6u1EY3OgJmZDQ8OOGZmVhcOOGZmVhcOOGZmVhcOOGZmVhcOOGZmVhcOOGYtSNIMSR7zYEOKA45ZHyRNkvQzSQskvSRpoaRbJM1odN5qRdIoSadL+pWkJZJC0omNzpc1Pwccs15I2huYA7wDuAg4BvgesBw4pXE5q7ltgFOBSUBXg/NiLWSjRmfAbAj7Kim4TI6IZ/I7JI1rTJbq4ilgfEQslLQDML+x2bFW4RKOWe92Ah4qBhuAiFhYei/pw5KukdQt6WVJT0g6S9Im+WMkXZRVy42X1CnpBUlPSTou27+zpGslvZid69OF4/fLqrcOz6q8FkpaIel6SRMquSFJ0yTdlF17efb+XYV7ezl/f2bV4oBj1rsFwO6S3tpPun8CXgX+E/gCcDPwZeB/y6QdAfwaeBo4AZgHnC3pX4DrgUey7U8B35O0R5lzfAX4O+DbwNnAO4GbJL2mr0xK+jhwLbCKVCV4Cqn67EZJ7+znHs02XET45ZdfZV7A/qSH8yrgbtID/iBgZCHdpmWOPQVYDbwut+0iIIB/y20bBTyfpf2n3PZxZEEst22/7PjFwOhCPgP4Rm7bjPTfe83nzYBngYuKeSdVmf2ml9/BDtm5T2z038Ov5n+5hGPWi4i4EZgKXA28GTgeuAZYLOkfc+lWAEgaIWlLSdsAtwIC3lbm1P+TO/ZF4CFSUJuV274Q+DOpWq/oxxGxtJDPh4AP9HE704CtgEskbVN6kQLODcBUSe19HG+2wdxpwKwPEXEHcLCkNuAtwAdJgeeHkp6IiBslTQRmkkogHYVTjC58XhkRTxW2LQMWRcTKMtu3KpOtx8tsewx4Tx+38qbs53V9pBkNLOljv9kGccAxq0BErCJ1Ee6SdDvwG+CTku4FbgJWkHq1/RHoAcaTqtCKtQire7nEql62q1x2KkyXV8rHp4DuXtIs6+ccZhvEAcds4O7Jfo4jlSrGAPtFxC2lBJKm1fD6byqzbWfgiT6O+WP2c0lE3FD9LJn1z204Zr2QtL+kcv9HSm0lj7K2ZLKmhJEdc1wNs/YPktZU1UnaH5hI6v3Wm2uBpcApxe7a2TnGVD2XZgUu4Zj17j+AUZJ+QequPILUCeAI4K/AuaQeZn8FfiTpP4GVpC7Lo2qYr8XA7ZJ+QGp3OTbbdnZvB0TEC5KOBn4C3C/pEmAR8DpSKW058P5SeknHZOcuBbb3SCo9L/4zIlz9ZgPmgGPWuy8DhwIHAv8MbAIsBC4BzoiIBQCSDiI97E8HXgSuAP4beKBG+ToLmJDlbyvgduDzUWaAal5E/FxSN6mt6UukHmpPkbp8/08h+ZeB1+c+H5C9AC7G7T02CIrwhLJmzUDSfqQOCh+PiJ82ODtmA+Y2HDMzqwsHHDMzq4u6BZzSglCF16LcfmVpFkrqkXRzNqAuf46tJM2StCx7zcr31snSTMrWK+nJJkA8VVJ/YxTMzKzG6t1pYB5pNHZJfrDbCaQR3J/K0p0KXC9pQkS8kKW5FNie1JsmSA2ds4APAUjagjQB4q3A20kNqxeReuD02oPHrBlExM30P8DTbMiqd8B5NSIWFTdmJZBjgTMj4ops25GkGXU/AVwgaVfgfcA+2XQjZNO335YFpXnA4aSeN0dGRA/wYHbccZLOiX56SGyzzTaxww47VOtezcyGhXvvvfeZiOh3LFe9A84bsm6Zr5C6Yn41Iv4E7AiMJTfPU0T0SLoVeBdwATCF1OX0jtz5bieVXt5FKhVNAW7Lgk3JbOBrpFlv+1xIaocddmDOnDkbcn9mZsOOpL5muVijnp0G7iZVl70fOIoUYO7I1vAYm6VZXDhmcW7fWNK0HGtKKdn7pwtpyp2DXJp1SDpa0hxJc5Ys8byFZma1UrcSTkSsM+2GpLuAPwFHAneVkhUOU2Fbb5MW9pVGvWwv5etC4EKAyZMne1CSmVmNNKxbdG4dkJ1JU2zA+qWQbVlbQlkEbJvvcZa9H1NIU+4csH7Jx8zM6qhhAUfSSGAX0tQa80nBYlph/1TWttncSZqfakruNFNIKxnm00zNji2ZRpqOZEHVb8LMzCpWz3E435b0bkk7ZuunX04KFj/K2mLOBU6UdIik3UjdmV8kdYUmIh4hzXh7gaS9JE0hdSa4JuuhRpZ2BXCRpN0kHQKcCPTbQ83MzGqrnr3UXkeaqXYb0qqCdwF7RUSpd8O3SKslnk+akPBu4IDcGBxI3Z7/g7W92a4CjintjIhl2Tok5wNzgOdI42/OqdE9mZkNCZ1d3cycPY+FS3sYN7qD6QdO4OA9xjc6W+vw5J05kydPDneLNrNm09nVzUlXzqVn5dqx9B3tbXzzkEl1CTqS7o2Iyf2l81xqZmZNbubseesEG4CelauYOXteL0c0hgOOmVmTW7i0Z0DbG8UBx8ysyY0b3TGg7Y3igGNm1uSmHziBjva2dbZ1tLcx/cAJDcpReV5i2sysyZU6Bgz1XmoOOGZmTaa3LtBDLcAUOeCYmTWRYhfo7qU9nHTlXIAhH3DchmNm1kSapQt0OQ44ZmZNpFm6QJfjgGNm1kSapQt0OQ44ZmZNpFm6QJfjTgNmZk2kWbpAl+OAY2bWZJqhC3Q5rlIzM7O6cMAxM7O6cMAxM7O6cMAxM7O6cKcBM7MhqhmWjR4IBxwzsyGomedM642r1MzMhqBmnjOtNy7hmJkNEfkqtOglTTPMmdYbBxwzsyGgWIXWm2aYM603rlIzMxsCylWhFTXLnGm9cQnHzGwI6KuqTOBeamZmVh3jRnfQXSbojB/dwe0n7g+kare9z7yxabtJO+CYmQ0B0w+csF4bjoD37DKGzq5uTr/6IZ5bsXLNvmbsJu02HDOzIeDgPcZz6J7jUW5bAJfd82emX37/OsGmpNm6SbuEY2Y2RNz06JL1ukOvXN1bB+mkmbpJu4RjZjZEDCZ4bNnRXoOc1IYDjpnZEDF604EHj+WvvEpnV3cNclN9DjhmZkNAZ1c3y8q00/Rn5apomnachgUcSV+VFJK+m9smSTMkLZTUI+lmSRMLx20laZakZdlrlqTRhTSTJN2SnaNb0qmS8m1xZmZDyszZ81jdx/62Ph5hzdKO05CAI2kv4CjggcKuE4Djgc8DbweeBq6XtHkuzaXA24D3A+/L3s/KnXsL4HpgcXaOLwDTgeNqcS9mZtXQX9BYFdFr0GmW6W7qHnAkbQlcAvwz8Fxuu4BjgTMj4oqIeBA4Etgc+ESWZldSkDk6Iu6IiDuBTwMflFSa7+FwYFPgyIh4MCKuAM4CjnMpx8yGqkqCxqpYv8daM01304gSzoXA5RFxY2H7jsBY4LrShojoAW4F3pVtmgK8CNyRO+52YHkhzW3ZsSWzgXHADsXMSDpa0hxJc5YsWTLYezIz2yDTD5xA+4jKvhO3SYg0C8E3D5nUNAM/6zoOR9JRwBuBI8rsHpv9XFzYvhgYn0uzJGJtmI+IkPR07vixwF/KnKO0b35+R0RcSAqCTJ48ue8O72ZmNTRq5EZlB3gWrY5g/pkH1SFH1VW3gJNVeX0DmBoRr/SRtPjQV2FbuaDQXxr1st3MrGFK6990L+1Z7yHWl2ZpsymqZ5XaFGAb4EFJr0p6FXg38Nns/V+zdGMLx23L2hLKImDbfFtM9n5MIU25c8D6pSczs4YorX9TmrCz0mDTTG02RfUMOJ3AJGD33GsO8NPs/WOkYDGtdICkkcBU1rbZ3AmMIgWvkinAZoU0U7NjS6YBC4EF1bwhM7PBqmT9m6LRHe1N1WZTVLcqtYhYCizNb5O0HHg265GGpHOBkyU9SgpAp5A6CVyaneMRSdcCF2TtQQIuAK6JiNLIp0uB04CLJH0deBNwInB6vu3HzKyRBjp2ZnRHOzM+PLFpgw0Mvck7vwV0AOcDWwF3AwdExAu5NIcD/8Ha3mxXAceUdkbEMknTsnPMIXW9Phs4p+a5NzOrUG/r3/Rmac/KpluOoEj+0r/W5MmTY86cOY3OhpkNA6U2nIFWq+UXZBsqJN0bEZP7S+e51MzMGuDgPcbzzUMm9TllTTnNMo1NOQ44ZmYNUOoSXW72gL40a5doGHptOGZmLW+w1WntI9S0XaLBJRwzs7obTJdoAJp8NkgHHDOzOhtsO0wzrX1TjgOOmVmdbUg7TDN3Gui1DUfSfCqcbSEi3lC1HJmZtbjpB05g+uX3s3LVwIeltGqnge/m3o8iLWB2D2nqGEhTyryDNKjSzMwGYhBDIJt5HjXoI+BExJpAIuki4KyI+EY+jaSTgImYmVnFZs6ex8rVA4s4bVJTz6MGlbfhHAL8rMz2nwMfrl52zMxa30DbYTra2zj7o29t6mADlQec5cB+ZbbvB6yoVmbMzIaDLTvaK07bCiWbkkoHfn4HOF/SZOCubNtewJHAjBrky8ysJZ3SOZelPf2v6glpoOfMw5q/ZFNSUcCJiG9JWgB8EfhotvkR4MiIKFfVZmZmBad0zuXiu56sKG0rLEdQVPHUNllgcXAxMxuEzq5uLqkw2AAtF2xgAAM/JY2U9HeSviJpdLZtJ0lb1y57ZmatYebseQPqCX3SlXPp7OquWX4aoaKAI+mNwKPA94AzgFKQ+VfSomlmZtaHgfZM61m5qqmnsSmn0hLOuaQVNl8L5H9rVwHvqXamzMxazWBmCGjmaWzKqbQN513AXhGxSusuFvQkMK7quTIzayGdXd0sf/nVAR/XzNPYlDOQ9XDKdRzfHlhWpbyYmbWcwa590+zT2JRTaZXadaS51EpC0hbA6cAvq54rM7MWMdi1b1plsGdepSWc44CbJM0DRgKXAW8EFrN2XI6ZmRUMph1m/OiOlgs2UPnAz4WSdgc+DryNVDK6ELgkIlqrVcvMrIrGje6ge4BBp9Wq0koqCjiS9gXuiIgfAj/Mbd9I0r4RcWutMmhm1ow6u7qZOXvegIMN0JKlG6i8Decm1o69ydsy22dmZpnOrm6m//z+QQUbgMO/f2f/iZpQpQFHlF8u6DWkmaTNzCwz46qHBrzeTd7tf3y25WYZgH6q1CRdlb0N4GJJL+d2twG7AXfUKG9mZk2p0tmg+zJz9ryWq1rrrw3nr9lPAc+x7iwDrwC/Bb5fg3yZmQ1rrTbLAPQTcCLiHwGypQlmRoQXWzMz60dvbRAD0WqzDEDlbTiXk8bdrEPSWyS9ubpZMjNrbhsabERrdo2udODnhcD5wAOF7W8GjgH2qWamzMyaTakbdDWqwg7fa/uWa7+Byks4bwHuKbP9d8CkSk4g6XOSHpD0fPa6U9JBuf2SNEPSQkk9km6WNLFwjq0kzZK0LHvNKq3Nk0szSdIt2Tm6JZ2qwoyjZmbVVJovrXtpzwaXbj651/Z8/eCKHqtNp9KAs4o05qZoK1LprxJ/Ab5CmqlgMnAj0CnpLdn+E4Djgc8DbweeBq6XtHnuHJdmx78feF/2flZpZza/2/WkKXfeDnwBmM6688CZmVXVYOdLKxrd0d6ywQYqDzi3ACdLaittkLQRcDJQ0SwDEfF/EfHriPhDRDwWEScDLwBTshLIscCZEXFFRDwIHAlsDnwiu96upCBzdETcERF3Ap8GPiipVNl5OLApcGREPBgRVwBnAce5lGNmtVKtHmUffOt2VTnPUFVpwDkBmAr8IavGmgU8Tmq7mT7Qi0pqk/T3wCjSOJ4dgbGkWakByOZou5W0Fg/AFOBF1h33cztp4Gk+zW2F+d1mk9bs2aGXvBwtaY6kOUuWLBnorZiZVa1H2U2PtvYzqKKAExHzSO04l5KmuHkNcAnw1oh4pNKLZe0rLwIvk5ar/khEzCUFG0hVYXmLc/vGAksiYk0Vafb+6UKacucgl6Z4bxdGxOSImDxmzJhKb8XMbI3pB06go71tnW2DqVJpxbE3eRUvwBYRT5Gq0DbEPGB3YDRwKPAjSfvlL1NIX+zOXq49rr806mW7mdmg5XuljRvdwaF7juemR5ewcGkPW3a08/xLK4kBPnVacexNXqVVaqXSyXcl/UrSdtm2gyXtUek5IuKVrA1nTkScBNwHfAlYlCUplkK2ZW0JZRGwbb4tJns/ppCm3Dlg/ZKPmdmgFHuldS/t4Yp7u5l+4ATmn3kQr7y6ioFOpdY+Qi059iavooAj6QBSF+jxwHuBUhjeCThtA6+/CTCfFCym5a45ktRuVGqzuZPU5jMld/wUYLNCmqnZsSXTgIXAgg3Ip5nZGuV6pfWsXMXM2fM4pXMuK1auHvg5D3trS469yau0hPM14LiI+AhpDrWSm4F3VHICSWdKmipph6y09E1gP9IibgGcC5wo6RBJuwEXkToJXAqQtRVdC1wgaS9JU4ALgGuyNiaytCuAiyTtJukQ4ETgnHzbj5nZhuitraV7aQ+X3PXkgM+31abtLR9soPI2nInAr8psf5by6+SUMxa4OPu5jDRrwfsjYna2/1ukktP5pPE9dwMHRMQLuXMcDvwHa3uzXUWa6QCAiFgmaVp2jjmkCUfPBs6pMI9mZv3qaxXPwXyzHS5fhysNOM+RqtMWFLa/jTSgs18R8al+9gcwI3v1luZZ4JP9nGcusG8leTIzG4zpB07gpCvnVmWwJ8CyKixn0AwqrVK7FJgp6XWkAL6RpHcD3wZ+XKvMmZkNRQfvMZ5D96xeFVir904rqTTgnEJq2H+C1HD/MGlqmt8CZ9Qma2ZmQ1e1BmkOh95pJRVVqUXESuBwSf9GqkYbAXRFxOO1zJyZ2VDVWxvOgA2jSbcq7RZ9qqRNI+JPEXF5RPwsIh6X1CHp1Fpn0sxsqBlRpUCxclUwc/a8/hO2gEqr1E4jVaUVbcqGjcMxM2tKAx3Y2ZdWn9KmpNKA09uKqXuQukabmdkgDZdOA3224Uh6gRRoAviTpHzQaQNGkibhNDMbVjraR9AziBkF1j9PmzsNZI4hlW5+SJq4c1lu3yvAgmxdGjOzYaOzq5tXq1CnJuCbh0waFrMMQD8BJyJ+BCBpPnBH1lvNzGxYmzl7HitXbXjAOXyv7YdNsIHKu0XfIum1ko4gTdj5bxHxjKS9gYURMb+muTQza6DiUgTV6hLd6guuFVUUcCTtCfyGNPhzIjATeIY0E/ObyJaBNjNrFaUgUwwu3Ut7eu1FNVDDpXdaSaW91L4NnBcRe5BW6yyZDexd9VyZmTVQfr2bcqrVI3q49E4rqTTg7An8qMz2p4DXVi87ZmaNV269m2obTr3TSioNOD2kJQOKdgGerl52zMwar9ZVXW3SsOqdVlJpwPk/4DRJm2SfQ9IOwFnAFTXIl5lZw9SyqqujvY2zP9r6q3uWU2nA+TJpobUlpOlsfgv8AVhKmknazKxlTD9wAh3tbTU596F7jh+WwQYq7xb9PLCPpP1ZO1v07yPihlpmzsysEQ7eYzxznniWS+56smodBEquuLebya/felgGnUpX/AQgIm4krYNjZtbSbnp0SdWDDUDPylXMnD3PASdvIMsORMS/Vyc7ZmZDQy07Dgy38TclfZVwDit8fj2p/WZh9nkcsAJYADjgmFlLKA34rEXppmS4jb8p6TXgRMSk0ntJ/wj8A3BkRDyZbdse+F/gklpn0sysHkoDPms5Bmc4jr8pqbSX2qnAsaVgA5C9Px4vwGZmLaLWAz6H2+zQRZV2GngtUK4MOBLYpnrZMTOrv97mTaum9hFi5mHDc/xNSaUlnOuB70vaS1Jb9toLuCDbZ2bWlPqbN61ahnuwgcoDzr8AfwbuAF7KXrcD3cBRtcmamVnt1WPetPGjO4Z9sIHKB34uAT4gaWdgV1JV5CMR8VgtM2dmVmv16KI8XDsJFA104OfjwOM1youZWd2N3rSd51bUdjFjl26SSqvUzMxaUtRywA1pZmhLHHDMbFhb1lPb0s2qWke0JuKAY2bDWq1H/Y8fprMKlFO3gCPpJEm/k/S8pCWSrpa0WyGNJM2QtFBSj6SbJU0spNlK0ixJy7LXLEmjC2kmSbolO0e3pFMll2vNbH21btAvd/7Orm72PvNGdjzxl+x95o10dnXXNA9DRUUBR9KbJU3IfZ4m6eIsiFS6aMR+wH8B7wL2B14FbpC0dS7NCaTZCz4PvJ20muj1kjbPpbmUtETC+4H3Ze9n5fK2BWls0OLsHF8ApgPHVZhPMxsmOru6Of3qh2p6jWKHgfy4nwC6l/Zw0pVzh0XQqbSE8wNgDwBJryOtALo18Dng65WcICIOjIj/jYgHI2IucAQwBtg7O6+AY4EzI+KKiHgQOBLYHPhElmZXUpA5OiLuiIg7gU8DH8wFxMNJk4wemV3rCtLKpMe5lGNmJaUHfy17qJXrMFBu3E9pyYJWV2nA2RX4ffb+MODuiPgAKWh8fJDX3jy7/nPZ5x2BscB1pQQR0QPcSioVAUwBXiQNQC25HVheSHNbdmzJbNLs1jsMMq9m1mLqMeDz4+/8m/W29TbuZzgsWVBpwGkDXsnevxf4Vfb+j6R51gbjPOA+4M7s89js5+JCusW5fWOBJRFru31k758upCl3jvw11pB0tKQ5kuYsWbJkMPdhZk2oHg/4rx88ab1tvXVSGA5LFlQacB4E/lXSVFLAuTbbPh54ZqAXlXQOsA9waEQUv2IU+xCqsK1cH8P+0qiX7UTEhRExOSImjxkzpt+8m1lraFTvtOkHTqCjfd2m7+GyZEGlAecrpDnTbgZ+krXBAHwYuGcgF5T0HVI13P4R8afcrkXZz2IpZFvWllAWAdvm22Ky92MKacqdA9Yv+ZjZMFPqIVbryTp7CyAH7zGebx4yifGjOxApMA2XJQsqnUvtVkljgC0i4rncrgtIq35WRNJ5wN8D+0XEo4Xd80nBYhrwuyz9SGAqqZcZpOq3UaR2mlI7zhRgs9znO4GzJI2MiJeybdNIK5UuqDSvZtZ66rHAWklfAeTgPcYPiwBTVPE4nKzqq03SOyVtkm1bEBFPV3K8pPOBfySVbp6TNDZ7jcrOFcC5wImSDsnG6FxE6iRwaZbmEVJ13gXZUglTSEHvmogodfG4lBQEL5K0m6RDgBOBc/JtP2Y2/My46qG6BJvRHe01v0YzqnQczuaSfk5qnL+D1HaDpO9JmlHhtT5L6pn2G+Cp3OvLuTTfAs4BzgfmANsBB0TEC7k0hwP3k3qzzc7eH1HaGRHLSCWacdk5zgfOzs5rZsNUZ1c3S2s8jU3JjA9P7D/RMFTpbNFnkR7gbwN+m9t+DXAGMKO/E0REv2NgshLIjL7OFxHPAp/s5zxzgX37u56ZDR+1HuBZ8sm9th+W1WWVqDTgfBj4SETcJylfLfUI8IbqZ8vMrHo6u7prvgRBSbmu0JZU2oazFfDXMts3B2pfIWpmtgG+euUDjc6CUXnA+R2plFNSKuV8mnVH/ZuZDSmdXd2sWLm60dkwKq9S+yowO5u5eSPSvGQTgXfgthIzG6JO6ZzLxXc9WbfreSmCvlVUwomIO0hzlW1Mms7mvaRxLVMi4vd9HWtm1gj1DjZQ+6UOml2lJZxSz68ja5gXM7Oq+cndf67r9bbatN290/pRccABkDSONE3MOiUjl3LMbKip99LOB71lu7perxlVFHAk7QFcDOzC2okwS4I0m7SZ2bB106Oebb4/lZZwLgT+TJrAcyHlZ2w2MxsSGrF65nBYz2ZDVRpw3gzsERGP1TIzZmbVMOOq+swqkDcc1rPZUJWOw5lLmcXLzMyGmnrOmVYi3EOtEpUGnK8C35L0t5JeK2nr/KuWGTQzG4hjL7uv7tcM+l6OwJJKq9RuyH5ex7rtN6WVNt1pwMwa7p1nXN+Q63rAZ2UqDTjvqWkuzMw20Cmdc1n8wit1v277CLk6rUKVrvh5S60zYmY2WI2YVQDSQmszPjzR1WkVqnjgp6RJpMk6dwL+KSKeknQw8EREdNUqg2Zm/WlEsDn3Y7s70AxQpSt+HkCaMXo8sD9QqrDcCTitNlkzM+vftHNurvs1d952MwebQai0l9rXgOMi4iNAvpL0ZtKM0WZmDfH408vrer0tNmnj+uP2q+s1W0WlAWci8Ksy258F3C3azBqi3jMK7LztZjxw+vvqes1WUmnAeY5UnVb0NuAv1cuOmVnl6r2Sp0s2G6bSgHMpMFPS60jjbjaS9G7g28CPa5U5M7O+eCXP5lJpwDkFmA88AYwCHgZuBH4LnFGbrJmZ9a7e1Wlbbdpe1+u1okrH4awEDpd0KrAHKVB1RcTjtcycmVlvpv+8vlPYnPahiXW9Xisa0AJsEfFH0hLTZmZ119nVzczZ8+huwFIA7ga94SpdgO2HvewK4CXgD8BlEbGwWhkzMyvp7Orm9Ksf4rkV9Z0FusRzpVVHpSWcMcBUYDXwYLZtN9LknfcChwD/LmlqRNR/qlYza1mdXd2cdOVcelauasj1O9rbPFdalVTaaeB24NfA6yJi34jYF3gdaWzOdcDrgV8CZ9ckl2Y2bM2cPa9hwQbgm4dMcnValVQacL4I/HtErChtyN6fAXwpIl4BzgJ2r34WzWw4a/TSzQ421VNpwBkFbFdm+9hsH8DzDLATgplZfxq5dLPbbqqr0oDzC+AHkg6TtIOk10s6DPgBcGWW5h3AY7XIpJkNX9MPnEB7m+p+XbfdVF+lAeczwGzgYlK36D9l768FPpuleQQ4qq+TSNpX0lWSuiWFpE8V9kvSDEkLJfVIulnSxEKarSTNkrQse82SNLqQZpKkW7JzdEs6VVL9/8Wa2aB1dnWz95k3cuxl97FyVfR/QBVttnGb225qoNKBnyuAz0g6nrQkgYA/RMTyXJpKeqeNIvVy+zHlp8Q5ATge+BQwDzgVuF7ShIh4IUtzKbA98H5St+z/AWbGPIVyAAAUsklEQVQBHwKQtAVwPXAr8HZgAnARsBx3ajBrCo3smfbJvbbn6wdPqvt1h4OBDvxcDgx6tryI+BXZrNOSLsrvy0ogxwJnRsQV2bYjgaeBTwAXSNoVeB+wT0TckaX5NHBbFpTmAYcDmwJHRkQP8GB23HGSzomI+n5VMrMBa1TPtAVnHlT3aw4nlVap1cOOpE4I15U2ZAHjVuBd2aYpwIvAHbnjbieVXvJpbsuOLZkNjAN2qEXGzay6GjGTgNXeUAo4Y7OfiwvbF+f2jQWW5Esp2funC2nKnSN/jTUkHS1pjqQ5S5Ys2YDsm1m1tDWgydU90mpvKHZjLlZ5qbCtXJVYf2nUy3Yi4kLgQoDJkye7us2sAUpzpC1c2sPI9hGsqnPNd3ub3COtDoZSwFmU/RwL/Dm3fVvWllAWAdtKUqmUk7X9jCmkKZZkts1+Fks+ZtZgxQ4CPXVe42arTds57UMT3SOtDoZSldp8UrCYVtogaSRpDrdSm82dpJ5uU3LHTQE2K6SZmh1bMg1YCCyoRcbNbPAaOXXN+NEddJ16gINNndQ14EgaJWl3Sbtn194++7x9VmI5FzhR0iGSdiN1Z36R1BWaiHiENPbnAkl7SZoCXABck/VQI0u7ArhI0m6SDgFOBNxDzWwIauTUNY2eNme4qXcJZzLQlb06gNOz9/+e7f8WcA5wPjCHNJ3OAbkxOJC6Pd9P6s02O3t/RGlnRCwjlWjGZec4nzT+5pxa3ZSZDd7I9sZVtDRy2pzhqK5tOBFxM2sb8MvtD2BG9uotzbPAJ/u5zlxg38Hk0czqp7Oru+5tNiWeuqb+hlKnATMbRk7pnMvFdz1Z12uWurOOH93B9AMnuO2mzhxwzKzm8t2ex43uYIfXdHD7H5+tax4cZBrPAcfMaqqzq5vpP7+flatTn53upT11nUnAc6MNHUOpW7SZtaCTrnxgTbCpt802bnOwGUIccMysZhrZKaBthDjjIw42Q4mr1MysJjq7ujn+Z/fX9ZqbbDSCV15dzTi31wxJDjhmVnWl6WrqOSfa3jttzSVHTek/oTWMA46ZVU1nVzenX/0Qz61YWdfrnvux3V2aaQIOOGZWFZ1d3Uy//P66Lwc9fnSHg02TcKcBM6uKmbPn1T3YeLaA5uISjplVRb3G1nS0j+Clle4Y0IwccMxswEozB3Qv7Vlv9cNaEXC4B3E2NQccMxuQ4swB9Qg2bRJnf/StLs00OQccM6tIvlRTTx3tbXzzkEkONi3AAcfM+nVK51wuuevJupRm8rz8c2txwDGzPjViGQEHmtbkgGNmazSq2gxgdEc7Mz7sINPKHHDMDGjcwE3PEjB8OOCYDXOHf//Oui+GVuJZAoYXBxyzYaizq5sZVz3E0p76znmW51kChh8HHLNhoLOrmxMuv59X6lxdlrfZxm20t41gWc9KzxIwTDngmLWwaefczONPL29oHrbatJ2uUw9oaB5saHDAMWsRjWyL6U17mzjtQxMbnQ0bIhxwzJpQWk3zPhpYQ9arEYIIXG1m63HAMRviGjHwcrA8DY31xQGnSk6/+iEeXvh8o7NhTe6BvyylZ+XqRmdjUDZuG8G40SP5yT1P8pN7miNA2lpvHrdFzas/HXDMGuDu+UOrrWVDjRD8zdYdbDNqk0ZnxYYwB5wqccOolXR2dXPsZfc1OhtVs8lGI3j51dW0SayKWPMzb3XASytXc9mnpzQol9YMHHDMBmCHE3/Z6CzUTW8Lnu3Yy+9gYQPmX7Pm4oBjw9JwChwD0SaxOqLPHmbjRneUndxz3OiOemTRmljLBhxJnwWmA9sBDwHHRsRtjc2VleOH/9CxOoL5Zx7UZ5rpB07gpCvn0rNy1ZptnqbGKtGSAUfSx4DzgM8Cv81+/lrSmyNiyHSfKfegXdDPf/YNPU9/af3wH94qKaWUSj0zZ89j4dIej7exiiliCI4c20CS7gYeiIijctseBy6PiJN6O27y5MkxZ86cemSxzwf7QILOQM7TX1oHm+HNY2hssCTdGxGT+0s3oh6ZqSdJGwN7AtcVdl0HvKv+OTIbGvbeaWvGZyWYNmmdn+NHdzjYWM21YpXaNkAbsLiwfTHwt8XEko4GjgbYfvvta545s2pQNn1MOeNHd/CeXcZw06NLXOVlQ0orBpyS4n9HldlGRFwIXAipSq0O+TIDoKN9BCPb23huxco1Y1scLKyVtWLAeQZYBYwtbN+W9Us9Zn0SnoTSrFpaudPA/RFxdG7bY8AVQ6XTALiX2mDsvdPWXHKUR7ObDSWVdhpo1YDzMWAWqTv07cBngH8GJkbEE70dV++AY2bWCioNOK1YpUZEXCbpNcAppIGfDwIf6CvYmJlZbbVkwAGIiP8C/qvR+TAzs6TlxuGYmdnQ5IBjZmZ14YBjZmZ10ZK91AZL0hKgER0LtiGNH2olvqfm0Ir3BK15X0P5nl4fEWP6S+SAMwRImlNJl8Jm4ntqDq14T9Ca99UK9+QqNTMzqwsHHDMzqwsHnKHhwkZnoAZ8T82hFe8JWvO+mv6e3IZjZmZ14RKOmZnVhQOOmZnVhQOOmZnVhQNOA0n6rKT5kl6SdK+kqY3O02BJOknS7yQ9L2mJpKsl7dbofFWTpK9KCknfbXReNpSk7ST9KPtbvSTpYUnvbnS+BktSm6Sv5f4/zZf0dUlNNUGxpH0lXSWpO/u39qnCfkmaIWmhpB5JN0ua2KDsDpgDToNka/acB3wD2AO4A/i1pO0bmrHB2480O/e7gP2BV4EbJG3dyExVi6S9gKOABxqdlw0laTRpnSgBBwG7Ap8Hnm5kvjbQV4DPAV8AdgG+mH3udcHFIWoUaTmVLwI9ZfafABxP+nu9nfQ3u17S5nXL4QZwL7UGyVYlfSAijsptexy4vK9VSZuFpFHAMuDgiLi60fnZEJK2BH5PCjinAg9GxDGNzdXgSfoG8O6I2LvReakWSdcAf42II3PbfgS8JiI+2LicDZ6kF4FjIuKi7LOAhcB3I+KMbFsHKeh8OSIuaFReK+USTgNI2hjYE7iusOs6UgmhFWxO+vf1XKMzUgUXkr4I3NjojFTJwcDdki6T9LSk+yQdkz3QmtVvgfdI2gVA0ptJJe1fNTRX1bUjMJbccyMieoBbaZLnRlPVb7aQbYA2YHFh+2Lgb+ufnZo4D7gPuLPRGdkQko4C3ggc0ei8VNEbSMuvfwc4E9gd+M9sX7O2T51F+pLzsKRVpGfbGdlCjK1ibPaz3HNjfJ3zMigOOI1VrM9UmW1NR9I5wD7APhGxqtH5GSxJE0htbFMj4pVG56eKRgBzclW3XZJ2JrV5NGvA+RjwD8AngIdIQfQ8SfMj4gcNzVn1Ne1zw1VqjfEMsIq131hKtmX9by9NRdJ3gI8D+0fEnxqdnw00hVQafVDSq5JeBd4NfDb7vEljszdoTwEPF7Y9AjRrhxWAmcC3I+KnETE3ImYB59B8nQb6sij72bTPDQecBsi+Ld8LTCvsmkbqrdaUJJ1H+oa5f0Q82uj8VEEnMIn0bbn0mgP8NHvfrKWe24EJhW1vojFrQVXLpqQvcXmraK1n3HxS0Fnz3JA0EphKkzw3XKXWOOcAsyTdQ3oAfAYYB3yvobkaJEnnk9o5Dgaek1T6FvZiRLzYuJwNXkQsBZbmt0laDjwbEQ82JldV8R3gDkknA5eRuuV/AfhqQ3O1Ya4GTpQ0n1SltgdwHPDjhuZqgLLenW/MPo4Atpe0O+nf3JOSzgVOlvQo8BhwCvAicGlDMjxQEeFXg16khtsFwMukEs++jc7TBtxL9PKa0ei8Vfk+byZ1S214XjbwPg4C7gdeIj24vkA2TKIZX6QOA+eSSmk9wJ9I7W8jG523Ad7Hfr38P7oo2y9gBqla9CXgFmC3Rue70pfH4ZiZWV20Uv2mmZkNYQ44ZmZWFw44ZmZWFw44ZmZWFw44ZmZWFw44ZmZWFw441vQkjcwWq2rKaegrJWmX7D6bZmE7STdI2uDBzJLGZDNb/032ecj+LiRNlvRENguA5TjgDDPZf9K+Xhc1Oo8DFREvAdsB11d6jKT3Zfc7qnY5syo6FbgiIv5cy4tI+qik6yQ9k/372KtMmk0l/bekv0p6UdKVkrYr7Y+IOaRF1D5fy7w2Iwec4We73OuoMtu+WO4gSe11yd0AZWsLERGLIuLlRuenGiSNkNRWx+ttXK9rDYakLYBPAfWY9XkUaW2d6X2kOZ80U8NhpMlcxwKdhfWE/hf4nCQ/Y/MaPdWBX417AX+X/gmst30X0nQah5GmzngJ+BfgtaS5t7qBFaRvcYcXjr2LNFfXTOBZ0mSD3yA3bQppKvkHSVOQ/BW4ibQyY2n//yNNkvkSaWbtTmCjbN8i0pxfPyatKHoxMDLL7wcL+f8YaT2el0jza72nsD//+l62773APcBy0jxqdwET+vgdLgJOJk3ouZy0IuMXCmm2Jj0slwDPAzcCu+f2fya7z/9HmsX5VeCNffxdyt5XlmYT0sNuQfb7fQz4UuH3/1PgcuDfsvz+Ods+EvgW8GfSdEt/AD6TO+69wO+yfU9ladtz+zfP/h7Ls/1fBm4o/W5z1zg7+ze0HLibNNlrX/9OPwks7OV3sdsA8rcFac6xUv6OL+Yvl/Z12fn3KmzfJvv7HJrb9sYs7btz2zYFVpKW6Gj4//Wh8mp4Bvxq4B+//4DzR9JknDuSJhbdIXt47U5axOtz2X++fXLH3kUKBKcAOwOHk2bt/Ui2//XZMcdk55sEfJos4GTXe5U0X9SuwFtJ69W3Z/sXkQLBl4Cdsv/svQWcJ4BDsvNcQJrkcFvS4ncfz9K8gfQNdQvSw/oFUoB8Q3bcEcDOffwOF2X3ewJpxuVjsgfNB7L9baQA1glMzn4n3yIF4zFZms9kD8nfkpZEmABs1sffpex9ZWk2A07LrrVjdp8vkPtiQAo4L5AC00Syhzbwi+zcB2f3/97Scdnf6iXSQm27ZGmWkBY5K533h9nx783+rv9HCrD5gHMFcBtpvaSdSBNsvgzs2sfv+ALgql5+F7sNIH8XkeZY2x/YDbgy+9sNJOB8INu+eWH7H4GTCtvuA05u9P/zofRqeAb8auAfv/+A87kKztFJbjJLUsC5qZDmtlIa0lK4q4GxvZzvXrKJCnvZvwj4eWFbbwHn+FyaNtK3/lOyz+/L0ozKpRmXbXvnAH6Hi4CrC9suBm7I3n+AFFw2LqR5lKwkRAo4AUzs51r93lcvx50LXJP7/FNSCSP/7X9Sdu79ejnH2aTSVL6k9BlSSbedVIorfvMfTQqGpdLjm0lfPl5bOPe1wDl95P9a4IJefhe7DTB/B+f2b0kKvAMJOP8ErCiT/g7gvMK2XwHfr/Tf0nB4eXkC68uc/AdJG5Gqj/6OtKTtxqRSwa8Lxz1Q+LyQVLKAVOVxGzBP0nWkhv4rIuKvWR34W0kPyIrz1Yc1y1tHxCpJvyM99MqKiIWSfgrcLOkG4Dek4NZd6XVyn0/O3u9JerA9u24VPyNJ3/BLeiLioX6us971yt2XpM+T2jy2BzpID9x5hXM8EBErc5/3IJXMbuvlmrsCd0T2JM38Njv/jsBWpOCXz9tSSY/k0u9Jajf+Y+F3sQmplNObDlLppS+V5u+eXP6WZdP8D1SU2VZu1c2e7PqWccCxviwvfD6ZVI12LOnb5HLSN8viypcrC5+D9J+diFgp6T2kqqMDgH8FzpS0N+lb/2DyVTUR8XFJbwMOBA4FviHpoIi4aZCnHAH8hVTNVLQs937FIM+/DklHAmeRqqruJn2D/1KZ6xd/h6JvfS1jHBUcD+l3sZIU3Irn6utv+gwpYFQjf72lqdQiYFNJm0fEC7nt5Vbd3JpUrWYZ96CwgdgH+EVEXBoR95Pqw9800JNExOqIuD0iTiN9630OOCz7dno/5R/Og7GmS2vW62syaSllWLta53q9wSLi9xHxzYiYSnpoH1HpdXKfS9f5Pamq7uWI+EPhtWRgt7P+9crc1z7AbRHxvYjoiog/kNqN+vN7Ukloai/7Hwb2LvTE2of0LX4BqXPCqkLetiSVPIrX2KbM7+KpPvLWRR8l0wHm7x25/G1BqpobiHuy8+RX3dyJ1OZVXHVzIumeLeMSjg3EY8BBkqaQGu6PIz1MK16aWNJU0oPgeuBp4O2k7tgPZ0nOAH6Wrdz4M1JAeB+pfrxYcurPFyT9iVRy+iLpW+iF2b4F2c8PZlV7K0j19v8AXEOqBtyZ9KD7ZT/XebekL5MayacBfw98JNv3K9JD5ypJXyH9DseR2naujoi7B3hP/d3XY8BhkqaR/i5HAO8k9crqVUTMlXQV8CNJXyQF/u2B8RFxKakx/nPAeZL+i9Sx4WvAd7K/y7OSZgFnS1pK+tt+jdRel7/GFcAlko4nffvfhtSI/3BEXN1L9mYDp0raIiKe7yVNJfm7OMvfMlKHgtOz/K0p9Uh6DfA3rK0C3lnSS6Reck9HxDPZfX5H0nOkThHnkaqKb82dZxdSCeeGXvI7PDW6Ecmvxr3ov9PAboXt2wBXkRqCF5OCww+Aa3Np7gK+XTjup8Dl2ftJpAfI06xdbfJLhfSHkr7Vvkx6MPyCdbtFH1NI31ungb8nfSN9mRTQ3ls47uvZfawmLe09ntQJYmF2zILsHtv6+B2WukVfTqoWWlTmfrYkjd1YSCpZPQlcArw+2/8Z4JkK/l793lf2u/gR6QvBc6QeXl8HHi339yicv4PUpb10/48Dn87tfy+p/ezl7D776na8iNS7sNgtepMsP/Oz38VT2d/3rf3cexfwz339G60wfz8hfbl4ivSF6XZSUCqlKXXgKL5OzKXZNPv38mx2r78Ativk93Sgs9H/x4fayyt+WsvJvl0+AkyKiAdrfK1FwNcj4ru1vM5wJ+lgUonlLVGlh5akDlL72qkRcX41zpmdd1PSGKYPR5p1wDKuUjOzIS8iOiWVxoP112uwLEnvIPVYm0MqdZ5MalO6vFr5zOwA/JuDzfoccMysKUTEdzbwFCJV872JVJ3XBUyNiGLvsg0SEQ+ztk3SclylZmZmdeFu0WZmVhcOOGZmVhcOOGZmVhcOOGZmVhcOOGZmVhf/H6hF8WRE622MAAAAAElFTkSuQmCC\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'], axis=1).A[:,0]\n",
    "    D[s]['genes_detected'] = np.sum(D[s]['E']>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)\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample1 8786 / 737280\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",
    "    if np.sum(ix) > 0:\n",
    "        print s, np.sum(ix), '/', D[s]['E'].shape[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---  Sample1 ---\n",
      "Pre-filter: 737280 barcodes\n",
      "Post-filter: 8786 barcodes\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",
    "    D[s] = filter_dict(D[s], tmpfilt)\n",
    "    print 'Post-filter: %i barcodes' %D[s]['E'].shape[0]"
   ]
  },
  {
   "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",
    "\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": [
      "8786 8758\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAGPCAYAAAAwSojkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xu8JVV95/3PF1GRi8KEZlriaJOJMggk6tMorTQak1YY8jgGn4wKKiSjGI0Qg8pgNKQ1MwZRGDQPEyGZ2KYjMaNMfBQNFy9AP4CYNhrAIDoKXmi6OQq2NLQXzG/+qDr0dnNOr92nz63P+bxfr3rt2qtWrb2qXofeX1at2pWqQpIkaXt2m+sOSJKk+c/AIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJJ2OUlWJ/GecGkWGRikRSzJ4Un+Z5Lbk/wwyYYkVydZPdd9mylJ9k7ytiSfTDKWpJKcOdf9kuY7A4O0SCV5FrAeeDqwBngd8D7gPuCtc9ezGbc/cBZwOPDFOe6LtMvYfa47IGnO/AFdOFheVd8d3JDkwLnp0qy4E/j5qtqQZBlw29x2R9o1OMIgLV7/FvjycFgAqKoN4+tJXpDk0iR3JPlRkm8meWeSRw7uk2RNf1nj55N8NMm9Se5Mcnq//YlJLkuypW/r1UP7P6e/PHBif8lgQ5L7k1yZ5OBRDijJqiSf7T/7vn79mUPH9qPB45M0GgODtHjdDjwlyS836v028ADwp8BpwFXAG4H3T1B3N+DvgbuAM4BbgXOTvBK4ErilL78TeF+Sp07Qxn8G/h/g3cC5wDOAzyb5ue11MslLgcuAn9JdUnkr3eWHzyR5RuMYJTV4SUJavM6h+xL/xyTrgXXAZ4FPV9UPB+qdUFX3D7y/MMnXgLcnOaOqvjOw7eHAh6vqjwGSrAU2ABcBr6yqv+zLPwp8iy6MnDrUr38NHFxV3+/rfhb4NPAGussoD5FkL+ACYG1VnTxQfiHwZeAdwK+OdFYkTcgRBmmRqqrPACuBjwNPpvtCvhTYlOS3BurdD5BktySPSbI/cA0Q4GkTNP0XA/tuofvC/imwdqB8A/Btussiw/5qPCwM9PPLwL/fzuGsAvYDPphk//EF2BP4FLAyycO3s7+kBkcYpEWsqq4DXpjkYcAvAb9OFxz+Msk3q+ozSQ4F3gU8B3jUUBP7Dr3/SVXdOVS2GdhYVT+ZoHy/Cbr1tQnKvgr8ynYO5Un96xXbqbMvMLad7ZK2w8Agiar6Kd0thl9Mci3dJYCXJfkC3WWK++kuB3wd2Ar8PN2tmMOjlP8yyUf8dJLyTNSdEesNGu/HycAdk9TZ3GhD0nYYGCQN+3z/eiDd/9UvAZ5TVVePV0iyagY//0kTlD0R+OZ29vl6/zpWVZ+a/i5Jcg6DtEgleW6Sif4NGJ8r8BW2jQw8+H/4/T6nz2DXXpHkwUsdSZ4LHEp398VkLgO+D7x1+HbPvo0l095LaZFxhEFavN4L7J3k7+hud9yNbhLjy4HvAecDP+jXP5DkT4Gf0N3yuPcM9msTcG2S/0E37+D1fdm5k+1QVfcmOQX4G+CfknwQ2Ag8jm6U5D7g2PH6SV7Xtz0eTH4lyfi/h39aVV6+kIYYGKTF643Ai4DnA/8JeCTdLZAfBP5rVd0OkOQ4ui/rtwFbgEuAPwNunKF+vRM4uO/ffsC1wKkT/cDUoKr6cJI76OZa/D7dHRJ3AjcwcOdG743AEwbeP69fAP4a5ztID5EqH/gmae4leQ7dBMuXVtWH5rg7koY4h0GSJDUZGCRJUpOBQZIkNTmHQZIkNTnCIEmSmrytcsD+++9fy5Ytm+tuSJI0a77whS98t6qaP25mYBiwbNky1q9fP9fdkCRp1iTZ3s+uP8hLEpIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJrf/v7MbpEkzSmfVqn5beNNc90DSRKOMEiSpBEYGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOPt55hy878xJT2u/3s46a5J5IkTd2sjjAkeWySDyQZS/LDJP+c5NkD25NkdZINSbYmuSrJoUNt7JdkbZLN/bI2yb5DdQ5PcnXfxh1JzkqS2TpOSZIWmlkLDP2X+rVAgOOAQ4BTgbsGqp0BvKEvP6LfdmWSfQbqXAw8DTgWOKZfXzvwOY8GrgQ29W2cBrwJOH0mjkuSpMVgNi9JnAHcWVWvGCi7bXylHwF4PXB2VV3Sl51EFxpOAC5McghdSDiqqq7r67waWJfk4Kq6FTgR2BM4qaq2Ajf3+52e5Lyqqhk/UkmSFpjZvCTxQuCGJH+b5K4kX0ryuoFLBQcBS4Erxnfov/CvAZ7ZF60AtgDXDbR7LXDfUJ11/b7jLgcOBJZN7yFJkrQ4zGZg+AXgtcA3gOcD7wHOBn633760f900tN+mgW1LgbHBUYJ+/a6hOhO1MfgZD0pySpL1SdaPjY3t6DFJkrQozGZg2A34x6p6c1V9sareD7yXbYFh3PAlgwyVTXRJoVUnk5RTVRdV1fKqWr5kyZLWMUiStCjNZmC4E/jnobJbgMf36xv71+FRgAPYNkKwEThg8I6Hfn3JUJ2J2oCHjjxIkqQRzGZguBY4eKjsScA3+/Xb6L7sV41vTLIHsJJtcxauB/amm6cwbgWw11Cdlf2+41YBG4Dbd/YgJElajGYzMPw34Mgkb0nyi0l+k+6WxwvgwbkI5wNnJjk+yWHAGrpJjhf3dW4BLqO7Y+LIJCuAC4FL+zsk6OveD6xJcliS44EzAe+QkCRpimbttsqq+ockLwTeAfwh8K3+9b8PVDsHeBRdiNgPuAF4XlXdO1DnRLq5D+N3U3wMeN3A52xOsqpvYz1wD3AucN4MHJYkSYvCrP40dFV9Apj0t5L7EYDV/TJZnbuBlzU+5ybg6Cl1UpIkPYQPn5IkSU0GBkmS1GRgkCRJTQYGSZLUZGCQJElNBgZJktRkYJAkSU0GBkmS1GRgkCRJTQYGSZLUZGCQJElNBgZJktRkYJAkSU0GBkmS1GRgkCRJTQYGSZLUZGCQJElNBgZJktRkYJAkSU0GBkmS1GRgkCRJTQYGSZLUZGCQJElNBgZJktRkYJAkSU0GBkmS1GRgkCRJTQYGSZLUZGCQJElNBgZJktRkYJAkSU0GBkmS1GRgkCRJTQYGSZLUZGCQJElNBgZJktRkYJAkSU2zFhiSrE5SQ8vGge3p62xIsjXJVUkOHWpjvyRrk2zul7VJ9h2qc3iSq/s27khyVpLM1nFKkrQQzfYIw63AYweWwwe2nQG8ATgVOAK4C7gyyT4DdS4GngYcCxzTr68d35jk0cCVwKa+jdOANwGnz8zhSJK0OOw+y5/3QFVtHC7sRwBeD5xdVZf0ZSfRhYYTgAuTHEIXEo6qquv6Oq8G1iU5uKpuBU4E9gROqqqtwM39fqcnOa+qahaOUZKkBWe2Rxh+ob9McFuSDyX5hb78IGApcMV4xf4L/xrgmX3RCmALcN1Ae9cC9w3VWdfvO+5y4EBg2TQfiyRJi8ZsBoYbgJPpLie8ii4gXJfk5/p16C4lDNo0sG0pMDY4StCv3zVUZ6I2GKjzM5KckmR9kvVjY2M7ekySJC0Ks3ZJoqr+fvB9ks8B3wBOAj43Xm1otwyVTXRJoVUnk5SP9+si4CKA5cuXe8lCkqQJzNltlVW1Bfgy8ERgfF7D8CjAAWwbIdgIHDB4x0O/vmSozkRtwENHHiRJ0ojmLDAk2QP4d8CdwG10X/arhravZNucheuBvenmKYxbAew1VGdlv++4VcAG4PZpPwhJkhaJ2fwdhncneXaSg5I8A/gI3Zf9B/q5COcDZyY5PslhwBq6SY4XA1TVLcBldHdMHJlkBXAhcGl/hwR93fuBNUkOS3I8cCbgHRKSJO2E2byt8nHA3wD7A2N08xaOrKpv9tvPAR4FXADsRzdJ8nlVde9AGycC72Xb3RQfA143vrGqNidZ1bexHrgHOBc4b4aOSZKkRWE2Jz2+pLG9gNX9Mlmdu4GXNdq5CTh6x3soSZIm47MkJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1zVlgSPIHSSrJ/ztQliSrk2xIsjXJVUkOHdpvvyRrk2zul7VJ9h2qc3iSq/s27khyVpLM1rFJkrTQzElgSHIk8CrgxqFNZwBvAE4FjgDuAq5Mss9AnYuBpwHHAsf062sH2n40cCWwqW/jNOBNwOkzcSySJC0Gsx4YkjwG+CDwn4B7BsoDvB44u6ouqaqbgZOAfYAT+jqH0IWEU6rquqq6Hng18OtJDu6bOhHYEzipqm6uqkuAdwKnO8ogSdLUzMUIw0XAR6rqM0PlBwFLgSvGC6pqK3AN8My+aAWwBbhuYL9rgfuG6qzr9x13OXAgsGx6DkGSpMVlVgNDklcBvwj84QSbl/avm4bKNw1sWwqMVVWNb+zX7xqqM1Ebg58x2KdTkqxPsn5sbGzUQ5EkaVGZtcDQXzJ4B3BiVf14O1Vr6H2Gyoa3j1Ink5RTVRdV1fKqWr5kyZLtdEuSpMVrNkcYVgD7AzcneSDJA8Czgdf269/r6w2PAhzAthGCjcABg3MR+vUlQ3UmagMeOvIgSZJGMKXAkORRSX4tyRN2YLePAocDTxlY1gMf6te/Svdlv2rgc/YAVrJtzsL1wN504WPcCmCvoTor+33HrQI2ALfvQH8lSVJvpMCQZE2S1/brjwA+Tzc58dYkx47SRlV9v79r4cGFbrLi3f37As4HzkxyfJLDgDV0kxwv7tu4BbgMuDDJkUlWABcCl1bVrf1HXQzcD6xJcliS44EzgfMG5z5IkqTRjTrC8Hzgc/36C+hudVwKrO6X6XIOcB5wAd3ow2OB51XVvQN1TgT+iS6wXN6vv3x8Y1VtphtROLBv4wLg3L5dSZI0BbuPWG8/ujsRoPsdhEuq6q4kHwLeMtUPr6rnDL0vGiGkqu4GXtZo9ybg6Kn2S5Ik/axRRxg2AocleRjdaMOn+vK9gZ/MRMckSdL8MeoIw18Cf0s3cfCnwKf78mcAX5mBfkmSpHlkpMBQVW9P8mXg8cCHB35H4QG6n12WJEkL2KgjDPTPZBgu+8D0dkeSJM1HkwaG/nbEkVTV/5qe7kiSpPloeyMMHxmxjQIeNg19kSRJ89SkgaGq5uJJlpIkaR4yFEiSpCbnMEiSpCbnMEiSpCbnMEiSpCZDgSRJaho5MCQ5NsknktyS5N/0Za9M8qsz1z1JkjQfjBQYkpwI/E/gq8Ay4OH9pocBZ8xIzyRJ0rwx6gjDGcCrqur36Z4fMe5zwFOmvVeSJGleGTUwPBG4foLyLcCjp687kiRpPho1MGwAnjRB+dHA16evO5IkaT4aNTBcBLw3ybP69/8myUnAOcCfzUjPJEnSvDHS462r6pwkjwGuBPYAPgv8CHh3VV0wg/2TJEnzwEiBAaCq3pLkvwJPphuZ+Oeq2jJjPZMkSfPGSIEhyVJg96r6DrB+oPxxwE+qatMM9U+SJM0Do85hWAscO0H58/ttkiRpARs1MBwBXDNB+Tpg+fR1R5IkzUejBobdgUdOUL7HJOWSJGkBGTUw3AC8ZoLy3wX+Yfq6I0mS5qNR75J4C/CZJL8MfLovey7wVODXZqJjkiRp/hhphKGqPgesAL4BHA+8CLgNWFFV181c9yRJ0nywI7/D8E/Ay2awL5IkaZ4adQ6DJElaxAwMkiSpycAgSZKaDAySJKnJwCBJkpp2KjAk+Q9JXjFdnZEkSfPTzo4wnA28fzo6IkmS5q+Rf4dhEr8KPHw6OiJJkuavnQoMVbVhujoiSZLmr5EvSSR5ZJLfTvLuJO9KcnKSkZ9UmeR3k9yY5Af9cn2S4wa2J8nqJBuSbE1yVZJDh9rYL8naJJv7ZW2SfYfqHJ7k6r6NO5KclSSj9lOSJD3USIEhyZOBrwLnAc8AjgTOB76a5JARP+s7wH8GngYsBz4DfDTJL/XbzwDeAJwKHAHcBVyZZJ+BNi7u9z8WOKZfXzvQz0cDVwKb+jZOA94EnD5iHyVJ0gRGvSTxHuBLwMur6gfw4JfzX9MFh+e3Gqiq/2+o6C1JXgOsSHIT8Hrg7Kq6pG//JLrQcAJwYR9MjgGOGn/gVZJXA+uSHFxVtwInAnsCJ1XVVuDmfr/Tk5xXVTXi8UqSpAGjXpJ4FvAH42EBoF9/C3DUjn5okocleQmwN3AdcBCwFLhioP2twDXAM/uiFcCWvv64a4H7huqs6/cddzlwILBsR/spSZI6owaGHwL7TlD+mH7bSPr5BVuAHwHvA36jqm6iCwvQXUoYtGlg21JgbHCUoF+/a6jORG0wUGe4T6ckWZ9k/djY2KiHIknSojJqYPg48OdJntWPDjwsyVHAhcDHduDzbgWeQjcH4s+ADyQ5bGD78CWDDJVNdEmhVSeTlHeFVRdV1fKqWr5kyZJG9yVJWpxGDQy/B3wNWEc3ovBD4Gq6iZCvH/XDqurHVfW/q2p9Vb2Zbl7E7wMb+yrDowAHsG2EYCNwwOAdD/36kqE6E7UBDx15kCRJIxopMFTV96vqPwAHA8cDLwIOrqrfqKrNO/n5jwRuo/uyXzW+IckewEq2zVm4nm7Ow4qB/VcAew3VWdnvO24VsAG4fSf6KUnSorZDP9xUVV+jG2nYYUnOBj4BfBvYh+7uh+cAx1VVJTmf7s6Jr9CNXLyVbpLjxf1n35LkMro7Jl5Fd6nhQuDS/g4J+rp/BKxJ8l+AJwFnAm/zDglJkqZu5MCQ5MV0PwV9AEMjE1X1ghGaWEp3G+ZSYDNwI3BsVV3ebz8HeBRwAbAfcAPwvKq6d6CNE4H3su1uio8Brxvox+Ykq/o21gP3AOfS/X6EJEmaopECQ5J30c1V+Czd8P4O/996VZ3c2F7A6n6ZrM7dwMsa7dwEHL2j/ZMkSZMbdYThFcBLq+ojM9kZSZI0P416l8RudHc0SJKkRWjUwHARjUsBkiRp4Rr1ksS+wAn9hMIbgZ8Mbqyq06a7Y5Ikaf4YNTA8mW2XJP7d0DZvV5QkaYEbKTBU1a/MdEckSdL8NeocBkmStIgZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNRkYJElSk4FBkiQ1GRgkSVKTgUGSJDUZGCRJUpOBQZIkNe0+1x3QxJad+Ykp7Xf72cdNc08kSXKEQZIkjcDAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaZi0wJHlzkn9I8oMkY0k+nuSwoTpJsjrJhiRbk1yV5NChOvslWZtkc7+sTbLvUJ3Dk1zdt3FHkrOSZDaOU5KkhWg2RxieA/x34JnAc4EHgE8l+VcDdc4A3gCcChwB3AVcmWSfgToXA08DjgWO6dfXjm9M8mjgSmBT38ZpwJuA02fioCRJWgx2n60PqqrnD75P8nJgM/As4OP9CMDrgbOr6pK+zkl0oeEE4MIkh9CFhKOq6rq+zquBdUkOrqpbgROBPYGTqmorcHO/3+lJzquqmo3jlSRpIZnLOQz79J9/T//+IGApcMV4hf4L/xq6UQmAFcAW4LqBdq4F7huqs67fd9zlwIHAsmk9AkmSFom5DAzvAb4EXN+/X9q/bhqqt2lg21JgbHCUoF+/a6jORG0MfsaDkpySZH2S9WNjY1M5DkmSFrw5CQxJzgOOAl5UVT8d2jx8ySBDZRNdUmjVySTlVNVFVbW8qpYvWbKk2XdJkhajWQ8MSf4b8FLguVX1jYFNG/vX4VGAA9g2QrAROGDwjod+fclQnYnagIeOPEiSpBHMamBI8h66CYzPraqvDG2+je7LftVA/T2AlWybs3A9sDfdPIVxK4C9huqs7PcdtwrYANw+LQciSdIiM5u/w3AB8Ft0owv3JFnaL3vDg3MRzgfOTHJ8/xsNa+gmOV7c17kFuIzujokjk6wALgQu7e+QoK97P7AmyWFJjgfOBLxDQpKkKZq12yqB1/avnx4qfxuwul8/B3gUcAGwH3AD8Lyquneg/onAe9l2N8XHgNeNb6yqzUlW9W2sp7sL41zgvOk6EEmSFpvZ/B2G5i8t9iMAq9kWICaqczfwskY7NwFH71gPJUnSZHyWhCRJajIwSJKkJgODJElqMjBIkqQmA4MkSWoyMEiSpCYDgyRJajIwSJKkJgODJElqMjBIkqQmA4MkSWoyMEiSpCYDgyRJajIwSJKkJgODJElqMjBIkqQmA4MkSWoyMEiSpCYDgyRJajIwSJKkJgODJElqMjBIkqQmA4MkSWoyMEiSpCYDgyRJajIwSJKkJgODJElqMjBIkqQmA4MkSWoyMEiSpCYDgyRJajIwSJKkJgODJElqMjBIkqQmA4MkSWoyMEiSpCYDgyRJaprVwJDk6CQfS3JHkkpy8tD2JFmdZEOSrUmuSnLoUJ39kqxNsrlf1ibZd6jO4Umu7tu4I8lZSTILhyhJ0oI02yMMewM3A78HbJ1g+xnAG4BTgSOAu4Ark+wzUOdi4GnAscAx/fra8Y1JHg1cCWzq2zgNeBNw+jQfiyRJi8bus/lhVfVJ4JMASdYMbutHAF4PnF1Vl/RlJ9GFhhOAC5McQhcSjqqq6/o6rwbWJTm4qm4FTgT2BE6qqq3Azf1+pyc5r6pqFg5VkqQFZT7NYTgIWApcMV7Qf+FfAzyzL1oBbAGuG9jvWuC+oTrr+n3HXQ4cCCybiY5LkrTQzafAsLR/3TRUvmlg21JgbHCUoF+/a6jORG0MfsaDkpySZH2S9WNjYzvRfUmSFq75FBjGDV8yyFDZRJcUWnUySTlVdVFVLa+q5UuWLNnRvkqStCjMp8CwsX8dHgU4gG0jBBuBAwbveOjXlwzVmagNeOjIgyRJGsF8Cgy30X3ZrxovSLIHsJJtcxaup7vTYsXAfiuAvYbqrOz3HbcK2ADcPhMdlyRpoZvt32HYO8lTkjyl/+zH9+8f389FOB84M8nxSQ4D1tBNcrwYoKpuAS6ju2PiyCQrgAuBS/s7JOjr3g+sSXJYkuOBMwHvkJAkaYpme4RhOfDFfnkU8LZ+/e399nOA84ALgPXAY4HnVdW9A22cCPwT3d0Ul/frLx/fWFWb6UYUDuzbuAA4t29XkiRNwWz/DsNVbJuAONH2Alb3y2R17gZe1vicm4Cjp9JHSZL0UPNpDoMkSZqnDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmgwMkiSpycAgSZKaDAySJKnJwCBJkpoMDJIkqcnAIEmSmnaf6w5oei078xNT2u/2s4+b5p5IkhYSRxgkSVKTgUGSJDUt2MCQ5LVJbkvywyRfSLJyrvskSdKuakEGhiQvBt4DvAN4KnAd8PdJHj+nHZMkaRe1UCc9ng6sqao/79+fmuQY4DXAm+euW/PXVCdLTpWTLCVp17LgAkOSRwD/F/DuoU1XAM+c/R5pIqMGlA894nsAvKSvb9CQpLmx4AIDsD/wMGDTUPkm4NeGKyc5BTilf7slya0z0J/vTnObi8aK7mV/+PXvAuSdc9mbXZp/hzvPc7jzPIc7bybO4RNGqbQQA8O4GnqfCcqoqouAi2aqE0nWV9XymWp/MfAc7jzP4c7zHO48z+HOm8tzuBAnPX4X+CmwdKj8AB466iBJkkaw4AJDVf0Y+AKwamjTKrq7JSRJ0g5aqJckzgPWJvk8cC3wO8CBwPvmoC8zdrljEfEc7jzP4c7zHO48z+HOm7NzmKqHXNZfEJK8FjgDeCxwM/D7VXXN3PZKkqRd04INDJIkafosuDkMkiRp+hkYdsKOPq8iybP7ej9M8o0kvzNbfZ3PduQ8Jjk+yRVJxpLcm+SGJC+Yzf7OR1N9dkqSo5I8kOTmme7jfDeF/54fkeTt/T4/SvKtJKfNVn/noymcwxOSfCnJ/Uk2JvnrJMN3uC0aSY5O8rEkdySpJCePsM/hSa5OsrXf76wkmYn+GRimaEefV5HkIOCTfb2nAn8C/GmSF81Oj+enKTz349nAZ4Dj+vqfBP5uMT9cbKrPTkmyH/BXwKdnvJPz3BTP4d8Ax9D98NvBwG8CN85wV+etKfyb+CxgLfAB4FDghcCTgQ/OSofnp73p5tz9HrC1VTnJo4Er6X4y4AjgNOBNdI9HmH5V5TKFBbgB+POhsq8BfzJJ/XcCXxsq+wvg+rk+ll3pPE7SxueBc+f6WHa1cwj8L+CPgNXAzXN9HLvSOQSeB2wG9p/rvs+XZQrn8I3AN4fKfgvYMtfHMh8WYAtwcqPOa4AfAI8aKHsrcAf9HMXpXBxhmIKB51VcMbRpe8+rWDFB/cuB5UkePr093DVM8TxOZB/gnunq165kquewv4toKfBfZq53u4YpnsMXAv8AnJ7kO0m+luS9Sfaewa7OW1M8h9cCj03yf6ezP/ASulFDjWYFsK6qBkcjLqf7GYFl0/1hBoap2d7zKia7/rZ0kvq79+0tRlM5jz8jye8Cj6Mb2lyMdvgcJjmcbmThxKr66cx2b5cwlb/DXwCOAn4ZeBHwOrrLE2tmpovz3g6fw6q6Hngp3SX1ctBbAAAFS0lEQVSIHwNjdD/hf9LMdXPBmex7ZXzbtDIw7JyRnlfRqD9R+WKzo+exq9TN/3gX3RffN2eiY7uQkc5hkkcCHwLeWFW3zUbHdiE78ne4W7/thKq6oaoupwsNL0ryr2ewj/PdyOcwyZOB9wJ/TDc6cQzdl9yFM9nBBWjWvlcW6i89zrSpPK9i4yT1HwC+N62923VM+bkffVhYC7yiqj42M93bJezoOXws3cSy9yd5f1+2G5AkDwD/vqqGh5UXuqn8Hd4J3FFVmwfKbulfH7+d/RaqqZzDNwOfr6p39e9vTHIfsC7JW6rq2zPT1QVlsu8VmIG/QUcYpqCm9ryK63no47VXAeur6ifT28NdwxTPI0n+I/DXdBOCPjJzPZz/pnAO7wAOB54ysLwP+N/9+qJ73soU/w6vBQ4cmrPwpP510Y12TfEc7kkXMgaNv5+R2wIXoOuBlUn2GChbBWwAbp/2T5vrmaC76gK8mO662yuBQ+huJ9oCPKHf/lfAXw3UPwi4Dzi/r//Kfv8XzfWx7GLn8SXAT+huO1o6sPyruT6WXeUcTrD/arxLYkf/DvcGvg18mO6WwGfR3Q734bk+ll3oHJ7c/7f8Gro5Ic+im0j6hbk+ljk8h3uzLcjfD5zVrz++3/4nwKcH6j+GbpThQ8BhwPF0d028YUb6N9cnaFdegNfSpbgf0aXrowe2XQVcNVT/2cA/9vVvA35nro9hPiw7ch779zXBctVs93s+LTv6tzi076IPDFM5h3S/vXBF/w/7HcAFwD5zfRy72Dk8Ffhyfw7vBC4GHjfXxzGH5+85k/z7tqbfvga4fWifw4FrgB/25/CPmIFbKqvKZ0lIkqQ25zBIkqQmA4MkSWoyMEiSpCYDgyRJajIwSJKkJgODJElqMjBImlNJTknyrST/kmT1XPdH0sQMDJImlGRNkkryFxNsO6ffdmn//vb+/WTLVZN8xn50P3j0LuDngXfP4CFJ2gkGBknb823gxUn2Gi9IsjvwcuBbA/WOoHuw1WPpnjoI8PSBsuMnaf8JdA/Bu7Sq7qyqLcMVkuyexGcLSHPMwCBpe24Evgb8x4Gy4+h+hvaq8YKqGquqjVW1kW1PX32wrKruHm44ycnAF/u33+hHIpYlWZ3k5iQnJ/k63c8M75XkmCTrktyT5O4klyc5ZKjNA5N8MMn3ktyf5EtJfmWazoW0qBkYJLX8D+C3B97/NvB+ut+43xl/y0NHI8YfaXwQcALwm8Av0wWUvege3vZ0ut/c3wx8PMkjAPpRkKuBZcBv0P3G/tt3so+ServPdQckzXsXA+9O8kTgXrov+VPZyS/jqtqa5GdGIwD6qw+PAF5eVZsGdrlkcP8kv0X3ZL6nA/8/XcBYCqyoqu/21b6+M32UtI2BQdJ2VdU9Sf6ObmTh+3RPHPzWDE8r+M5QWCDJvwX+GHgGsIRuhHQ34PF9lacCNw6EBUnTyMAgaRR/CXwA2AKcNQufd98EZR+ne4z0q/vXB4B/phuNAHBipDSDnMMgaRSfBn4M7A98dLY/PMnPAYcA76iqT1XVLcA+/Oz/9Pwj8EtJ9p/t/kmLgYFBUlNVFfBLwEFV9aM56MI9wHeBVyX5xSTPBt5HN8ow7mLgLuCjSVYmOSjJC8bvkkjy9CRfSfL0We+9tAAYGCSNpKruraofzNFn/wvwYrrQcjPdjz39Id0tl+N17gOeTXe54uPAl4G3se1ujj2Bg/tXSTso3f84SJIkTc4RBkmS1GRgkCRJTQYGSZLUZGCQJElNBgZJktRkYJAkSU0GBkmS1GRgkCRJTQYGSZLU9H8A3fOdTpLiNpEAAAAASUVORK5CYII=\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",
    "\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'], 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: 8786 barcodes\n",
      "Post-filter: 8758 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]"
   ]
  },
  {
   "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": [
      "18232 / 33694\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",
    "    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": [],
   "source": [
    "# normalize by total counts\n",
    "E = E[:,keep_genes]\n",
    "E = tot_counts_norm_sparse(E)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(8758, 18232)\n"
     ]
    }
   ],
   "source": [
    "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": "markdown",
   "metadata": {},
   "source": [
    "## Save base directory files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "main_spring_dir=os.getcwd()\n",
    "\n",
    "if not os.path.exists(main_spring_dir):\n",
    "    os.makedirs(main_spring_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Option to save the barcode information for later use\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": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# option to save mitochondrial percentages for later use\n",
    "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": 22,
   "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": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving hdf5 file for fast gene loading...\n",
      "FO538757.3_ENSG00000279928 0 / 18232\n",
      "FAM132B_ENSG00000178752 3000 / 18232\n",
      "SUPT3H_ENSG00000196284 6000 / 18232\n",
      "MAPKAP1_ENSG00000119487 9000 / 18232\n",
      "CCDC168_ENSG00000175820 12000 / 18232\n",
      "KCNJ2_ENSG00000123700 15000 / 18232\n",
      "MIR99AHG_ENSG00000215386 18000 / 18232\n",
      "Saving hdf5 file for fast cell loading...\n",
      "0 / 8758\n",
      "3000 / 8758\n",
      "6000 / 8758\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": 24,
   "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": 25,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Sample1']\n",
      "1612367428.28\n",
      "________________ /FullDataset_v1\n",
      "Filtering genes\n",
      "Using 937 genes\n",
      "Running PCA\n",
      "Building kNN graph\n",
      "kNN built in 7.14478 sec\n",
      "Saving SPRING files to /home/I0369218/precision-immunology/CHAMBERLAIN.Mat/10X_PBMCs_Signac_GitHub/FullDataset_v1\n",
      "Running ForceAtlas2\n",
      "213.905304193\n",
      "214.041133165\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAGPCAYAAABWJglCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xt8XGWdP/DPdyaZTNs0F5LSTKBNuAkigUoKaxEVLSJV3HW1P3XddEVdSaGw4gWLi1giosUrKoUE71Ivu8u6iIXK/SqIJlAIKOXatDBJaYYmaZqkaWae3x8zZ5hMzm2u58yZz/v1mleTmTNnnpmk+Z7neb7P9xGlFIiIiMj7fE43gIiIiIqDQZ+IiKhMMOgTERGVCQZ9IiKiMsGgT0REVCYY9ImIiMoEgz5RCRKRK0REpd13n4jc51CTiKgEMOgTFYiIDIhIV+Lr+0REGdzeUsA2pL7OjIi8JiJ9IvJ9ETk+x3N/WUQ+kK+25qoY7RGRU0Vkk4j8VUQOJD7XpkK+JlE+VTjdACIvEpE2AEsB3Jpy9yCAL+oc/nyBm3MPgJ8BEAC1AE4C8HEA60RkvVLqO1me98sAfgvg5ry0MnfFaM97AZwH4CkAzwI4oYCvRZR3DPpEhfE+AK8C+GvKfWNKqc0OtOW59NcVkUsB/AHAt0Xk70qp2xxoVym6HsDVSqlJEbkCDPpUYji8T1QY5wDYqjKocy0i5yaGi1vT7m9N3H9uvhqnlIoA+CiAGQCXpbxWQES+KiK9IrJXRCZF5FER+ce0NikAVQA+njJ9cF/isUNE5Nsi8qSI7BORcRG5V0TeqvOeP5wYKh8TkVER6ReRr6QdU5M4347EkPpOEdkoIlV22mNFRLYkpmLE4LGd2mNKqd1KqUk75yVyI/b0ifJMRA4B8BYA16Q95BeRxrT7ppVSY8Vp2WxKqZ0icj+Ad4pITaIdNQDWAvgvAD8FMA/AvwK4WURWKaVuTzx9TeLxPwO4IXHf7sS/RwJYDeAmxKcu6gD8O4C7RWS5UuopABCRMxEfjr8HwJcARAEcC+DtWhtFZB6A+wAckXidFwAsA/B5AMcB0Obwzdpj5beIj8ycBuBPKa9dD+AsAN/P5OKNyM0Y9Iny72wAMQB3pN1/NIA9affdiviogFOeArASQCuAJwHsBbBEKXVAO0BEfghgG4AvALgdAJRSm0XkxwBe1Jmy6AdwlFIqmnKOGwA8A+AzAD6duPt9APYBeE/qsWk+i3hwb1dK/T3lfE8B2CQib1NKPWjRHis3A5hEfOTjTyn3fwhAJeIXBUSewKBPlH/nAHhIpwe/C8An0+5LvwgotvHEvwsBIBF8o0B8qB9ANeLTgPcD+IidE6ZdMAQBLEA8ifAvANpTDh1LPHYWgK0Gp/sw4oF4T9ooyV2Jf98F4EE77TJp77iI3ApgtYh8RikVSzz0EQDPK6X6cjk/kZsw6BPlkYj4AbwHwFU6D08ope7Sud9J1Yl/92l3iMi/I97DfiPiwVpja4hbRHyIr1I4D/Fh+VQvpXx9PeJB/TYRCSMeyH8H4JaU4fQ3IL7awOji6FA7bbLht4hPSZwB4B4RORTAOwF8I0/nJ3IFBn2i/FoB4BDMXqpnl1FQ9WffHEsnIN6zfwkARORfAPwIwC0ArkZ8BcIMgE8A+JjNc64H8HUAv0B8GV0k8RpfAnCUdpBSakhETgLwbsSnRM4G8G8AtorI+xKB34f4nL9R8H3Z7hu1cCviIw8fTbze/0P8c+fQPnkKgz5Rfp0D4AWl1PYsnrs38W9d2v2tObXIgIgsBfAOAH9WSmk9/Y8AeBHAB1KT10TkEzqnMLpI+QiA+5RS56a9XtecEyg1jXjAvTWRIf8NxC8atKS6FwAstDlCknWynVJqSkR+D+CDIrIO8eD/lFLq6WzPSeRGXLJHlF/vA7Aly+dqRXremXb/uuyboy+xwuA3iF/4fy3lIS2hTlKOPRLAP+ucZj+Aep37o5g9LQAROQ3xUZDU+xpSv09cZDye+FY7728BnJK+ZDDx/HkiUp1yl1F77PoNgAYA5wJ4K9jLJw9iT58oTxI95xMAfC6b5yul/iYiDwG4KhGUdwP4R8SnC3JxjIh0IB6IaxCfI/9/iM/nf1Yp9ceUY28B8EEAt4jILQAOA3AB4tXnTko7by+AM0XkC4gPs7+qlLoncY4rROSXiCfZHYP4/P7f8HoOAQD8OJGcdzfiSY6HAbgQ8cqF9yeO+TbiF1K/E5HNiCcDViG+tO/DiFfI+7NFe+y6E/GpiO8kPqs5QV9EWhBfHgi8vrTwYhEZBzCglLoxg9cjKj6lFG+88ZaHG+LBcR+AgM5j9wF4xsY5WhDPZJ8AMAxgE4A3IT50fW7KcVcg0TlOe4370u5TKbco4lMIjyFeQ+B4gzZcgvgQ/xTiS/o6DF7vWMTnv8cT578vcX8A8XyAlxFfCvdXxOfrfw5gR8rzPwTgjwCGABwAsDNxzBFprzM/8frPJI6LJM75FQB1Vu3J8GfYnXjuXw0ePyPtM029Zfx6vPFW7JsoxZoTRPmQWPZ1QCn1QafbQkSkh8P7RPlzP2YXdyEichX29InIs2xsezutlHqtKI0hcgEGfSLyrMRGPGbuV0qdUYy2ELkBh/eJyMvebfH4XovHiTyFPX0iIqIy4cmefmNjo2ptbXW6GUREREXR19c3rJRaZHWcJ4N+a2srent7nW4GERFRUYjIgJ3jWIaXiIioTDDoExERlQkGfSIiojLBoE9ERFQmGPSJiIjKBIM+ERFRmWDQJyIiKhMM+kRERGWCQZ+IiKhMMOgTERGVCQZ9IiKiMsGgT0REVCYY9ImIPKanpwfNzc3w+Xxobm5GT0+P000il/DkLntEROWqp6cHF198MaampgAAg4ODuPjiiwEAnZ2dTjaNXIA9fSIiD+nq6koGfM3U1BS6urocapG3lPooCnv6REQeMjQ0lNH9ZJ8XRlHY0yci8pCmpqaM7if7vDCKwqBPROQhGzZsQDAYnHVfMBjEhg0bHGqRd3hhFIVBn4jIQzo7O3HNNdcgFApBRBAKhXDNNdeUzPCzm3lhFIVBn4jIYzo7OxEOhxGLxRAOhxnw88QLoygM+kRERDZ4YRRFlFJOtyHvli9frnp7e51uBhERUVGISJ9SarnVcezpExERlQkGfSLytFIvpkKUTyzOQ0Se5YViKkT5xJ4+EXmWF4qpuAVHTLyBPX0i8iwvFFNxA46YeAd7+kTkWV4opuIGHDHxDgZ9IvIsLxRTcQOjkZHBwUEO9ZcYBn0i8iwvFFNxA7OREW2on4G/NLA4DxERmUqf09cTCoUQDoeL2CpKxeI8RESUF6kjJkaYHFkaGPSJiMiStomPUeBncmRpYNAnIiLbmByZuyu3/A03PPCCI69dMkFfRP5TRJSIXOt0W4iIyhWTI3Pz4p5x/OxPL2H32AFHXr8kivOIyFsAfBrAk063hYio3HV2djLIZ+mau55DVYUf559xlCOv7/qevojUAvgVgE8B2Otwc4iIiLKyfWgf/vBkGOe+tRWN1VWOtMH1QR/ADQBuUkrd43RDiIiIsvW9O59FdaACnW8/0rE2uDroi8inARwN4HIbx54nIr0i0rtnz57CN46IiMimp14ZxR+fHsInTz8CdfMDjrXDtUFfRI4F8HUA/6qUmrY6Xil1g1JquVJq+aJFiwrfQCIiIpu+e+ezqJ1XiU+97QhH2+HaoA9gBYBGAE+JyIyIzAB4B4ALEt87MyFCRESUgb6BvbjnmVfR+Y4jUROsdLQtbs7evxlAei3dnwF4DvERAMvePxERkdO+d+ezaKwO4NzTWp1uint7+kqpEaXUU6k3APsBvJb43nubBhBRxnp6etDc3Ayfz8cd38h1HnkhgoeeH8badxyF+QHn+9nOt4CIKEvpG8FoO74B4DpycpxSCt+5YzuaaoLoeEuL080B4OKevh6l1BlKqQudbgcRuUNXV9ecnd+mpqbQ1dXlUIuIXnf/s3vQO7AXF77raAQr/U43B0CJBX0iolRGO7txxzdyWryX/ywOr5+HDy9f4nRzkhj0iahkGe3sxh3fyGl3/G03+l8ZxWdWHoNAhXtCrXtaQkSUIbfv+MYkw/IUjSl8945ncWTjAvzzmw9zujmzMOgTUcly845vWpLh4OAglFLJJMNSDfy8gLFvy5NhbN+9Dxe/+w2o8LsrzIoXV74tX75c9famL/EnIiqMnp4edHV1YWhoCE1NTdiwYQO6urowODg459hQKIRwOOxAK7OXvkoCiI+ouOUCy01mojG8+3sPoKrCh9v+423w+aQorysifUqp5VbHuesShIhcjz2+2Yx69HoBHyjNJEOukrDvd4+9gpeG9+Oz735D0QJ+JrhOn4hs47r4uYwCos/nQywWm3N8KSYZcpWEPQdmovj+3c/hpMNrcdbxi51uji729InINvb45jIKfLFYzNVJhpngKgl7fvuXXXhlZBJfeM+xEHFfLx9g0CeiDLDHN5dR4NOSCt2YZJgpt6+ScIPJ6Siuvfd5/MMRh+D0oxudbo4hBn0iso09vrnMAmJnZyfC4TBisRjC4XBJBnzA3ask3OIXj+zAnn0HcImLe/kAgz4RZYA9vrnKJSB65QKmEMamDqL7/hdwxrGLsLz1EKebY4qJfERkm/aHPn15WrkHgM7OzrL/DMrZTx58CSMTB/GFs451uimW2NMnooywx0cAl25q9u6fxk8eegmrTmjCCYfVOt0cS+zpExFRRrh083Xd97+AiekZfO7db3C6Kbawp09ERBnh0s243WNT+PnDO/CBNx+GYxYvdLo5tjDoExE5pFSHyLl0M+6H9zyHmFL47Jml0csHGPSJiBxRyhvy2F26WaoXNXbsjEzgt3/ZhY+eshRLDpnvdHNsY9AnInJAKQ+R21m6WcoXNXZcc9ezqPALLnrX0U43JSMM+kREDijlIXI7tQlK+aLGyrO79+H/tr2Cj5/WikNrgtZPcBFurUtE5IDm5mbPbL2rx+fzQS++iIjuRkSlpPPGXjz8fAQPfPGdqF8QcLo5ALi1LhG5nJfne+3wenVDr5ZsfmLXCG5/ejc+/fYjXRPwM8GgT0RFZ3e+18sXBl4v3+vVi5pv3v4MDlkQwCdPP8LppmRHKeW5W3t7uyIi9wqFQgrAnFsoFEoe093drYLB4KzHg8Gg6u7udrDllInu7m4VCoWUiKhQKFTyP7uHntujWtZvUT9+8EWnmzIHgF5lIz5yTp+Iis7OfK/X57yptCil8IHrHsaesSnc84UzEKz0O92kWTinT0SuZWe+t5Sz28l7bn96N57YNYKLz3yD6wJ+Jhj0iajo7Mz3ejURjEpPNKbwnTu246hFC/DBkw9zujk5YdAnoqKzk8RWyolgXk5ALEf/9/greO7VcXz+rGNR4S/tsMk5fSJyrZ6eHnR1dWFoaAhNTU3YsGGD67Pb03egA+IXK17KzC8nB2aieNe378chCwK45cK3QkScbpIuzukTUUnR6x13dnYiHA4jFoshHA6XRND0ciW6cvSbR3filZFJXPKeY10b8DNR4XQDiIi8tD87ExC9Y/zADH54z/NYcWQD3nZMo9PNyQv29InIcV7qHTMB0Tt++tBLiOyfxhfP9kYvH2DQJyIX8FLvuJQTEOl1kfEDuOGBF3H2m5rw5qX1Tjcnbxj0iUqQ17LDvdQ79np53XJx3X0vYGJ6Bl94zxucbkpeMegTlRgv7lPutd5xKSYg0uteGZnEjY8MYHX74Tj60IVONyevXBv0RWSdiDwpImOJ2yMi8j6n20XkNC/Mf6ePVABg75hc43t3PgsIcPGZ3urlA+7O3n8ZwHoAzyF+cfJxADeLSLtS6klHW0bkoFKf/zbK1L/mmmtYU58c99zuffjdYy/jU6cfgea6eU43J+9c29NXSv1eKbVVKfW8UupZpdRlAPYBWOF024icVOrz314YqSDv+ubt27EgUIELzjja6aYUhGuDfioR8YvIRwFUA3jY6fYQ6YlEIujv70dfXx/6+/sRiUQK8jqlPv9d6iMV5F29O17DnX/bjbVnHIX6BQGnm1MQrg76ItImIuMADgDoBvDPSql+g2PPE5FeEends2dPUdtJFIlEMDAwgOnpaQDA9PQ0BgYGChL4Sz07vNRHKsiblFLYuPUZHLqwCp986xFON6dgXB30AWwHsAzAWwBcD+AXInKC3oFKqRuUUsuVUssXLVpUzDYSIRwOz9kfXilVsDnqUs4OL/WRCq/y2jLQTN3191fRO7AXF5/5BswLlO7WuVZcHfSVUtOJOf1epdSXAGwD8Fmn20WUTuvh272/nJX6SIUXeXEZaCZmojF884/P4MjGBfjw8sOdbk5BldQueyJyD4CwUqrD7DjuskfF1t/frxvgA4EA2traHGgRkX3Nzc0YHBycc38oFCqLFRX//ddd+OL/PonujpNx9gkhp5uTlZLfZU9ENorI20SkNTG3/w0AZwD4lcNNI5qjubl5Tm1uEUmuQSdys3JOrpw6GMV373wWy5bU4T1v8n5eiWuDPoAmAJsRn9e/G8ApAFYppbY62ioiHQ0NDWhpaUEgEM/4DQQCaGlpQUNDg8Mto1LjxNx6OSdX/vzhHRgam8Klq47zzKY6Zlwb9JVS5yqlWpRSVUqpQ5VSZyqlbne6XURGGhoa0NbWhvb2drS1tTHg51G5JJk5NbdersmVe/dPY9O9z+Ndxx2KtxxZHv9fXRv0iYiA8koyc6pwUbkmV1577/PYf2AGl646zummFE1JJfLZxUQ+Imf09PSgq6sLQ0NDaGpqwoYNG3IOHOWUZObz+eYs/QTi+SGxWMyBFnnXrtcm8K7v3IcPvvlwXL36RKebk7OST+QjotJSqB55OSWZGc2hK6U8Pa3hhG/fsR1+n+Cz7/bepjpmGPSJKC8KNTRdTklmenPrGi9PaxRb/8uj+P22MP799CPRVKv/eXsVgz4R5UWheuR2ksy8kuiXOreuhxsT5U4pha/f9nccsiCAzncc6XRzio5Bn4jyolA9cqskM68l+mkllo2Wj3lxWqOY7nt2Dx55MYLPrDwGC4OVTjen6Bj0iWwo1g56payQy77M9hrw6la9bp/WKMXRlZloDN+47e9oaZiPfzl1qdPNcQSDPpGFYu6gV8qcWvbl1UQ/N6+dL9XRlZv6Xsazu8dx6dnHIVBRnuGvPN81UQaKvYNeKUvtkW/YsAFdXV0F7wm6vUecLTevnS/F0ZX9B2bwnTufxfKWepx9Qmn/buSCQZ/IQjF20CvFoVIzxewJurlHnCu3bqFciqMrNzzwIvbsO4DL3vfGsii3a4RBn8iCVk/f7v2ZKtWhUjPF7Am6uUfsVaU2urJ7bAo3PPAizjkxhDcvrXe6OY5iRT4iC9qcfur/FRHJ24Y6Xqw4x8py3qZdqKZe2AWDQddebH3xpidw8+Nh3P35d2DJIfOdbk5BsCIfUZ4Uege9UhwqtVKsnqDbp0WK2b5ivlYpja78fXAM/9P3Mj5+WotnA35GlFKeu7W3tyuiUhEKhRSAObdQKOR007LW3d2tgsHgrPcTDAZVd3d3Sb1GLorZPrd/Fk6JxWKq48d/Videcbsa2T/tdHMKCkCvshEfHQ/Qhbgx6FMp8eof7O7ubhUKhZSIqFAolPf34/aLpWK2z+2fhVPu+ftu1bJ+i/rJgy863ZSCsxv0ObxP5LBSGirNRKEzz90+LWLUjsHBwbwPwbv9s8hFttMWB6MxfO3Wv+HIxgVYs6KlwK0sHQz6RC7g1qVZbp4zd3sGuVk7lMrvKg23fxbZymVly68f3YkX9uzHl977RlT6GeqS7AwHlNqNw/tEuXP7tEMptk/vlo8heLd/FtnKdtpiZP+0OqnrdvUvNzyiYrFYkVrrLHBOn4hyYfQHt66urqBz9ZkodN5ArlLbZxT0RSTvr+XGzyIbRp+b1Wf21T88rVov3aKefmW0SC11nt2gz3X6RKTLaK19Ojevz3aTYtRj6OnpQVdXF4aGhtDU1IQNGzaU9M8lm8/speH9OOt79+NDJx+OjR86sdBNdA2u0yeinNidD3Z7zXW3KHS5YC9WdszmM/vGbX9HwO/D5856Q6GbV5rsDAeU2o3D+0S5szsnjTwOUXtdoYbgu7u7lc/n8+SyvUw+swef3aNa1m9R197zXBFb6A7g8D6H94lylT5cPDExgdHR0TnHlXLJ4FKnVxI3lVHpY69NBcxEY3jvDx7E5MEo7vzsOxCs9DvdpKLi8D4R5ayzsxMbNmxAU1MThoaGICKoqKiYdYxXdrQrVXqbG6XSm6bRmwq48MILUV9fX9DlmYVcAvrrv+zEs7vHcdl7jy+7gJ8RO8MBpXbj8D5RZoyGUPWG+CsrK1Vtba2nssQz4bYsebOVAUbL9oxWZth5brYKuazwtfED6sQrymuJXjpwyR4R6eno6EjO//p8PrVixQrDP8Ys7zqbG9fDG/2MfD6fYbvMLhQK9XMu5O/SV27uV0dcukX9fbB8luilsxv0OadPVEbWrFmDzZs32zo2FAphaGgIen8jynWLXDdug5zNNrdG7yNdPn/OhdpuefvQPrz3Bw/iY6cuxZUfOCGXJpY0zukT0Ry//vWvbR+rJXjpybS8q9PlfPP1+k7WuNfeg4jA7/dDRNDc3AwAGe/doLcUTk8+y/gWolTw8PAwvvjbRzHPD5x9+AwikUjW5yobdoYDSu3G4X0ifbAxpIuUYdd8DGc7PSSez9d3arrDbPlktu8lNTehtrZWVVZWltSc/vDwsLr29w+plvVb1Fd/c7/q7e1VfX19anh4OG9tLiXgnD4RpTNay20WSHJNXHM6LyCfr+/UBYxV4l2+6ven/pw7OjrynrCYzyTIvzy2TbVfcZt6+9f/qB79y19Vb2+v6u3tVU8++WTO7SxFDPpENEdHR4du0FixYoXpH+Nc/lhnWz89X/L9+k5k71sl3uX7s3R6dMaOz//8XtWyfov65daHkwFfu5UjBn0i0pWevd/R0WF6fK4BwEs9fTsKcVFQjJ6+nddzy6qNx597WR31pS1qzaa75gR89vTNb0zkIyozN954I6LRKJRSiEajuPHGG02P1yv+MjU1hfXr19tKjit0zXkrxXz9QtW/N0u8K8R7cTJh0UokEsHXbv07KnyCNScunPVYanIjGbBzZVBqN/b0ifLH7ppus96/0wVtivX6hewhp9ZN0EZqCvVe3NzT/8nWR1XL+i3q8l/dP6eXX65JfEp5YHgfwJcA/BXAGIA9AP4A4AQ7z2XQJ8ofO9Xb3BQUMpHvi4Fi5S8U+iLGrXP6k9Mz6h++ulWd9rWt6pGU5L1ynsvX2A36bh7ePwPAdQBOA/AuADMA7hKRQ5xsFFG5sbumG3DH8K9dhRiKL8Ra9HTF2EK3s7Mz47X/xdB9/wsY2h/Fp95cg0qfzHosEAg41KoSY+fKwA03ANUAogDeb3Use/pUqoo1DJ3p66QfX1dXZ9nTd3pI30omQ9h2kx+L0UN289B7Ib20Z1wdc9lt6ryfPaL6+vpm9fDLeX2+BqU+vD+noUAo8ct9utWxDPpUioo1pFqMgjtGxWTsrBQo1oWC3aF4o2WORu8l09URhWq3l8RiMdXx4z+rN33lj2r36KQaHh5WTz75ZDJbv9wDvlLeDPr/DeBxAH6Dx88D0Augd+nSpfn7JImKJN89OKMAmq/XMQvQZnkAZsl+xZxHtvs5GBU08vl8jryHTH5+bh9tsWvLE2HVsn6L+ulDLzrdFNfyVNAH8F0AYQBH2jmePX3KhVO9iHz24MyCTzF6imYZ/0YXF06sp7cK0N3d3YbvA0Be34PdAG33wsKtyXiZGpucVqdedad67/cfUAdnok43x7U8E/QBfA/AIIDj7D6HQZ+yNTw87Nh8YT6Dntm5ChVcU4OWWblfo4sLJ4atjQJtd3e3Yd6CWU8/2/egF6ArKipUXV2d7kWAnQsEr8z9d93ytGq9dIt6fOdep5viap4I+gC+D2A3gDdm8jwGfdJk2mvXjnWiylc+e2ZmwaoQPUCzDWGc6OnnMqxt973ozdVn+x5qa2stXy/Tn1GuF1FumBrof3lEHXHpFvWfvyvPKnuZKPmgD2AT4mv03wWgKeVWbfVcBv3ylB7gd+zYYbvXnvpco1sxZPOHVu85Zj3t7u7uvP9Bt7uW36qATz4uRgpVNjj1Nm/ePN3zZfPaVlMI2V4A5TrV4PTUwMGZqDrnBw+q9ivvVCP7p4v2uqXKC0Hf6Bf/CqvnMuiXH71heaNbeq/dznP7+vocemfmjP4457PHaIfVHH62SwOd2P421wqEVu8h/XE7vfxMe+na62QbuN0wNfCjB15QLeu3qC1PhIv2mqXMbtCX+LHesnz5ctXb2+t0M6hAIpEIwuEwpqenEQgE0NzcnPzervb29uTX/f39tp7b2tqKhoaGrNpcKM3NzRgcHJxzv8/nQywWM3xeKBRCOBwueDtSX2/Dhg1FKe7i8/mg93dNREw/E43Ve0mV6eeoFdZJ38vArtTX6+npQVdXF4aGhtDU1KT7+do5Rk+un2Gudr02gbO+9wBOO6oBP/74coiI9ZPKnIj0KaWWWx3n5op8RHNEIhEMDAwkg/T09DR27NiRUcBPr9xl97k7duxAX18f+vv7EYlE7Dc6TU9Pj62NauwwqoAXi8VMq+jZDWp2WVXtGxwcxNq1a5Mbopi951w/H6PqdyJi65z5rEDY09OD+vp6iAhEBOeff37WAb+yshITExPw+Xyoq6vDRRddVLCqfMWoLGhEKYUv3/wURICvfuAEBvw8Y9CnkhIOh3V7IHbp7cKVafnO6elpDAwMZBX4811C1eyP8OrVqy3bki9a2Vafz/pPitl7zuXz0S4WjC5oYrGYrXPqlaCtra3VPdbsQqKnpwcXXXQRRkZGkvdl+7srIlBKYXR0NPnvwYMHZx0zNTWFrq6uWa+f7Wfp5M6ItzwRxv3P7sEXzjoWh9XNK/jrlRsO71NJ6evrs32siKChoQFjY2OzpgLSh+i10YNM/y8EAgG0tbVl9ByjoJTtcHtPTw/Wrl2r+1goFDLt0QeDwbzXUzcaFtaj956z/XzMhs2NpjoyGSq3Myyf/nlmMk1QV1cHALM1ebHlAAAgAElEQVQuEDSVlZVYsGCB7mPpUoffc/1dy3ZqIBd790/jzO/ej8Pr5+F3F7wVfh97+XbZHd5n0KeSMTAwgOHhYdvHZzIHH4lEsGvXLkSj0YzalJobYEch5kpzGf4s9tx+Kr33nO3nYxbghoaGTM+pF9D1LohSg6BRe1I/T7sXQNprdXV16b6H2tpajI2N2TqXndcv1rx8Ni7+7ePY8uQgbrnwdBzfXON0c0oK5/TJUzIN+IFAIKOku4aGBixbtgytra22h/uz2dWrEHOloVBI9347FwP53hUvk/lwvfec7edj9D60Xqqe2tpaNDc3Y+3atXN68OlD5UB82D8cDienCazaYTQlAMQDcvrudUbvYWxszNbvR/rwu5Pz8tm462+7cfO2MC5459EM+AWUVdAXEZ+IzM93Y4iMZBLwtXn7SCSC/v7+jJLvGhoa0NbWhvb2dtOgnpobkMnrZDpXaiepbcOGDaioqJhzv52eYb4DQOp8OADDOX6j95ztXLJZgNM7Z2VlJcbHx01HJcwuiOwEVLPP/7rrrkMsFkM4HEZnZyd6enoML9KM3kNFRQXq6uoMt751cl4+U6MTB/Gf/9eP45oW4sJ3Hu10c7zNzrq+xC/vKgC/BDCA+N72UQD7ATwI4DIAzXbPVegb1+l7y44dOyzX0D/++OOzqu4ZPSe1MI9VtT6j9fuPP/548thsyvbmu8a6UsqyZKzerZjFVjJZf59tgaL0z6qyslLV1tYqEVG1tbWzStraWRtvtibdzs/GbL2/1bn0zpmvwk1u9Pn/3qaO/NKt6sldI043pWQhX+v0ReQDAL4JYCGA2wD8BfHNbyYBHALgBABvA7ACwM8BXK6U2mPzmqMgOKfvHZFIBDt27MjoOdXV1RgfH9d9zO/3Y9myZbrJeyKClpaWWdMCejUBUh83WuOfTZJfukwSsTJJoNPOUax188WSOudeW1uL8fFxzMzMJB9Pnae3+rzsJDlaJbrZ/fmZ1Vq47rrrPPUz0nPv9lfxiZ/9FReccRS+ePZxTjenZOUtkU9EHgVwJYDblFKG2R8ichiAzwDYrZT6TobtzSsGfe+wWzgnE+3t7XkL1marCTJN8ktnFpi0BDUt2BglgaWfrxyCCGAdcM0SDvN1QWQ3QbAUE+7yZd/UQZz1vQdQXVWBLf9xOqoq/E43qWTlLZFPKfUPSqktZgE/cdwrSqkvOh3wyVvyHfCtzpvp6xnN+2eS5Gc0b28235669nrt2rXYu3cvKisrTV9HKVUWAR8wT+wDjOe7u7u7k/PsudJb7683elBqCXf51PWHv2H32BS+ufpEBvwisZXIJyJMpaSiGxgYyPs5/f74Hxa7wdoqSa+5uXlOApZeASAjZgVUMsmEn5qamjWUrafUg0gmlfqsAqndgJyr1Ix/o4sJOwl3+azi6Ba3Pz2Em/pexgVnHI03L613ujllw272/l4ROVT7RkSuE5HGArWJCJFIJKOMfbuWLFkCwF6w1iv5m16Jr6GhAS0tLQgEAti6dSvOOeccLF++HG1tbbb+MHd1dRkuF9MLTGbMpuoqKytdmbVtV6bV5Yyy3ScnJ5OBE4BlQNZrR2rwXbNmTc7B2OoCRKvsl/reL7roopII/EYXzXv2HcCXftePNzXX4D9WHuNwK8uLreI8IhID0KSUejXx/RiAZUqpFwvcvqxwTr/0PfbYYzmV2wVmV2Lz+/1YsmRJwZL07M7f6rUxk/ncTIrfpOro6MCNN96Y8fPcor6+XrcinVlxodREu5qaGkxMTMwqXZtpRcJsqvLZZZYUaPTe6+rqsHfv3oxep5iMkmWXLl2KS299CQ88N4xbLzodxyxe6GArvSOvFfl0gv4+ACcx6FMhZFqIx4zf70c0GjUswWsmkyS9bEueZvq8bHdpy3flvWIyKzWc6855mSQ32r3gysfOe6kXD2ZFlnK9MC4ko4vm+3ZO44ePvobLzzkenzr9CAda5k2FqMjn3t8u8pR8DutrZXW1ofmBgQFs27YNfX196Ovrw7Zt2wyL6WSSpGeVOGZEr7BORUWF6VB8VVWV7v1mwSHflfeKKb0yXiq7eQpmuxHa3YTG7meY6WdtNsVTyvQC/tD4DH7UtxenLK3BJ05rLX6jKKOg/w0R+YGI/ABAAMAG7fuU+4lyksuWtVaUUhgeHp5VXz8ajWLHjh26r1tTMzd/1ShJL5cMbL3cAj1aj3B0dHTOY7W1tbpV+TJph1uZBVE7eQpm1e4A/QCrlzhn9zM0K7+rx+qC0eh8mb5OsaVfHM/EFK55dBQ+Af69LYi9e19zqGXlzW7QfwDAUQDaEreHASxN+b4N8SI9RDnZtWtXxs/Jx37b6cOxkUhE90KgoaFBd4rAKNN+YmJiVi8yPZhceumlc7ZIPXjwINauXQu/3481a9Yk79frEQLx4eT58+fPOY9GywQv1QzwXC5YtAslqymA1MBrlDS4cuVKW6spxsfHM/psrS4Yr776at3RoKuvvtr2awDFXwGQniz7m6fG8dxrB3F+ey0a5/tKdrqp1NkK+kqpM5RS77S4vavQjSVvi0QiGe9yB+RnXjN1KFKrAqh33uHhYd2LAS0DW9siVTM6OpocPtYLJmbbpcZiMWzevDkZ+M16hGa94aqqKjz00ENZ763uNLPevNUQuNGFUrrUwGs03H733XfPyrJP/1lrZmZmMhqat1qy19nZiWuvvXZWz766utr2+QHz1Q+FuhjQVrYAwBO7D+Dm7ftx5hHzcNqS+HstVA0OMsd1+uQahViXn4lt27bhscceMy37u3XrVpxwwgnw+Xyoq6tDXV0dfD4fQqEQ7rrrLoyNjc15jjZ8bDcApfv1r38NwLxHaNYbHh0dxebNmz05b2w1f25nfj19TbzZxZW27v766683/VkODg5mFExT8zTq6up0VwAcOHAg+fXIyEhGF21GFzLr16/Py8Wg0dK8hoYGTMT8+MFfRnF4TQU+uez1UJLNLpWUO7vZ+1EAoZTs/esAfEUplf+F1HnA7P3Sk8+M/ULZunUrrrrqqqwCt4jkNCKhlDLN8gaQVVZ/KZR6raur081jALJfHaEtlcy2Zr6dTP5gMGi5hNPuUs9sV4ekv1+7MlmBYLaPRX39Iej40Z/Qu3MUG1c2oKW2ctbjmaymIXP5zt5PnzTtAMDeP+WN2wM+AGzatCmrgA/Ee+NG28xa0Z5nVsQlfUvbTNrlZj09PYYBH9Af+k/tYU9MTMwpTxwMBudsbZt+Tr25+5UrVya/tjOCYGdkxW7mfrarQzSZ/pwzWYEQDofnXFAopRAOh/GjB1/Ewy+N4nPvbMUxixYAiPfwGfCdk91fobkXAURZc2JYP5uhxd27d2f1WtrwcbY96o997GPJr83KumqPGc01G7ULcG+ZV7Pph7q6ujkBO33uenR0FEop033n03V2dmL16tVz7r/pppts7YtgJj2YGo0WDA4O2lo5YLcdRnkDRr8rmbw/o7n5x18Zxzdv3473tjWhc+XxaGtrQ3t7O9ra2hjwHcR1+uQ4J3r52SQRLV68OKvX0oJMNkusVqxYkVElvZ6eHt1thSsrK9HR0aE7SmBV4rYQFwSp56yvr0/mRqSf36zHuXHjxjltW79+/Zye88zMDObNm4dYLJbckdDqvdx9991z7kvtgZvtixAMBjF//nzT964xG/2xWjmQnotgxmiUaOPGjbrvY3Bw0PbPWu8C+rXJKL775xG0NMzH1R86MS8rbCg/uE6fyKZ169bZ3gBHo/2x6+np0U3yM7NixQo8/PDDho/rBeOuri7djXf8fj/uvvtu3TKvZkPMmda8tyP9nCMjI8keefr5jXqcdXV1eOihh7B27dpZbTOaChgaGsrovVgNp6dPp2jBWwumRjseKqVmvaad0R+9lQOhUAirV6+2dQGj0RslMpsW0nZw1GpTGJ0/fWneTEzhO38ewVQU6O5ox8Kg+e6PVFx2E/nug3VPX7ll2R4T+UpHoRP4qqurdXu+2dq6dSs2bdqE3bt3Y+HChRARjI2NWW52o5Sy3AUvlVWt/GzL8WpSk8XMkrxS9y9IlUtZXztJcNr5jRLdVq9ejc2bN9t+TaP3kfpadtqYr8Q57Tx2S/umJ1xmu9eDGau2mJ0/dR+LX/Tvxy3P7MP3P7oM/7TssKzaQpnLayIf1+lTIdgJ+NpWuNnKNOBrO+WdcsopOOecc7B161b4fL5kT2bVqlXYsmULvvrVr2L+/PkYGxvD4sWLTYfuDx48mFHAB4DTTz/d9PFsl/9pUoeqzeZvjQJlJole6SMSdoKcXo86dVhab/jdjFmPOr0wT319vW4bMxlOt5oT115TrwyznfMVonSv1c90amoKV1xxhe420w0NDWhra0O4IoRbntmHc09rZcB3qWwT+YhyZqeHn02xHjN6QT31sauuugpDQ0NQSmFoaAhXXXVVMvCbHbd//37DId1sWP3xzkct/dTAk+m0hd1Er56eHlx44YWzhtTtUEolh5T1hqXN3n9dXV3yIsHOigntvWhb2OoVTKqtrc2oF231maZ+fnbmu1NXDgC5Z/NbtcmIlsyqt810/8ujuOSmJ9DeUo//fO8bs24HFZbl/wgR+bKILLBzMhF5q4i8P/dmEeWfWVAH9JfkTU1N4Yc//OGsiw+942ZmZjB//vysl+Wls/rjnY+ldk1NTck8gKmpqVnz0mYqKytNe7ypPfu1a9dmPMqhMZtzN3v/GzduTF4kWE1fpm5u1NXVZVjKeP78+YYBXy+3Qhuh0BsBSh0xMHvNVKkrB4Dc9nowYufiLzWZVVuWBwCvjk3h07/sRcOCKnR3tCNQwf6kW9n5yRwNYKeI3CAi54hI8rdKRIIicrKI/IeI/AXALwG4d4Nn8iS/32/aW9J695dffrluUN+0aRMA4yV56fcbHTc2Nobrrrtuzh/OysrKOUO4wWAQV155pe6mPoD1H2+zoGtVoU97/ZUrVyYT24D4ELgWkMwCv1UATE2Ws2P+/PmG0yNGm+FMTEzoHt/R0TGrbVafQyAQSB5vdqFl9JhZcmBnZydGRkbQ3d2tu2rC6jVTpX8OVqV7s2FV6yEYDGLdunWz7puensbUwSg+/ctejE4exI/+bTkWLdTfBZJcQilleUN8Q50exAN6FMAMgMnE1zEAvQDOAxCwc75C39rb2xW5244dO1Rvb29Otx07diTP9+STT+oec+WVV6pgMKgQT0TVvYmI6u3tVU1NTbqPNzU1qSuvvFI1NTUpEVE+n8/wuCeffHLWsdpztfsAJJ/f1NSkVq9ePad9wWBQdXd3W36GHR0dc9oQDAaTr2f2nru7u1UoFNJ9LBQKqe7ubtPnh0IhJSLJYzVG57Q6l1JKiYjhz8fsPWvHpLZd093drSoqKkxf307btTams/oMjT6nbD6v1M9Be29W57dzjJHU52q/x+n/v5544gl10a8fUy3rt6it/YO2z035B6BX2YiPtrL3NSLiA3AigBYA8wAMA9imXFaOl9n77pbvjP1AIGC47v6cc86x7E3V1NTgnnvu0S2zGwwGcc4552DLli2miXPBYBCXXXYZVq1aZXiM2fkfeugh7N69W7csLDA7OzoQCKC5uRkNDQ3J4fmhoSEsXrwYp59+evJcgP5mRFrmuFGGuZYpXl9fb7ohUOp7sLMSwIj2ekZJfnV1ddi7dy96enqwdu3ajNqjzdObDaFr7TU6tqKiAtdee63u6IbZ+9Urw7t69epZSydXrlyJ//qv/7I1xJ/Jiok1a9bgV7/61Zy2ZZvhb1Rq986wH9c/9DIuec+xWPfOozM6J+WX3ez9jIJ+qWDQd7e+vr6ivdYpp5xiGYRqa2uT2eBbt27Ft7/97eR675qaGoiIaSlY7Rxf+MIXTIO+2QVIU1MTLrroIjQ0NCSDuHYBsHr1asPa5lpls/7+fvz+97+33BugoqIC1dXVGB0dNay7b7Zczoj2HLM6+VbPtQq6XV1dthMBNVZ7HtTW1uLqq6/G+vXrk+1OfY72uFGQNKvtb2cNfjAYhIhgcnLS9LjKykr88Ic/tBWs16xZY7qcMdvllukXnr2vVaBr64v45zcfhu9++CQW4HGY3aBvvVZk9kmXGjykAEwppfZkcj4qP8Uuubt48WLLnn560ZzU3czsFtQZHR3FVVddBQCGgd+sjO/Q0BA2bNgAEUkGPW1++K677kr23hcvXox169Zh1apVCIfDyaDf3NxsuDeA1hutqanBxMREsveuFwzTt3QFMOsixKxsbHNzs2GP1Sj4VlZWYmJiAj6fD01NTfD7/XPOoW1Vm01mulnAr6iowPvf/35ceOGFs5INU5+jV1kvdXSltrYWFRUVs56f3sM3Y/c4qwvX1DZZHZtthn9DQ0Py9+3e7a/ia7f34vSjG1lxr8RkOrwfg3mRnjEAPwPwRaVUdim7ecCevjtp+9QXk52d8ZqamrBlyxYA9qYDzKSeK7UN3/rWtzKuyGfnta644opkcDYbag6FQpicnNQdrjfbcS6d3XX26bRRhfSAOT4+bju7327v2a6Ojg7cfffdlu8ndXhfb/SjsrISfr8/eZ9Wz97O1EgmjHromRZqyqWwEgA8sWsEH73hzzhy0QL89ry3sOKeS+R7lz3NvwB4GcCXAbw7cfsygJ0APgngCgBrAFye4Xl1icjbReQWEXlFRJSInJuP85Izih3wgXiv+7LLLjPM4k7PSM52Ux2j52/duhVdXV15D/hAvMdmp2QtEO+JGwWhTAJpNmv6U9vW2dmJDRs2oKmpCSMjIxkt58v3FsA33XSTrQuYmZkZXHzxxQD0i+IcPHhw1n0jIyMYHx/PqG6DnWONLkYzKdSUa4b/juH9+OTP/4rGhQH87BOnMOCXoEyD/vkAPquU+oZS6p7E7RsAPg/gk0qp7wP4D8QvDvKhGsBTAD6D+GoBKlHbt2937LW1Knq9vb248sor0dTUBBFBU1PTnOQ7o011ampqZj3PaIlZTU3NrOI/X//617Nep25H+iYwdqq76VE2N9rJdgvflStXoqenB3V1dcl6+dnKVy2ETCoaTk1Noaenx/YokFa3IXWpXkdHh+EFk1LKcojc6KLO7meZa5ne8MgkOn7yKBSAX3ziVBy6MPOLP3JepsP7kwBOVEo9l3b/GxDP4p8vIq0A/qaUMt5mKpuGiowDuFAp9XOrYzm87z7FTN7LhVGGffrFwcaNG3HTTTfNeX4uQ9AigoqKCluZ3OnPi8ViOO200/DII49k9dqpQqEQNmzYYFnb3WioX29OO9P9B2pra00TAjP5nNPn3LOlXehkcsESCoWwcuXKZMZ+TU2N4V4NZrkARgHb7ooGwDovwMzusSl8uOcRvDY+jV99+h9w4uH2tm+m4ilI9r6IPAPgD0qpS9Lu/xaA9yuljhORUwD8n1Lq8EwbbfHaDPolqtCb6gD5ne9N3VRHWwZ3xx13JIfoa2troZQqyJB9KqvMc40WWOxsQFNbWwsRyXq+OXU+2GjTl2AwmPN8djYBVo+I4Prrr88q89/oXFZLAAshGAziwIEDc/IuMtmw5/rrr7fs5estDY0FqvHRGx7B0OgUbvz3f0DLgpju8lFyVqGC/vsA/C+AFwD8FfGkvlMAHAXgQ0qp20TkAgBHK6U+l1XLjV/bNOiLyHmIFwjC0qVL24udJU7GSqWXr0ebky/kEL0Rv99vufeA1gO84IILLC96tHXiN910U9ab9ejt9pa+xPD888/PqVepvc7111+f006Cqe21s17finbBM2/evJzalKts6yI0NTXhqaeeMgzQemvx900rXPmnMYRHp/HLT52KIxcqy+Wj5IyCJPIppW4FcAyAWwDUAKhLfH2sUuq2xDHX5Tvg22zbDUqp5Uqp5YsWLSr2y5NHbdq0yZGAD1hvNpRa0tUq4NfV1aGqqgqbN2/OKWAppZK7DjY3NwPAnM1w8rEvgIjg/PPPR1VVFerq6pLz4tm0V2un3vI7u7QEuJ6eHkcDPmB/h8R0u3fvnrNJTqpwODwrmI9MRXHFfRHs2juFn5y7HKe0HjLnGGB2DX5yv4wzYpRSu5RSX1JKfVAp9c9Kqf9USu0sROOo9JX6iEuu2fyFIiKYmJjA+eefj+bmZtMksI6ODkxNTWVcNMeI9kd/cHAQa9euxZo1a2Y9nm2Gfypts5zR0VGMjIxAKYWJiYnkcrhMaO3M5f1rF1e5bF2bT9nskLh48eJkgI5EIujv75+1TW5qVcs9E1F8+d7XMDgexaVvrcNpRzUCgGHlS6P7yX0yTvUVkcUA1gE4HvHh/acBXK+UcudfR3JMJBIp+Fx+odkp7uMELSAC8aBmNMy7YsUK3H333QXtnW7evBmbN29OJgCmF/XJV9VPrYqgE7Rg75bfhdQevp1iQKlLU7VtcbWfi/a9Np308tgMvvrAa5iaUfjK2+txYnN18jxGJa8DgUA+3hYVQUY9fRF5K4DnAXwM8SV0UwA6ADwnIivy3TgRqRaRZSKyLNHWpYnvjSoDkouUei8fANatW5f1MrhiisVimDdvXnI5m8/nQ0dHBx5++OGcE9jsGhwcxIUXXphc2qcN+2fTOzfiVNlwbTnjvHnzLI+trKxMJkzma3lhqtSphosvvtgyadLn881ZfaI3RA8AL47M4PJ7I5iJAV894xC8cVFVcnoEgO6oUupUD7lfpol8jwDoB7BWKRVL3OcD0A3gBKXUaXltnMgZAO7VeegXSqlzjZ7H7H3nFSNjv1jsVPVzC62OQGqPu9i949ra2lmBKB9JdKUoFArl/YIrdTTFTua+nY2gNI8PHcD3Hh3D/ArgK2+vR2vD/OQSw9RMfQDM3nehQmXvTwJYppTannb/cQAeV0pZXwYXAYO+s5wot1sM2lK+TIZ48106NhN2l/wVQvrr9vT02FphUE6CwSCqqqps5xqkb/5jlrkvIrP2abCy9fn9+Om2fTiuqQY/OXc5QrXzDHfWY6a+OxWqDO8ogCN07j8CQH4LTVPJ8mLAB16v7JdJxrSTQc5OwNeGnws9GtDZ2enYBYgb1dXVYfXq1Yaf+7x58+asNhgdHbVVdjkUCiEWi2HLli2WAT8aU/jx42P48eP7cPqR9fiftSsQqo333Zip702ZBv3fAviJiPyriBwhIq0i0gHgRwB+k//mUanxwjy+lXXr1tnOmM7H8rVCSR2FyGdA1gJZahnf+vr6sg762jLHUCiE7u7uZEVHo/n4yclJTExMzLk/vexy+u+hlrAXiUQsk+vGp2P4xp/2YuvzE+hYHsLPPrUCC6pez19hpr43ZZqh9EUAAuCniecKgGkA1wO4NL9No1Jkdx7fTuEZt9J6T1rVPqN96ZuamrBu3TrX5gMUahRCr458vnecKzWxWAyhUAhDQ0Po6urC5ORk1r8T2vRSZ2cnpqamsHHjxlnbLp999tnYsWMHGhsbEYlEdC+2nn/tIL7zyAhem4riGx9sw7+cOjc3mpn63pTRnH7ySSLzEa/CJwCeV0rNvSR1EOf0nWNVfa+1tRXA64lAXmBVr79QW+tSeUqtrNff32/4/0hLskv9v6aUwtYXJvCLbftQP8+Hz6+ow0ffrb/winP6pcXunL5lT19EbrFxDABAKfWPtlpHnrRt2zbLY9L/iHjBqlWr8MQTT8zagGdqagrf+ta3ko+vWrUKb3vb2zA5yc0iKXva8L02jWZ24Tw9PY2GhobkxcHe8Sl0943h4Zen0B6qwkWn1qKhWn+aSqvBn/p/lZn63mBneF+/ZiNRGjvD9V4L+EC8p3/zzTfPuX9sbCw5/7pq1SpXDvHnwsnVAaUil9Ub2l4Jd91116zh+1WrViUT6oyG4IH4FJo2EvDEqwdx7aN7MXIgho62avzTsQvgE0FNTc2c5xn18BnwvcEy6CulPlGMhlBps9PL9yqz+vwzMzPYtGkTVq1a5drqftliwLeWbcBvamrCRRddhM7OTsPVMNPT0/D7/YbniEaj2Dd1EL94Yh/uemkSh9dUYP1b63H0IZXJYyKRCKqrq2cFc7OsfQb90uf+UmNUEko1KS8frOrza4+vW7cOl19+eTGaRCUqvZiO1fJXs/93veEp/PjxMUQmYvjAsQvwkTdVI+CfnWCpF8yZte9tDPqUs1LeOjcfrHrwwWAQp556KgvTkC6tyE4mxXTMDI3P4Kfb9qFv8AAOX+jH1951CI5tMM64Tw/mzNr3NgZ9ohytW7cOXV1dhkP8TN4jI5mUybUyNRPD/z2zH7/fvh9+n+DfTlyI9x4zH5U+88JL6cG8ubnZcE6fSh+DPuWk3Hv5wOvr9t26Hp/cSavjkGvAPxhVuPOlCfzv3/Zj5EAMpy8J4t9OWoiGecbz/Rq9YK4N9bO+vjcx6FPWIhEu7CDKRL4CPRAvoXv/zkn899Pj2DMRw/GNlbjktDoc12h/GD51zb22TI+B3tsY9ClrXq2xn41Nmzaxl1/GampqMDMzkyydq22mMzY2lre5es3UTAz37JjEH56dwKv7oziqvgJr22tx0uJAxnsopAb81CH96enpZC0ABn5vYdCnrLCXP5tVBj+VjpqaGkxMTMzK0aioqMCGDRsAvF5+Od/B3MreqSj++PwE/vjCBManFY5tqMS5Jy3Eqc1VWW2YlDqXz2V65YNBn7LCXv5sXluDX66CwSAuueQSAMbBvVhBHgCiSuGJoWnc9dIEesMHEFPAKc1V+KdjF2Q0jJ8ufS6fy/TKB4M+ZWz79u1ON8F1rDbWqampwVlnnYUtW7ZwGqCAfD4fPvjBD+KOO+5I7nWgVQ7Um0/funWrK4J7up2jB/GnXVO4b2ASwxMx1FT5cM4bFuDMI+aheWFuf7b15uu5TK98ZLXhjttxw53CYsa+PrMAondMTU0NpqenTZf0HXnkkRgeHk4GMLOyq04SEXzoQx/CSSedZLq50Lx585Lvd968eQgEAhgbG0sOqR88eDB57Pz58zE5OTnns0qtxk0AAB4JSURBVNy4cSN+97vfIRaLQUQQDAYxNTVV9OH2fFJK4ZV9UTzy8hT+tGsKu8Zm4APQtjiAM4+Yj1MOq7JcemfGaqMcbq5T+uxuuMOgTxnZvn07xsfHnW6Gp9i5WCDvmY4qPL1nGn2DB/DY4AHs3h+vrvfGxkq8dUkQbzk8iPqg9bI7K3Yz8Zm9X9oY9Bn0C4K9fKLsHIwqPPvaQTy9ZxpPvzqNZyPTmI4BAT/QdmgVTg5V4ZRQFRrm5x7o/X4/lixZwqBdRvK2tS6RhnP5RPZElcLu8Sief+0gnkvcXho5iJkYIABa6ypw1lHzceLiKpxwaABV/uyH7jU+nw9Lly5loCdTDPpkG4f1iWaLKYXIZAzhfTN4ZWwGA6Px287RGRyIxkdRq/yCo+or8N6j5+P4RQG8sTGA6oAvr+3w+/1YtmxZXs9J3sSgT7awl0/lSCmF/QcVIpNRDE9EsXt/FK8mbrvHowiPz2A6ZaO76oCgtbYSZx45Dy21FTiqvhJLairgzyEJz45y3uWSMsOgT7awl09eMR1V2H8whvFphbEDMYxPxzB2IH4bmdJuUYxMxRCZjCV77JqADzh0gR+HLqjACYcGcNjCCjQv9OOwhRWoC/qyKpSTKy6tI7sY9MkSk/eoGJRSiCpgJqYwEwMOav9GFWZiCtPR+H3T0dQbcCDx9dSMwoGZ+L9TMzFMzihMHlSYnFGYOBjDxEGF/dMxTJvscDy/UlBX5UNd0IfWukq0N/vQMM8fv8334dAFftRVORPYjXAHPMoEg76F2/oH8bew/prjcjE4tM/pJthSjIUoVi9h1Yb0h41Wz6jk4/YfS79fQaV8nf5Y/GuV0obZ96X+qxBTr9/3+tfx+7XvYypxXOKYaCz+fTTxWFQBsVg8yS2aeDz+b/zrmTz8/AI+oKpCUFUhmF/hw7xKwYJKQeP8SiyoFCwI+JL/VlcKaqp8WFjlQ00g/m8gDwl1haQtpePSOsoWg76Fu//+Km7e9orTzXBUKS3rzMefbMtzWByQ/vCcwyX9W9F/WOYebvSY3v0CACIpX89+TOT1+yXRDu2+9H99En9Mm5r2SeK+xDG+xK3SJ/AlnuP3+ZL3+0RQ4Yt/7Zf4MRU+gd8H+BNfV/gAv0/iXwtQ4RdU+uLnrPAJAn6g0i8I+AWVie+rKgRVifuq/FLwuXOnaQGeQZ6yxaBv4TsfPgnf+fBJTjfDMRzaJ3IHv9/PYE85y++6EfIUZuwTFZdZrgAz9CkfGPTJEDP2iYrL5/MZZuIzQ5/ygUGfdLGXT1R80WgUzc3Nc3r8zNCnfOGcPuliL5+o+AKBQHLenhn6VAgM+jRHJBJxuglEnlBdXY3Gxkbs2rUrOSfv9/tRX1+PSCQyZytbrTfPDH0qFAZ9mmPHjh1ON4GopDQ2NqK6utqwd64XwM2OJyoUBn2a5amnnnK6CUQlQS9QZxK02ZsnJ7g+6IvIBQAuARAC8DSAi5VSDzrbKu86cOCA000gyhufL14yN3W5W2qw3r59u27+SiAQQE1NzZwheCDeq29paSl424kKwdVBX0Q+AuD7AC4A8FDi360icrxSaqejjfMgZuyTE3w+H2KxWDIYA8DAwMCc+W6fz4doNGo4FB6JRDIeLj/22GNNH+cQPHmNq4M+gM8B+LlS6keJ7y8SkbMBnA/gS841y5uYsU/54vf7AZgXlLEKopkG20IMl3MInrzGtUFfRAIA2gF8O+2hOwCcVvwWeRsz9ikbqcE9nz1wBluiwnBt0AfQCMAPYHfa/bsBnJl+sIicB+A8AFi6dGnBG+c1zNinTGQyr80ATuQebg76mvQt3kTnPiilbgBwAwAsX768dLaFc4Ft27Y53QQqAZzTJip9bi7DOwwgCqAp7f5DMbf3TzngRh7e5ff70draali33e/3m27yogkEAmhra2PAJypxrg36SqlpAH0A3p320LsBPFz8FnnTwMCA002gAgkEAli2bBkaGhoM67kvWbIELS0tyYsCvYsA1n0n8g63D+9/F8CNIvIXAH8CsBZAM4BuR1vlIcPDw043gQogPVBb1XNP7cFnk3hHRKXB1UFfKfVfItIA4MuIF+d5CsB7lVLsnuYBq+95k1GgtptQx8Q7Iu9yddAHAKXUdQCuc7odXsTqe94iImhpaWHAJiJDrg/6VBhcl+8tHIYnIjsY9MsU1+WXFr/fj1gsNqc0LXv2RJQJBv0yxIz90qJl2QOZl6YlIkrFoF+GmLHvfoFAwPa+7EREdjHolxnO5bufVgiHiCjfGPTLDOfy3U1bX8+18kRUCAz6ZYS9fPfx+/3w+/2zgjswez/56enpZB5GQ0MDLwiIKGsM+mWEvfzi0ubljWgJeukBu7+/f1aWPgAopRAOhwGYXxAQEZlh0C8TrL5XXK2trckgnNoz15j10I0uFKanpxEOhw0vCBj0icgKg36ZYPU9Y1VVVbY/H21tPPD68jm/3w8gvluhXjDPtKyt0QiB2ciB2YgCEZGGQb8M2FmXLyJzepDlwirgF3v5XHNz86whfOD1BL/0EYPUNhIRWWHQLwN21uWffPLJ2LZtG6LRaBFaVDqcWD5ntSOe0QUBEZEVBn2P27Ztm+1jGfBnczKYGk0JWF0QEBGZYdD3ODuBvL29HYB1tnkpycd0hVvr2nPrWyLKls/pBlDhZFpjv6ampkAtKb5cA77f78eOHTvQ39/P+gZE5Bns6XuYnbn8QCCAvr4+BAKBsh3eFxH4fD5Eo9HkbnbaZ8F18ETkJQz6HmW3d6oN53tlWD8bqcP4/f39cy5+uA6eiLyCw/sexep79qUGc66DJyIvY9D3oHKpviciEJG8ntNovTvXwRORFzDoe5Cd6nKlHsQCgQBaWlpyTthrbGyc9X1zc/OcCwmugycir+Ccvsds377d1nF6Vd9KQWpNewCGFerSiQgWLFiA8fHx5H2NjY3JkroaroMnIi9j0PeY1KBmRFuXD9gPmk7z+/2or69HOBzGjh07ksHY6OKluroa09PTWQVuroMnIq9i0C9jWnDr6+tzuimm2tvbEYlEdLeUbWlpQUtLC3vmREQ2MOh7iJ3gndrLj0Qi2LlzZyGblDMt98BsS9m2tjYGeSIiGxj0PSLTqnHpPWc3Sk2g41I6IqLcMXvfI+ysy0/t5e/atcvVAV/Lztd68FxKR0SUO/b0PcBOxn7qMrRIJOLakrsiorvRjdke80REZA+DvgfYydg/+eSTk1+Hw+FCNseQ1c53Zkl4XEpHRJQ7Bv0yVMx5cC3Qa0HaaIlgIBBAW1ub6bm4lI6IKDec0y9xmWbsA/E178VQXV2NysrKWfex4h0RkXMY9EuYnYz9fNemt6uxsRH79++ftYuftkVtS0tLMgEvPWGPiIgKh8P7JcxOxn7qXL7GThKfiKChoQHDw8MZt8vv92NsbIzr6omIXIY9/RKl9ZrN5DKM39LSkvHaf82SJUu4rp6IyIUY9EuUnR74smXL5txnJ5BXV1djfHzc1ev4iYgoc64N+iJynojcKyIjIqJEpNXpNrmF3cCtZ9euXabPExE0NjZmNayvcWpJIBERmXNt0AcwH8AdAK5wuB2uY2cu/9hjj51zn52iPNrmNbnQ1tHrYQU9IiLnuDboK6WuUUp9A8BDTrfFTezM5ev18rVa+1YaGhos592tVgRoa/K5NI+IyF1cG/RJn51hd71evp1a+42NjQDMe+ONjY1oaWkxfFwL7A0NDVyaR0TkMp5Zsici5wE4DwCWLl3qcGsKw85cvha409lZpqcFc70699q5tax+vZK6fr8fS5YsSQZ2vQp6kUiEpXSJiBxS1J6+iHwtkZRndjsjm3MrpW5QSi1XSi1ftGhRnlvuDnbm8s164Xbp9dJbW1uT59bb2x6IB32zAK5NMaQX7Ml2aSAREWWm2D39awBstjhmZzEaUmoee+wxy2NaW1sNH7Oz2U2q9F56JBJBf3+/6Xy/VS6A3sWCVrCHvX0iosIratBXSg0DyH4tWBmzs2beKHBGIhHT52vz8EZD7wMDA7ZyCawy81mwh4jIWa6d0xeRJgBNAN6QuOt4EakDsFMp9ZpzLSsdqUHciBbcAcyax9eG3sfHx22v2a+pqTF9PBAIGO6wR0REhefm7P21AB4H8KvE97cmvv9Hx1rkkGx20kufPzei1cE3GnrPpEjP2NiY6eNcxkdE5CzX9vSVUleAhXmwfft2y2P0MvaNku2M5GOI3eoc2vQDs/eJiJzh2qBPcePj45bH6GXs2wni2oY8+cqetzNMr7eMj4iIisPNw/tlz04v36jGvtUOeyKCJUuWmFbq0+rwW1Xg047lMD0Rkbuxp+9idnr5RjX2Y7GY6fO06njbtm0znAbQjqmurp4zJA9wmJ6IqNQw6Jcwox641Xx+IBBIBmizSn1mlfVSHyciotLA4X2XspOxf/LJJ+vebzafz2F4IqLyxaDvQnYS68zm2c3m89M3vTE61iongIiISg+DvgvZqbFv1MuPRCK6Q/YigtbW1jlD8kuWLNE9j9H9RERUujinX4KMMvaB+Hy+Hp/PZzovz6Q8IiLvY9B3maeeesryGL2MfY3RfL5Vwh6DPBGR9zHou0gkEsGBAwdMjzHr5ZvlAmiFc7ifPRFR+WLQdxE7c/lGvXyrIjvaLnp6m+oAXH5HRFQOmMjnEkYBO5VejX2N2dp8LWPfbD97IiLyPgZ9l7Czm51ejX2N2dp8rRfP/eyJiMobg74H2JnLT//a6BgiIvIuBn0XsLOxTnt7u+FjZsPzqdX3uJ89EVF5YyKfC1htrGM2lw/YG9pP/ZrZ+0RE5YlB32F2EvjM5vKB+PC8XuDXG7bnmnwiovLF4X0HDQwMWCbwtba2Wp6Hw/ZERGQHe/oOspOxb6dXzmF7IiKyg0Hfxazm8lNx2J6IiKxweN8hjz32mOnjjY2NlnP5REREmWDQd0BqKVwjDPhERJRvDPoOsDOXT0RElG8M+i6UyVw+ERGRXQz6RbZt2zbTx0WEQ/tERFQQDPpFFo1GTR8/+eSTi9QSIiIqNwz6LsKNb4iIqJAY9F2EFfSIiKiQGPSLzO/3694vIiyuQ0REBcWgX2TLli2bE/j9fj/n8omIqOBYhtcBy5Ytc7oJRERUhhj0iyASiXAzHCIichyDfoFFIpFZZXenp6cxMDAAwN4OekRERPnCOf0CC4fDc+rsK6UQDocdahEREZUrVwZ9ETlERH4oIs+IyKSI7BKR60Wk5LrG09PTGd1PRERUKK4M+gCaARwG4IsA2gB0AHg7gN842ahsGBXcYSEeIiIqNlfO6SulngLwwZS7nheRSwBsEZEapdSYQ03LSCQS0S27KyIsxENEREXnyqBvoAbAAQATTjfEjvQEPo3f78eSJUuYxEdEREXn1uH9WUSkDsCVAH6klJoxOOY8EekVkd49e/YUt4E69BL4gHjQZ8AnIiInFDXoi8jXRERZ3M5Ie84CAH8A8Aric/y6lFI3KKWWK6WWL1q0qLBvxAYm8BERkdsUe3j/GgCbLY7ZqX0hItUAbkt8e45SaqpQDcu3QCCgG+CZwEdERE4patBXSg0DGLZzrIgsBLAVgAA4Wyk1Xsi25Vtzc/OcOX0m8BERkZNcmciXCPh3IJ689wEACxLD/ADwmlLK9WPk2rw9y+8SEZFbuDLoA2gH8JbE18+mPfZOAPcVtTVZamhoYJAnIiLXcGXQV0rdh/iwPhEREeVJSSzZIyIiotwx6BMREZUJBn0iIqIywaBPRERUJhj0iYiIygSDPhERUZlg0CciIioTDPpERERlgkGfiIioTIjenu+lTkT2ABhIuasRNjf6KVFef3+A998j31/p8/p75PtztxallOW+8p4M+ulEpFcptdzpdhSK198f4P33yPdX+rz+Hvn+vIHD+0RERGWCQZ+IiKhMlEvQv8HpBhSY198f4P33yPdX+rz+Hvn+PKAs5vSJiIiofHr6REREZY9Bn4iIqEyUfNAXkfNE5F6R/9/emUffNV1x/POtMYhqEFMRLMsYNfOrKVFpQ0eJVZSiFjVVS03L0kZUV1OtElVDVVFRkhrWMg8RghIhLFNIJCuGICTRVvyQid0/9nncvLz5997vvfze/qx113v33H3O3fue+84+59x939H/JJmkfnnH+0n6u6QZkj5JnyMk9SpT7vBUXnZ7t5G2lNClpI1J5iuSRkn6IG2jJK1RplwlO99J12a8pG0aZUclpPrKv+657YwS+QYUybNld+pfKela5+s6uoJ8QyW9LGlB+jywO/StBkl9JF0qaUq6r2ZKukLSmmXyHVWkDlfuLt1LIelESa9Jmi/pGUl7lZHfJ8nNT+3O8d2lazVIOlvS05LmSZoj6U5J25bJU+x3Ori79K6UWtpySf0lPZLu37clDZOk7tK5kSzzTh9YBXgAGF7k+JbAcsAJwDbAycARwCUVlD0VWC+z9e+irrVSzkaAG4Edgf2Bwen7qDLlngmchl+TXYDZwFhJvbuob1eYyZLXfD3gRMCAWyrIv01e3mmNUbMuXMuSuh5XSlhSBzAG+Cewffq8WdJuDdazWtYHNsDvr/7A4cDewE0V5P2YvPo3s/kN0rNiJB2Mtxm/A3YAngDulbRREflNgHuS3A7ACOBSSUO7R+OqGABcDnwd2BdYDDwoqU8FeQezZH091CAdu0rFbbmk1YGxwHt4u/hz4Azgl41Xsxswsx6xATvjjqFfBbInAu+XkRkOvNRsuyqxEdgqpe+RSdszpW1RpCwBs4BzMmm9gA+B45pta56uY4EHysgMSPau1Wx9K7RpPPCXKvOMAcbmpT0I3NRseyrQ/QDgM2D1EjJHAZ3N1rWIbhOBv+WlTQNGFJG/AJiWl3Y1MKHZtlRg62rAp8B3S8j0S7+3nZutbwX2VNWW4wPEeUCvTNqvgLdJwe/L8tYTRvq1sDrw3wrkNk1TO69JGi1p00YrViMdQCc+qsjxOPAR3nsvxCbAuvgMAgBm9gnwaIk83U4aMX2Dyl+nmSRplqRxkgY2ULV6cIikuZImS7qwghmWDjL1lbifFqqvEqwOLMBH8qXoJekNSW9JukvSDt2gW0kkrQjsxNLX/gGKX/tidbWzpBXqq2Hd6Y3PAlfSRt4mabakxyUd1GC9ukI1bXkH8FhqD3Pcj89g9Wukkt1B2zn9NB13Oj6dVYqJ+Mhjf+BY3EE+Ue65ZJNYF5hjqUsKkL7PTseK5QGfwsryXok8zeBY/P+wby8jNwvvoQ8FhuDTeeMk7d1Y9WrmRuAwYCBwPq73bWXyrEvr19dSpNiS8/GR8uISolOBo4HvA4cC84HHJW3eeC1Lshb+iLCaa1+srpZP5bUylwDPARNKyHTi7egP8VmcccAYSYc3Xr2qqbYtL1Z3uWPLNMs3W4FCSPotcE4ZsYFmNr7KctfBe2xjgYtLyZrZvXl5nwRmAEcCF1Vz3iK61NvGQn+4oCLppfJVkqdqarFX0vL4j/U6M1tUKqOZTcWdRo4J8oDH0/HZi4ZTjY1mlp25eFHSDGCipB3N7NkS+bulvgpRYx2uCtyJT42eWSqjmU0g42gkPYE7n5Px56rNptprX0i+UHrLIOki/NHgnmb2aTE5M5sL/CmTNEnSWngd39BYLaujxrZ8mau7SmlJpw+MpPyN82Y1BUpaFw8yeQn4cXZUXAlm1ilpMlCvUUc9bXwX6CtJObtSpOnaLN1jzeYB77nOzKT3LZGnK9Ri73fxoJurazznROCQGvPWQlfqdBL+HHVzoJjTf5elRxqNqq9CVGWfpNXwYDaA71iVAXlm9qmkSdTvN1crc/G6qebaF6urxcD7ddWuTki6GP+9DDSzGTUUMRH4SX21qj8VtOXF6g6677fWMFrS6adeZN2WOJS0HvAwMBk4tMwUY7EyVsbfBHi4HjrV2cYJePBNB1881+8AVmXJ5/xZXsNv7kHA0/C5jXvhkap1pUZ7jwUeMbNXazzt9vi0f7fQxTrtj08hl9J3Al5ff8ykDaJ4HdeVauxL8Qn34iOkwWbWWe35Usd1O+D5avPWEzNbKOkZ/FrfnDk0CLi1SLYJwA/y0gYBk8rNWjUDSZfgDn+AmU2psZhu/b3VSgVt+QTgAkkrZzqqg4B3gNcbr2GDaXYkYVc3vEe2PfAjfOrlgLTfJx1fH3gVj5beMMnntuUy5UwBfpbZvxDYBw942w24C4/o3LjVbEwy9wIvArvjDv9F4M7M8Q2SjQdm0s5KNg0BtgVG4zd27xao143w0dVhRY5fD1yf2T8Fb2Q3x1/bG5Gu1ZBm21JA982AYfjbGP1Sfb6Cj/Cz9+Q4MtHheNDYYuBsvNE6G1gE7NZsm/Ls6403nLnRVPY3t2IJ+84FvgVsmu7va5J9u7aATQcDC4Fj8LdlLsGfa29c5H7cBA+kHZnkj0n5hzbblgK2XZbagX3z6mq1jMwIYFxm/8jUHm0FbIE/RlsInNpsewrYV7ItL2Dbl/EB0ejULg5J8qc125a6XI9mK1CHCh2eGvf87ah0/Kgix5d49S3tD8/s5xzgQvx55K3A1q1oY5Lpg0+9zkvbDcAameP9CuRRKnsWHjT1CLBts+s06XYe8B9g5SLHxwPjM/tnAtOBT1K+x4ADmm1HEd03TNf6fTyifTruRPrkyb2OxzNk0w7CO28L8Y5CK3ZqBpT4zQ0oZh8eZ/NGuiaz8fibjmbbk9HvxKTzAuAZYO9i92NK2wfvyC3AZ9aOb7YNRewqVlfDMzLXAa9n9o8EXsY7NvPwx1OHN9uWIvaVbMvzbUtp/fFYoPmpfTyXHvC6npnFgjtBEARB0C603St7QRAEQdCuhNMPgiAIgjYhnH4QBEEQtAnh9IMgCIKgTQinHwRBEARtQjj9IAiCIGgTwukHwTKOpOsk3dWgcofVu9xGIKmvpDmSvtpsXYKglQmnHwTBUkjqj692N7LZulSCmc3G/xXvvGbrEgStTDj9IAgKcTJwq5nNa7YiVXAtcJikPs1WJAhalXD6QdCDkLSSpJGS3pM0X9KTkvbMk/m2pKnp+KOSDpFkaSliJC2Hr5N+R16+1yUNS9P+H0qaKelgSWtIGi2pU9I0Sd/My7e1pLtTntmSbkqrXuaO7yLpAUlzJc2T9G9JHXllmKSfSrpZ0keSZuSv3W5mL+F/tzqk61cyCHom4fSDoGfxB3xxmKOBHfCFl+5LK00iaSPgNuBu4GvAn1OeLNvhi45MKlD+KcBTwI7Av4B/ADfiS+huj/9f+Q1pJbPcCpeP4kta7wrsh68IeYekXPvTGxiFr/C4K/AccE9anz3LMOD2pPcY4BpJG+fJPIX/530QBAUIpx8EPQRJqwInAGeZ2d1m9gpwPL4G+ElJ7ARgBr5i2FQzuwW4Mq+ojfEFV94tcJr7zexyM5uGL0KyEjDdzK43s+nA+cDa+OpkufM9b2ZnmdkrZvYCcASwC77KIGb2kJmNSsen4I8W5gOD8849ysxuSOf5Nb7i4F55Mu/gi0sFQVCAcPpB0HPYDFgBeDyXYGaf4svcbp2StgSetiVX2pqYV04vYJGZfVbgHC9kyu4EPsZnE3K8lz77ps+dgL3T1H+npE5gZkbfXOT9XyW9KukD4MOUf6MS514MzMmcJ8cnSf8gCAqwfLMVCIKgbih9Flo60zIy5ZbWnAusKGkVM/s479iiAuUuytuHLwYUX8IfJZxe4Dy5DsI/gHWAU/li6dpxwIoVnDt/4NIH7wwEQVCAGOkHQc9hOr5m+OeBeykorwNf+xzgFXxqPcuuefvPpc+t6TrPAtsAb5jZ9LztwySzJ3BpeiQxGR/pr1fj+bZN5wyCoADh9IOgh2BmHwFXAL+XdICkrdL+OsDlSexKYDNJF0raQtIQ4LhcEamcObjjXCLqv0Yuw4MCx0jaTdKmkvaTdJWk3knmVeDwFOW/CzAa77xUhaRV8McJ99VB7yDokYTTD4KexVl4VP21+Ih9O2Cwmc0CMLM3gKHA94Dn8Sn13B/azM+UcxVwWFeVMbN3gD2Az3BnPBnvCCxIG/ibBqsBz+AO/xp8mr9avg+8aWaPdU3rIOi5aMl4niAI2g1JvwB+A3wlF7wnaSVgCnDEsuJEJT0FjDSzG5utSxC0KhHIFwRthqSTgKfxgLfd8dffrstG65vZAklH4oFxLY+kvsAtwE3N1iUIWpkY6QdBmyHpYvwf99YE3sKn1H9jZlU/Rw+CYNkinH4QBEEQtAkRyBcEQRAEbUI4/SAIgiBoE8LpB0EQBEGbEE4/CIIgCNqEcPpBEARB0CaE0w+CIAiCNuH/94nMUZi7HIkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "merge_setup = {'/FullDataset_v1':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, \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"
   ]
  }
 ],
 "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
}