notebooks/GraphISS_pipeline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Graph-ISS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0. Specify source code paths and create results folder (if not exists already)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "\n",
    "n_threads = 5 # Intra tile parallelization, each imaging round running on a different thread\n",
    "graph_iss_path = \"\" # Enter absolute path to graph-ISS repository\n",
    "prePro_path = graph_iss_path+\"/prePro_pipeline/\"\n",
    "graph_decoding_path = graph_iss_path+\"/pgm_pipeline/\"\n",
    "\n",
    "dataset_folder = \"\" # Enter absolute path to dataset folder\n",
    "tagList_path = \"\" # Enter absolute path to tagList\n",
    "tagList = pd.read_csv(tagList_path, sep = \",\", usecols = [0,1], header = None, names = [\"Seq\",\"Gene\"])\n",
    "\n",
    "results_folder = \"\" # Enter absolute path to execution folder\n",
    "if not os.path.exists(results_folder):\n",
    "    os.makedirs(results_folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Create array of intput images\n",
    "The dataset folder should have the following data structure:\n",
    "```\n",
    "dataset_folder\n",
    "├── Base1 (1st imaging round)\n",
    "│   ├── chan0.tif (general stain/anchor channel)\n",
    "│   ├── chan1.tif (nuclei channels/DAPI)\n",
    "│   ├── chan2.tif (T channel)\n",
    "│   ├── chan3.tif (G channel)\n",
    "│   ├── chan4.tif (C channel)\n",
    "│   └── chan5.tif (A channel)\n",
    "├── Base2 (2nd imaging round)\n",
    "│   ├── chan0.tif (general stain/anchor channel)\n",
    "│   ├── chan1.tif (nuclei channels/DAPI)\n",
    "│   ├── chan2.tif (T channel)\n",
    "│   ├── chan3.tif (G channel)\n",
    "│   ├── chan4.tif (C channel)\n",
    "│   └── chan5.tif (A channel)\n",
    "├── Base3 (3rd imaging round)\n",
    "│   ├── chan0.tif (general stain/anchor channel)\n",
    "│   ├── chan1.tif (nuclei channels/DAPI)\n",
    "│   ├── chan2.tif (T channel)\n",
    "│   ├── chan3.tif (G channel)\n",
    "│   ├── chan4.tif (C channel)\n",
    "│   └── chan5.tif (A channel)\n",
    "├── Base4 (4th imaging round)\n",
    "│   ├── chan0.tif (general stain/anchor channel)\n",
    "│   ├── chan1.tif (nuclei channels/DAPI)\n",
    "│   ├── chan2.tif (T channel)\n",
    "│   ├── chan3.tif (G channel)\n",
    "│   ├── chan4.tif (C channel)\n",
    "│   └── chan5.tif (A channel)\n",
    "└── BaseN (n imaging round)\n",
    "    ├── chan0.tif (general stain/anchor channel)\n",
    "    ├── chan1.tif (nuclei channels/DAPI)\n",
    "    ├── chan2.tif (T channel)\n",
    "    ├── chan3.tif (G channel)\n",
    "    ├── chan4.tif (C channel)\n",
    "    └── chan5.tif (A channel)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "tif_files = []\n",
    "# r=root, d=directories, f = files\n",
    "for r, d, f in os.walk(dataset_folder):\n",
    "    for file in f:\n",
    "        if '.tif' in file:\n",
    "            tif_files.append(os.path.join(r, file))\n",
    "            \n",
    "tif_files.sort()\n",
    "\n",
    "img_CSV = results_folder+'img_CSV/'\n",
    "# save into file\n",
    "if not os.path.exists(img_CSV):\n",
    "    os.makedirs(img_CSV)\n",
    "    \n",
    "    pd.DataFrame({\"Key\":np.arange(len(tif_files)), \"File\":tif_files}).to_csv(img_CSV+'/out.csv', index=False, header=True, sep='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Estimate normalization values for each image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" norm.py\n",
    "Compute normalization intervals for each channel and round.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    sys.argv[1] : input csv file contaning an array of input images\n",
    "    sys.argv[2] : upper percentile value of 98th percentile distri-\n",
    "        bution of image patches for signal level estimation.\n",
    "    sys.argv[3] : number of running threads. Each thread run over \n",
    "        an image.\n",
    "    sys.argv[4] : output csv file where to store computed intervals\n",
    "    sys.argv[5] : number of random patches used to estimate norma-\n",
    "        lization intervals\n",
    "\"\"\"\n",
    "\n",
    "norm_folder = results_folder+\"norm/\"\n",
    "\n",
    "if not os.path.exists(norm_folder):\n",
    "    os.makedirs(norm_folder)\n",
    "    \n",
    "    %%capture log\n",
    "    %run $prePro_path'/norm.py' $img_CSV'/out.csv' 99 $n_threads $norm_folder'/out.h5' 50000\n",
    "    %store log.stderr > $norm_folder'log_strerr'\n",
    "    %store log.stdout > $norm_folder'log_stdout'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Register whole slides"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" mainReg.py\n",
    "Coarse registration.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    sys.argv[1] : Path to pre-processing library\n",
    "    sys.argv[2] : Output folder where to store aligned images\n",
    "    sys.argv[3] : Input csv file contaning an array of input images\n",
    "    sys.argv[4] : Number of resolution levels used for registration\n",
    "    sys.argv[5] : Type of registration (valid arguments: \"translati-\n",
    "        on\", \"rigid\", or \"affine\")\n",
    "    sys.argv[6] : String flag that enables sequencing rounds registra-\n",
    "        tion using maximum projected images of general stain and\n",
    "        nuclei channel if enabled. Genaral stain images only are used\n",
    "        if disabled. (valid arguments: \"Nuclei\": feature enabled,\n",
    "        other string: feature disabled)\n",
    "    sys.argv[7] : Type of registration procedure. (valid arguments:\n",
    "        \"DO1\" : if only the general stain of the first sequencing round\n",
    "        is available\n",
    "        \"DO\" : if a general stain image is available for each sequencing\n",
    "        round\n",
    "    sys.argv[8] : highest resolution level for multiresolution image reg-\n",
    "        istration\n",
    "    sys.argv[9] : lower resolution level for multiresolution image regi-\n",
    "        stration\n",
    "    sys.argv[10] : flag to enable BSpline registration after rigid\n",
    "\"\"\"\n",
    "\n",
    "registration_folder = results_folder+\"registration/\"\n",
    "\n",
    "if not os.path.exists(registration_folder):\n",
    "    os.makedirs(registration_folder)\n",
    "    os.makedirs(registration_folder+'/folder1')\n",
    "    \n",
    "    # Redirect stdout, stderr to file\n",
    "    f_out =open(registration_folder+'/log_strout','w')\n",
    "    f_err =open(registration_folder+'/log_strerr','w')\n",
    "    orig_stdout = sys.stdout\n",
    "    orig_stderr = sys.stderr\n",
    "    sys.stdout = f_out\n",
    "    sys.stderr = f_err\n",
    "    # Run registration\n",
    "    %run $prePro_path'/mainReg.py' $prePro_path $registration_folder'/folder1' $img_CSV'/out.csv' 8 'translation,rigid' 'Nuclei' 'DO' 512 4 'noBspline'\n",
    "    # Close log files\n",
    "    f_out.close()\n",
    "    f_err.close()\n",
    "    # Restore stdout, stderr\n",
    "    sys.stdout = orig_stdout\n",
    "    sys.stderr = orig_stderr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialize ipython cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, 1, 2, 3]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import ipyparallel as ipp\n",
    "rc = ipp.Client()\n",
    "rc[:].push(dict(prePro_path=prePro_path, graph_decoding_path=graph_decoding_path, results_folder=results_folder, norm_folder=norm_folder, n_threads=n_threads, tagList=tagList))\n",
    "rc.ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Progress bar\n",
    "import sys\n",
    "import time\n",
    "import ipyparallel\n",
    "import ipywidgets\n",
    "from IPython.display import clear_output, display\n",
    "\n",
    "def progress_bar(ar):\n",
    "    # We create a progress bar.\n",
    "    w = ipywidgets.IntProgress()\n",
    "    # The maximum value is the number of tasks.\n",
    "    w.max = len(ar.msg_ids)\n",
    "    label = ipywidgets.HTML()\n",
    "    box = ipywidgets.VBox(children=[label, w])\n",
    "    # We display the widget in the output area.\n",
    "    display(box)\n",
    "    # Repeat:\n",
    "    while not ar.ready():\n",
    "        # Update the widget's value with the\n",
    "        # number of tasks that have finished\n",
    "        # so far.\n",
    "        w.value = ar.progress\n",
    "        label.value = u'{index} / {size}'.format(\n",
    "                        index=ar.progress,\n",
    "                        size=len(ar)\n",
    "                    )\n",
    "    w.value = w.max"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Image tiling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AsyncResult: _push>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tile_size_x = 1024\n",
    "tile_size_y = 1024\n",
    "\n",
    "alltiles = results_folder + \"alltiles/\"\n",
    "register_tiles = alltiles + \"registered/\"\n",
    "\n",
    "rc[:].push(dict(alltiles=alltiles, registration_folder=registration_folder, register_tiles=register_tiles, tile_size_x=tile_size_x, tile_size_y=tile_size_y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def runTiling(reg_tif):\n",
    "    import os\n",
    "    import sys\n",
    "    \n",
    "    reg_tif_name = os.path.splitext(reg_tif)[0]\n",
    "    os.makedirs(register_tiles+reg_tif_name)\n",
    "    \n",
    "    # Redirect stdout, stderr to file\n",
    "    f_out =open(register_tiles+reg_tif_name+'/log_strout','w')\n",
    "    f_err =open(register_tiles+reg_tif_name+'/log_strerr','w')\n",
    "    orig_stdout = sys.stdout\n",
    "    orig_stderr = sys.stderr\n",
    "    sys.stdout = f_out\n",
    "    sys.stderr = f_err\n",
    "    # Run bioformat command line tool bfconvert\n",
    "    !bfconvert -tilex $tile_size_x -tiley $tile_size_y $registration_folder'/folder1/'$reg_tif $register_tiles$reg_tif_name/%x_%y.tif    \n",
    "    # Close log files\n",
    "    f_out.close()\n",
    "    f_err.close()\n",
    "    # Restore stdout, stderr\n",
    "    sys.stdout = orig_stdout\n",
    "    sys.stderr = orig_stderr\n",
    "    \n",
    "    return 0\n",
    "\n",
    "if not os.path.exists(alltiles):\n",
    "    os.makedirs(alltiles)\n",
    "    if not os.path.exists(register_tiles):\n",
    "        os.makedirs(register_tiles)\n",
    "            \n",
    "        # Run parallel\n",
    "        view =  rc.load_balanced_view()\n",
    "        reg_tifs = os.listdir(registration_folder+'/folder1')\n",
    "        ar = view.map(runTiling, reg_tifs)\n",
    "        progress_bar(ar)            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Tiles processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8fa813d50f004a88b1557170df0e300c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HTML(value=''), IntProgress(value=0, max=738)))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def runTiles(tile_tif):\n",
    "    import sys\n",
    "    import os\n",
    "    import glob\n",
    "    import pandas as pd    \n",
    "    tile = os.path.splitext(tile_tif)[0]\n",
    "    # Collect tile channels and rounds\n",
    "    if not os.path.exists(alltiles+tile):\n",
    "        os.makedirs(alltiles+tile)\n",
    "        tifs = glob.glob(register_tiles+'/*/'+tile_tif)\n",
    "        tifs.sort()\n",
    "        pd.DataFrame({'Key':[tile]*len(tifs), \"File\":tifs}).to_csv(alltiles+tile+'/out.csv', index=False, header=True, sep='\\t')\n",
    "            \n",
    "        \n",
    "    # Run pre-processing\n",
    "    prePro_folder = results_folder+tile+'_prePro/'\n",
    "    if not os.path.exists(prePro_folder):\n",
    "        os.makedirs(prePro_folder)\n",
    "        os.makedirs(prePro_folder+\"/folder1\")\n",
    "        \n",
    "        \"\"\" main.py\n",
    "            Performs tile registration, tile normalization using\n",
    "            precomputed values, tophat filtering, signal candida\n",
    "            -te detection and signal merging.\n",
    "\n",
    "            Parameters\n",
    "            ----------\n",
    "            sys.argv[1] : Path to pre-processing library\n",
    "            sys.argv[2] : Input csv file contaning an array of in-\n",
    "                put images\n",
    "            sys.argv[3] : Output folder where to store aligned im-\n",
    "                ages\n",
    "            sys.argv[4] : Output dataframe of signal candidate dete-\n",
    "                ction after merging\n",
    "            sys.argv[5] : Output dataframe of signal candidate dete-\n",
    "                ction before merging\n",
    "            sys.argv[6] : Number of thread used for parallelizing\n",
    "                sequencing rounds\n",
    "            sys.argv[7] : h maxima threshold\n",
    "            sys.argv[8] : csv file containing image normalization\n",
    "                intervals\n",
    "            sys.argv[9] : Output hdf5 object storing normalized im-\n",
    "                ages\n",
    "            sys.argv[10] : Output hdf5 object storing signal candi-\n",
    "                dates binary masks\n",
    "            sys.argv[11] : Type of registration procedure. (valid\n",
    "                arguments: \"DO1\" : if only the general stain of\n",
    "                the first sequencing round is available, \"DO\" : if\n",
    "                a general stain image is available for each sequen-\n",
    "                cing round\n",
    "        \"\"\"\n",
    "        tile_path = alltiles+tile+'/out.csv'\n",
    "        \n",
    "        # Redirect stdout, stderr to file\n",
    "        f_out =open(prePro_folder+'/log_strout','w')\n",
    "        f_err =open(prePro_folder+'/log_strerr','w')\n",
    "        orig_stdout = sys.stdout\n",
    "        orig_stderr = sys.stderr\n",
    "        sys.stdout = f_out\n",
    "        sys.stderr = f_err\n",
    "        # Run prePro\n",
    "        %run $prePro_path'/main.py' $prePro_path $tile_path $prePro_folder\"/folder1\" $prePro_folder\"/out1\" $prePro_folder\"/out2\" $n_threads 0.05 $norm_folder'/out.h5' $prePro_folder\"/out3\" $prePro_folder\"/candidates_max.h5\" \"DO\"\n",
    "        # Close log files\n",
    "        f_out.close()\n",
    "        f_err.close()\n",
    "        # Restore stdout, stderr\n",
    "        sys.stdout = orig_stdout\n",
    "        sys.stderr = orig_stderr\n",
    "        \n",
    "    # Predict probabilities of signal candidate detections\n",
    "    predProb_folder = results_folder+tile+'_predProb/'\n",
    "    if not os.path.exists(predProb_folder):\n",
    "        os.makedirs(predProb_folder)\n",
    "        \n",
    "        \"\"\" get_proba_DNN.py\n",
    "            Signal candidate probability predictions.\n",
    "\n",
    "            Parameters\n",
    "            ----------\n",
    "            sys.argv[1] : Path to signal probability prediction library\n",
    "            sys.argv[2] : Input dataframe of signal candidate detection\n",
    "                after merging\n",
    "            sys.argv[3] : Output dataframe of signal probability predictions\n",
    "            sys.argv[4] : (Optional) probability threshold for signal candid-\n",
    "                    date predictions\n",
    "        \"\"\"\n",
    "        # Redirect stdout, stderr to file\n",
    "        f_out =open(predProb_folder+'/log_strout','w')\n",
    "        f_err =open(predProb_folder+'/log_strerr','w')\n",
    "        orig_stdout = sys.stdout\n",
    "        orig_stderr = sys.stderr\n",
    "        sys.stdout = f_out\n",
    "        sys.stderr = f_err\n",
    "        %run $prePro_path'/network/get_proba_DNN.py' $prePro_path'/network' $prePro_folder\"/out1\" $predProb_folder\"/out1\"\n",
    "        # Close log files\n",
    "        f_out.close()\n",
    "        f_err.close()\n",
    "        # Restore stdout, stderr\n",
    "        sys.stdout = orig_stdout\n",
    "        sys.stderr = orig_stderr\n",
    "    \n",
    "    \n",
    "    # Graph-based decoding\n",
    "    GM_folder = results_folder+tile+'_GM/'\n",
    "    if not os.path.exists(GM_folder):\n",
    "        os.makedirs(GM_folder)\n",
    "\n",
    "        \"\"\" main.py\n",
    "            Graph-based decoding of signal candidates.\n",
    "\n",
    "            Parameters\n",
    "            ----------\n",
    "            sys.argv[1] : Path to decoding library\n",
    "            sys.argv[2] : Dataframe of signal probability predictions\n",
    "            sys.argv[3] : Output file stroring dataframe of decoded\n",
    "                barcodes\n",
    "            sys.argv[4] : Output csv file of decoded barcodes\n",
    "            sys.argv[5] : Number of threads for parallelization\n",
    "            sys.argv[6] : d_th parameter for building graph connected\n",
    "                components\n",
    "            sys.argv[7] : d_max parameter for building graph connected\n",
    "                components\n",
    "            sys.argv[8] : Output list of sets of signals detections contribu-\n",
    "                ting to decoded barcodes\n",
    "            sys.argv[9] : enable barcodes decoding guided by the list\n",
    "                of targeted barcodes (valid arguments, \"prior\": feature\n",
    "                enabled, other string: feature disabled)\n",
    "            sys.argv[10] : csv file of targeted barcode sequences\n",
    "            sys.argv[11] : Dataframe of signal candidate detection befo-\n",
    "                re merging\n",
    "            sys.argv[12] : hdf5 object storing normalized images\n",
    "            sys.argv[13] : Path to signal probability prediction library\n",
    "        \"\"\"\n",
    "        # Redirect stdout, stderr to file\n",
    "        f_out =open(GM_folder+'/log_strout','w')\n",
    "        f_err =open(GM_folder+'/log_strerr','w')\n",
    "        orig_stdout = sys.stdout\n",
    "        orig_stderr = sys.stderr\n",
    "        sys.stdout = f_out\n",
    "        sys.stderr = f_err\n",
    "        %run $graph_decoding_path'/main.py' $graph_decoding_path $predProb_folder\"/out1\" $GM_folder'/out2' $GM_folder'/out1' $n_threads 3 4 $GM_folder'/out3' blind $tagList $prePro_folder\"/out2\" $prePro_folder\"/out3\" $prePro_path'/network'\n",
    "        # Close log files\n",
    "        f_out.close()\n",
    "        f_err.close()\n",
    "        # Restore stdout, stderr\n",
    "        sys.stdout = orig_stdout\n",
    "        sys.stderr = orig_stderr\n",
    "        \n",
    "    # Add tile offset to decoded barcode coordinates\n",
    "    barcodes_folder = results_folder+'/barcodes/'\n",
    "    if not os.path.exists(barcodes_folder):\n",
    "        os.makedirs(barcodes_folder)\n",
    "\n",
    "    df = pd.read_csv(GM_folder+'/out1')\n",
    "    x, y = tile.split('_')\n",
    "    df.loc[:,'global_X_pos'] = df.loc[:,'global_X_pos'] + int(x) * tile_size_x\n",
    "    df.loc[:,'global_Y_pos'] = df.loc[:,'global_Y_pos'] + int(y) * tile_size_y\n",
    "    df.to_csv(barcodes_folder+tile+'_barcodes_withOffset.csv', index=False, header=True, sep=',')\n",
    "    return 0\n",
    "        \n",
    "# Run parallel\n",
    "view =  rc.load_balanced_view()\n",
    "tiles = os.listdir(register_tiles+os.listdir(register_tiles)[0])\n",
    "ar = view.map(runTiles, tiles)\n",
    "progress_bar(ar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Output file (barcodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "be9d0fd9b2064a9c92d174eb17f9f725",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=732), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "from tqdm.notebook import tqdm\n",
    "# Concatenate results in a single file\n",
    "barcodes_folder = results_folder+'/barcodes/'\n",
    "\n",
    "for i, tile_barcodes in enumerate(tqdm(os.listdir(barcodes_folder))):\n",
    "    if i==0:\n",
    "        barcodes_df = pd.read_csv(barcodes_folder+tile_barcodes)\n",
    "    else:\n",
    "        df = pd.read_csv(barcodes_folder+tile_barcodes)\n",
    "        barcodes_df = pd.concat([barcodes_df, df])\n",
    "\n",
    "# Add corresponding gene names to decoded barcodes\n",
    "d = pd.Series(tagList.Gene.values,index=tagList.Seq).to_dict()\n",
    "barcodes_df[\"Gene\"] = barcodes_df['letters'].map(d)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Auxiliary function for computing final quality score of decoded barcodes\n",
    "from sklearn.metrics import auc\n",
    "\n",
    "def runROC(exp_df,unexp_df):\n",
    "    ppv = []\n",
    "    tnr = []\n",
    "    tpr = []\n",
    "    fpr = []\n",
    "    n_bins = 100\n",
    "    fp_minQ = unexp_df.Q.min()\n",
    "    fp_maxQ = unexp_df.Q.max()\n",
    "    tp_minQ = exp_df.Q.min()\n",
    "    tp_maxQ = exp_df.Q.max()\n",
    "    Q_min = np.amin([fp_minQ,tp_minQ])\n",
    "    Q_max = np.amax([fp_maxQ,tp_maxQ])\n",
    "    step = (Q_max-Q_min)/n_bins\n",
    "    Q=Q_min\n",
    "    for n in range(n_bins):\n",
    "        tp = len(exp_df[exp_df.Q>=Q])\n",
    "        fp = len(unexp_df[unexp_df.Q>=Q])\n",
    "        tn = len(unexp_df[unexp_df.Q<Q])\n",
    "        fn = len(exp_df[exp_df.Q<Q])\n",
    "        ppv.append(tp/(tp+fp))\n",
    "        tnr.append(tn/(tn+fp))\n",
    "        tpr.append(tp/(tp+fn))\n",
    "        fpr.append(fp/(fp+tn))\n",
    "        Q = Q + step\n",
    "    return fpr, tpr, auc(fpr, tpr)\n",
    "\n",
    "# Find best quality threshold that miximize expected vs unexpected separation\n",
    "def find_d1_param(barcodes_df, tagList_df):\n",
    "    r = np.linspace(2,5,20) # 20 linearly spaced numbers\n",
    "    auc_list=[]\n",
    "    for v in tqdm(r):\n",
    "        def T_quality(x):\n",
    "            return np.clip(1-np.log(1+x)/v,0,1)\n",
    "        barcodes_df[\"Q\"]=barcodes_df.seq_quality_min*barcodes_df.max_dist.apply(T_quality)\n",
    "\n",
    "        exp_df = barcodes_df[barcodes_df.letters.isin(tagList_df.Seq)]\n",
    "        unexp_df = barcodes_df[~barcodes_df.letters.isin(tagList_df.Seq)]\n",
    "        unexp_tagList = unexp_df.letters.unique()\n",
    "        auc_list.append(runROC(exp_df,unexp_df)[-1])\n",
    "\n",
    "    return r[np.argmax(auc_list)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "002967260fa245aaaef582985c98fb02",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=20), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "d1 = find_d1_param(barcodes_df,tagList)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_quality(x):\n",
    "    return np.clip(1-np.log(1+x)/d1,0,1)\n",
    "barcodes_df[\"quality\"]=barcodes_df.seq_quality_min*barcodes_df.max_dist.apply(compute_quality)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save barcode file\n",
    "barcodes_df.to_csv(results_folder+'/barcodes.csv', index=False, header=True, sep=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot quality of decoded barcodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Number of reads')"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEGCAYAAACkQqisAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXjU5bn/8fdNAHFhU0GpqMFzQEGWGIMgyCYKKCpC4adIFXDvwVaxvwq2VlzaaltRD4haPLJoVRSoSBVbQQgIqCwaUEQWNdpUKgjKKh5C7vPHfDMGmEkmQyaTyXxe15UrM893mfsxmDvP+jV3R0REJB41kh2AiIikLiURERGJm5KIiIjETUlERETipiQiIiJxq5nsACrb8ccf75mZmckOQ0QkZaxcufJrd28U6VjaJZHMzExWrFiR7DBERFKGmX0e7Zi6s0REJG5KIiIiEjclERERiVvajYlEsm/fPgoKCti7d2+yQ5EUUqdOHZo2bUqtWrWSHYpI0iiJAAUFBdStW5fMzEzMLNnhSApwd7Zu3UpBQQHNmjVLdjgiSaPuLGDv3r0cd9xxSiASMzPjuOOOU+tV0p6SSEAJRMpL/2ZElEREROQwaEwkgkfmrq/Q+428sEWF3q+i5Ofns3TpUq666qpyXTds2DAuueQSBg4ceEB59+7deeihh8jJyanIMGOSn5/PJZdcwocffljpny2SztQSSWP5+fk8//zzyQ4DgMLCwmSHIBXgkbnrw1+SHpREqoi//OUvnHPOOWRlZXHTTTexf/9+Pv/8c5o3b87XX39NUVERXbp04Y033iA/P58zzjiDoUOH0rZtWwYOHMiePXsAWLlyJd26dePss8+md+/ebNq0CYCNGzdywQUX0K5dO7Kzs/nkk08YPXo0b731FllZWTzyyCPs37+fX/7yl7Rv3562bdvy5z//GQjNRLrlllto1aoVffv2ZfPmzaXWo1OnTrRu3Zply5YBsGzZMjp16sRZZ51Fp06dWLduHQBTpkxh0KBBXHrppfTq1QuAP/7xj7Rp04Z27doxevRoAPLy8ujYsSNt27alf//+fPPNN+G6tmvXjnPPPZcJEyaEY4hWj02bNtG1a1eysrJo3bo1b731VoX9/ETSlZJIFbB27VpefPFFlixZQl5eHhkZGTz33HOceuqpjBo1iptvvpmxY8fSqlWr8C/bdevWceONN7J69Wrq1avH448/zr59+/jZz37GjBkzWLlyJddeey2//vWvARgyZAgjRoxg1apVLF26lCZNmvDggw/SpUsX8vLyGDlyJE8//TT169dn+fLlLF++nKeeeorPPvuMl19+mXXr1vHBBx/w1FNPsXTp0qh12b17N0uXLuXxxx/n2muvBeCMM85g0aJFvP/++9x333386le/Cp//9ttvM3XqVObPn8/rr7/OrFmzePfdd1m1ahV33HEHANdccw1/+MMfWL16NW3atOHee+8FYPjw4YwbN4633377gBii1eP555+nd+/e5OXlsWrVKrKysiruhyiSpjQmUgW8+eabrFy5kvbt2wPw3Xff0bhxYwCuv/56pk+fzpNPPkleXl74mpNPPpnOnTsD8JOf/IRx48bRp08fPvzwQy688EIg9Bd5kyZN2LlzJ//617/o378/EFokF8kbb7zB6tWrmTFjBgDbt29nw4YNLFq0iMGDB5ORkcGPfvQjzj///Kh1GTx4MABdu3Zlx44dfPvtt+zcuZOhQ4eyYcMGzIx9+/aFz7/wwgs59thjAZg3bx7Dhw/nqKOOAuDYY49l+/btfPvtt3Tr1g2AoUOHMmjQoEPKr776al5//fVS69G+fXuuvfZa9u3bx+WXX64kUolKdm9V1TFCiY+SSBXg7gwdOpQHHnjgkGN79uyhoKAAgF27dlG3bl3g0OmlZoa7c+aZZx7yl/mOHTtijmP8+PH07t37gPI5c+bEPJ01Uly/+c1v6NGjBy+//DL5+fl07949fPzoo48+4PNj/ZzSzo1WD4BFixbx2muvcfXVV/PLX/6Sa665JqbPk4qjhFK9qDurCujZsyczZswIjzVs27aNzz8P7bw8atQohgwZwn333ccNN9wQvuaLL74IJ4sXXniB8847j9NPP50tW7aEy/ft28eaNWuoV68eTZs2ZdasWQB8//337Nmzh7p167Jz587wPXv37s0TTzwRbimsX7+e3bt307VrV6ZNm8b+/fvZtGkTCxYsiFqXF198EYDFixdTv3596tevz/bt2znppJOA0DhINL169WLSpEnh8Z1t27ZRv359GjZsGB6/ePbZZ+nWrRsNGjSgfv36LF68GIDnnnuuzHp8/vnnNG7cmBtuuIHrrruO9957L/oPRWKmwfT0ppZIBJX911GrVq347W9/S69evSgqKqJWrVpMmDCB/Px8li9fzpIlS8jIyGDmzJlMnjyZHj160LJlS6ZOncpNN91E8+bN+elPf0rt2rWZMWMGP//5z9m+fTuFhYXcdtttnHnmmTz77LPcdNNN3H333dSqVYvp06fTtm1batasSbt27Rg2bBi33nor+fn5ZGdn4+40atSIWbNm0b9/f+bPn0+bNm1o0aJFuAspkoYNG9KpUyd27NjBpEmTALjjjjsYOnQoDz/8cKldYX369CEvL4+cnBxq167NxRdfzO9//3umTp3KzTffzJ49ezjttNOYPHkyAJMnT+baa6/lqKOOOqDVcf3110esR25uLn/605+oVasWxxxzDM8880wF/QRF0pe5e7JjqFQ5OTl+8EOp1q5dS8uWLZMUUflpTUTVkWr/dhIhWguk5B9jsZwjVZeZrXT3iAvA1J0lIiJxU3dWCsrMzFQrRJIqlvEPjZGkB7VEREQkbkoiIiISNyURERGJm5KIiIjETQPrkSw4dOX4YelxZ8XeL4lyc3OpXbs2nTp1Ktd1mZmZrFixguOPP/6A8mOOOYZdu3ZVZIgxmzJlCitWrOCxxx5LyueLVAdKIlIuubm5HHPMMeVOIolQWFhIzZr6J1xZEjHbSlugpD51Z1UB+fn5tG7dOvz+oYce4p577gFCD3oaNWoU55xzDi1atAhv/xFtu/OXX36ZCy64AHdn06ZNtGjRgn//+99MmTKFfv360adPH04//fTwTrgQeRt6gL///e9kZ2fTrl07evbsSX5+Pk8++SSPPPIIWVlZvPXWW2zZsoUf//jHtG/fnvbt27NkyRIAtm7dSq9evTjrrLO46aabKG1R6y9+8Quys7Pp2bMnW7ZsAeCpp56iffv2tGvXjh//+MfhrVCGDRvG7bffTo8ePRg1ahS7du1i+PDhtGnThrZt2zJz5kwgtBVMmzZtaN26NaNGjQp/1uTJk8Or7otjBaLWY+HChWRlZZGVlcVZZ511wDYxIqIkkhIKCwtZtmwZjz76aPiXf7Ttzvv378+JJ57IhAkTuOGGG7j33ns58cQTgdBzPZ577jny8vKYPn06K1asiLoN/ZYtW7jhhhuYOXMmq1atYvr06WRmZnLzzTczcuRI8vLy6NKlC7feeisjR45k+fLlzJw5k+uvvx6Ae++9l/POO4/333+fyy67jC+++CJi3Xbv3k12djbvvfce3bp1C9dvwIABLF++nFWrVtGyZUuefvrp8DXr169n3rx5jB07lvvvv5/69evzwQcfsHr1as4//3y+/PJLRo0axfz588nLy2P58uXMmjWLTZs2MWbMGJYsWcLcuXP56KOPwveMVo+HHnqICRMmkJeXx1tvvcWRRx5Z8T9gkRSmvoAUMGDAAADOPvts8vPzgejbnTdr1ozx48fTunVrOnbsGN6aHULbrh933HHhey5evJiaNWtG3Ib+nXfeoWvXrjRr1gwgvF37webNm3fAL+MdO3awc+dOFi1axF//+lcA+vbtS8OGDSNeX6NGDa644gogtKV9cV0//PBD7rrrLr799lt27dp1wN5YgwYNIiMjI/z506ZNCx9r2LAhixYtonv37jRq1AgIPUtl0aJFAAeUX3HFFaxfv77UenTu3Jnbb7+dIUOGMGDAAJo2bRqxHiLpSkmkCqhZsyZFRUXh93v37j3g+BFHHAFARkZG+DGypW13/q9//YsaNWrw1VdfUVRURI0aoQZntO3jI21DP3v27Ji2ZS8qKuLtt9+O+Bd6rNu6R7pm2LBhzJo1i3bt2jFlyhRyc3PD55S1fXxpXWfRYopWj9GjR9O3b1/mzJlDx44dmTdvHmeccUZ5qyVSbak7qwo44YQT2Lx5M1u3buX777/n1VdfLfOaaNudFxYWMnz4cJ5//nlatmzJww8/HL5m7ty5bNu2je+++45Zs2bRuXPnqNvQn3vuuSxcuJDPPvssXA4csn18r169DpjdVPzgrK5du4a3Z3/99dfDj7Q9WFFRUbg19fzzz3PeeecBsHPnTpo0acK+ffsO2Ob9YAd//jfffEOHDh1YuHAhX3/9Nfv37+eFF16gW7dudOjQgdzcXLZu3cq+ffuYPn16mfX45JNPaNOmDaNGjSInJ4ePP/44aiwi6UgtkUgqeUpurVq1uPvuu+nQoQPNmjWL6S/daNudjx07li5dutClSxeysrJo3749ffv2BeC8887j6quvZuPGjVx11VXk5IQ25Yy0DX3Hjh2ZOHEiAwYMoKioiMaNGzN37lwuvfRSBg4cyCuvvML48eMZN24cI0aMoG3bthQWFtK1a1eefPJJxowZw+DBg8nOzqZbt26ccsopEetx9NFHs2bNGs4++2zq168ffh7J/fffT4cOHTj11FNp06ZN1AHtu+66ixEjRtC6dWsyMjIYM2YMAwYM4IEHHqBHjx64OxdffDH9+vUD4J577uHcc8+lSZMmZGdnhycRRKvHo48+yoIFC8jIyKBVq1ZcdNFF5fvhilRz2gqe9NjOW2siEiMd/u0US/SGipriW3UldSt4M8sws/fN7NXgfTMze9fMNpjZi2ZWOyg/Ini/MTieWeIedwbl68ysd4nyPkHZRjMbnei6iIjIgSpjTORWYG2J938AHnH35sA3wHVB+XXAN+7+n8AjwXmYWSvgSuBMoA/weJCYMoAJwEVAK2BwcK5EMGzYMLVCRKTCJXRMxMyaAn2B3wG3W2hqzPnAVcEpU4F7gCeAfsFrgBnAY8H5/YBp7v498JmZbQTOCc7b6O6fBp81LTj3h3ma5RBplo9IadKhK1jPBJGyJLol8ihwB1A8f/U44Ft3LwzeFwAnBa9PAv4JEBzfHpwfLj/ommjlhzCzG81shZmtKF4RXVKdOnXYunVrWvxSkIrh7mzdupU6deokOxSRpEpYS8TMLgE2u/tKM+teXBzhVC/jWLTySAkwYhZw94nARAgNrB98vGnTphQUFBApwYhEU6dOHS0+lLSXyO6szsBlZnYxUAeoR6hl0sDMagatjabAl8H5BcDJQIGZ1QTqA9tKlBcreU208nKpVatWeGW2iIjELmHdWe5+p7s3dfdMQgPj8919CLAAGBicNhR4JXg9O3hPcHy+h/qXZgNXBrO3mgHNgWXAcqB5MNurdvAZsxNVHxEROVQyFhuOAqaZ2W+B94HinfWeBp4NBs63EUoKuPsaM3uJ0IB5ITDC3fcDmNktwD+ADGCSu6+p1JqIiKS5Skki7p4L5AavP+WH2VUlz9kLDIpy/e8IzfA6uHwOMKcCQxURkXLQ3lkiIhI3JREREYmbNmAUkTAtLpTyUktERETipiQiIiJxU3eWiFQJJbvStC186lBLRERE4qYkIiIicVN3lkia04wsORxqiYiISNyUREREJG7qzhJJQ+rCkoqiloiIiMRNSUREROKmJCIiInFTEhERkbgpiYiISNyUREREJG5KIiIiEjclERERiZsWG4qkiVRaYKht4VOHWiIiIhI3JREREYmbkoiIiMRNYyIikjI0VlL1qCUiIiJxKzOJmNnRZlYjeN3CzC4zs1qJD01ERKq6WFoii4A6ZnYS8CYwHJiSyKBERCQ1xJJEzN33AAOA8e7eH2iV2LBERCQVxJREzOxcYAjwWlCmAXkREYkpGdwG3Am87O5rzOw0YEFiwxIRCUmllfbpqMwk4u4LgYUl3n8K/DyRQYlIxdAvYEm0qEnEzP4GeLTj7n5ZQiISEZGUUVpL5KHg+wDgROAvwfvBQH4CYxIRkRQRNYkE3ViY2f3u3rXEob+Z2aKERyYiIlVeLLOzGgWD6QCYWTOgUeJCEhGRVBFLEhkJ5JpZrpnlEpqZdVtZF5lZHTNbZmarzGyNmd0blDczs3fNbIOZvWhmtYPyI4L3G4PjmSXudWdQvs7Mepco7xOUbTSz0eWquYiIHLZYZmf93cyaA2cERR+7+/cx3Pt74Hx33xVsk7LYzF4HbgcecfdpZvYkcB3wRPD9G3f/TzO7EvgDcIWZtQKuBM4EfgTMM7PindcmABcCBcByM5vt7h/FWHcRETlMsW7A2Bw4HWhH6Bf7NWVd4CG7gre1gi8HzgdmBOVTgcuD1/2C9wTHe5qZBeXT3P17d/8M2AicE3xtdPdP3f1/gWnBuSIiUkli2YBxDDA++OoB/BGIaXqvmWWYWR6wGZgLfAJ86+6FwSkFwEnB65OAfwIEx7cDx5UsP+iaaOWR4rjRzFaY2YotW7bEErqIiMQglpbIQKAn8G93H06oNXJELDd39/3ungU0JdRyaBnptOC7RTlW3vJIcUx09xx3z2nUSHMCREQqSixJ5Dt3LwIKzaweoVbFaWVccwB3/xbIBToCDcyseCymKfBl8LoAOBkgOF4f2Fay/KBropWLiEgliSWJrDCzBsBTwErgPWBZWReZWaPgOszsSOACYC2h2V0Dg9OGAq8Er2cH7wmOz3d3D8qvDGZvNSM0PrMMWA40D2Z71SY0+D47hvqIiEgFiWV21n8FL580s78D9dx9dQz3bgJMNbMMQsnqJXd/1cw+AqaZ2W+B94Gng/OfBp41s42EWiBXBp+/xsxeAj4CCoER7r4fwMxuAf4BZACT3H1NTLUWEZEKUWYSCWZIDQFOc/f7zOwUMzvH3UttjQSJ5qwI5Z8SGh85uHwvMCjKvX4H/C5C+RxgTll1EBGRxIhlK/jHgSJCU3PvA3YCM4H2CYxLROKknXulMsWSRDq4e7aZvQ/g7t8UrzIXEZH0FsvA+r5gXMMhNGBOqGUiIiJpLpaWyDjgZaCxmf2O0MypuxIalYhIGUp22428sEUpZ0oixTI76zkzW0lowaEBl7v72oRHJiKHWvDAD6973Jm8OEQCpSYRM6sBrHb31sDHlROSiBygZOKIVq6EIklS6phIsFJ9lZmdUknxiIhICollTKQJsMbMlgG7iwv1jHUREYklidyb8ChEpMJ0/GJi+PU7p9yYxEgkHcQysL6wMgIRkcMQbdxEJMFiaYmISGVTUpAUEeuTDUVERA4RNYmY2ZvB9z9UXjgiIpJKSuvOamJm3YDLzGwaBz1J0N3fS2hkIulGXViSgkpLIncDowk9MfDhg445oV19RaQKePvTrckOQdJU1CTi7jOAGWb2G3e/vxJjEhGRFBHLFN/7zewyoGtQlOvuryY2LBERSQVlzs4ysweAWwk9nvYj4NagTERE0lws60T6AlnBPlqY2VRCz0bXjm8iVZxWr0uixbrYsAGwLXhdP0GxiKQfzciSFBdLEnkAeN/MFhCa5tsVtUJERITYBtZfMLNcoD2hJDLK3f+d6MBERKTqi6k7y903AbMTHItIekhSF5bGRyQRtHeWiIjETbv4iqQorVKXqqA8z1gXkWpCXVtSUfSMdRERiZuesS4iInHTM9ZFRCRuMT1j3cxOBZq7+zwzOwrISHxoItWIVqZLNRXLBow3ADOAPwdFJwGzEhmUiIikhljWiYwAOgM7ANx9A9A4kUGJiEhqiCWJfO/u/1v8xsxqEnqyoYiIpLlYBtYXmtmvgCPN7ELgv4C/JTYsEYlECwylqokliYwGrgM+AG4C5gD/k8igRKTyaOGhHI5YZmcVBQ+iepdQN9Y6d1d3loiIxDQ7qy/wCTAOeAzYaGYXxXDdyWa2wMzWmtkaM7s1KD/WzOaa2Ybge8Og3MxsnJltNLPVZpZd4l5Dg/M3mNnQEuVnm9kHwTXjzMzK/59ARETiFUt31ligh7tvBDCz/wBeA14v47pC4Bfu/p6Z1QVWmtlcYBjwprs/aGajCXWXjQIuApoHXx2AJ4AOZnYsMAbIIdQSWmlms939m+CcG4F3CHWz9YkhLpHKobUhkgZimZ21uTiBBD4FNpd1kbtvcvf3gtc7gbWE1pj0A6YGp00FLg9e9wOe8ZB3gAZm1gToDcx1921B4pgL9AmO1XP3t4PutWdK3EtERCpB1JaImQ0IXq4xsznAS4RaAoOA5eX5EDPLBM4iNK5yQvCQK9x9k5kVrzk5CfhnicsKgrLSygsilIuISCUprTvr0hKvvwK6Ba+3AA1j/QAzOwaYCdzm7jtKGbaIdMDjKI8Uw42Eur045RRtSCwiUlGiJhF3H364NzezWoQSyHPu/teg+CszaxK0QprwQ9dYAXByicubAl8G5d0PKs8NyptGOP8Q7j4RmAiQk5OjmWUi1cwjc9eHX4+8sEUSI0k/ZQ6sm1kz4GdAZsnzy9oKPpgp9TSw1t0fLnFoNjAUeDD4/kqJ8lvMbBqhgfXtQaL5B/D74llcQC/gTnffZmY7zawjoW6ya4DxZdVHJNVU5gJDrRmR8opldtYsQsngb0BROe7dGbga+MDM8oKyXxFKHi+Z2XXAF4TGWCA0u+piYCOwBxgOECSL+/lhHOY+d98WvP4pMAU4ktCsLM3MEhGpRLEkkb3uPq68N3b3xUQetwDoGeF8J7TZY6R7TQImRShfAejRvVJ1aFqvpJlYksh/m9kY4A3g++LC4um7IiKSvmJJIm0IdUudzw/dWR68FxGRNBZLEukPnFZyO3gRERGIbcX6KqBBogMREZHUE0tL5ATgYzNbzoFjIqVO8RWR1FZyui9Urym/WldScWJJImMSHoWIiKSkWJ4nsrAyAhERqQjRWhkly6XixLJifSc/7ElVG6gF7Hb3eokMTCRlJGBtiB6DK6kilpZI3ZLvzexy4JyERSQiIikjltlZB3D3WWiNiIiIEFt31oASb2vwwxMGRUQkzcUyO6vkc0UKgXxCTyEUEZE0F8uYyGE/V0REUl8qbhOvGVmJV9rjce8u5Tp39/sTEI+IiKSQ0loiuyOUHQ1cBxwHKIlI+tKW7yJA6Y/HHVv82szqArcSelDUNGBstOtEJD5aGyKpqNQxETM7FrgdGAJMBbLd/ZvKCExERKq+0sZE/gQMACYCbdx9V6VFJSIiKaG0lsgvCO3aexfwa7Pwk26N0MC6tj0RSVOpOFNLEqO0MZFyr2YXEZH0okQhIiJxi2XFuoiApvWKRKCWiIiIxE0tERE5LBpkT29KIiJJpAWGkuqURESkwqhVkn40JiIiInFTS0SkNJqRJVIqtURERCRuSiIiIhI3JREREYmbkoiIiMRNSUREROKm2VkilUwLDKU6UUtERETippaIyMG0NkQkZmqJiIhI3BKWRMxskpltNrMPS5Qda2ZzzWxD8L1hUG5mNs7MNprZajPLLnHN0OD8DWY2tET52Wb2QXDNOCvx/F4RSb6OX0wMf0n1lciWyBSgz0Flo4E33b058GbwHuAioHnwdSPwBISSDjAG6ACcA4wpTjzBOTeWuO7gzxIRkQRL2JiIuy8ys8yDivsB3YPXU4FcYFRQ/oy7O/COmTUwsybBuXPdfRuAmc0F+phZLlDP3d8Oyp8BLgdeT1R9REQieWTu+vDrkRe2SGIkyVHZA+snuPsmAHffZGaNg/KTgH+WOK8gKCutvCBCeURmdiOhVgunnHLKYVZBpPw0rVeqq6oyOyvSeIbHUR6Ru08EJgLk5OREPU/SlGZjJZyeM1J9VXYS+crMmgStkCbA5qC8ADi5xHlNgS+D8u4HlecG5U0jnC8iknAlu7CiladL11ZlT/GdDRTPsBoKvFKi/JpgllZHYHvQ7fUPoJeZNQwG1HsB/wiO7TSzjsGsrGtK3EtERCpJwloiZvYCoVbE8WZWQGiW1YPAS2Z2HfAFMCg4fQ5wMbAR2AMMB3D3bWZ2P7A8OO++4kF24KeEZoAdSWhAXYPqIpIQ0VoektjZWYOjHOoZ4VwHRkS5zyRgUoTyFUDrw4lRRCRR0qVrq6oMrItUrkoYTNeMLEkHSiIiUqk0U6t6URIREYlA4yCxURIREUmw6jw+ol18RUQkbmqJiFQgDaantsrowqpurRIlEUkf2t5EpMIpiYhI0qT7TK3q0CpREhGRtKZZWIdHSUREpApI1VaJkohUb1qZLpJQSiIiIlVMKrVKlEREpEpI90H2RKiMZKTFhiIiEje1RKT60ThIylOrJHUoiYiIVGFVfXxE3VkiIhI3tUSketCWJiJJoSQiEiONgySHxkd+UBW7ttSdJSIicVMSERGRuKk7S1KXxkHSjrq2fhCta6uyN5RUEhEphcZBJBUkcydiJRFJLWp9iFQpSiIiB1HrIzWoa6tq0MC6iIjETS0RqfoS3IWllodI/NQSERGRuKklIlWTBtClHKKNj5S3XMpPSUTSkrqwqq+SCSKWcjk8SiJSdWjsQyTlKIlIcilxiKQ0DayLiEjc1BKRyqfWh1QhFTlWko6D9EoiUjkSkDiULKSqScdZX0oikjhKHJLGYmnhVIdEk/JJxMz6AP8NZAD/4+4PJjmk9KNkIZK2UjqJmFkGMAG4ECgAlpvZbHf/KLmRVSNKECIJE621UlEtlAPv/1CF3PNgKZ1EgHOAje7+KYCZTQP6AYlJIrH8Qu1xZ0I+ulxi/MVf8pf5uacdV+Y5IlI5UmlhZKonkZOAf5Z4XwB0OPgkM7sRKE7tu8xsXZyfdzzwdemn/CrOW1dZMdS52km3OqdbfSEd63z92MOp86nRDqR6ErEIZX5IgftE4LBTu5mtcPecw71PKlGdq790qy+ozhUp1RcbFgAnl3jfFPgySbGIiKSdVE8iy4HmZtbMzGoDVwKzkxyTiEjaSOnuLHcvNLNbgH8QmuI7yd3XJPAjU2e0q+KoztVfutUXVOcKY+6HDCGIiIjEJNW7s0REJImUREREJG5KIiccHhEAAARfSURBVDEwsz5mts7MNprZ6GTHUxnMbJKZbTazD5MdS2Uws5PNbIGZrTWzNWZ2a7JjSjQzq2Nmy8xsVVDne5MdU2Uxswwze9/MXk12LJXBzPLN7AMzyzOzFRV6b42JlC7YWmU9JbZWAQZX961VzKwrsAt4xt1bJzueRDOzJkATd3/PzOoCK4HLq/PP2cwMONrdd5lZLWAxcKu7v5Pk0BLOzG4HcoB67n5JsuNJNDPLB3LcvcIXWKolUrbw1iru/r9A8dYq1Zq7LwK2JTuOyuLum9z9veD1TmAtoR0Rqi0P2RW8rRV8Vfu/Ks2sKdAX+J9kx1IdKImULdLWKtX6l0u6M7NM4Czg3eRGknhBt04esBmY6+7Vvs7Ao8AdQFGyA6lEDrxhZiuDbaAqjJJI2WLaWkWqBzM7BpgJ3ObuO5IdT6K5+353zyK028M5Zlatuy7N7BJgs7uvTHYslayzu2cDFwEjgu7qCqEkUjZtrZImgnGBmcBz7v7XZMdTmdz9WyAX6JPkUBKtM3BZMEYwDTjfzP6S3JASz92/DL5vBl4m1E1fIZREyqatVdJAMMj8NLDW3R9OdjyVwcwamVmD4PWRwAXAx8mNKrHc/U53b+rumYT+X57v7j9JclgJZWZHB5NFMLOjgV5Ahc26VBIpg7sXAsVbq6wFXkrw1ipVgpm9ALwNnG5mBWZ2XbJjSrDOwNWE/jLNC74uTnZQCdYEWGBmqwn9sTTX3dNiymuaOQFYbGargGXAa+7+94q6uab4iohI3NQSERGRuCmJiIhI3JREREQkbkoiIiISNyURERGJm5KISAUxs6Zm9oqZbTCzT83sMTM7Is575ZpZTvB6jpk1CL7+q2KjFjk8SiIiFSBYrPhXYJa7NweaA0cCfzzce7v7xcGK8gaAkohUKUoiIhXjfGCvu0+G0J5UwEjgGjO7xcweKz7RzF41s+7B6yfMbEVpz/MIngVxPPAg8B/BQsg/mdmzZtavxHnPmdlliauiyKFqJjsAkWriTELPIAlz9x3BHk2l/X/2a3ffFjy35k0za+vuq6OcOxpoHWyYiJl1I5SoXjGz+kAnYOhh1kOkXNQSEakYRuTdnSPtAl3S/zOz94D3CSWiVrF+oLsvBP7TzBoDg4GZwTY9IpVGSUSkYqwh9KS8MDOrR2jfoq0c+P9aneB4M+D/Az3dvS3wWvGxcngWGAIMBybHFbnIYVASEakYbwJHmdk1EH6s8ljgMeAzIMvMapjZyfywDXc9YDew3cxOIPSsh9LsBOoeVDYFuA0gHTYGlapHSUSkAnhoJ9P+wEAz20Co9VHk7r8DlhBKJB8ADwHFj+FdRagbaw0wKTivtM/YCiwxsw/N7E9B2VeEdpdWK0SSQrv4iiSAmXUCXgAGJPIpemZ2FKHklO3u2xP1OSLRqCUikgDuvtTdT01wAil+iNR4JRBJFrVEREQkbmqJiIhI3JREREQkbkoiIiISNyURERGJm5KIiIjE7f8ANnPyr/mHaUAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "exp_df = barcodes_df[barcodes_df.letters.isin(tagList.Seq)]\n",
    "unexp_df = barcodes_df[~barcodes_df.letters.isin(tagList.Seq)]\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(exp_df.quality,bins=100,alpha=0.5,label='expected barcodes')\n",
    "plt.hist(unexp_df.quality,bins=100,alpha=0.5,label='unexpected barcodes')\n",
    "plt.legend()\n",
    "plt.xlabel(\"Quality\")\n",
    "plt.ylabel(\"Number of reads\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}