notebooks/scran_preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "empirical-relaxation",
   "metadata": {},
   "source": [
    "# plug-in_scapy_scran"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "greek-smart",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T14:59:57.822022Z",
     "start_time": "2020-04-10T14:59:57.819335Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Thomas Walle\n",
      "\n",
      "2021-05-01 19:12:28.476726\n"
     ]
    }
   ],
   "source": [
    "import datetime\n",
    "\n",
    "print(\"Thomas Walle\\n\")\n",
    "print(datetime.datetime.now())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "intimate-saturn",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T15:00:01.576328Z",
     "start_time": "2020-04-10T15:00:01.570091Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/wallet/miniconda3/envs/2020_run_scran/lib/python3.8'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Which python\n",
    "import sys\n",
    "\n",
    "sys.path[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "young-thesis",
   "metadata": {},
   "source": [
    "## Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "discrete-regular",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:32:52.846019Z",
     "start_time": "2020-04-10T16:32:52.840928Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Packages\n",
    "import scipy\n",
    "from scipy import io\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scanpy as sc\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy import sparse\n",
    "from matplotlib import rcParams\n",
    "\n",
    "%matplotlib inline\n",
    "sns.set_style(\"white\")\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "opened-concern",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T01:04:24.616212Z",
     "start_time": "2020-04-10T01:04:24.613407Z"
    }
   },
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "crazy-miller",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T15:00:08.011647Z",
     "start_time": "2020-04-10T15:00:07.985655Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Functions\n",
    "\n",
    "# Count distribution\n",
    "def total_molecules_per_cell(adat, count=str):\n",
    "\n",
    "    # Plot distribution of total molecules per cell\n",
    "    a = adat.copy()\n",
    "\n",
    "    print(\"Median total count of UMI is:\", np.median(a.obs[count]))\n",
    "\n",
    "    # Per cell (Total molecules per cell)\n",
    "    fig, (ax1, ax2, ax3) = plt.subplots(nrows=3,\n",
    "                                        ncols=1,\n",
    "                                        sharex=True,\n",
    "                                        figsize=(15, 5))\n",
    "\n",
    "    a.obs[count][a.obs[count] != 0].plot(kind=\"hist\",\n",
    "                                         bins=100,\n",
    "                                         ax=ax1,\n",
    "                                         color=\"r\")\n",
    "    a.obs[count][a.obs[count] != 0].plot(kind=\"kde\", ax=ax2, color=\"r\")\n",
    "    a.obs[count][a.obs[count] != 0].plot(kind=\"box\",\n",
    "                                         vert=False,\n",
    "                                         ax=ax3,\n",
    "                                         color=\"r\")\n",
    "\n",
    "    x = ax3.set_xlabel(\"Molecule counts \"+ count +\" per cell\")\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "# Randomize cells for plotting UMAP\n",
    "def randomize_cells(anndata_object):\n",
    "\n",
    "    \"Randomize cells for plotting UMAP\"\n",
    "\n",
    "    index_list = np.arange(anndata_object.shape[0])\n",
    "    np.random.shuffle(index_list)\n",
    "    anndata_object = anndata_object[index_list]\n",
    "\n",
    "    return anndata_object\n",
    "\n",
    "\n",
    "# Percentage of variance explained\n",
    "def pcs_variance(annadata_object, var_explained):\n",
    "\n",
    "    for idx, var in enumerate(annadata_object.uns[\"pca\"][\"variance_ratio\"]):\n",
    "\n",
    "        var_exp = np.sum(annadata_object.uns[\"pca\"][\"variance_ratio\"][:idx])\n",
    "\n",
    "        if var_exp > var_explained:\n",
    "\n",
    "            return idx\n",
    "\n",
    "\n",
    "# Plot function\n",
    "def plot_cumsum(pc_var, pc, cumvar, color=\"r\"):\n",
    "    \n",
    "    plt.figure()\n",
    "    plt.plot(np.cumsum(pc_var), c=\"k\", zorder=0)\n",
    "    plt.scatter(pc, cumvar, c=color, zorder=1)\n",
    "    plt.xlabel('Number of Components')\n",
    "    plt.ylabel('Variance (%)')  #for each component\n",
    "    plt.title('Dataset Explained Variance')\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "# kneepoint\n",
    "def kneepoint(vec, color=\"r\"):\n",
    "    curve = [1 - x for x in vec]\n",
    "    nPoints = len(curve)\n",
    "    allCoord = np.vstack((range(nPoints), curve)).T\n",
    "    np.array([range(nPoints), curve])\n",
    "    firstPoint = allCoord[0]\n",
    "    lineVec = allCoord[-1] - allCoord[0]\n",
    "    lineVecNorm = lineVec / np.sqrt(np.sum(lineVec**2))\n",
    "    vecFromFirst = allCoord - firstPoint\n",
    "    scalarProduct = np.sum(vecFromFirst *\n",
    "                           numpy.matlib.repmat(lineVecNorm, nPoints, 1),\n",
    "                           axis=1)\n",
    "    vecFromFirstParallel = np.outer(scalarProduct, lineVecNorm)\n",
    "    vecToLine = vecFromFirst - vecFromFirstParallel\n",
    "    distToLine = np.sqrt(np.sum(vecToLine**2, axis=1))\n",
    "    idxOfBestPoint = np.argmax(distToLine)\n",
    "\n",
    "    cumvar = np.cumsum(vec[:idxOfBestPoint + 1])[-1]\n",
    "\n",
    "    print(\"Kneepoint happens at PC:\")\n",
    "\n",
    "    print(\"PC\", idxOfBestPoint)\n",
    "\n",
    "    print(\"The cumulativate variation explaind is:\")\n",
    "\n",
    "    print(cumvar)\n",
    "\n",
    "    plot_cumsum(vec, idxOfBestPoint, cumvar, color)\n",
    "\n",
    "    return idxOfBestPoint\n",
    "\n",
    "\n",
    "# Treatment colors\n",
    "treatment_colors = [\"red\", \"green\", \"yellow\", \"blue\"]\n",
    "\n",
    "\n",
    "# Select PCs\n",
    "def select_pcs(adata_obj, color=\"r\"):\n",
    "\n",
    "    # Slope function\n",
    "    def slope(x1, y1, x2, y2):\n",
    "        return (y2 - y1) / (x2 - x1)\n",
    "\n",
    "    pc_v = adata_obj.uns[\"pca\"][\"variance_ratio\"]\n",
    "\n",
    "    max_slope = slope(1, np.cumsum(pc_v[0])[-1], 2, np.cumsum(pc_v[:2])[-1])\n",
    "\n",
    "    for i, v in enumerate(pc_v):\n",
    "\n",
    "        x1 = i\n",
    "        x2 = x1 + 1\n",
    "\n",
    "        if x1 == 0:\n",
    "            y1 = np.cumsum(pc_v[x1])[-1]\n",
    "        else:\n",
    "            y1 = np.cumsum(pc_v[:x1 + 1])[-1]\n",
    "\n",
    "        y2 = np.cumsum(pc_v[:x2 + 1])[-1]\n",
    "\n",
    "        s = slope(x1, y1, x2, y2)\n",
    "\n",
    "        if s <= max_slope / 10:\n",
    "\n",
    "            print(\"Change of slope no less than 1/10 of PC1 to PC2 is in:\")\n",
    "\n",
    "            print(\"PC\", x1)\n",
    "\n",
    "            print(\"The cumulativate variation explaind is:\")\n",
    "\n",
    "            print(y1)\n",
    "\n",
    "            break\n",
    "\n",
    "    plot_cumsum(pc_v, x1, y1, color)\n",
    "\n",
    "\n",
    "# Genes to potentially remove\n",
    "def genes_to_keep_remove(adata_obj):\n",
    "\n",
    "    # Ribosomal genes\n",
    "    riboL = [gene for gene in adata_obj.var.index if \"RPL\" == gene[0:3]]\n",
    "    riboS = [gene for gene in adata_obj.var.index if \"RPS\" == gene[0:3]]\n",
    "\n",
    "    ribo = list(set(riboL + riboS))\n",
    "\n",
    "    ribosomal_genes_df = pd.read_csv(\n",
    "        \"/home/ec2-user/software/gene_sets/qc/ribosomal_genes.csv\", sep=\",\")\n",
    "\n",
    "    ribos = ribo + list(ribosomal_genes_df[\"ribosomal_genes\"])\n",
    "\n",
    "    # Malat1\n",
    "    malat1 = [\"MALAT1\"]\n",
    "\n",
    "    # Mitochondrial genes\n",
    "    mitos = [gene for gene in adata_obj.var.index if \"MT-\" == gene[0:3]]\n",
    "    mito_ribos = [gene for gene in adata_obj.var.index if \"MRPL\" == gene[0:3]]\n",
    "\n",
    "    # Actin genes\n",
    "    actns = [gene for gene in adata_obj.var.index if \"ACTN\" == gene[0:4]]\n",
    "\n",
    "    remove_genes = {\n",
    "        \"malat1\": malat1,\n",
    "        \"mitos\": mitos,\n",
    "        \"mito_ribos\": mito_ribos,\n",
    "        \"actns\": actns,\n",
    "        \"ribos\": ribos\n",
    "    }\n",
    "\n",
    "    genes_to_remove = set()\n",
    "    for geneset in remove_genes.values():\n",
    "\n",
    "        genes_to_remove = genes_to_remove | set(geneset)\n",
    "\n",
    "    keep_genes = list(set(adata_obj.var.index) - genes_to_remove)\n",
    "\n",
    "    return (keep_genes, remove_genes)\n",
    "\n",
    "\n",
    "# 64 color palette\n",
    "godsnot_64 = [\n",
    "    # \"#000000\",  # remove the black, as often, we have black colored annotation\n",
    "    \"#FFFF00\",\n",
    "    \"#1CE6FF\",\n",
    "    \"#FF34FF\",\n",
    "    \"#FF4A46\",\n",
    "    \"#008941\",\n",
    "    \"#006FA6\",\n",
    "    \"#A30059\",\n",
    "    \"#FFDBE5\",\n",
    "    \"#7A4900\",\n",
    "    \"#0000A6\",\n",
    "    \"#63FFAC\",\n",
    "    \"#B79762\",\n",
    "    \"#004D43\",\n",
    "    \"#8FB0FF\",\n",
    "    \"#997D87\",\n",
    "    \"#5A0007\",\n",
    "    \"#809693\",\n",
    "    \"#6A3A4C\",\n",
    "    \"#1B4400\",\n",
    "    \"#4FC601\",\n",
    "    \"#3B5DFF\",\n",
    "    \"#4A3B53\",\n",
    "    \"#FF2F80\",\n",
    "    \"#61615A\",\n",
    "    \"#BA0900\",\n",
    "    \"#6B7900\",\n",
    "    \"#00C2A0\",\n",
    "    \"#FFAA92\",\n",
    "    \"#FF90C9\",\n",
    "    \"#B903AA\",\n",
    "    \"#D16100\",\n",
    "    \"#DDEFFF\",\n",
    "    \"#000035\",\n",
    "    \"#7B4F4B\",\n",
    "    \"#A1C299\",\n",
    "    \"#300018\",\n",
    "    \"#0AA6D8\",\n",
    "    \"#013349\",\n",
    "    \"#00846F\",\n",
    "    \"#372101\",\n",
    "    \"#FFB500\",\n",
    "    \"#C2FFED\",\n",
    "    \"#A079BF\",\n",
    "    \"#CC0744\",\n",
    "    \"#C0B9B2\",\n",
    "    \"#C2FF99\",\n",
    "    \"#001E09\",\n",
    "    \"#00489C\",\n",
    "    \"#6F0062\",\n",
    "    \"#0CBD66\",\n",
    "    \"#EEC3FF\",\n",
    "    \"#456D75\",\n",
    "    \"#B77B68\",\n",
    "    \"#7A87A1\",\n",
    "    \"#788D66\",\n",
    "    \"#885578\",\n",
    "    \"#FAD09F\",\n",
    "    \"#FF8A9A\",\n",
    "    \"#D157A0\",\n",
    "    \"#BEC459\",\n",
    "    \"#456648\",\n",
    "    \"#0086ED\",\n",
    "    \"#886F4C\",\n",
    "    \"#34362D\",\n",
    "    \"#B4A8BD\",\n",
    "    \"#00A6AA\",\n",
    "    \"#452C2C\",\n",
    "    \"#636375\",\n",
    "    \"#A3C8C9\",\n",
    "    \"#FF913F\",\n",
    "    \"#938A81\",\n",
    "    \"#575329\",\n",
    "    \"#00FECF\",\n",
    "    \"#B05B6F\",\n",
    "    \"#8CD0FF\",\n",
    "    \"#3B9700\",\n",
    "    \"#04F757\",\n",
    "    \"#C8A1A1\",\n",
    "    \"#1E6E00\",\n",
    "    \"#7900D7\",\n",
    "    \"#A77500\",\n",
    "    \"#6367A9\",\n",
    "    \"#A05837\",\n",
    "    \"#6B002C\",\n",
    "    \"#772600\",\n",
    "    \"#D790FF\",\n",
    "    \"#9B9700\",\n",
    "    \"#549E79\",\n",
    "    \"#FFF69F\",\n",
    "    \"#201625\",\n",
    "    \"#72418F\",\n",
    "    \"#BC23FF\",\n",
    "    \"#99ADC0\",\n",
    "    \"#3A2465\",\n",
    "    \"#922329\",\n",
    "    \"#5B4534\",\n",
    "    \"#FDE8DC\",\n",
    "    \"#404E55\",\n",
    "    \"#0089A3\",\n",
    "    \"#CB7E98\",\n",
    "    \"#A4E804\",\n",
    "    \"#324E72\"\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "integral-worship",
   "metadata": {},
   "source": [
    "## Load data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "liked-exhaust",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T15:00:15.167659Z",
     "start_time": "2020-04-10T15:00:14.830097Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(109655, 22639)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Read raw adata\n",
    "PATH = #insert path to data here\n",
    "FILE = #insert path to data here\n",
    "adata = sc.read_h5ad(PATH + FILE)\n",
    "samples_to_exclude = []\n",
    "adata.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "exceptional-youth",
   "metadata": {},
   "source": [
    "## Dealing with negative size factors\n",
    "\n",
    "It is possible for the deconvolution algorithm to yield negative estimates for the size factors. These values are obviously nonsensical and `computeSumFactors` will raise a warning if they are encountered. Negative estimates are mostly commonly generated from low quality cells with few expressed features, such that most genes still have zero counts even after pooling. They may also occur if insufficient filtering of low-abundance genes was performed.\n",
    "\n",
    "To avoid negative size factors, the best solution is to increase the stringency of the filtering.\n",
    "\n",
    "If only a few negative size factors are present, they are likely to correspond to a few low-quality cells with few expressed features. Such cells are difficult to normalize reliably under any approach, and can be removed by increasing the stringency of the quality control.\n",
    "\n",
    "If many negative size factors are present, it is probably due to insufficient filtering of low-abundance genes. This results in many zero counts and pooled size factors of zero, and can be fixed by filtering out more genes with a higher `min.mean` - see “Gene selection” below.\n",
    "\n",
    "Another approach is to increase in the number of `sizes` to improve the precision of the estimates. This reduces the chance of obtaining negative size factors due to estimation error, for cells where the true size factors are very small.\n",
    "\n",
    "As a last resort, `positive=TRUE` is set by default, which uses `cleanSizeFactors` to coerce any negative estimates to positive values. This ensures that, at the very least, downstream analysis is possible even if the size factors for affected cells are not accurate. Users can skip this step by setting `positive=FALSE` to perform their own diagnostics or coercions.\n",
    "\n",
    "REF: https://rdrr.io/bioc/scran_final/man/computeSumFactors.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "public-remains",
   "metadata": {},
   "source": [
    "`scTranform` filters out genes that are expressed in no more than 5 cells by default. However, when applying `scran`, negative factors were present in data if the number of cells was limited to 5. By limiting the cells to 10, no negative factors were found. Thus, for consistency I remove genes that are expressed in no more than 10 cells. Similarly, `scanpy` tutorial suggests 3 cells as a threshold as a basic filter, but for the reasons mentioned above, I think 10 is a reasonable balance.\n",
    "\n",
    "REF: https://icb-scanpy-tutorials.readthedocs-hosted.com/en/latest/pbmc3k.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sexual-mexico",
   "metadata": {},
   "source": [
    "#define samples to normalize (if necessary)\n",
    "a_length=len(adata)\n",
    "print('length of original adata is', a_length)\n",
    "\n",
    "adata = adata[~adata.obs['sample_id'].isin(samples_to_exclude)]\n",
    "b_length=len(adata)\n",
    "print('length of filtered adata is', b_length)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "christian-investment",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T15:38:54.146933Z",
     "start_time": "2020-04-10T15:38:54.024303Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(109655, 19984)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "### Remove genes expressed in less tha 10 cells\n",
    "sc.pp.filter_genes(adata, min_cells=10)\n",
    "\n",
    "adata.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "pacific-netscape",
   "metadata": {},
   "outputs": [],
   "source": [
    "#add raw .X\n",
    "adata1 = adata.raw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "stainless-slovenia",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AnnData object with n_obs × n_vars = 109655 × 19984\n",
       "    obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'batch', 'cohort', 'total_molecules', 'ID', 'age_rounded_up', 'postmenopause', 'TNM_version', 'T_prefix', 'T', 'T_suffix', 'N_prefix', 'N', 'N_suffix', 'M_prefix', 'M', 'M_suffix', 'HER2', 'ER', 'PR', 'CTX', 'KI67_precent_range_prePD1', 'KI67_precent_range_postPD1', 'TCR_clonotype_expansion', 'Bassez_nCount_RNA', 'Bassez_nFeature_RNA', 'Bassez_patient_id', 'Bassez_timepoint', 'Bassez_expansion', 'Bassez_BC_type', 'Bassez_cellType', 'Bassez_cellSubType', 'n_genes', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_100_genes', 'pct_counts_in_top_200_genes', 'pct_counts_in_top_500_genes', 'percent_mito', 'n_counts', 'pheno_k10', 'pheno_k20', 'pheno_k30', 'pheno_k40', 'pheno_k50', 'pheno_k60', 'pheno_k70', 'pheno_k80', 'pheno_k90', 'pheno_k100', 'doublets', 'major_subsets'\n",
       "    var: 'features-cohort_1', 'features-cohort_2', 'n_cells', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts'\n",
       "    uns: 'doublets_colors', 'major_subsets_colors', 'neighbors', 'pca', 'umap'\n",
       "    obsm: 'X_pca', 'X_tsne', 'X_umap'\n",
       "    varm: 'PCs'\n",
       "    layers: 'imputed'\n",
       "    obsp: 'connectivities', 'distances'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "shaped-asthma",
   "metadata": {},
   "outputs": [],
   "source": [
    "#filter genes\n",
    "adata1 = adata1[adata.obs_names,adata.var_names]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "separate-least",
   "metadata": {},
   "outputs": [],
   "source": [
    "#add .X to adata\n",
    "adata1.obs = adata.obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "engaged-organization",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.X = adata1.X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "vulnerable-ghost",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  (0, 1500)\t1.0\n",
      "  (0, 5284)\t1.0\n",
      "  (0, 9558)\t1.0\n",
      "  (0, 18297)\t3.0\n",
      "  (0, 13119)\t1.0\n",
      "  (0, 13705)\t1.0\n",
      "  (0, 1813)\t1.0\n",
      "  (0, 8688)\t1.0\n",
      "  (0, 10062)\t1.0\n",
      "  (0, 18252)\t1.0\n",
      "  (0, 5552)\t1.0\n",
      "  (0, 16341)\t1.0\n",
      "  (0, 1682)\t1.0\n",
      "  (0, 13585)\t1.0\n",
      "  (0, 9005)\t1.0\n",
      "  (0, 19776)\t1.0\n",
      "  (0, 865)\t1.0\n",
      "  (0, 14846)\t1.0\n",
      "  (0, 3711)\t1.0\n",
      "  (0, 12114)\t1.0\n",
      "  (0, 19706)\t1.0\n",
      "  (0, 9011)\t1.0\n",
      "  (0, 5855)\t1.0\n",
      "  (0, 3393)\t1.0\n",
      "  (0, 12798)\t1.0\n",
      "  :\t:\n",
      "  (109654, 17431)\t1.0\n",
      "  (109654, 8494)\t1.0\n",
      "  (109654, 16589)\t2.0\n",
      "  (109654, 3088)\t3.0\n",
      "  (109654, 19775)\t1.0\n",
      "  (109654, 1367)\t2.0\n",
      "  (109654, 5491)\t1.0\n",
      "  (109654, 9242)\t1.0\n",
      "  (109654, 5079)\t1.0\n",
      "  (109654, 12899)\t1.0\n",
      "  (109654, 7057)\t1.0\n",
      "  (109654, 3435)\t5.0\n",
      "  (109654, 8660)\t5.0\n",
      "  (109654, 12001)\t1.0\n",
      "  (109654, 5967)\t1.0\n",
      "  (109654, 3150)\t1.0\n",
      "  (109654, 17819)\t1.0\n",
      "  (109654, 11514)\t1.0\n",
      "  (109654, 1335)\t1.0\n",
      "  (109654, 1195)\t1.0\n",
      "  (109654, 13716)\t1.0\n",
      "  (109654, 9336)\t1.0\n",
      "  (109654, 19047)\t1.0\n",
      "  (109654, 8489)\t1.0\n",
      "  (109654, 12324)\t1.0\n"
     ]
    }
   ],
   "source": [
    "print(adata.X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "extreme-ranch",
   "metadata": {},
   "source": [
    "## Prepare scran input files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "delayed-columbus",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:30:17.244660Z",
     "start_time": "2020-04-10T16:30:17.217220Z"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate total molecules per cell\n",
    "adata1 = adata.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "prostate-seeking",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:30:26.708552Z",
     "start_time": "2020-04-10T16:30:26.488598Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "### QC metrics\n",
    "sc.pp.calculate_qc_metrics(\n",
    "    adata1,\n",
    "    expr_type=\"counts\",\n",
    "    var_type=\"genes\",\n",
    "    inplace=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "tribal-float",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:30:43.711855Z",
     "start_time": "2020-04-10T16:30:42.856552Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Median total count of UMI is: 2072.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5cAAAE+CAYAAAAZGeTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deVzVVeL/8fdlExIVsQAHEWzSybXl625pYWAFKOIybY6Z1WSiOZqVOVo5tmurk+bUaIv6KxO1pLLEXGjM0haltLRiAtKrCW4oq+f3xx1uIIvg3UBez8fjPu79nM/5nHM+93Kv9+35fD7XYowxAgAAAADAAV6eHgAAAAAAoOEjXAIAAAAAHEa4BAAAAAA4jHAJAAAAAHAY4RIAAAAA4DDCJQAAAADAYT6eHkBD0qtXL4WHh3t6GAAAAADgETk5Odq6dWuV6wiXdRAeHq6UlBRPDwMAAAAAPCIpKanadRwWCwAAAABwGOESAAAAAOAwwiWAc0NBQc3LAAAAcCnOuQRwbvD3lyyW35eN8dxYAAAAGiFmLgEAAAAADiNcAgAAAAAcRrgEAAAAADiMcAkAAAAAcBjhEgAAAADgMMIlAAAAAMBhhEsAAAAAgMMIlwAAAAAAhxEuAZybCgpqXgYAAIBTuSVc/vDDD+7oBgB+5+8vWSy/3/z9PT0iAACAc5qPOzp56KGHVFxcrKFDhyohIUHNmzd3R7cAAAAAADdxS7hctmyZMjMztWLFCg0bNkzdunVTUlKS+vXr547uAQAAAAAu5pZwKUlRUVGaNGmSunTpotmzZ+u7776TMUaTJ09WbGysu4YBAAAAAHABt4TL3bt3KyUlRRs3blTfvn21YMECde7cWVarVTfccAPhEgAAAAAaOLeEy3/84x8aMWKEJk+eLP9yF9UIDQ3VPffc444hADgXFRRwoR4AAIB6wi3hcuHChfL395e3t7ck6dSpUyosLFRAQIASExPdMQQA56KyK8JKkjGeHQsAAEAj55afIhkzZowKyv3G3MmTJzVmzBh3dA0AAAAAcAO3hMvCwkI1bdrUvty0aVOdPHmyxm327dunUaNG6brrrlNcXJxee+01SdLhw4c1ZswYxcbGasyYMTpy5IgkyRij2bNnKyYmRgkJCfr222/tba1cuVKxsbGKjY3VypUr7eUZGRlKSEhQTEyMZs+eLcPMBwAAAACcFbeEy4CAgAphLyMjo8K5l1Xx9vbWAw88oA8++EBvvfWWli5dqr1792rhwoXq06ePPvroI/Xp00cLFy6UJG3atEmZmZn66KOP9I9//EMPP/ywJFsYnTdvnt5++20tX75c8+bNswfShx9+WLNmzdJHH32kzMxMbdq0yTVPAAAAAACc49xyzuWDDz6oe+65RyEhIZKkgwcP6tlnn61xm5CQEHv9wMBAXXjhhbJarUpLS9Mbb7whSUpMTNSoUaM0depUpaWlKTExURaLRZdeeqmOHj2qAwcO6PPPP1e/fv0UFBQkSerXr582b96snj176vjx47rsssvsbaWlpWnAgAGuehoAAAAA4JzllnDZrVs3ffDBB/r5559ljNGFF14oX1/fWm+fnZ2tXbt26ZJLLtGhQ4fsoTMkJES5ubmSJKvVqrCwMPs2YWFhslqtlcpDQ0OrLC+rDwAAAACoO7eES0nauXOncnJyVFpaql27dklSra4Um5+fr4kTJ+rBBx9UYGBgtfWqOl/SYrHUuRwAAAAAUHduCZdTp05VVlaWLr74YvvPkVgsljOGy+LiYk2cOFEJCQmKjY2VJLVq1UoHDhxQSEiIDhw4oODgYEm2mcf9+/fbt92/f79CQkIUFhamzz//3F5utVrVs2fPausDAAAAAOrOLeEyIyND77//fp1mBo0xmj59ui688MIKP1sSHR2tVatW6c4779SqVas0cOBAe/mbb76puLg4ffPNN2rWrJlCQkJ0xRVX6JlnnrFfxCc9PV2TJ09WUFCQmjZtqq+//lqXXHKJVq1apVGjRjl3xwEAAACgkXBLuGzfvr0OHjxYp5nB7du3a/Xq1erQoYOGDBkiSZo8ebLuvPNOTZo0Se+8845at26t559/XpI0YMAAbdy4UTExMQoICNBjjz0mSQoKCtLdd9+t4cOHS5LGjx9vv7jPww8/rGnTpqmgoED9+/dX//79nbnbAAAAANBoWIwbftxx1KhR2r17t7p161bhQj4LFixwdddOlZSUpJSUFE8PA0B5ZUdEGPP74+qWAQAA4JCaMpFbZi4nTJjgjm4AAAAAAB7ilnDZs2dP5eTk6L///a/69u2rkydPqrS01B1dAwAAAADcwMsdnbz99tuaOHGiZs6cKcl2xdbx48e7o2sAAAAAgBu4JVwuWbJEy5Yts/9OZVRUlHJzc93RNQAAAADADdwSLv38/OTn52dfLikpcUe3AAAAAAA3ccs5lz169NCCBQtUUFCgTz/9VEuXLlV0dLQ7ugYAAAAAuIFbZi7vvfdeBQcHq0OHDnrrrbc0YMAATZo0yR1dAwAAAADcwC0zl15eXho5cqRGjhzpju4AAAAAAG7mlnAZHR0tS/kfM/+ftLQ0d3QPAAAAAHAxt4TLFStW2B8XFRXpgw8+0JEjR9zRNQAAAADADdxyzmXLli3tt9DQUN1666367LPP3NE1AAAAAMAN3DJz+e2339ofnzp1ShkZGcrPz3dH1wAAAAAAN3BLuHziiSd+79DHR+Hh4Xruuefc0TUAAAAAwA3cEi7feOMNd3QDAAAAAPAQt4TLRYsW1bh+zJgx7hgGAAAAAMBF3BIuMzIytHPnTkVHR0uSPvnkE3Xv3l2tW7d2R/cAAAAAABdzS7jMy8tTSkqKAgMDJUnJycm655579Oijj7qjewAAAACAi7nlp0h+/fVX+fn52Zf9/PyUk5Pjjq4BAAAAAG7glpnLIUOGaPjw4YqJiZHFYtHHH3+sxMREd3QNADYFBZK/f/XLAAAAcIhbwuW4cePUv39/bdu2TZL0+OOPq1OnTu7oGgBs/P0li+X3ZWM8NxYAAIBzkFsOi5WkkydPKjAwUKNHj1ZYWJiysrLc1TUAAAAAwMXcEi7nzZunV155RQsXLpQkFRcXa+rUqe7oGgAAAADgBm4Jlx9//LHmz5+vgIAASVJoaKjy8/Pd0TUAAAAAwA3cEi59fX1lsVhk+d/5TidOnHBHtwAAAAAAN3HLBX2uu+46zZw5U0ePHtXbb7+tFStWaOTIke7oGgAAAADgBm4Jl2PHjtWnn36qpk2b6ueff9bEiRPVr18/d3QNAAAAAHADl4fL0tJSjR07VosXLyZQAgAAAMA5yuXnXHp7e8vf31/Hjh1zdVcAAAAAAA9xy2GxTZo0UUJCgvr27avzzjvPXv73v//dHd0DAAAAAFzMLeHyqquu0lVXXeWOrgAAAAAAHuDScPnrr7/qD3/4g4YOHerKbgAAAAAAHubScy7Hjx9vfzxhwgRXdgUAAAAA8CCXhktjjP1xVlaWK7sCAAAAAHiQS8OlxWKp8jEAAAAA4Nzi0nMud+/ercsvv1zGGBUWFuryyy+XZJvRtFgs+vLLL13ZPQAAAADATVwaLnft2uXK5gHg7BUUSP7+1S8DAACgTtzyUyQAUO/4+0vlD9cvd444AAAA6s6l51wCAAAAABoHwiUAAAAAwGGESwAAAACAwwiXAAAAAACHES4BQLJdLbaqxwAAAKgVrhYLAFLFq8dy5VgAAIA6Y+YSAE53+swlM5kAAABn1OjD5aZNmzRo0CDFxMRo4cKFnh4OgPqgbBaz7HY6wiYAAEAljTpclpaWatasWXrllVeUmpqqNWvWaO/evZ4eFoD65kxh8+TJisuETwAA0Ag16nC5Y8cORUZGKiIiQn5+foqLi1NaWpqnhwWgvjs9bAYEVB8+Tw+eZ1ommAIAgAaqUYdLq9WqsLAw+3JoaKisVqsHRwTgnFA+fJ4ePM+0fLry4bOuQbWmZUIvAABwMosxjfeyiB988IHS09P16KOPSpJWrVqlnTt3asaMGVXW79Wrl8LDw905RAAAAACoN3JycrR169Yq1zXqnyIJCwvT/v377ctWq1UhISHV1q/uSQQAAACAxq5RHxbbtWtXZWZmKisrS0VFRUpNTVV0dLSnhwUAAAAADU6jnrn08fHRzJkzdfvtt6u0tFTDhg1T+/btPT0sAAAAAGhwGvU5lwAAAAAA52jUh8UCAAAAAJyDcAkAAAAAcBjhEgAAAADgMMIlAAAAAMBhhEsAAAAAgMMIlwAAAAAAhxEuAQAAAAAOI1wCAAAAABxGuAQAAAAAOIxwCQAAAABwGOESAAAAAOAwH08PoCHp1auXwsPDPT0MAAAAAPCInJwcbd26tcp1hMs6CA8PV0pKiqeHAQAAAAAekZSUVO06DosFcG4pKpJKSjw9CgAAgEaHcAng3HHkiNSpk9SqlbRhg6dHAwAA0KgQLgGcOx56SPrxR+noUWnaNMkYT48IAACg0SBcAjg3FBVJb7wh3XCD9MIL0mefSd9+6+lRAQAANBqESwDnho8+knJzpVtukRISbGUbN3p2TAAAAI0I4RLAueGjj6TzzpNiYqTISKlNG2nTJk+PCgAAoNEgXAI4N2zYIPXrJ/n5SRaLdOWVUnq6p0cFAADQaBAuATR8hw5JO3dKV131e9lll0m//mo7VBYAAAAuR7gE0PBt3267793797LOnW33XNQHAADALQiXABq+b76x3V966e9lhEsAAAC3IlwCaPi+/lqKiJCCg38va9tWCgyUvvvOc+MCAABoRAiXABq+r7+uOGsp2S7qc/HF0vffe2ZMAAAAjQzhEkDDVlIi7dkjdepUeV27dlJmptuHBAAA0BgRLgE0bFlZUnGx1L595XVRUbZweeqUu0cFAADQ6BAuATRse/fa7i+6qPK6du2koiJp/373jgkAAKARIlwCaNh+/NF2/8c/Vl4XFWW7//lntw0HAACgsSJcAmjY9u6V/P2lP/yh8rp27Wz3nHcJAADgcoRLAA3bTz/ZQqRXFR9nkZG2e2YuAQAAXI5wCaBhy862/aZlVQICpLAwZi4BAADcoF6Fy02bNmnQoEGKiYnRwoULK60vKirSpEmTFBMToxEjRig7O9u+7uWXX1ZMTIwGDRqkzZs3n7HNKVOmaNCgQYqPj9e0adNUXFzs2p0D4BrZ2VJ4ePXro6KYuQQAAHCDehMuS0tLNWvWLL3yyitKTU3VmjVrtLfsKpD/s3z5cjVv3lwff/yxbr31Vs2ZM0eStHfvXqWmpio1NVWvvPKKHnnkEZWWltbY5uDBg/Xhhx/qvffeU2FhoZYvX+72fQbgoJISyWqtOVzyW5cAAABu4bJwOWHCBG3YsEGnavn7cjt27FBkZKQiIiLk5+enuLg4paWlVaizfv16DR06VJI0aNAgbdmyRcYYpaWlKS4uTn5+foqIiFBkZKR27NhRY5sDBgyQxWKRxWJRt27dZLVanfsEAHC9/fttv2F5ppnLX36RSkvdNiwAAIDGyGXh8sYbb9R7772n2NhYzZkzRz+W/VxANaxWq8LCwuzLoaGhlQKf1WpV69atJUk+Pj5q1qyZ8vLyqt22Nm0WFxdr9erVuvLKK896XwF4SE6O7b5Nm+rrtGtnm+EsqwsAAACX8HFVw3379lXfvn117NgxrVmzRrfddptat26tESNGaPDgwfL19a1Q3xhTqQ2LxVKrOtWVVzVrenqbjzzyiLp3767u3bvXar8A1CNlgfFMM5eS7bzL6i78AwAAAIe59JzLvLw8paSkaPny5erYsaP+8pe/6LvvvtNtt91WqW5YWJj2799vX7ZarQoJCalUZ9++fZKkkpISHTt2TEFBQdVue6Y2582bp9zcXE2bNs1p+wzAjcrCZVW/cVmm7OdI/vtf148HAACgEXNZuExOTtbNN9+sgoICLViwQAsWLND111+vGTNmKD8/v1L9rl27KjMzU1lZWSoqKlJqaqqio6Mr1ImOjtbKlSslSWvXrlXv3r1lsVgUHR2t1NRUFRUVKSsrS5mZmerWrVuNbS5fvlzp6el65pln5FXV7+MBqP8OHLD9vuX551dfp2y2knAJAADgUi47LHbEiBEaMGBAhbKioiL5+fkpJSWl8kB8fDRz5kzdfvvtKi0t1bBhw9S+fXs9//zz6tKliwYOHKjhw4dr6tSpiomJUYsWLfTss89Kktq3b6/rrrtO119/vby9vTVz5kx5e3tLUpVtStJDDz2kP/zhD/rzn/8sSYqJiVFycrKrng4ArnDwoNSqlS1gVsff3/Zbl4RLAAAAl7KYqk5YdIKhQ4faZxlrKmtIkpKSqgzGADxk2DBp927p229rrte7txQYKK1b555xAQAAnKNqykROn7k8ePCgrFarCgoK9N1339kvtnP8+HGdPHnS2d0BaMx++0264IIz14uMlL780vXjAQAAaMScHi7T09OVkpKi/fv36/HHH7eXN23aVJMnT3Z2dwAas4MHpc6dz1wvMlJatcr2m5icYw0AAOASTg+XQ4cO1dChQ7V27VoNGjTI2c0DwO8OHqz5Yj5loqKkoiLJapX+91u5AAAAcC6nh8vVq1dryJAhysnJ0aJFiyqtHzNmjLO7BNAYlZZKhw7V/rBYScrMJFwCAAC4iNPDZdl5lSdOnHB20wDwu9xcyZi6hcv//lfq08e14wIAAGiknB4ub7jhBkniZz0AuNbBg7b7uoZLAAAAuITLrmzx1FNP6fjx4youLtbo0aPVq1cvrV692lXdAWhs6hIumzWTgoMJlwAAAC7ksnD56aefKjAwUBs2bFBYWJjWrl2rV1991VXdAWhsfvvNdl+bcCnZZi8zM102HAAAgMbOZeGypKREkrRx40bFxcUpKCjIVV0BaIzKZi5rc7VYyRYumbkEAABwGZeFy6uvvlrXXnutMjIy1KdPH+Xm5qpJkyau6g5AY3O24dIY140JAACgEXP6BX3K3HvvvbrjjjsUGBgob29vBQQE6KWXXnJVdwAam4MHpRYtJD+/2tWPjJTy821XmW3VyrVjAwAAaIRcFi4l6ccff1ROTo5KS0vtZYmJia7sEkBjcfBg7c+3lKSoKNt9ZibhEgAAwAVcFi6nTp2qrKwsXXzxxfL29pYkWSwWwiUA56hruCz/cyT/93+uGRMAAEAj5rJwmZGRoffff18Wi8VVXQBozH777ffAWBtldbliLAAAgEu47II+7du318GyC24AgLPVdeYyONh2juaPP7puTAAAAI2Yy2Yu8/LyFBcXp27dusnX19devmDBAld1CaCxMMYWLmt7pVhJslikDh2kH35w3bgAAAAaMZeFywkTJriqaQCN3dGjUnFx3WYuJVu43LzZNWMCAABo5Fx2WGzPnj0VHh6ukpIS9ezZU127dlWnTp1c1R2AxqTskPuzCZe//CKdPOn8MQEAADRyLguXb7/9tiZOnKiZM2dKkqxWq8aPH++q7gA0Jo6ES0nau9e54wEAAIDrwuWSJUu0bNkyBQYGSpKioqKUm5vrqu4ANCa//Wa7r2u4vPhi2/133zl3PAAAAHBduPTz85Ofn599uaSkxFVdAWhsznbmsmNHydtbyshw/pgAAAAaOZdd0KdHjx5asGCBCgoK9Omnn2rp0qWKjo52VXcAGpOycFmXq8VKUpMmtkNjd+50/pgAAAAaOZfNXN57770KDg5Whw4d9NZbb2nAgAGaNGmSq7oD0JgcPCgFBEhNm9Z92y5dCJcAAAAu4LKZSy8vL11zzTW65pprFBwc7KpuADRGBw/W/ZDYMt26ScuX237OpHlz544LAACgEXN6uDTGaN68eXrzzTfty15eXrrllluUnJzs7O4ANEaOhMsePWz327dLV1/tvDEBAAA0ck4/LPa1117Tl19+qXfeeUdbt27V559/ruXLl+urr77S4sWLnd0dgMbIapVCQ89u27Jw+fnnzhsPAAAAnB8uV61apblz5yoiIsJeFhERoaefflqrVq1ydncAGqP9+6WwsLPbNjhYuugiaetW544JAACgkXN6uCwpKanyHMvg4GB+jgSA406dss1cnm24lKS+faX0dFtbAAAAcAqnh0tfX9+zWgcAtXLokFRa6li4jI62nbfJ710CAAA4jdMv6LN7925dfvnllcqNMSoqKnJ2dwAam/37bfeOhMuyC/mkpdmuHgsAAACHOT1c7tq1y9lNAsDvnBEu27aVOnWS3n1X+tvfnDMuAACARs7ph8UCgEs5I1xKUlKStGmT7fBYAAAAOIxwCaBh2bfPdu9ouBwxwnZBn6VLHR8TAAAACJcAGphffpFatJCaNXOsnW7dpJ49pZdfloxxztgAAAAaMcIlgIbl55+ldu2c01ZysrRrl7R6tXPaAwAAaMQIlwAaFmeGyxtvlDp0kKZNkwoLndMmAABAI0W4BNBwGCNlZkpRUc5pz8dHeu45afduaeZM57QJAADQSBEuATQcBw5IJ086b+ZSkq67TvrrX6WnnpLmzXNeuwAAAI2M03/nEgBc5uefbffODJeS9MILtqvQTphg62P2bCkgwLl9AAAAnOOYuQTQcHzzje2+c2fntuvnJ61YIY0bJz3zjNS1q/Tmm1JxsXP7AQAAOIcRLgE0HNu3Sy1bOu+cy/J8fKSXXpLWrZP8/aVRo6SLLrIdLpuT4/z+AAAAzjGESwANx7Zt0v/9n2SxuK6PgQOlHTuk996zhdj775ciIqSrr5Yee0zaulUqKXFd/wAAAA1UvQqXmzZt0qBBgxQTE6OFCxdWWl9UVKRJkyYpJiZGI0aMUHZ2tn3dyy+/rJiYGA0aNEibN28+Y5tZWVkaMWKEYmNjNWnSJBUVFbl25wA45uhRKSND6t7d9X15eUnx8dLGjdIPP0gzZki5udL06VLv3lJwsO1CQDNnSikp0qef2uodOiSVlrp+fAAAAPVQvbmgT2lpqWbNmqVFixYpNDRUw4cPV3R0tC666CJ7neXLl6t58+b6+OOPlZqaqjlz5ui5557T3r17lZqaqtTUVFmtVo0ZM0Zr166VpGrbnDNnjm699VbFxcVp5syZeuedd3TTTTd5avcBnMmKFbZzIIcMcW+/7dtLjzxiux08KG3YIKWlSf/5j/TRR9KpU5W3adHCFkCDg22H8bZsWfHx6beydc2a2YItAABAA1RvwuWOHTsUGRmpiIgISVJcXJzS0tIqhMv169crOTlZkjRo0CDNmjVLxhilpaUpLi5Ofn5+ioiIUGRkpHbs2CFJVbb5xz/+UZ999pnmzp0rSRo6dKjmzZtHuATqq6Ii6cUXbedA9urluXFccIE0YoTtJkn5+dL339tC58GDttnN3FwpL892f+iQ7XFWlu0+L6/miwR5eUlBQbagGRRkO/fT319q0sR2q+rx2awvX+bnJ3l72/o+/QYAAFAH9SZcWq1WhYWF2ZdDQ0PtAbF8ndatW0uSfHx81KxZM+Xl5clqteqSSy6psK3VapWkKtvMy8tT8+bN5ePjY69TVr/BOXlSWrZMOnGidvWNqVv7dalP2/V7LA217ZIS2wzhV19J77zj2vMt66ppU+nyy2tf3xhbIC0LmuVvZaG07Hb4sFRYaHuPlz0uKLDdn/7YVaoKnGW36gJpXep4qg2Lpfq/o4ZeXh/HxHPhvnJntwXAszp0kGJiPD2KOqk34dJU8WXTctoHXnV1qis/VcXhaqe3eabyeu+bb6Q77qj60DzgXNGunfTqq9KwYZ4eiWMsFikw0Hb73xEVDjPGNrNbVeiszeOiItvnR9mttLTicm3XOXPbkhLX91vdc9mQy892GwBA/XTRRdKePZ4eRZ3Um3AZFham/fv325etVqtCQkIq1dm3b5/CwsJUUlKiY8eOKSgoqMZtqypv2bKljh49qpKSEvn4+Gj//v2V+moweve2zXjU5ff46hqk61K/obZdn8ZC25X5+tat7cbEYvn9EFfAWepbcHZH0G7o5c5uC4DnNWvm6RHUWb0Jl127dlVmZqaysrIUGhqq1NRU+zmRZaKjo7Vy5UpddtllWrt2rXr37i2LxaLo6GhNmTJFY8aMkdVqVWZmprp16yZjTJVtWiwW9erVS2vXrlVcXJxWrlyp6OhoD+25E7Ro4ekRAADOJRxCCQA4C/UmXPr4+GjmzJm6/fbbVVpaqmHDhql9+/Z6/vnn1aVLFw0cOFDDhw/X1KlTFRMToxYtWujZZ5+VJLVv317XXXedrr/+enl7e2vmzJny9vaWpCrblKSpU6fqb3/7m5577jl17NhRI8ou0AEAAAAAqDOLqeqERVSpV69eCg8P9/QwAAAAAMAjcnJytHXr1irXES4BAAAAAA7jh8wAAAAAAA4jXAIAAAAAHEa4BAAAAAA4jHAJAAAAAHAY4RIAAAAA4DDCJQAAAADAYYRLAAAAAIDDCJcAAAAAAIcRLgEAAAAADiNcAgAAAAAcRrgEAAAAADiMcAkAAAAAcBjhEgAAAADgMMIlAAAAAMBhPp4eQEPSq1cvhYeHe3oYAAAAAOAROTk52rp1a5XrCJd1EB4erpSUFE8PAwAAAAA8Iikpqdp1HBYLAAAAAHAY4RIAAAAA4DDCJQAAAADAYYRLAAAAAIDDCJcAAAAAAIcRLgEAAAAADiNcAgAAAAAcRrgEAAAAADiMcAkAAAAAcBjhEgAAAADgMMIlAAAAAMBhhEsAAAAAgMMIlwAAAAAAhxEuAQAAAAAOI1wCAAAAABxGuAQAAAAAOIxwCQAAAABwGOESAAAAAOAwwiUAAAAAwGGESwAAAACAwwiXAAAAAACHES4BAAAAAA4jXAIAAAAAHEa4BAAAAAA4jHAJAAAAAHAY4RIAAAAA4DDCJQAAAADAYYRLAAAAAIDDCJcAAAAAAIcRLgEAAAAADiNcAgAAAAAcRrgEAAAAADiMcAng3BEcLFksNd+Cgz09SgAAgHOSj6cHAABOk5cnGfP7ssVScbmsDAAAAE7HzCUAAAAAwGE1hsujR49qyZIlNTaQnZ2t995774wdZWdnKz4+vm6jc6La7AuABqZLF5c1bbFYuHFrVLeIiAi1atXqjHWWLVumZcuWqUuXLvL29lZERIQCAwPtdXx9fdWtWzd5eXlV2DYwMFATJkxQly5d5CNS/I8AAByxSURBVOXlJV9f3wrbeHl5qUuXLlq2bJn9fThhwgT5+/vLYrHI399fEyZMkKRK/Zcfd0REhL0fb29ve5vltylfFhERUWkcERERioiIqFC3vKracpf63LcrxubJ/QVwFkwNsrKyTFxcXE1VzGeffWbuvPPOGuvUti1Xckb/Q4cOddJoADjF6R9hZ1qurqxSFXHj1mhu3t7e9ntvb28TFRVVYb2vr6/9PigoyDRr1sxccMEFZv369eb11183AQEBRpK58cYbzfvvv2+aNGli37Zp06bmtttuMxaLxV6WmJhoLrjgAtOyZUvTokUL06JFCxMcHGzOP/98M336dNOuXTuzdOlSk5ycbHx8fMzcuXNNfn6+mTt3rvHx8TGxsbGmXbt29v6DgoKMl5eXeeCBB8xHH31kX54+fbopKioy69evNxdccIF9zOXLmjdvbsLCwsy0adNMWFiYCQ4ONoGBgaZFixamdevW5vXXXzfr16+3j8kYY5YuXWrvv6yt8utdqT737YqxeXJ/AVSvpkxU47esSZMmma5du5rBgwebJ554wjzxxBMmLi7OxMfHm9TUVGOMMSNGjDCXX365GTx4sFm0aJHJysoyN954o0lMTDSJiYlm+/btxpgzh7uSkhLzxBNPmPj4eBMfH29ef/11Y4wx//nPf8yQIUNMfHy8eeCBB0xhYaExxpirr77aHDp0yBhjzI4dO8wtt9xijDHmhRdeMA888IC55ZZbTHR0tHnttdeq3Ber1WpuuukmM3jwYBMXF2e++OILh55IAB4gwiU3buVvXl5exsvLq8b1koy/v7992dfX114+d+5cI8m0bNnShIWF2cPmuHHjjCQTFRVlfH19TVRUlDHGmM6dOxtfX18zbtw407lzZ2OMsYdRSWb9+vXGGGNvV5Jp0qSJiYqKMuvXrzdRUVEVHnfu3NmsX7/edO7c2TRp0sTMnTu3wntz7ty5xmKx2Nvt3LmziYqKMnPnzrX3HxUVZUJDQ+3LZWVlYy5f5uvra+9v/fr1Zv369fb9Kys3xlR4XFa3vPLrXak+9+2KsXlyfwFU76zDZflA+OGHH5pbb73VlJSUmIMHD5oBAwYYq9VaaebyxIkTpqCgwBhjzM8//2zv/EzhcsmSJSY5OdkUFxcbY4zJy8szBQUFpn///uann34yxhgzdepUs2jRImNMzeHyz3/+syksLDSHDh0yPXv2NEVFRZX6f/XVV81LL71kjLEF22PHjtX0VBhjCJdAvSNVvp2+vjbbnHbzdEDgxs1Tt/z8fPvjssBpsVjM4cOH7Y/L7o0x9jqHDx82Xl5e/3uL/d5eUVGRMcZUaLds+6KiImOxWIyXl5cpKiqyB+Oyx2XjKa+snbJ2vby8jMViMfn5+fb+y8rKlo0xlZbLj72sv6KiIlNUVFRhfGXblH9cVre88utdqT737YqxeXJ/AVSvpkxU6wv6bN++XXFxcfL29tb555+vHj16aOfOnZXqlZSU6O9//7sSEhJ0zz336Mcff6xV+1u2bNENN9wgHx/bBWyDgoL0888/q02bNmrXrp0kaejQodq2bdsZ2xowYID8/PwUHBys4OBgHTp0qFKdrl27KiUlRS+++KJ++OEHBQYG1mqcAOqZ8rHwbLap6gY0UF5eXvLyqv6f9rJ1/v7+9uWy8wwlacGCBZKkli1bKiQkRN7e3oqMjNS0adMkSZGRkfL19VVkZKQkqWPHjvL19dW0adPUsWNHSZKvr6+9v/T09ArtSlKTJk0UGRmp9PR0RUZGqm3btkpPT1fbtm3VsWNHpaenq2PHjmrSpEmF7crasVgs9nY7duyoyMhILViwwN5/27ZtFRISYl8uG3fbtm0rtNW2bVv5+vra+0tPT1d6erp9/8rKy/aj7HFZ3fLKr3el+ty3K8bmyf0FcHZqHS5NLb9wLV68WOeff75Wr16tFStWqLi4uNbtWyyWSmXV8fb2tq8vLCyssM7Pz69CvZKSkkrb9+jRQ2+++aZCQ0N13333adWqVbUaJwAA9dWpU6d06tSpGtdLUkFBQYX6p06dkre3t+677z5FRUUpLy9P+/fvV2lpqXJycjR//nz5+vrq8OHD8vf3V35+vj755BPdf//98vHx0fz589WtWzd98MEHFcJtQkKCxo4dq3vvvddedt111yk/P1/Dhg1TXl6e8vLyNHz4cB0/flyJiYkaO3aspk+frjvuuEP333+/nnnmGZ04cULPPPOM7r//fsXExGjs2LH2/g8fPqypU6cqISFBH3/8sQ4fPqyDBw8qMTFRxcXF+uSTT5Sfn28fc/mygIAA3XTTTRo8eLBuuukmDR8+XE2aNFFeXp5uvvlm3X///frkk0/sY5Kk6dOn2/sva6v8eleqz327Ymye3F8AZ6mmKc/c3Fxz1VVXGWOMWbt2rbnttttMSUmJOXTokLnqqqvMgQMHzM6dO83NN99s3+bRRx81r776qjHGmHfeecd06NDBGHPmw2KXLl1qJkyYUOmw2AEDBpjMzExjjDH333+/Wbx4sTHGmNGjR5sNGzbY+yx/WOwrr7xibzcuLs5kZWVV2BdjjMnOzrb3tWjRIjN79uyangpjDIfFAvXO6efdnP6RVtVHXM0fe+Wqef4QRW7c3Hlr06aNCQ4OPmOdpUuXmqVLl5rOnTsbLy8v06ZNG9O0aVN7HR8fH9O1a9cKF/GRbBf3SU5ONp07dzYWi8X4+PhU2MZisZjOnTtXuFhLcnKy/QJBTZo0McnJycYYU6n/8uNu06aNvR8vLy97m+W3KV/Wpk2bSuNo06aNadOmTYW65VXVlrvU575dMTZP7i+AqtWUiSzG1DwlOWXKFH3//fe68sorJUmbN2+WxWLRuHHjdP3116u4uFi333678vLylJSUpKuuukoTJkxQQECAevXqpTfffFNfffWVsrOzddddd2nNmjVV9lNSUqKnn35amzdvlo+Pj0aOHKlbbrlFW7Zs0ZNPPqnS0lJ16dJFjzzyiPz8/LRt2zZNnz5drVq10iWXXKKMjAy98cYbevHFF3Xeeedp7NixkqT4+HgtWLBAbdq0qbAvHTp00KuvviofHx+dd955evLJJxUREVHTU6GkpCSlpKTUWAeAB1ksqnBY6+nL1ZUBAACgVmrKRGcMl/gd4RKo5wiXAAAALlVTJvJx81gAwLVOO3e70nLLlu4bCwAAQCPi9nC5efNmzZkzp0JZmzZt9M9//tPdQwFwrmFGEgAAwGPcHi6vvPJK+/mbAAAAAIBzQ61/igQAAAAAgOoQLgEAAAAADiNcAgAAAAAcRrgEAAAAADiMcAkAAAAAcBjhEgAAAADgMMIlAAAAAMBhhEsAAAAAgMMIlwAAAAAAhxEuAQAAAAAOI1wCAAAAABxGuAQAAAAAOIxwCQAAAABwGOESAAAAAOAwwiUAAAAAwGGESwAAAACAwwiXAAAAAACHES4BAAAAAA4jXAIAAAAAHEa4BAAAAAA4jHAJAAAAAHAY4RIAAAAA4DDCJQAAAADAYYRLAAAAAIDDCJcAAAAAAIcRLgEAAAAADiNcAgAAAAAcRrgEAAAAADiMcAkAAAAAcBjhEgAAAADgMMIlAAAAAMBhPp4eQEOSk5OjpKQkTw8DAAAAADwiJyen2nUWY4xx41gAAAAAAOcgDosFAAAAADiMcAkAAAAAcBjhEgAAAADgMMIlAAAAAMBhhEsAAAAAgMMIl0704osv6sorr9SQIUM0ZMgQbdy40b7u5ZdfVkxMjAYNGqTNmzfbyzdt2qRBgwYpJiZGCxcutJdnZWVpxIgRio2N1aRJk1RUVCRJKioq0qRJkxQTE6MRI0YoOzvbfTvYwD355JO69tprlZCQoPHjx+vo0aOSpOzsbHXr1s3+us2cOdO+TUZGhhISEhQTE6PZs2er7OLKhw8f1pgxYxQbG6sxY8boyJEjkiRjjGbPnq2YmBglJCTo22+/df+OnuOqe8/A+fbt26dRo0bpuuuuU1xcnF577TVJ7vmsQ+1ER0crISFBQ4YMsf9U1tl8Pq1cuVKxsbGKjY3VypUr7eXVfQbizH766Sf7e2TIkCG6/PLLtXjxYt4/HjZt2jT16dNH8fHx9jJ3vGeq6wMVVfX68P2tgTFwmhdeeMG88sorlcr37NljEhISTGFhofnll1/MwIEDTUlJiSkpKTEDBw40v/zyiyksLDQJCQlmz549xhhjJk6caNasWWOMMWbGjBlmyZIlxhhj3nzzTTNjxgxjjDFr1qwx99xzj5v2ruHbvHmzKS4uNsYY89RTT5mnnnrKGGNMVlaWiYuLq3KbYcOGmS+//NKcOnXKjB071mzYsMEYY8yTTz5pXn75ZWOMMS+//LK9rQ0bNpixY8eaU6dOma+++soMHz7c1bvVqNT0noHzWa1Wk5GRYYwx5tixYyY2Ntbs2bPHLZ91qJ2rr77aHDp0qEJZXT+f8vLyTHR0tMnLyzOHDx820dHR5vDhw8aY6j8DUTclJSWmb9++Jjs7m/ePh33++ecmIyOjwr/77njPVNcHKqrq9eH7W8PCzKUbpKWlKS4uTn5+foqIiFBkZKR27NihHTt2KDIyUhEREfLz81NcXJzS0tJkjNFnn32mQYMGSZKGDh2qtLQ0SdL69es1dOhQSdKgQYO0ZcsW/ie5lq644gr5+PhIki699FLt37+/xvoHDhzQ8ePHddlll8lisSgxMdH+OqSlpSkxMVGSlJiYqHXr1lUot1gsuvTSS3X06FEdOHDAhXvVuFT3noFrhISEqHPnzpKkwMBAXXjhhbJardXWd+ZnHc5eXT+f0tPT1a9fPwUFBalFixbq16+fNm/eXONnIOpmy5YtioiIUHh4eLV1eP+4R48ePdSiRYsKZe54z1TXByqq6vXh+1vDQrh0siVLlighIUHTpk2zT7VbrVaFhYXZ64SGhspqtVZbnpeXp+bNm9vfSGFhYfYvdFarVa1bt5Yk+fj4qFmzZsrLy3PX7p0zVqxYof79+9uXs7OzlZiYqFtuuUXbtm2TVPl1K/86HDp0SCEhIZJsX8Bzc3PPuA0cV917Bq6XnZ2tXbt26ZJLLpHk+s861N7YsWOVlJSkt956S1LdP59q+7rx+py91NTUCof58f6pX9zxnqmuD9QN39/qPx9PD6ChufXWW/Xbb79VKp80aZJuvPFG3X333bJYLHr++ef1xBNP6PHHH69yZtFisejUqVNVllelrLy6tmBT0+tzzTXXSJLmz58vb29vDR48WJLtw+WTTz5Ry5YtlZGRofHjxys1NfWsnmteH9fi+fWM/Px8TZw4UQ8++KACAwPd8lmH2lm2bJlCQ0N16NAhjRkzRhdeeGG1dat7fepajropKirS+vXrNWXKFEni/dOA8J6pX/j+1jAQLuto8eLFtao3YsQI3XXXXZJs//tRfgrfarXa/9ekqvKWLVvq6NGjKikpkY+Pj/bv32+vHxYWpn379iksLEwlJSU6duyYgoKCnLR3Dd+ZXp+VK1dqw4YNWrx4sf1Dw8/PT35+fpKkLl26qG3btvr5558rvW7lX4dWrVrpwIEDCgkJ0YEDBxQcHCyp8mtdfhs4rqb3ElyjuLhYEydOVEJCgmJjYyVJ559/vn29qz7rUDuhoaGSbJ9JMTEx2rFjR50/n8LCwvT555/by61Wq3r27MnnmZNs2rRJnTt3tr9veP/UP+54z1TXB2qH728NB4fFOlH5Y7PXrVun9u3bS7JdzS81NVVFRUXKyspSZmamunXrpq5duyozM1NZWVkqKipSamqqoqOjZbFY1KtXL61du1aS7Q0VHR1tb6vsqmRr165V7969+Z+VWtq0aZP+9a9/af78+QoICLCX5+bmqrS0VJLsr09ERIRCQkLUtGlTff311zLGaNWqVRo4cKAk2+uwatUqSaqy3Bijr7/+Ws2aNePDyYmqe8/ANYwxmj59ui688EKNGTPGXu6Ozzqc2YkTJ3T8+HH7408//VTt27ev8+fTFVdcofT0dB05ckRHjhxRenq6rrjiiho/A1F7qampiouLsy/z/ql/3PGeqa4PnBnf3xoWi+FqME4zdepU7d69W5IUHh6uWbNm2f8w58+frxUrVsjb21sPPvigBgwYIEnauHGjHnvsMZWWlmrYsGEaN26cJNub5G9/+5uOHDmijh07as6cOfLz81NhYaGmTp2qXbt2qUWLFnr22WcVERHhmR1uYGJiYlRUVGSf6b3kkks0a9YsrV27Vi+88IK8vb3l7e2tCRMm2P+B3rlzp6ZNm6aCggL1799fM2bMkMViUV5eniZNmqR9+/apdevWev755xUUFCRjjGbNmqXNmzcrICBAjz32mLp27erJ3T7nVPeegfNt27ZNN998szp06CAvL9v/RU6ePFlr1qxx+WcdziwrK0vjx4+XJJWWlio+Pl7jxo07q8+nd955Ry+//LIk6a677tKwYcMkVf8ZiNo5efKkrrrqKq1bt07NmjWT5J7vCqje5MmT9fnnnysvL0+tWrXShAkTdM0117j8PVPd+xIVVfX6LFy4kO9vDQjhEgAAAADgMA6LBQAAAAA4jHAJAAAAAHAY4RIAAAAA4DDCJQAAAADAYYRLAAAAAIDDCJcAAAAAAIcRLgEADvnTn/6kqVOn2pdLSkrUu3dv/fWvf61xu5SUFM2aNcupY3nxxRf16quvOrVNR2RnZ+u9996r0zYLFiw4Y52jR49qyZIltWrvsssuq1P/zlab/akvyv9N1re/JQBoCAiXAACHnHfeedqzZ48KCgokSZ9++qlCQ0M9PKr6IScnR2vWrKnTNmU/zF6To0ePatmyZWc7LLeqzf44S0lJidv6AgBURrgEADisf//+2rBhgyQpNTVVcXFx9nWHDx/W3XffrYSEBI0cOVK7d++utH1ubq4mTJigYcOGadiwYdq+fbskKT8/X9OmTVNCQoISEhK0du1aSRVn4z788EM98MADldr85ZdfNHbsWCUlJemmm27Sjz/+WKlOde2vWbNGCQkJio+P19NPP22vX12/DzzwgGbPnq0bbrhBAwcO1IcffihJmjt3rrZt26YhQ4Zo8eLF2rNnj4YPH64hQ4YoISFBmZmZFcYzZ84cFRQUaMiQIZoyZYokadGiRYqPj1d8fLwWL15sb/eXX37RkCFD9OSTTyo/P1+jR4/W0KFDlZCQoHXr1lXxKlXtX//6lxISEjR48GDNmTNHkrRr1y6NHDlSCQkJGj9+vI4cOSJJGjVqlHbu3CnJ9ppFR0dLss34JScna+zYsYqNjdVTTz1V5f6cOHFCd955pwYPHqz4+Hi9//77lcYzatQoPfroo7rhhhsUHx+vHTt2SJJOnDihadOmadiwYUpMTLTvY0pKiiZOnKi77rpLt912W6X2Vq1aZd+/shn26v7eAACO8fH0AAAADd/111+vl156SVdffbW+//77Cl/YX3zxRXXq1EkvvfSStmzZovvvv1+rV6+usP2jjz6q0aNHq3v37vr11181duxYffDBB3rppZcUGBhoP7S0LOTUxowZM/TII48oKipK33zzjR555BG9/vrrFepU1b7VatWcOXOUkpKi5s2b67bbbtO6det0zTXX1NjfgQMHtHTpUv30008aN26crr32Wk2ZMkX//ve/7bN3//jHP/SXv/xFgwcPVlFRkU6dOlWhjXvvvVdLliyxPz8ZGRlKSUnR22+/LWOMRo4cqZ49e2rKlCnas2ePvV5JSYn++c9/KjAwULm5ufrzn/+sgQMHymKx1DjmjRs3Ki0tTW+//bYCAgJ0+PBhSdJ9992nGTNmqGfPnnr++ec1b948TZ8+vca2du3apVWrVsnPz0/XXnutRo0aVWl/1q5dq5CQEC1cuFCSdOzYsSrbOnnypP7f//t/+uKLL/Tggw9qzZo1WrBggXr37q3HH39cR48e1YgRI9S3b19J0tdff613331XQUFBFdrZs2eP5s+fr2XLlik4ONi+f9X9vQEAHEO4BAA47OKLL1Z2drbWrFmjAQMGVFi3fft2vfjii5KkPn366PDhw5VCxX/+8x/t3bvXvnz8+HEdP35cW7Zs0TPPPGMvb9GiRa3Gk5+fr6+++kr33HOPvayoqKhSvara/+KLL9SzZ08FBwdLkhISEvTFF1+cMVxec8018vLy0kUXXaTffvutyjqXXnqpFixYoP379ys2NlZRUVE1trl9+3Zdc801Ou+88yRJMTEx2rZtm33GsIwxRs8884y++OILeXl5yWq16rffftMFF1xQY/tbtmxRUlKSAgICJElBQUE6duyYjh07pp49e0qShg4dWuF5rE6fPn3UrFkzSdIf//hH5eTkqHXr1hXqdOjQQU8++aSefvppXX311erevXuVbZXNfPfo0UPHjx/X0aNHlZ6ervXr1+vf//63JKmwsFD79u2TJPXr169SsJSkzz77TNdee639tSyrU93fGwDAMYRLAIBTREdH66mnntLrr79unyGSbMHndKfPqJ06dUpvvfWW/P39K5QbY844+1ZYWFipzBij5s2bV5ohraremdqvbb9+fn5n3CYhIUGXXHKJNmzYoLFjx2r27Nnq06dPjeOrjffee0+5ublKSUmRr6+voqOjq3xeqmq/Lvvv7e1tH9PpYb38/nt7e6u0tLTS9u3atVNKSoo2btyouXPnql+/fkpOTq5U7/QxlS2/8MILuvDCCyus++abb+zh+HTVPX/V/b0BABzDOZcAAKcYPny47r77bv3pT3+qUN6jRw+9++67kqStW7eqZcuWCgwMrFDniiuu0Jtvvmlf3rVrlyTbjFT58rLDYs8//3z9+OOPOnXqVJXnFwYGBqpNmzb2Qx2NMVWe61lV+926ddMXX3yh3NxclZaWKjU1VT169KhVv6dr2rSp8vPz7ctZWVmKiIjQX/7yF0VHR+v777+vtI2Pj4+Ki4sl2Z67devW6eTJkzpx4oTWrVun7t27V2r32LFjatWqlXx9ffXZZ58pJyfnjGMr2/8VK1bo5MmTkmznxzZr1kzNmzfXtm3bJEmrV6+27394eLgyMjIkyX5e6ZmU3x+r1aqAgAANGTJEY8eO1XfffVflNmXnYm7btk3NmjVTs2bN7H8jZYGxum3L69Onjz788EPl5eXZ90+q/u8NAOAYZi4BAE4RFham0aNHVypPTk62XzQnICBATzzxRKU606dP16xZs5SQkKDS0lJ1795ds2bN0rhx4zRr1izFx8fLy8tLycnJio2N1ZQpU/TXv/5VrVu3Vvv27XXixIlKbT799NN6+OGHNX/+fJWUlOj666/XxRdfXKFOde1PnjxZo0ePljFG/fv3tx8SW5t+y/vTn/4kb29vDR48WElJSSosLNS7774rHx8fnX/++Ro/fnylbUaOHKnBgwerU6dOmjt3rpKSkjRixAhJtgDfqVMnSdLll1+u+Ph4XXnllbrjjjs0btw4JSUlqWPHjpVm96rTv39/7d69W8OGDZOvr68GDBigyZMn68knn9RDDz2kkydPKiIiQo8//rgk6bbbbtOkSZP07rvvqlevXrXqo/z+JCYm6qmnnpKXl5d8fHz08MMPV7lNixYtdMMNN+j48eN67LHHJEl33323HnvsMQ0ePFjGGIWHh5/xSrTt27fXXXfdpVGjRsnLy0udOnXSE088Ue3fGwDAMRZT22NuAAAAXGzUqFG677771LVrV08PBQBQRxwWCwAAAABwGDOXAACc477//nvdd999Fcr8/Py0fPlyD40IAHAuIlwCAAAAABzGYbEAAAAAAIcRLgEAAAAADiNcAgAAAAAcRrgEAAAAADiMcAkAAAAAcNj/B8QxFf+noXgrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "total_molecules_per_cell(adata1, count=\"total_counts\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "educated-wrist",
   "metadata": {},
   "source": [
    "### Sparse gene matrix\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "sorted-sheet",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:31:26.108833Z",
     "start_time": "2020-04-10T16:31:25.794298Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<109655x19984 sparse matrix of type '<class 'numpy.float32'>'\n",
       "\twith 146938269 stored elements in Compressed Sparse Row format>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mtx_file = sparse.csr_matrix(adata1.X)\n",
    "\n",
    "mtx_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "further-patch",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write sparse matrix\n",
    "PATH = #insert out path here\n",
    "PATH += \"input/\"\n",
    "FILE = \"counts_v1.mtx\"\n",
    "\n",
    "scipy.io.mmwrite(PATH + FILE, mtx_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "committed-acrylic",
   "metadata": {},
   "source": [
    "### Gene names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "sticky-lease",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:32:56.521777Z",
     "start_time": "2020-04-10T16:32:56.516194Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0              UBTF\n",
       "1             KCNC3\n",
       "2        AL157838.1\n",
       "3            MINDY3\n",
       "4             CDC5L\n",
       "            ...    \n",
       "19979         AZGP1\n",
       "19980         TSHZ3\n",
       "19981         HMGA2\n",
       "19982       SLC35A1\n",
       "19983        ZC3H7B\n",
       "Length: 19984, dtype: object"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g_file = pd.Series(adata1.var_names)\n",
    "\n",
    "g_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "green-disclaimer",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write gene names\n",
    "PATH = #out path here\n",
    "PATH += \"input/\"\n",
    "FILE = \"genes.csv\"\n",
    "\n",
    "g_file.to_csv(PATH + FILE, index=False, header=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "brown-payment",
   "metadata": {},
   "source": [
    "### Cell names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "healthy-empire",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:33:34.915287Z",
     "start_time": "2020-04-10T16:33:34.911705Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['JUQ024_AAACCTGCAACAACCT-1', 'JUQ024_AAACCTGCAAGAAGAG-1',\n",
       "       'JUQ024_AAACCTGGTCTCCACT-1', 'JUQ024_AAACCTGTCAACGAAA-1',\n",
       "       'JUQ024_AAACGGGCACAGAGGT-1', 'JUQ024_AAACGGGTCAACGAAA-1',\n",
       "       'JUQ024_AAACGGGTCCGAACGC-1', 'JUQ024_AAACGGGTCGTAGGAG-1',\n",
       "       'JUQ024_AAAGATGAGCAGGCTA-1', 'JUQ024_AAAGATGCAACTGCTA-1',\n",
       "       ...\n",
       "       'JUQ124_TTCGAAGCAGTCGTGC-1', 'JUQ124_TTCTACAAGCTATGCT-1',\n",
       "       'JUQ124_TTCTCCTGTACCAGTT-1', 'JUQ124_TTCTCCTGTACCGTAT-1',\n",
       "       'JUQ124_TTCTTAGGTGGTTTCA-1', 'JUQ124_TTCTTAGTCACCATAG-1',\n",
       "       'JUQ124_TTTACTGCACCAGATT-1', 'JUQ124_TTTACTGCACGGTGTC-1',\n",
       "       'JUQ124_TTTACTGGTAGCTAAA-1', 'JUQ124_TTTATGCAGCTAGTCT-1'],\n",
       "      dtype='object', length=109655)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata1.obs_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "worse-baker",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:33:46.368004Z",
     "start_time": "2020-04-10T16:33:46.363225Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0         JUQ024_AAACCTGCAACAACCT-1\n",
       "1         JUQ024_AAACCTGCAAGAAGAG-1\n",
       "2         JUQ024_AAACCTGGTCTCCACT-1\n",
       "3         JUQ024_AAACCTGTCAACGAAA-1\n",
       "4         JUQ024_AAACGGGCACAGAGGT-1\n",
       "                    ...            \n",
       "109650    JUQ124_TTCTTAGTCACCATAG-1\n",
       "109651    JUQ124_TTTACTGCACCAGATT-1\n",
       "109652    JUQ124_TTTACTGCACGGTGTC-1\n",
       "109653    JUQ124_TTTACTGGTAGCTAAA-1\n",
       "109654    JUQ124_TTTATGCAGCTAGTCT-1\n",
       "Length: 109655, dtype: object"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bc_file = pd.Series(adata1.obs_names)\n",
    "\n",
    "bc_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "first-doctrine",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write gene names\n",
    "PATH = #out path here\n",
    "PATH += \"input/\"\n",
    "FILE = \"cells.csv\"\n",
    "\n",
    "bc_file.to_csv(PATH + FILE, index=False, header=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "auburn-worst",
   "metadata": {},
   "source": [
    "## Run scran\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "interpreted-hunger",
   "metadata": {},
   "source": [
    "On the terminal do:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "listed-mattress",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:35:25.353142Z",
     "start_time": "2020-04-10T16:35:25.348887Z"
    }
   },
   "source": [
    "```\n",
    "PATH_I=~/input/\n",
    "PATH_O=~/output/\n",
    "PATH_T=~/something_here\n",
    "\n",
    "m_file=combo_fixed_counts_v1.mtx\n",
    "g_file=genes.csv\n",
    "bc_file=cells.csv\n",
    "ofile=sct_combo_fixed_counts_v1.csv\n",
    "tfile=run_scran.R\n",
    "\n",
    "Rscript $PATH_T$tfile $PATH_I$m_file $PATH_I$g_file $PATH_I$bc_file $PATH_O$ofile\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "multiple-robinson",
   "metadata": {},
   "source": [
    "## Load scran output\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "economic-retro",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:37:58.722752Z",
     "start_time": "2020-04-10T16:37:57.863995Z"
    }
   },
   "outputs": [],
   "source": [
    "import anndata2ri\n",
    "\n",
    "anndata2ri.activate()\n",
    "\n",
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "committed-button",
   "metadata": {},
   "source": [
    "Install following R packages only first time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "configured-symbol",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:50:38.402986Z",
     "start_time": "2020-04-10T16:50:38.395640Z"
    }
   },
   "outputs": [],
   "source": [
    "#%%R\n",
    "\n",
    "#paste(\"R\", \"T\", sep = \"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "remarkable-enterprise",
   "metadata": {},
   "outputs": [],
   "source": [
    "#%%R\n",
    "\n",
    "#BiocManager::install(\"S4Vectors\")\n",
    "#BiocManager::install(\"SummarizedExperiment\")\n",
    "#BiocManager::install(\"SingleCellExperiment\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "distributed-thickness",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:55:09.106054Z",
     "start_time": "2020-04-10T16:55:00.366063Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "class: SingleCellExperiment \n",
      "dim: 19984 105773 \n",
      "metadata(0):\n",
      "assays(2): X logcounts\n",
      "rownames(19984): UBTF KCNC3 ... SLC35A1 ZC3H7B\n",
      "rowData names(1): value\n",
      "colnames(105773): JUQ024_AAACCTGCAACAACCT-1 JUQ024_AAACCTGCAAGAAGAG-1\n",
      "  ... JUQ124_TTTACTGGTAGCTAAA-1 JUQ124_TTTATGCAGCTAGTCT-1\n",
      "colData names(0):\n",
      "reducedDimNames(0):\n",
      "spikeNames(0):\n",
      "altExpNames(0):\n"
     ]
    }
   ],
   "source": [
    "%%R -o scr_matrix\n",
    "\n",
    "library(Matrix)\n",
    "\n",
    "PATH1 <- #data path here\n",
    "PATH2 <- \"/output/\"\n",
    "\n",
    "PATH <- paste(PATH1, PATH2, sep = \"\")\n",
    "\n",
    "FILE <- \"sct_counts_v1.rds\"\n",
    "\n",
    "FILE <- paste(PATH, FILE, sep = \"\")\n",
    "\n",
    "scr_matrix <- readRDS(FILE)\n",
    "\n",
    "scr_matrix <- as(scr_matrix, 'SingleCellExperiment')\n",
    "\n",
    "scr_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "indie-vintage",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:55:39.308901Z",
     "start_time": "2020-04-10T16:55:39.304345Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,\n",
       "        0.        ],\n",
       "       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,\n",
       "        0.        ],\n",
       "       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,\n",
       "        0.        ],\n",
       "       ...,\n",
       "       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,\n",
       "        0.        ],\n",
       "       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,\n",
       "        0.        ],\n",
       "       [1.94247272, 0.        , 0.        , ..., 0.        , 0.        ,\n",
       "        0.        ]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Add log counts from scran\n",
    "scr_matrix.X = scr_matrix.layers[\"logcounts\"]\n",
    "\n",
    "scr_matrix.X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "configured-extra",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:56:13.423235Z",
     "start_time": "2020-04-10T16:56:13.328419Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AnnData object with n_obs × n_vars = 105773 × 19984"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata1scr = sc.AnnData(X=np.array(scr_matrix.X),\n",
    "                       obs={'obs_names': scr_matrix.obs_names},\n",
    "                       var={'var_names': scr_matrix.var_names})\n",
    "\n",
    "adata1scr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "shaped-skirt",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:56:28.989169Z",
     "start_time": "2020-04-10T16:56:28.772346Z"
    }
   },
   "outputs": [],
   "source": [
    "### QC metrics\n",
    "sc.pp.calculate_qc_metrics(adata1scr,\n",
    "                           expr_type=\"counts\",\n",
    "                           var_type=\"genes\",\n",
    "                           inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "lesbian-asthma",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:56:38.855895Z",
     "start_time": "2020-04-10T16:56:38.818688Z"
    }
   },
   "outputs": [],
   "source": [
    "a = adata1scr.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "individual-greece",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-04-10T16:56:51.457580Z",
     "start_time": "2020-04-10T16:56:50.688792Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Median total count of UMI is: 1795.3435\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5cAAAE+CAYAAAAZGeTqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/S0lEQVR4nO3deVyVZf7/8fdhUxPN0ACHHB1Ly738uaSOOKK4AYq4ZE1mZlnmOpqVWY5ZmpY6lS3IVLYZVsZgSS6JqWllaotSWloxIROHEjdQVu/fH/eXk+zggXOzvJ6Px/2Qc597+Ry4wfM+13Vfl80wDEMAAAAAADjBzeoCAAAAAAA1H+ESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE7zsLqAmqRnz54KCAiwugwAAAAAsERycrL27t1b7HOEywoICAhQTEyM1WUAAAAAgCUiIiJKfI5usQAAAAAApxEuAQAAAABOI1wCAFCczMzSHwMAgAK45xIAgOLUry/ZbH88NgzragEAoAag5RIAgPK4uOWSVkwAAIogXAIAUB75LZkXt2bmI2wCAFAzw2VeXp7Cw8N19913S5JOnTqliRMnatCgQZo4caJOnz7t2Hb16tUKDg7W4MGD9cknnzjWJyQkKCwsTMHBwXr88cdl0N0JAOq2igTEi4OmzWY+BgCgjquR4fL111/X1Vdf7XgcFRWlXr16aevWrerVq5eioqIkSceOHVNcXJzi4uL00ksv6dFHH1VeXp4kaeHChVq0aJG2bt2qxMRE7dq1y5LXAgCoJgoHRgAAUCE1LlympKRox44dGj16tGNdfHy8wsPDJUnh4eHatm2bY31ISIi8vLzUokULtWzZUgcPHlRqaqrS09N1ww03yGazKTw8XPHx8Va8HAAAAACoFWpcuFyyZInmzp0rN7c/Sj9x4oR8fX0lSb6+vkpLS5Mk2e12+fv7O7bz8/OT3W4vst7f3192u91FrwAAUC1wnyQAAJWqRoXLjz/+WD4+PurYsWO5ti/uPkqbzVbiegBAHUI3WAAAKlWNmufyyy+/1Pbt27Vr1y5lZWUpPT1d9913n5o2barU1FT5+voqNTVVPj4+kswWyZSUFMf+drtdvr6+RdanpKQ4Wj4BALVUZmbVDbxT+NhVeS4AAKqpGtVyOWfOHO3atUvbt2/XypUrdeONN2r58uUKCgpSbGysJCk2NlYDBgyQJAUFBSkuLk7Z2dlKSkpSYmKiOnfuLF9fXzVs2FBff/21DMMosA8AoJaqypZKRo8FAKBmtVyWZPLkyZo1a5bWr1+v5s2b65lnnpEktWnTRkOHDtWwYcPk7u6uBQsWyN3dXZI5Wuy8efOUmZmpwMBABQYGWvkSAAAAAKBGsxkWTvD4ww8/qG3btladvsIiIiIUExNjdRkAgEt1cYulYVz64/JsCwBALVRaJrK05fKf//yncnJyNHLkSIWFhalx48ZWlgMAAAAAuESWhsvo6GglJibqvffe06hRo9S5c2dFRESoT58+VpYFAAAAAKggy++5bNWqlWbNmqWOHTvq8ccf13fffSfDMDR79mwNGjTI6vIAAAAAAOVgabg8cuSIYmJitHPnTvXu3VuRkZHq0KGD7Ha7xo0bR7gEADiHKUEAAHAZS8PlY489pjFjxmj27Nmqf9F//n5+fpo5c6aFlQEAaoX8KUIk1w6yw7yXAIA6yNJwGRUVpfr16zumB7lw4YKysrLUoEEDhYeHW1kaAACX7uJQKzF6LACgTnCz8uQTJ05UZmam4/H58+c1ceJECysCAAAAAFwKS8NlVlaWGjZs6HjcsGFDnT9/3sKKAAAAAACXwtJw2aBBA3377beOxwkJCQXuvQQAoEIu6g0DAABcy9J7Lh966CHNnDlTvr6+kqTffvtN//rXv6wsCQBQk3GvIwAAlrE0XHbu3FmbNm3Szz//LMMw1Lp1a3l6elpZEgAAAADgElgaLiXp0KFDSk5OVl5eng4fPixJjBQLAAAAADWMpeFy7ty5SkpK0nXXXeeYjsRmsxEuAQC1C/NeAgDqAEvDZUJCgj788EPZLr4/BgCA8qopIY17QQEAdYClo8W2adNGv/32m5UlAABqsvzQlr8AAADLWNpyefLkSYWEhKhz584FBvKJjIy0sCoAAAAAQEVZGi6nT59u5ekBAAAAAJXE0nDZo0cPJScn67///a969+6t8+fPKy8vz8qSAADVWU25xxIAgDrI0nsu33nnHc2YMUMLFiyQJNntdk2dOtXKkgAA1VltuccyM7P4rwEAqMEsDZdr165VdHS0vL29JUmtWrVSWlqalSUBAFD1Lg7JtMQCAGoJS8Oll5eXvLy8HI9zc3MtrAYAAAAAcKksveeye/fuioyMVGZmpvbs2aO33npLQUFBVpYEAKhOuMcSAIAaw9KWy/vuu08+Pj5q27at3n77bfXr10+zZs2ysiQAQHVSW+6xBACgDrC05dLNzU1jx47V2LFjrSwDAAAAAOAkS8NlUFCQbMV8Eh0fH29BNQAAAACAS2VpuHzvvfccX2dnZ2vTpk06ffp0qfv8+uuvuv/++/X77787Wj4nTJigU6dO6R//+IeSk5MVEBCgp59+WpdffrkkafXq1Vq/fr3c3Nz08MMPq2/fvpKkhIQEzZs3T5mZmerXr5/mz59fbNgFALgQ91kCAFAjWXrP5RVXXOFY/Pz8dPvtt+vzzz8vdR93d3c9+OCD2rRpk95++2299dZbOnbsmKKiotSrVy9t3bpVvXr1UlRUlCTp2LFjiouLU1xcnF566SU9+uijysvLkyQtXLhQixYt0tatW5WYmKhdu3ZV+WsGAJTh4vss64LC81wy7yUAoIayNFx+++23juXQoUOKjo5WRkZGqfv4+vqqQ4cOkiRvb2+1bt1adrtd8fHxCg8PlySFh4dr27ZtkswutiEhIfLy8lKLFi3UsmVLHTx4UKmpqUpPT9cNN9wgm82m8PBwuuMCAFyv8KBFtNoCAGooS7vFLl261PG1h4eHoztreR0/flyHDx9Wly5ddOLECfn6+koyA2haWpokyW63q0uXLo59/Pz8ZLfb5eHhIX9/f8d6f39/2e12J18RAKDC6AYLAECtYGm4fOONNy5534yMDM2YMUMPPfSQvL29S9zOMIwi62w2W4nrAQAult9yl6+Yv88AAKD6szRcrlmzptTnJ06cWOz6nJwczZgxQ2FhYRo0aJAkqWnTpkpNTZWvr69SU1Pl4+MjyWyRTElJcexrt9vl6+tbZH1KSoqj5RMAUIVoqQQAoFay9J7LhIQERUdHy263y263a926dTp27JgyMjJKvPfSMAzNnz9frVu3LhA+g4KCFBsbK0mKjY3VgAEDHOvj4uKUnZ2tpKQkJSYmqnPnzvL19VXDhg319ddfyzCMAvsAAKpQ4XsMAQBArWBpy+XJkycVExPj6NY6bdo0zZw5U4sXLy5xnwMHDmjDhg1q27atRowYIUmaPXu2Jk+erFmzZmn9+vVq3ry5nnnmGUlSmzZtNHToUA0bNkzu7u5asGCB3N3dJZmjxeZPRRIYGKjAwMAqfsUAAAAAUDtZGi7/97//ycvLy/HYy8tLycnJpe7TrVs3ff/998U+99prrxW7fsqUKZoyZUqR9Z06ddLGjRsrUDEAAAAAoDiWhssRI0Zo9OjRCg4Ols1m00cffeSYTgQAUEtwj2XFFP5+8f0DANQQlobLKVOmKDAwUPv375ckPfHEE2rfvr2VJQEAKhujwVYM3y8AQA1l6YA+knT+/Hl5e3trwoQJ8vf3V1JSktUlAQCckZlpdQUAAMAClobL5557Ti+99JKioqIkmVOMzJ0718qSAADOYjRYAADqJEvD5UcffaQXX3xRDRo0kCT5+fmVOAUJAKAao7USAIA6z9Jw6enpKZvNJtv/fbJ97tw5K8sBAFyqi1srAQBAnWTpgD5Dhw7VggULdObMGb3zzjt67733NHbsWCtLAgCgemH0WABADWFZuDQMQ8OGDdNPP/2khg0b6ueff9aMGTPUp08fq0oCAJQXAcd1GD0WAFBDWBYubTabpk6dqpiYGAIlAFR3hcMkgQcAABRi6T2XXbp00cGDB60sAQBQHowAW30UHjyJwZQAANWEpfdc7t27V+vWrVNAQIBjxFhJ+uCDDyysCgCAaoxWYwBANWVJuPzf//6nP/3pT/r3v/9txekBAIUV7vZ6/rx00Yd+AAAAZbGkW+zUqVMlSQEBAVq6dKkCAgIKLAAAFyvc7bVBA7rBAgCACrEkXBoXdeFJSkqyogQAAGqHi++5tPL+y8LnPn+++K+Le8x9owBQK1gSLm0XfQpu4xNxAHCNi9/AF35zj5rr4lbnwtPDlDX4jzOBsPDj0lq/C7eEF37MtDYAUCtYcs/lkSNH1LVrVxmGoaysLHXt2lWS2aJps9n05ZdfWlEWANQupU0fYhgMClMblTVlTEmBMN/F10Vx10hZjyur7sL3/DKvKgDUCJaEy8OHD1txWgBwjYvfCBd+U1zWm+jS3lRXdF9GFa17yvqZV9drorSQm/8YAFDtWTrPJQDUSKV1JZQKdg8srKyBc0obSMeZfYGajLk9AaBGIFwCQEVVZGTVwts6cy6grirt94igCQDVhiXdYgEAAC5Z4fuHAQDVAi2XAFAWWkaA6osuswBQbdByCQBlqa6DoADg9xMAqhFaLgGgMFo+gJqLlkwAsAwtlwBQGC0hQM3F7y8AWIaWSwCQaN0Aaquypg7idx8AKk2dD5e7du3S4MGDFRwcrKioKKvLAeAqhd9QMuUHUDtVZOoggiYAOKVOh8u8vDwtWrRIL730kuLi4rRx40YdO3bM6rIAVIXSwiSBEqi7SvtgibAJABVSp8PlwYMH1bJlS7Vo0UJeXl4KCQlRfHy81WUBqAyESQAVVdbficJdauliCwAF1Olwabfb5e/v73js5+cnu91uYUUALhlhEkBlK6tLbWldbAmiAOogm2HU3WHUNm3apN27d2vx4sWSpNjYWB06dEiPPPJIsdv37NlTAQEBriwRAAAAAKqN5ORk7d27t9jn6vRUJP7+/kpJSXE8ttvt8vX1LXH7kr6JAAAAAFDX1elusZ06dVJiYqKSkpKUnZ2tuLg4BQUFWV0WAAAAANQ4dbrl0sPDQwsWLNCdd96pvLw8jRo1Sm3atLG6LAAAAACocer0PZcAAAAAgMpRp7vFAgAAAAAqB+ESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4jXAJAAAAAHAa4RIAAAAA4DTCJQAAAADAaR5WF1CT9OzZUwEBAVaXAQAAAACWSE5O1t69e4t9jnBZAQEBAYqJibG6DAAAAACwRERERInP0S0WAIDSpKdLZ85IFy5YXQkAANUa4RIAgIvl5UkffijdfLN01VVSo0bS5ZdLPj5SSIi0caNkGFZXCQBAtUO3WAAA8n3wgTRvnvTtt1LTptKQIVLHjpKHh3TsmBk6w8Kk3r2ltWulVq2srhgAgGqDcAkAQFqaNH269NZb0rXXStHRUkSE5OVVcLucHOn116XZs6XrrzfDaN++lpQMAEB1Q7dYAEDd9vnnUqdO0jvvSI8+Kh06JI0bVzRYSpKnpzRpkvTVV1Lz5tLQodKuXa6vGQCAaohwCQCou95+W/rb36T69aW9e6UFC8wAWZbWraWPP5ZatJBGjJCOHq3yUgEAqO4IlwCAuumZZ8wWym7dzGDZtWvF9vf3lzZtktzdpfBw6ezZKikTAICagnAJAKh7Vq2SZs0y76uMj5eaNbu047RqZbZ+Hjki3XdfZVYIAECNQ7gEANQtq1dLM2ZII0dK69ZJ9eo5d7wBA8xgGRUlbd5cOTUCAFADES4BAHXHxo3Svfea81WuW1e++yvL49FHpQ4dpLvukjIyKueYAADUMIRLAEDd8NVX5j2WN9xgdmUtbjTYS1W/vtkievy4tHhx5R0XAIAahHAJAKj97HYpLEzy8THnpmzYsPLP0aePdNtt0vLl0g8/VP7xAQCo5giXAIDaLTfXbLFMSzODZfPmVXeuZcukBg2kmTMlw6i68wAAUA0RLgEAtdvDD0s7dpjdVrt0qdpz+ftLCxeaA/ts2VK15wIAoJohXAIAaq/YWLM18Z57pPHjXXPOqVOl1q2lBx6Q8vJcc04AAKoBwiUAoHY6dkyaMEHq1k16+mnXndfLyxzU5+BBae1a150XAACLES4BALXPuXPSqFGSh4e0fr3zc1lW1Nix0v/7f9Ijj0iZma49NwAAFiFcAgBqF8Mwu8EeOiS99ZbUsqXra3BzM7vj/vKL9Pzzrj8/AAAWIFwCAGqXqCjpjTfMgXUGD7aujgEDzPMvXiydOmVdHQAAuAjhEgBQe+zbJ82YIQ0dao4Sa7Vly8xguXSp1ZUAAFDlCJcAgNrhxAlp9GhzHss33jC7plqtSxfp73+XnnlGOn7c6moAAKhSLv+fd/r06dqxY4cuXLjg6lMDAGqr3Fxp3DgpJcUcwKdpU6sr+sNjj0kXLkiPPmp1JQAAVCmXh8ubb75ZH3zwgQYNGqTly5frxx9/LHOfXbt2afDgwQoODlZUVFSR5w3D0OOPP67g4GCFhYXp22+/LXPfTZs2KSQkRNddd50OHTpUOS8OAGCNmTOlbdukyEhz6pHqpFUracoU6ZVXpCNHrK4GAIAq4/Jw2bt3b61YsUL/+c9/FBAQoDvuuEPjxo3Te++9p5ycnCLb5+XladGiRXrppZcUFxenjRs36tixYwW22bVrlxITE7V161Y99thjWrhwYZn7tm3bVqtWrVL37t2r/DUDAKrQc89JL7wgzZ0rTZxodTXFmz9fatjQ/BcAgFrKkhtSTp48qZiYGL377rtq166dbrvtNn333Xe64447imx78OBBtWzZUi1atJCXl5dCQkIUHx9fYJv4+HiFh4fLZrPp+uuv15kzZ5SamlrqvldffbVat27tktcLAKgiW7aYrZbDh0tPPGF1NSW78krpvvukmBhp716rqwEAoEq4PFxOmzZNf//735WZmanIyEhFRkZq2LBheuSRR5SRkVFke7vdLn9/f8djPz8/2e32Urfx9/eX3W4v174AgBpq3z5p1CipUydp7VrJ3d3qiko3e7bk6ys9+KA5FycAALWMh6tPOGbMGPXr16/AuuzsbHl5eSkmJqbI9kYx/wHbbLZybVOefQEANdCRI+Z0I1deKW3aJHl7W11R2by9pUcekaZPN1tchwyxuiIAACqVy1sun3766SLrbrrpphK39/f3V0pKiuOx3W6Xr69vqdukpKTI19e3XPsCAGqY48elQYPMlsqPPjKnHqkpJk+WWrc2Wy8ZNR0AUMu4LFz+9ttvSkhIUGZmpr777jt9++23+vbbb7V3716dP3++xP06deqkxMREJSUlKTs7W3FxcQoKCiqwTVBQkGJjY2UYhr7++ms1atRIvr6+5doXAFCD/PqrNHCgdPq0tHmzdM01VldUMV5e5tQk33wjrVtndTUAAFQql3WL3b17t2JiYpSSkqInLhp0oWHDhpo9e3aJ+3l4eGjBggW68847lZeXp1GjRqlNmzaKjo6WZE5t0q9fP+3cuVPBwcFq0KCBlixZUuq+kvTRRx/pscceU1pamu6++261a9dOL7/8chV+BwAATklJkYKCzJbLTZukG26wuqJLM26c9NRT0sMPS6NHm4ETAIBawGYUd2NiFdqyZYsGDx7sylNWmoiIiGLvCwUAVLGUFKl/fykpyQyWfftaXZFzNm827xl99lnzHkwAAGqI0jKRy1ouN2zYoBEjRig5OVlr1qwp8vzE6jo3GQDAWvktlklJ0ocf1vxgKUmDB0t/+5vZRfb226VGjayuCAAAp7nsnsv8+yrPnTunjIyMIgsAAEUkJpph8r//NYNlYKDVFVUOm01aulT67Tdp5UqrqwEAoFK4rOVy3Lhxksx5LgEAKFNCgjkq7Pnz5qiwvXtbXVHl6tlTioiQli+Xpkwx58AEAKAGc/lUJE8++aTS09OVk5OjCRMmqGfPntqwYYOrywAAVGeffmq2WNps0ief1L5gmW/xYuncOfNfAABqOJeHyz179sjb21s7duyQv7+/tmzZwiitAIA/vP++Od3IlVdKe/ZIHTtaXVHVue466Y47pBdflH7+2epqAABwisvDZW5uriRp586dCgkJUZMmTVxdAgCgOjIM6cknpfBwqUMHafduqVUrq6uqegsXSu7u0oIFVlcCAIBTXB4u+/fvryFDhighIUG9evVSWlqa6tWr5+oyAADVSVaWNHGi9MAD0pgx0s6ddecexIAAacYMae1a6eBBq6sBAOCSuXyeS0k6ffq0vL295e7urvPnzys9PV1XXnmlq8uoMOa5BIAqkJJiBsrdu81WvAULzHst65KTJ6XWrc17S+PirK4GAIASVYt5Li/2448/Kjk5WXl5eY514eHhVpQCALDSJ59IN90knTolrVtnfl0XXXGF9OCD5rJrV+2ZcgUAUKe4PFzOnTtXSUlJuu666+Tu7i5JstlshEsAqEsMQ1qxwgxTrVtLW7ZInTpZXZW1ZsyQVq2SZs+WvvhCcnP5nSsAADjF5eEyISFBH374oWx1rcsTAMB0+rQ5QmpMjDnP45o1UuPGVldlvQYNzAGN/v5383syaZLVFQEAUCEu/1i0TZs2+u2331x9WgBAdfDNN1K3btKGDWbL5fr1BMuL3Xyz1KePNG+e2VUYAIAaxOUtlydPnlRISIg6d+4sT09Px/rIyEhXlwIAcJULF6RnnjG7wTZtKn38sdS3r9VVVT82m/Tss2YAf/RR6V//sroiAADKzeXhcvr06a4+JQDASr/+Kt1+u7R1qzR8uPTyy1KzZlZXVX117SpNnmyGzL//3QyaAADUAC7vFtujRw8FBAQoNzdXPXr0UKdOndS+fXtXlwEAcIUPPpA6dzZHhX3xRSk2lmBZHkuXSv7+5n2X2dlWVwMAQLm4PFy+8847mjFjhhYsWCBJstvtmjp1qqvLAABUpbQ0acIEs6XyqqukAweke+6pe/NXXqomTcwwfvCgtGyZ1dUAAFAuLg+Xa9euVXR0tLy9vSVJrVq1UlpamqvLAABUlZgYqX176a23pIcflj7/XGrXzuqqap7hw6Vx46THHpO+/NLqagAAKJPLw6WXl5e8vLwcj3Nzc11dAgCgKtjt0tix0qhR0p/+JO3bZwajevWsrqzmWrVK8vOTxowxp3ABAKAac3m47N69uyIjI5WZmak9e/Zo5syZCgoKcnUZAIDKkp0trVwptW0rvf++tGSJtHevdP31VldW8zVrJr39tvTLL+bcoIZhdUUAAJTI5eHyvvvuk4+Pj9q2bau3335b/fr106xZs1xdBgCgMmzaZA7YM2eOOT/jN9+YczReNNUUnNS7t3nfZUyM9MQTVlcDAECJXD4ViZubmwYOHKiBAwfKx8fH1acHAFSG3bulhQul+HipTRspLk4aNszqqmqvf/zDHBRp/nypeXNp4kSrKwIAoAiXhUvDMPTcc8/pzTffdDx2c3PTrbfeqmnTprmqDACAMy4Olb6+ZnfYqVOli+6lRxWw2aQ1a6TffpPuuktq2tQc8AcAgGrEZd1iX3vtNX355Zdav3699u7dqy+++ELvvvuuvvrqK7366quuKgMAUFE5OdK6dVKvXlLfvtKhQ9KKFdLPP5stagRL1/Dykt57T+ra1Rw06f8+rAUAoLpwWbiMjY3VihUr1KJFC8e6Fi1a6KmnnlJsbKyrygAAlNfvv5uD8/zlL9LNN5uPn33WDJWzZ0uXXWZ1hXVPo0bSRx+ZIX/8eGn5cgb5AQBUGy7rFpubm1vsPZY+Pj5MRwIA1YVhmF1fX37ZHKU0M1MaOFBavVoaOlRyc/k4cCjs8svNgZRuu02aO9ecR/Tf/5auuMLqygAAdZzL3iV4ljJyYGnPAQBcICXFHJH0uuukwEBzZNIJE6SEBLOlLCSEYFmd1KsnRUdLTz4pbdggtWtndpO9cMHqygAAdZjLWi6PHDmirl27FllvGIays7NdVQYAIF9KivSf/0jr10s7dpjBpG9f6aGHpNGjpYYNra4QpXFzM1suBw6U7rnH7Cb7xBPmfbBjx0qNG1tdIQCgjnFZuDx8+LCrTgUAKE5OjrRvn7Rtm9kauWeP2Q322mvNQDl+vNS2rdVVoqJuuEH69FPzQ4LHHzdHk502zRyAqX9/qUMH6eqrJX9/ydvbvFfWzc382eflSbm55pKRIZ09K6Wnl7xkZJRvycoyz+HhIbm7m/96e5uj3DZtKgUEmFPYXHON1LGj1KKFOSIuAKBGc/k8lwAAFzAM6fhxc27E/fvNZc8eMyDYbOaIowsWSGPGSO3b88a+pnN3l266yWyx3LtXevdd6eOPzWljCg/4Y7OZwS8vr+LncXMzW7QbNjRDav7XDRtKzZr98XW9emZL+MXh9exZ6cQJ6YcfzA84zp7947h+flL37lKfPtKgQdL119MNGwBqIMIlAFjNMMw339nZ5r95eeVbcnKktDRzFNfff5d+/VU6duyP5cwZ8/ju7mbr0K23ml0o+/eXihlgDbWAzSbdeKO5SOY1cOyY9OOPZrA7e9ZcLlwo2KqY37JYeGnU6I+vGzaU6tevnA8iDMOcs/PoUembb6QvvjBD8caN0rx5ZlANDpZCQ82FLr4AUCMQLgGgMuTmmm/ijx41WwyTk803z/lv5vOX9HRzBNbMTLPrYP5SGdNJuLub04Zcc43Uu7fZ3bVbN6lLF6lBA+ePj5qncWOzlbqYMQ8sZbNJvr7m0qePdO+95vqUFLNVc+tWc4mONuf3HDTIvA94+HBGxQWAaoxwCQCXIidH2rVLio+Xtm83W18yM/943s3NbH1p3Nhs/WnUSPrTn8wWoPr1zW6DhRcvrz9aktzdy148PMwWyCuvNM/l42OuB2oqf3+zhf3WW83W1c8/N+8lfe89s1XT01MKC5MmTpSGDDF/BwAA1QZ/lQGgIr76SoqKMt/w/v67+ea2Rw9p6lSpc2dzKo+rrjLvISPoAZfOzc1sge/dW1qxwrxvODranHIlJsYMouPHS3fcYf7eAQAsR7gEgLJcuCBt3iwtX24OktKggdk976abzHsYGzWyukKgdrPZzAF/unc352ONi5NefVX617+kp56S/vpXc5Tc0aPNgYYAAJZgKDYAKElmpvTyy+ZgOCEh5iiXTz4p/e9/0rp10siRBEvA1Tw9pfBwKTbWvL952TLJbpcmTDC7nk+bZnZTBwC4HOESAAo7ccKcL7BVK+nOO837Id94Q/rpJ3PS+iZNrK4QgGR2P7//fun7781eBSEh0ksvmVOZ9Ogh/fvfBac8AQBUKcIlAOT77jvpnnvMCd0fecScnH7bNunLL80BRry8rK4QQHFsNulvf5PWrjV7Fjz9tHTunDR5stmaOXGitGGDuQ4AUGUIlwDqtpwc6YMPzKkOOnSQXntNuuUW6dAhadMmacCAypnXD4Br+PhIM2eav8OffiqNGSP95z9mV9qmTc3WzaeeMufVzMmxuloAqFUY0AdA3WMY5qTtb75p3jv5++9m68bixWZLR7NmVlcIwFk2m9Srl7lERppTB23YIH30kfThh+Y2DRua82z26mXOCdu9u9nVFgBwSWpEuNy1a5cWL16sCxcuaMyYMZo8eXKB5w3D0OLFi7Vz507Vr19fS5cuVYcOHUrd99SpU/rHP/6h5ORkBQQE6Omnn9bll1/u8tcGwEUyM6UdO8xRJjdulBITzXspR4wwu7wOGWIOFAKg9vHyMkd2HjjQfJySIn3yibRzpxk6H3vMHBVaMqcS6tbtj6VjR/PDJ3owAECZqn24zMvL06JFi7RmzRr5+flp9OjRCgoK0jXXXOPYZteuXUpMTNTWrVv1zTffaOHChXr33XdL3TcqKkq9evXS5MmTFRUVpaioKM2dO9fCVwqgUuXlmd3idu2S4uPNeyfPnTOnERk40LynctQoiQ+VgLrH39/sLjtmjPk4PV36+mtp3z5zPs39+83RaPM1amTOpdmunfnvNddIf/6zufj5mXNyAgCqf7g8ePCgWrZsqRYtWkiSQkJCFB8fXyBcxsfHKzw8XDabTddff73OnDmj1NRUJScnl7hvfHy83njjDUlSeHi4xo8fT7gEaqqsLOm//5UOHjTfIH71lbRnj3T6tPl8q1bS7bdLoaHmoB8NGlhXK4Dqx9vbnCvzr3/9Y92pU+bfku++k44ckQ4fNj+kev31gvt6epqtnVddJV15pdmt/uKlaVOpcWMzoDZq9MfX9eq59CUCgCtU+3Bpt9vl7+/veOzn56eDBw+Wuo2/v7/sdnup+544cUK+vr6SJF9fX6WlpVXly6g6589L0dGlj4BnGKUfw+rnq0MN1f356lBDdXg+M9MMjKdPS2fOSCdPSr/8YnZxy9/f3d1sWRg7VgoMNJc//7n0YwNAYU2aSP37m8vFzpwxu9X/8ou5JCWZ/yYnm1Oi7Nlj3sedl1f68T09iw+d9eubzxW3uLmZ3XMLLyWtd2bbkrYD4DrXXisFB1tdRYVU+3BpFPOG01boj1tJ25Rn3xrvm2+ku+76414RoCYr6/ezXj2zG2v+0qSJNHSo1LKl2TrZsaPUvr355gwAqkLjxlLnzuZSEsMwPwT7/Xdz3twzZ8z5Ns+e/ePr4tadOGF+iJaTU/xy4YJ57OKWkp4DUHNdc4109KjVVVRItQ+X/v7+SklJcTy22+2OFseStklJSZGvr69ycnJK3Ldp06ZKTU2Vr6+vUlNT5ePjU8WvpIrceKOUllb2cOplvWm3+vnqUEN1f7461FAZrwEAajubzfzwq0kT882h1cobRMvaDoBreXtbXUGFVftw2alTJyUmJiopKUl+fn6Ki4vTihUrCmwTFBSkN998UyEhIfrmm2/UqFEj+fr6ysfHp8R9g4KCFBsbq8mTJys2NlYDBgyw4uVVDgYkAQAAJaFLKwAXqfbh0sPDQwsWLNCdd96pvLw8jRo1Sm3atFF0dLQk6eabb1a/fv20c+dOBQcHq0GDBlqyZEmp+0rS5MmTNWvWLK1fv17NmzfXM888Y9lrBAAAAICazmYUd2MiitWzZ08FBARYXQYAAAAAWCI5OVl79+4t9jnCJQAAAADAacz6CwAAAABwGuESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4jXAJAAAAAHAa4RIAAAAA4DTCJQAAAADAaYRLAAAAAIDTCJcAAAAAAKcRLgEAAAAATvOwuoCapGfPngoICLC6DAAAAACwRHJysvbu3Vvsc4TLCggICFBMTIzVZQAAAACAJSIiIkp8jm6xAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4jXAJAAAAAHAa4RIAAAAA4DTCJQAAAADAaYRLAAAAAIDTCJcAAAAAAKcRLgEAAAAATiNcAgAAAACcRrgEAAAAADiNcAkAAAAAcBrhEgAAAADgNMIlAAAAAMBphEsAAAAAgNMIlwAAAAAApxEuAQAAAABOI1wCAAAAAJxGuAQAAAAAOI1wCQAAAABwGuESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4jXAJAAAAAHAa4RIAAAAA4DTCJQAAAADAaYRLAAAAAIDTPKwuAACAAnx8pJMnra6ial1xhZSWZnUVAABUKsIlAKB6OXlSMgxrzm2zuebcNlvVnwMAABejWywAAAAAwGmlhsszZ85o7dq1pR7g+PHj+uCDD8o80fHjxxUaGlqx6ipReV4LgJph8ODBcnNzk81mq1YLANey+ne+KpdL+RvXtGlTTZ8+XR07dpS7u7tatGihFi1ayM3NTZ6eno7tPD095ebmpo4dOyo6OlqSFB0d7djv4vUVNX36dNWvX182m03169fX9OnTnfoZV1ZdAFyjzHBZ1i9xcnKyNm7cWKlFVYXyvBYA1d/gwYO1detWBQQEWF1KEQRMwHVq+++bcQnds9PS0vT8888rPDxcr776qnJzc5Weni5vb281bty4wNKsWTOFh4dr/vz5mj59uubPn69Vq1YpMzNTq1at0vz58yv8vmn69OmKjIzUkiVLlJGRoSVLligyMvKSA2Z0dHSl1AXAhYxSzJo1y+jUqZMxfPhwY+nSpcbSpUuNkJAQIzQ01IiLizMMwzDGjBljdO3a1Rg+fLixZs0aIykpybj55puN8PBwIzw83Dhw4IBhGIaRlJRkhISElHiu3NxcY+nSpUZoaKgRGhpqvP7664ZhGMann35qjBgxwggNDTUefPBBIysryzAMw+jfv79x4sQJwzAM4+DBg8att95qGIZhPPvss8aDDz5o3HrrrUZQUJDx2muvFfta7Ha7ccsttxjDhw83QkJCjH379pX2rTAMwzBGjhxZ5jYAqpbNZjOmTJli2Gw2Q1K1W1AJrPw+uurcXCtOs/p3vaoWd3f3Cm3v6elpTJkypcD+HTp0MDp06GBs377daNWqleHp6en4ulWrVo6v87epV6+esX379gLf3+3btxsdOnSo0M+kXr16xooVKwqsW7FihVGvXr1L+hnn1+dsXQAqV2mZqNT/3S4OhJs3bzZuv/12Izc31/jtt9+Mfv36GXa73fj888+NyZMnO/Y5d+6ckZmZaRiGYfz888+Ok5cVLteuXWtMmzbNyMnJMQzDME6ePGlkZmYagYGBxk8//WQYhmHMnTvXWLNmjWEYpYfLm266ycjKyjJOnDhh9OjRw8jOzi5y/pdfftl44YUXDMMwg+3Zs2dL+1YYhkG4BKoDScapU6csfwNYYrhkqZzFugvMdedhcWqx+ne9Oi2F/ya6ubkZbm5uRnZ2tuODuPyv89dfvE3+8xfL36Zil7WMjIyMAusyMjIM6dJ+r/Lrc7YuAJWrtExU7gF9Dhw4oJCQELm7u6tZs2bq3r27Dh06VGS73NxcPfzwwwoLC9PMmTP1448/luv4n332mcaNGycPD3MA2yZNmujnn3/WVVddpb/85S+SpJEjR2r//v1lHqtfv37y8vKSj4+PfHx8dOLEiSLbdOrUSTExMVq1apV++OEHeXt7l6tOANay2WyaN29e9e0SZ/lb7lqw1BVWf59r+lJLubu7V2h7T09PzZs3r8D+7dq1U7t27bR79261bNlSnp6ejq///Oc/a/fu3frzn//s2KZevXravXt3gePu3r1b7dq1q1At9erVU2RkZIF1kZGRqlevXoWOky+/PmfrAuA65Q6XRjn/kL/66qtq1qyZNmzYoPfee085OTnlPn7hN4ulndPd3d3xfFZWVoHnvLy8CmyXm5tbZP/u3bvrzTfflJ+fn+6//37FxsaWq04A1goODtaLL75YLe+5BABn5eXlVWj7nJwcvfjii47HFy5cUHh4uB544AHdcsstOnXqlOrXr6/Ro0crLS1NaWlpGj16tDIyMhQeHq5Jkybprrvu0qRJk/Txxx8rJydHH3/8sSZNmqT58+dXqJa77rpLDzzwgFauXKlz585p5cqVeuCBB3TXXXdV6Dj55s+fXyl1AXCdUue5bNiwoTIyMiSZYeztt9/WyJEjdfr0ae3fv1/333+/7Ha7YxtJOnv2rPz9/eXm5qb//Oc/5f4j2adPH61bt049evSQh4eHTp06pdatWys5OVn//e9/1bJlS23YsEHdu3eXJAUEBCghIUH9+vXT1q1byzz+xa9FMgci8vPz09ixY3Xu3Dl9++23Cg8PL1etAKyzZcsWDR48WB999JHVpRRR3g/hADivuA+laxObzVbhvyk+Pj665ZZbFBsbq8OHD+tPf/qTvL29lZycXODDdg8PD+Xl5Sk2NlaLFy/WzTffrN69e2v69Ok6fPiw2rVr51hfEatWrZIkPfTQQ5ozZ47q1aune+65x7G+ovLP72xdAFyn1HB5xRVXqGvXrgoNDVXfvn3Vtm1bjRgxQjabTXPnztWVV16pJk2ayN3dXcOHD1dERIRuueUWTZ8+XZs3b1bPnj112WWXlauQMWPGKDExUcOHD5eHh4fGjh2rW2+9VU888YRmzpypvLw8dezY0fEHZdq0aZo/f75Wr16tLl26lHn84l7Lyy+/LA8PD1122WVatmxZueoEYL0tW7ZYXQKAaoAPdCrPzTffXCmhbdWqVZccJotTWXUBcA2bwV/mcouIiFBMTIzVZQBA7WazWXdPnavObeVrBADACaVlolJbLgEAsISV3R1dce4rrqj6cwAA4GIuD5effPKJli9fXmDdVVddpeeff97VpQAAqiNa9AAAqJFcHi779u2rvn37uvq0AAAAAIAqVO6pSAAAAAAAKAnhEgAAAADgNMIlAAAAAMBphEsAAAAAgNMIlwAAAAAApxEuAQAAAABOI1wCAAAAAJxGuAQAAAAAOI1wCQAAAABwGuESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE4jXAIAAAAAnEa4BAAAAAA4jXAJAAAAAHAa4RIAAAAA4DTCJQAAAADAaYRLAAAAAIDTCJcAAAAAAKcRLgEAAAAATiNcAgAAAACcRrgEAAAAADiNcAkAAAAAcBrhEgAAAADgNMIlAAAAAMBphEsAAAAAgNMIlwAAAAAApxEuAQAAAABOI1wCAAAAAJzmYXUBNUlycrIiIiKsLgMAAAAALJGcnFziczbDMAwX1gIAAAAAqIXoFgsAAAAAcBrhEgAAAADgNMIlAAAAAMBphEsAAAAAgNMIlwAAAAAApxEugXLatWuXBg8erODgYEVFRVldDlxs3rx56tWrl0JDQx3rTp06pYkTJ2rQoEGaOHGiTp8+7Xhu9erVCg4O1uDBg/XJJ5841ickJCgsLEzBwcF6/PHHxYDdtc+vv/6q8ePHa+jQoQoJCdFrr70miesFBWVlZWn06NEaPny4QkJC9Oyzz0riOkHJ8vLyFB4errvvvlsS1wqKFxQUpLCwMI0YMcIxhaJLrxUDQJlyc3ONAQMGGL/88ouRlZVlhIWFGUePHrW6LLjQF198YSQkJBghISGOdcuWLTNWr15tGIZhrF692njyyScNwzCMo0ePGmFhYUZWVpbxyy+/GAMGDDByc3MNwzCMUaNGGV9++aVx4cIFY9KkScaOHTtc/2JQpex2u5GQkGAYhmGcPXvWGDRokHH06FGuFxRw4cIFIz093TAMw8jOzjZGjx5tfPXVV1wnKNErr7xizJ4925g8ebJhGPwfhOL179/fOHHiRIF1rrxWaLkEyuHgwYNq2bKlWrRoIS8vL4WEhCg+Pt7qsuBC3bt31+WXX15gXXx8vMLDwyVJ4eHh2rZtm2N9SEiIvLy81KJFC7Vs2VIHDx5Uamqq0tPTdcMNN8hmsyk8PJzrqBby9fVVhw4dJEne3t5q3bq17HY71wsKsNlsatiwoSQpNzdXubm5stlsXCcoVkpKinbs2KHRo0c71nGtoLxcea0QLoFysNvt8vf3dzz28/OT3W63sCJUBydOnJCvr68kM1CkpaVJKvl6Kbze39+f66iWO378uA4fPqwuXbpwvaCIvLw8jRgxQr1791bv3r25TlCiJUuWaO7cuXJz++OtO9cKSjJp0iRFRETo7bffluTaa8Wjsl4EUJsZxfQzt9lsFlSCmqCk64XrqG7JyMjQjBkz9NBDD8nb27vE7bhe6i53d3dt2LBBZ86c0dSpU/XDDz+UuC3XSd318ccfy8fHRx07dtTevXvL3J5rpW6Ljo6Wn5+fTpw4oYkTJ6p169YlblsV1wrhEigHf39/paSkOB7b7XbHJ0Cou5o2barU1FT5+voqNTVVPj4+kkq+XgqvT0lJ4TqqpXJycjRjxgyFhYVp0KBBkrheULLGjRurZ8+e+uSTT7hOUMSXX36p7du3a9euXcrKylJ6erruu+8+rhUUy8/PT5L5f05wcLAOHjzo0muFbrFAOXTq1EmJiYlKSkpSdna24uLiFBQUZHVZsFhQUJBiY2MlSbGxsRowYIBjfVxcnLKzs5WUlKTExER17txZvr6+atiwob7++msZhlFgH9QehmFo/vz5at26tSZOnOhYz/WCi6WlpenMmTOSpMzMTH366adq3bo11wmKmDNnjnbt2qXt27dr5cqVuvHGG7V8+XKuFRRx7tw5paenO77es2eP2rRp49JrhZZLoBw8PDy0YMEC3XnnncrLy9OoUaPUpk0bq8uCC82ePVtffPGFTp48qcDAQE2fPl2TJ0/WrFmztH79ejVv3lzPPPOMJKlNmzYaOnSohg0bJnd3dy1YsEDu7u6SpIULF2revHnKzMxUYGCgAgMDrXxZqAIHDhzQhg0b1LZtW40YMUKSef1wveBiqampevDBB5WXlyfDMDRkyBD1799f119/PdcJyoW/KSjsxIkTmjp1qiTznu7Q0FAFBgaqU6dOLrtWbEZxnWoBAAAAAKgAusUCAAAAAJxGuAQAAAAAOI1wCQAAAABwGuESAAAAAOA0wiUAAAAAwGmESwAAAACA0wiXAACnXHvttZo7d67jcW5urm688Ubdfffdpe4XExOjRYsWVWotq1at0ssvv1ypx3TG8ePH9cEHH1Ron8jIyDK3OXPmjNauXVuu491www0VOn9lK8/rqS4uviar27UEADUB4RIA4JTLLrtMR48eVWZmpiRpz5498vPzs7iq6iE5OVkbN26s0D6rV68uc5szZ84oOjr6UstyqfK8nsqSm5vrsnMBAIoiXAIAnBYYGKgdO3ZIkuLi4hQSEuJ47tSpU7r33nsVFhamsWPH6siRI0X2T0tL0/Tp0zVq1CiNGjVKBw4ckCRlZGRo3rx5CgsLU1hYmLZs2SKpYGvc5s2b9eCDDxY55i+//KJJkyYpIiJCt9xyi3788cci25R0/I0bNyosLEyhoaF66qmnHNuXdN4HH3xQjz/+uMaNG6cBAwZo8+bNkqQVK1Zo//79GjFihF599VUdPXpUo0eP1ogRIxQWFqbExMQC9SxfvlyZmZkaMWKE5syZI0las2aNQkNDFRoaqldffdVx3F9++UUjRozQsmXLlJGRoQkTJmjkyJEKCwvTtm3bivkpFe/f//63wsLCNHz4cC1fvlySdPjwYY0dO1ZhYWGaOnWqTp8+LUkaP368Dh06JMn8mQUFBUkyW/ymTZumSZMmadCgQXryySeLfT3nzp3T5MmTNXz4cIWGhurDDz8sUs/48eO1ePFijRs3TqGhoTp48KAk6dy5c5o3b55GjRql8PBwx2uMiYnRjBkzdM899+iOO+4ocrzY2FjH68tvYS/pegMAOMfD6gIAADXfsGHD9MILL6h///76/vvvC7xhX7Vqldq3b68XXnhBn332mR544AFt2LChwP6LFy/WhAkT1K1bN/3vf//TpEmTtGnTJr3wwgvy9vZ2dC3NDznl8cgjj+jRRx9Vq1at9M033+jRRx/V66+/XmCb4o5vt9u1fPlyxcTEqHHjxrrjjju0bds2DRw4sNTzpaam6q233tJPP/2kKVOmaMiQIZozZ45eeeUVR+vdY489pttuu03Dhw9Xdna2Lly4UOAY9913n9auXev4/iQkJCgmJkbvvPOODMPQ2LFj1aNHD82ZM0dHjx51bJebm6vnn39e3t7eSktL00033aQBAwbIZrOVWvPOnTsVHx+vd955Rw0aNNCpU6ckSffff78eeeQR9ejRQ88884yee+45zZ8/v9RjHT58WLGxsfLy8tKQIUM0fvz4Iq9ny5Yt8vX1VVRUlCTp7NmzxR7r/PnzWrdunfbt26eHHnpIGzduVGRkpG688UY98cQTOnPmjMaMGaPevXtLkr7++mu9//77atKkSYHjHD16VC+++KKio6Pl4+PjeH0lXW8AAOcQLgEATrvuuut0/Phxbdy4Uf369Svw3IEDB7Rq1SpJUq9evXTq1KkioeLTTz/VsWPHHI/T09OVnp6uzz77TCtXrnSsv/zyy8tVT0ZGhr766ivNnDnTsS47O7vIdsUdf9++ferRo4d8fHwkSWFhYdq3b1+Z4XLgwIFyc3PTNddco99//73Yba6//npFRkYqJSVFgwYNUqtWrUo95oEDBzRw4EBddtllkqTg4GDt37/f0WKYzzAMrVy5Uvv27ZObm5vsdrt+//13XXnllaUe/7PPPlNERIQaNGggSWrSpInOnj2rs2fPqkePHpKkkSNHFvg+lqRXr15q1KiRJOnqq69WcnKymjdvXmCbtm3batmyZXrqqafUv39/devWrdhj5bd8d+/eXenp6Tpz5ox2796t7du365VXXpEkZWVl6ddff5Uk9enTp0iwlKTPP/9cQ4YMcfws87cp6XoDADiHcAkAqBRBQUF68skn9frrrztaiCQz+BRWuEXtwoULevvtt1W/fv0C6w3DKLP1LSsrq8g6wzDUuHHjIi2kxW1X1vHLe14vL68y9wkLC1OXLl20Y8cOTZo0SY8//rh69epVan3l8cEHHygtLU0xMTHy9PRUUFBQsd+X4o5fkdfv7u7uqKlwWL/49bu7uysvL6/I/n/5y18UExOjnTt3asWKFerTp4+mTZtWZLvCNeU/fvbZZ9W6desCz33zzTeOcFxYSd+/kq43AIBzuOcSAFApRo8erXvvvVfXXnttgfXdu3fX+++/L0nau3evrrjiCnl7exfY5q9//avefPNNx+PDhw9LMlukLl6f3y22WbNm+vHHH3XhwoVi7y/09vbWVVdd5ejqaBhGsfd6Fnf8zp07a9++fUpLS1NeXp7i4uLUvXv3cp23sIYNGyojI8PxOCkpSS1atNBtt92moKAgff/990X28fDwUE5OjiTze7dt2zadP39e586d07Zt29StW7cixz179qyaNm0qT09Pff7550pOTi6ztvzX/9577+n8+fOSzPtjGzVqpMaNG2v//v2SpA0bNjhef0BAgBISEiTJcV9pWS5+PXa7XQ0aNNCIESM0adIkfffdd8Xuk38v5v79+9WoUSM1atTIcY3kB8aS9r1Yr169tHnzZp08edLx+qSSrzcAgHNouQQAVAp/f39NmDChyPpp06Y5Bs1p0KCBli5dWmSb+fPna9GiRQoLC1NeXp66deumRYsWacqUKVq0aJFCQ0Pl5uamadOmadCgQZozZ47uvvtuNW/eXG3atNG5c+eKHPOpp57SwoUL9eKLLyo3N1fDhg3TddddV2Cbko4/e/ZsTZgwQYZhKDAw0NEltjznvdi1114rd3d3DR8+XBEREcrKytL7778vDw8PNWvWTFOnTi2yz9ixYzV8+HC1b99eK1asUEREhMaMGSPJDPDt27eXJHXt2lWhoaHq27ev7rrrLk2ZMkURERFq165dkda9kgQGBurIkSMaNWqUPD091a9fP82ePVvLli3TP//5T50/f14tWrTQE088IUm64447NGvWLL3//vvq2bNnuc5x8esJDw/Xk08+KTc3N3l4eGjhwoXF7nP55Zdr3LhxSk9P15IlSyRJ9957r5YsWaLhw4fLMAwFBASUORJtmzZtdM8992j8+PFyc3NT+/bttXTp0hKvNwCAc2xGefvcAAAAVLHx48fr/vvvV6dOnawuBQBQQXSLBQAAAAA4jZZLAABque+//173339/gXVeXl569913LaoIAFAbES4BAAAAAE6jWywAAAAAwGmESwAAAACA0wiXAAAAAACnES4BAAAAAE4jXAIAAAAAnPb/Ae8HxXsE8Ei9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "total_molecules_per_cell(a, count=\"total_counts\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "based-camping",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "saved to: /data/peer/wallet/Jupyter/2021_SPADE/2021_Bassez_BRCA/data/BRCA-X-TIL-X-Bassez_2021-X-cohort1_2_raw_filtered_clustered_drops_annotated_nodrops_log1p_clustered_leukocytes_scran_v3_210501.h5ad\n"
     ]
    }
   ],
   "source": [
    "### Save\n",
    "PATH = #h5ad path\n",
    "FILE = #h5ad path\n",
    "adata1scr.write(PATH + FILE)\n",
    "print('saved to:', PATH + FILE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pressed-radius",
   "metadata": {},
   "source": [
    "# End"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "experimental-personal",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dressed-workstation",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}