SAM-parser.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Goal of this analysis is to:\n",
    "#1 Group reads based on barcode and perform filtering and quality control on barcode groups\n",
    "#2 Determine Amplicon Sequencing Bias to help adjust primer ratios\n",
    "#3 Generate Tables of Read Mappings for Each Barcode\n",
    "\n",
    "#Begin with imports\n",
    "import timeit\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "pd.set_option(\"display.precision\", 2) #floats displayed 3 digits\n",
    "cpallet = [\"#6CC860\",\"#8CC860\",\"#C4E54C\",\"#F7DD48\",\"#F5A232\",\"#FC8609\",\"#FE6315\"]*3\n",
    "import matplotlib.pyplot as plot\n",
    "plot.rcParams['svg.fonttype'] = 'none'\n",
    "import scipy.stats as stats\n",
    "from labellines import labelLine, labelLines\n",
    "\n",
    "header = [] #global var to keep track of the header of the SAM file\n",
    "\n",
    "#def some classes and basic functions for dealing with the SAM files\n",
    "class Read (object):\n",
    "    '''this reads class is specific to SAM files, contains a bunch of attribtues:\n",
    "ID(barcode), Sequence, Mapping position, Mapping quality, CIGAR string, and FLAG and align'''\n",
    "    def __init__ (self, ID = \"default\", seq = \"none\", MapPos = -1, MapQ = -1, CIGAR = \"none\", FLAG = \"none\", align = \"none\"):\n",
    "        self.ID = ID\n",
    "        self.seq = seq\n",
    "        self.MapPos = MapPos\n",
    "        self.MapQ = MapQ\n",
    "        self.CIGAR = CIGAR\n",
    "        self.FLAG = FLAG\n",
    "        self.align = align\n",
    "        self.species = \"N/A\"\n",
    "        return\n",
    "        \n",
    "    def __str__ (self):\n",
    "        return \"ID: %s \\t Seq: %s \\t MapPos: %d \\t CIGAR: %s \\t FLAG: %s \\t align: %s\" % (self.ID, self.seq, self.MapPos, self.CIGAR, self.FLAG, self.align)\n",
    "\n",
    "class Bargroup (object):\n",
    "    def __init__(self, bar = \"default\", reads = [], species = \"undet\", num_genes = -1):\n",
    "        self.bar = bar\n",
    "        self.reads = reads\n",
    "        self.species = species\n",
    "        self.num_genes = num_genes\n",
    "        return\n",
    "        \n",
    "\n",
    "#this function takes the SAM file as input and parses it into a dictionary of barcodes and reads\n",
    "def in_put1 (infile, limit = 0):      \n",
    "    bargroups = {}\n",
    "    count = 0\n",
    "    \n",
    "    for line in infile:  #each line is a new read\n",
    "        if line[0] == \"@\":        #ignore the line if it is a header\n",
    "            header.append(line)\n",
    "            continue       \n",
    "        \n",
    "        count += 1        \n",
    "\n",
    "        split = line.split(\"\\t\")    #split line into separate words\n",
    "        #now assign each segment of the split into a read\n",
    "        try:\n",
    "            #read = Read (ID = split[0][:15], align = split[2], seq = split[9], MapPos = int(split[3]), MapQ = split[4], CIGAR = split[5], FLAG = split[1])\n",
    "            read = Read (ID = split[0], align = split[2], seq = split[9], MapPos = int(split[3]), MapQ = int(split[4]), CIGAR = split[5]) \n",
    "\n",
    "        except:\n",
    "            print( split[3])\n",
    "            print( \"read could not be stored\")\n",
    "        \n",
    "        if read.ID[:15] in bargroups:\n",
    "            bargroups[read.ID[:15]].reads.append(read)\n",
    "        else:\n",
    "            bargroups[read.ID[:15]] = Bargroup(bar = read.ID, reads = [read])\n",
    "        \n",
    "        if limit != 0:\n",
    "            if count >= limit:\n",
    "                break\n",
    "    print( \"%d reads inputted\" % count)\n",
    "    return bargroups\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/flan3@ad.wisc.edu/SICMA-analysis/scripts/DoTA-seq/Bfragilis-paper\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['21-2-BF', '21-8-BF']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Set paths and get all the SAM files to be analyzed\n",
    "print(os.getcwd())\n",
    "samplename = \"BFragilis-summary\"\n",
    "os.chdir(\"/mnt/scratch/freeman/BF\")\n",
    "files = os.listdir(os.getcwd())\n",
    "files = [f for f in files if f[-4:] == \".sam\" and f.find(\"Unde\") == -1]\n",
    "files.sort()\n",
    "\n",
    "prefix = [f[:-4] for f in files]\n",
    "prefix #prefix is used to keep track of samples and for naming files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2586193 reads inputted\n",
      "113639 reads inputted\n",
      "it took 16.609799654921517 seconds\n"
     ]
    }
   ],
   "source": [
    "#import the reads into memory usin the in_put function\n",
    "tic = timeit.default_timer()\n",
    "\n",
    "libraries = []\n",
    "for f in files:\n",
    "    filename = f\n",
    "    with open(filename, \"r\") as infile:\n",
    "        libraries.append(in_put1(infile))\n",
    "toc = timeit.default_timer()\n",
    "\n",
    "print(f\"it took {toc-tic} seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "21-2-BF.sam the file contains 149522 bargroups\n",
      "21-8-BF.sam the file contains 11578 bargroups\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAF+CAYAAADay/61AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5wcZZn3/883k0AMSBATfDgHMQq4CkIW1PVZENyVgCuw4kqEgKAJqOADK6vw8wD+fFRkiauIGAhkIwkLKKAcDLAr4bC6oATkaEAzHCSghIOAcsjxev7oGug0Mz1V3VVdXT3f9+s1r0xX31V9VU9fubqq7rpvRQRmZmZVMarsAMzMzLJw4TIzs0px4TIzs0px4TIzs0px4TIzs0px4TIzs0px4eoikh6S9D5J/5+kc5NlkySFpNFlx2dWFc6l3ubC1YUi4usR8Yk0bSUdLuk2Sc9JWibptPrElHSMpMWSVkiaN8y21pd0nqSHJf1Z0q8lTR1mnYckvSjpL5L+JOmnkraqe36epJXJ8wM/H0mzb2btKiuXkvaTJC1M8uKPks5sVjSdS+m5cFXfOOA4YAKwO7A3cELd848B/xeYm2Jbo4FHgD2A8cCXgB9KmjTMev8QERsCmwGPA99teP60iNiw7ufiFLGYdVqeuQRwFrCcWl7sTC2vPjXMOs6lFFy4upCkUyQtaFh8pKTHJP1B0mcHFkbE9yPivyNiZUQ8ClwA/E3d85dFxE+Ap4Z73Yh4PiJOiYiHImJtRFwFPAjsmibuiHgJuATYMU17s6KVlUuJbYEfRsRLEfFH4BrgrWlWdC4158JVHe8FJgN/D5wo6X1DtPtb4N48XlDSG4A3p92epHHAR4Bb8nh9s4J0Kpe+AxwsaZykLYCp1IrXsJxLzblwVcdXkiOiu4F/B6Y1NpB0BDAFOL3dF5M0hto3zh9ExH3DNP+JpGeA54C/A/614fkTJD2T/DzZbmxmbepULt1I7QjrOWAZsBj4yTDrOJdScOGqjkfqfn8Y2Lz+SUkHAKcCUyMi1Qda0tV1F3kPqVs+CpgPrASOGa49cEBEbAysn7S/UdL/qnv+9IjYOPmZkG53zQpTeC4lOXQtcBmwAbXrZq8DvjlY+7pNOZdScOGqjq3qft+a2oViACTtA8yhdmH37rQbjIipdRd5L0i2JeA84A3AhyJiVbP2DdtbExGXAWuA92TcP7NO6UQubZK8zpkRsSIinqJ2dLfvEO0bt+dcasKFqzq+lJwrfytwBHAxgKS9qJ3S+1BE/KpxJUmjJY0F+oA+SWOHuY/l+8AO1BL3xSwBqmZ/at8sl2RZ16yDCs+l5EjtQeCTyXobA4cDd6YJ0LnUnAtXddwILAWuo3a64D+T5V+i1nV9Yd2ph6vr1vsi8CJwInBo8vsXB3sBSdsAR1HruvvHIU5lDOZKSX+hdl7+a8DhEZFLBxGzAhSeS4l/BPYBnkhebzVw/DCxOZdSkCeSNDOzKvERl5mZVYoLl5mZVYoLl5mZVYoLl5mZVYoLl5mZVUrl5qWZMGFCTJo0qewwzIZ02223PRkRE8uOYzjOJet2Q+VS5QrXpEmTWLx4cdlhmA1J0sNlx5CGc8m63VC55FOFZmZWKS5cZmZWKS5cZmZWKS5cZmZWKS5cZmZWKS5cZmZWKS5cZmZWKS5cZmZWKS5cZmZWKS5cZmZWKS5cZmZWKS5cZjnoX97PpxZ8io2O3QgmsGvZ8ZhVVZpcqtwgu2Z561/ez/EXH89Vd11FEOs8N3bMWLbaeCsefOpBVq9dDcCYvjH805R/4isf/ArbbbodV999NQfNPohVa1axas2qMnbBrPL6l/fzmQs/w8J7Fg7btrDCJWku8AFgeUT81SDPC/gOsC/wAvCxiLi9qHhs5Opf3s+XL/8yP7rtRy8Xlj71MaZvDC+tfqnpui+teonfPfG7dZatWrOKC355AZfcdglnffQsjr3oWF5Y+UJh8TuXrNddfffVHHjWgaxYvSJV+yKPuOYBZwLnD/H8VGBy8rM78P3kX7O2LFqyiCPnHcnDTw89u8iaWMOa1Wvaep0Vq1cwc/5ManWjUPNwLlmP6l/ezwHfO4CVa1amXqewwhURN0ma1KTJ/sD5ERHALZI2lrRZRPyhqJisNy1asojp503nsWcf6/hrr4k1NJxdzJ1zyXrZIecekqloQbnXuLYAHql7vCxZ9qpkkzQTmAmw9dZbdyS4kWLPPfcE4IYbbig1jrQW3LyAGefPGPYU3wjjXBphqpa3Qznhhyfwywd/mXm9MgvXYOdXBv3uGhHnAOcATJkypeDvt9YtTrr0JE695tSyw6gC55JVznEXHcd3rvtOS+uWWbiWAVvVPd4S6Py5HusKHzjjA/z07p+WHUZVOZesUhYtWdRy0YJyC9cVwDGSLqJ2IflZn5MfGWZdO4sTLjmh7DByM3rU6Je7ypfEuWSVcvSCo9tav8ju8BcCewITJC0DTgbGAETEbGAhte67S6l14T2iqFisO+x8ys7c+eidZYeRq9eMeQ1nTjuTYy86trD7uJxL1mt+t/x3wzdqoshehdOGeT6ATxf1+tZdxhw1puyjktytP3p9Lv3kpUx921T2eMse/Nt//Rvzb5nPczyX6+s4l8zW5SGfrHA7n7JzZYvW2DFjmbzpZEaPeuU73np963HI7odw71fuZerbpgKw3abbceYhZ/Lsd5+FJ7mtrHjNul3/8v62t+Ehn6xwVTs9uOXGW/KDI3/AXjvsVXYoZj3n+B8e3/Y2XLisUJpR+KgSmY3SKE770Gl89v2fLTsUsxHnqruuansbLlw2Ivgoyqw71C7JtseFywqzaMmi0l77fdu/j9nTZ7PdptuVFoOZrSvz/wlrGHSIHBcuK8ze39q78NeYsOEELp55sY+kzCrgw2d/ONsKz/P7wRa7cFmlHPW3RzF7+uyywzCzjBYtWcTTzz+dbaWV/HmwxS5c1rU2XH9DLv/05T6aMusB08+bnql9n/pYw+BTD7lwWdc4cZ8T+caHvlF2GGaWs/7l/ZmnHZq22zQWzFkw6HO+Adm6houWWW868KwDM69zygdPGfI5Fy4rxAaf2qDsEMysCyxasoi7H7070zrv3PadTXsEu3BZIV5Y9UKm9jHHU0OZ9aJWRoJf8InBTxEOcOEyM7PCZB0JfhSjhr3/0oXLzMwK0cqAuh/d/aPDtnHhMjOzQpx8xcmZ2o/SqKadMl5u12I8ZmZmTV1868WZ2s+ZPifVMG0uXGZmVois8/Ad+b+PTNXOhctK5x6FZr0njwkjh+LCZWZmudvvjP0ytd9/p/1Tt3Xhstx14+SRZtY5i5Ys4v7H78+0zqx/mpW6rQuXmZnlqpWbjrPMnefCZWZmucp603Gf+jK1d+EyM7NSTdttWqb2LlxmZlaqNDcd13PhMjOz3GTtBv/1A7+e6foWuHBZyQ7+64PLDsHMcpR1mKeT9j0p82u4cFmpLpx5YdkhmFmOsg7z1AoXLjMzy03WYZ5a4cJlZmaV4sJluVpwc/OZS83MBmQZ5qmeC5flavrc6WWHYGYlWbRkUab2WYZ5qufCZWZmucg61FPWbvADXLjMzCwXWYd6apULl5mZVYoLl5mZddw2m2zT8rouXGZm1rasQz3N/djcll/LhctKE3Oi7BDMLCfH//D4TO332mGvll/LhcvMzNp21Z1Xdey1XLjMzKxtQefOoLhwmZlZW7LeeNwuFy4zM2vLkT84MlP7Vod6GuDCZWZmbXn4qYcztW91qKcBLlxmZtayrN3gofWhnga4cFluRs3wx8lspMnaDb6dG48H+H8ay00nexWZWXfI2g2+nRuPB7hwmZlZy7J+YW3nxuMBLlxmZlYphRYuSftIul/SUkknDvL8eElXSrpT0r2SjigyHrMqch5Zr2i3G/yAwgqXpD7ge8BUYEdgmqQdG5p9GvhNROwE7AnMkrReUTGZVY3zyLpZ1h6F7XaDH1DkEdduwNKIeCAiVgIXAY3lNoDXShKwIfA0sLrAmMyqxnlkXStrj8J2u8EPKLJwbQE8Uvd4WbKs3pnADsBjwN3A/4mItY0bkjRT0mJJi5944omi4rUO8sjwqeWWR+Bcsnx1cmDdekUWLg2yrPF/q/cDdwCbAzsDZ0ra6FUrRZwTEVMiYsrEiRPzj9Sse+WWR+BcsnyVdQtMkYVrGbBV3eMtqX0jrHcEcFnULAUeBLYvMCazqnEemTUosnDdCkyWtG1yofhg4IqGNr8H9gaQ9AbgLcADBcZkVjXOI+sJeYyYMWB0bltqEBGrJR0DXAv0AXMj4l5JRyfPzwa+CsyTdDe1UyKfj4gni4rJrGqcR9atsvYozGPEjAGFFS6AiFgILGxYNrvu98eAvy8yBrOqcx5ZN8raozCPETMGeOQMMzPLrKweheDCZWZmLShzUG0XLjMzqxQXLsvF6BmFXi41swrLs0chuHBZTtawpuwQzKxL5dmjEFy4zMysYHn2KAQXLjMzy2jBzQtKfX0XLjMzy+STF3yy1Nd34TIzs0z+suIvpb6+C5eZmVWKC5d1nOfiMhs58u4KDy5cZmZWoLy7woMLl5mZZZC1R2HeXeHBhcvMzDIou0chuHCZmVkGZfcoBBcuMzOrGBcuMzMrRBE9CsGFy8zMUupf3p+pfRE9CsGFy8zMUjr+h8dnal9Ej0Jw4TIzs5SuvPPKskMAYMjZ/yTt0mzFiLg9/3DMeo9zyXrBoiWLyg7hZc2mrZ2V/DsWmALcCQh4O/BL4D3FhmbWM5xLVnlH/uDIskN42ZCnCiPivRHxXuBhYJeImBIRuwLvAJZ2KkCzqnMuWS94+KmHM7V/3/bvKyiSdNe4to+IuwceRMQ9wM6FRWSVM+2caWWHUBXOJRsxZk+fXdi2m50qHLBE0rnAAiCAQ4ElhUVklXPRrReVHUJVOJesklq5vrXdptsVEElNmsJ1BPBJ4P8kj28Cvl9YRGa9y7lklXTY3MMytS/yNCGkKFwR8ZKk2cDCiLi/0GjMephzyarq0WcezdS+yNOEkOIal6QPAncA1ySPd5Z0RaFRmfUg55KNFEWeJoR0nTNOBnYDngGIiDuASQXGZNarnEtWOd10/9aANIVrdUQ8W3gkZr3PuWSVk/X61v477V9QJK9I0znjHkkfBfokTQY+A/xPsWFZr4o5UXYIZXIuWeVkvb41659mDd+oTWmOuI4F3gqsAC4EngOOKzIosx7lXLKeV/T1LUjXq/AF4AvAFyT1ARtExEuFR2bWY5xLZvlI06vwPyRtJGkD4F7gfkn/UnxoZr3FuWS9ruj7twakOVW4Y0Q8BxwALAS2BqYXGpVZb3IuWU8r+v6tAWkK1xhJY6gl2+URsYracDVmlo1zySola1f4TlzfgnSF62zgIWAD4CZJ21C7qGxm2TiXrFKydoXvlGELV0ScERFbRMS+ERHA74H3Fh+aWW9xLlnVZO0K3ynNZkD+54ZFATwJ/DwiHiw0KrMe4lwyy1ezI67XNvxsRG321qslHdyB2Mx6hXPJLEdDHnFFxFcGWy5pE+BngCdhMkvBuWRVlLVjRieGehqQpnPGOiLiaUAFxGIVpBn+KLTKuWTdLGvHjE4M9TQgc+GStBfwpwJiMRtRnEvWzbJ2zOhUV3ho3jnjbl59j8kmwGNAd/aRNOtCziWzfDUbq/ADDY8DeCoini8wHrNe5Fwyy1GzzhkPdzIQs17lXLKq6V/en6l9p8YoHJD5GpeZmfW2o+Yflal9p8YoHFBo4ZK0j6T7JS2VdOIQbfaUdIekeyXdWGQ8Vq7x648vO4RKch5Zp11333WZ2neyYwakmI8rmYLhxYhYK+nNwPbA1ckAoc3W6wO+B/wdsAy4VdIVEfGbujYbA2cB+0TE7yVt2sa+WJd75sxnyg6hVK3kkvPI7NXSHHHdBIyVtAVwHXAEMC/FersBSyPigYhYSe0my8Y71D4KXBYRvweIiOVpAzeroFZyyXlk1iBN4VIyc+s/At+NiAOBHVOstwXwSN3jZcmyem8GXifpBkm3SRq0a7CkmZIWS1r8xBNPpHhps67USi7llkfgXLLhZR0xY/PxmxcUydBSFS5J7wIOAX6aLBv2FCODjwjQeC/LaGBXYD/g/cCXklMo664UcU5ETImIKRMnTkzx0mZdqZVcyi2PwLlkw8s6Ysb8j88vKJKhpSlAxwEnAT+OiHslvRG4PsV6y4Ct6h5vSe2Gy8Y2Tyb3szwv6SZgJ+C3KbZvVjWt5JLzyDoq64gZe+2wV0GRDG3YwhURNwI31j1+APhMim3fCkyWtC3wKHAwtXPx9S4HzpQ0GlgP2B34t3Shm1VLi7nkPDJr0GzIpytpMq14RHyw2YYjYrWkY4BrgT5gbvIt8+jk+dkRsUTSNcBdwFrg3Ii4p4X9MOta7eSS88js1ZodcZ2e/PuPwP8CFiSPp1GbfnxYEbEQWNiwbHbD438F/jXN9swqqq1cch5Zpyy4ecHwjep0esSMAc2GfLoRQNJXI+Jv6566MjmHbmYpOJesKmbOn5mpfadHzBiQplfhxOQiMgDJuXZ3RzLLzrlkXe3FVS9mat/pETMGpOlVeDxwg6QHkseTgGwDWVlP8iSSmTmXzHKQplfhNZImUxueBuC+iFhRbFhmvce5ZN0s6/WtsaPHFhTJ8NIccQFMBt4CjAV2kkREnF9cWGY9y7lkXenoBUdnaj/nsDkFRTK8NIPsngzsSW1omoXAVODngJPNLAPnknWz51dmm9f00HcdWlAkw0vTOeMgYG/gjxFxBLU78tcvNCqz3uRcMstBmsL1YkSsBVZL2ghYDrxxmHXM7NWcS2Y5SHONa3Ey388c4DbgL8CvCo3Kek6ZF3K7iHPJulL/8v5M7cu68XhAml6Fn0p+nZ0MK7NRRNxVbFjWa178frb7Q3qRc8m61VHzs92VUdaNxwOGPVWomkMlfTkiHgKekbRb8aGZ9RbnknWr6+67LlP7sm48HpDmGtdZwLuojasG8GdqU4mbWTbOJbMcpLnGtXtE7CLp1wAR8SdJ6xUcl1kvci6Z5SDNEdcqSX0k0zJImkht6gQzy8a5ZJU3bsy4skNIVbjOAH4MbCrpa9RumPx6oVGZ9SbnknWdRUsWZWp/9vSzC4okvTS9Ci+QdBu1GycFHBARSwqPzKzHOJesG3303MYJtZsrc8SMAanGKoyI+4D7ACRtLOkLEfG1QiMz60HOJes2jz/3eNkhZDbkqUJJW0k6R9JVkj4haZykWcBvgU07F6JZtTmXzPLV7IjrfOBG4FJgH+AW4F7g7RHxxw7EZl1s51N2LjuEKnEumeWoWeHaJCJOSX6/VtLjwF97/iADuPPRO8sOoUqcS9aVsnbM2Hz85gVFkk3Ta1ySXkftIjLAH4FxkjYAiIinC47NrGc4l6wbZe2YMf/j8wuKJJtmhWs8tYFA6+dnvz35N/Co1mZpOZesK2XtmLHXDnsVFEk2QxauiJjUwTjMepZzySxfaW5ANjOzHpN1KpNu4sJlZjYCTT9veqb2737juwuKJDsXLitczImyQzCzBjc/cHOm9ud//PyCIskuVeGS9B5JRyS/T5S0bbFhmfUm55JVVdlzcNVLM5HkycDngZOSRWOABUUGZdaLnEtm+UhzxHUg8EHgeYCIeAx4bZFBmfUo55J1hawdM8aOHltQJK1JU7hWRkTwyhxCGxQbklnPci5ZV8jaMWPOYXMKiqQ1aQrXDyWdDWwsaQbwM6C79sKsGpxL1hWydszohqlM6qWZj+t0SX8HPAe8BfhyRPxX4ZGZ9Rjnklk+hi1cko4HfuQEswF9M/rKDqGSnEtm+UhzqnAjaiNa/7ekT0t6Q9FBWXdby9qyQ6gq55KVbta1szK177aOGZCicEXEVyLircCngc2BGyX9rPDIzHqMc8m6wecv/Xym9t3WMQOyjZyxnNp0DE/hWVvN2uFcstKsiTWZ2ndbxwxIdwPyJyXdAFwHTABmRMTbiw7MrNc4l8zyMWznDGAb4LiIuKPoYMx6nHPJKqWP7uyINWThkrRRRDwHnJY83qT+ec/aapaOc8m6RdaOGd886JsFRdKeZkdc/wF8gNrMrcG6s7d61lZLxSPDA84l6xKfu/Rzmdp/9v2fLSiS9jSbAfkDyb8evdqsDc4l6xZrozduZUnTOeO6NMvMrDnnkpWpyjMeN2p2jWssMA6YIOl1vHJ6YyNq96CYWQrOJesGH5794UztR3XxPMPNrnEdBRxHLbFu45Vkew74XsFxmfUS55KV7teP/DpT+9MOOq2gSNrX7BrXd4DvSDo2Ir7bwZisi5106UnDN7J1OJesirq1YwakGx3+u5L+CtgRGFu3/PwiA7PudOo1p5YdQmU5l8zykWZ0+JOBPakl20JgKvBzwMlmloFzycqy4OYFmdp348C69dJcfTsI2Bv4Y0QcAewErJ9m45L2kXS/pKWSTmzS7q8lrZF0UKqozaqppVxyHlm7PnH+JzK178aBdeulKVwvRsRaYLWkjagNEDrsDZOS+qhdeJ5K7RvmNEk7DtHum8C1WQI3q6DMueQ8sjysWL0iU/tuHFi3XprCtVjSxtSmGL8NuB34VYr1dgOWRsQDEbESuAjYf5B2xwKXUktis17WSi45j8wapOmc8ank19mSrgE2ioi7Umx7C+CRusfLgN3rG0jaAjgQ2Av466E2JGkmMBNg6623TvHSZt2nxVzKLY+Sts6lESbr9a0qaHYD8i7NnouI24fZtgZZ1jhw3beBz0fEGmmw5slKEecA5wBMmTLFg99VhMcprGkzl3LLI3AujURZr29tumH3TxHX7Iir2TDCQe3bXTPLgK3qHm8JPNbQZgpwUZJsE4B9Ja2OiJ8Ms22zKmknl5xH1pas17cunHlhQZHkp9kNyO9tc9u3ApMlbQs8ChwMfLThNV4edFTSPOAqJ5v1mjZzyXlkHbXXDsMdk5QvzX1chw22fLibJiNitaRjqPVy6gPmRsS9ko5Onp/dQrxmldVKLjmPrB2LliwqO4RCpJkBuf5i71hq96HcToqbJiNiIbUbLeuXDZpoEfGxFLGYVVlLueQ8slZ95OyPZGq/ybhNhm/UBdL0Kjy2/rGk8cD8wiKyrqUZzS/8W3POJeu0J59/MlP7Hx39o4IiyVcr49a/AEzOOxCzEci5ZF2lCte3IN01rit5pfvtKGp37/+wyKDMepFzyTqpF+/fGpDmGtfpdb+vBh6OiGUFxWPWy5xL1jFZ799ar2+9giLJX5prXDcCJGOrjU5+3yQini44NrOe4lyyTsp6/9Z5h59XUCT5S3OqcCbwVeBFYC21O/mDFAPtmtkrnEvWzbp9YN16aU4V/gvw1ojI1j3FzBo5l6wj+pf3lx1CodL0Kuyn1vvJLDWPUzgo55J1xF6nZ+sdWJX7twakOeI6CfgfSb8EXj5pGhGfKSwqs97kXLLC9S/v5/d/+n2mdapy/9aANIXrbGARcDe18/Jm1hrnkhXuw7M/nHmdqty/NSBN4VodEf9ceCRmvc+5ZIX79SO/LjuEwqW5xnW9pJmSNpO0ycBP4ZFZV/FwT7lwLlnXGdXSAErlSnPENTCFwkl1y9yF1yw755IV6qRLTxq+UYPTDjqtgEiKleYG5G2Ha2Nmw3MuWdFOvebUzOt89v2fLSCSYhU2H5eZrcu5ZN3m9INOH75RFyp0Pi4zW4dzyQrTymnCKh5tgefjMusY55IVqZXThFXl+bgsdx41IzXnkpXmqL89quwQWub5uGxY7gqfD+eSdZPZ02eXHULLPB+XWec4l8xyMGThkvQm4A0DcwjVLf/fktaPiN4eftgsJ84lK9qsa2dlav/Wzd5aUCSd0ewa17eBPw+y/MXkOTNLx7lkhfrcJZ/L1P7yYy4vKJLOaFa4JkXEXY0LI2IxMKmwiMx6j3PJCrU245jN2226XUGRdEazwjW2yXOvyTsQsx7mXLLCLLh5QdkhdFyzwnWrpBmNCyV9HLituJDMeo5zyQpz+NzDyw6h45r1KjwO+LGkQ3gluaYA6wEHFh2YVZPv4RqUc8kKsWjJosynCcePHV9QNJ0zZOGKiMeBd0t6L/BXyeKfRsSijkRm1iOcS1aUA8/K/r3nsk9dVkAknZVmyKfrges7EIt1Id98nB/nkuXtuZeey7xO1WY7Hkz1ZhAzMzMeeOKBzOuMGTWmgEg6z4XLzKyCHnn6kczrzP3Y3AIi6TwXLjOzivnNY79pab1D33VozpGUw4XLzKxC+pf388Sfn8i8XlUnjRyMC5eZWYXs8+19WlqvqpNGDsaFy4aUtUfhxA0mFhSJmQ1Y+sTSzOv0wr1b9Vy4LDfLv7287BDMelqrwzv1wr1b9Vy4zMwqYvrc6ZnXGTNqTE/cu1XPhcvMrAJOuvSklta75rhrco6kfC5cNqit/2XrskMwszqnXnNq5nWmbDOl5462wIXLhvDIM9lvbjSz7nLrF28tO4RCuHBZLkZr2GEvzaxFbzv5bWWH0FVcuCwXq85ZVXYIZj3pHV95B/c8dk/m9XrphuNGLlz2Kh4R3qw7zLp2FncsuyPzekI9dcNxIxcuM7MudcIlJ7S03s/++Wc5R9JdXLjMzLpQ//L+ltftxZ6E9Vy4rG0xJ8oOwazn7Pa13coOoWu5cNk6+mb0lR2C2Yi3aMkinn7h6ZbWPXGfE3OOpvu4cNk61rK27BDMRrSr776avb+1d0vrjh87nm986Bs5R9R9Ci1ckvaRdL+kpZJe9TVA0iGS7kp+/kfSTkXGY1ZFzqORo395P/uesW9L647WaJ757jM5R9SdCitckvqA7wFTgR2BaZJ2bGj2ILBHRLwd+CpwTlHx2PBa6Qbv61vFch6NLNt/afuW1x1J91IWecS1G7A0Ih6IiJXARcD+9Q0i4n8i4k/Jw1uALQuMx6yKnEcjxKxrZ7F67eqW1t3vbfvlHE13K7JwbQHUD3i3LFk2lI8DVw/2hKSZkhZLWvzEE9mnrDarsNzyCJxL3azVe7YArvrMVTlG0v2KLFyDnXca9LySpPdSS7jPD/Z8RJwTEVMiYsrEiZ5l10aU3PIInEvdato501ped6QdbQEUOTLqMmCrusdbAo81NpL0duBcYGpEPFVgPNaEr291LefRCOctOqgAAA/oSURBVHDRrRe1vO5IO9qCYo+4bgUmS9pW0nrAwcAV9Q0kbQ1cBkyPiN8WGItZVTmPelw7Y4POP3J+jpFUR2FHXBGxWtIxwLVAHzA3Iu6VdHTy/Gzgy8DrgbMkAayOiClFxWSDu/H+G2Fy2VHYYJxHva2dorXpRpty6LsOzTGa6ih0EqWIWAgsbFg2u+73TwCfKDIGK4ZPE3aO86g3bfCpDVped/y48eyw2Q45RlMtHjljhLvx/hvLDsFsxJkxbwYvrHqh5fV33mrnHKOpHheuEWz0DM9abNZp/cv7OfcX57a8/ki9rlXPhWsEW8Oaltbb7LWb5RyJ2cjxpi+8qa31R+p1rXouXCNUOxeFH/vWq3pjm1kK4z89vq31fbRV48I1ArVTtMaNGZdjJGYjw9Hzj0YzxHMrn2t5GztutqOPthK+yDHCtFO0AJ4/6/mcIjEbGUbPHM2aaO20/IAJG0zg3v//3pwiqj4fcVlq1/3zdWWHYFYZ7/zaO9EMtV20xo8dzxPf9riS9XzENYK0e7S11w575RSJWW9rN9fqjZQ5trJw4Roh2k0k33Bs1twHzvgAP737p7lu03k3OBeuEaDdouWeTGbNjZoxihh80P6WuWgNzde4elwepyzck8ns1RYtWYRmCM1Q7kVr+ze0PhPySOAjrh6WR9Hytz6zdc26dhafu/RzrI21hb3Gkv+7pLBt9wIXrh7lomWWr9d/5vU8/eLThb+O8254Llw9KM8eTWYj2Ws//Vr+svIvHXs9F610XLh6TF5FywlkI80OX9yB+x6/r5TX3vPNe3L9v1xfymtXkQtXD3HRMktn1rWz+MKPv8CKNStKjWP+kfPd+akFLlw9YNGSRez9rb1z2ZaLlvWq9Y9en5VrVpYdBgCTJ07mt1//bdlhVJYLV8XleT3LRcuqbuNjNubZFc+WHUZTPspqnwtXReXdAcNFy7pZp3r0FWmnLXbijlPuKDuMnuDCVTFF9Bjc4y175L5Ns0bNPrsDX5x6sUfsJuM24anvPFV2GD3FhasiikpoFy1rR5piNFy7NM9X0fp96/O1A7/GZ9//2bJD6TkuXF2uFxPausNtD9+2zucr6+niNMUo5sSI+wyPX388z5zpEd2L5LEKu9DA+GdFJ7yva1m9LJ+3tG1HUtHa7237EXPCRasDfMTVJYqYEqEZFy0bzMBRkqWz+fjNuelzN7HdptuVHcqI4sLVBTr5rdT/KZm1ZtyYcTx/1vNlh2G4cJWq06dRXLTMsnvj699I/6n9ZYdhdVy4OqjM8/0uWmbD81FVNbhwdYALlln32u9t+3HVZ64qOwzLwIUrR93Wg8pFy7rBen3rlT5G4O6TdueWL9xSagyWHxeuNnVbsQIXLOsuK2avSJUnrdzz5etPI5MLVxPdWJSaccGyTklbZAY+kzEnhhydfb2+9Vgxe8U67c2aceGiegWqkZPdyjBc8Wr8XA4UJ7N29XzhGi6xqly0XLCsbP4MWhl6unD14sCe/o/CzEa6yhWutAODVrEoDcXFyszsFZUfZLeXCtRgXLTMzNZV+cIFvVm8Yk64aJmZDaJypwp7lYuUmVk6Llwlc8EyM8vGhatDXKDMzPLhwlUgFyszs/y5cLXIRcnMrBwuXMNwgTIz6y49W7iyDgJqZmbV0BP3cQ1luKLkomVmVj09e8Q1wMXJzKy39PQRl5mZ9Z5CC5ekfSTdL2mppBMHeV6Szkiev0vSLkXGY1ZFziOzdRVWuCT1Ad8DpgI7AtMk7djQbCowOfmZCXy/qHjMqsh5ZPZqRR5x7QYsjYgHImIlcBGwf0Ob/YHzo+YWYGNJmxUYk1nVOI/MGhTZOWML4JG6x8uA3VO02QL4Q30jSTOpfZOEDfMOc2S74YYbyg7Bmsstj8C51CtGet4WecQ12E1UjV380rQhIs6JiCkRMYWxucRmVhW55RE4l6w3FFm4lgFb1T3eEnishTZmI5nzyKxBkYXrVmCypG0lrQccDFzR0OYK4LCkV9Q7gWcj4lWnN8xGMOeRWYPCrnFFxGpJxwDXAn3A3Ii4V9LRyfOzgYXAvsBS4AXgiJZeyzcZW4/qZB6Bc8mqQRHV+qBqooIDXnnsRLNuI+m2iJhSdhzDcS5Ztxsqlyo35NOu2+zK4jmLyw7DrPKcS1ZVHvLJzMwqxYXLzMwqxYXLzMwqxYXLzMwqxYXLzMwqxYXLzMwqxYXLzMwqxYXLzMwqxYXLzMwqxYXLzMwqxYXLzMwqxYXLzMwqpXqjw0t/Bu5vcfXxwLMtthlseeOy+sfD/T4BeDJV1OljTNMmj/2oX1bkfuS1D/WPO/G32CYiJra4vY5pM5cGk+ZzmaVtO3//wZYN9blu5+8/mCzvQ5r2I/l9GDyXIqJSP8DiNtY9p9U2gy1vXFb/eLjfq7wfDcsK24+89qFJ7B35W3TrTzv73e57UfRneLjPRF6f4Tw+E0XmQS+9D/U/I+1U4ZVttBlseeOyKzP+3qqy9yOPfUiznbz2of5xGX+LkSLLe1H0Z3iwZXn/7YeSddtF5sFgy6r6PrysiqcKF0cFJukbjveje/TCPrRipO53I78PNVV6H6p4xHVO2QHkxPvRPXphH1oxUve7kd+Hmsq8D5U74jIzs5GtikdcZmY2grlwmZlZpbhwmZlZpVS+cEnaQNIPJM2RdEjZ8bRK0hslnSfpkrJjaZWkA5K/w+WS/r7seFolaQdJsyVdIumTZcfTCb2SR3nohVzMQ1fnc543nOV449pcYDlwT8Pyfajd6b8UODFZNh34h+T3i8uOvdX9qHvukrLjzmEfXgecV3bsOezHqG7bj6L2uZvzqMS/f1flYonvQ9flc7cecc2j9ka+TFIf8D1gKrAjME3SjsCWwCNJszUdjDGNeaTfj241j+z78MXk+W4yjwz7IemDwM+B6zobZq7m0Rt5lId5VD8X8zCPHsjnrixcEXET8HTD4t2ApRHxQESsBC4C9geWUUs66LL9ybgfXSnLPqjmm8DVEXF7p2NtJuvfIiKuiIh3A5U9bdYreZSHXsjFPPRKPlfpA7oFr3wjhFqibQFcBnxI0vepxvA7g+6HpNdLmg28Q9JJ5YSW2lB/i2OB9wEHSTq6jMAyGupvsaekMySdDSwsJ7TC9Eoe5aEXcjEPlcvn0WUHkIEGWRYR8TxwRKeDacNQ+/EU0FUfjiaG2oczgDM6HUwbhtqPG4AbOhtKx/RKHuWhF3IxD5XL5yodcS0Dtqp7vCXwWEmxtKMX9qMX9gF6Zz+yGIn7PBS/FzWVex+qVLhuBSZL2lbSesDBwBUlx9SKXtiPXtgH6J39yGIk7vNQ/F7UVO99KLtb4xBdNi8E/gCsovZt4OPJ8n2B3wL9wBfKjnMk7Ecv7EMv7Yf32e+F34fwILtmZlYtVTpVaGZm5sJlZmbV4sJlZmaV4sJlZmaV4sJlZmaV4sJlZmaV4sIFSFoj6Q5Jd0q6XdK7y46pW0h6SNKEFtf9mKQzB1l+iqQTWtzmqGQcwXsk3S3pVknbJs8tlLRxK9u1fDiXhuZcyk+Vxios0osRsTOApPcD3wD2SLOiJAGKiLUp2o6OiNVtRdqiMl87Zx8BNgfeHhFrJW0JPA8QEfuWGpmBc6lKKptLPuJ6tY2APwFI2lDSdck3x7sl7Z8snyRpiaSzgNuBrSR9XNJvJd2QzBp6ZtJ2nqRvSboe+KakTST9RNJdkm6R9Pak3TrfnJJvQZOSn/tUm532LtVm5R2XtDlV0m+S5ac37kiyzXMk/SdwfrKt/0725+Vvw8lo6Dck275P0gXJfyL123qNpGskzUgeHyrpV8m367NVm9MHSUck78ONwN80eZ93krRI0u/qtjl/4D1OHl+g2rxY9TYD/jDwn1tELIuIgb/XQ5ImSDo6iesOSQ8m7z2S/l7Szcm+/0jShk3is/Y5l5xLxSh76I5u+KE2cd4dwH3As8CuyfLRwEbJ7xOozQ4qYBKwFnhn8tzmwEPAJsAY4L+BM5Pn5gFXAX3J4+8CJye/7wXckfx+CnBCXUz3JK8zCQjgb5Llc4ETkte6H14e/WTjQfbrFOA24DXJ43HA2OT3ycDi5Pc9k/3ektqXmZuB9yTPPZTE8DPgsGTZDtSmvhiTPD4LOIxaIvwemAisB/xi4H0YJK47gdck7+sjyXu4B/CTpM144EFgdMO6WyYx3QHMAt5R99xDwIS6xwN/i39IXucmYIPkuc8DXy77s9drPziX9sS5VPiPj7hqXoyInSNie2qzg56ffEsS8HVJd1H7sG0BvCFZ5+GIuCX5fTfgxoh4OiJWAT9q2P6PImJgVtn3APMBImIR8HpJ44eJ75GI+EXy+4JkG88BLwHnSvpH4IUh1r0iIl5Mfh8DzJF0dxJj/Synv4raN6611D7Ik+qeuxz494g4P3m8N7ArcKukO5LHbwR2B26IiCeiNiHdxU326fKIeDEingSuB3aLiBuBN0naFJgGXBoNp2QiYhnwFuAkav/hXSdp7yFe4zvAooi4Enhnsr+/SGI+HNimSXzWGueSc6lwvsbVICJuVu0C6kRqA09OpPatcZWkh4CxSdPn61YbbD6besO1DWA16566HdvwfEOYsVrSbtQ+6AcDx1D71tnstY8HHgd2Sl7rpbrnVtT9voZ1Pxu/AKZK+o+ofcUS8IOIWGeSPUkHDBLrUF61T8m/86nNOnwwcOSgK0asAK4Grpb0OHAAcF1DLB+jlkzHDCwC/isipqWMz9rkXAKcS4XwEVcDSdsDfcBT1A6xlyeJ9l6G/lbxK2APSa+TNBr4UJOXuIlkOnhJewJPRsRz1A7Nd0mW7wJsW7fO1pLelfw+Dfh5ck55fEQsBI4Ddk6xe+N55Zz29GQ/0/gytffjrOTxddRmRd00iXcTSdsAvwT2VG0G2THAh5tsc39JYyW9ntrplVuT5fOS/SEi7m1cSdIukjZPfh8FvB14uKHNrtROAR0ar1zovwX4G0lvStqMk/TmlPtvLXAuDcq5lAMfcdW8Jjnkhdq3icMjYo2kC4ArJS3mlfP2rxIRj0r6OrUP22PAb6id5x7MKcC/J6dMXqB2mA1wKXBYEset1KYYGLAEOFy1qeR/B3yfWuJcLmlsEvPxKfbzLOBSSR+mdkrh+WHa1zsOmCvptIj4nKQvAv+ZfOBXAZ+OiFsknULtvP4fqF1sHyqhfwX8FNga+GpEPAYQEY9LWgL8ZIj1NqV2imb9uu00dhM+htp1i+uT6+KLI+ITyTfHC+vW/SLrvs/WPufS8JxLbfK0JjmRtGFE/CX5lvhjYG5E/DiH7U4CroqIv2p3W1WgWi+vu4FdImKo/7CshzmX8tHLueRThfk5JfmGdw+1HjxDfcuxIUh6H7Vv4t/ttUSzTJxLber1XPIRl5mZVYqPuMzMrFJcuMzMrFJcuMzMrFJcuMzMrFJcuMzMrFJcuMzMrFL+H2mWeUsPhL7nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "##### QC all libraries for how many reads per bargroup, decide which bargroups to keep based on number of reads they contain\n",
    "\n",
    "#vlines are the cutoffs for removing barcode groups that are too small (too few reads)\n",
    "#adjust vlines to match the inflection points of each library\n",
    "vlines = [50, 50, 80, 80, 80, 80]\n",
    "\n",
    "fig, axs = plot.subplots(1, len(libraries), figsize = [3*len(libraries), 5])\n",
    "fig.tight_layout()    \n",
    "\n",
    "#A jackpottocurve is a cumulative histogram of the barcode grouped reads binned by group size (# of reads per group)\n",
    "def jackpottocurve(fig, ax, barsize, vline):         \n",
    "        #prepare data to plot the jackpottocurve\n",
    "        X = np.array(barsize) #make an np array of the data so we can use the various np functions\n",
    "        X.sort()\n",
    "        X_lorenz = X.cumsum() / X.sum() #turns it into a cumsum percentage (for y-axis)\n",
    "        X_lorenz = np.insert(X_lorenz, 0, 0) #insert the 0,0 point\n",
    "        X = np.insert(X, 0, 0) #insert 0,0 point for X as well\n",
    "    \n",
    "        ## scatter plot of the jackpottocurve (modificaiton of the lorenz curve)\n",
    "        ax.scatter(X, X_lorenz, \n",
    "                   marker='o', color='darkgreen', s=100) #x axis is cumulative barcode number, yaxis is cumulative total # reads\n",
    "        ax.set_title(f\"lib{prefix[i]}\")\n",
    "        ax.set_xlim((1,max(barsize)))\n",
    "        ax.set_xscale(\"log\", basex=10)\n",
    "        ax.set_xlabel(\"Bargroups ranked by Size\")\n",
    "        ax.set_ylabel(\"Cumulative Reads Used\")\n",
    "        ax.vlines(vline, 0, 1)\n",
    "        return\n",
    "\n",
    "\n",
    "for i in range(len(libraries)):\n",
    "    #confirm how many bargroups there are\n",
    "    print(f\"{files[i]} the file contains {len(libraries[i])} bargroups\")\n",
    "    \n",
    "    #let's now plot the bargroup sizes\n",
    "    #first obtain the list of bargroup sizes\n",
    "    barsize = []\n",
    "    lib = libraries[i]\n",
    "    for b in lib:\n",
    "        barsize.append(len(lib[b].reads))\n",
    "    #plot them on the figure\n",
    "    ax = axs[i]\n",
    "    jackpottocurve(fig, ax, barsize, vlines[i])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "library 21-2-BF.sam\n",
      "130697 groups removed, 18825 groups kept, cutoff is 50\n",
      "library 21-8-BF.sam\n",
      "10838 groups removed, 740 groups kept, cutoff is 50\n"
     ]
    }
   ],
   "source": [
    "#now let's filter out the barcodes with specified cutoff # of reads (vlines)\n",
    "def filter_bar(groups, cutoff):\n",
    "    filtered_groups = {}\n",
    "    deleted = 0\n",
    "    for key in groups:\n",
    "        if len(groups[key].reads) >= cutoff:\n",
    "            filtered_groups[key] = groups[key]\n",
    "        else:\n",
    "            deleted += 1\n",
    "    print(f\"{deleted} groups removed, {len(filtered_groups)} groups kept, cutoff is {cutoff}\")\n",
    "    return filtered_groups\n",
    "\n",
    "filtered_libraries = []\n",
    "for i in range(len(libraries)):\n",
    "    print(f\"library {files[i]}\")\n",
    "    filtered_libraries.append(filter_bar(libraries[i], vlines[i]))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>PSA_OFF</th>\n",
       "      <th>PSA_ON</th>\n",
       "      <th>PSB_OFF</th>\n",
       "      <th>PSB_ON</th>\n",
       "      <th>PSD_OFF</th>\n",
       "      <th>PSD_ON</th>\n",
       "      <th>PSE_OFF</th>\n",
       "      <th>PSE_ON</th>\n",
       "      <th>PSF_OFF</th>\n",
       "      <th>PSF_ON</th>\n",
       "      <th>PSG_OFF</th>\n",
       "      <th>PSG_ON</th>\n",
       "      <th>PSH_OFF</th>\n",
       "      <th>PSH_ON</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>TGAAAAATGTTAGTG</th>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>44.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>32.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>58.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TCACATCAAGTAAAT</th>\n",
       "      <td>0.0</td>\n",
       "      <td>21.0</td>\n",
       "      <td>19.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TGATCTGCGTTAGAA</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>21.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ATATTGGGCCAAGGG</th>\n",
       "      <td>0.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>39.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>23.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AATACTAATTCGCAT</th>\n",
       "      <td>0.0</td>\n",
       "      <td>23.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>22.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>17.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>42.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ATTGAAGCGACCTAA</th>\n",
       "      <td>0.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>24.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>14.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TGGTTCGAACATCCT</th>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>15.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AATTGAACTATTTTT</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>15.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AACTGAAAATATACT</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>15.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>CAAGAGAACTGCGAC</th>\n",
       "      <td>0.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>18.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>18825 rows × 14 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                 PSA_OFF  PSA_ON  PSB_OFF  PSB_ON  PSD_OFF  PSD_ON  PSE_OFF  \\\n",
       "TGAAAAATGTTAGTG      0.0     2.0      9.0     0.0      0.0    44.0      0.0   \n",
       "TCACATCAAGTAAAT      0.0    21.0     19.0     0.0      0.0    17.0      0.0   \n",
       "TGATCTGCGTTAGAA      0.0     0.0     14.0     0.0      0.0    21.0      0.0   \n",
       "ATATTGGGCCAAGGG      0.0    18.0     24.0     0.0      0.0    39.0      0.0   \n",
       "AATACTAATTCGCAT      0.0    23.0     27.0     0.0      0.0    24.0      0.0   \n",
       "...                  ...     ...      ...     ...      ...     ...      ...   \n",
       "ATTGAAGCGACCTAA      0.0     7.0     11.0     0.0      0.0    24.0      0.0   \n",
       "TGGTTCGAACATCCT      0.0     4.0     12.0     0.0      0.0    15.0      0.0   \n",
       "AATTGAACTATTTTT      0.0     0.0     15.0     0.0      0.0    12.0      0.0   \n",
       "AACTGAAAATATACT      0.0     0.0      2.0     0.0      0.0    14.0      0.0   \n",
       "CAAGAGAACTGCGAC      0.0     6.0      5.0     0.0      0.0     6.0      0.0   \n",
       "\n",
       "                 PSE_ON  PSF_OFF  PSF_ON  PSG_OFF  PSG_ON  PSH_OFF  PSH_ON  \n",
       "TGAAAAATGTTAGTG    17.0      3.0     0.0     32.0     0.0      0.0    58.0  \n",
       "TCACATCAAGTAAAT     9.0     17.0     0.0     14.0     0.0     24.0     0.0  \n",
       "TGATCTGCGTTAGAA    11.0     17.0     0.0     12.0     0.0     16.0     0.0  \n",
       "ATATTGGGCCAAGGG    23.0     18.0     0.0     25.0     0.0      8.0     0.0  \n",
       "AATACTAATTCGCAT     8.0     22.0     0.0     17.0     0.0      0.0    42.0  \n",
       "...                 ...      ...     ...      ...     ...      ...     ...  \n",
       "ATTGAAGCGACCTAA     5.0      0.0     0.0     16.0     0.0      0.0    14.0  \n",
       "TGGTTCGAACATCCT     3.0      7.0     0.0      4.0     0.0      0.0    11.0  \n",
       "AATTGAACTATTTTT     6.0      0.0     0.0      9.0     0.0      0.0    11.0  \n",
       "AACTGAAAATATACT    14.0      5.0     0.0     15.0     0.0     14.0     0.0  \n",
       "CAAGAGAACTGCGAC     2.0     18.0     0.0      5.0     0.0      0.0    10.0  \n",
       "\n",
       "[18825 rows x 14 columns]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#turn the grouped reads into a pandas dataframe\n",
    "\n",
    "#datastructure as such: rows = each barcode, columns = # reads for each target\n",
    "dfs = []\n",
    "for i in range(len(filtered_libraries)):\n",
    "    groups = filtered_libraries[i]\n",
    "    rowlist = []\n",
    "    keylist = []\n",
    "    mappednames = set()\n",
    "    for b in groups:\n",
    "        keylist.append(b)\n",
    "        row = {}\n",
    "        for r in groups[b].reads:\n",
    "            if r.align in row:\n",
    "                row[r.align] += 1\n",
    "            else:\n",
    "                row[r.align] = 1\n",
    "        rowlist.append(row)\n",
    "        \n",
    "    df = pd.DataFrame(rowlist, index = keylist)\n",
    "    df = df.fillna(0) #replace NaN values with 0s\n",
    "    #fix column headers\n",
    "    #sort column header\n",
    "    df = df.reindex(sorted(df.columns), axis=1) #reindex based on sorted new columns\n",
    "    #fix names\n",
    "    cols = df.columns.tolist()\n",
    "    cols = [c if c.find(\"|\") == -1 else c.split(\"|\")[5] for c in cols]\n",
    "    df.columns = cols    \n",
    "    #add to the dfs list\n",
    "    dfs.append(df)    \n",
    "\n",
    "libraries = None #clear the libraries variable to free up memory\n",
    "    \n",
    "dfs[0]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export dataframe to excel sheet\n",
    "for i in range(len(prefix)):\n",
    "    with pd.ExcelWriter(f'{prefix[i]}.xls') as writer:  \n",
    "        dfs[i].to_excel(writer)"
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}