n_to_d_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import statments\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.pyplot import figure\n",
    "import random\n",
    "\n",
    "# load data\n",
    "norm_df = pd.read_json('./data/segments_df.json') # use segments_df, not pretty_segments, which has had putative Ns changed to Ds\n",
    "norm_df = norm_df[norm_df.pretty]\n",
    "deam_df = pd.read_json('./data/n_to_d_segments_df.json') # this is already filtered for \"pretty\" segments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate and print the percent putative deamidation of Asn in different conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normal conditions, % N-> D:10\n",
      " Basic conditions, % N-> D:96\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/rb/4nfyc06110b944pkssbl1yrw0000gn/T/ipykernel_628/263080559.py:2: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  deam_nn = len(deam_df[deam_df.aa == 'N'][deam_df['max'] < 1.3])\n",
      "/var/folders/rb/4nfyc06110b944pkssbl1yrw0000gn/T/ipykernel_628/263080559.py:3: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  deam_nd = len(deam_df[deam_df.aa == 'N'][deam_df['max'] >= 1.3])\n",
      "/var/folders/rb/4nfyc06110b944pkssbl1yrw0000gn/T/ipykernel_628/263080559.py:7: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  norm_nn = len(norm_df[norm_df.aa == 'N'][norm_df['max'] < 1.3])\n",
      "/var/folders/rb/4nfyc06110b944pkssbl1yrw0000gn/T/ipykernel_628/263080559.py:8: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  norm_nd = len(norm_df[norm_df.aa == 'N'][norm_df['max'] >= 1.3])\n"
     ]
    }
   ],
   "source": [
    "deam_n = len(deam_df[deam_df.aa == 'N'])\n",
    "deam_nn = len(deam_df[deam_df.aa == 'N'][deam_df['max'] < 1.3])\n",
    "deam_nd = len(deam_df[deam_df.aa == 'N'][deam_df['max'] >= 1.3])\n",
    "\n",
    "\n",
    "norm_n = len(norm_df[norm_df.aa == 'N'])\n",
    "norm_nn = len(norm_df[norm_df.aa == 'N'][norm_df['max'] < 1.3])\n",
    "norm_nd = len(norm_df[norm_df.aa == 'N'][norm_df['max'] >= 1.3])\n",
    "\n",
    "print(f\"Normal conditions, % N-> D:{round(norm_nd/(norm_n)*100)}\\n Basic conditions, % N-> D:{round(deam_nd/deam_n*100)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the percent putative deamidation of Asn in different conditions for each run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/rb/4nfyc06110b944pkssbl1yrw0000gn/T/ipykernel_628/2814218114.py:7: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  deam_n = len(deam_df[deam_df.run == run][deam_df.aa == 'N'])\n",
      "/var/folders/rb/4nfyc06110b944pkssbl1yrw0000gn/T/ipykernel_628/2814218114.py:8: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  deam_nn = len(deam_df[deam_df.run == run][deam_df.aa == 'N'][deam_df['max'] < 1.3])\n",
      "/var/folders/rb/4nfyc06110b944pkssbl1yrw0000gn/T/ipykernel_628/2814218114.py:9: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  deam_nd = len(deam_df[deam_df.run == run][deam_df.aa == 'N'][deam_df['max'] >= 1.3])\n",
      "/var/folders/rb/4nfyc06110b944pkssbl1yrw0000gn/T/ipykernel_628/2814218114.py:15: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  norm_n = len(norm_df[norm_df.run == run][norm_df.aa == 'N'])\n",
      "/var/folders/rb/4nfyc06110b944pkssbl1yrw0000gn/T/ipykernel_628/2814218114.py:16: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  norm_nn = len(norm_df[norm_df.run == run][norm_df.aa == 'N'][norm_df['max'] < 1.3])\n",
      "/var/folders/rb/4nfyc06110b944pkssbl1yrw0000gn/T/ipykernel_628/2814218114.py:17: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "  norm_nd = len(norm_df[norm_df.run == run][norm_df.aa == 'N'][norm_df['max'] >= 1.3])\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAFyCAYAAADvSk1ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4dUlEQVR4nO3deViU9f4//ufNjqyyKKtogKIppqRyNBRIPeZHxRBPaR3Xcj2fr5hJWobivp60rjKto9iHoFTQitRKhcQCcy200FIRQtzQBEYWGe7fH/6Y4zg3MCOz83xc11xXc9/vuefF3YzPed/L+y2IoiiCiIjoERaGLoCIiIwTA4KIiCQxIIiISBIDgoiIJDEgiIhIEgOCiIgkMSCIiEgSA4KIiCQxIJogiiLKy8vBewmJqDViQDShoqICLi4uqKioMHQpRER6x4AgIiJJDAgiIpLEgCAiIkkMCCIiksSAICIiSQwIIiKSxIAgIiJJDAgiIpJk8ICoqKhAQkIChg4dCk9PTwiCgCVLlki2PXXqFAYPHgxHR0e4uroiNjYWly5dkmz73nvvISQkBLa2tujUqROSkpJw//59Hf4lRETmxeABUVZWhq1bt6KmpgajR49utF1BQQEiIyNRW1uLnTt3Ytu2bbhw4QIiIiJw8+ZNpbYrVqzAnDlzEBsbi2+++QazZs3CypUrMXv2bB3/NaZHLpcjOzsbaWlpyM7Ohlwu12g9ET1glt8V0cDq6+vF+vp6URRF8ebNmyIAcfHixSrtxo4dK3p4eIh3795VLCssLBStra3FhIQExbJbt26JdnZ24rRp05Rev2LFClEQBPHcuXNq13b37l0RgNJ7mpP09HTRz89PBKB4+Pn5ienp6WqtJ6IHzPW7YvCAeFhjAXH//n3R3t5enD59usprhg4dKgYHByuep6SkiADE3NxcpXZXr14VAYgrVqxQux5zDoj09HRREASlDzQAURAEURAEcf78+U2uN/UPPpG2NPddMuXvisEPManj4sWLqKqqQmhoqMq60NBQ/PHHH6iurgYAnD17FgDQo0cPpXbe3t7w8PBQrG/N5HI55syZIzlKbcOyf//7302uj4+PN48uNFELqPNdMuXvikkERFlZGQDAzc1NZZ2bmxtEUcSdO3cUbW1tbeHg4CDZtmFbUmpqalBeXq70MEc5OTn4888/G10vimKTH2hRFFFcXIycnBxdlEdkMtT5Lpnyd8UkAqKBIAhqrVO33aNWrVoFFxcXxcPf3//xCjVypaWlRrUdIlOl7nfAVL8rJhEQ7u7uACD56//27dsQBAGurq6KttXV1bh3755kW6leSIOFCxfi7t27ikdxcbF2/gAj4+3tbVTbITJV6n4HTPW7YhIBERgYCHt7e+Tn56usy8/PR1BQEOzs7AD899zDo22vXbuGW7duoXv37o2+j62tLZydnZUepkgmk0EQBAiCAJlMprI+IiICfn5+jfamBEGApaVlo9sXBAH+/v6IiIjQWs1Epkid75Ipf1dMIiCsrKwwcuRIZGRkKM3uVlRUhKysLMTGxiqWDRs2DHZ2dkhOTlbaRnJyMgRBaPJei9bC0tISmzZtAqB6yK3h+WuvvaYIGan1GzdubDJEiMzVwz/Aqqurm/0umfR3xVCXTz1s37594q5du8Rt27aJAMSxY8eKu3btEnft2iXKZDJRFEXxt99+Ex0dHcWBAweK+/btEzMyMsTu3buLPj4+4o0bN5S2t3z5clEQBPHNN98Us7OzxXXr1om2trbiq6++qlFdpnqZa2VlpeJSu8rKykbbSV277e/v3+R9EA+vJ2qNpL5f5vpdEURR4vosPevYsSOuXLkiue7y5cvo2LEjAODkyZN44403kJubCysrK0RHR2P9+vUIDAxUed27776L999/H4WFhfDy8sLkyZPx1ltvwdraWu26ysvL4eLigrt375rU4SaZTAZHR0cAQGVlpeQVXQ3kcjlycnJQWloKb29vREREKP3aaW49UWvT2PfLHL8rRhEQxqo1BAQRaaY1fb9M4hwEERHpHwOCiIgkMSCIiEgSA4KIiCQxIIiISJKVoQsgItLUq+EJGr+mXqzHtYo/UVUrg72NA7yc/GAhaP4b+b68VvHfsyPfgrWljcbb0JaP8tbqdPsMCBMnde01ESm7fPsC8q4cgqy2UrHMwcYR4QHPopNbZwNWZtwYECYsIyMDc+bMURpu2M/PD2vWrDFgVUTG5fLtCzj0+xcqy2W1lTj0+xd4Njim0ZCQ6nW0JgwIE5WRkYG4uDiViUpKSkrw0ksvGagqIuNSL9Yj78qhJtvkXTmMgLZBKoebGut19PEfpJNajRFPUpsgdWaxIiLgWsWfSv/AS5HVVuBahfKkPw29jkdfK6utRPbFr7Vep7FiQJig5maxIqIHqmpVh7tvrp06vY7WggFhgkx1dioifbO3UW+cpIfbqdPraC0YECbIVGenItI3Lyc/ONg4NtnGwcZJ6eSzur2O1oABYYKam8WKiB6wECwQHvBsk23CA6KVTlCr2+toDRgQJkidGeGI6IFObp3xbHCMSk/CwcZJ8hJXdXsdU/rOM+hNcvrAgDBRsbGx2L17N3x9fZWW+/n54dNPPzVQVUTGqZNbZ7zw1HQM7/oCogJHYHjXF/DCU9Mk7394nF6HueJ9ECYsNjYWMTExKndSV1dXG7o0IqNjIVjAx7mDWm0beh2q90E4ITwgutXcfc2AMHGWlpaIjIw0dBlEZqeTW2cEtA3SyvhNpooBQUTUCE16Heao9UQhERFphD0IEyU1iqulpaWhyyIiM8KAMEGNjeK6adMmxMbGGrAyIjInDAgdefrAIp1s984PBbi0PF1l+Z9//okxY8bgiUVj4Bz2hGL5M98lwdLOsNdqnxi23KDvT0SPh+cgTIgor0fxh9822aZ4y3cQ5fV6qoiIzBkDwoRUnivG/VsVTba5f7Mcst9K9FQREZkzBoQJuX9bvREm1W1HRNQUBoQJsXZrenwYTdsRETWFAWFCHJ/0h7WHU5NtrD2d4dDVt8k2RETqYECYEMHSAv4zhjbZxn/6EAiW/N9KRC3Hf0lMTNsBIXhi0RiVnoS1pzOeWDQGbQeEGKgyIjI3vA/CBLUdEALX8M4Prmq6XQlrN0c4PunPngMRaRUDwkQJlhZwCg0wdBlEZMb4k5OIiCQxIIiISBIDgoiIJDEgiIhIEgOCiIgkMSCIiEgSA4KIiCTxPggjIMrredMbERkdBoSB3fmhAMUffqs0z4O1hxP8ZwzlsBlEZFD8mWpADdOHPjoJ0P1bFbi0PB13figwUGVERAwIg+H0oURk7BgQBqLu9KGV54r1VBERkTIGhIFw+lAiMnYmFRCnT5/G6NGj4ePjgzZt2iAkJARLly7FvXv3lNqdOnUKgwcPhqOjI1xdXREbG4tLly4ZqGppnD6UiIydyQTEr7/+iv79+6OwsBAbN25EZmYmXnzxRSxduhTjxo1TtCsoKEBkZCRqa2uxc+dObNu2DRcuXEBERARu3rxpwL9AmbrThzo+6a+nioiIlJnMZa6pqamorq5Geno6AgMDAQDR0dEoLS3F1q1bcefOHbRt2xaJiYmwtbVFZmYmnJ2dAQBhYWEIDg7G+vXrsWbNGkP+GQoN04deWp7eaBtOH0pEhmQy//pYW1sDAFxcXJSWu7q6wsLCAjY2Nqirq0NmZibGjBmjCAcACAgIQFRUFPbs2aPXmpujq+lDLe1sELb/LYTtfwuWdjbaKJWIWiGT6UFMnDgRGzduxMyZM7FmzRp4enri+++/x5YtWzB79mw4ODjg/PnzqKqqQmhoqMrrQ0ND8d1336G6uhp2dnYG+AukcfpQIjJWJhMQHTt2RG5uLp5//nnFISYA+H//7/9h48aNAICysjIAgJubm8rr3dzcIIoi7ty5A29vb8n3qKmpQU1NjeJ5eXm5Fv+CxnH6UCIyRiYTEIWFhRg5ciTat2+P3bt3w9PTE8eOHcPy5ctRWVmJ//znP4q2giA0up2m1q1atQpJSUlarZuIyFS1KCBu3ryJqqoqleUdOnRoyWYlLViwAOXl5Thz5gwcHBwAAAMHDoSHhwemTJmCCRMmwMvLC8B/exIPu337NgRBgKura6PvsXDhQrz22muK5+Xl5fD351VERNQ6aRwQFRUVmDt3LtLS0lBdXS3ZRi6Xt7iwR505cwbdunVThEODPn36AADOnj2LAQMGwN7eHvn5+Sqvz8/PR1BQUJPnH2xtbWFra6vdwomITJTGAREfH4/U1FRMnToVoaGhevsH1cfHB2fPnkVlZSUcHf9781hubi4AwM/PD1ZWVhg5ciQyMjKwdu1aODk9uDqoqKgIWVlZmDt3rl5qJSIyBxoHxNdff43Vq1djzpw5uqinUfHx8Rg9ejSGDBmCuXPnwsPDA3l5eVi1ahW6deuG5557DgCQlJSEPn36YMSIEViwYAGqq6uRmJgIDw8PzJs3T681ExGZMo2vpayurkaPHj10UUuTRo0ahUOHDsHZ2Rlz5szBiBEjsGPHDkyfPh1HjhyBjc2D6/1DQkKQnZ0Na2trxMXFYdKkSQgKCsKRI0fg6emp97qJiEyVxj2I4cOHIycnB9HR0bqop0lRUVGIiopqtl1YWBgOHjyoh4qIiMyXxgGxaNEixMXFwcnJCSNHjoS7u7tKG6n7EIiIyLRoHBDdu3cHAMyfPx/z58+XbKOLq5iIiEi/NA6IxMTEJm82IyIi86BxQCxZskQHZZC2iPJ6jutERFrRojupq6urFcNsG9MAeK3VnR8KUPzht0pTmVp7OMF/xtDHHhmWiFqvx/pp+eOPPyIiIgJOTk7w8/ODk5MTBg0apLhpjfTvzg8FuLQ8XWWe6/u3KnBpeTru/FBgoMqIyFRp3IPIy8tDdHQ0XF1dMW3aNPj4+KCkpAQZGRmIjo5GdnY2+vXrp4taqRGivB7FH37bZJviLd/BNbwzDzcRkdoe6yR1aGgosrKylMZFWrduHaKiopCYmIhvvvlGq0VS0yrPFav0HB51/2Y5Ks8Vc1hxIlKbxj8n8/LykJCQoDJonoODA+bPn8/DTAZw/3alVtsREQGPERByubzRAfrs7Ox4D4QBWLs5Nt9Ig3ZERMBjBETPnj2xefNmyXVbtmxBz549W1wUacbxSX+Vea0fZe3pDMcnObcFEalP43MQCxYswOjRo9GrVy+8/PLL8Pb2RmlpKVJTU3HmzBns3btXB2VSUwRLC/jPGIpLy9MbbeM/fQhPUBORRjQOiFGjRiElJQUJCQlKQ234+voiJSUFI0eO1GqBpJ62A0LwxKIxqvdBeDrDf/oQ3gdBRBp7rBvlxo8fj3HjxuH8+fMoKyuDu7s7unTpwiE4DKztgBC4hnfmndREpBWPfSe1IAgICeGvUmMjWFrwUlYi0gq1AuLIkSPo3bs3HB0dceTIkWbbDxw4sMWFERGRYakVEJGRkcjLy0Pfvn0RGRnZ6KEkURQhCAIvdSUiMgNqBURWVha6desGADh8+DDPNRARtQJqBcSgQYMU/x0ZGamrWoiIyIhofHlLdHQ0CgqkRwa9cOGCQeaqJiIi7dM4ILKzs1FeXi65rqKiAt9//32LiyIiIsPT6gXypaWlaNOmjTY3SUREBqLWOYgvvvgCX3zxheL5smXL4OnpqdSmqqoK2dnZ6NWrl3YrJCIig1ArIH799Vfs2rULwIMb5A4fPgwLC+XOh62tLXr06IFNmzZpv0oiItI7tQJi4cKFWLhwIQDAwsICWVlZ6Nu3r04LIyIiw9J4qI36+npd1EFEREaGo7gREZEktXoQTzzxBPbs2YOePXuiU6dOTd5JLQgCLl68qLUCiYjIMNS+k9rZ2Vnx3xxqg4jI/KkVENu3b1f8d3Jysq5qISIiI8JzEEREJEnjgDh8+LDinggAuH79OoYPHw4vLy9MmDAB1dXVWi2QiIgMQ+OASExMxK+//qp4npCQgJycHPTv3x+7d+/GunXrtFogEREZhsYBceHCBfTu3RsAUFdXhz179mDNmjXIyMjA0qVLkZaWpvUiiYhI/zQOiPLycri6ugIATp48CZlMhlGjRgEA+vbti6KiIq0WSEREhqFxQLRr1w6///47AODgwYMICAiAn58fgAfDfVtbW2u3QiIiMgiNh9oYNmwY3nzzTZw7dw7JycmYOHGiYl1BQQE6duyozfqIiMhANA6IlStXoqioCB999BH69u2LRYsWKdalpqaif//+Wi2QiIgMQ+OA8PDwwIEDByTXZWVlwc7OrsVFERGR4bXoRrnq6mqUlpYq7n1wdnaGjY2NVgojIiLDeqyA+PHHHxEREQEnJyf4+fnByckJgwYNQm5urrbrIyIiA9H4EFNeXh6io6Ph6uqKadOmwcfHByUlJcjIyEB0dDSys7PRr18/XdRKRER6pHFAJCYmIjQ0FFlZWXBwcFAsX7duHaKiopCYmIhvvvlGq0USEZH+aXyIKS8vDwkJCUrhAAAODg6YP38+DzMREZkJjQNCLpfD1tZWcp2dnR3kcnmLiyIiIsPTOCB69uyJzZs3S67bsmULevbs2eKiiIjI8DQOiAULFuDAgQPo1asXNmzYgNTUVGzYsAFhYWHYv38/Fi5cqIs6FY4ePYrhw4ejbdu2sLe3R3BwMJYtW6bU5tSpUxg8eDAcHR3h6uqK2NhYXLp0Sad1ERGZG41PUo8aNQopKSlISEjA/PnzFct9fX2RkpKCkSNHarXAh6WmpuKf//wn/vGPf+CTTz6Bo6MjLl68iKtXryraFBQUIDIyEk899RR27tyJ6upqJCYmIiIiAmfOnIGnp6fO6iMiMicaBwQAjB8/HuPGjcP58+dRVlYGd3d3dOnSRadzVZeUlGDatGmYPn06PvjgA8XyqKgopXaJiYmwtbVFZmamYh7tsLAwBAcHY/369VizZo3OaiQiMiePfSe1IAgICQnBgAEDEBISotNwAICPP/4YMpkMb7zxRqNt6urqkJmZiTFjxijCAQACAgIQFRWFPXv26LRGIiJz8lgBUVhYiOnTp6Nz585wd3dH586dMX36dFy+fFnb9SkcOXIEbm5uKCgowFNPPQUrKyu0a9cOM2bMQHl5OQDg4sWLqKqqQmhoqMrrQ0ND8ccff3BKVCIiNWkcEGfOnEGvXr2QnJwMX19fDB06FL6+vkhOTkavXr1w5swZHZT54BDTvXv3MHbsWLzwwgs4ePAg5s+fj08++QTDhw+HKIooKysDALi5uam83s3NDaIo4s6dO42+R01NDcrLy5UeREStlcbnIOLj4+Hp6YmDBw+iQ4cOiuVXrlzBkCFDMHfuXGRlZWm1SACor69HdXU1Fi9ejAULFgAAIiMjYWNjg/j4eBw6dAht2rQBgCYPdzW1btWqVUhKStJu4UREJkrjHsRPP/2EpKQkpXAAHhznX7JkCY4dO6a14h7m7u4OAPj73/+utPy5554D8ODS1oY2DT2Jh92+fRuCICimS5WycOFC3L17V/EoLi7WUvVERKZH44BwcXGBi4uL5DpXV1elk8PaJHVeAQBEUQQAWFhYIDAwEPb29sjPz1dpl5+fj6CgoCbnq7C1tYWzs7PSg4iotdI4IMaPH4+PP/5Yct1HH32EcePGtbgoKWPGjAEA7N+/X2n5vn37AADh4eGwsrLCyJEjkZGRgYqKCkWboqIiZGVlITY2Vie1ERGZI43PQfTu3Ru7d+9G3759MW7cOHh5eeHatWtIS0vDjRs3MHbsWGRkZCjaa+sf5aFDh2LkyJFYunQp6uvrER4ejhMnTiApKQkjRozAM888AwBISkpCnz59MGLECCxYsEBxo5yHhwfmzZunlVqIiFoDQWw4RqMmC4umOx2CICgO+wiCoNXB+6qqqpCUlITU1FSUlpbCx8cHL730EhYvXqw0gODJkyfxxhtvIDc3F1ZWVoiOjsb69esRGBio0fuVl5fDxcUFd+/e1fhw09MHFjXfqJU4MWy5oUsgM/NqeIKhSzAKH+Wt1en2Ne5B6OIKJXXZ29tj9erVWL16dZPtwsLCcPDgQT1VRURknjQOiEGDBumiDiIiMjKPPdQGERGZNwYEERFJYkAQEZEkBgQREUliQBARkSQGBBERSVLrMteH74xWB4e0ICIyfWoFRFxcnNId0lIahtEWBAF1dXXaqY6IiAxGrYBo7u7p6upqvPvuu9i/fz8sLS21UhgRERmWWgHR1N3TO3bswOLFi1FUVISYmBisWrVKa8UREZHhPPZJ6n379iE0NBSTJ0+Gv78/jh49ij179iAkJESb9RERkYE81oxykZGRGDlyJOrr67F3717k5OSgf//+uqiPiIgMRO2A+P333zF27FiEh4fj4sWL2Lp1K3755ReMGjVKl/UREZGBqHUOYtasWfjPf/4DBwcHrF69GnPmzFGaf4GIiMyPWgHx4YcfQhAEdOvWDUePHsXRo0cbbSsIAr744gutFUhERIahVkB06NABgiCgpKQEJSUlTbZtuB+CiIhMm1oBUVhYqOMyiIjI2HAsJiIiksSAICIiSQwIIiKSxIAgIiJJDAgiIpLEgCAiIkkMCCIikqTWfRCPSklJQWpqKq5cuYKqqiqldYIg4OLFi1opjoiIDEfjgFizZg0WLlyIbt26oWfPnhyTiYjITGkcEFu3bsXs2bPx3nvv6aIeIiIyEhqfg7h27Rqef/55XdRCRERGROOACAsL4zkGIqJWQOOA+Pe//40NGzbg5MmTuqiHiIiMhMbnICZPnoyysjL07dsXXl5ecHd3V1ovCAJ+/vlnrRVIRESGoXFAuLu7w8PDQxe1EBGREdE4ILKzs3VQBhERGRveSU1ERJLU6kEUFRXB29sb1tbWKCoqarZ9hw4dWlwYEREZlloB0alTJ+Tm5qJv377o2LFjs/NOy+VyrRRHRESGo1ZAbNu2DYGBgYr/bi4giIjI9KkVEBMnTlT896RJk3RVCxERGRGepCYiIkkMCCIiksSAICIiSQwIIiKSxIAgIiJJDAgiIpKkcUAcPnwYu3btUjy/fv06hg8fDi8vL0yYMAHV1dVaLZCIiAxD44BITEzEr7/+qniekJCAnJwc9O/fH7t378a6deu0WmBTPv74YwiCAEdHR5V1p06dwuDBg+Ho6AhXV1fExsbi0qVLequNiMjUaRwQFy5cQO/evQEAdXV12LNnD9asWYOMjAwsXboUaWlpWi9SSklJCV5//XX4+PiorCsoKEBkZCRqa2uxc+dObNu2DRcuXEBERARu3rypl/qIiEydxgFRXl4OV1dXAMDJkychk8kwatQoAEDfvn3VGsxPG2bMmIGBAwdiyJAhKusSExNha2uLzMxMDB8+HLGxsfj6669x8+ZNrF+/Xi/1ERGZOo0Dol27dvj9998BAAcPHkRAQAD8/PwAABUVFbC2ttZuhRJSUlLw/fff44MPPlBZV1dXh8zMTIwZMwbOzs6K5QEBAYiKisKePXt0Xh8RkTnQeMKgYcOG4c0338S5c+eQnJysNE5TQUEBOnbsqM36VNy4cQPx8fFYvXq1IpgedvHiRVRVVSE0NFRlXWhoKL777jtUV1fDzs5Op3USEZk6jQNi5cqVKCoqwkcffYS+ffti0aJFinWpqano37+/Vgt81KxZs9ClSxfMnDlTcn1ZWRkAwM3NTWWdm5sbRFHEnTt34O3trbK+pqYGNTU1iufl5eVaqpqIyPRoHBAeHh44cOCA5LqsrCyd/jJPT0/HV199hdOnTzc75HhT6xtbt2rVKiQlJbWoRiIic6HVG+WcnZ3x008/aXOTCpWVlZg9ezb+93//Fz4+Pvjrr7/w119/oba2FgDw119/QSaTwd3dHcB/exIPu337NgRBUJxkf9TChQtx9+5dxaO4uFgnfwsRkSnQWkAcO3YMQ4cOxaBBg7S1SSW3bt3C9evXsWHDBrRt21bxSEtLg0wmQ9u2bfHSSy8hMDAQ9vb2yM/PV9lGfn4+goKCGu3l2NrawtnZWelBRNRaqR0Qn332GaKiotCtWzc8//zzOHPmDIAHJ4VHjRqF/v374+jRo3j99dd1UqiXlxeysrJUHn//+99hZ2eHrKwsLF++HFZWVhg5ciQyMjJQUVGheH1RURGysrIQGxurk/qIiMyNWucgPvvsM4wfPx4A4OnpiczMTGRlZSEtLQ0vvPACZDIZJkyYgGXLlkleWaQNdnZ2iIyMVFmenJwMS0tLpXVJSUno06cPRowYgQULFqC6uhqJiYnw8PDAvHnzdFIfEZG5UasH8d5776F79+4oLCzE9evXcevWLQwaNAjPP/88bG1t8f3332P79u06CwdNhYSEIDs7G9bW1oiLi8OkSZMQFBSEI0eOwNPT09DlERGZBLV6EGfPnsWWLVvQoUMHAICLiwvWr1+Pr776CqtWrcIzzzyj0yKbkpycjOTkZJXlYWFhOHjwoP4LIiIyE2r1ICoqKtCpUyelZQ3Pe/Toof2qiIjI4NQ+Sf3ovQMNz/UxtAYREemf2jfKbdiwAe3bt1c8F0URgiBg3bp1Ssf1BUHApk2btFslERHpndoB8fAkQQ/7/PPPlZ4zIIiIzINaAVFfX6/rOoiIyMiodQ6Ck+wQEbU+agWEr68v4uLisH//foiiqOuaiIjICKgVEGPHjsX+/fsxYsQI+Pv7Y9GiRbh48aKuayMiIgNSKyA+/fRTlJaW4v3334evry9WrlyJzp07IyoqCikpKaiurtZ1nUREpGdq3wfh7OyMGTNm4NixYzh37hzmzp2LgoICTJgwAV5eXpg5cyaOHz+uy1qJiEiPHmu4765du2L9+vX4888/sXfvXkRGRmLbtm0IDw+XnOqTiIhMT4vmg7C0tMSoUaOwZcsW/Otf/wIAnDt3TiuFERGRYWk85WgDuVyOL7/8Etu3b8eBAwdQV1eH0NBQTJ06VZv1ERGRgWgcEOfOncO2bduQkpKCW7duwdnZGa+88gqmTp2KsLAwXdRIREQGoFZAlJeXIzU1Fdu2bcPJkycBAAMHDsTUqVMRFxfX6BSeRERkutQKCC8vL9TU1MDb2xsLFizAlClTEBgYqOvaiIjIgNQKiGHDhmHq1Kl47rnnYGHRovPaRERkItQKiIyMDF3XQURERobdASIiksSAICIiSQwIIiKSxIAgIiJJDAgiIpLUooCora1FUlIS+vXrh759+2Lx4sWoqanRVm1ERGRAjz0WEwDMmzcPOTk5mDp1KioqKvDee+/hxo0b2Lx5s7bqIyIiA1ErIMrKyuDu7q6yPD09HadPn0b79u0BAF26dMGMGTMYEEREZkCtQ0xdunTBJ598orLc2toaMplM8Vwmk8HGxkZ71RERkcGo1YPYvHkz4uPjsWPHDmzdulUxDtOkSZMQHR2NsWPHQiaT4dNPP8Xs2bN1WjAREemHWj2IsWPH4rfffkNwcDB69uyJlStXoq6uDkuWLMHbb7+NwsJC3LhxA+vXr8fy5ct1XTMREemB2iepnZ2d8eGHH2LChAmYPn060tLSsHXrVkydOpWTBBERmSGNL3Pt378/Tp8+jRdeeAGDBw/GrFmzUFFRoYvaiIjIgNQOiMuXL2Pr1q1499138fPPP2PRokU4c+YMzp8/j5CQEKSnp+uyTiIi0jO1AuLLL79E165dsWHDBuzYsQPh4eFYunQpgoODcejQIaxYsQIzZszAqFGj8Oeff+q6ZiK9kMvlyM7ORlpaGrKzsyGXyw1dEpFeqRUQCxYswLx583D+/HmcPHkSu3btwrJlyxSHliZNmoSCggK4urqiW7duOi2YSB8yMjLQsWNHREVFYfz48YiKikLHjh05Nwq1KmoFRElJCaKiohTPIyMjIZfLce3aNcUyd3d3fPLJJ9i7d6/WiyTSp4yMDMTFxan0hktKShAXF8eQoFZDrYAYMGAAkpKS8NNPP+Hs2bOYO3cuvL298cQTT6i0jY6O1nqRRPoil8sxZ84ciKKosq5hWXx8PA83UaugVkBs3boV1tbW+Nvf/obQ0FAcO3YMu3fvhqWlpa7rI9KrnJycJs+jiaKI4uJi5OTk6LEqIsNQ6z4IPz8/HD58GFVVVaitrYWLi4uu6yIyiNLSUq22IzJlGo3mam9vD3t7e13VQmRw3t7eWm1HZMo4YRDRQyIiIuDn5wdBECTXC4IAf39/RERE6LkyIv1jQBA9xNLSEps2bQIAlZBoeL5x40aef6NWgQFB9IjY2Fjs3r0bvr6+Ssv9/Pywe/duxMbGGqgyIv1q0YxyROYqNjYWMTExyMnJQWlpKby9vREREcGeA7UqDAiiRlhaWiIyMtLQZRAZDA8xERGRJJMJiMOHD2PKlCkICQmBg4MDfH19ERMTg5MnT6q0PXXqFAYPHgxHR0e4uroiNjYWly5dMkDVRESmy2QCYvPmzSgsLMScOXOwb98+bNq0CTdu3EB4eDgOHz6saFdQUIDIyEjU1tZi586d2LZtGy5cuICIiAjcvHnTgH8BEZFpMZlzEO+//z7atWuntGzYsGEICgrCypUrFWNAJSYmwtbWFpmZmXB2dgYAhIWFITg4GOvXr8eaNWv0XjsRkSkymR7Eo+EAAI6OjujWrRuKi4sBAHV1dcjMzMSYMWMU4QAAAQEBiIqKwp49e/RWLxGRqTOZgJBy9+5dnDp1Ck8++SQA4OLFi6iqqkJoaKhK29DQUPzxxx+orq7Wd5lERCbJZA4xSZk9ezZkMhneeustAEBZWRkAwM3NTaWtm5sbRFHEnTt3Gh1Hp6amBjU1NYrn5eXlOqiaiMg0mGwP4u2338ann36Kd955B2FhYUrrGhtHp7l1q1atgouLi+Lh7++vtXqJiEyNSQZEUlISli9fjhUrVuBf//qXYrm7uzuA//YkHnb79m0IggBXV9dGt7tw4ULcvXtX8Wg4t0FE1BqZ3CGmpKQkLFmyBEuWLMGbb76ptC4wMBD29vbIz89XeV1+fj6CgoJgZ2fX6LZtbW1ha2ur9ZqJiEyRSfUgli1bhiVLlmDRokVYvHixynorKyuMHDkSGRkZqKioUCwvKipCVlYWB1kzEXK5HNnZ2UhLS0N2dnaLpveUyWQQBAGCIEAmk2mxSiLzZzI9iA0bNiAxMRHDhg3D//zP/yAvL09pfXh4OIAHPYw+ffpgxIgRWLBgAaqrq5GYmAgPDw/MmzfPEKWTBjIyMjBnzhylaT/9/PywadMmBjyRnplMQHz11VcAgAMHDuDAgQMq6xsmlA8JCUF2djbeeOMNxMXFwcrKCtHR0Vi/fj08PT31WjNpJiMjA3FxcYr/lw1KSkoQFxfHobaJ9MxkAiI7O1vttmFhYTh48KDuiiGtk8vlmDNnjko4AA/CXxAExMfHIyYmhkNuE+mJSZ2DIPOVk5OjdFjpUaIoori4GDk5OXqsiqh1Y0CQUSgtLdVqOyJqOQYEGYXG7m5/3HZE1HIMCDIKERER8PPza/ROd0EQ4O/vj4iICD1XRtR6MSDIKFhaWmLTpk0AVIdDaXi+ceNGnqAm0iMGBBmN2NhY7N69G76+vkrL/fz8eIkrkQGYzGWu1DrExsYiJiYGOTk5KC0thbe3NyIiIthzIDIABgQZHUtLS0RGRhq6DKJWj4eYiIhIEgOCiIgkMSCIiEgSA4KIiCQxIIiISBIDgoiIJPEyVzIZcrmc90cQ6REDgkwCZ5oj0j8eYiKj1zDT3KPzRTTMNJeRkWGgyojMGwOCjFpzM80BQHx8PORyub5LIzJ7DAgyapxpjshwGBBk1DjTHJHhMCDIqHGmOSLDYUCQUeNMc0SGw4Ago8aZ5ogMh/dBkM69Gp7Q4m1EB41C3pVDkNVWKpa1sXZEeEA09q/Nw/61eZKvuy+vVfz37Mi3YG1p0+JaHtdHeWsN9t5Ej4MBQSahk1tnBLQNwrWKP1FVK4O9jQO8nPxgIbATTKQrDAgyGRaCBXycOxi6DKJWgz+/iIhIEnsQRNQoDpDYujEgiEgSB0gkBgQZpXqxniekDahhgMRHx8BqGCBx9+7dLQoJ9kxMAwOCjM7l2xdULml1sHFEeMCz6OTW2YCVtQ7NDZAoCALi4+MRExPzWP+os2diOviTjIzK5dsXcOj3L5TCAQBktZU49PsXuHz7goEqaz10OUAih243LQwIMhr1Yj3yrhxqsk3elcOoF+v1VFHrpKsBEjl0u+lhQJDRuFbxp0rP4VGy2gpcq2j81y21nK4GSOTQ7aaHAUFGo6pWptV29HiaGyARAHx9fTUeIJFDt5seBgQZDXsbB622o8fT1ACJDdauXavxCWoO3W56GBBkNLyc/OBg49hkGwcbJ3g5+empotYrNjYWu3fvhq+vr+T6mJgYjbfJodtNDwOCjIaFYIHwgGebbBMeEM37IfQkNjYWhYWFyMrKQmpqKvbt29ei7XHodtPDbxoZlU5unfFscIxKT8LBxgnPBsfwPgg9s7S0RGRkJMaNG4eBAwe2eHuN9Uz8/PxafPMdaR9vlCOjw6G9pT19YJFB319e/d+5NZ75LgmWdo85t0YboN2H/0Sbc8W4f7sS1m6OcHzSHystT2HlgVNqbaLX470zaYgBQUaJQ3ubN8HSAk6hAYYug5rRun+SERFRoxgQREQkiQFBRESSGBBERCSJAUFERJLMNiAqKysRHx8PHx8f2NnZ4amnnsJnn31m6LKIiEyG2V7mGhsbi+PHj2P16tXo3LkzUlNTMW7cONTX12P8+PGGLo+0RF8zz5nyDHeivB6Vj9xzIFiaRu1kWGYZEPv27cN3332nCAUAiIqKwpUrVzB//ny88MILvJ3fDOhr5jlTnuHuzg8FKP7wW9y/VaFYZu3hBP8ZQ9F2QIhG27K0s0HY/re0XSIZMbP8GbFnzx44Ojpi7NixSssnT56Mq1ev4tixYwaqjLRF3ZnnrC1t8Eq/+Xil33xYW2p+568pz3B354cCXFqerhQOAHD/VgUuLU/HnR8KDFQZmQqzDIizZ8+ia9eusLJS7iCFhoYq1pPp0tfMc6Y8w50or0fxh9822aZ4y3cQ5cZXOxkPszzEVFZWhieeeEJluZubm2K9lJqaGtTU1Cie3717FwBQXl6ucQ1yWU3zjVqJ2rra5htpoFTNmeeK/7oM7xYMDa7t93mcz9HDNPlMVZwrUuk5POr+zXLcPXkRTk+a3pAm2v5MmaqWfKacnJyanBQKACCaoeDgYHHYsGEqy69evSoCEFetWiX5usWLF4sA+OCDDz7M/nH37t1m/y01yx6Eu7u7ZC/h9u3bAP7bk3jUwoUL8dprryme19fX4/bt23B3d28+aY1MeXk5/P39UVxcDGdnZ0OXY1DcF9rDffmAOewHJyenZtuYZUD06NEDaWlpqKurUzoPkZ+fDwDo3r275OtsbW1ha2urtMzV1VVndeqDs7OzyX6AtY37Qnu4Lx8w9/1gliepn3/+eVRWViI9PV1p+Y4dO+Dj44N+/foZqDIiItNhlj2I5557DkOGDMHMmTNRXl6OoKAgpKWl4cCBA0hJSeE9EEREajDLgACAjIwMvPXWW0hMTMTt27cREhKCtLQ0vPjii4YuTS9sbW2xePFilUNmrRH3hfZwXz7QWvaDIIqiaOgiiIjI+JjlOQgiImo5BgQREUliQBARkSQGhBlJTk6GIAiws7PDlStXVNZHRkY2eg9IS96v4WFlZQU/Pz9MnjwZJSUlGm/v3r17WLJkCbKzsx+7pqtXr2LJkiU4c+aMyrolS5YY5IbHR/eTnZ0dvLy8EBUVhVWrVuHGjRt6r0ld/Ey17s8UA8IM1dTUYNGiRXp7v+3btyM3NxffffcdXn31VaSlpSEiIgIymUyj7dy7dw9JSUkt/jInJSVJfplfeeUV5ObmPva2W+rh/fT+++/jqaeewpo1a9C1a1ccPHjQYHWpg5+p1vmZYkCYoWHDhiE1NRU///yzXt6ve/fuCA8PR1RUFBYvXoyEhARcvnwZe/fu1cv7q8vPzw/h4eEGe/+G/RQREYExY8bgnXfewS+//AIHBwfExsbi+vXrBqutOfxMSTP3zxQDwgwlJCTA3d0db7zxhkHev+EL03BIIjIyEpGRkSrtJk2ahI4dOwIACgsL4enpCQBISkpSdJ0nTZoEAPjjjz8wefJkBAcHo02bNvD19cXIkSMVw6cAQHZ2Nvr06QPgwdwfDdtYsmQJANXDAaNHj0ZAQADq61WHvO7Xrx969+6teC6KIj744AM89dRTsLe3R9u2bREXF4dLly493k76/3Xo0AEbNmxARUUFtmzZ0qJt6RI/U63zM8WAMENOTk5YtGgRvvnmGxw+fFjv7//HH38AgOLLqQ5vb28cOHAAADB16lTk5uYiNzcXb7/9NoAH3Xx3d3esXr0aBw4cwPvvvw8rKyv069cP58+fBwD07t0b27dvBwAsWrRIsY1XXnlF8j2nTJmCoqIilX1UUFCAn376CZMnT1Ysmz59OuLj4zF48GDs3bsXH3zwAc6dO4f+/fu3+Ffa8OHDYWlpiSNHjrRoO7rEz1Qr/Uy1aFxtMirbt28XAYjHjx8Xa2pqxCeeeEJ8+umnxfr6elEURXHQoEHik08+qfX3y8vLE+/fvy9WVFSImZmZoqenp+jk5CReu3ZN8b6DBg1Sef3EiRPFgIAAxfObN2+KAMTFixc3+951dXVibW2tGBwcLM6dO1ex/Pjx4yIAcfv27SqvaRjOvcH9+/fF9u3bi+PHj1dql5CQINrY2Ii3bt0SRVEUc3NzRQDihg0blNoVFxeL9vb2YkJCQpO1Pvz/pTHt27cXu3bt2uR2DIGfqdb9mWIPwkzZ2Nhg+fLlOHHiBHbu3KnT9woPD4e1tTWcnJwwYsQIeHl5Yf/+/Wjfvr3W3qOurg4rV65Et27dYGNjAysrK9jY2OD333/Hb7/99ljbtLKywssvv4yMjAzF5FByuRz/93//h5iYGLi7uwMAMjMzIQgCXn75ZdTV1SkeXl5e6NmzZ4tOgDYQTWBAA36mmmdunykGhBl78cUX0bt3b7z11lu4f/++zt7nk08+wfHjx3H69GlcvXoVv/zyCwYMGKDV93jttdfw9ttvY/To0fjqq69w7NgxHD9+HD179kRVVdVjb3fKlCmorq7GZ599BgD45ptvUFpaqnQo4Pr16xBFEe3bt4e1tbXSIy8vD7du3WrR3yaTyVBWVgYfH58WbUcf+Jlqnjl9psx2sD4CBEHAmjVrMGTIEGzdulVn79O1a1c8/fTTja63s7NT/Jp6mCZfgpSUFEyYMAErV65U2UZL5uzo1q0b+vbti+3bt2P69OnYvn07fHx8MHToUEUbDw8PCIKAnJwcycHZWjpg29dffw25XC550tXY8DPVPHP6TLEHYeYGDx6MIUOGYOnSpaisbHp+ZV3p2LEjLly4oDTfd1lZGX788Ueldg1fCqlfb4IgqHxpvv76a5Wbp5raRmMmT56MY8eO4ejRo/jqq68wceJEpSHhR4wYAVEUUVJSgqefflrl0aNHD7Xf61FFRUV4/fXX4eLigunTpz/2dvSJn6nmmctnij2IVmDNmjUICwvDjRs38OSTT+r9/f/5z39iy5YtePnll/Hqq6+irKwMa9euVZmJy8nJCQEBAfjiiy/w7LPPws3NDR4eHujYsSNGjBiB5ORkhISEIDQ0FCdPnsS6devg5+entI3AwEDY29vj008/RdeuXeHo6AgfH58mu9rjxo3Da6+9hnHjxqGmpkZxGWSDAQMGYNq0aZg8eTJOnDiBgQMHwsHBAaWlpTh69Ch69OiBmTNnNrsfzp49qzjWfOPGDeTk5GD79u2wtLTEnj17NLpCx9D4mWoln6kWneImo9LUlQ3jx48XAejkipOmrqRosGPHDrFr166inZ2d2K1bN/Hzzz9XueJEFEXx4MGDYq9evURbW1sRgDhx4kRRFEXxzp074tSpU8V27dqJbdq0EZ955hkxJydH8mqWtLQ0MSQkRLS2tla6guXRK04e1rB/BgwY0OjfsG3bNrFfv36ig4ODaG9vLwYGBooTJkwQT5w40eTf3rCfGh42NjZiu3btxEGDBokrV64Ub9y40eTrDYmfqQda62eK80EQEZEknoMgIiJJDAgiIpLEgCAiIkkMCCIiksSAaOVu3bqFOXPmoGPHjrC1tUX79u3x3HPP4fbt24YuTe+4L7TjlVdeQffu3eHq6gp7e3t07twZ8+fPb/HdwabGHPYD74Noxa5evYqIiAhYWVnh7bffRnBwMG7duoWsrCzU1tYaujy94r7QHplMhmnTpiEoKAh2dnY4ceIEVqxYgX379uH06dOwsbExdIl6YQ77gZe5tmKjR4/GiRMnkJ+fj7Zt2xq6HIPivtCtzZs3Y9asWTh06BCio6MNXY7BmNp+4CEmI9YwGcm5c+cwbtw4uLi4oH379pgyZYrkODSaKCwsxJdffolXX33VJP5B5L7QHl3uy8Y03NFrZWU8By24H5rHgDABY8aMQefOnZGeno4FCxYgNTUVc+fOVayXy+VKQwY39nh4lqucnByIoggfHx+MGzcOjo6OsLOzQ2RkpEHn2G0O94X26GJfPqyurg4ymQw//PAD3n77bTzzzDNaH5FVG7gfmqCV+7FJJxpu41+7dq3S8lmzZol2dnZKk7bgodvuG3s0DDEgiqK4atUqEYDo7OwsxsTEiAcOHBDT09PF0NBQ0c7OTvz555/1+ac2i/tCe3S5Lxs0TIjT8Bg+fLhYXl6ujz9PbdwPzTONfk4rN2rUKKXnoaGhqK6uxo0bN9C+fXts2bIFFRUVzW7Hw8ND8d8Nv3b8/PyQnp6uGGnyb3/7G4KCgrB27VqkpKRo8a/QDu4L7dHFvmzQo0cPHD9+HPfu3cOZM2ewevVqDBkyBIcPH0abNm209jdoA/dD4xgQJqBhFqoGjw4/HBQUpNbsURYW/z2i2LDNwYMHKw1D7O3tjZ49e+LUqVMtrlsXuC+0Rxf7soGDg4NiPoeBAweiX79+CA8Px5YtW5QO3xgD7ofG8RyEGXj22WdVZqWSekyZMkXxmtDQ0Ea3J4qi5IfdFHBfaM/j7MvGPP3007CwsMCFCxf0ULl2teb9wB6EGXicLnC/fv3g5+eHb7/9FnK5XPHL+erVq/j5558xfvx4ndWrS9wX2tOSQyuP+v7771FfX4+goCBtlKZXrXk/MCDMQJcuXTR+jYWFBd555x384x//QExMDGbOnAmZTIZly5bBxsYGCxcu1EGlusd9oT2Psy8zMzPx0UcfYdSoUQgICMD9+/dx4sQJbNy4EUFBQXjllVd0UKluteb9wIBoxeLi4rBnzx6sWLECcXFxsLW1xaBBg/D5558jMDDQ0OXpFfeFdgQFBcHGxgbLli3D9evXATyYHnTq1KlYsGABXFxcDFyhfpjLfuCd1EREJKl1nn0jIqJmMSCIiEgSA4KIiCQxIIiISBIDgoiIJDEgiIhIEgOCiIgkMSCIiEgSA4KIiCQxIIiISBIDgoiIJDEgiIhI0v8HX6+22UnSZH4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p_deam_nn = []\n",
    "p_deam_nd = []\n",
    "p_norm_nn = []\n",
    "p_norm_nd = []\n",
    "\n",
    "for run in deam_df.run.unique():\n",
    "    deam_n = len(deam_df[deam_df.run == run][deam_df.aa == 'N'])\n",
    "    deam_nn = len(deam_df[deam_df.run == run][deam_df.aa == 'N'][deam_df['max'] < 1.3])\n",
    "    deam_nd = len(deam_df[deam_df.run == run][deam_df.aa == 'N'][deam_df['max'] >= 1.3])\n",
    "    # print(f\"{run} : {deam_nd/deam_n}, {deam_nn/deam_n}\")\n",
    "    p_deam_nn.append(deam_nn/deam_n*100)\n",
    "    p_deam_nd.append(deam_nd/deam_n*100)\n",
    "\n",
    "for run in norm_df[norm_df.aa == 'N'].run.unique():\n",
    "    norm_n = len(norm_df[norm_df.run == run][norm_df.aa == 'N'])\n",
    "    norm_nn = len(norm_df[norm_df.run == run][norm_df.aa == 'N'][norm_df['max'] < 1.3])\n",
    "    norm_nd = len(norm_df[norm_df.run == run][norm_df.aa == 'N'][norm_df['max'] >= 1.3])\n",
    "    # print(f\"{run} : {norm_nd/norm_n}, {norm_nn/norm_n}\")\n",
    "    p_norm_nn.append(norm_nn/norm_n*100)\n",
    "    p_norm_nd.append(norm_nd/norm_n*100)\n",
    "p_deam_nn, p_deam_nd, p_norm_nn, p_norm_nd\n",
    "\n",
    "bar_xs = [0,0.8,1.8,2.6]\n",
    "figure(figsize=(4,4))\n",
    "\n",
    "conditions = [p_norm_nn, p_norm_nd, p_deam_nn, p_deam_nd]\n",
    "bar_ys = [np.mean(x) for x in conditions]\n",
    "bar_errs = [np.std(x) for x in conditions]\n",
    "\n",
    "plt.bar(bar_xs, bar_ys, color = ['#33B77B','#643775','#33B77B', '#643775'])\n",
    "plt.xticks(bar_xs, \n",
    "        labels=[l + f\"\\nn={len(c)}\" for l, c in zip(['N', 'Putative D', 'N', 'Putative D'], conditions)],\n",
    "        fontsize=12)\n",
    "plt.yticks(fontsize=12)\n",
    "plt.errorbar(bar_xs, bar_ys, yerr=bar_errs, ls='none', color='black')\n",
    "plt.ylabel(\"% VRs in N position\", fontsize=12)\n",
    "ax = plt.gca()\n",
    "ax.spines[['right', 'top']].set_visible(False)\n",
    "new_rc_params = {'text.usetex': False,\n",
    "\"svg.fonttype\": 'none',\n",
    "\"font.family\":'arial'\n",
    "}\n",
    "\n",
    "for bar_x, ys in zip(bar_xs, conditions):\n",
    "    for y in ys:\n",
    "        plt.scatter(bar_x + random.random()/2-0.25, y, color = 'black')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.13 ('base')",
   "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.9.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "b01409a2fbaa3f35eaee56680c62c394066017bc60fbb802363dadef564d4717"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}