repeatmasker_to_bedtools-Copychr1libwhole_vslibchr1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pa\n",
    "import numpy as np\n",
    "import itertools\n",
    "import csv\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En tete a lire a chaque fois : les ID ne sont pas uniques, ça vient peut_etre des etoilllllllles je vais voir ce que je peux faire au pire j'en recrée.\n",
    "\n",
    "edit : \n",
    "non, ça vient du fait que les TE sont morcelés, et du coup pour le meme te on garde la meme ID meme si les séquences sont discontinues. du coup je fais quoi je sais pas.\n",
    "\n",
    "edit lundi 27 : \n",
    "j'ai recrée de nouvelle ID, ça devrait suffire pour faire fonctionner tout ça. \n",
    "(ça fonctionne effectivement, mais en l'état l'annotation n'est pas de bonne qualité. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "names_col_data = ['SW_score',\"perc_div\",\"perc_del\",\"perc_ins\", \"query_seq\", \"pos_query_begin\", \"pos_query_end\", \"pos_query_left\",\"strand\", \"matching_repeat\", \"repeat_class_family\", \"pos_in_repeat_begin\", \"pos_in_repeat_end\", \"pos_in_repeat_left\", \"TE_A_ID\",\"star\"]\n",
    "data = pa.read_csv(\"/home/tkastylevsky/FASTA_files/repeatmasker/gallus_gallus_proper/gallus_chr1_lib_whole_genome/chr1.fa.out\",header = None, skiprows = 2, sep = '\\s+', names = names_col_data)\n",
    "data['strand'] = data['strand'].replace('C','-')\n",
    "data.insert(1,\"A_ID\",list(range(1,len(data)+1)))\n",
    "\n",
    "data_split = data.repeat_class_family.str.split('/',expand = True)\n",
    "data = data.join(data_split)\n",
    "data.columns = ['SW_score','A_ID',\"perc_div\",\"perc_del\",\"perc_ins\", \"query_seq\", \"pos_query_begin\", \"pos_query_end\", \"pos_query_left\",\"strand\", \"matching_repeat\", \"repeat_class_family\", \"pos_in_repeat_begin\", \"pos_in_repeat_end\", \"pos_in_repeat_left\", \"TE_A_ID\",\"star\",\"repeat_class\",\"repeat_family\"]\n",
    "data = data.drop('repeat_class_family', axis = 1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_intersect = data\n",
    "data_intersect = data_intersect.drop([\"SW_score\",\"perc_del\",\"perc_ins\", \"pos_query_left\", \"pos_in_repeat_begin\", \"pos_in_repeat_end\", \"pos_in_repeat_left\", \"star\"], axis = 1 )\n",
    "data_intersect.insert(1, \"feature\", \"Similarity\")\n",
    "data_intersect.insert(1, \"source\", \"RepeatMasker\")\n",
    "data_intersect.insert(1, \"phase\", \".\")\n",
    "\n",
    "data_intersect = data_intersect[['query_seq','source','feature','pos_query_begin','pos_query_end', 'perc_div', 'strand','phase', 'A_ID']]\n",
    "data_intersect.to_csv('/home/tkastylevsky/FASTA_files/repeatmasker/gallus_gallus_proper/gallus_chr1_lib_whole_genome/chr1.fa.out_1.gff', sep = '\\t', header = False, index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12350\n"
     ]
    }
   ],
   "source": [
    "#nombre de hit qui se superposent à un meilleur hit au sein d'un meme out. \n",
    "print(data.star.count())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ici, la ref c'est mon annotation précédente\n",
    "\n",
    "\n",
    "\n",
    "names_col_ref = ['SW_score',\"perc_div\",\"perc_del\",\"perc_ins\", \"query_seq\", \"pos_query_begin\", \"pos_query_end\", \"pos_query_left\",\"strand\", \"matching_repeat\", \"repeat_class_family\", \"pos_in_repeat_begin\", \"pos_in_repeat_end\", \"pos_in_repeat_left\", \"TE_B_ID\",\"star\"]\n",
    "ref = pa.read_csv(\"/home/tkastylevsky/FASTA_files/repeatmasker/gallus_gallus_proper/gallus_proper_chr1/chr1.fa.preTueJan281534092020.RMoutput/chr1.fa.out\",header = None, skiprows = 2, sep = '\\s+', names = names_col_ref)\n",
    "ref['strand'] = ref['strand'].replace('C','-')\n",
    "ref_chr1 = ref[ref['query_seq']=='chr1']\n",
    "ref = ref_chr1\n",
    "ref.insert(1,\"B_ID\",list(range(1,len(ref)+1)))\n",
    "ref_split = ref.repeat_class_family.str.split('/',expand = True)\n",
    "ref = ref.join(ref_split)\n",
    "ref.columns = ['SW_score','B_ID',\"perc_div\",\"perc_del\",\"perc_ins\", \"query_seq\", \"pos_query_begin\", \"pos_query_end\", \"pos_query_left\",\"strand\", \"matching_repeat\", \"repeat_class_family\", \"pos_in_repeat_begin\", \"pos_in_repeat_end\", \"pos_in_repeat_left\", \"TE_B_ID\",\"star\",\"repeat_class\",\"repeat_family\"]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_intersect = ref\n",
    "ref_intersect = ref_intersect.drop([\"SW_score\",\"perc_del\",\"perc_ins\", \"pos_query_left\",\"repeat_class_family\", \"pos_in_repeat_begin\", \"pos_in_repeat_end\", \"pos_in_repeat_left\", \"star\"], axis = 1 )\n",
    "ref_intersect.insert(1, \"feature\", \"Similarity\")\n",
    "ref_intersect.insert(1, \"source\", \"RepeatMasker\")\n",
    "ref_intersect.insert(1, \"phase\", \".\")\n",
    "\n",
    "ref_intersect =ref_intersect[['query_seq','source','feature','pos_query_begin','pos_query_end', 'perc_div', 'strand','phase', 'B_ID']]\n",
    "ref_intersect.to_csv('/home/tkastylevsky/FASTA_files/repeatmasker/gallus_gallus_proper/gallus_proper_chr1/chr1.fa.out_1.gff', sep = '\\t', header = False, index  = False)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour passer a la suite, il faut faire un détour par bedtools pour récupérer les intersections(ça m'a semblé plus pratique)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#results A match_raw combien de fois A match dans B, sans tenir compte du type ni de la longueur des matchs. \n",
    "\n",
    "results = pa.read_csv(\"/home/tkastylevsky/results/results_bedtools/chk_galGal6_chr1_libwholegenome_compared_to_chr1_libchr1.txt\",header = None, sep = '\\s+')\n",
    "results = results.drop([0,1,2,5,6,7,9,10,11,14,15,16], axis = 1)\n",
    "names_col_2 = ['A_start','A_end','A_ID','B_start','B_end','B_ID','intersect']\n",
    "results.columns = names_col_2\n",
    "\n",
    "results['intersect_A'] = results['intersect']/(results['A_end']-results['A_start'])\n",
    "results['intersect_B'] = results['intersect']/(results['B_end']-results['B_start'])\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Deux petits graphes qui résument les pourcentages d'identité entre les matchs. J'ai ai fais un en logscale "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAD4CAYAAAAgs6s2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAW6klEQVR4nO3dfbRddX3n8ffHRARaEZToMAk0OKaOkdUuMYN0nGmtUQy0JXQWdMKaDtFmzOhgax/WjGC7BpfKWtgnKmt8wpISqBWRtpKpMJmIUGdmCRKKw6MMd8CBWxiJBtGKQoPf+eP8rnNMzr33JHefc7k379daZ2Xv7/7tvX+/3Es+7MeTqkKSpC49Z747IElafAwXSVLnDBdJUucMF0lS5wwXSVLnls53B54tjj766Fq5cuV8d0OSFpTbbrvt61W1bO+64dKsXLmSnTt3znc3JGlBSfJ/BtU9LSZJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSercyMIlyZYkjyW5q6/2e0m+kuSOJH+Z5Mi+ZecnmUhyX5I39dXXtdpEkvP66scnuSXJ/Uk+leSQVn9em59oy1eOaoyStBCtPO+zP/QZhVEeuVwOrNurtgM4oap+AvhfwPkASVYDG4BXtnU+nGRJkiXAh4BTgdXA2a0twAeAi6tqFfA4sKnVNwGPV9XLgItbO0nSGI0sXKrqC8DuvWr/tar2tNmbgRVtej1wVVU9VVUPAhPASe0zUVUPVNXTwFXA+iQBXg9c09bfCpzRt62tbfoaYG1rL0kak/m85vIrwPVtejnwcN+yyVabrv4i4Jt9QTVV/6FtteVPtPb7SLI5yc4kO3ft2jXnAUmSeuYlXJL8NrAH+MRUaUCzOoD6TNvat1h1aVWtqao1y5bt88ZoSdIBGvsr95NsBH4eWFtVU//oTwLH9jVbATzSpgfVvw4cmWRpOzrpbz+1rckkS4EXsNfpOUnSaI31yCXJOuBdwOlV9WTfom3Ahnan1/HAKuBLwK3AqnZn2CH0Lvpva6F0I3BmW38jcG3ftja26TOBz/eFmCRpDEZ25JLkk8DrgKOTTAIX0Ls77HnAjnaN/eaqeltV3Z3kauAeeqfLzq2qZ9p23gFsB5YAW6rq7raLdwFXJXk/cDtwWatfBlyZZILeEcuGUY1RkjTYyMKlqs4eUL5sQG2q/YXAhQPq1wHXDag/QO9usr3r3wPO2q/OSpI65RP6kqTOGS6SpM4ZLpKkzhkukqTOGS6SpM4ZLpKkzhkukqTOGS6SpM4ZLpKkzhkukqTOGS6SpM4ZLpKkzhkukqTOGS6SpM4ZLpKkzhkukqTOGS6SpM4ZLpKkzhkukqTOGS6SpM4ZLpKkzhkukqTOGS6SpM4ZLpKkzhkukqTOjSxckmxJ8liSu/pqL0yyI8n97c+jWj1JLkkykeSOJCf2rbOxtb8/yca++quT3NnWuSRJZtqHJGl8Rnnkcjmwbq/aecANVbUKuKHNA5wKrGqfzcBHoBcUwAXAa4CTgAv6wuIjre3Ueutm2YckaUxGFi5V9QVg917l9cDWNr0VOKOvfkX13AwcmeQY4E3AjqraXVWPAzuAdW3ZEVX1xaoq4Iq9tjVoH5KkMRn3NZeXVNWjAO3PF7f6cuDhvnaTrTZTfXJAfaZ97CPJ5iQ7k+zctWvXAQ9KkvTDni0X9DOgVgdQ3y9VdWlVramqNcuWLdvf1SVJ0xh3uHytndKi/flYq08Cx/a1WwE8Mkt9xYD6TPuQJI3JuMNlGzB1x9dG4Nq++jntrrGTgSfaKa3twClJjmoX8k8Btrdl305ycrtL7Jy9tjVoH5KkMVk6qg0n+STwOuDoJJP07vq6CLg6ySbgIeCs1vw64DRgAngSeAtAVe1O8j7g1tbuvVU1dZPA2+ndkXYYcH37MMM+JEljMrJwqaqzp1m0dkDbAs6dZjtbgC0D6juBEwbUvzFoH5Kk8Xm2XNCXJC0ihoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlz8xIuSX4jyd1J7kryySSHJjk+yS1J7k/yqSSHtLbPa/MTbfnKvu2c3+r3JXlTX31dq00kOW/8I5Skg9vYwyXJcuDXgDVVdQKwBNgAfAC4uKpWAY8Dm9oqm4DHq+plwMWtHUlWt/VeCawDPpxkSZIlwIeAU4HVwNmtrSRpTObrtNhS4LAkS4HDgUeB1wPXtOVbgTPa9Po2T1u+Nkla/aqqeqqqHgQmgJPaZ6KqHqiqp4GrWltJ0piMPVyq6m+B3wceohcqTwC3Ad+sqj2t2SSwvE0vBx5u6+5p7V/UX99rnenq+0iyOcnOJDt37do198FJkoD5OS12FL0jieOBfwj8CL1TWHurqVWmWba/9X2LVZdW1ZqqWrNs2bLZui5JGtJ8nBZ7A/BgVe2qqr8H/gL4p8CR7TQZwArgkTY9CRwL0Ja/ANjdX99rnenqkqQxGSpckpzQ4T4fAk5Ocni7drIWuAe4ETiztdkIXNumt7V52vLPV1W1+oZ2N9nxwCrgS8CtwKp299kh9C76b+uw/5KkWSydvQkAH23/UF8O/FlVffNAd1hVtyS5BvgbYA9wO3Ap8FngqiTvb7XL2iqXAVcmmaB3xLKhbefuJFfTC6Y9wLlV9QxAkncA2+ndibalqu4+0P5KkvbfUOFSVf8sySrgV4CdSb4E/ElV7TiQnVbVBcAFe5UfoHen195tvwecNc12LgQuHFC/DrjuQPomSZq7oa+5VNX9wO8A7wJ+BrgkyVeS/ItRdU6StDANe83lJ5JcDNxL73mUX6iqV7Tpi0fYP0nSAjTsNZf/BHwceHdVfXeqWFWPJPmdkfRMkrRgDRsupwHf7btg/hzg0Kp6sqquHFnvJEkL0rDXXD4HHNY3f3irSZK0j2HD5dCq+rupmTZ9+Gi6JEla6IYNl+8kOXFqJsmrge/O0F6SdBAb9prLrwOfTjL1GpVjgH85mi5Jkha6YR+ivDXJPwZeTu/FkF9p7wWTJGkfwx65APwTYGVb51VJqKorRtIrSdKCNlS4JLkS+EfAl4FnWrkAw0WStI9hj1zWAKvb24glSZrRsHeL3QX8g1F2RJK0eAx75HI0cE97G/JTU8WqOn0kvZIkLWjDhst7RtkJSdLiMuytyH+d5MeAVVX1uSSH0/siLkmS9jHsK/ffClwDfKyVlgOfGVWnJEkL27AX9M8FXgt8C37wxWEvHlWnJEkL27Dh8lRVPT01k2QpvedcJEnax7Dh8tdJ3g0cluSNwKeB/zy6bkmSFrJhw+U8YBdwJ/BvgesAv4FSkjTQsHeLfZ/e1xx/fLTdkSQtBsO+W+xBBlxjqaqXdt4jSdKCtz/vFptyKHAW8MLuuyNJWgyGuuZSVd/o+/xtVf0R8PoR902StEAN+xDliX2fNUneBjz/QHea5Mgk1yT5SpJ7k/xUkhcm2ZHk/vbnUa1tklySZCLJHXt93fLG1v7+JBv76q9Ocmdb55IkOdC+SpL237Cnxf6gb3oP8FXgl+aw3w8C/6WqzkxyCHA48G7ghqq6KMl59O5QexdwKrCqfV4DfAR4TZIXAhfQO2VXwG1JtlXV463NZuBmene2rQOun0N/JUn7Ydi7xX62qx0mOQL4aeDNbdtPA08nWQ+8rjXbCtxEL1zWA1e075K5uR31HNPa7qiq3W27O4B1SW4CjqiqL7b6FcAZGC6SNDbD3i32mzMtr6o/3I99vpTeMzN/kuQngduAdwIvqapH2/YeTTL1epnlwMN960+22kz1yQH1fSTZTO8Ih+OOO24/hiBJmsmwD1GuAd7O///H+23AanrXXfb32stS4ETgI1X1KuA79E6BTWfQ9ZI6gPq+xapLq2pNVa1ZtmzZzL2WJA1tf74s7MSq+jZAkvcAn66qf3MA+5wEJqvqljZ/Db1w+VqSY9pRyzHAY33tj+1bfwXwSKu/bq/6Ta2+YkB7SdKYDHvkchzwdN/808DKA9lhVf1f4OEkL2+ltcA9wDZg6o6vjcC1bXobcE67a+xk4Il2+mw7cEqSo9qdZacA29uybyc5ud0ldk7ftiRJYzDskcuVwJeS/CW9U0y/CFwxh/3+KvCJdqfYA8Bb6AXd1Uk2AQ/Re1ATend7nQZMAE+2tlTV7iTvA25t7d47dXGf3im8y4HD6F3I92K+JI3RsHeLXZjkeuCft9Jbqur2A91pVX2ZH37qf8raAW2L3vfJDNrOFmDLgPpO4IQD7Z8kaW6GPS0GvWdRvlVVHwQmkxw/oj5Jkha4YZ/Qv4DeMyfnt9JzgT8dVackSQvbsEcuvwicTu+2YarqEebw+hdJ0uI2bLg83a59FECSHxldlyRJC92w4XJ1ko8BRyZ5K/A5/OIwSdI0hr1b7PeTvBH4FvBy4D9W1Y6R9kyStGDNGi5JltB7OPENgIEiSZrVrKfFquoZ4MkkLxhDfyRJi8CwT+h/D7izvdb+O1PFqvq1kfRKkrSgDRsun20fSZJmNWO4JDmuqh6qqq3j6pAkaeGb7ZrLZ6Ymkvz5iPsiSVokZguX/i/eeukoOyJJWjxmC5eaZlqSpGnNdkH/J5N8i94RzGFtmjZfVXXESHsnSVqQZgyXqloyro5IkhaP/fk+F0mShmK4SJI6Z7hIkjpnuEiSOme4SJI6Z7hIkjpnuEiSOme4SJI6Z7hIkjo3b+GSZEmS25P8VZs/PsktSe5P8qkkh7T689r8RFu+sm8b57f6fUne1Fdf12oTSc4b99gk6WA3n0cu7wTu7Zv/AHBxVa0CHgc2tfom4PGqehlwcWtHktXABuCVwDrgwy2wlgAfAk4FVgNnt7aSpDGZl3BJsgL4OeCP23yA1wPXtCZbgTPa9Po2T1u+trVfD1xVVU9V1YPABHBS+0xU1QNV9TRwVWsrSRqT+Tpy+SPgPwDfb/MvAr5ZVXva/CSwvE0vBx4GaMufaO1/UN9rnenq+0iyOcnOJDt37do11zFJkpqxh0uSnwceq6rb+ssDmtYsy/a3vm+x6tKqWlNVa5YtWzZDryVJ+2O273MZhdcCpyc5DTgUOILekcyRSZa2o5MVwCOt/SRwLDCZZCnwAmB3X31K/zrT1SVJYzD2I5eqOr+qVlTVSnoX5D9fVf8KuBE4szXbCFzbpre1edryz1dVtfqGdjfZ8cAq4EvArcCqdvfZIW0f28YwNElSMx9HLtN5F3BVkvcDtwOXtfplwJVJJugdsWwAqKq7k1wN3APsAc6tqmcAkrwD2A4sAbZU1d1jHYkkHeTmNVyq6ibgpjb9AL07vfZu8z3grGnWvxC4cED9OuC6DrsqSdoPPqEvSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSerc2MMlybFJbkxyb5K7k7yz1V+YZEeS+9ufR7V6klySZCLJHUlO7NvWxtb+/iQb++qvTnJnW+eSJBn3OCXpYDYfRy57gN+qqlcAJwPnJlkNnAfcUFWrgBvaPMCpwKr22Qx8BHphBFwAvAY4CbhgKpBam819660bw7gkSc3Yw6WqHq2qv2nT3wbuBZYD64GtrdlW4Iw2vR64onpuBo5McgzwJmBHVe2uqseBHcC6tuyIqvpiVRVwRd+2JEljMK/XXJKsBF4F3AK8pKoehV4AAS9uzZYDD/etNtlqM9UnB9QH7X9zkp1Jdu7atWuuw5EkNfMWLkl+FPhz4Ner6lszNR1QqwOo71usurSq1lTVmmXLls3WZUnSkOYlXJI8l16wfKKq/qKVv9ZOadH+fKzVJ4Fj+1ZfATwyS33FgLokaUzm426xAJcB91bVH/Yt2gZM3fG1Ebi2r35Ou2vsZOCJdtpsO3BKkqPahfxTgO1t2beTnNz2dU7ftiRJY7B0Hvb5WuBfA3cm+XKrvRu4CLg6ySbgIeCstuw64DRgAngSeAtAVe1O8j7g1tbuvVW1u02/HbgcOAy4vn0kSWMy9nCpqv/O4OsiAGsHtC/g3Gm2tQXYMqC+EzhhDt2UJM2BT+hLkjpnuEiSOme4SJI6Z7hIkjpnuEiSOme4SJI6Z7hIkjpnuEiSOme4SJI6Z7hIkjpnuEiSOme4SJI6Z7hIkjpnuEiSOme4SJI6Z7hIkjpnuEiSOme4SJI6Z7hIkjpnuEiSOme4SJI6t3S+OyBJGr2V5312rPszXDrQ/0P76kU/N489kaRnB0+LSZI655GLJC1S4z4V1m/RhkuSdcAHgSXAH1fVRfPcJUkaifkMkeksynBJsgT4EPBGYBK4Ncm2qrpn1Pse5ofsdRnp4DXdNdpnY0DMxaIMF+AkYKKqHgBIchWwHhh5uAxjsf0S9ZvuPxYDVVMW8+///lrMfxeLNVyWAw/3zU8Cr9m7UZLNwOY2+3dJ7jvA/R0NfP0A112oBo45HxjceLr6AuPP+eBw0I05H5jTmH9sUHGxhksG1GqfQtWlwKVz3lmys6rWzHU7C4ljPjg45oPDKMa8WG9FngSO7ZtfATwyT32RpIPOYg2XW4FVSY5PcgiwAdg2z32SpIPGojwtVlV7krwD2E7vVuQtVXX3CHc551NrC5BjPjg45oND52NO1T6XIiRJmpPFelpMkjSPDBdJUucMl/2QZF2S+5JMJDlvwPLnJflUW35LkpXj72W3hhjzbya5J8kdSW5IMvCe94VktjH3tTszSSVZ0LetDjPeJL/Ufs53J/mzcfexa0P8Xh+X5MYkt7ff7dPmo59dSrIlyWNJ7ppmeZJc0v5O7khy4px2WFV+hvjQuzHgfwMvBQ4B/ieweq82/w74aJveAHxqvvs9hjH/LHB4m377wTDm1u75wBeAm4E1893vEf+MVwG3A0e1+RfPd7/HMOZLgbe36dXAV+e73x2M+6eBE4G7pll+GnA9vecETwZumcv+PHIZ3g9eKVNVTwNTr5Tptx7Y2qavAdYmGfRA50Ix65ir6saqerLN3kzvmaKFbJifM8D7gN8FvjfOzo3AMON9K/ChqnocoKoeG3MfuzbMmAs4ok2/gEXwnFxVfQHYPUOT9cAV1XMzcGSSYw50f4bL8Aa9Umb5dG2qag/wBPCisfRuNIYZc79N9P7PZyGbdcxJXgUcW1V/Nc6OjcgwP+MfB348yf9IcnN74/hCNsyY3wP8cpJJ4DrgV8fTtXm1v/+9z2hRPucyIsO8Umao184sIEOPJ8kvA2uAnxlpj0ZvxjEneQ5wMfDmcXVoxIb5GS+ld2rsdfSOTP9bkhOq6psj7tuoDDPms4HLq+oPkvwUcGUb8/dH37150+m/Xx65DG+YV8r8oE2SpfQOp2c6DH22G+o1OkneAPw2cHpVPTWmvo3KbGN+PnACcFOSr9I7N71tAV/UH/b3+tqq+vuqehC4j17YLFTDjHkTcDVAVX0ROJTeCy0Xs05fm2W4DG+YV8psAza26TOBz1e7UrZAzTrmdoroY/SCZaGfi4dZxlxVT1TV0VW1sqpW0rvOdHpV7Zyf7s7ZML/Xn6F34wZJjqZ3muyBsfayW8OM+SFgLUCSV9ALl11j7eX4bQPOaXeNnQw8UVWPHujGPC02pJrmlTJJ3gvsrKptwGX0Dp8n6B2xbJi/Hs/dkGP+PeBHgU+3exceqqrT563TczTkmBeNIce7HTglyT3AM8C/r6pvzF+v52bIMf8W8PEkv0Hv1NCbF/j/KJLkk/RObR7driVdADwXoKo+Su/a0mnABPAk8JY57W+B/31Jkp6FPC0mSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSerc/wOt4F4+a4oVbAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot = results.intersect_A.plot.hist(bins = 100)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD4CAYAAAAD6PrjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAARzElEQVR4nO3de4yld13H8fenraUUcVEL0bRdtri1siEa61i8i4q6XLZVRGzVKKZ2Fa3GSwz1EkWNSeMNJTQpqzYVDK0Fte7KkioqlpgiLeClFxvXWulYY0F0UUFq4esf5/TJcZzZeaYzv/PMc877lUxynufcvr+d2fM5v8vzPKkqJEkCOG3oAiRJu4ehIEnqGAqSpI6hIEnqGAqSpM4ZQxewHeecc07t27dv6DIkaVTe/e53f6Cqnr7efaMOhX379nHXXXcNXYYkjUqSf9zoPoePJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1Bn1wWuStEz2XfOW7vaD176oyXvYU5AkdXZNKCR5XpJ3JLk+yfOGrkeSllHTUEhyQ5JHkty9Zv/BJPcnOZHkmunuAv4TOAtYbVmXJGl9rXsKNwIHZ3ckOR24DngBcAC4IskB4B1V9QLglcBPN65LkrSOpqFQVbcDH1yz+xLgRFU9UFWPAjcDl1XVx6f3/xvwpJZ1SZLWN8Tqo3OBh2a2V4HnJnkJ8LXA04DXbvTkJIeBwwB79+5tWKYkLZ8hQiHr7Kuq+l3gdzd7clUdAY4ArKys1A7XJklLbYjVR6vA+TPb5wEPD1CHJGmNIULhTuDCJBckORO4HDi6lRdIcijJkZMnTzYpUJKWVeslqTcBdwAXJVlNcmVVPQZcDdwG3AfcUlX3bOV1q+pYVR3es2fPzhctSUus6ZxCVV2xwf7jwPGW7y1J2rpdc0TzVjh8JEltjDIUHD6SpDZGGQqSpDYMBUlSZ5Sh4JyCJLUxylBwTkGS2hhlKEiS2jAUJEmdUYaCcwqS1MYoQ8E5BUlqY5ShIElqw1CQJHUMBUlSZ5Sh4ESzJLUxylBwolmS2hhlKEiS2jAUJEkdQ0GS1DEUJEmdUYaCq48kqY1RhoKrjySpjVGGgiSpDUNBktQxFCRJHUNBktQxFCRJHUNBktQZZSh4nIIktTHKUPA4BUlqY5ShIElqw1CQJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSZ5Sh4GkuJKmNUYaCp7mQpDZGGQqSpDYMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSZ1eFQpKnJHl3khcPXYskLaOmoZDkhiSPJLl7zf6DSe5PciLJNTN3vRK4pWVNkqSNte4p3AgcnN2R5HTgOuAFwAHgiiQHkjwfuBf4l8Y1SZI2cEbLF6+q25PsW7P7EuBEVT0AkORm4DLgE4GnMAmKjyQ5XlUfX/uaSQ4DhwH27t3brnhJWkJNQ2ED5wIPzWyvAs+tqqsBkrwc+MB6gQBQVUeAIwArKyvVtlRJWi5DhELW2dd9uFfVjfMrRZI0a4jVR6vA+TPb5wEPb+UFkhxKcuTkyZM7WpgkLbshQuFO4MIkFyQ5E7gcOLqVF6iqY1V1eM+ePU0KlKRl1XpJ6k3AHcBFSVaTXFlVjwFXA7cB9wG3VNU9LeuQJPXTa04hyXOq6u7NH/l/VdUVG+w/Dhzf6uvN1HMIOLR///4n+hKSpHX07Slcn+RdSb4nydOaVtSDw0eS1EavUKiqLwG+hckE8V1J3pjkq5tWJkmau95zClX1d8BPMDkVxZcDr0nyt0le0qo4SdJ89QqFJJ+d5NVMJoa/EjhUVc+e3n51w/o2qsclqZLUQN+ewmuB9wCfU1XfW1XvAaiqh5n0HubKOQVJaqPvEc0vBD5SVR8DSHIacFZVfbiq3tCsOknSXPXtKbwNePLM9tnTfZKkBdI3FM6qqv98fGN6++w2JW3OOQVJaqNvKPxXkosf30jyecBH2pS0OecUJKmNvnMKPwC8KcnjJ677dOCb2pQkSRpKr1CoqjuTfBZwEZNTX/9tVf1P08okSXO3lespfD6wb/qcz01CVb2+SVWSpEH0PSHeG4DPAP4S+Nh0dwGDhIInxJOkNvr2FFaAA1W1Ky5/WVXHgGMrKytXDV2LJC2SvquP7gY+rWUhkqTh9e0pnAPcm+RdwEcf31lVlzapSpI0iL6h8KqWRUiSdoe+S1L/LMkzgQur6m1JzgZOb1uaJGne+p46+yrgzcDrprvOBW5tVVSPejzNhSQ10Hei+XuBLwY+BN0Fd57RqqjNeJoLSWqj75zCR6vq0SQAJDmDyXEKkqSG9l3zlrm+X9+ewp8l+THgydNrM78JONauLEnSEPqGwjXA+4G/Ab4LOM4AV1yTJLXVd/XRx4Ffm/5IkhZU33Mf/QPrzCFU1bN2vCJJ0mC2cu6jx50FfCPwKTtfzvzMTt48eO2LBqxEknaPXnMKVfWvMz//VFW/Anxl49o25HEKktRG3+Gji2c2T2PSc3hqk4p68CypktRG3+GjX5q5/RjwIPCyHa9GkjSovquPvqJ1IZKk4fUdPvqhU91fVb+8M+VIkoa0ldVHnw8cnW4fAm4HHmpRlCRpGFu5yM7FVfUfAEleBbypqr6zVWGSpPnrGwp7gUdnth8F9u14NZK05OZ9Ary1+obCG4B3Jfk9Jkc2fz3w+mZVSZIG0Xf10c8leSvwpdNd31FV721XliRpCH3PkgpwNvChqvpVYDXJBY1qkiQNpO/lOH8KeCXwo9NdnwD8VquiJEnD6NtT+HrgUuC/AKrqYQY8zYXnPpKkNvqGwqNVVUxPn53kKe1K2pzXaJakNvqGwi1JXgc8LclVwNvwgjuStHD6rj76xem1mT8EXAT8ZFX9UdPKJGlJDH1swqxNQyHJ6cBtVfV8wCCQpAW26fBRVX0M+HASB/AlacH1PaL5v4G/SfJHTFcgAVTV9zepSpI0iL6h8JbpjyRpgZ0yFJLsrar3VdVvzqsgSdJwNusp3ApcDJDkd6rqG9qXJEmLbzetOJq12URzZm4/q2UhkqThbRYKtcFtSdIC2mz46HOSfIhJj+HJ09tMt6uqPqlpdZKkuTplKFTV6fMqZEizY3sPXvuiASuRpGH1XZIqSdqm3Tq5PGsrF9lpKsmzk1yf5M1JXjF0PZK0jJqGQpIbkjyS5O41+w8muT/JiSTXAFTVfVX13cDLgJWWdUmS1td6+OhG4LXA6x/fMT3B3nXAVwOrwJ1JjlbVvUkuBa6ZPkeSRm8MQ0azmvYUqup24INrdl8CnKiqB6rqUeBm4LLp449W1RcB39KyLknS+oaYaD4XeGhmexV4bpLnAS8BngQc3+jJSQ4DhwH27t3brkpJWkJDhELW2VdV9Xbg7Zs9uaqOAEcAVlZWPKBO0q4ztiGjWUOEwipw/sz2ecDDA9SxLo9ZkLTMhgiFO4ELk1wA/BNwOfDNW3mBJIeAQ/v3729QniRt3Zh7B7NaL0m9CbgDuCjJapIrq+ox4GrgNuA+4Jaqumcrr1tVx6rq8J49XgxOknZS055CVV2xwf7jnGIyWZI0jFGe5mJew0fOL0iCjT8LFmXIaNYoQ6GqjgHHVlZWrhq6FknLZRGDYNauOfeRJGl4hoIkqTPK4SOXpEpqbdGHiTYyyp6CS1IlqY1RhoIkqY1RDh8NweWp0uJb1iGjWaPsKSQ5lOTIyZMnhy5FkhbKKHsKHqcgaafYO/i/RhkKQ1v7R+RwkqRFYShIWjr2DjZmKOwAJ6ElLYpRhsJuPnjNgJA0ZqNcfeTBa5LUxih7CmNhr0HaPZxH6MdQGIBhIc2HQbB1hsKc+McpaQwMBUmjZ+9754wyFHbz6qOt8o9Z0m7i6iNJUmeUPQVJ2ojzd9szyp6CJKkNewq7yEbfcJxrkDQv9hQkSR17CiPgCiVJ82JPQZLUGWVPYZGOU9gqew2SWhplKHg5zgkDQtJOG2UoSFo8fb7k+EWoPUNhQfifRdJOMBQk7WobHb/jkcttGAqSdh0/8IfjklRJUseewgJyfkHSE2UoLBHDQtJmDAVJW+KXi8VmKCw4J+zUkgGxeEYZCst8mgtpt/Lgs8UwylDwNBfS7mYPdbxGGQpqx29yWk+LD3mDY3cyFJaUH/6S1mMoqBcvFbp8/Ca/nAwF7di5ZdY+3sAYhhO+2g5DQc3Yu5DGx1DQoPzGKu0uhoK2xXFnabEYCpI6hrwMBY3Wbh562s21SadiKGjuWh8ItWgrbnaqbnsB6sNQ0KgsywfbPNu5LP+m6sdQ0ELrcwzGmHoNs7baO5L6MBS0a4zlg3qrdfY9qG87H+B++GunGAralYb6kBtLMEmt7KpQSPJ1wIuAZwDXVdUfDlySltg8AsJv+NptmodCkhuAFwOPVNVzZvYfBH4VOB349aq6tqpuBW5N8snALwKGgkZjHh/whohaO20O73EjcHB2R5LTgeuAFwAHgCuSHJh5yE9M75ckzVHznkJV3Z5k35rdlwAnquoBgCQ3A5cluQ+4FnhrVb1nvddLchg4DLB3795WZWuJ9Pn27Td0LYt59BTWcy7w0Mz26nTf9wHPB16a5LvXe2JVHamqlapaefrTn96+UklaIkNNNGedfVVVrwFeM+9itFh2+4Ff9jq0mw3VU1gFzp/ZPg94uO+TkxxKcuTkyZM7XpgkLbOhQuFO4MIkFyQ5E7gcONr3yVV1rKoO79mzp1mBkrSMmodCkpuAO4CLkqwmubKqHgOuBm4D7gNuqap7WtciSTq1eaw+umKD/ceB40/kNZMcAg7t379/O6VJktYYavhoWxw+kqQ2RhkKkqQ2DAVJUmeUoeCSVElqY5Sh4JyCJLWRqhq6hicsyfuBf3yCTz8H+MAOljMGtnk52OblsJ02P7Oq1j1P0KhDYTuS3FVVK0PXMU+2eTnY5uXQqs2jHD6SJLVhKEiSOsscCkeGLmAAtnk52Obl0KTNSzunIEn6/5a5pyBJWsNQkCR1Fj4UkhxMcn+SE0muWef+JyX57en9f7HO9aRHp0ebfyjJvUn+OskfJ3nmEHXupM3aPPO4lyapJKNfvtinzUleNv1d35PkjfOucaf1+Nvem+RPk7x3+vf9wiHq3ClJbkjySJK7N7g/SV4z/ff46yQXb/tNq2phf4DTgb8HngWcCfwVcGDNY74HuH56+3Lgt4euew5t/grg7OntVyxDm6ePeypwO/BOYGXouufwe74QeC/wydPtZwxd9xzafAR4xfT2AeDBoeveZpu/DLgYuHuD+18IvJXJJY6/APiL7b7novcULgFOVNUDVfUocDNw2ZrHXAb85vT2m4GvSrLeNaTHYtM2V9WfVtWHp5vvZHI51DHr83sG+Fng54H/nmdxjfRp81XAdVX1bwBV9cica9xpfdpcwCdNb+9hC5f53Y2q6nbgg6d4yGXA62vincDTknz6dt5z0UPhXOChme3V6b51H1OTK8KdBD51LtW10afNs65k8k1jzDZtc5LPBc6vqj+YZ2EN9fk9fybwmUn+PMk7kxycW3Vt9Gnzq4BvTbLK5CJe3zef0gaz1f/vm2p+5bWBrfeNf+0a3D6PGZPe7UnyrcAK8OVNK2rvlG1OchrwauDl8ypoDvr8ns9gMoT0PCa9wXckeU5V/Xvj2lrp0+YrgBur6peSfCHwhmmbP96+vEHs+OfXovcUVoHzZ7bP4/93J7vHJDmDSZfzVN213a5Pm0nyfODHgUur6qNzqq2Vzdr8VOA5wNuTPMhk7PXoyCeb+/5t/35V/U9V/QNwP5OQGKs+bb4SuAWgqu4AzmJy4rhF1ev/+1YseijcCVyY5IIkZzKZSD665jFHgW+f3n4p8Cc1ncEZqU3bPB1KeR2TQBj7ODNs0uaqOllV51TVvqrax2Qe5dKqumuYcndEn7/tW5ksKiDJOUyGkx6Ya5U7q0+b3wd8FUCSZzMJhffPtcr5Ogp823QV0hcAJ6vqn7fzggs9fFRVjyW5GriNycqFG6rqniQ/A9xVVUeB32DSxTzBpIdw+XAVb1/PNv8C8InAm6Zz6u+rqksHK3qberZ5ofRs823A1yS5F/gY8CNV9a/DVb09Pdv8w8CvJflBJsMoLx/zl7wkNzEZ/jtnOk/yU8AnAFTV9UzmTV4InAA+DHzHtt9zxP9ekqQdtujDR5KkLTAUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1PlfDtm64kNsINYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot = results.intersect_A.plot.hist(bins = 100)\n",
    "plot = plot.set_yscale('log')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAD4CAYAAAAgs6s2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAW6UlEQVR4nO3dfbRddZ3f8ffHRARm5EGJ1knAYM3YiazpElNkajt1jGLAKaFrwTSsmSHSjOk46DyuVnC6GpfKKnYeGFn1CSVjoI6IjJV0hKYRcWy75CGI5Uktt2DhDlSiQbSiMMFv/zi/6zom5+ae3LvPudyb92uts+7e3/3be/9+3CQf9t6/c06qCkmSuvSs+e6AJGnxMVwkSZ0zXCRJnTNcJEmdM1wkSZ1bOt8deKY47rjjauXKlfPdDUlaUG6//fZvVdWyfeuGS7Ny5Up27do1392QpAUlyf8ZVPe2mCSpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXO+Q1+SDjErL/zsT6x/45I3dn4Or1wkSZ0zXCRJnRtZuCTZmuTRJHf31f4oydeS3JnkPyU5pm/bRUkmknw9yRv66utabSLJhX31E5PckuS+JJ9MclirP6etT7TtK0c1RknSYKO8cvkYsG6f2k7gpKr6eeB/ARcBJFkNbABe3vb5QJIlSZYA7wdOB1YD57a2AO8FLq2qVcBjwKZW3wQ8VlUvBS5t7SRJYzSycKmqLwJ79qn916ra21ZvBla05fXA1VX1ZFU9AEwAp7TXRFXdX1VPAVcD65MEeC1wbdt/G3BW37G2teVrgbWtvSRpTObzmcu/AG5oy8uBh/q2TbbadPXnA9/pC6qp+k8cq21/vLXfT5LNSXYl2bV79+45D0iS1DMv4ZLkD4G9wMenSgOa1SzqBzrW/sWqy6tqTVWtWbZsvy9SkyTN0tjf55JkI/DLwNqqmvpHfxI4vq/ZCuDhtjyo/i3gmCRL29VJf/upY00mWQoczT635yRJozXWK5ck64C3A2dW1RN9m7YDG9pMrxOBVcCtwG3AqjYz7DB6D/23t1C6CTi77b8RuK7vWBvb8tnA5/tCTJI0BiO7cknyCeA1wHFJJoEt9GaHPQfY2Z6x31xVv1lV9yS5BriX3u2yC6rq6XactwI7gCXA1qq6p53i7cDVSd4D3AFc0epXAFclmaB3xbJhVGOUJA02snCpqnMHlK8YUJtqfzFw8YD69cD1A+r305tNtm/9h8A5B9VZSVKnfIe+JKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpc4aLJKlzhoskqXOGiySpcyMLlyRbkzya5O6+2vOS7ExyX/t5bKsnyWVJJpLcmeTkvn02tvb3JdnYV39lkrvaPpclyYHOIUkan1FeuXwMWLdP7ULgxqpaBdzY1gFOB1a112bgg9ALCmAL8CrgFGBLX1h8sLWd2m/dDOeQJI3JyMKlqr4I7NmnvB7Y1pa3AWf11a+snpuBY5K8CHgDsLOq9lTVY8BOYF3bdlRVfamqCrhyn2MNOockaUzG/czlhVX1CED7+YJWXw481NdustUOVJ8cUD/QOSRJY/JMeaCfAbWaRf3gTppsTrIrya7du3cf7O6SpGmMO1y+2W5p0X4+2uqTwPF97VYAD89QXzGgfqBz7KeqLq+qNVW1ZtmyZbMelCTpJ407XLYDUzO+NgLX9dXPa7PGTgUeb7e0dgCnJTm2Pcg/DdjRtn0vyaltlth5+xxr0DkkSWOydFQHTvIJ4DXAcUkm6c36ugS4Jskm4EHgnNb8euAMYAJ4AjgfoKr2JHk3cFtr966qmpok8BZ6M9KOAG5oLw5wDknSmIwsXKrq3Gk2rR3QtoALpjnOVmDrgPou4KQB9W8POockaXyeKQ/0JUmLiOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknq3LyES5LfS3JPkruTfCLJ4UlOTHJLkvuSfDLJYa3tc9r6RNu+su84F7X615O8oa++rtUmklw4/hFK0qFt7OGSZDnw28CaqjoJWAJsAN4LXFpVq4DHgE1tl03AY1X1UuDS1o4kq9t+LwfWAR9IsiTJEuD9wOnAauDc1laSNCbzdVtsKXBEkqXAkcAjwGuBa9v2bcBZbXl9W6dtX5skrX51VT1ZVQ8AE8Ap7TVRVfdX1VPA1a2tJGlMxh4uVfU3wB8DD9ILlceB24HvVNXe1mwSWN6WlwMPtX33tvbP76/vs8909f0k2ZxkV5Jdu3fvnvvgJEnAkOGS5KSuTpjkWHpXEicCPwP8FL1bWPuqqV2m2Xaw9f2LVZdX1ZqqWrNs2bKZui5JGtKwVy4fSnJrkt9Kcswcz/k64IGq2l1Vfwt8GviHwDHtNhnACuDhtjwJHA/Qth8N7Omv77PPdHVJ0pgMFS5V9Y+AX6X3j/auJH+R5PWzPOeDwKlJjmzPTtYC9wI3AWe3NhuB69ry9rZO2/75qqpW39Bmk50IrAJuBW4DVrXZZ4fRe+i/fZZ9lSTNwtKZm/RU1X1J/g2wC7gMeEULh3dU1acP4ji3JLkW+DKwF7gDuBz4LHB1kve02hVtlyuAq5JM0Lti2dCOc0+Sa+gF017ggqp6GiDJW4Ed9Gaiba2qe4btnyRp7oYKlyQ/D5wPvBHYCfzTqvpykp8BvkTv1tbQqmoLsGWf8v30Znrt2/aHwDnTHOdi4OIB9euB6w+mT5Kk7gx75fIfgI/Qu0r5wVSxqh5uVzOSJP3YsOFyBvCDvttOzwIOr6onquqqkfVOkrQgDTtb7HPAEX3rR7aaJEn7GTZcDq+q/ze10paPHE2XJEkL3bDh8v0kJ0+tJHkl8IMDtJckHcKGfebyu8Cnkky9GfFFwD8fTZckSQvdUOFSVbcl+XvAy+h9vMrX2rvrJUnaz9BvogT+AbCy7fOKJFTVlSPplSRpQRv2TZRXAX8X+ArwdCsXYLhIkvYz7JXLGmB1+0wvSZIOaNjZYncDf2eUHZEkLR7DXrkcB9yb5FbgyaliVZ05kl5Jkha0YcPlnaPshCRpcRl2KvJfJ3kxsKqqPpfkSHofZy9J0n6G/ZrjNwPXAh9upeXAZ0bVKUnSwjbsA/0LgFcD34XeF4cBLxhVpyRJC9uw4fJkVT01tdK+y95pyZKkgYYNl79O8g7giCSvBz4F/OfRdUuStJANGy4XAruBu4B/Se8rhP0GSknSQMPOFvsRva85/shouyNJWgyG/WyxBxjwjKWqXtJ5jyRJC97BfLbYlMOBc4Dndd8dSdJiMNQzl6r6dt/rb6rqz4DXjrhvkqQFatjbYif3rT6L3pXMc0fSI0nSgjfsbLE/6Xv9O+CVwK/M9qRJjklybZKvJflqkl9I8rwkO5Pc134e29omyWVJJpLc2R90STa29vcl2dhXf2WSu9o+lyXJbPsqSTp4w84W+6WOz/s+4L9U1dlJDgOOBN4B3FhVlyS5kN7057cDpwOr2utVwAeBVyV5HrCF3lVUAbcn2V5Vj7U2m4Gb6U2bXgfc0PEYJEnTGPa22O8faHtV/emwJ0xyFPCLwJvavk8BTyVZD7ymNdsGfIFeuKwHrmxfVHZzu+p5UWu7s6r2tOPuBNYl+QJwVFV9qdWvBM7CcJGksRn2ttga4C30PrByOfCbwGp6z10O9tnLS+i9IfPPk9yR5KNJfgp4YVU9AtB+Tn122XLgob79J/v6MV19ckBdkjQmB/NlYSdX1fcAkrwT+FRV/cYsz3ky8LaquiXJ++jdApvOoOclNYv6/gdONtO7fcYJJ5xwoD5Lkg7CsFcuJwBP9a0/Bayc5TkngcmquqWtX0svbL7ZbnfRfj7a1/74vv1XAA/PUF8xoL6fqrq8qtZU1Zply5bNcjiSpH0NGy5XAbcmeWeSLcAtwJWzOWFV/V/goSQva6W1wL3AdmBqxtdG4Lq2vB04r80aOxV4vN022wGcluTYNrPsNGBH2/a9JKe2WWLn9R1LkjQGw84WuzjJDcA/bqXzq+qOOZz3bcDH20yx+4Hz6QXdNUk2AQ/S+xQA6M32OgOYAJ5obamqPUneDdzW2r1r6uE+vedDHwOOoPcg34f5kjRGwz5zgd504e9W1Z8nWZbkxKp6YDYnraqv8JMfKTNl7YC2Re/LygYdZyuwdUB9F3DSbPomSZq7Yb/meAu9acEXtdKzgf84qk5Jkha2YZ+5/DPgTOD7AFX1MH78iyRpGsOGy1Pt9lQBtPelSJI00LDhck2SDwPHJHkz8Dn84jBJ0jSGnS32x0leD3wXeBnwb6tq50h7JklasGYMlyRL6L1/5HWAgSJJmtGMt8Wq6mngiSRHj6E/kqRFYNj3ufwQuKt98vD3p4pV9dsj6ZUkaUEbNlw+216SJM3ogOGS5ISqerCqto2rQ5KkhW+mZy6fmVpI8pcj7oskaZGYKVz6vxvlJaPsiCRp8ZgpXGqaZUmSpjXTA/2/n+S79K5gjmjLtPWqqqNG2jtJ0oJ0wHCpqiXj6ogkafEY9rPFJEkamuEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknqnOEiSeqc4SJJ6pzhIknq3LyFS5IlSe5I8ldt/cQktyS5L8knkxzW6s9p6xNt+8q+Y1zU6l9P8oa++rpWm0hy4bjHJkmHuvm8cvkd4Kt96+8FLq2qVcBjwKZW3wQ8VlUvBS5t7UiyGtgAvBxYB3ygBdYS4P3A6cBq4NzWVpI0JvMSLklWAG8EPtrWA7wWuLY12Qac1ZbXt3Xa9rWt/Xrg6qp6sqoeACaAU9proqrur6qngKtbW0nSmMzXlcufAf8a+FFbfz7wnara29YngeVteTnwEEDb/nhr/+P6PvtMV99Pks1JdiXZtXv37rmOSZLUjD1ckvwy8GhV3d5fHtC0Zth2sPX9i1WXV9WaqlqzbNmyA/RaknQwZvqysFF4NXBmkjOAw4Gj6F3JHJNkabs6WQE83NpPAscDk0mWAkcDe/rqU/r3ma4uSRqDsV+5VNVFVbWiqlbSeyD/+ar6VeAm4OzWbCNwXVve3tZp2z9fVdXqG9psshOBVcCtwG3Aqjb77LB2ju1jGJokqZmPK5fpvB24Osl7gDuAK1r9CuCqJBP0rlg2AFTVPUmuAe4F9gIXVNXTAEneCuwAlgBbq+qesY5Ekg5x8xouVfUF4Att+X56M732bfND4Jxp9r8YuHhA/Xrg+g67Kkk6CL5DX5LUOcNFktQ5w0WS1DnDRZLUOcNFktQ5w0WS1DnDRZLUOcNFktQ5w0WS1DnDRZLUOcNFktQ5w0WS1DnDRZLUOcNFktQ5w0WS1DnDRZLUOcNFktQ5w0WS1DnDRZLUOcNFktQ5w0WS1DnDRZLUOcNFktQ5w0WS1Lmxh0uS45PclOSrSe5J8jut/rwkO5Pc134e2+pJclmSiSR3Jjm571gbW/v7kmzsq78yyV1tn8uSZNzjlKRD2XxcuewF/qCqfg44FbggyWrgQuDGqloF3NjWAU4HVrXXZuCD0AsjYAvwKuAUYMtUILU2m/v2WzeGcUmSmrGHS1U9UlVfbsvfA74KLAfWA9tas23AWW15PXBl9dwMHJPkRcAbgJ1VtaeqHgN2AuvatqOq6ktVVcCVfceSJI3BvD5zSbISeAVwC/DCqnoEegEEvKA1Ww481LfbZKsdqD45oD7o/JuT7Eqya/fu3XMdjiSpmbdwSfLTwF8Cv1tV3z1Q0wG1mkV9/2LV5VW1pqrWLFu2bKYuS5KGNC/hkuTZ9ILl41X16Vb+ZrulRfv5aKtPAsf37b4CeHiG+ooBdUnSmMzHbLEAVwBfrao/7du0HZia8bURuK6vfl6bNXYq8Hi7bbYDOC3Jse1B/mnAjrbte0lObec6r+9YkqQxWDoP53w18OvAXUm+0mrvAC4BrkmyCXgQOKdtux44A5gAngDOB6iqPUneDdzW2r2rqva05bcAHwOOAG5oL0nSmIw9XKrqvzP4uQjA2gHtC7hgmmNtBbYOqO8CTppDNyVJc+A79CVJnTNcJEmdM1wkSZ0zXCRJnTNcJEmdM1wkSZ0zXCRJnTNcJEmdM1wkSZ0zXCRJnTNcJEmdM1wkSZ0zXCRJnTNcJEmdM1wkSZ0zXCRJnTNcJEmdM1wkSZ0zXCRJnTNcJEmdM1wkSZ0zXCRJnVs63x2QJI3eygs/O9bzGS4d6P+lfeOSN85jTyTpmcFwkaRFatxXK/0WbbgkWQe8D1gCfLSqLpnnLknSyM1noPRblOGSZAnwfuD1wCRwW5LtVXXvqM893S/W22WSYPrb6MPUF5JFGS7AKcBEVd0PkORqYD0w8nCZzkL9A9KVYf6yTBfAPtNavBbbP6gHa7pxLobxp6rmuw+dS3I2sK6qfqOt/zrwqqp66z7tNgOb2+rLgK/P8pTHAd+a5b4LlWM+NDjmQ8Ncxvziqlq2b3GxXrlkQG2/FK2qy4HL53yyZFdVrZnrcRYSx3xocMyHhlGMebG+iXISOL5vfQXw8Dz1RZIOOYs1XG4DViU5MclhwAZg+zz3SZIOGYvytlhV7U3yVmAHvanIW6vqnhGecs631hYgx3xocMyHhs7HvCgf6EuS5tdivS0mSZpHhoskqXOGy0FIsi7J15NMJLlwwPbnJPlk235LkpXj72W3hhjz7ye5N8mdSW5M8uL56GeXZhpzX7uzk1SSBT1tdZjxJvmV9nu+J8lfjLuPXRviz/UJSW5Kckf7s33GfPSzS0m2Jnk0yd3TbE+Sy9p/kzuTnDynE1aVryFe9CYG/G/gJcBhwP8EVu/T5reAD7XlDcAn57vfYxjzLwFHtuW3HApjbu2eC3wRuBlYM9/9HvHveBVwB3BsW3/BfPd7DGO+HHhLW14NfGO++93BuH8ROBm4e5rtZwA30Huf4KnALXM5n1cuw/vxR8pU1VPA1EfK9FsPbGvL1wJrkwx6Q+dCMeOYq+qmqnqird5M7z1FC9kwv2eAdwP/HvjhODs3AsOM983A+6vqMYCqenTMfezaMGMu4Ki2fDSL4H1yVfVFYM8BmqwHrqyem4FjkrxotuczXIa3HHiob32y1Qa2qaq9wOPA88fSu9EYZsz9NtH7P5+FbMYxJ3kFcHxV/dU4OzYiw/yOfxb42ST/I8nN7RPHF7JhxvxO4NeSTALXA28bT9fm1cH+fT+gRfk+lxEZ5iNlhvrYmQVk6PEk+TVgDfBPRtqj0TvgmJM8C7gUeNO4OjRiw/yOl9K7NfYaelem/y3JSVX1nRH3bVSGGfO5wMeq6k+S/AJwVRvzj0bfvXnT6b9fXrkMb5iPlPlxmyRL6V1OH+gy9JluqI/RSfI64A+BM6vqyTH1bVRmGvNzgZOALyT5Br1709sX8EP9Yf9cX1dVf1tVD9D7gNdVY+rfKAwz5k3ANQBV9SXgcHof7riYdfqxWYbL8Ib5SJntwMa2fDbw+WpPyhaoGcfcbhF9mF6wLPR78TDDmKvq8ao6rqpWVtVKes+ZzqyqXfPT3Tkb5s/1Z+hN3CDJcfRuk90/1l52a5gxPwisBUjyc/TCZfdYezl+24Hz2qyxU4HHq+qR2R7M22JDqmk+UibJu4BdVbUduILe5fMEvSuWDfPX47kbcsx/BPw08Kk2d+HBqjpz3jo9R0OOedEYcrw7gNOS3As8Dfyrqvr2/PV6boYc8x8AH0nye/RuDb1pgf+PIkk+Qe/W5nHtWdIW4NkAVfUhes+WzgAmgCeA8+d0vgX+30uS9AzkbTFJUucMF0lS5wwXSVLnDBdJUucMF0lS5wwXSVLnDBdJUuf+P79fag9fCjUBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot = results.intersect_B.plot.hist(bins = 100)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD4CAYAAAAD6PrjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAARr0lEQVR4nO3de4xk513m8e8Tm1wcwgA7RiDHk3EYx2QULYq3CSBuBgKaJEwMIQR7QbtBXs+SxSAuQjEXQQAhhWsgiiVnApaJUS5OADNDJvISluCwclhPkl3wBWsHY+LBCE+4TIAEjJMff1TNoRi6p09P91unTtX3I7VU51TV6d/bl3rO+77nkqpCkiSAJw1dgCRpcRgKkqSOoSBJ6hgKkqSOoSBJ6lw4dAHbsXv37tq7d+/QZUjSqHzgAx/4SFVdvN5zow6FvXv3cvz48aHLkKRRSfJnGz3n8JEkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6hoIkqWMoSJI6oz55TZJWyd4b39U9fvi1L2nyPewpSJI6hoIkqbMwoZDkqiTvS3JzkquGrkeSVlHTUEhyS5LHktx71voDSR5MciLJjdPVBfw98FTgZMu6JEnra91TuBU4MLsiyQXATcCLgP3AtUn2A++rqhcBrwZ+rHFdkqR1NA2FqroL+OuzVr8AOFFVD1XV48DbgKur6pPT5/8GeMpG20xyKMnxJMdPnTrVpG5JWlVDzClcAjwys3wSuCTJy5K8EbgNeMNGb66qw1W1VlVrF1+87j0iJEnnaYjzFLLOuqqqXwd+fd7FSJL+1RA9hZPApTPLzwQe3coGkhxMcvj06dM7WpgkrbohQuEe4PIklyV5MnANcGQrG6iqo1V1aNeuXU0KlKRV1fqQ1LcCdwNXJDmZ5LqqegK4AbgTeAC4varua1mHJKmfpnMKVXXtBuuPAcfOd7tJDgIH9+3bd76bkCStY2HOaN4Kh48kqY1RhoIkqQ1DQZLUGWUoeEiqJLUxylBwTkGS2hhlKEiS2jAUJEmdUYaCcwqS1MYoQ8E5BUlqY5ShIElqw1CQJHUMBUlSx1CQJHVGGQoefSRJbYwyFDz6SJLaGGUoSJLaMBQkSR1DQZLUMRQkSZ1RhoJHH0lSG6MMBY8+kqQ2RhkKkqQ2DAVJUsdQkCR1DAVJUsdQkCR1DAVJUmeUoeB5CpLUxihDwfMUJKmNUYaCJKkNQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEmdUYaCl7mQpDZGGQpe5kKS2hhlKEiS2jAUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEmdhQqFJE9P8oEkXzd0LZK0ipqGQpJbkjyW5N6z1h9I8mCSE0lunHnq1cDtLWuSJG2sdU/hVuDA7IokFwA3AS8C9gPXJtmf5IXA/cBfNq5JkrSBC1tuvKruSrL3rNUvAE5U1UMASd4GXA18KvB0JkHx8STHquqTZ28zySHgEMCePXvaFS9JK6hpKGzgEuCRmeWTwBdW1Q0ASV4JfGS9QACoqsPAYYC1tbVqW6okrZYhQiHrrOs+3Kvq1vmVIkmaNcTRRyeBS2eWnwk8upUNJDmY5PDp06d3tDBJWnVDhMI9wOVJLkvyZOAa4MhWNlBVR6vq0K5du5oUKEmrqlcoJHne+Ww8yVuBu4ErkpxMcl1VPQHcANwJPADcXlX3nc/2JUk7q++cws3TvfpbgbdU1d/2eVNVXbvB+mPAsZ7fW5I0J716ClX1pcC3MJkLOJ7kLUm+pmll5+CcgiS10XtOoar+P/DDTM46/grg9Un+OMnLWhV3jlqcU5CkBvrOKfzHJK9jMgfwVcDBqnru9PHrGtYnSZqjvj2FNwAfBD6/qr6jqj4IUFWPMuk9zJXDR5LURt9QeDGTCeaPAyR5UpKLAKrqtlbFbcThI0lqo28ovAd42szyRdN1kqQl0jcUnlpVf39mYfr4ojYlSZKG0jcU/iHJlWcWkvwn4ONtSpIkDaXvyWvfDbwjyZlrFH0O8M1tStpckoPAwX379g1VgiQtpV6hUFX3JPk84AomVzn946r656aVnbueo8DRtbW164eqQZKW0VYunf0FwN7pe56fhKp6c5OqJEmD6BUKSW4DPhf4v8AnpqsLMBQkaYn07SmsAfuryjudSdIS63v00b3AZ7csZCs8o1mS2ugbCruB+5PcmeTIma+WhZ2LZzRLUht9h49e07IISdJi6HtI6u8leRZweVW9Z3rdowvaliZJmre+l86+Hngn8MbpqkuAO1oVJUkaRt85he8AvgT4KHQ33PmsVkVJkobRNxT+qaoeP7OQ5EIm5ykMwqOPJKmNvqHwe0l+EHja9N7M7wCOtivr3Dz6SJLa6BsKNwKngD8C/jtwjAHuuCZJaqvv0UefBN40/ZIkLam+1z76U9aZQ6iqZ+94RZKkwWzl2kdnPBX4JuAzd74cSdKsvTe+a67fr9ecQlX91czXn1fVLwBf1bg2SdKc9R0+unJm8UlMeg7PaFLRnMym78OvfcmAlUjS4ug7fPRzM4+fAB4GXrHj1fTk7TglqY2+Rx99ZetCtsLbcUpSG32Hj773XM9X1c/vTDmSpCFt5eijLwDO3EPhIHAX8EiLoiRJw+gbCruBK6vq7wCSvAZ4R1X9t1aFSZLmr+9lLvYAj88sPw7s3fFqJEmD6ttTuA34P0l+g8mZzd8AvLlZVZKkQfQ9+ugnk7wb+LLpqm+rqg+1K0uSNIS+w0cAFwEfrapfBE4muaxRTZKkgfS9HeePAq8GfmC66lOAX21VlCRpGH3nFL4BeD7wQYCqejTJqC9zIUmLaN4XwDtb3+Gjx6uqmF4+O8nT25W0OW/HKUlt9A2F25O8Efj0JNcD72HAG+54O05JaqPv0Uc/O70380eBK4AfqarfblqZJGnuNg2FJBcAd1bVCwGDQJKW2KbDR1X1CeBjSRyrkaQl1/foo38E/ijJbwP/cGZlVX1Xk6okSYPoGwrvmn5JkpbYOUMhyZ6q+nBV/cq8CpIkDWeznsIdwJUASX6tqr6xfUmStFqGPmFt1mYTzZl5/OyWhUiShrdZKNQGjyVJS2iz4aPPT/JRJj2Gp00fM12uqvq0ptVJkubqnKFQVRfMqxBJ0vD6HpK61GYneR5+7UsGrESShrWVm+xIkpacPQVJGsAiHYY6a2F6Ckmem+TmJO9M8qqh65GkVdQ0FJLckuSxJPeetf5AkgeTnEhyI0BVPVBV3w68AlhrWZckaX2tewq3AgdmV0wvxX0T8CJgP3Btkv3T514K/D7wO43rkiSto+mcQlXdlWTvWatfAJyoqocAkrwNuBq4v6qOAEeSvAt4S8vaJGneFnUeYdYQE82XAI/MLJ8EvjDJVcDLgKcAxzZ6c5JDwCGAPXv2tKtSklbQEKGQddZVVb0XeO9mb66qw8BhgLW1tR2/9IbnLEhaZUMcfXQSuHRm+ZnAowPUIUk6yxA9hXuAy5NcBvw5cA3wn7eygSQHgYP79u1rUJ4k7ZwxzCPMan1I6luBu4ErkpxMcl1VPQHcANwJPADcXlX3bWW7VXW0qg7t2uVtoyVpJ7U++ujaDdYf4xyTyZI0ZmPrHcwa5WUuHD6SNE+rdADKKEOhqo4CR9fW1q4fuhZJq2sZw2KUoSBJrY15CGg7DIVzWMa9AEnbs+xhMcpQcE5B0qJZlrBYmEtnb4WHpEpSG6MMBUlSG4aCJKnjnEJPTjpLWgWjDAXPU5DUwrJMFm+Hw0eSpM4oewpDO3tvwuEkabzsHfxb9hQkSZ1RhkKSg0kOnz59euhSJGmpjDIUPHlNktpwTkHSynEeYWOj7ClIktqwp7ADPLFNWnz2DvqxpyBJ6oyypzD2S2dvtMdiL0PaWfYOtm6UoeBlLiRtxCDYnlGGwiJzfkGaP4Ng5zinIEnq2FNoaKt7L/YyJA3NnoIkqWNPYUHZa5A0BHsKkqTOKHsKYz9PoTV7GZLO1yhDYdXOU+jzIe8heVo1/s23McpQkLSaDIL2DIWRcWhIUkuGwoi51yRppxkKkhaOPeLhGAqSFsJGPV97xPNlKKwQL9ktaTOGwpJzL0vSVhgKkubGuYLF52UuJEkdewpqxr1CaXxGGQpe+0gajpddWW6jDIVVu/bRPG31H94egM7Yzk2ltDhGGQraWX3+Of0HllaDoaBBDdXrsLezPn8uMhSkkVqED3B7kMvHUNCG5j2stAgfcqvID3bNMhQ0F37g/6ux/CwMi9XkyWuSpI49Bc3dPPdAt7pXPpa9eKkVQ0Gjsp0rvS7zcMh2wk+aZShoIS3yHnufYFrk+qVzMRSkHrazZ70oe+WLUocWm6GglTemD0vPPldrhoIW3lY/CMc0XDPWurW8DAVpA0Pucbu3r6F4noIkqbNQoZDk65O8KclvJvnaoeuRpFXTfPgoyS3A1wGPVdXzZtYfAH4RuAD4pap6bVXdAdyR5DOAnwX+Z+v6pJ2y3SEfh4y0COYxp3Ar8AbgzWdWJLkAuAn4GuAkcE+SI1V1//QlPzx9XtoyP1yl89d8+Kiq7gL++qzVLwBOVNVDVfU48Dbg6kz8FPDuqvrgettLcijJ8STHT5061bZ4SVoxQ80pXAI8MrN8crruO4EXAi9P8u3rvbGqDlfVWlWtXXzxxe0rlaQVMtQhqVlnXVXV64HXz7sYSdLEUD2Fk8ClM8vPBB7t++YkB5McPn369I4XJkmrbKhQuAe4PMllSZ4MXAMc6fvmqjpaVYd27drVrEBJWkXNQyHJW4G7gSuSnExyXVU9AdwA3Ak8ANxeVfe1rkWSdG7N5xSq6toN1h8Djp3PNpMcBA7u27dvO6VJks6yUGc09+XwkSS1McpQkCS1kaoauobzluQU8Gfn+fbdwEd2sJwxsM2rwTavhu20+VlVte6JXqMOhe1Icryq1oauY55s82qwzauhVZsdPpIkdQwFSVJnlUPh8NAFDMA2rwbbvBqatHll5xQkSf/eKvcUJElnMRQkSZ2lD4UkB5I8mOREkhvXef4pSd4+ff4Pkuydf5U7q0ebvzfJ/Un+MMnvJHnWEHXupM3aPPO6lyepJKM/fLFPm5O8Yvq7vi/JW+Zd407r8be9J8nvJvnQ9O/7xUPUuVOS3JLksST3bvB8krx++vP4wyRXbvubVtXSfjG5//OfAM8Gngz8P2D/Wa/5H8DN08fXAG8fuu45tPkrgYumj1+1Cm2evu4ZwF3A+4G1oeuew+/5cuBDwGdMlz9r6Lrn0ObDwKumj/cDDw9d9zbb/OXAlcC9Gzz/YuDdTO5R80XAH2z3ey57T2Hd236e9ZqrgV+ZPn4n8NVJ1rsJ0Fhs2uaq+t2q+th08f1M7mcxZn1+zwA/Afw08I/zLK6RPm2+Hripqv4GoKoem3ONO61Pmwv4tOnjXWzhPi2LqNa/nfGsq4E318T7gU9P8jnb+Z7LHgob3fZz3dfU5JLep4H/MJfq2ujT5lnXMdnTGLNN25zk+cClVfVb8yysoT6/5+cAz0nyv5O8P8mBuVXXRp82vwb41iQnmVyF+TvnU9pgtvr/vqmhbsc5L+ve9vM8XjMmvduT5FuBNeArmlbU3jnbnORJwOuAV86roDno83u+kMkQ0lVMeoPvS/K8qvrbxrW10qfN1wK3VtXPJfli4LZpmz/ZvrxB7Pjn17L3FPrc9rN7TZILmXQ5z9VdW3S9bnWa5IXADwEvrap/mlNtrWzW5mcAzwPem+RhJmOvR0Y+2dz3b/s3q+qfq+pPgQeZhMRY9WnzdcDtAFV1N/BUJheOW1bburXxepY9FPrc9vMI8F+nj18O/K+azuCM1KZtng6lvJFJIIx9nBk2aXNVna6q3VW1t6r2MplHeWlVHR+m3B3R52/7DiYHFZBkN5PhpIfmWuXO6tPmDwNfDZDkuUxC4dRcq5yvI8B/mR6F9EXA6ar6i+1scKmHj6rqiSRnbvt5AXBLVd2X5MeB41V1BPhlJl3ME0x6CNcMV/H29WzzzwCfCrxjOqf+4ap66WBFb1PPNi+Vnm2+E/jaJPcDnwC+v6r+ariqt6dnm78PeFOS72EyjPLKMe/kTW9nfBWwezpP8qPApwBU1c1M5k1eDJwAPgZ827a/54h/XpKkHbbsw0eSpC0wFCRJHUNBktQxFCRJHUNBktQxFCRJHUNBktT5F5yVmpNQJuu3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot = results.intersect_B.plot.hist(bins = 100)\n",
    "plot.set_yscale('log')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On va regarder le pourcentage de hits pour nos données et la ref dans le bedtools intersect. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "ID_A = pa.DataFrame(list(range(1,len(data)+1)))\n",
    "ID_B = pa.DataFrame(list(range(1,len(ref)+1)))\n",
    "ID_A.columns = ['ID']\n",
    "ID_B.columns = ['ID']\n",
    "\n",
    "A_presence = ID_A[ID_A.ID.isin(results.A_ID)]\n",
    "B_presence = ID_B[ID_B.ID.isin(results.B_ID)]\n",
    "A_absence = ID_A[~ID_A.ID.isin(results.A_ID)]\n",
    "B_absence = ID_B[~ID_B.ID.isin(results.B_ID)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0              5\n",
       "1              6\n",
       "2              9\n",
       "3             10\n",
       "4             12\n",
       "           ...  \n",
       "163495    163549\n",
       "163496    163550\n",
       "163497    163551\n",
       "163498    163553\n",
       "163499    163553\n",
       "Name: A_ID, Length: 163500, dtype: int64"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "on va juste regarder rapidement comment ça change lorsque l'on demande a bedtools d'etre plus stringeant sur ses intersections. (on accepte que les hits a 80% reciproques)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "149210 VS 132561\n",
      "147892 VS 132519\n",
      "14349 VS 30998\n",
      "15536 VS 30909\n"
     ]
    }
   ],
   "source": [
    "results_best = pa.read_csv(\"/home/tkastylevsky/results/results_bedtools/chk_galGal6_chr1_libwholegenome_compared_to_chr1_libchr1_80.txt\",header = None, sep = '\\s+')\n",
    "results_best = results_best.drop([0,1,2,5,6,7,9,10,11,14,15,16], axis = 1)\n",
    "results_best.columns = names_col_2\n",
    "\n",
    "A_presence_best = ID_A[ID_A.ID.isin(results_best.A_ID)]\n",
    "B_presence_best = ID_B[ID_B.ID.isin(results_best.B_ID)]\n",
    "A_absence_best = ID_A[~ID_A.ID.isin(results_best.A_ID)]\n",
    "B_absence_best = ID_B[~ID_B.ID.isin(results_best.B_ID)]\n",
    "print(len(A_presence), 'VS', len(A_presence_best))\n",
    "print(len(B_presence), 'VS', len(B_presence_best))\n",
    "print(len(A_absence), 'VS', len(A_absence_best))\n",
    "print(len(B_absence), 'VS', len(B_absence_best))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.912270189962032 VS 0.8104781760710202\n",
      "0.9049367305480089 VS 0.810870842205742\n"
     ]
    }
   ],
   "source": [
    "#des petits pourcentages sur les totaux: \n",
    "perc_hit_A = len(A_presence)/(len(A_presence)+len(A_absence))\n",
    "perc_hit_B = len(B_presence)/(len(B_presence)+len(B_absence))\n",
    "\n",
    "perc_hit_A_best = len(A_presence_best)/(len(A_presence_best)+len(A_absence_best))\n",
    "perc_hit_B_best = len(B_presence_best)/(len(B_presence_best)+len(B_absence_best))\n",
    "\n",
    "print(perc_hit_A, 'VS', perc_hit_A_best )\n",
    "print(perc_hit_B, 'VS', perc_hit_B_best )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "on s'attaque au type matching ! yay ! \n",
    "\n",
    "edit : \n",
    "\n",
    "lundi 27 janvier,  il y a un soucis dans le results_type, qui vient sans doute du merge et de la façon dont il est fait (il est beaucoup trop petit)\n",
    "\n",
    "\n",
    "edit : c'est corrigé. il fallait faire un \"inner\" et pas un \"outer\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_type = results.drop(['intersect'], axis = 1)\n",
    "results_type = pa.merge(results_type, data, on = 'A_ID', how = 'inner')\n",
    "results_type = results_type[['A_ID','B_ID','intersect_A','intersect_B', 'matching_repeat','repeat_class','repeat_family','TE_A_ID']]\n",
    "results_type.columns = ['A_ID','B_ID','intersect_A','intersect_B', 'matching_repeat_A', 'repeat_class_A','repeat_family_A','TE_A_ID']\n",
    "ref_type = ref[['B_ID','matching_repeat', 'repeat_class','repeat_family','TE_B_ID']]\n",
    "ref_type.columns = ['B_ID','matching_repeat_B','repeat_class_B','repeat_family_B','TE_B_ID']\n",
    "\n",
    "\n",
    "results_type = pa.merge(results_type, ref_type, on = 'B_ID', how = 'inner')\n",
    "\n",
    "results_type.to_csv('/home/tkastylevsky/results/bedtools_analysis/chk_galgal6_chr1_libwhole_vs_chr1_libchr1_repeat_fam_intersect.csv',sep = '\\t', index = False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour la suite, voir le code R bedtools display, pour se faire une idée générale de la qualité de l'annotation (pas bonne). c'est trèèèèèèèèèèèsssssss long d'essayer de le faire sur Python. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "ans = results_type.groupby('A_ID').apply(lambda grp: grp.nlargest(1,'intersect_A'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['A_ID', 'B_ID', 'intersect_A', 'intersect_B', 'matching_repeat_A',\n",
       "       'repeat_class_A', 'repeat_family_A', 'TE_A_ID', 'matching_repeat_B',\n",
       "       'repeat_class_B', 'repeat_family_B', 'TE_B_ID'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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
}