LARRY_for_10X.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c669e9ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gzip, time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import defaultdict\n",
    "\n",
    "def is_valid(bc):\n",
    "    return (bc[4:6]=='TG'\n",
    "            and bc[10:12]=='CA' \n",
    "            and bc[16:18]=='AC' \n",
    "            and bc[22:24]=='GA' \n",
    "            and bc[28:30]=='GT'\n",
    "            and bc[34:36]=='AG')\n",
    "\n",
    "def in_filtered_list(cell_bc, filtered_cell_barcodes):\n",
    "    num_Ns = sum([c=='N' for c in cell_bc])\n",
    "    if num_Ns > 1: return False\n",
    "    elif num_Ns == 1: return np.any([cell_bc.replace('N',c) in filtered_cell_barcodes for c in 'ACTG'])\n",
    "    else: return cell_bc in filtered_cell_barcodes\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2736e0ad",
   "metadata": {},
   "source": [
    "## Set input params\n",
    "\n",
    "N_READS: Cutoff for the min number of reads supporting a given UMI/GFP-BC pair in a particular cell\n",
    "\n",
    "N_UMIS: Cutoff for the min number of UMIs supporting a given cell/GFP-BC pair\n",
    "\n",
    "N_HAMMING: Cutoff for the max hamming distance to collapse two GFP-BCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "0114a73f",
   "metadata": {},
   "outputs": [],
   "source": [
    "larry_prefix = 'GTTGCTAGGAGAGACCATATG'\n",
    "N_READS = 10\n",
    "N_UMIS = 3\n",
    "N_HAMMING = 3\n",
    "\n",
    "output_prefix = 'Larry_day-0'\n",
    "\n",
    "sample_paths = {\n",
    "    'Larry_day-0': {\n",
    "        'filtered_cell_barcodes_path':'barcodes/barcodes.tsv.gz',\n",
    "        'larry_barcode_fastq_paths':['barcodes/Larry_system_S7_L001_{}_001.fastq.gz',\n",
    "            'barcodes/Larry_system_S6_L001_{}_001.fastq.gz',\n",
    "            'barcodes/Larry_system_S8_L001_{}_001.fastq.gz',\n",
    "            'barcodes/Larry_system_S5_L001_{}_001.fastq.gz'\n",
    "        ]\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e3d2694",
   "metadata": {},
   "source": [
    "### Count all (sample, cell-BC, UMI, LARYY-BC) combinations with a valid cell-BC and LARRY-BC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0684cde6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 1000000 lines in 14.133716106414795 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 2000000 lines in 28.340556859970093 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 3000000 lines in 41.80684494972229 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 4000000 lines in 56.10518288612366 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 5000000 lines in 70.00225687026978 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 6000000 lines in 84.42418193817139 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 7000000 lines in 98.42103576660156 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 8000000 lines in 112.3191339969635 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 9000000 lines in 126.36743497848511 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 10000000 lines in 140.66051578521729 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 11000000 lines in 154.63344597816467 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 12000000 lines in 168.5036759376526 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 13000000 lines in 183.1948699951172 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 14000000 lines in 197.05925583839417 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 15000000 lines in 210.66526794433594 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 16000000 lines in 224.2083830833435 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 17000000 lines in 237.74522995948792 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 18000000 lines in 251.82353687286377 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 19000000 lines in 266.00370478630066 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 20000000 lines in 280.11109471321106 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 21000000 lines in 293.83448600769043 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 22000000 lines in 307.7327370643616 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 23000000 lines in 321.73655581474304 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 24000000 lines in 335.1471040248871 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 25000000 lines in 348.52025079727173 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 26000000 lines in 362.19200897216797 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 27000000 lines in 375.7206280231476 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 28000000 lines in 389.0827898979187 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 29000000 lines in 402.4177248477936 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 30000000 lines in 415.7678940296173 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 31000000 lines in 429.09841895103455 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 32000000 lines in 442.69291710853577 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 33000000 lines in 456.2546420097351 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 34000000 lines in 469.8486239910126 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 35000000 lines in 483.20930004119873 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 36000000 lines in 497.09227991104126 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 37000000 lines in 511.26597785949707 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 38000000 lines in 525.0689489841461 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 39000000 lines in 538.7632620334625 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 40000000 lines in 552.4738030433655 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 41000000 lines in 566.0128998756409 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 42000000 lines in 579.754903793335 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 43000000 lines in 593.2382838726044 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 44000000 lines in 606.7555840015411 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 45000000 lines in 620.5737419128418 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 46000000 lines in 634.0089199542999 seconds\n",
      "barcodes/Larry_system_S7_L001_{}_001.fastq.gz: Processed 47000000 lines in 647.7209160327911 seconds\n",
      "ERROR extracting barcodes/Larry_system_S7_L001_{}_001.fastq.gz\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 1000000 lines in 13.88553500175476 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 2000000 lines in 27.69385600090027 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 3000000 lines in 41.50436305999756 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 4000000 lines in 53.67356896400452 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 5000000 lines in 67.31585311889648 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 6000000 lines in 81.04626774787903 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 7000000 lines in 94.7165687084198 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 8000000 lines in 108.51787495613098 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 9000000 lines in 122.1899938583374 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 10000000 lines in 136.52827286720276 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 11000000 lines in 150.21489000320435 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 12000000 lines in 164.1510248184204 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 13000000 lines in 177.9083650112152 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 14000000 lines in 192.69267773628235 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 15000000 lines in 206.40870785713196 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 16000000 lines in 220.16033291816711 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 17000000 lines in 233.85000872612 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 18000000 lines in 247.57693696022034 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 19000000 lines in 261.41951179504395 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 20000000 lines in 275.90805101394653 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 21000000 lines in 289.8791286945343 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 22000000 lines in 303.8281989097595 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 23000000 lines in 318.7013359069824 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 24000000 lines in 332.69948077201843 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 25000000 lines in 346.7188878059387 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 26000000 lines in 360.87361693382263 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 27000000 lines in 376.67424988746643 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 28000000 lines in 391.59569478034973 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 29000000 lines in 407.68695974349976 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 30000000 lines in 424.7530617713928 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 31000000 lines in 443.1041269302368 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 32000000 lines in 460.38357496261597 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 33000000 lines in 478.8512759208679 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 34000000 lines in 496.6092767715454 seconds\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 35000000 lines in 512.647481918335 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 36000000 lines in 529.0447018146515 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 37000000 lines in 545.037318944931 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 38000000 lines in 565.0010149478912 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 39000000 lines in 584.8846368789673 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 40000000 lines in 601.7119328975677 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 41000000 lines in 617.5741739273071 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 42000000 lines in 633.5346648693085 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 43000000 lines in 649.2436039447784 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 44000000 lines in 665.4648399353027 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 45000000 lines in 681.6242208480835 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 46000000 lines in 697.1213748455048 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 47000000 lines in 714.4043698310852 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 48000000 lines in 730.1609930992126 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 49000000 lines in 746.1193459033966 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 50000000 lines in 761.6569538116455 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 51000000 lines in 777.5359599590302 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 52000000 lines in 805.7587578296661 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 53000000 lines in 823.692549943924 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 54000000 lines in 840.0294888019562 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 55000000 lines in 854.6428778171539 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 56000000 lines in 869.2348408699036 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 57000000 lines in 883.8786189556122 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 58000000 lines in 898.3419768810272 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 59000000 lines in 912.593948841095 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 60000000 lines in 927.2461559772491 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 61000000 lines in 941.8562350273132 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 62000000 lines in 956.1678850650787 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 63000000 lines in 970.8963198661804 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 64000000 lines in 984.9767999649048 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 65000000 lines in 998.7180027961731 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 66000000 lines in 1012.5251970291138 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 67000000 lines in 1026.34472990036 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 68000000 lines in 1040.2899868488312 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 69000000 lines in 1054.3572187423706 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 70000000 lines in 1068.3529028892517 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 71000000 lines in 1082.4897179603577 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 72000000 lines in 1096.944359779358 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 73000000 lines in 1111.0051720142365 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 74000000 lines in 1124.98140001297 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 75000000 lines in 1138.8252680301666 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 76000000 lines in 1152.9142227172852 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 77000000 lines in 1166.74977684021 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 78000000 lines in 1180.6564447879791 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 79000000 lines in 1194.3947467803955 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 80000000 lines in 1208.9594690799713 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 81000000 lines in 1222.9998469352722 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 82000000 lines in 1236.7506430149078 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 83000000 lines in 1250.6549189090729 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 84000000 lines in 1264.5151867866516 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 85000000 lines in 1278.4629907608032 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 86000000 lines in 1292.3115179538727 seconds\n",
      "barcodes/Larry_system_S6_L001_{}_001.fastq.gz: Processed 87000000 lines in 1305.6937410831451 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 1000000 lines in 14.18682599067688 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 2000000 lines in 28.37282395362854 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 3000000 lines in 41.837624073028564 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 4000000 lines in 55.03623604774475 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 5000000 lines in 69.46619009971619 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 6000000 lines in 84.50206804275513 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 7000000 lines in 99.09303498268127 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 8000000 lines in 113.48776912689209 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 9000000 lines in 127.80174708366394 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 10000000 lines in 142.6965479850769 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 11000000 lines in 157.57042717933655 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 12000000 lines in 172.85713815689087 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 13000000 lines in 189.38178205490112 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 14000000 lines in 205.6254801750183 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 15000000 lines in 221.13742208480835 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 16000000 lines in 235.65522003173828 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 17000000 lines in 249.83126497268677 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 18000000 lines in 264.209037065506 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 19000000 lines in 278.1631369590759 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 20000000 lines in 292.032821893692 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 21000000 lines in 306.4205811023712 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 22000000 lines in 321.13605403900146 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 23000000 lines in 336.88095211982727 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 24000000 lines in 351.4104640483856 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 25000000 lines in 366.34513092041016 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 26000000 lines in 381.13604402542114 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 27000000 lines in 395.3523061275482 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 28000000 lines in 409.89400005340576 seconds\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 29000000 lines in 424.8110420703888 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 30000000 lines in 439.10650420188904 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 31000000 lines in 454.5688261985779 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 32000000 lines in 468.9592549800873 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 33000000 lines in 483.204017162323 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 34000000 lines in 497.52818608283997 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 35000000 lines in 511.8300552368164 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 36000000 lines in 526.5843069553375 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 37000000 lines in 540.3441829681396 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 38000000 lines in 554.0650930404663 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 39000000 lines in 567.8844411373138 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 40000000 lines in 581.8078720569611 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 41000000 lines in 595.5759360790253 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 42000000 lines in 609.5067930221558 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 43000000 lines in 623.676787853241 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 44000000 lines in 637.7840728759766 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 45000000 lines in 651.8482689857483 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 46000000 lines in 665.9295611381531 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 47000000 lines in 680.6526830196381 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 48000000 lines in 695.1109099388123 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 49000000 lines in 709.7385041713715 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 50000000 lines in 725.3905930519104 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 51000000 lines in 740.390248298645 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 52000000 lines in 755.3926050662994 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 53000000 lines in 770.5125832557678 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 54000000 lines in 785.6627900600433 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 55000000 lines in 800.5316441059113 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 56000000 lines in 815.614470243454 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 57000000 lines in 830.5427691936493 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 58000000 lines in 846.8937010765076 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 59000000 lines in 862.0512800216675 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 60000000 lines in 876.7979929447174 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 61000000 lines in 891.6465630531311 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 62000000 lines in 906.6070692539215 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 63000000 lines in 921.5935161113739 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 64000000 lines in 937.0257730484009 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 65000000 lines in 952.2529389858246 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 66000000 lines in 967.3366770744324 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 67000000 lines in 982.8540992736816 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 68000000 lines in 998.4399371147156 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 69000000 lines in 1014.5440740585327 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 70000000 lines in 1031.3018951416016 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 71000000 lines in 1047.05939412117 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 72000000 lines in 1061.982393026352 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 73000000 lines in 1077.019031047821 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 74000000 lines in 1091.9734032154083 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 75000000 lines in 1107.2295632362366 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 76000000 lines in 1122.1756830215454 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 77000000 lines in 1137.2514879703522 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 78000000 lines in 1151.9921391010284 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 79000000 lines in 1166.9366550445557 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 80000000 lines in 1181.9059331417084 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 81000000 lines in 1196.9203219413757 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 82000000 lines in 1214.1221771240234 seconds\n",
      "barcodes/Larry_system_S8_L001_{}_001.fastq.gz: Processed 83000000 lines in 1229.7130711078644 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 1000000 lines in 14.508250951766968 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 2000000 lines in 28.992111921310425 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 3000000 lines in 42.038073778152466 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 4000000 lines in 56.871898889541626 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 5000000 lines in 72.43606615066528 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 6000000 lines in 87.37192392349243 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 7000000 lines in 105.40743708610535 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 8000000 lines in 121.54121804237366 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 9000000 lines in 137.66126799583435 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 10000000 lines in 152.01716589927673 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 11000000 lines in 166.56402587890625 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 12000000 lines in 182.25203585624695 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 13000000 lines in 196.7080729007721 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 14000000 lines in 211.15145182609558 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 15000000 lines in 225.7959599494934 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 16000000 lines in 240.37956500053406 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 17000000 lines in 256.0026068687439 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 18000000 lines in 271.62835001945496 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 19000000 lines in 287.02884793281555 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 20000000 lines in 303.3858890533447 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 21000000 lines in 319.74830389022827 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 22000000 lines in 336.0916759967804 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 23000000 lines in 353.0307321548462 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 24000000 lines in 370.38219594955444 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 25000000 lines in 388.8426549434662 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 26000000 lines in 408.07146096229553 seconds\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 27000000 lines in 424.79282307624817 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 28000000 lines in 441.04706478118896 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 29000000 lines in 456.2985291481018 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 30000000 lines in 471.19148898124695 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 31000000 lines in 486.21552205085754 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 32000000 lines in 501.31154894828796 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 33000000 lines in 516.3139588832855 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 34000000 lines in 532.7527208328247 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 35000000 lines in 549.0961720943451 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 36000000 lines in 565.3879539966583 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 37000000 lines in 581.2682580947876 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 38000000 lines in 598.1516900062561 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 39000000 lines in 614.1828348636627 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 40000000 lines in 630.4159729480743 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 41000000 lines in 646.7817170619965 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 42000000 lines in 663.0252590179443 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 43000000 lines in 679.4278240203857 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 44000000 lines in 696.3770658969879 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 45000000 lines in 713.3322508335114 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 46000000 lines in 730.2148427963257 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 47000000 lines in 746.543224811554 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 48000000 lines in 762.5293970108032 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 49000000 lines in 778.6613440513611 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 50000000 lines in 793.6079409122467 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 51000000 lines in 809.2759220600128 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 52000000 lines in 825.5032279491425 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 53000000 lines in 842.4204549789429 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 54000000 lines in 858.8511719703674 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 55000000 lines in 875.1376819610596 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 56000000 lines in 889.0220649242401 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 57000000 lines in 903.8127629756927 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 58000000 lines in 917.7523877620697 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 59000000 lines in 931.5333399772644 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 60000000 lines in 945.6778991222382 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 61000000 lines in 959.902575969696 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 62000000 lines in 973.9063320159912 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 63000000 lines in 988.1499300003052 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 64000000 lines in 1002.3028290271759 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 65000000 lines in 1016.7270500659943 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 66000000 lines in 1031.8210020065308 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 67000000 lines in 1047.8906118869781 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 68000000 lines in 1062.0173258781433 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 69000000 lines in 1076.0531599521637 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 70000000 lines in 1089.8325350284576 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 71000000 lines in 1103.489984035492 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 72000000 lines in 1117.128420829773 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 73000000 lines in 1131.5499949455261 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 74000000 lines in 1146.2258341312408 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 75000000 lines in 1160.1735899448395 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 76000000 lines in 1174.033961057663 seconds\n",
      "barcodes/Larry_system_S5_L001_{}_001.fastq.gz: Processed 77000000 lines in 1187.6943490505219 seconds\n"
     ]
    }
   ],
   "source": [
    "counts = {}\n",
    "\n",
    "for sample,paths in sample_paths.items():\n",
    "    filtered_cell_barcodes = gzip.open(paths['filtered_cell_barcodes_path']).read().decode('utf-8').split('\\n')\n",
    "    filtered_cell_barcodes_set = set(filtered_cell_barcodes)\n",
    "\n",
    "    for fastq_path in paths['larry_barcode_fastq_paths']:\n",
    "        R1 = gzip.open(fastq_path.format('R1'))\n",
    "        R2 = gzip.open(fastq_path.format('R2'))\n",
    "        counter = 0\n",
    "        start_time = time.time()\n",
    "        while True:\n",
    "            counter += 1\n",
    "            if counter % 1000000 == 0: print(fastq_path+ ': Processed {} lines in {} seconds'.format(counter, time.time()-start_time))\n",
    "            try:\n",
    "                r1_line = R1.readline().decode('utf-8')\n",
    "                r2_line = R2.readline().decode('utf-8')\n",
    "            except:\n",
    "                print('ERROR extracting {}'.format(fastq_path))\n",
    "                break\n",
    "            if r2_line == '': break\n",
    "            if r2_line[0] in '@+': continue\n",
    "            if larry_prefix in r2_line:\n",
    "                larry_bc = r2_line.split(larry_prefix)[1][:40]\n",
    "                cell_bc = r1_line[:16]+'-1'\n",
    "                umi = r1_line[16:24]\n",
    "                if is_valid(larry_bc) and in_filtered_list(cell_bc, filtered_cell_barcodes):\n",
    "                    combo = (sample, cell_bc, umi, larry_bc)\n",
    "                    if combo in counts:\n",
    "                        counts[combo] += 1\n",
    "                    else:\n",
    "                        counts[combo] = 1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "266fba58",
   "metadata": {},
   "source": [
    "## Set parameters for clonal analsysis\n",
    "\n",
    "N_READS: Cutoff for the min number of reads supporting a given UMI/GFP-BC pair in a particular cell\n",
    "\n",
    "N_UMIS: Cutoff for the min number of UMIs supporting a given cell/GFP-BC pair\n",
    "\n",
    "N_HAMMING: Cutoff for the max hamming distance to collapse two GFP-BCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "30762cef",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Retaining 1506741 out of 5280870 (Sample,Cell-BC,UMI,GFP-BC) combinations\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAD4CAYAAAAdIcpQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWY0lEQVR4nO3de5BV1b3g8e9PYpsKreHli5cYmxCJJhnpoJaJReXhhbn2xcq9MbQxxhsDZbQnRZlKZCoEIyEV54aqlKhjLhAKSKGiqcTQc1EzkwQ1jomihSOMITYYQoPKS7TaJKi45g+OPU3bh727z+k+fbq/n6pT9F77tc5i0z/WY68VKSUkSTqW4yqdAUlS/2ewkCRlMlhIkjIZLCRJmQwWkqRM76l0BgBGjRqVJkyYUOlsDApbt24FYNKkSRXOiaRSPfXUU/tSSif3xb36RbCYMGECGzdurHQ2BoVp06YBsGHDhormQ1LpImJHX92ros1QEdEQEUtfffXVSmZDkpShosEipdScUprz/ve/v5LZkCRlsGYhScpkzUKSlMmhs5KkTDZDSZIy2QwlScpkM5QkKVO/eCmvOybM+4+jtv98yz9WKCeSNHjYZyFJymSfhSQpk30WkqRMBgtJUiaDhSQpkx3ckqRMdnBLfejOO+/k1FNPpba2lv379/PYY48xceJEamtruf/++yudPakom6HUpQkTJnDKKafw+uuvt6ctX768ffGkY4kIhg4dSm1tLWPGjOGGG27g8OHD7funTZvGe9/7Xmpra9s/DQ0NR13jhRde4LjjjuNrX/vaMa8/cuRIPv3pT7N27dqjjtmyZQuXXHIJI0aMYNiwYUyZMoX169d3sxSObcOGDYwdOzb38W+++SY33HADv/rVr2hra2PkyJEsWLCApqYm2trauOyyy8qaP6mcDBYq6vDhw9x66609OveZZ56hra2Nhx9+mLVr17JixYqj9t9+++20tbW1f5qbm4/av3r1aoYPH87atWs5dOhQ0etv3bqVq6++mqamJm6++eb2/Q0NDXz2s5/lpZdeYs+ePSxZsoSTTjqpR9+lXF5++WX+/ve/8+EPf7g9bceOHUdtS/2VwUJFffOb32Tx4sUcPHiwx9eoq6vjoosuYtOmTbnPSSmxevVqFi1axPHHH/+uQNLRqFGj+NKXvsSdd97JD37wA/bv38++fft44YUXmD17NjU1NdTU1HDRRRfxiU98ouh1li1bxtlnn82JJ57I5MmTefrpp4EjtZiWlpb2466++mrmz5/P66+/zowZM9i9e3d77Wj37t0cOnSIuXPnMnr0aEaPHs3cuXM5dOgQf/rTn9rXPR82bBif+tSnOOuss9i+fTsNDQ3U1tZ2GRSl/sJgoaLq6+uZNm0aixcv7vE1/vjHP/Loo49SV1eX+5zf/e53tLa2MmvWLC6//HJWrVqVec7MmTN56623eOKJJxg5ciR1dXVceeWV3H///bz88svHPPe+++7ju9/9LqtXr+a1115j3bp1jBw58pjnDB06lAceeIDRo0e3145Gjx7N97//fX7/+9+zadMmnnnmGZ544gkWLVrEBz/4QbZs2QLAwYMH+c1vfsO2bdsYP348zc3NtLW1ccIJJ+QuI6mvlT1YRMRxEfH9iLgtIr5c7uurby1cuJDbbruNvXv3duu88847j6FDh3L22Wczbdo0rrvuuqP2f/3rX2fYsGHtn+985zvt+1atWsWMGTMYPnw4V1xxBQ8++CB79uw55v2OP/54Ro0axYEDB4gIfvvb3zJhwgS+8Y1vcPrpp3PxxRfz/PPPd3nu8uXL+da3vsXHP/5xIoK6ujrOOOOMbn3fd6xZs4YFCxZwyimncPLJJ3PTTTfx05/+tEfXkvqTXMEiIlZExJ6I2NwpfXpEbI2IloiYV0ieCYwF3gRay5td9bVzzjmHSy+9lFtuuaVb5z399NO0tbWxdu1a/vCHPxzVUQ6wZMkSDh482P753ve+B8Df/vY37rvvPr74xS8CcOGFFzJ+/HjuuuuuY97vzTffZO/evYwYMQKAsWPHcvvtt7Nt2zZ27NjB0KFDueqqq7o8d+fOnZx11lnd+n7F7N69+6hAc8YZZ7B79+6yXFuqpLw1i5XA9I4JETEEuAOYAUwGGiNiMjAJ+N8ppRuAdw9lUdW5+eabWbZsGbt27erWeRHB5ZdfzoUXXsjChQtznfOLX/yC1157jeuuu47TTjuN0047jV27dmU2Rf3yl7/kPe95D1OnTn3XvnHjxnH99dezefPmLs48sn/btm1d7nvf+97HX//61/btl156qf3niHjX8aNHj2bHjh3t23/5y18YPXr0MfMuVYNcwSKl9AhwoFPyVKAlpbQ9pfQGcA9HahWtwCuFYw5TRETMiYiNEbGxu00c6lt1dXV84QtfYMmSJT06f968eSxbtuyoX7TFrFq1iq985Ss8++yzbNq0iU2bNvHYY4/xzDPP8Oyzz77r+AMHDrBmzRquv/56brzxRkaOHMkrr7zCTTfdREtLC2+//Tb79u1jxYoVXHDBBV3e86tf/SqLFy/mqaeeIqVES0tL+y/8j33sY9x1110cPnyYBx98kIcffrj9vFNPPZX9+/fT8aXSxsZGFi1axN69e9m3bx8LFy7kyiuv7G6RSf1OKX0WY4CdHbZbC2k/B/4hIm4DHil2ckppaUqpPqVUf/LJJ5eQDfWFBQsWvKspKa9zzz2Xiy++mB/+8IftaU1NTUe9ZzFlyhR27drFr3/9a+bOndteqzjttNOYMmUK06dPP6p28dGPfpTa2lrq6upYvnw5P/rRj9prLzU1Nfz5z3/mM5/5DCeddBLnnHMOJ5xwAitXruwyf5///Of59re/zRVXXMGJJ57IZZddxoEDR/5vdOutt9Lc3MywYcNYs2bNUe9CfOhDH6KxsZEPfOADDBs2jN27dzN//nzq6+v5yEc+wrnnnst5553H/Pnze1RuUn8SKaV8B0ZMAP5HSumcwva/ANNTSl8tbH8JOD+l1JT75hENQENdXd3sYp2Pnbn4UWneealuw4YNFc2HpNJFxFMppfq+uFcpNYtdwLgO22MLaZKkAaaUYPEkMDEizoyIGmAWsK47F3BuqOr06KOPHtWE1PEjaWDKtQZ3RNwNTANGRUQrcFNK6ScR0QQ8BAwBVqSUtnTn5h2aobqXa1XUJz/5Sdra2iqdDUl9KFewSCk1FklfD/R4draUUjPQXF9fP7un15Ak9T7Xs5AkZXI9C0lSJicSlCRlshlKkpTJZihJUiaboSRJmWyGkiRlshlKkpQp10t5A13nyQnBCQolqSOboSRJmWyGkiRlcjSUJCmTwUKSlMkOblWUgwuk6mDNQpKUqaI1Cxc/UldcZ13qfyoaLFz8aPDpqtlJUv9nM5QkKZMd3OpV1iSkgcFgoX7PPgyp8myGkiRlKnvNIiKmAd8DtgD3pJQ2lPseGtx8N0Pqe7mCRUSsAC4F9qSUzumQPh24FRgCLE8p3QIkoA14L9Ba9hyr3/CXtjR45G2GWglM75gQEUOAO4AZwGSgMSImA4+mlGYANwI3ly+rkqRKyRUsUkqPAAc6JU8FWlJK21NKbwD3ADNTSm8X9r8CnFDsmhExJyI2RsTGvXv39iDrkqS+UkqfxRhgZ4ftVuD8iPgc8A/AMOD2YienlJZGxItAQ01NzZQS8qF+xKGy0sBU9g7ulNLPgZ/nPLZq3+C2vV7SYFLK0NldwLgO22MLabm5Up4kVYdSgsWTwMSIODMiaoBZwLruXMCV8iSpOuQKFhFxN/A4MCkiWiPimpTSW0AT8BDwHHBvSmlLd25uzUKSqkOuPouUUmOR9PXA+p7evJr7LAYjO6+lwaui031Ys5Ck6uB6FgIc3SXp2KxZSJIyWbPoQ/7vXVK1cj0LDQiueSH1LpuhJEmZKhosfClPkqqDzVD9jM0pkvojl1WVJGWqaM0iIhqAhrq6ukpmo1/L89Z0b9U+fGNb0jscOjsAVDKgSBoc7LMYJN4JKC9t33/UtiTlYZ+FJCmTNQsNSI4qk8rLl/IkSZl8KU+SlMk+C0lSJoOFJCmTwUKSlMlgIUnK1CvBIiKGRsTGiLi0N64vSepbuYJFRKyIiD0RsblT+vSI2BoRLRExr8OuG4F7y5lRSVLl5K1ZrASmd0yIiCHAHcAMYDLQGBGTI+KzwP8F9pQxn5KkCsr1BndK6ZGImNApeSrQklLaDhAR9wAzgVpgKEcCyN8iYn1K6e3O14yIOcAcgPHjx/f4C0iSel8p032MAXZ22G4Fzk8pNQFExNXAvq4CBUBKaSmwFKC+vj6VkA9JUi/rtbmhUkors45xPQtJqg6ljIbaBYzrsD22kCZJGmBKCRZPAhMj4syIqAFmAeu6cwHnhpKk6pB36OzdwOPApIhojYhrUkpvAU3AQ8BzwL0ppS3dubmzzkpSdcg7GqqxSPp6YH1Pb+6yquorXa0M6BoXUn6uZyFJyuR6FpKkTNYsJEmZrFlIkjI5RbkkKZPNUJKkTDZDSZIy2QwlScpksJAkZbLPQpKUyT4LSVImm6EkSZkMFpKkTAYLSVKmXltWVervOk9b7pTlUnGOhpIkZXI0lCQpk30WkqRMBgtJUiaDhSQpU9mDRUScHRE/joifRcTXyn19SVLfyxUsImJFROyJiM2d0qdHxNaIaImIeQAppedSStcClwMXlT/LkqS+lrdmsRKY3jEhIoYAdwAzgMlAY0RMLuz7J+A/gPVly6kkqWJyBYuU0iPAgU7JU4GWlNL2lNIbwD3AzMLx61JKM4AvljOzkqTKKOUN7jHAzg7brcD5ETEN+BxwAseoWUTEHGAOwPjx40vIhlQevtEtFVf26T5SShuADTmOWxoRLwINNTU1U8qdD0lS+ZQyGmoXMK7D9thCWm6+wS1J1aGUYPEkMDEizoyIGmAWsK47F3BuKEmqDnmHzt4NPA5MiojWiLgmpfQW0AQ8BDwH3JtS2tJ7WZUkVUquPouUUmOR9PWUMDw2pdQMNNfX18/u6TUkSb3PKcolSZmcolySlMmJBCVJmWyGkiRlshlKkpTJZihJUiaboSRJmWyGkiRlKvtEgtJA0XkWWnAmWg1e9llIkjLZZyFJymSfhSQpk81QkqRMBgtJUiaDhSQpkx3ckqRMdnBLkjLZDCVJyuQb3FI3dH6r2ze6NVhYs5AkZTJYSJIy9UozVERcBvwjcBLwk5TSr3rjPpKkvpG7ZhERKyJiT0Rs7pQ+PSK2RkRLRMwDSCndn1KaDVwLfKG8WZYk9bXuNEOtBKZ3TIiIIcAdwAxgMtAYEZM7HDK/sF+SVMVyN0OllB6JiAmdkqcCLSml7QARcQ8wMyKeA24BHkgpPd3V9SJiDjAHYPz48T3IulR5jo7SYFFqn8UYYGeH7VbgfOC/AJ8B3h8RdSmlH3c+MaW0NCJeBBpqamqmlJgPSVIv6pXRUCmlJSmlKSmla7sKFB2O8w1uSaoCpQaLXcC4DttjC2m5ODeUJFWHUpuhngQmRsSZHAkSs4Ar8p6cUmoGmuvr62eXmA+p37JfQwNBd4bO3g08DkyKiNaIuCal9BbQBDwEPAfcm1La0o1rWrOQpCrQndFQjUXS1wPre3JzaxaSVB0qOpFgRDQADXV1dZXMhlQ2nZucpIHC9SwkSZlcKU+SlMmahSQpk1OUS5Iy2cEt9THfu1A1shlKkpTJZihJUiaDhSQpk0NnJUmZ7LOQJGWyGUqSlMlgIUnKVNH3LCR1Pfmg716ov7GDW5KUyQ5uSVIm+ywkSZkMFpKkTAYLSVImg4UkKVPZg0VEfCAifhIRPyv3tSVJlZHrPYuIWAFcCuxJKZ3TIX06cCswBFieUrolpbQduMZgIfWca16ov8lbs1gJTO+YEBFDgDuAGcBkoDEiJpc1d5KkfiFXsEgpPQIc6JQ8FWhJKW1PKb0B3APMLHP+JEn9QCl9FmOAnR22W4ExETEyIn4M/KeI+K/FTo6IORGxMSI27t27t4RsSJJ6W9nnhkop7QeuzXHc0oh4EWioqamZUu58SJLKp5SaxS5gXIftsYU0SdIAU0qweBKYGBFnRkQNMAtY150LODeUJFWHXMEiIu4GHgcmRURrRFyTUnoLaAIeAp4D7k0pbenOzZ11VpKqQ64+i5RSY5H09cD6nt48pdQMNNfX18/u6TUkSb3P9SwkSZlcz0KSlMmJBCVJmSq6BndENAANdXV1lcyG1O85V5QqzWYoSVImm6EkSZkcDSVJymQzlCQpk81QkqRMjoaSBghHTKk32QwlScpkM5QkKZPBQpKUyWAhScpksJAkZXI0lDSIOYJKeTkaSpKUyWYoSVImg4UkKZPBQpKUyWAhScpU9tFQETEU+O/AG8CGlNKact9DktS3ctUsImJFROyJiM2d0qdHxNaIaImIeYXkzwE/SynNBv6pzPmVJFVA3maolcD0jgkRMQS4A5gBTAYaI2IyMBbYWTjscHmyKUmqpFzNUCmlRyJiQqfkqUBLSmk7QETcA8wEWjkSMDZxjGAUEXOAOQDjx4/vbr6lQa3zy3RSbyulg3sM/78GAUeCxBjg58A/R8SdQHOxk1NKS4GbgadrampKyIYkqbeVvYM7pfQ68K85j20Gmuvr62eXOx+SpPIppWaxCxjXYXtsIS23iGiIiKWvvvpqCdmQJPW2UoLFk8DEiDgzImqAWcC67lzAuaEkqTrkHTp7N/A4MCkiWiPimpTSW0AT8BDwHHBvSmlLd25uzUKSqkPe0VCNRdLXA+t7enP7LCSpOlR0ug9rFpJUHVzPQpKUyZqFJClTpJQqnQciYi+wI+fho4B9vZidwcAyLJ1lWDrLsHSTUkon9sWNKroG9ztSSifnPTYiNqaU6nszPwOdZVg6y7B0lmHpImJjX93L9SwkSZkMFpKkTNUYLJZWOgMDgGVYOsuwdJZh6fqsDPtFB7ckqX+rxpqFJKmPGSwkSZn6ZbAotuZ3h/0REUsKa3//n4g4r6/z2B91VW4RMSIi/mdEPF/4c3iRc79cOOb5iPhy3+W68rpTbnmfvYiYEhHPFo5bEhHRV9+nr5Sr3PI8e3mf42rQ2+WW59nr0e/QlFK/+wAXA+cBm4vs/8/AA0AAFwB/qHSe+8Onq3ID/g2YV/h5HvDfujhvBLC98Ofwws/DK/19+mO55X32gCcK+6Nw/IxKf8/+WG55n708z3G1fHq73PI8ez35HVrxgjtGgU44RrD4d6Cxw/ZW4PRK57k/fDqXW8eyAU4HtnZxTiPw78XKdzB88pZbnmevcPwfi5XvQPqUWm55n708z3E1fXqr3PI+ez35Hdovm6FyKLb+t97t1JTSi4WfXwJO7eIYy/PdipVbnrIaU0g/1jEDVXfLLe+zl+c5rmblKre8z163/81Xa7BQD6Qj/4VwrHQ3WW4901vlNtD/Pvrr96vWYFHy+t+DyMsRcTpA4c89XRxjeb5bsXLLU1a7CunHOmag6m655X328jzH1axc5Zb32ev2v/lqDRbrgKsKPfoXAK92qMLpaOuAd0ZKfBn4ZRfHPARcEhHDC6MwLimkDWbFyi3z2StsvxYRFxRGolxF1+U+EHW33PI+e3me42pWlnLrxrPX/d+hle7oKdL5czfwIvAmR9rSrgGuBa4t7A/gDmAb8CxQX+k894dPkXIbCfwaeB74X8CIwrH1wPIO534FaCl8/rXS36Ufl1vRZw/Y1OHnemBz4bjbKcyWMJA+ZSy3Lp89YPk7xxW7bjV++qDcunz2Sv0d6nQfkqRM1doMJUnqQwYLSVImg4UkKZPBQpKUyWAhScpksJAkZTJYSJIy/T+4zDufaEy2YAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "num_reads = [v for k,v in counts.items()]\n",
    "plt.hist(np.log(num_reads)/np.log(10), bins=50)\n",
    "plt.axvline(np.log(N_READS)/np.log(10),c='k')\n",
    "plt.xticks(range(5),np.logspace(0,4,5))\n",
    "plt.text(np.log(N_READS)/np.log(10)*1.1,10**6,'N_READS cutoff', fontsize=12)\n",
    "plt.yscale('log')\n",
    "\n",
    "counts_filtered = {k:v for k,v in counts.items() if v >= N_READS}\n",
    "print('Retaining '+repr(len(counts_filtered))+ ' out of '+repr(len(counts))+' (Sample,Cell-BC,UMI,GFP-BC) combinations')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74776937",
   "metadata": {},
   "source": [
    "## Collapse GFP-BCs by hamming distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "12108f1c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mapped 500 out of 6910 barcodes\n",
      "Mapped 1000 out of 6910 barcodes\n",
      "Mapped 1500 out of 6910 barcodes\n",
      "Mapped 2000 out of 6910 barcodes\n",
      "Mapped 2500 out of 6910 barcodes\n",
      "Mapped 3000 out of 6910 barcodes\n",
      "Mapped 3500 out of 6910 barcodes\n",
      "Mapped 4000 out of 6910 barcodes\n",
      "Mapped 4500 out of 6910 barcodes\n",
      "Mapped 5000 out of 6910 barcodes\n",
      "Mapped 5500 out of 6910 barcodes\n",
      "Mapped 6000 out of 6910 barcodes\n",
      "Mapped 6500 out of 6910 barcodes\n",
      "\n",
      "Collapsed 2404 barcodes\n"
     ]
    }
   ],
   "source": [
    "def hamming(bc1,bc2): return np.sum([x1 != x2 for x1,x2 in zip(bc1,bc2)])\n",
    "\n",
    "all_gfp_bcs = sorted(set([k[3] for k in counts_filtered]))\n",
    "good_gfp_bcs = []\n",
    "bc_map = {}\n",
    "for i,bc1 in enumerate(all_gfp_bcs):\n",
    "    if i > 0 and i % 500 == 0: print('Mapped '+repr(i)+' out of '+repr(len(all_gfp_bcs))+' barcodes')\n",
    "    mapped = False\n",
    "    for bc2 in good_gfp_bcs:\n",
    "        if hamming(bc1,bc2) <= N_HAMMING:\n",
    "            mapped = True\n",
    "            bc_map[bc1] = bc2\n",
    "            break\n",
    "    if not mapped:\n",
    "        good_gfp_bcs.append(bc1)\n",
    "\n",
    "print('\\nCollapsed '+repr(len(bc_map))+' barcodes')\n",
    "for bc in good_gfp_bcs: bc_map[bc] = bc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "405f016e",
   "metadata": {},
   "source": [
    "## Filter GFP-barcodes by UMI\n",
    "Index (UMI, GFP-BC) combinations to cells, then filter out BCs supported by fewer than N_UMIS UMIs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "a2e9c2b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Final annotation has 3836 clones in 5675 cells\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAskElEQVR4nO3deXxU9dn38c+VnewJCUlIAoQtYQdFXEBcsBapFdtahGrFpXqr1do+ba297a3cPnbxXmq1j8VaRdS24kKrtqLWtYpiJbJD2BICSQghbIEkQLbr+eOcxCFMYCCTTCZzvV+veeXM7yxzzWT55pzfOb8jqooxxpjQExboAowxxgSGBYAxxoQoCwBjjAlRFgDGGBOiLACMMSZERQS6gFORlpamgwYNCnQZxhgTVD7//PM9qprevj2oAmDQoEEUFhYGugxjjAkqIrLdW7sdAjLGmBBlAWCMMSHKAsAYY0KUBYAxxoQoCwBjjAlRPgWAiEwXkU0islVE7vEyf6CIvCsia0TkAxHJ8Zg3V0S2uI+5Hu1nishad5uPioj45y0ZY4zxxUkDQETCgceAy4CRwBwRGdlusf8BnlXVscADwC/ddVOB+4GzgUnA/SKS4q4zH7gZGOY+pnf63RhjjPGZL3sAk4Ctqlqiqg3AImBmu2VGAu+50+97zP8y8Laq7lPV/cDbwHQRyQISVfVTdcajfha4snNvpWOvrqrgj596PQ02KIgItoNkjPE3XwIgGyjzeF7utnlaDXzdnf4akCAifU+wbrY7faJtAiAit4hIoYgUVldX+1Du8ZasreSppdtOa11jjOmt/NUJ/CPgAhFZCVwAVADN/tiwqj6hqhNVdWJ6+nFXMvukIDOR0r11HG7wS0nGGNMr+BIAFUCux/Mct62Nqu5U1a+r6gTgXrftwAnWrXCnO9ymPxVkJqAKm6sOddVLGGNM0PElAJYDw0QkT0SigNnAa54LiEiaiLRu66fAAnf6LeBSEUlxO38vBd5S1UrgoIic4579cx3wqh/ej1cFWYkAbNplAWCMMa1OGgCq2gTcgfPHvAh4UVXXi8gDInKFu9iFwCYR2QxkAD93190H/F+cEFkOPOC2AdwOPAlsBYqBN/z1ptobkBpLTGQYGy0AjDGmjU+jgarqEmBJu7b7PKZfBl7uYN0FfLFH4NleCIw+lWJPV3iYkJ+RwMZdB7vj5YwxJiiEzJXA+ZkJbNx1COesU2OMMSEUAInsq2uguvZooEsxxpgeIWQCYERmAmAdwcYY0ypkAiDfDYCNlRYAxhgDIRQAfeOjSU+ItjOBjDHGFTIBAM4FYZuq7EwgY4yBEAuA/IwENlfV0tTcEuhSjDEm4EIqAAqyEmloaqF0b32gSzHGmIALrQBo7Qi2C8KMMSa0AmBov3jCxE4FNcYYCLEAiIkMJy8tzs4EMsYYQiwAwOkHsENAxhgTigGQkUDZvsPUHm0KdCnGGBNQIRcA+TYkhDHGACEYACPs5jDGGAOEYABkJ/chLiqcTdYPYIwJcSEXAGFhQn5mAkW2B2CMCXE+BYCITBeRTSKyVUTu8TJ/gIi8LyIrRWSNiMxw268RkVUejxYRGe/O+8DdZuu8fn59ZyeQn5nIJrs5jDEmxJ00AEQkHHgMuAwYCcwRkZHtFvsZzr2CJ+DcNP53AKr6J1Udr6rjgW8D21R1lcd617TOV9XdnX43PirITKDmcCO7Dh7prpc0xpgex5c9gEnAVlUtUdUGYBEws90yCiS600nATi/bmeOuG3BfDAlhh4GMMaHLlwDIBso8npe7bZ7mAdeKSDnOzePv9LKdq4Hn27U97R7++Q8REW8vLiK3iEihiBRWV1f7UO7JFWTamUDGGOOvTuA5wEJVzQFmAM+JSNu2ReRsoF5V13msc42qjgHOdx/f9rZhVX1CVSeq6sT09HS/FJsUG0lWUgwbK+1MIGNM6PIlACqAXI/nOW6bp5uAFwFUdRkQA6R5zJ9Nu//+VbXC/XoI+DPOoaZuk5+ZYIeAjDEhzZcAWA4ME5E8EYnC+WP+WrtldgDTAERkBE4AVLvPw4BZeBz/F5EIEUlzpyOBy4F1dKOCzESKq2tptJvDGGNC1EkDQFWbgDuAt4AinLN91ovIAyJyhbvYD4GbRWQ1zn/61+sX51hOBcpUtcRjs9HAWyKyBliFs0fxB3+8IV8VZCbQ2KyUVNd158saY0yPEeHLQqq6BKdz17PtPo/pDcDkDtb9ADinXVsdcOYp1upX+R43h2mdNsaYUBJyVwK3GpIeT0SYWD+AMSZkhWwAREWEMSQ93k4FNcaErJANAICCrAQ7FdQYE7JCOgDyMxPYWXOEmsONgS7FGGO6XUgHQOuQEJur7DCQMSb0hHgAOENC2GGg4PPxxx8zduxYoqKiuPDCCztsM8Z0LKQDICsphoSYiKA+E+jCCy/kjjvuOK594cKFxMfHtz2fN28eIsK0adOOW3b+/PmICKNHj+5w/a7wwQcfICLs2bPnlNe96667GDduHMXFxfzlL3/psM0Y07GQDgARYURmYlAHwKnIzMxk6dKllJaWHtP+1FNPMWDAgMAUdZq2bt3KxRdfTG5uLqmpqR22GWM6FtIBAE5H8OYQuTlM3759+cpXvsLTTz/d1rZmzRo2btzIVVdddcrba2ho4N///d8ZOHAg0dHRDB48mEcffRTw/t99aWkpIkJhYSGlpaVcdNFFAKSnpyMiXH/99QAcPXqU73//+2RkZBATE8M555zD0qVLj9lGTU0NN954IyLCwoULvbYZY04s5AOgICuBQ0ebqDhwONCldIubbrqJZ555hpYWZwykp556ilmzZpGQcOpXQ8+dO5dnn32WX//61xQVFfHUU0+RnJzs07q5ubksXrwYgPXr11NZWckjjzwCwN13380LL7zAggULWLlyJWPGjGH69OlUVlaSm5tLZWUlsbGx/OY3v6GyspJvfvObx7VdffXVp/x+jAk1Pg0F0Zu13Rym8hA5KbEBrqbrTZ8+ncbGRt59912mTp3KH//4R1599VXeeeedU9rOli1bWLRoEW+88QbTp08HYPDgwT6vHx4e3naYpl+/fqSlOYPH1tXVMX/+fJ588km+8pWvAPD444/z3nvv8dhjj/Hggw+SmZmJiJCUlERmZiYAcXFxx7UZY04s5PcAhmc4AbApRE4FDQ8PZ+7cuSxYsIBXXnmF9PR0pkyZcsrbWblyJWFhYW2HcfyluLiYxsZGJk/+Ymip8PBwzj33XDZs2ODX1zIm1IX8HkBCTCQ5KX0oCtJTQRMTE6mpqTmu/cCBAyQlJXld54YbbmDs2LGUlpZyww03dEldYWHO/xaefSuNjZ274K6Dm8YZY05TyO8BgHMYKFjHBMrPz2fFihXHdWKvWLGC/Px8r+sMGzaMSZMmUVhYyNy5c0/rdcePH09LSwvvv/++1/mtd2+rrKxsa1u1atUxy0RFRQHQ3Nzc1jZkyBCioqL4+OOP29qam5tZtmwZI0eOPK1ajTHeWQDgXBBWsqeOo03NJ1+4h7ntttsoKSnhzjvvZPXq1WzatImHH36Y559/nh//+McdrvfGG2+we/fu0z5ePnz4cGbNmsV3vvMdFi9ezLZt2/joo4947rnnABg6dCi5ubnMmzePzZs3849//IMHH3zwmG0MHDgQEeH111+nurqa2tpa4uLiuO222/jJT37CkiVLKCoq4rbbbqOqqorbb7/9tGo1xnhnAYBzKmhzi7J1d22gSzllgwcP5sMPP2TLli1ceumlTJo0iUWLFvHSSy9x2WWXdbhebGwsKSkpnXrtZ599lm9961t873vfo6CggOuvv77tcFRkZCSLFi2ipKSEcePGcf/99/OLX/zimPWzs7P5z//8T+69914yMjLaLmh76KGHuPrqq7nhhhsYP348a9as4c033yQrK6tT9RpjjiXBdP77xIkTtbCw0O/b3VJ1iC89/CG/njWOr5+R4/ftd1brse9g+l4ZY3oOEflcVSe2b/dpD0BEpovIJhHZKiL3eJk/QETeF5GVIrJGRGa47YNE5LCIrHIfj3usc6aIrHW3+agEsIcvLy2OqPCwkLki2BhjwIcAEJFw4DHgMmAkMEdE2vfG/QznXsETcG4a/zuPecWqOt593OrRPh+4GRjmPqaf/tvonIjwMIb2i7cAMMaEFF/2ACYBW1W1RFUbgEXAzHbLKJDoTicBO0+0QRHJAhJV9VP35vHPAleeSuH+ZjeHMcaEGl8CIBso83he7rZ5mgdcKyLlODePv9NjXp57aOifInK+xzbLT7JNAETkFhEpFJHC6upqH8o9PQWZCew+dJT9dQ1d9hrGGNOT+OssoDnAQlXNAWYAz4lIGFAJDHAPDf0f4M8ikniC7RxHVZ9Q1YmqOrH13PKu0HZvADsMZIwJEb4EQAWQ6/E8x23zdBPwIoCqLgNigDRVPaqqe932z4FiYLi7vufpNt622a3axgTaZYeBjDGhwZcAWA4ME5E8EYnC6eR9rd0yO4BpACIyAicAqkUk3e1ERkQG43T2lqhqJXBQRM5xz/65DnjVL+/oNKUnRJMSGxm0VwQbY8ypOulYQKraJCJ3AG8B4cACVV0vIg8Ahar6GvBD4A8i8gOcDuHrVVVFZCrwgIg0Ai3Araq6z9307cBCoA/whvsIGBGhIDORIgsAY0yI8GkwOFVdgtO569l2n8f0BmCyl/UWA4s72GYhMNrbvEDJz0zgxcIyWlqUsDAbeMwY07vZUBAeRmQlUN/QTNn++kCXYowxXc4CwEO+eyZQUaUdBjLG9H4WAB6GZ8QjgnUEG2NCggWAh9ioCAamxtqpoMaYkGAB0E5+EN8cxhhjToUFQDsFmYmU7q3jcEPw3RzGGGNOhQVAOwWZCbQobNltewHGmN7NAqCdgiwbE8gYExosANoZkBpLTGQYG+1UUGNML2cB0E54mDA8I4FNVXYmkDGmd7MA8KLAzgQyxoQACwAv8jMT2VPbQPWho4EuxRhjuowFgBcj3HsD2F6AMaY3swDwIt9uDmOMCQEWAF70jY8mPSHaTgU1xvRqFgAdsI5gY0xvZwHQgfyMBDZXHaK5RQNdijHGdAkLgA4UZCVytKmF0r11gS7FGGO6hE8BICLTRWSTiGwVkXu8zB8gIu+LyEoRWSMiM9z2L4nI5yKy1v16scc6H7jbXOU++vnvbXVeQWtHsF0RbIzppU4aACISDjwGXAaMBOaIyMh2i/0MeFFVJwCzgd+57XuAr6rqGGAu8Fy79a5R1fHuY3cn3offDe0XT5jAJjsTyBjTS/myBzAJ2KqqJaraACwCZrZbRoFEdzoJ2AmgqitVdafbvh7oIyLRnS+768VEhpOXFmdnAhljei1fAiAbKPN4Xu62eZoHXCsi5cAS4E4v2/kGsEJVPS+vfdo9/PMfIiLeXlxEbhGRQhEprK6u9qFc/ynITLQAMMb0Wv7qBJ4DLFTVHGAG8JyItG1bREYBDwH/5rHONe6hofPdx7e9bVhVn1DViao6MT093U/l+qYgM4Ed++qpO9rUra9rjDHdwZcAqAByPZ7nuG2ebgJeBFDVZUAMkAYgIjnAX4HrVLW4dQVVrXC/HgL+jHOoqUdpvSJ4U5XtBRhjeh9fAmA5MExE8kQkCqeT97V2y+wApgGIyAicAKgWkWTgdeAeVf24dWERiRCR1oCIBC4H1nXyvfjdCPfmMHZBmDGmNzppAKhqE3AH8BZQhHO2z3oReUBErnAX+yFws4isBp4HrldVddcbCtzX7nTPaOAtEVkDrMLZo/iDn99bp2Un9yEuKtwCwBjTK0X4spCqLsHp3PVsu89jegMw2ct6DwIPdrDZM30vMzDCwoT8zASKKu1UUGNM72NXAp9EfmYim6oO4ezQGGNM72EBcBIFmQkcqG+k6qDdHMYY07tYAJxEgd0bwBjTS1kAnERBpnMmkF0QZozpbSwATiIpNpKspBg7E8gY0+tYAPhgVP9Elm7dw/66hkCXYowxfmMB4IO7pg3nQH0DP355jZ0NZIzpNSwAfDAmJ4l7LhvBO0VVPPNJaaDLMcYYv7AA8NGNkwcxraAfv1iykXUVNYEuxxhjOs0CwEciwn9/cxypcVHc+fxKam2EUGNMkLMAOAWpcVE8Mns82/fW8R+vrLP+AGNMULMAOEVnD+7LXdOG89eVFSxe0X5UbGOMCR4WAKfhjouHcs7gVP7jlXUUV9cGuhxjjDktFgCnITxMeGT2BPpEhfPdP63gSGNzoEsyxphTZgFwmjISY/ifb45l465D/GJJUaDLMcaYU2YB0AkXF2TwnSl5PLtsO2+uqwx0OcYYc0osADrp7ukFjM1J4u6X11C+vz7Q5RhjjM98CgARmS4im0Rkq4jc42X+ABF5X0RWisgaEZnhMe+n7nqbROTLvm4zWERFhPHbORNoUfje8ytpbG4JdEnGGOOTkwaAiIQDjwGXASOBOSIyst1iP8O5V/AEnJvG/85dd6T7fBQwHfidiIT7uM2gMbBvHL/8+hhW7DjAw29vDnQ5xhjjE1/2ACYBW1W1RFUbgEXAzHbLKJDoTicBO93pmcAiVT2qqtuAre72fNlmUPnquP7MPiuX+f8s5qMt1YEuxxhjTsqXAMgGyjyel7ttnuYB14pIOc7N4+88ybq+bBMAEblFRApFpLC6umf/Yb3/q6MYmh7PD15YTfUhu4WkMaZn81cn8BxgoarmADOA50TEL9tW1SdUdaKqTkxPT/fHJrtMn6hw/t+3zuDQkUb+z4uraGmxoSKMMT2XL3+kK4Bcj+c5bpunm4AXAVR1GRADpJ1gXV+2GZTyMxO4/6uj+GjLHh7/sDjQ5RhjTId8CYDlwDARyRORKJxO3dfaLbMDmAYgIiNwAqDaXW62iESLSB4wDPjMx20GrTmTcvnKmCz+9x+b+Xz7/kCXY4wxXp00AFS1CbgDeAsowjnbZ72IPCAiV7iL/RC4WURWA88D16tjPc6ewQbgTeC7qtrc0Tb9/eYCRUT45TfGkJUUw/eeX0lNfWOgSzLGmONIMA1pPHHiRC0sLAx0GT5bVXaAq+Z/wiUjMph/7RmIyGltp3W9YPpeGWN6DhH5XFUntm+3K4G70PjcZO6ens+b63fZrSSNMT2OBUAX+86UwVyUn868v23gV29spMmuFDbG9BAWAF0sLEyYf+2ZzJk0gMf/Wcx1Cz5jT61dI2CMCTwLgG4QExnOL78+hv+6aiyF2/fz1d8uZeUOOzvIGBNYFgDdaNbEXP5y23mEhwmzfr+MP3663Tp2jTEBYwHQzUZnJ/H3O6dw3pA0fvbKOn740moON9gdxYwx3c8CIACSY6N4+vqzuGvaMP66soKvz/+E7XvrAl2WMSbEWAAESFiY8IMvDWfB3LOo2F/PV3+7lPc2VgW6LGNMCLEACLCLCvrx9zvPJycllhsXFvLrtzfTbIPIGWO6gQVADzCgbyx/uf08vnFGDo++u4UbFi5nf11DoMsyxvRyFgA9RExkOP/zzbH8/Guj+bR4L5f/dilry2sCXZYxphezAOhBRIRrzh7Ii7eei6ryjcc/4YXlOwJdljGml7IA6IHG5ybztzuncNagFH6yeG2gyzHG9FIWAD1U3/honr3xbG6/cEhb27Y9dqqoMcZ/LAB6sPAw4e7pBW3PL3/0I15d1StunGaM6QEsAILIiKxE7lq0ih+/tJr6hqZAl2OMCXIWAEFk0S3ncOfFQ3l5RTmX/3YpG3YeDHRJxpgg5lMAiMh0EdkkIltF5B4v8x8WkVXuY7OIHHDbL/JoXyUiR0TkSnfeQhHZ5jFvvB/fV68UER7GDy/N5083nc2hI01c+buPeW5ZqQ0oZ4w5LSe9JaSIhAObgS8B5Tg3dJ+jqhs6WP5OYIKq3tiuPRXYCuSoar2ILAT+rqov+1pssN0S0l+83RJyT+1RfvTSaj7YVM2XR2XwX98YR1JsZKBKNMb0YJ25JeQkYKuqlqhqA7AImHmC5efg3Bi+vauAN1S13peCzYmlxUezYO5Z3DtjBO8W7WbGox9RWLov0GUZY4KILwGQDZR5PC93244jIgOBPOA9L7Nnc3ww/FxE1riHkKI72OYtIlIoIoXV1dU+lBs6wsKEm6cO5mX3HgNXP/Epj72/1cYSMsb4xN+dwLOBl1X1mAHuRSQLGAO85dH8U6AAOAtIBX7ibYOq+oSqTlTVienp6X4ut3cYn5vM3783hctGZ/Lfb23iugX/YvfBI4EuyxjTw/kSABVArsfzHLfNG2//5QPMAv6qqo2tDapaqY6jwNM4h5rMaUqMieS3cybw0DfG8Pn2/Vz2yEd8sGl3oMsyxvRgvgTAcmCYiOSJSBTOH/nX2i8kIgVACrDMyzaO6xdw9woQp4fzSmDdKVVujiMiXH3WAP52xxTS4qO5/unl/HJJEQ1NLYEuzRjTA500AFS1CbgD5/BNEfCiqq4XkQdE5AqPRWcDi7TdaUUiMghnD+Kf7Tb9JxFZC6wF0oAHT/tdmGMMy0jg1Tsmc83ZA/j9hyV88/fL2LHX+t6NMcc66WmgPYmdBnrq36slayv5yeI1HG1s4WsTsrl56mCG9ov3d4nGmB6so9NAIwJRjOk+M8ZkMS43mcc/KObFwjJeKCzjkhEZ3HrBYCYOSg10ecaYALI9gCDQmT0AT3trj/Lssu08u6yU/fWNnDEgmVumDuHSkRmEhYk/SjXG9EAd7QFYAAQBfwVAq/qGJl4qLOfJpSWU7TvM4LQ4bp46mK9NyCYmMtwvr2GM6TksAIKYvwOgVVNzC2+s28UTH5awtqKGtPhobpg8iGvPHmjDShjTi1gABLGuCoBWqsqy4r08/mEJH26uJjYqnNlnDeCm8/PITu7TJa9pjOk+FgBBrKsDwNOGnQf5w0cl/G31ThT46tgsbpk6hJH9E7v8tY0xXcMCIIh1ZwC0qjhwmAVLt7Hosx3UNTRz/rA0br1gCOcN6dtWjzEmOFgABLFABECrmvpG/viv7Tz9cSl7ao8yqn8i/3bBEGaMziQi3O4nZEwwsAAIYoEMgFZHGpt5ZWUFT3xUQkl1HTkpffjOlDxmnZVLbJRdTmJMT2YBEMR6QgC0amlR3imq4okPSyjcvp/k2EiuO2cg1503iLR4ryN6G2MCzAIgiPWkAPD0+fZ9/P6fJbxdVEVUeBhXnZnDzecPZlBaXKBLM8Z4sAAIYj01AFoVV9fy5EclLP68gsaWFqaPyuSWqYOZMCAl0KUZY7AACGo9PQBa7T50hGc+KeW5Zds5eKSJSXmp/NvUwVyU38+GmjAmgCwAgliwBECr2qNNvLC8jKc+KmFnzRGG9Yvn5qmDmTm+P9ERNtSEMd3NAiCIBVsAtGpsbuH1NZU8/s9iNu46RFp8FNeeM5BrzxloHcbGdCMLgCAWrAHQSlX5pHgvTy3dxnsbdxMVEcaV4/tz45Q8CjLtCmNjuprdD8AEjIgweWgak4emUVxdy9Mfb2Px5xW8WFjO5KF9uWlKHhcOt34CY7qbT3sAIjIdeAQIB55U1V+1m/8wcJH7NBbop6rJ7rxmnNs+AuxQ1Svc9jxgEdAX+Bz4tqo2nKgO2wMIzj0Abw7UN/D8Z2U880kpuw4eYXB6HDdMzuMbZ2TbhWXG+NlpHwISkXBgM/AloBznJvFzVHVDB8vfCUxQ1Rvd57Wqetw9CEXkReAvqrpIRB4HVqvq/BPVYgHQewKgVWNzC0vWVrJg6TZWl9eQ1CeSOZMGMPe8gWQl2UikxvhDRwHgy2Auk4Ctqlri/oe+CJh5guXnAM+fpBgBLgZedpueAa70oRbTy0SGhzFzfDavfHcyL996LucN6csTHxZz/kPv873nV7K67ECgSzSm1/JlXzsbKPN4Xg6c7W1BERkI5AHveTTHiEgh0AT8SlVfwTnsc0BVmzy2mX1qpZveRESYOCiViYNSKdtXzzOflPLC8jJeW72TMwemcN25A7mooB+JMXajGmP8xd8HW2cDL6tqs0fbQFWtEJHBwHsishao8XWDInILcAvAgAED/Fqs6ZlyU2P52eUjueuSYbxUWM7Tn2zjrkWriAwXzs7ryyUj+jFtRAa5qbGBLtWYoOZLH8C5wDxV/bL7/KcAqvpLL8uuBL6rqp90sK2FwN+BxUA1kKmqTe1foyPWB9D7+gB80dyirNixn3eKqnhnQxXF1XUAFGQmMG1EPy4ZkcG4nGQ7i8iYDnSmEzgCpxN4GlCB0wn8LVVd3265AuBNIE/djYpIClCvqkdFJA1YBsxU1Q0i8hKw2KMTeI2q/u5EtVgAhGYAtLdtTx3vFlXxTlEVy0v309yipMVHM62gH5eMzGDK0DT6RNkVx8a06tSFYCIyA/gNzmmgC1T15yLyAFCoqq+5y8wDYlT1Ho/1zgN+D7TgdDj/RlWfcucNxulQTgVWAteq6tET1WEBYAHQ3oH6Bj7YVM07RVX8c1M1h442ER0RxpShaUwbkcG0Ef3ISIwJdJnGBJRdCRzELAB809DUwvLSfby9wdk7KN9/GICxOUlclN+Ps/NSmTAgxfYOTMixAAhiFgCnTlXZXFXr9BsUVbGq7ACqEBEmjM5O4uy8VM4alMrEQSkkx0YFulxjupQFQBCzAOi8msONrNi+n89K97F82z7WlNfQ0NwCQH5GAmflpXDWoFQm5aXaBWim17EACGIWAP53pLGZ1WUHWF66j39t28eK7fupa3DOXs5J6cOkvFQmDUrlrLxUBqfFtX0PjAlGFgBBzAKg6zU1t1BUeahtD2F56T721jlDU6XFRzFhQAojsxIZkZXIyKxEclP7WCiYoGEBEMQsALqfqlJcXcdyNxBWlR1g2946Wr8F8dERFGQmMMINhRFZCRRkJloHs+mRLACCmAVAz1Df0MSmXYcoqjxEUeVBiioPsnHXIWqPOiOaiEBe37i2QGgNh6ykGNtbMAFl9wMwppNioyKYMCDlmJvdt7Qo5fsPs8ENhKLKg6ypOMDrayvblkmOjWRMdhIXFzhXLdsQFqansD2AIGB7AMHn0JFGNu76Yk/hs2372oawyM9wh7AYmcF4G8LCdAM7BBTELAB6h9I9dW3XJXgOYXFxQTqXjMhgyrA0uxmO6RIWAEHMAqD3qalv5IPNu3mnaDcfbNrNoSPOEBaTh6ZxiQ1hYfzMAiCIWQD0bo3NLSzfto+33b2Dsn1fDGExrSCDS0b2Y2RWonUkm9NmARDELABCh6qyZXctb2+o4t2iKla6Q1ikJ0QzNjuJUf0TGdk/idHZiWQn27UIxjd2FpAxQUBEGJ6RwPCMBL570VD21B7lvY27WVa8l/U7a3h/025a3P8DkvpEMqp/IqPdYBjVP5G8tHjCrVPZ+Mj2AIKA7QGYVocbmtm46yDrdx5k/c4a1u90rkVoaHLGNeoTGc6IrARGuXsJo/onMSwjnugIu0AtlNkhoCBmAWBOpLG5ha27a78IhYqDbKg82HaBWmS4MLRfAqPdvYTR2UmMyEokLtoOAIQKC4AgZgFgTlVLi7JjXz3r3L2EdRU1bNh5sG18o9arlke5h49G93e+psTZ0Ni9kfUBGBNCwsKEQWlxDEqL4/Kx/QHnH4iqg0dZv7OGdRXO3sKK7fv52+qdbev1T4ppC4XWw0iZiTaURW9lAWBMiBARMpNiyEyKYdqIjLb2/XUNbKh09hJaDyO9U1TVNvBdalxUWyC0HkIamBprVzD3Aj4FgIhMBx7BuSfwk6r6q3bzHwYucp/GAv1UNVlExgPzgUSgGfi5qr7grrMQuACocde7XlVXdebNGGNOXUpcFJOHpjF5aFpbW93RprbO5nUVzh7DU0tLaGx2UiE+OqKts7k1HIZlxBMZHhaot2FOw0kDQETCgceALwHlwHIReU1VN7Quo6o/8Fj+TmCC+7QeuE5Vt4hIf+BzEXlLVQ+483+sqi/7560YY/wlLjqCMwemcubA1La2hqYWNlcdYsPOg219Cy8WllHv3kgnKjyM/MwEJxDcw0gjbIjsHs2XPYBJwFZVLQEQkUXATGBDB8vPAe4HUNXNrY2qulNEdgPpwIFO1GyMCYCoiDBGZycxOjuJWeQC0NyibNtTx/qdNW3B8Ob6XSxaXgZAmMCQ9HhG9k8kPzOBgswE8jMT6W9DZPcIvgRANlDm8bwcONvbgiIyEMgD3vMybxIQBRR7NP9cRO4D3gXuUdWjXta7BbgFYMCAAT6Ua4zpLuFhwtB+8QztF8/M8dmA09m8s+ZIW5/Chp01fLZtH6+u+qKzOSE6guGZCW2hMDzD+Zoca2chdSd/dwLPBl5W1WbPRhHJAp4D5qpqi9v8U2AXTig8AfwEeKD9BlX1CXc+EydOtPMgjenhRITs5D5kJ/fhy6My29prDjeyueoQm3Z98fj76p38+V9NbctkJEa3hUF+ZiL5GQkMy4gnJtIOI3UFXwKgAtz9PUeO2+bNbOC7ng0ikgi8Dtyrqp+2tqtq6x0zjorI08CPfC3aGBN8kvpEctagVM4a9EW/QuupqRt3HXRCwQ2IZ5Ztb7u6OUxgQGosua2PlFhyU/uQmxJLTkofUuOi7HDSafIlAJYDw0QkD+cP/2zgW+0XEpECIAVY5tEWBfwVeLZ9Z6+IZKlqpTjfuSuBdaf7Jowxwcnz1NQL8/u1tTc1t7B9Xz2bdh1i465DFO+upWx/PevWVrK/vvGYbcRFhZPjhoLzNZbclD5tgRFvVzx36KSfjKo2icgdwFs4p4EuUNX1IvIAUKiqr7mLzgYW6bGXq84CpgJ9ReR6t631dM8/iUg6IMAq4FY/vB9jTC8QER7GkPR4hqTHM2NM1jHzao82Ubav3nnsP0z5/nrK9jlflxXvpa7hmCPQpMRGMiA1lsHp8QxJj3O22y+egX1jQ36MJBsKIgjYUBDG+EZV2V/f6IaDEwxl++vZsbeekupadtYcaVu29dBSayC0hUN6fK8bEsOGgjDG9HoiQmpcFKlxUYzLTT5uft3RJrbtqaO4upbi3bUUVzvTH23d09bnAM7Vz56BMLSfcyprv4ToXtXfYAFgjAkZcdERbdcyeGpuUXYeOMzWdsHwTlFV2zUNAGnxUc4NeTyGxhgQxMNiWAAYY0JeeJi0dRpf5NEZDXCgvoEtu2tZX1HDup3O8BhPfFhCk3tnnoToCEa0DrXdP4lR2YkMTY8nIgiGxbAAMMaYE0iOjTru9NWjTc1s3lXrjKzqDovx/Gc7ONLoHEaKigijIDPhmAH0RmQl9LhOZwsAY4w5RdER4YzJSWJMzheHkppblJLqL27Ms67iIK+v2cnzn+0AnBvzjMxKZFxuMuNykhmXm8zgtLiAHj6ys4CCgJ0FZExwUlXK9x9mbUUNq8sPsLrsAGvLa9pOVU2IjmBsblJbIIzPTSYjMcbvddhZQMYY081EvuhbaL2eoblFKa6uZVWZEwiryw8c06eQmRjDuNwkJxBykhmTk0RCTGSX1GcBYIwx3Sg8TBie4QyAN2uiM8rOkcZm1u882BYIq8sO8Nb6KsC5feeQ9Hgev/YMhvZL8GstFgDGGBNgMZHhnDkwhTMHprS17a9rYE1FDavLDrCq7AD9uuDQkAVAELBj/8aEnpS4KC4Yns4Fw9O77DV6/omqxhhjuoQFgDHGhCgLAGOMCVEWAMYYE6IsAIwxJkRZABhjTIiyADDGmBBlAWCMMSEqqAaDE5FqYPtprp4G7PFjOf5m9XWO1dc5Vl/n9PT6BqrqcVeUBVUAdIaIFHobDa+nsPo6x+rrHKuvc3p6fR2xQ0DGGBOiLACMMSZEhVIAPBHoAk7C6uscq69zrL7O6en1eRUyfQDGGGOOFUp7AMYYYzxYABhjTIjqdQEgItNFZJOIbBWRe7zMjxaRF9z5/xKRQd1YW66IvC8iG0RkvYjc5WWZC0WkRkRWuY/7uqs+9/VLRWSt+9qFXuaLiDzqfn5rROSMbqwt3+NzWSUiB0Xk++2W6dbPT0QWiMhuEVnn0ZYqIm+LyBb3a0oH6851l9kiInO7sb7/FpGN7vfvryKS3MG6J/xZ6ML65olIhcf3cEYH657wd70L63vBo7ZSEVnVwbpd/vl1mqr2mgcQDhQDg4EoYDUwst0ytwOPu9OzgRe6sb4s4Ax3OgHY7KW+C4G/B/AzLAXSTjB/BvAGIMA5wL8C+L3ehXOBS8A+P2AqcAawzqPtv4B73Ol7gIe8rJcKlLhfU9zplG6q71Igwp1+yFt9vvwsdGF984Af+fD9P+HvelfV127+/wL3Berz6+yjt+0BTAK2qmqJqjYAi4CZ7ZaZCTzjTr8MTBMR6Y7iVLVSVVe404eAIiC7O17bj2YCz6rjUyBZRLICUMc0oFhVT/fKcL9Q1Q+Bfe2aPX/GngGu9LLql4G3VXWfqu4H3gamd0d9qvoPVW1yn34K5Pj7dX3VwefnC19+1zvtRPW5fzdmAc/7+3W7S28LgGygzON5Ocf/gW1bxv0lqAH6dkt1HtxDTxOAf3mZfa6IrBaRN0RkVPdWhgL/EJHPReQWL/N9+Yy7w2w6/sUL5OcHkKGqle70LiDDyzI95XO8EWePzpuT/Sx0pTvcQ1QLOjiE1hM+v/OBKlXd0sH8QH5+PultARAURCQeWAx8X1UPtpu9Auewxjjgt8Ar3VzeFFU9A7gM+K6ITO3m1z8pEYkCrgBe8jI70J/fMdQ5FtAjz7UWkXuBJuBPHSwSqJ+F+cAQYDxQiXOYpSeaw4n/++/xv0u9LQAqgFyP5zlum9dlRCQCSAL2dkt1zmtG4vzx/5Oq/qX9fFU9qKq17vQSIFJE0rqrPlWtcL/uBv6Ks6vtyZfPuKtdBqxQ1ar2MwL9+bmqWg+LuV93e1kmoJ+jiFwPXA5c44bUcXz4WegSqlqlqs2q2gL8oYPXDfTnFwF8HXiho2UC9fmdit4WAMuBYSKS5/6XOBt4rd0yrwGtZ1xcBbzX0S+Av7nHDJ8CilT11x0sk9naJyEik3C+R90SUCISJyIJrdM4nYXr2i32GnCdezbQOUCNx+GO7tLhf16B/Pw8eP6MzQVe9bLMW8ClIpLiHuK41G3rciIyHbgbuEJV6ztYxpefha6qz7NP6WsdvK4vv+td6RJgo6qWe5sZyM/vlAS6F9rfD5yzVDbjnCFwr9v2AM4PO0AMzqGDrcBnwOBurG0KzuGANcAq9zEDuBW41V3mDmA9zlkNnwLndWN9g93XXe3W0Pr5edYnwGPu57sWmNjN3984nD/oSR5tAfv8cIKoEmjEOQ59E06f0rvAFuAdINVddiLwpMe6N7o/h1uBG7qxvq04x89bfwZbz4rrDyw50c9CN9X3nPuztQbnj3pW+/rc58f9rndHfW77wtafOY9lu/3z6+zDhoIwxpgQ1dsOARljjPGRBYAxxoQoCwBjjAlRFgDGGBOiLACMMSZEWQAYY0yIsgAwxpgQ9f8B8VBsjy9yp3UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cell_data = {}\n",
    "for sample,paths in sample_paths.items():\n",
    "    filtered_cell_barcodes = gzip.open(paths['filtered_cell_barcodes_path']).read().decode('utf-8').split('\\n')\n",
    "    for cell_bc in filtered_cell_barcodes:\n",
    "        cell_data[(sample,cell_bc)] = {}\n",
    "\n",
    "for sample,cell_bc,umi,larry_bc in counts_filtered.keys():\n",
    "    if (sample,cell_bc) in cell_data:\n",
    "        if not larry_bc in cell_data[(sample,cell_bc)]:\n",
    "            cell_data[(sample,cell_bc)][larry_bc] = 0\n",
    "        cell_data[(sample,cell_bc)][larry_bc] += 1\n",
    "\n",
    "num_cells_with_barcode = np.zeros(20)\n",
    "for larry_bc_counts in cell_data.values():\n",
    "    if len(larry_bc_counts)>0:\n",
    "        num_cells_with_barcode[:np.min([20,np.max(list(larry_bc_counts.values()))])] += 1\n",
    "efficiency = num_cells_with_barcode / len(cell_data)\n",
    "plt.plot(range(20),efficiency)\n",
    "plt.plot([N_UMIS,N_UMIS],[np.min(efficiency),np.max(efficiency)],'-k',linewidth=2)\n",
    "plt.text(N_UMIS*1.1,np.max(efficiency)*.95,'UMI cutoff',fontsize=14)\n",
    "\n",
    "final_BCs = {}\n",
    "for k,larry_bc_counts in cell_data.items():\n",
    "    final_BCs[k] = '-'.join(sorted([k for k,v in larry_bc_counts.items() if v >= N_UMIS]))\n",
    "print('\\nFinal annotation has '+repr(len(set(final_BCs.values())))+' clones in '+repr(len([k for k,v in final_BCs.items() if len(v)>0]))+' cells')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd790104",
   "metadata": {},
   "source": [
    "## Write output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "a1516f21",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "594635"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output = []\n",
    "for sample,paths in sample_paths.items():\n",
    "    filtered_cell_barcodes = gzip.open(paths['filtered_cell_barcodes_path']).read().decode('utf-8').split('\\n')\n",
    "    for cell_bc in filtered_cell_barcodes:\n",
    "        output.append(sample+','+cell_bc+','+final_BCs[(sample,cell_bc)])\n",
    "open(output_prefix+'.larry_clones.csv','w').write('\\n'.join(output))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df70b660",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "base"
  },
  "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}