results/data_analysis/Table_1_and_2_reviewed.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import stats as st\n",
    "import seaborn as sns\n",
    "\n",
    "from IPython.core.pylabtools import figsize\n",
    "\n",
    "import numpy.random as r\n",
    "from pylab import *\n",
    "from matplotlib.gridspec import GridSpec\n",
    "\n",
    "import sys\n",
    "sys.path.insert(0, '../../utils')\n",
    "import splicing_utils as spu\n",
    "import single_cell_plots as scp\n",
    "from single_cell_plots import *\n",
    "\n",
    "plt.rcParams[\"axes.edgecolor\"] = \"black\"\n",
    "plt.rcParams[\"axes.linewidth\"] = 1\n",
    "plt.rcParams[\"axes.facecolor\"] = 'white'\n",
    "\n",
    "import matplotlib as mpl\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "mpl.rcParams[\"mathtext.fontset\"] = \"stix\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = '/mnt/c/Users/ferna/Desktop/SingleCell/data/'\n",
    "%run -i '../../utils/load_data.py'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "from scipy.stats import spearmanr\n",
    "import rpy2\n",
    "import rpy2.robjects.packages as rpackages\n",
    "import rpy2.robjects as robjects\n",
    "import rpy2.robjects.numpy2ri as rpyn\n",
    "from statsmodels.stats.multitest import multipletests\n",
    "dt = rpy2.robjects.packages.importr('diptest')\n",
    "\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.cluster import AgglomerativeClustering\n",
    "\n",
    "from scipy.stats import hypergeom\n",
    "def hyper_test(M, n, N, k):\n",
    "    \n",
    "    hpd = hypergeom(M, n, N)\n",
    "    p_depleted = hpd.cdf(k)\n",
    "    p_enriched = hpd.sf(k-1)\n",
    "    \n",
    "    return p_depleted, p_enriched"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpl.rcParams[\"mathtext.fontset\"] = \"stix\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table 1: Number of skipped exons events per dataset\n",
    "\n",
    "We consider an alternative exon as observed if the following two premises are true:\n",
    "* At least one informative junction read is observed in at least 10% of the cells.\n",
    "* The observations amount to an average $\\hat{\\Psi}$ between 0.05 and 0.95."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chen dataset\n",
      "Total observed exons:\n",
      "3276\n",
      "Mean reads per event\n",
      "3.1\n"
     ]
    }
   ],
   "source": [
    "print('Chen dataset')\n",
    "print('Total observed exons:')\n",
    "print(np.sum((chen_PSI.mean(axis=1) <= 0.95) & (chen_PSI.mean(axis=1) >= 0.05) & (chen_PSI.isna().mean(axis=1) <= 0.9)))\n",
    "print('Mean reads per event')\n",
    "print(round(chen_read_counts.loc[(chen_PSI.mean(axis=1) <= 0.95) & (chen_PSI.mean(axis=1) >= 0.05) & (chen_PSI.isna().mean(axis=1) <= 0.9)].mean(axis=1).mean(),1) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lescroart dataset\n",
      "Total observed exons:\n",
      "3007\n",
      "Mean reads per event\n",
      "3.2\n"
     ]
    }
   ],
   "source": [
    "print('Lescroart dataset')\n",
    "print('Total observed exons:')\n",
    "print(np.sum((lescroart_PSI.mean(axis=1) <= 0.95) & (lescroart_PSI.mean(axis=1) >= 0.05) & (lescroart_PSI.isna().mean(axis=1) <= 0.9)))\n",
    "print('Mean reads per event')\n",
    "print(round(lescroart_read_counts.loc[(lescroart_PSI.mean(axis=1) <= 0.95) & (lescroart_PSI.mean(axis=1) >= 0.05) & (lescroart_PSI.isna().mean(axis=1) <= 0.9)].mean(axis=1).mean(), 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trapnell dataset\n",
      "Total observed exons:\n",
      "4457\n",
      "Mean reads per event\n",
      "14.4\n"
     ]
    }
   ],
   "source": [
    "print('Trapnell dataset')\n",
    "print('Total observed exons:')\n",
    "print(np.sum((trapnell_PSI.mean(axis=1) <= 0.95) & (trapnell_PSI.mean(axis=1) >= 0.05) & (trapnell_PSI.isna().mean(axis=1) <= 0.9)))\n",
    "print('Mean reads per event')\n",
    "print(round(trapnell_read_counts.loc[(trapnell_PSI.mean(axis=1) <= 0.95) & (trapnell_PSI.mean(axis=1) >= 0.05) & (trapnell_PSI.isna().mean(axis=1) <= 0.9)].mean(axis=1).mean(), 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Song dataset\n",
      "Total observed exons:\n",
      "5355\n",
      "Mean reads per event\n",
      "88.8\n"
     ]
    }
   ],
   "source": [
    "print('Song dataset')\n",
    "print('Total observed exons:')\n",
    "print(np.sum((song_PSI.mean(axis=1) <= 0.95) & (song_PSI.mean(axis=1) >= 0.05) & (song_PSI.isna().mean(axis=1) <= 0.9)))\n",
    "print('Mean reads per event')\n",
    "print(round(song_read_counts.loc[(song_PSI.mean(axis=1) <= 0.95) & (song_PSI.mean(axis=1) >= 0.05) & (song_PSI.isna().mean(axis=1) <= 0.9)].mean(axis=1).mean(), 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fletcher dataset\n",
      "Total observed exons:\n",
      "684\n",
      "Mean reads per event\n",
      "2.7\n"
     ]
    }
   ],
   "source": [
    "print('Fletcher dataset')\n",
    "print('Total observed exons:')\n",
    "print(np.sum((das_PSI.mean(axis=1) <= 0.95) & (das_PSI.mean(axis=1) >= 0.05) & (das_PSI.isna().mean(axis=1) <= 0.9)))\n",
    "print('Mean reads per event')\n",
    "print(round(das_read_counts.loc[(das_PSI.mean(axis=1) <= 0.95) & (das_PSI.mean(axis=1) >= 0.05) & (das_PSI.isna().mean(axis=1) <= 0.9)].mean(axis=1).mean(), 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shalek dataset\n",
      "Total observed exons:\n",
      "380\n",
      "Mean reads per event\n",
      "213.4\n"
     ]
    }
   ],
   "source": [
    "print('Shalek dataset')\n",
    "print('Total observed exons:')\n",
    "print(np.sum((shalek_PSI.mean(axis=1) <= 0.95) & (shalek_PSI.mean(axis=1) >= 0.05) & (shalek_PSI.isna().mean(axis=1) <= 0.9)))\n",
    "print('Mean reads per event')\n",
    "print(round(shalek_read_counts.loc[(shalek_PSI.mean(axis=1) <= 0.95) & (shalek_PSI.mean(axis=1) >= 0.05) & (shalek_PSI.isna().mean(axis=1) <= 0.9)].mean(axis=1).mean(), 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extent of bimodality among exons\n",
    "\n",
    "In this notebook, we quantify how many events qualify as biomodal. For this, we use the following ad-hoc quartile definition of bimodality:\n",
    "An intermediate exon ($0.2 \\leq \\mu (\\hat{\\Psi}) \\leq 0.8)$ is bimodal if the following are true:\n",
    "\n",
    "1. The first quartile of $\\hat{\\Psi}$ is equal or less than 0.25.\n",
    "\n",
    "2. The third quartile of $\\hat{\\Psi}$ is equal or greater than 0.75.\n",
    "\n",
    "### Shalek bimodal events\n",
    "\n",
    "We check some of the exons reviewed in the Shalek et al., 2013 paper that describes bimodal splicing. None of the exons we checked are bimodal by the quartile definition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "shalek_int_genes, shalek_int_exons = spu.get_int_events(shalek_PSI, shalek_mrna_counts, 0.05)\n",
    "shalek_int_exons = [x for x in shalek_int_exons if x in mrna_per_event_shalek.index]\n",
    "shalek_PSI_filtered, shalek_PSI_mrna_filtered, shalek_good_exons, mrna_filtered, reads_filtered = filter_psi(shalek_PSI, \n",
    "                                            shalek_int_exons, mrna_per_event_shalek, shalek_coverage_tab['SJ_coverage'], \n",
    "                                                                        shalek_read_counts, 10,\n",
    "                                                                       cell_min=0.5)\n",
    "\n",
    "good_cells = shalek_PSI_filtered.columns[shalek_PSI_filtered.isna().mean() <= 0.5]\n",
    "shalek_PSI_good = shalek_PSI_filtered[good_cells]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "shalek_paper_bin_exon = ['Acpp_AE', 'Clec7a_other_1', 'Irgm1_1', 'Irf7_1', 'Clec4n_2', 'Sat1_nmdSE_1', 'Zfp207', \n",
    "                         'Abi1_7', 'Srsf7_nmdSE_1', 'Psmg4_1']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Clec7a_other_1\n",
      "0.0\n",
      "0.13411938098747234\n",
      "Irgm1_1\n",
      "0.3618864292589028\n",
      "0.5\n",
      "Clec4n_2\n",
      "0.5853566009104704\n",
      "0.6755955139493929\n",
      "Sat1_nmdSE_1\n",
      "0.0014144271570014145\n",
      "0.11048371917937136\n",
      "Abi1_7\n",
      "0.0995260663507109\n",
      "0.6094674556213018\n",
      "Srsf7_nmdSE_1\n",
      "0.0\n",
      "0.1562189054726368\n"
     ]
    }
   ],
   "source": [
    "for event in [x for x in shalek_paper_bin_exon if x in shalek_PSI_filtered.index]:\n",
    "    print(event)\n",
    "    print(shalek_PSI.loc[event, mrna_per_event_shalek.columns].quantile(0.25))\n",
    "    print(shalek_PSI.loc[event, mrna_per_event_shalek.columns].quantile(0.75))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is using the MISO calculations of $\\hat{\\Psi}$ from the Shalek paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Irgm1_1\n",
      "0.19\n",
      "0.98\n",
      "Irf7_1\n",
      "0.11\n",
      "0.8899999999999999\n",
      "Clec4n_2\n",
      "0.515\n",
      "0.665\n",
      "Zfp207\n",
      "0.0\n",
      "0.20500000000000002\n",
      "Abi1_7\n",
      "0.395\n",
      "0.9\n"
     ]
    }
   ],
   "source": [
    "sra_meta = pd.read_csv(data_dir + 'shalek/shalek.meta.tab', sep='\\t', index_col=0)\n",
    "shalek_PSI_paper = pd.read_csv(data_dir + 'shalek/shalek.psi_paper.csv', index_col = 0)\n",
    "shalek_counts_paper = pd.read_csv(data_dir + 'shalek/shalek.expression_paper.csv', index_col = 0)\n",
    "shalek_pca = pd.read_csv(data_dir + 'shalek/shalek.pca.tab', index_col = 0, sep='\\t')\n",
    "shalek_cells = shalek_PSI_paper.columns[1:18]\n",
    "\n",
    "for event in shalek_paper_bin_exon:\n",
    "    gen = event.split('_')[0]\n",
    "    if len(shalek_PSI_paper.loc[shalek_PSI_paper.gene==gen, shalek_cells].index) >= 1:\n",
    "        print(event)\n",
    "        print(shalek_PSI_paper.loc[shalek_PSI_paper.loc[shalek_PSI_paper.gene==gen, shalek_cells].index[0],\n",
    "                                         shalek_cells].quantile(0.25))\n",
    "        print(shalek_PSI_paper.loc[shalek_PSI_paper.loc[shalek_PSI_paper.gene==gen, shalek_cells].index[0],\n",
    "                                         shalek_cells].quantile(0.75))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Song bimodal events\n",
    "\n",
    "Here we compare the modality of observations as determined in the Song et al., 2017 paper, versus the quartile definition of bimodality. We also see the proportion of exons defined as bimodal (either by Song et al., or by the quartile definition that are selected by the 10 mRNA filter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "song_modalities = pd.read_csv(data_dir + 'song/song_event_modalities.tab', sep='\\t', index_col = 0)\n",
    "song_outrigger = pd.read_csv(data_dir + 'song/song_outrigger_psi.tab', index_col = 0)\n",
    "outrigger_iPSC = [x for x in list(song_outrigger.index) if x[0] == 'P']\n",
    "outrigger_NPC = [x for x in list(song_outrigger.index) if ((x[0] == 'C') or (x[0] == 'N'))]\n",
    "outrigger_MN = [x for x in list(song_outrigger.index) if x[0] == 'M']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "song_modalities_original = pd.read_csv('../../../GSE85908_modalities_tidy.csv', index_col=0)\n",
    "song_skipped = [x for x in song_modalities_original.index if len(song_modalities_original.loc[x, \n",
    "                                                                                'event_id'].split('@exon:')) == 3]\n",
    "song_modalities_original = song_modalities_original.loc[song_skipped]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iPSC\n",
      "3223 out of 11071 bimodal.\n",
      "29.11%\n",
      "NPC\n",
      "2245 out of 8222 bimodal.\n",
      "27.3%\n",
      "MN\n",
      "1269 out of 5206 bimodal.\n",
      "24.38%\n"
     ]
    }
   ],
   "source": [
    "ipsc_observations_all = len(song_modalities_original.loc[song_modalities_original.phenotype == 'iPSC'])\n",
    "ipsc_bimodal_all = np.sum(song_modalities_original.loc[song_modalities_original.phenotype == 'iPSC'].modality=='bimodal')\n",
    "\n",
    "npc_observations_all = len(song_modalities_original.loc[song_modalities_original.phenotype == 'NPC'])\n",
    "npc_bimodal_all = np.sum(song_modalities_original.loc[song_modalities_original.phenotype == 'NPC'].modality=='bimodal')\n",
    "\n",
    "mn_observations_all = len(song_modalities_original.loc[song_modalities_original.phenotype == 'MN'])\n",
    "mn_bimodal_all = np.sum(song_modalities_original.loc[song_modalities_original.phenotype == 'MN'].modality=='bimodal')\n",
    "\n",
    "print('iPSC')\n",
    "print(str(ipsc_bimodal_all)+' out of '+str(ipsc_observations_all)+' bimodal.')\n",
    "print(str(round(100*(ipsc_bimodal_all/ipsc_observations_all), 2)) + '%')\n",
    "\n",
    "print('NPC')\n",
    "print(str(npc_bimodal_all)+' out of '+str(npc_observations_all)+' bimodal.')\n",
    "print(str(round(100*(npc_bimodal_all/npc_observations_all), 2)) + '%')\n",
    "\n",
    "print('MN')\n",
    "print(str(mn_bimodal_all)+' out of '+str(mn_observations_all)+' bimodal.')\n",
    "print(str(round(100*(mn_bimodal_all/mn_observations_all), 2)) + '%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_proportions(PSI_tab, subpop, mrna_counts, mrna_per_event, read_counts, coverage_tab, psi_int = 0.1, mrna_min = 10,\n",
    "                   outrigger_tab = '', outrigger_modalities = '', outrigger_subpop = '', \n",
    "                    outrigger_cell_type = '', cell_min=0.5):\n",
    "    \n",
    "    PSI_filtered = process_subpop(subpop, PSI_tab, mrna_counts, \n",
    "                                    mrna_per_event, read_counts, coverage_tab['SJ_coverage'], psi_int, \n",
    "                                  mrna_min, cell_min=cell_min)\n",
    "    \n",
    "    PSI_observed = (PSI_tab[subpop].isna().mean(axis = 1) <= (1-cell_min))\n",
    "    PSI_int = (np.abs(0.5-PSI_tab[subpop].mean(axis = 1)) <= (0.5-psi_int))\n",
    "    PSI_unfiltered = PSI_tab.loc[(PSI_observed & PSI_int), subpop]\n",
    "    \n",
    "    filtered_bimodal_by_song = 0\n",
    "    unfiltered_bimodal_by_song = 0\n",
    "    \n",
    "    assert np.all([x in PSI_unfiltered.index for x in PSI_filtered[0].index])\n",
    "    \n",
    "    suma = len(PSI_unfiltered.index)\n",
    "    suma_f = len(PSI_filtered[0].index)\n",
    "    suma_unfiltered = 0\n",
    "    suma_filtered_all = 0\n",
    "    suma_filtered = 0\n",
    "    suma_outrigger = 0\n",
    "    suma_outrigger_total = 0\n",
    "\n",
    "    for evento in PSI_unfiltered.index:\n",
    "\n",
    "        q25_unfiltered = PSI_unfiltered.loc[evento].quantile(0.25)\n",
    "        q75_unfiltered = PSI_unfiltered.loc[evento].quantile(0.75)\n",
    "        if (q25_unfiltered <= 0.25) and (q75_unfiltered >= 0.75):\n",
    "            suma_unfiltered += 1\n",
    "            \n",
    "        if len(outrigger_modalities) >= 1:\n",
    "            if evento in outrigger_modalities.index:\n",
    "                outrigger_event = outrigger_modalities.loc[evento, 'song_annotation']\n",
    "                if outrigger_modalities.loc[evento, outrigger_cell_type] == 'bimodal':\n",
    "                    unfiltered_bimodal_by_song += 1\n",
    "\n",
    "        if evento in PSI_filtered[0].index:\n",
    "\n",
    "            q25_filtered_all = PSI_tab.loc[evento, subpop].quantile(0.25)\n",
    "            q75_filtered_all = PSI_tab.loc[evento, subpop].quantile(0.75)        \n",
    "            if (q25_filtered_all <= 0.25) and (q75_filtered_all >= 0.75):\n",
    "                suma_filtered_all += 1\n",
    "\n",
    "            q25_filtered = PSI_filtered[0].loc[evento].quantile(0.25)\n",
    "            q75_filtered = PSI_filtered[0].loc[evento].quantile(0.75)\n",
    "            if (q25_filtered <= 0.25) and (q75_filtered >= 0.75):\n",
    "                suma_filtered += 1\n",
    "                \n",
    "                \n",
    "            if len(outrigger_modalities) >= 1:\n",
    "                if evento in outrigger_modalities.index:\n",
    "                    outrigger_event = outrigger_modalities.loc[evento, 'song_annotation']\n",
    "                    if outrigger_modalities.loc[evento, outrigger_cell_type] == 'bimodal':\n",
    "                        filtered_bimodal_by_song += 1\n",
    "                    \n",
    "\n",
    "        if len(outrigger_modalities) >= 1:\n",
    "            if evento in outrigger_modalities.index:\n",
    "                suma_outrigger_total += 1\n",
    "                outrigger_event = outrigger_modalities.loc[evento, 'song_annotation']\n",
    "\n",
    "                q25_outrigger = outrigger_tab.loc[outrigger_subpop, outrigger_event].quantile(0.25)\n",
    "                q75_outrigger = outrigger_tab.loc[outrigger_subpop, outrigger_event].quantile(0.75)\n",
    "\n",
    "                if (q25_outrigger <= 0.25) and (q75_outrigger >= 0.75):\n",
    "                    suma_outrigger += 1\n",
    "\n",
    "\n",
    "    p_deplete = hyper_test(suma, suma_unfiltered, suma_f, suma_filtered)[0]\n",
    "    \n",
    "    \n",
    "    print('Total intermediate exons: ' + str(suma))\n",
    "    print('Total intermediate exons that are bimodal by quartile definition: ' + str(suma_unfiltered))\n",
    "    print('#################')\n",
    "    if len(outrigger_modalities) >= 1:\n",
    "        print('Intermediate exons observed in Song et al.: ' + str(suma_outrigger_total))\n",
    "        print('Song et al. exons that are bimodal by Song et al. definition: ' + str(unfiltered_bimodal_by_song))\n",
    "        print(str(round(100*unfiltered_bimodal_by_song/suma_outrigger_total))+'%')\n",
    "        print('Song et al. exons that are bimodal by quartile definition: ' + str(suma_outrigger))\n",
    "        print(str(round(100*suma_outrigger/suma_outrigger_total))+'%')\n",
    "        print('#################')\n",
    "    print('Selected exons: ' + str(suma_f))\n",
    "    print('Selected exons that are bimodal by Song et al. definition: ' + str(filtered_bimodal_by_song))\n",
    "    print(str(round(100*filtered_bimodal_by_song/suma_f))+'%')\n",
    "    print('Selected exons that are bimodal by quartile definition: ' + str(suma_filtered))\n",
    "    print(str(round(100*suma_filtered/suma_f))+'%')\n",
    "    print('#################')\n",
    "    print('Depletion of quartile bimodality in selected exon (p-value): ' + str(p_deplete))\n",
    "                    \n",
    "    \n",
    "    return (suma, suma_f, suma_unfiltered, suma_filtered_all, suma_filtered, suma_outrigger, suma_outrigger_total, \n",
    "            p_deplete, filtered_bimodal_by_song, unfiltered_bimodal_by_song) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Percent bimodal, from the supplementary data from Song et al., 2017\n",
    "\n",
    "Modalities as reported in **GSE85908_modalities_tidy.csv.gz**; see process_data.ipynb for the code that we used to map the exon IDs. We only use skipped exons that are reported both by Song et al and by us.\n",
    "\n",
    "##### Percent of shared skipped exons that are reported as bimodal in the Song et al. paper\n",
    "\n",
    "(All shared exons)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total exons in iPSC: 8895\n",
      "Bimodal exons in iPSC: 2835\n",
      "Percentage bimodal in iPSC: 31.87%\n"
     ]
    }
   ],
   "source": [
    "print('Total exons in iPSC: ' + str((len(song_modalities) - (song_modalities.iPSC == '-').sum())))\n",
    "print('Bimodal exons in iPSC: ' + str((song_modalities.iPSC == 'bimodal').sum()))\n",
    "print('Percentage bimodal in iPSC: ' + str(round(((song_modalities.iPSC == 'bimodal').sum() / (len(song_modalities) - (song_modalities.iPSC == '-').sum())) * 100, 2))+'%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total exons in NPC: 6539\n",
      "Bimodal exons in NPC: 1981\n",
      "Percentage bimodal in NPC: 30.3%\n"
     ]
    }
   ],
   "source": [
    "print('Total exons in NPC: ' + str((len(song_modalities) - (song_modalities.NPC == '-').sum())))\n",
    "print('Bimodal exons in NPC: ' + str((song_modalities.NPC == 'bimodal').sum()))\n",
    "print('Percentage bimodal in NPC: ' + str(round(((song_modalities.NPC == 'bimodal').sum() / (len(song_modalities) - (song_modalities.NPC == '-').sum())) * 100, 2))+'%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total exons in MN: 4156\n",
      "Bimodal exons in MN: 1131\n",
      "Percentage bimodal in MN: 27.21%\n"
     ]
    }
   ],
   "source": [
    "print('Total exons in MN: ' + str((len(song_modalities) - (song_modalities.MN == '-').sum())))\n",
    "print('Bimodal exons in MN: ' + str((song_modalities.MN == 'bimodal').sum()))\n",
    "print('Percentage bimodal in MN: ' + str(round(((song_modalities.MN == 'bimodal').sum() / (len(song_modalities) - (song_modalities.MN == '-').sum())) * 100, 2))+'%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Percent bimodal in intermediate exons\n",
    "\n",
    "We break it down as: 1) Percent bimodal in intermediate exons according to Song et al. 2) Bimodal according to the quartile definition. 3) Percent bimodal after filtering according to Song et al. 4) Bimodal after filtering according to the quartile definition.\n",
    "\n",
    "Here we do the analysis in the cell type as labeled in the original paper. The reason is that: 1) The modality definitions in Song et al. are assigned to cell types. 2) We were unable to match cells specifically, as their ID is not matched to the accession run ID in the supplementary materials of Song et al.\n",
    "\n",
    "We observe that:\n",
    "* The percent of exons that are labeled as bimodal by Song et al is larger for intermediate exons than for all exons.\n",
    "* The percent of selected intermediate exons that are bimodal decrease after filtering for both definitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total intermediate exons: 616\n",
      "Total intermediate exons that are bimodal by quartile definition: 269\n",
      "#################\n",
      "Intermediate exons observed in Song et al.: 579\n",
      "Song et al. exons that are bimodal by Song et al. definition: 455\n",
      "79%\n",
      "Song et al. exons that are bimodal by quartile definition: 177\n",
      "31%\n",
      "#################\n",
      "Selected exons: 55\n",
      "Selected exons that are bimodal by Song et al. definition: 8\n",
      "15%\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "0%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 2.5476050557883045e-15\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(616, 55, 269, 0, 0, 177, 579, 2.5476050557883045e-15, 8, 455)"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_proportions(song_PSI, song_iPSC, song_mrna_counts, mrna_per_event_song, song_read_counts, song_coverage_tab, \n",
    "                psi_int = 0.2, mrna_min = 10,\n",
    "                   outrigger_tab = song_outrigger, outrigger_modalities = song_modalities, outrigger_subpop = outrigger_iPSC, \n",
    "                outrigger_cell_type = 'iPSC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total intermediate exons: 220\n",
      "Total intermediate exons that are bimodal by quartile definition: 105\n",
      "#################\n",
      "Intermediate exons observed in Song et al.: 203\n",
      "Song et al. exons that are bimodal by Song et al. definition: 176\n",
      "87%\n",
      "Song et al. exons that are bimodal by quartile definition: 71\n",
      "35%\n",
      "#################\n",
      "Selected exons: 28\n",
      "Selected exons that are bimodal by Song et al. definition: 14\n",
      "50%\n",
      "Selected exons that are bimodal by quartile definition: 2\n",
      "7%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 1.2111076400049347e-06\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(220, 28, 105, 5, 2, 71, 203, 1.2111076400049347e-06, 14, 176)"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_proportions(song_PSI, song_NPC, song_mrna_counts, mrna_per_event_song, song_read_counts, song_coverage_tab, \n",
    "                psi_int = 0.2, mrna_min = 10,\n",
    "                   outrigger_tab = song_outrigger, outrigger_modalities = song_modalities, outrigger_subpop = outrigger_NPC, \n",
    "                outrigger_cell_type = 'NPC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total intermediate exons: 188\n",
      "Total intermediate exons that are bimodal by quartile definition: 85\n",
      "#################\n",
      "Intermediate exons observed in Song et al.: 176\n",
      "Song et al. exons that are bimodal by Song et al. definition: 154\n",
      "88%\n",
      "Song et al. exons that are bimodal by quartile definition: 59\n",
      "34%\n",
      "#################\n",
      "Selected exons: 20\n",
      "Selected exons that are bimodal by Song et al. definition: 16\n",
      "80%\n",
      "Selected exons that are bimodal by quartile definition: 4\n",
      "20%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 0.01345020242302973\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(188, 20, 85, 8, 4, 59, 176, 0.01345020242302973, 16, 154)"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_proportions(song_PSI, song_MN, song_mrna_counts, mrna_per_event_song, song_read_counts, song_coverage_tab, \n",
    "                psi_int = 0.2, mrna_min = 10,\n",
    "                   outrigger_tab = song_outrigger, outrigger_modalities = song_modalities, outrigger_subpop = outrigger_MN, \n",
    "                outrigger_cell_type = 'MN')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### The issue with Song et al's definition of bimodality\n",
    "\n",
    "There are some events that are called bimodal, but they certainly look more like unimodal exons with inflation of 0 and 1 values; or at the very least, included/excluded modalities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEoCAYAAACZ5MzqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debgcVZnH8e+PTRSRIMiOBJVVUNAYQLYggmyCCyJRlLBlhEFRdBQRMOKoiBsiDBgFE4ZVYCKRnWEx4BiQfZOdCCFIRCCArCHv/HFOm6bTfW9X9Xpzf5/n6ae6q05XvVV9b79d55w6pYjAzMysjEV6HYCZmQ1dTiJmZlaak4iZmZXmJGJmZqU5iZiZWWlOImZmVpqTiFmXSApJ13Rw/ePyNsZ1ahv9pNH+SpohaUZvohp+nESsIUlrS/qppJslPSXp1Ty9XtKPJb2/znsm5H/skHRQg/VW/vn/c4Btv0nSM7ncmU3E+ilJl0qaneP8h6S7JZ0uae9ie153/VtKOl/S45JeztPLJe3U6ro7TdI1VZ9J5fGcpJskHS7pjXXes52kKZJmSXpF0tOS7pN0rqQvSVIv9sX6z2K9DsD6T/6COCo/FgFuBs4BngKWBt4DfBH4qqSDI+LEBquaIOn0iHi2RBifBpYBAviEpOUi4h8N4p0IHAC8CFwEPAwsBbwD+CgwBphcIobK+o8Avgs8CVwIPA4sD2yc131xk6taD3ihbBxNmAJMJ8VXz2RgBiBgNeATwPeA3SRtERGvAkg6PM+fC1wK3AssDqwJbA3sDvxXXm7DnJOI1XMUMAF4FBgbEX+sLSBpBeDLpC/6eh4A3gV8A/hWiRjGA/OAnwD/AewN/LROHJuTEshMYLOImFmzfHHSF30pkj5FSiD/C3wiIp6rs/6mRMQ9ZeNocv1zgDkDFJkUEddUXuTkeAswGvgMMFnSGsDRwLPAFhFxR/UKJC0CbAe81t7obahydZa9jqR3AEcArwA71ksgABExOyIOB45tsKpfALOAr0harWAMGwCbAlcCP8yxHNCg+OZ5en5tAslxvhoRVxTZflUci+TtvwB8pjaBVNZfYH0LtIlUVf+NkbS7pBskvZCrDc+WtGqB9RdqE4mIx4H/yS9H5+kmwKLA1bUJJL9nXkRcFgXGS5L0Vknfk3Rn3rc5km6TdIykpeqU/YGkv0h6MZe9UtL2zW6vQQxL5Gq4m3PV3Au57eQCSR9uZd3DnZOI1dqHdIZ6XkTcNVjhiGhUpfECcCTwRlLVSBHj83RSrsK6EFhX0pZ1ylaquNYuuI1mfJBUhXMx8LSknSV9Q9IhkjZr87YOAk4nVTedCNxJqtL7X0lvaPO2qlXaNipJoXI83yFp0ZZXLq1Jqg49HHgJOAk4lXTm+BXgbVVl1wBuAg4D/g6cTKpGXQ+4VFKjHxLNmAT8nFQtdxpwPDAN2BDYoYX1DnuuzrJalV/2V7VhXZOAQ4C9JP0sIm4d7A2SlgT2IlXLTKlazydIyeXamrdcmsvuKGkqcDbwZ+CBIr+WG/hAnj5B+iLcsCbWacDuEfH3FrcD6YvsA9W//nOHgrHAbsBv27CN15G0Mum4Alyfp9OBv5L29WpJk/KyeyKiTBXW6cAawOER8YOa7S8PPF81a3IuOzYizq4qNwK4Bjhe0tSIeKJIAJKWAfYkJahNavdD0nJF1mev5zMRq7VSnj5Wu0DSyFz9Uv34cqMVRcQ84Oukv7MfN7n9PYBlgXMi4sU87xLSF/nukpat2cZjwMeBB0mN6GcA9wFzcm+tvVr4Rb1Cnn6BdEb1YVLHgg2Ay4CtgHNLrrvW8XWqj36Vp6NrC5c0Ln9m35F0CnA3aR9vICVfIuKfwK7ArcCWwCmks6LnJP1B0kHNnhkp9d77YF7XD2uXR8STEfFSLvteUqP9+dUJJJd7Bvg2sCTwyeK7TZDOuF4mtbPVxlG3w4Y1x2ciVqu2eqPaSNI/c7W/Asc1WllEXCbpcmB7STtFxGA9mSpVFr+pWsdcSWcAhwKfI1VFVG/jaklrk86itib1mtoc+Eh+7C1pl4h4eZBt16okH5HOOG7Lr++S9HFSstpa0mYR8aeC6651Y515j+bpsnWWlVHd1fmfwP3A+cBPq9t2IuJ2YGNJo4BtgPeT2qi2yo/xkraJiKcH2d6meXpZ/kExkEr14DKSJtRZXqn2Wm+Q9SwgIp6V9HvSj4xbJZ1POqO9PiI62VtuWHASsVqPA+sCCzTo5p49ApC0GNBso/LXSL9Gj5V0WaNCktYDtiBVnUyvWfwbUhI5gJokkmObR/piuDavS6ReRJNJZxAHMkCya6DyJflQVQKpbO/FvC/7kc4UWk0iz9SZV2lvarltItumunfWYCLiRqqSm6TRpOP5XtKPiYZnodmIPF3grLaOSpXSdvnRyJubWFc9nyb1FPwM8J087yVJ5wFfK1pFZvO5OstqVXpjbduuFeZqmsnAu4F9ByhaaVBfVzUXxwGVqp4NJH2wiW1GRFxO6mkG8KESod+bp/W+4GF+klngYr2FUUTcABycXzZzPCvHrZkeZpWuyYdEhAZ47FM0bkhJPyImRMTawNtJ7W7X5el5ZdZpiZOI1ZpE+gW8ez4zaJcjSD22jiZdCPg6uZ79c6Q661NJdfG1j8pZTJFeOpVuuWWusJ5GOhZrSVqizvIN8nRGiXUPVUWOZ+Vs8iO5u3QzZev1wGuriHg0Is4gVXXeD2zhxvXynETsdSLiQeA/gSWASwb41T+iwfxG651FunBwJepXg3ySVKVxWUTsFxH71z6AT5Hq8vfIPW6QtIOkT9S76E/Sm6u2Na1IvDnmJ0ldTJchXYBZve7tSF9Cc0g9xBYKkkbn603qDYWyOKlKCJo4nhFxE/B/wEZV76te33K5N16l6uxa0ugEdc9WJW2odJFrIZLeJmmTOouWInWUmEu6FslKcJuI1XM06ZfmkcAfJd1E6sHzFCl5jCS1M0CxL+djSVVW76qzrFKV9etGb46I5ySdC4wjVUOcSGq/+RnpOo5rSb8s55KG9dg5x3s9cEKBOKsdSroA71uStiIdhzVIPcJeAw7IvYcWFquQ2p9OkHQdqQfXS8DKpG7IK5FGIzi6yfXtReqe+31Jn8zPBawFbE/6/Gbksp8hdS0/RdKXSJ/bM6TP8j2kM7/NgNkF92lVYLqkv5C6aj8KvAXYJe/P8fUuJLXmOInYAvL1FRMknUXq3roN6R98KVJ1xoOki8b+OyJuLrDe5yV9m3QR2b9IWovUq2o28PtBVvMrUhI5gJRETicN0bEdqcF3K1Lj6zOkxvxzgV9HRKlfmhExO/+KPYKUODYlHYOLgB/U6QAw1F1J+qy3J/XKGkVKxM8C95Au2Dux2S/diHhY0vtIXb0/RmpTeYmUOH5CVUKIiJm5W/AXSWemnyV1KvgbKZn9gvltY0XMIHUEGEP6W16e9IPoXtKFjWc3eqMNTq1fj2VmZsOV20TMzKw0JxEzMyvNbSI2bEgaSWpPacZxC1mDuVlHuE3Ehg1JY4Crmyy+ZkTM6Fw0ZguHYZdE8tXPZmZWQETUvcB0WFZnDbfEaWbWijQUXX1uWDczs9KcRMzMrDQnETMzK81JxMzMSnMSMTOz0pxEzMysNCcRMzMrzUnEzMxKcxIxM7PShuUV62ZmvTLysIt6st0Zx+zckfX6TMTMzEpzEjEzs9KcRMzMrDQnETMzK81JxMzMSnMSMTOz0pxEzMysNCcRMzMrzUnEzMxKcxIxM7PSnETMzKw0JxEzMyutpwMwSjoV2AWYHREb5HnnAOvkIiOAZyJiozrvnQE8B7wGzI2IUV0J2szM/qXXo/hOAk4ATqvMiIhPV55L+gkwZ4D3bxMRT3YsOjMzG1BPk0hETJM0st4ySQL2AD7UzZjMzKx5/dwmsiXwRETc32B5AJdLuknS+IFWJGm8pBsl3dj2KM3MhrFeV2cNZCxw1gDLN4+IWZJWAK6QdE9ETKtXMCImAhMBJEX7QzUzG5768kxE0mLAJ4BzGpWJiFl5OhuYAozuTnRmZlbRl0kE+DBwT0TMrLdQ0lKSlq48B7YH7uxifGZmRo+TiKSzgD8B60iaKWm/vGhPaqqyJK0i6eL8ckXgOkm3ATcAF0XEpd2K28zMkl73zhrbYP64OvNmATvl5w8B7+1ocGZmNqh+rc4yM7MhwEnEzMxKcxIxM7PSnETMzKw0JxEzMyvNScTMzEpzEjEzs9KcRMzMrDQnETMzK81JxMzMSnMSMTOz0pxEzMysNCcRMzMrzUnEzMxKcxIxM7PSnETMzKw0JxEzMyvNScTMzErr9T3WT5U0W9KdVfMmSHpM0q35sVOD9+4g6V5JD0g6rHtRm5lZRa/PRCYBO9SZ/7OI2Cg/Lq5dKGlR4ERgR2B9YKyk9TsaqZmZLaCnSSQipgFPlXjraOCBiHgoIl4BzgZ2a2twZmY2qF6fiTRysKTbc3XXsnWWrwo8WvV6Zp5Xl6Txkm6UdGO7AzUzG876MYmcBLwT2Ah4HPhJnTKqMy8arTAiJkbEqIgY1Z4QzcwM+jCJRMQTEfFaRMwDfkWquqo1E1i96vVqwKxuxGdmZvP1XRKRtHLVy48Dd9Yp9mdgLUlrSloC2BOY2o34zMxsvsV6uXFJZwFjgOUlzQS+DYyRtBGpemoG8G+57CrAryNip4iYK+lg4DJgUeDUiLirB7tgZjas9TSJRMTYOrNPaVB2FrBT1euLgQW6/5qZWff0XXWWmZkNHU4iZmZWWqEkIunnvjLczMwqip6JfBG4Q9I0SZ/NPaPMzGyYKppE9gCuBDYHTgNmSfqxpHXaHpmZmfW9QkkkIs6LiO1JV5QfC7wCHArcLekqSXtIWrwDcZqZWR8q1bAeETMi4pvA25l/drI1cBYwU9Ixkt7ZvjDNzKwftdQ7KyLmVp2dbEYaeuRtwNeBeyVdKOn9bYjTzMz6UMtdfCVtLelM4A+kkXT/DhwHXEe6OHC6pE+3uh0zM+s/pa5Yl/RWYG9gPLA2aVTdP5JG4D03Il7N5UYD/wNMAM5pQ7xmZtZHCiURSVuQxrL6JLAk8DzwS+CkiLijtnxE3CDpN8A32hCrmZn1maJnItPy9C7SWcdpEfH8IO95LD/MzGwhU7RN5Bxg64jYMCL+q4kEQkScHBFrlgvPzMz6WaEzkQaj7pqZ2TBVdOyst0naStLSDZa/JS9fvj3hmZlZPytanXUEcCEwr8Hy14DfA99sJSgzMxsaiiaR7YDLI+Kf9Rbm+ZcDH2k1MDMz639Fk8jqwIODlHkolzMzs4Vc0SQSwGDDvy9Buu/5oCSdKmm2pDur5v1I0j2Sbpc0RdKIBu+dIekOSbdKurHpPTAzs7YpmkTuZYCqKknKyx9ocn2TgB1q5l0BbBAR7wHuY+D2lW0iYqOIGNXk9szMrI2KJpHzgHUlnSDpjdUL8usTgHVocoiTiJgGPFUz7/KImJtfTgdWKxijmZl1SdEr1o8HxgIHAh+TNI10NfqqwFbAKsBtpAEY22FfGiekAC6XFMAvI2Jim7ZpZmZNKnqx4YuSxgD/RbqPyJ5Vi+cBZwIHR8SLrQYm6VvAXOCMBkU2j4hZklYArpB0Tz6zqbeu8aTBIs3MrI0Kj+IbEc8An5F0CPABYATwDHBDRDzZjqAk7Q3sAmwbEdEgjll5OlvSFGA088f2qi07EZiY1113fWZmVlypoeABIuLvwMVtjAUASTuQRv3dOiJeaFBmKWCRiHguP98eOLrdsZiZ2cBavilVKySdBfwJWEfSTEn7kRrnlyZVUd0q6eRcdhVJlaS1InCdpNuAG4CLIuLSHuyCmdmwVvhMJN+Qal9S9dGy1L8mJCJi28HW1WBAx1MalJ1FulMiEfEQ8N5mYzYzs84oelOqdYFrSPdR1wBF3e5gZjYMFK3O+jGwAvBD4B3A4hGxSJ1HU1esm5nZ0Fa0OmtLUvvD4Z0IxszMhpaiZyIC7u5EIGZmNvQUTSI3kYY1MTMzK5xEjgZ2yletm5nZMFe0TWR14ALSmFVnkc5MnqlXMCJOazE2MzPrc0WTyCRS910Bn8uP2u68yvOcRMzMFnJFk8g+HYnCzMyGpKKj+E7uVCBmZjb09HTsLDMzG9pKjeIr6W3AJ4H1gKUiYv+q+WsCd7TjniJmZtbfygzAuB/pDodLMr8Rff+8eEXSqLzjaTCQopmZLTwKVWdJ2o50c6f7gI8DJ1Uvj4g7gbuAj7UrQDMz619Fz0S+ATxOumHUs5I2rlPmdmCzliMzM7O+V7RhfRRwYUQ8O0CZmcBK5UMyM7OhomgSWQL45yBlRgCvlQvHzMyGkqJJZAbw/kHKbALcWyoaMzMbUoomkQuALSV9qt5CSfsA7wHObzUwMzPrf0WTyLHAI8BZks4hN6BLOji/ngjcD/yi2RVKOlXSbEl3Vs17q6QrJN2fp8s2eO/eucz9kvYuuC9mZtaiQkkkIp4GtgauAz4FbE+6VuT4/Pr/gG0jYrB2k2qTgB1q5h0GXBkRawFX5tevI+mtwLdJ1WejgW83SjZmZtYZhS82jIhHgDGS3kM6E1kOmANMj4ibSqxvmqSRNbN3A8bk55OBa0jdi6t9BLgiIp4CkHQFKRmdVTQGMzMrp9SwJwARcTvpmpBOWDEiHs/beVzSCnXKrAo8WvV6Zp63AEnjSVfRm5lZG5VOIn1AdebV3tskzYyYSGqvQVLdMmZmVlyhJCLpqCaLRkR8t0Q8FU9IWjmfhawMzK5TZibzq7wAViNVe5mZWZcUPROZMMCyyi/8yqCMrSSRqcDewDF5ekGdMpcB369qTN8e+GYL2zQzs4KKJpFtGswfAXwA+BJwEXBysyvM92ofAywvaSapx9UxwG/ziMGPkHp+IWkU8IWI2D8inpL0XeDPeVVHVxrZzcysOxTRviYCSRsCNwB7RkS9s4eekxTt3GczsyJGHnZRT7Y745idS79XEhFRrx26vXc2jIg7SFVPh7dzvWZm1p86cXvcR4ANOrBeMzPrM51IIpsAvjWumdkwULSL79sHWM/qwAHAFsBvW4zLzMyGgKK9s2bQ4IK+TKQBGL9WNiAzMxs6iiaR06ifROYBT5N6Zl0QES+3GpjZcDMUe+2YFUoiETGuQ3GYmdkQ1ImGdTMzGyacRMzMrLSivbOuKrmdiIhtS77XzMz6VNGG9TF5GjQeir3pIdrNzGxoK1qdtSRphN2HgX2ANYE35um+wEOkYU/eEBGLVD0WbWPMZmbWJ4omkSOBUcCoiJgcEX+NiJfzdBLz73d+ZJvjNDOzPlQ0iXwWOD8inqm3MA/Ffh6wV6uBmZlZ/yuaRFYBXhmkzKvAyuXCMTOzoaRoEpkJ7CZpiXoLJb0B2A14rNXAzMys/xVNIpOBdwFXSdpK0qIAkhaVtDVwJfAOYFJbozQzs75UtIvvMcD7gV2Bq4F5kp4C3kpKSCL13jqmnUGamVl/KnQmEhGvRsTHSA3nVwFzSAlkDuks5LMR8bGImNtKUJLWkXRr1eNZSV+uKTNG0pyqMke1sk0zMyuu6JkIABFxJnBmm2OpXv+9wEaQqspIbSxT6hS9NiJ26VQcZmY2sKEwdta2wIMR8ddeB2JmZq9XKolIeo+kYyRdIOl/q+aPlLSHpGXbFyJ7Amc1WLaZpNskXSLp3QPEO17SjZJubGNcZmbDXuHqLElHA4czPwFVj4u1COkL/8vAL1oNLncl3hX4Zp3FNwNrRMTzknYCfgesVW89ETERmJjX6XG8zMzapNCZiKQ9gSOAK0htFj+oXh4RDwE3kr7422FH4OaIeKJ2QUQ8GxHP5+cXA4tLWr5N2zUzsyYUrc76EvAAsFtE3E79q9f/QoMzghLG0qAqS9JKkpSfjybtyz/atF0zM2tC0eqsDYFJETHQ0CezgBXLh5RIehOwHfBvVfO+ABARJwO7AwdKmgu8COwZEa6qMjProqJJRMC8QcqsCLxULpz5IuIFYLmaeSdXPT8BOKHV7ZiZWXlFk8j9wAcbLczXdGwB3NVKUGY2PIw87KKebHfGMTv3ZLsLo6JtIr8F3ifpqw2Wf5M0tlbHLkQ0M7P+UfRM5DjgU8CxkvYgd++V9GNgS9INq6aTu9OamdnCrVASiYgXJW0D/Jx0g6rKbW8PJbWVnA4c3OrYWWZmNjQUvtgwIuYA4yQdCnyA1Pg9B7ghIv7e5vjMzKyPFUoikj4PPBERl+Vb4V7WmbDMzGwoKNqwfiqwQycCMTOzoadoEvlbifeYmdlCqmhCuBTYRpITiZmZFU4i3wKWBk7xYIdmZla0d9ZZpJ5Ynwf2lDSDVMVVO2ZVRMS2rYdnZmb9rGgSGVP1/A3AOvlRywMhmpkNAwNWZ0n6Uh5mHYCIWKTJx6IDrdfMzBYOg7WJHEdVl15Jr0k6srMhmZnZUDFYEnmJVG1VofwwMzMbNIk8DHxEUvVNptzeYWZmwOBJ5JfA+4BZkl7L8ybkaq2BHh6A0cxsGBiwd1ZEHC9pNrAzsAqwDfAIMKPzoZmZWb8btItvRJwNnA0gaR7wm4g4utOB5e3NAJ4DXgPmRsSomuUiDUu/E/ACMC4ibu5GbGZmVvw6ke8A13QgjoFsExFPNli2I7BWfmwCnJSnZmbWBUVvSvWdTgVS0m7AaRERwHRJIyStHBGP9zowM7PhoPBNqbosgMslBfDLiKi97e6qwKNVr2fmea9LIpLGA+M7GajZUDXysIt6HULXDcd97pR+TyKbR8QsSSsAV0i6JyKmVS2vd83KAl2Qc/KZCJATkpmZtUFfD+keEbPydDYwBRhdU2QmsHrV69WAWd2JzszM+jaJSFpK0tKV58D2wJ01xaYCn1eyKTDH7SFmZt3Tz9VZKwJTUi9eFgPOjIhLJX0BICJOBi4mde99gNTFd58exWpmNiz1bRKJiIeA99aZf3LV8wD+vZtxmZnZfH1bnWVmZv3PScTMzEpzEjEzs9KcRMzMrDQnETMzK81JxMzMSnMSMTOz0pxEzMysNCcRMzMrzUnEzMxKcxIxM7PSnETMzKw0JxEzMyvNScTMzEpzEjEzs9KcRMzMrLS+vSlVPxp52EU92e6MY3buyXbNzAbjMxEzMyutL5OIpNUlXS3pL5LuknRInTJjJM2RdGt+HNWLWM3MhrN+rc6aC3w1Im6WtDRwk6QrIuLumnLXRsQuPYjPzMzo0zORiHg8Im7Oz58D/gKs2tuozMysVl8mkWqSRgIbA9fXWbyZpNskXSLp3V0NzMzM+rY6CwBJbwbOB74cEc/WLL4ZWCMinpe0E/A7YK0G6xkPjO9osGZmw1DfnolIWpyUQM6IiP+pXR4Rz0bE8/n5xcDikpavt66ImBgRoyJiVEeDNjMbZvoyiUgScArwl4j4aYMyK+VySBpN2pd/dC9KMzPr1+qszYHPAXdIujXPOxx4O0BEnAzsDhwoaS7wIrBnREQvgjUzG676MolExHWABilzAnBCdyIyM7N6+rI6y8zMhgYnETMzK81JxMzMSnMSMTOz0pxEzMysNCcRMzMrzUnEzMxKcxIxM7PSnETMzKw0JxEzMyvNScTMzErry7GzzEYedlGvQzCzJvhMxMzMSnMSMTOz0pxEzMysNCcRMzMrzUnEzMxKcxIxM7PSnETMzKy0vk0iknaQdK+kByQdVmf5GySdk5dfL2lk96M0Mxve+jKJSFoUOBHYEVgfGCtp/Zpi+wFPR8S7gJ8BP+xulGZm1pdJBBgNPBARD0XEK8DZwG41ZXYDJufn5wHbSlIXYzQzG/b6ddiTVYFHq17PBDZpVCYi5kqaAywHPFm7MknjgfFVr9sdb0fJ51hm1qJOfY/0axKp9y0fJcqkmRETgYktByXdGBGjWl3PUOJ9XvgNt/0F73M79Wt11kxg9arXqwGzGpWRtBiwDPBUV6IzMzOgf5PIn4G1JK0paQlgT2BqTZmpwN75+e7AVRFR90zEzMw6oy+rs3Ibx8HAZcCiwKkRcZeko4EbI2IqcArw35IeIJ2B7NmF0FquEhuCvM8Lv+G2v+B9bhv5x7uZmZXVr9VZZmY2BDiJmJlZaU4idQy3IVea2N9DJd0t6XZJV0paoxdxttNg+1xVbndJIWnIdwdtZp8l7ZE/67skndntGNutib/tt0u6WtIt+e97p17E2S6STpU0W9KdDZZL0vH5eNwu6X0tbzQi/Kh6kBryHwTeASwB3AasX1PmIODk/HxP4Jxex93h/d0GeFN+fuBQ3t9m9zmXWxqYBkwHRvU67i58zmsBtwDL5tcr9DruLuzzRODA/Hx9YEav425xn7cC3gfc2WD5TsAlpOvsNgWub3WbPhNZ0HAbcmXQ/Y2IqyPihfxyOum6naGsmc8Y4LvAscBL3QyuQ5rZ5wOAEyPiaYCImN3lGNutmX0O4C35+TIseD3akBIR0xj4erndgNMimQ6MkLRyK9t0EllQvSFXVm1UJiLmApUhV4aiZva32n6kXzJD2aD7LGljYPWIuLCbgXVQM5/z2sDakv4oabqkHboWXWc0s88TgL0kzQQuBr7YndB6puj/+6D68jqRHmvrkCtDQNP7ImkvYBSwdUcj6rwB91nSIqSRocd1K6AuaOZzXoxUpTWGdLZ5raQNIuKZDsfWKc3s81hgUkT8RNJmpGvPNoiIeZ0Pryfa/t3lM5EFDbchV5rZXyR9GPgWsGtEvNyl2DplsH1eGtgAuEbSDFLd8dQh3rje7N/1BRHxakQ8DNxLSipDVTP7vB/wW4CI+BOwJLB8V6Lrjab+34twElnQcBtyZdD9zVU7vyQlkKFeTw6D7HNEzImI5SNiZESMJLUD7RoRN/Ym3LZo5u/6d6ROFEhanlS99VBXo2yvZvb5EWBbAEnrkZLI37saZXdNBT6fe2ltCsyJiMdbWaGrs2pE/w650hFN7u+PgDcD5+b+A49ExK49C7pFTe7zQqXJfb4M2F7S3cBrwH9ExD96F3VrmtznrwK/kvQVUrXOuCH8gxBJZ5GqI5fP7TzfBhYHiIiTSe0+OwEPAC8A+7S8zSF8vMzMrMdcnWVmZqU5iZiZWWlOImZmVpqTiJmZleYkYmZmpTmJmJlZaU4iZmZWmpOItZ2kkfkeHJMGmmcLL3/ew4eTiNkQImmUpN9IekjSi0h2x9UAAAdTSURBVJKelXSHpB9Jamk01prtOAlYU5xErO9IWkbStDxqMJL2k/Ro/lKrPB7JQ1ogaa9cfpneRt45eayjH5LGg9oLuAc4njQEzwvA14D7JO3euyhf5zFgPeCbvQ7EOstjZ1nfiYg5kp4HTpN0T0ScIulU4ELSuD8XRsRHASR9ADgNuCQi5vQu6o47Evg6MAPYJSLuql4o6ZPA6cDZkraLiKu7H+J8EfEqKdHZQs5nIsOcpNFK94t/TNLLkh6XdLmkPWrKbSLpPEl/k/RKPjP4paRVOhTabNK9Dz4KkAfFq4yie3NVuV1zuZZGXq2uvpG0rqTfSXpK0j8lXSdp+wHe29SxqdnG2vm4z5Y0T9KYgWIjJZFXSaMJ31VbJiLOB75CGmjwJKV7olTePyZvd0KD9c9QGvK+8noC8HB+uXfNGeC4ZvZloOqwfFZ1iNK93F/Kf3sn5DPQ18VSpnzV+zr6uVjiM5FhTNIBwEmkEVunAvcDK5BuPHUQ+T4LkvYBfgW8nMs9SrrPxP7ARyVtGhGPtDm8ysigq9eZVz1q6Gp15rViTeBPwJ2k4e9XBj4NXCLpMxFxTnXhksfmncD1wH3AGcAbgWcHiGkf0v/qbyPijgHK/ZqUbNYh3Tis7NnINcAI4BDSfcl/V7Xs1pqyRfcF4ETgQNJ9LCYCr5B+DIwmjTj7aovlu/W5GNDzG8v70ZsHsD7pn+8p4N11lq+Wp2uT/mkfAFatKfMhUgKaUjN/JOlLfdJA8waJ7/u5/JSqeRPyvAlV86bmed9r8XhU4gvgRzXLRuVj9TTwlqr5hY5NzTa+XyC2K/N7Dmii7Bm57BFV88bUHrea98wAZgz2GQ5wvBbYl0bvB7bM8+8FRlTNXwKYlpfNKFu+m5+LH+nh6qzh60DSr9vvRv3qkZlV5RYHDomIx2rKXEX6Ev+opKXbHN/v83QzpRsKLSDPH51ftute6HOAo6tnRLoZ1RmkX+cfr1pU9tg8AXynQEwr5+mjA5Z6fZlOVTPWKrovlZu5fS+qbrsbEa9QvxG+aHno3udiuDprONs0Ty8ZpNxmebp1bsSutQKpHn5t4KY2xUZE/EnSd0nVMz+V9PXq5ZKWAo4FVgSOjnRr03a4OSKeqzP/GtIX2sbA5Dyv7LG5LYrdYrhyX+xmquyKlG2HovuycZ5eV2fZdGBui+Whe5+L4SQynI3I08cGLAXL5el/DFLuza2Fs6CIOErS5cA4Uv3+CnnR3sCOwB3A5hHxf23c7BMN5v8tT6u7EZc9Nn+rW6qxx4F1gbc3UbbSRtTSLU8LKLovleO3wHGOiNck1d5JsWh56N7nYrh31nBWqRoY7AK1SrfZZSJCAzz+0IkgI+K6iNg/IjZh/hnA5IjYJM9vZwKBdGZTz0p5Wt2NuOyxKXqWUPkV/uGBCklalNT+AfDHqkXz8rTRj8ZWrq8pui+VhuoFjnOOf7ma2UXLQ/c+F8NJZDibnqc7Nlluyw7G0k/e16B9Z0ye3lI1r1vHZhKpMfjjkt49QLl9SW0h9wLVX5BP5+nqtW+Q9C7mn5VWey1PFy0a7CAqx2+LOss2ZcFEV7Q8DL+/2Z5yEhm+TiLVJx8paf3ahZIq1SInkHom/UzS2nXKLSGpW/+slS80DVgKyH3+/3VdQwHLAEfVrGsU8FnSL9wpVYu6cmwi4iFSb7XFgakNPq+PAT8nffkfFBHzqhbfQ/pFv5ukFare80bSVe/1PE36Zd5MFVoRp+Xpt1Q1wkDuJPH9NpSH/vqbXei5TWSYioi7JR0EnAzcIukC0nUiy5G6tD4HbBMR90jaFzgVuEvSpaR+9IuTvmC2JF3ot24Xwn53zXQglR9I9RpeBzIN2F/SJqQqocp1IosA/xYR/7puoMvHZgKwFHAocJuky4C78rY+CGwCvAiMzT2Q/iUiXpX0c1InhVskTSH9729HuvZiVu3GIuJ5SdcDW0o6I+/Xa8DUiLi97E5ExB8kTQTGk47Z+aQv/I+SkvQs5le/FS6f39NPf7MLv173Mfajtw9ST5bzSVeIv0L6p7wU2L2m3IakapW/ki7geor5F+R9qKbsSFq8TqRmffuS+vxH1eMB4N8HeM8tpF/fyza5jX/FRxrz6QLSr/EXSMnkIwO8t6lj08oxqFrHaFLb0MOkpPF83taPydf2NHifgMOAB/Pn/Aipd9ubqHOdSH7Pu0hdrf9B+qIOYFwz+zLQclJC/grpDOnl/Dd3Iuks8Dng1lbK9+JzGc4P5YNottCQNIL0xfeTiPj6YOXze0aSvpgnR8S4jgVnDUlai3TGcHZEjG13eesMt4nYwmhLUpXHT3sdiC1I0kqqGtsrz3sTcFx+OaWV8tZdbhOxhU5E/B5YstdxWENfBsZKuoZ0PctKwLaka1wuAc5tsbx1kZOImXXbFcB7ge2Bt5I6P9xH6il2XCxYx160vHWR20TMzKw0t4mYmVlpTiJmZlaak4iZmZXmJGJmZqU5iZiZWWlOImZmVpqTiJmZlfb/GIdIbF6DHkcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "song_modalities.loc['GNAS_6'].song_annotation\n",
    "\n",
    "plt.hist(song_outrigger.loc[outrigger_iPSC, \n",
    "                   'exon:chr20:57470667-57470739:+@exon:chr20:57473996-57474040:+@exon:chr20:57478586-57478640:+'].dropna())\n",
    "plt.title('GNAS_6 in iPS cells', fontsize=20)\n",
    "plt.ylabel('frequency', fontsize=20)\n",
    "plt.xlabel('cell $\\Psi$, per Outrigger', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEkCAYAAADeqh2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5gcVZ3/8fcnYZABlHBTYMJNRe5iELkY5KIiFxUigksAFdeFR/2x3ta4gKiI7oJmF5SVXRcRAUUBASMiGJSAKMolGBABo4BZyQTlGi4mQAjf3x+nGjqV6qnumZ6unp7P63nm6Zmq01Xfrpmpb9U5p85RRGBmZjaUCVUHYGZm3c/JwszMSjlZmJlZKScLMzMr5WRhZmalnCzMzKyUk4VZG0kKSdeN4vaPyvZx1Gjto5s0+rySFkhaUE1U45OTxTgn6TWSTpP0W0mPSlqWvd4k6T8kvb7gPSdl/8Ah6SMNtlv7J//SEPteXdLirNz3moj1UEk/lfRgFucjku6S9F1J72/tkxdu/02SLpX0gKRnsterJR0w0m2PNknX1f1Oal9PSrpV0gmS+gves4+kH0paJOlZSY9J+qOkH0j6qCRV8VmsO61SdQBWjexE8LnsawLwW+Ai4FHgpcBrgX8G/kXSsRFxZoNNnSTpuxHxxDDC+AdgLSCAgyWtGxGPNIj3LOBoYCnwE+DPwBrAK4F3AnsB5w0jhtr2TwS+CDwMXAE8AKwHTMm2fWWTm9oaWDLcOJrwQ+BGUnxFzgMWAAImAwcD/wYcJGn3iFgGIOmEbPlzwE+B+UAfsDmwJ3AI8N/ZejMni3Hsc8BJwP3A9Ii4IV9A0suBj5NO6EXuAV4N/CvwmWHEcAzwPPCfwAzg/cBpBXFMJSWKhcBuEbEwt76PdEIfFkmHkhLFz4GDI+LJgu03JSL+MNw4mtz+48DjQxQ5NyKuq/2QJcF5wM7A4cB5kjYFTgaeAHaPiDvqNyBpArAPsLy90dtY5mqocUjSK4ETgWeB/YsSBUBEPBgRJwBfabCp/wIWAZ+QNLnFGLYDdgWuAb6cxXJ0g+JTs9dL84kii3NZRPyslf3XxTEh2/8S4PB8oqhtv4XtrdRmUVdtt5ekQyTdLGlJVt13oaSBFrbfUptFRDwAXJb9uHP2ugswEbg2nyiy9zwfEbOjhbGAJK0j6d8k/T77bI9Lul3SqZLWKCh7iqS7JS3Nyl4j6W3N7q9BDKtm1We/zarUlmRtGz+S9NaRbNucLMarD5DuKi+JiDvLCkdEo6qIJcBngX5SlUYrjslez82qnq4AtpL0poKytaqp17S4j2a8kVT1ciXwmKS3S/pXSR+TtFub9/UR4LukaqIzgd+TquJ+Luklbd5XvVrbQ+3kXzuer5Q0ccQblzYnVWOeADwN/A9wDulO8BPA+nVlNwVuBY4DHgK+Qar+3Br4qaRGFwzNOBf4Gqk67XzgDOB6YHtgvxFs13A11HhVu1Kf04ZtnQt8DDhS0ukRcVvZGyStBhxJqk75Yd12DiYlkV/m3vLTrOz+ki4HLgRuAe5p5eq3gTdkr38jnfC2z8V6PXBIRDw0wv1AOmG9of5qPmvYnw4cBFzchn2sQNKGpOMKcFP2eiPwf6TPeq2kc7N1f4iI4VQ9fRfYFDghIk7J7X894Km6RedlZadHxIV15SYB1wFnSLo8Iv7WSgCS1gIOIyWiXfKfQ9K6rWzPVuY7i/Fpg+x1ML9C0mZZtUn918cbbSgingc+Tfpb+o8m9/8eYG3goohYmi27inTCPkTS2rl9DALvAu4lNWZfAPwReDzrHXXkCK6QX569foh0h/RWUgP/dsBsYA/gB8Pcdt4ZBdU+38xed84XHqajst/ZFyR9C7iL9BlvJiVZIuLvwIHAbcCbgG+R7nKelPQLSR9p9k5HqbfcG7NtfTm/PiIejoins7I7kBrPL61PFFm5xcDngdWAd7f+sQnSHdQzpHawfByFHSeseb6zGJ/y1RL1NiP909b7P+CrjTYWEbMlXQ28TdIBEVHWc6hW1fDtum08J+kC4JPAe0lVCPX7uFbSa0h3RXuSeilNBfbNvt4v6R0R8UzJvvNqSUakO4jbs5/vlPQuUlLaU9JuEfGbFredN7dg2f3Z69oF64ajvgvx34E/AZcCp9W3vUTE74ApknYC9gZeT2pD2iP7OkbS3hHxWMn+ds1eZ2cXDkOpVeutJemkgvW16qqtS7azkoh4QtKPSRcTt0m6lHSHelNEjGbvtHHDyWJ8egDYClipYTXrSSMASasAzTbufop0dfkVSbMbFZK0NbA7qcrjxtzqb5OSxdHkkkUW2/OkE8Avs22J1GvnPNIdwYcZIqk1UDsZ3leXKGr7W5p9lg+SrvxHmiwWFyyrtQeNuO0gs3d9b6gyETGXuiQmaWfS8dyBdNHQ8K4yMyl7XekutUCtKmif7KuRNZvYVpF/IPXMOxz4QrbsaUmXAJ9qtWrLVuRqqPGp1vvpLe3aYFa9ch6wLfCPQxStNWxvpdxDZECtimY7SW9sYp8REVeTenYBvHkYoc/PXotO5PBiMlnpobZeFBE3A8dmPzZzPGvHrZkeXbUuvx+LCA3x9YFW44aU3CPipIh4DbAJqV3sV9nrJcPZpr3IyWJ8Opd0RXtIdqXfLieSekidTHpgbgVZPfh7SXXK55DqyvNftbuSVnrF1Lq7DueJ4+tJx2ILSasWrN8ue10wjG2PVa0cz9rd4b5ZN+Rmyhb1eGuriLg/Ii4gVVH+Cdjdjdwj42QxDkXEvcCXgFWBq4a4ip/UYHmj7S4iPWC3AcXVF+8mVUXMjogPRsQ/5b+AQ0l17e/JerggaT9JBxc9HCdpzbp9Xd9KvFnMD5O6bq5FelCxftv7kE42j5N6ZPUESTtnz2sUDQHSR6rKgSaOZ0TcCvwaeF3d++q3t27W+61W5fVL0tP6hXefkrZXehi0JZLWl7RLwao1SB0WniM9y2PD5DaL8etk0pXjZ4EbJN1K6jHzKClJbEZqB4DWTsJfIVU1vbpgXa0K6uxGb46IJyX9ADiKVH1wJql95XTScxC/JF0pPkcazuLtWbw3AV9vIc56nyQ9qPYZSXuQjsOmpB5Yy4Gjs946vWIjUvvQ1yX9itRj6mlgQ1L33g1IT+ef3OT2jiR1e/13Se/OvhewBfA20u9vQVb2cFKX7W9J+ijp97aY9Lt8LelObjfgwRY/0wBwo6S7SV2g7wdeBrwj+zxnFD1wac1zshinsucTTpL0fVK30b1J/8hrkKoh7iU9XPWdiPhtC9t9StLnSQ9bvUDSFqReTA8CPy7ZzDdJyeJoUrL4Lmloin1IDa97kBpBF5Ma1X8AnB0Rw7pyjIgHs6vSE0kJYlfSMfgJcEpBQ/xYdw3pd/02Ui+onUgJ9wngD6QH285s9uQaEX+WtCOpC/U0UpvH06QE8Z/UnfgjYmHW3fafSXeaR5Aa9/9KSlr/xYttV61YQGqQ34v0t7we6cJnPukBwAsbvdGao5E/02RmZr3ObRZmZlbKycLMzEq5zcJ6iqTNSO0dzfhqjzVcm40at1lYT5G0F3Btk8U3j4gFoxeNWe/oyWSRPQ1sZmYtiojChzF7thqqF5Ogmdlo0hDTrruB28zMSjlZmJlZKScLMzMr5WRhZmalnCzMzKyUk4WZmZVysjAzs1JOFmZmVqrSh/IknUOanOTBiNiuYL1IY+sfQJqu86hW5lao0omz7uD7N93P8ggmSkzfZWO+NG37qsMCYNa8QWbOns+ixUvZaFI/M/bdkmlTmplCeezufyT7rH/vWv19SPDYkmVMlFgewUCT22smhqIywArL9t5qfa79w0MjPn6z5g3yhR/fyWNLlgEwqb+Pkw7c9oVt5df3901gtb6JLF6yrGFszcYy3N/HaP3tdPP/a7eodLiPbFayp4DzGySLA0iTpBxAmsnsaxFRNHVi/n1R5ec6cdYdfPfGv6y0/MhdN6n8D3DWvEGOv+wOli5b/sKy/r6JnHLw9h1JGFXsfyT7LHpvkbLtNRNDUZm+CQLBsuWN/577Jog1V1tlhZN4M59rxiW3r7Tdvgli5qE7ABSuz5fNx9bMcR3u72O0/naO+OZvuOHeR1da3g3/r50mqeFwH5VWQ0XE9aTZrBo5iJRIIputbJKkDTsT3fB9/6b7W1reSTNnz1/pxLd02XJmzp7fs/sfyT6L3lukbHvNxFBUZtnzMeQJu1bmsSXLCGBw8VKOv+wOZs0bHPI9M2fPL9zusueDmbPnN1xfFlszx3W4v4/R+NuZNW+wMFFAd/y/dpNub7MYIM2lW7MwW7YSScdImitpbkciG8LyBnc1jZZ30qLFS1ta3gv7H8k+W4lrqLLNxNCuY9DMCbQs1pHEUvbe4f4+RuNvZ6jj1A3/r92k25NF0e1Q4W8wIs6KiJ0iYqdRjqnUxAaDcTVa3kkbTepvaXkv7H8k+2wlrqHKNhNDO49B2Qm0LNaRxDJBGvLOZri/j9H42xnqOHXD/2s36fZksRDYuO7nycCiimJp2vRdNm5peSfN2HdL+vsmrrCsv2/iC42Vvbj/keyz6L1FyrbXTAxFZfomiL6JrZ+0yk6gM/bdsnC7fRPEjH23bLi+mdiWRwxZFTbc38do/O0MdZy64f+1m3R7srgceJ+SXYHHI+KBqoMq86Vp23Pkrpu8cGUyUeqaxrJpUwY45eDtGZjUj4CBSf0da9yuav8j2Wf+vZP6+1h79T7gxSvPZrbXTAxFZWYeugMzD9lhhWVH7rrJCvHkT9jNnECnTRlg5iE7vPBZyLY189AdmDZloHB9f98E1l69b6XYiq7Ah6oKG+7vYzT+dhpdDEx91Tpd8f/aTaruDfV9YC9gPeBvwOeBPoCI+EbWdfbrwH6krrMfiIjSNomqe0OZdVLVXaE3P+4nhXXDAv586ts7FsdwVX38uslQvaF6dqa8XvxcZt1o6qlzGCyo+x+Y1M8Nx725gohsuLq266yZjX1Vt4NZZ/TstKpm1hm1KhtX5fQ2V0OZmRngaigzMxshJwszMyvlZGFmZqWcLMzMrJSThZmZlXKyMDOzUk4WZmZWysnCzMxKOVmYmVkpD/dhNoZ4hFSriof7MBsjZs0b5PjL7lhhHur+vokdnY/EepuH+zDrATNnz18hUUBz822btYOThdkY0Wi+6LL5ts3awcnCbIxoNF902XzbZu3gZGE2RniSIauSe0OZjRGeZMiq5N5QZtYV3C24ekP1hnKyMLPKuVtwd3CyMLOuNvXUOQwW9OoamNTPDce9uYKIxpZ23ZUNlSzcZmFmlXO34OHL35UNLl7K8ZfdAdDWuzL3hjKzyrlb8PB16mFNJwszq1wnuwXPmjfI1FPnsPlxP2HqqXOYNW+w7fvopE7dlbkayqzLjMdeQZ3qFtypKptO2mhSf2F7T7vvytzAbdZF3CtodPViQ3o7/2Y8kKDZGOHBAkdXLzakT5sywCkHb8/ApH5ESnyjcXHhaiizLtKLJ7Nu0qkqm06bNmVg1O88fWdh1kXcK2h0eXyt4as8WUjaT9J8SfdIOq5g/SaSrpU0T9LvJB1QRZxmneCT2ejqVJVNL6q0gVvSROCPwD7AQuAWYHpE3FVX5ixgXkT8j6RtgCsjYrOS7bqB28as8dgbyrpDNz/BvTNwT0TcByDpQuAg4K66MgG8LPt+LWBRRyM067BO1D+btarqZDEA3F/380Jgl1yZk4CrJf0zsAbw1qINSToGOGYUYjQbNt8lWK+ous2i6HYnX380HTg3IiYDBwDfkbRS3BFxVkTsFBE7jUKcZi2r9X8fXLyU4MUHwMb6E8M2PlWdLBYCG9f9PJmVq5k+CFwMEBG/AVYD1utIdGYj4GcmrJdUXQ11C7CFpM2BQeAw4PBcmb8AbwHOlbQ1KVk81NEozYZhuM9MuOrKulGldxYR8RxwLDAbuBu4OCLulHSypAOzYv8CHC3pduD7wFHu6mRjwXCemXDVlXUrjw1lNkqGM2ZPp8cu8l2M1evmrrNmPWs4I6l2criPXhyB1UaPk4XZKGr1mYlOjl00VAO8k4XlVd0byszqdHK4Dw9aaK1wsjDrIp0cu8iDFlor3MBtNk55oiXLcwO3ma1ktKYydQ+r3uQ7CzNrSlESgBWTzd5brc+ltw76bmWMGurOwsnCzEoVVVn1TRAIli1/8X9NrDy4G4ztOa7HE1dDmY1j7agWKupmu+z5ldNCo0s097Aa+5wszHpYux68G+nJ3j2sxj53nTXrYe0a+baVk32+DsPTwvYGJwuzHtauB++KHhbsmyD6Jq6YGvr7JnLErpt4juse5Goosx7WruFDGnWzLVrmxNCb3BvKrIf5wTtrhXtDmY1To/XgnY0/vrMws67mJ8I7xw/lmdmY5Gq0zhoqWbg3lJl1rXZ1/bWRc7Iws67lOTe6h5OFmXUtz7nRPdwbynrOeGwQ7abP3M5YZuy7ZWGbhZ8I7zw3cFtPGY8Not30mU+cdQcX3PiXFQYUHGks3ZQIe517Q9m4MfXUOYVPLPfyENnd8plnzRvkExfd5iHKxzD3hrJxYzw2iHbLZ545e76HKO9hThbWU8Zjg2i3fOahEkIvH//xwsnCekrR6Ki93iDayc88a94gU0+dw+bH/YSpp85h1rzBF9Y1SgjKYrSxzb2hrKeMx7GQOvWZyyZSKuq5JOCIXTdpKRY3aHcnN3CbWVOaaUgf6Ym+m3p2jUceddZsnGrnVXozDenTpgyM6KQ+1PAeThbVcpuFWY+qXaUPLl5K8GK1UX07Qys60ZDeLT27bGWVJwtJ+0maL+keScc1KPMeSXdJulPS9zodo9lY1O5B+DrRkN4tPbtsZZUmC0kTgTOB/YFtgOmStsmV2QI4HpgaEdsCH+94oGZjULuv0qdNGeCUg7cf1fm1x2NvtrGi6jaLnYF7IuI+AEkXAgcBd9WVORo4MyIeA4iIBzsepdkY1K75t+uNtE2ime3D+OrNNlZUnSwGgPvrfl4I7JIr8xoASTcAE4GTIuKn+Q1JOgY4ZpTiNBtzxuogfKOdkGx4qk4WRV208n1eVwG2APYCJgO/lLRdRCxe4U0RZwFnQeo62/5Qzaox3B5Nvkq3dqo6WSwENq77eTKwqKDMjRGxDPizpPmk5HFLZ0I0q07Zg3BlfJVu7VJ1b6hbgC0kbS5pVeAw4PJcmVnA3gCS1iNVS93X0SjNKjKSHk1DDc1h1qpK7ywi4jlJxwKzSe0R50TEnZJOBuZGxOXZurdJugtYDsyIiEeqi9qsfcqqmBr1XBpcvJSpp85p+L6R3pGY5Xm4D7OKNDO0RaMhNvKafZ/nlbCheD4Lsy7UTBXTjH23pG9C4f/ukO/r5JPQru4aH5wszCrS7FhLa67WXG1x/fuGehK6nSf3dg8pYt3LycKsIs0ObbF4ybKWt9foSei9t1q/rSf3dg8pYt3LycKsIs0ObdHME9f59zUamuPaPzzU1pO7B/4bP6p+zsJs3Gr2obmiJ7H7Jog1V1uFxUuWNXxf0TMWn7jotsJYhntyH40hRaw7OVmYVaiZh+ba+SR2u0/uY3VIEWtdaddZSS+NiCc7FE9buOusWbHRmInO06D2jqG6zjaTLH4BPAE8RBq3aRFwRUTc1O5A28XJwqwxn9ytkZEmiyOB83hx0L/aG74QESe3Lco2crIwM2vdiJJFtoGdgANJA/2tSxpGfH3ShEQ3tjHWtnCyMOttvjsaHSO9s/hIRPx3blkf8FWgLyKOya37cET8zwhjHhEnC7PeNRrtLpaMdLiPT+cXZMOFzwC2Lyj/idbCMzNrnh8ErEYzXWc3kTQLmMfKExO9XNLnsu8FbAu8qo3xmZmtwA8CVqPZ5yzeSWqzyAvgpLrvxcoJxcysbfwgYDWaTRZXAXNZMREIeB+ppxSk+Si2Aw5qW3RmZjl+ELAazTRw/yUiNilY3g9cGxG75pb/KSK2aG+YrXEDt1lvc2+o0THS3lD/GBHn5JatApwGrFbQG+oDEfHtEcY8Ik4WZmata8dzFjuS2i02BtYGdgNeAbwpIn7dxljbwsnCzKx1QyWL0jYLSUew4hPctdcvdmOiMDOz9mt2bKingAdJDdwPAD/uxie3a3xnYWbWupG2WbwsIp4YlchGiZOFmVnrRtxmMdY4WZiZtW6kw32Ymdk452RhZmalPK2qWY/wg2o2mtxmYdYDPGy3tYPbLMx6nIftttHmZGHWAzxst402JwuzHtBoeG4P223t4mRh1gNm7Lsl/X0TV1jmYbutndwbyqwH1Bqx3RvKRkvlvaEk7Qd8jTR50tkRcWqDcocAPwDeEBFzS7bp3lBmZi3q2t5QkiYCZwL7A9sA0yVtU1DupcBHgZs6G6GZmUH1bRY7A/dExH0R8SxwIcXTsn4R+ArwdCeDMzOzpOpkMQDcX/fzwmzZCyRNATaOiCuG2pCkYyTNlTRkFZWZmbWu6gbuorqxFxobJE0ATgeOKttQRJwFnJW9zw0W1jYeRsOs+mSxkDRVa81kYFHdzy8FtgOukwSwAXC5pAPLGrnN2iE/jMbg4qUcf9kdAE4YNq5U2htK0irAH4G3AIPALcDhEXFng/LXAZ9ybyjrlKmnzmGw4CnogUn93HDcmyuIaGzzXVp3G9Ec3KMpIp6TdCwwm9R19pyIuFPSycDciLi8yvjMPIxG+/gubWyruhqKiLgSuDK37HMNyu7ViZjMajaa1F94Z9GJYTR67Sp8qMEOx/LnGi+q7g1l1tWqGkajdhU+uHgpwYtX4bPmDY7qfkeT79LGNicLsyFMmzLAKQdvz8CkfkRqq+jEHBG9OOS4Bzsc2yqvhjLrdtOmDHS8mqQXr8Jn7Ltl4QRNHuxwbPCdhVkX6sWr8Kru0qw9Kh9IcDQMt+tsrzUo2tjlaVKtCkN1nXWyyPif07qNL16s05wsmuCHr8xsvOvaIcq7SS82KJqZtYuTRaYXGxTNzNrFySLjOYzNzBrzcxYZz2FsZtaYG7jNzAxwA7eZmY2Qk4WZmZVysjAzs1JOFmZmVsrJwszMSjlZmJlZKScLMzMr5WRhZmalnCzMzKyUk4WZmZVysjAzs1JOFmZmVsrJwszMSjlZmJlZKScLMzMr5WRhZmalnCzMzKyUk4WZmZWqPFlI2k/SfEn3SDquYP0nJd0l6XeSrpG0aRVxmpmNZ5UmC0kTgTOB/YFtgOmStskVmwfsFBGvBS4BvtLZKM3MrOo7i52BeyLivoh4FrgQOKi+QERcGxFLsh9vBCZ3OEYzs3Gv6mQxANxf9/PCbFkjHwSuKloh6RhJcyXNbWN8ZmYGrFLx/lWwLAoLSkcCOwF7Fq2PiLOAs7KyhdswM7PhqTpZLAQ2rvt5MrAoX0jSW4HPAHtGxDMdis3MzDJVV0PdAmwhaXNJqwKHAZfXF5A0Bfhf4MCIeLCCGM3Mxr1Kk0VEPAccC8wG7gYujog7JZ0s6cCs2ExgTeAHkm6TdHmDzZmZ2ShRRO9V70uKXvxcZmajSRIRUdSWXHk1lJmZjQFOFmZmVsrJwszMSjlZmJlZKScLMzMr5WRhZmalnCzMzKyUk4WZmZVysjAzs1JOFmZmVsrJwszMSjlZmJlZKScLMzMr5WRhZmalnCzMzKyUk4WZmZVysjAzs1JOFmZmVsrJwszMSjlZmJlZqVWqDsDMmjdr3iAzZ89n0eKlbDSpnxn7bsm0KQNVh2XjgCKi6hjaTlL04uey8W3WvEGOv+wOli5b/sKy/r6JnHLw9k4Y1haSiAgVrXM1lNkYMXP2/BUSBcDSZcuZOXt+RRHZeOJkYTZGLFq8tKXlZu3kZGE2Rmw0qb+l5Wbt5GRhNkbM2HdL+vsmrrCsv28iM/bdsqKIbDxxbyizDhlpT6ZaWfeGsiq4N5RZB7gnk40F7g1lVjH3ZLKxzsnCrAPck8nGusqThaT9JM2XdI+k4wrWv0TSRdn6myRt1vkozUbGPZlsrKs0WUiaCJwJ7A9sA0yXtE2u2AeBxyLi1cDpwJc7G6XZyLknk411Vd9Z7AzcExH3RcSzwIXAQbkyBwHnZd9fArxFUmEDjFm3mjZlgFMO3p6BSf0IGJjU78ZtG1Oq7jo7ANxf9/NCYJdGZSLiOUmPA+sCD9cXknQMcMzohWo2MtOmDDg52JhVdbIoukPI93ltpgwRcRZwFqSusyMPzczMaqquhloIbFz382RgUaMyklYB1gIe7Uh0ZmYGVJ8sbgG2kLS5pFWBw4DLc2UuB96ffX8IMMdP3JmZdVal1VBZG8SxwGxgInBORNwp6WRgbkRcDnwL+I6ke0h3FIdVF7GZ2fjk4T7MzAzwcB9mZjZCThZmZlbKycLMzEpV/ZzFqPFD3mZm7dOTDdztImluROxUdRzNGEuxwtiK17GOjrEUK4yteEcjVldDmZlZKScLMzMr5WQxtLOqDqAFYylWGFvxOtbRMZZihbEVb9tjdZuFmZmV8p2FmZmVcrIwM7NSThY0NQ/4UZIeknRb9vVPFcV5jqQHJf2+wXpJOiP7HL+TtGOnY8zFUxbvXpIerzuun+t0jHWxbCzpWkl3S7pT0scKynTF8W0y1q44tpJWk3SzpNuzWL9QUOYlki7KjutNkjbrfKRNx9oV54K6eCZKmifpioJ17T2uETGuv0ij3d4LvBJYFbgd2CZX5ijg610Q6x7AjsDvG6w/ALiKNGHUrsBNXR7vXsAVVR/XLJYNgR2z718K/LHg76Arjm+TsXbFsc2O1ZrZ933ATcCuuTIfAb6RfX8YcFEXx9oV54K6eD4JfK/od93u4+o7i+bmAe8KEXE9Q0/8dBBwfiQ3ApMkbdiZ6FbWRLxdIyIeiIjfZt8/CdxNmtK3Xlcc3yZj7QrZsXoq+7Ev+8r3qjkIOC/7/hLgLapgCIYmY+0akiYDbwfOblCkrcfVyaJ4HvCif7x3Z1UPl0jauGB9N2j2s3ST3bLb/qskbVt1MADZ7foU0pVlva47vkPECl1ybLOqktuAB4GfRUTD4xoRzwGPA+t2NsqkiVihe84FXwU+DTzfYH1bj6uTRXNzfP8Y2CwiXgv8nBezdbdpar7yLvJbYNOI2AH4L2BWxfEgaU3gUuDjEfFEfnXBWyo7viWxds2xjYjlEfE60rTJO0vaLrGdW7gAAAumSURBVFeka45rE7F2xblA0juAByPi1qGKFSwb9nF1smhiHvCIeCQinsl+/Cbw+g7F1qpm5jTvGhHxRO22PyKuBPokrVdVPJL6SCffCyLisoIiXXN8y2LttmObxbEYuA7YL7fqheMqaRVgLSquvmwUaxedC6YCB0paQKo6f7Ok7+bKtPW4Olk0MQ94rl76QFIdcTe6HHhf1mtnV+DxiHig6qAakbRBrQ5V0s6kv8dHKopFpCl8746I0xoU64rj20ys3XJsJa0vaVL2fT/wVuAPuWKXA+/Pvj8EmBNZq2wnNRNrt5wLIuL4iJgcEZuRzllzIuLIXLG2HteeHaK8WdHcPOAflXQg8BwpMx9VRaySvk/q5bKepIXA50mNcETEN4ArST127gGWAB+oIs6aJuI9BPiwpOeApcBhVZwkMlOB9wJ3ZHXWACcAm0DXHd9mYu2WY7shcJ6kiaSEdXFEXJH7//oW8B1J95D+vw6rIM5mY+2Kc0Ejo3lcPdyHmZmVcjWUmZmVcrIwM7NSThZmZlbKycLMzEo5WZiZWSknC7McSedKiqpGPzXrRk4WtoJsbJyjJf1C0qOSlikNM/47SWdnfczry++VnVivqyhkM+uAcf9Qnr0oexjpCtIQB4uBn5CGDFgHeBVwOLAVuSfczaz3OVlYvemkRHE7sGdEPF6/UtLqwC5VBGZm1XI1lNV7Y/Z6bj5RAETEkoi4dqQ7qVVbZeMXnS1pUNJySUfVlVld0vFKs5H9XdJTkn4jaXrB9laVdKykKyX9n6Rnsiq0n0vaf4g43irpl9n2H5U0S9JWQ5Q/UNI1kh7I9rEoq677SIuf/x+y7Twq6WlJCyR9X9JOuXIvkXRcVgW4RNITWbzvyZXbLTumRYMf1srcncW8Tm75vtlxezhbf6+kmbUxknJlF2RfL5N0Wvb9MkknZes3kvQ5STdI+qukZ7Nj9D1JWzeIS5I+Jumu7FgMSvq6pLVq+2vwvulKswU+lr3vbkknSnpJo2NgI+M7C6tXG2juNR3Y1zrAjcBTwGWkMfn/BpCdqOaQ5mn4LXAO6cJmX+B7kraNiBNz2/oa8GvgZ8BDpHF+3glcKenoiFhhghhJhwAXAc9mrw8AuwO/AX6XD1bSMcD/An8lDVP9MPBy4LWkMaL+u+wDSxLwbdLgbg9nn/sh0ui1ewPzgblZ2VVJ45XtSRrM7kxgddKYTxdJel1EnAAQEb+RNB94h6R1I+KR3H53JlUfXhoRj9Yt/xzwBdK4QVeQ5nB4LfAp4ABJuxUMfb4q6XezDnA18ATw52zdHsBxwLWkEXGfArbIYj5Q0tSIuD23vTOBD5NG7z2L9Ps4kDQpWR+wrOA4fgv4R1IV6WWkKtNdgS+SJvjZJ5u/wdppJNPs+au3vkgn52dJJ+7vAAeT5kQY6j17kcbIv66F/UT2dT6wSsH6c7P1n84tXw34aRbf6+qWvwSYXLCdtYDfk06G/XXL1yQlxmXATrn3nF4X32Z1y28FngFeXrCf9Zr83Mdk270ZWCu3biKwYd3Px2dlr6w/RqQEtSBb98aC8scW7PfMbN0765btnS37NTApV/6obN3pueW1/f4cWKNgPy8HXlqwfAdS4rgqt/xN2fbm18dASkjXZ+sWNIjtsvrfabbupGzdx6r+X+rFr8oD8Fd3fQHvIV1lR93XI8AP6082deWHmywanXjXJY3oeUuD9+6Qvf8rTe7rk1n5PeqWHZEtO6+g/FqkK9WiZPF3YO0RHNs7su1OaaLsn0hJcauCdR/MtnNO3bLJwPL8cctOvI+Q7trqk84Ps21s22D/80iT69QvqyWLHYbx2S8Hngb66padnW3vfQXlpzZIFvNISX5SwXsmku7Ybh6t/4/x/OVqKFtBRFws6YekK8/dSXcbuwPTgGmSzgeOiuy/cwQWRMSDBcvfQPqnj1pdeE5f9rpCHbjStKEzSFUhG5LuQurVT3+6Y/b6i/zGI+JxpWG/98ytugD4T+BOSRdl770hIh4qiHElktYAtgP+FhHzSsq+FHg1MBgR+bkfIFUDQfrd1OJeKOkaYB9J20TEXdmqd5KqjE6PFatmdiOddA+VdGjBPlYF1i+o1nqagmq6utjfDnwI2AlYj5WrutcjXYzUx/+rgk3dSLpoqN/26qSLhYeBj6t4OulnyP1tWHs4WdhKImIZqT76anihS+27SW0H7yNdlY50ms6/NlhemyP4DdlXI2vWvlGaiGgO6e/5GtJV7BNk1VWkievrGz7Xyl7/1mxsEXGapIeBjwAfBT5OSmi/AGZExNwhYgWoNRgPlpSrj6/RxEq15flG6HOBfUhtIv+aLatNfpOf/nNd0vH6fEkstSq7mgcbXShI+iip7egxUtvRX0jzfgTpYmMHmvw9RMRySfnJmtYmTRW6fhNxW5s5WVipiFgOXCxpe+BE4M2MPFk0ujOp9cI6PSI+2eS2TgT6gb0j4rr6FZKOJyWLon28osH2NihaGBHnA+dnDfBvBN5FamidLWnrBndKNYuz14EhyuTjK4yDdOdUX67mh6QkeaSkE0h3FPsDt8fKDcuPAxMiYh1a0yhRrEJqLP8rsGPkZhCUtFvB22qN568A7suVn0hKaPXJtfZ550XEjlhHueusteLJ7LXw/r9NbibdEbyphfe8Gng0nygy+eokSD2sCtdJWot0N9JQRCyOiCsj4mjS1fw6ZfFGxN9Jje2vkDSlpOyTwL3AgKQtCorsnfsctfctBS4GNiJNCXoE6YIwf1cBqZpn7az6rh3WI93p/LogUazJi1V/9WrVcbsXrNuV3MVspDnF7wS2zXcBttHnZGEvyPqu7yNppb8LSRsAR2c/Xj9aMWRX5xcAO0n6bHbFmo/lVZI2r1u0AFhH0mtz5T5I6m6b9yNSVcnh+WcbSD1q1sq/QdJ+RbGQegBBqm4pc0b2+r9ZUqrf/gStOL/zOaSkPDO7yq6VWw/4bF2ZvHOz1/dlX8+Rjmfe6dnrNyVtlF8paY2seq9ZD5KOweuz5FDbTh+pamq9gvecn71+pv54ZN2G/73Bfk4jtaec0+BZkLUl+a5jFLgayurtAnwM+KukX/Fi//nNgbeTqnp+BFwyynEcS+qffzLw3iyWv5GumLcmtWVMr4vvq6Sk8CtJF5OqK3YiXbFeQurn/4KIeCp7buIi4JdZg3XtOYvtSMlwj1xMFwJPZ7EsIJ3I35TFciupO2mZs7N9vA/4k6QfkZ6z2IhUtXcOKVkB/AepCukg4HZJV5KesziUlKC+EhErNQxHxA1Kcy4fSuoM8OOi6rGIuEbSccApWSxXko7nmsCmpLuuX5Ge6C8VEc9LOoP0nMUd2WdblXQXtA7p2Yu9c+/5haSzSF2K75R0KanR/Z2k3+Ei0l1m/XvOkfR6UtvRvZJmk9pG1iH9ne5BepblQ83EbS2oujuWv7rnC9gY+H+kuu/5pDrlZ0kn0iuBI0n13PXveQupHvvqFvZT2tWWdKI5lvQcwOOkXi5/ITVgfxxYN1f+HaSqlSdJ7QNXk04cR2X7O6pgH/uQTohLSHcaPyI9vHYuK3ed/VB2XO7Lyj9Kqkb5NAXPFpR8tiNIvakeJ/Uu+jPp6n/HXLnVgBNI1VdLs8/2K2B6yfZP5MVuz+8uKbs7qepqUfa7fgi4jXQFn38GZQG5rqy59auQuirflcX7V9LzOpsWHdPsPROAT5AePHwmi+NM0t3dk8BtDfb1Dl58kPDZbF83A1+ioLuxv0b+pezAmw2L0vAb3wMuiIgjq47HekPWVvNH4MKIWGmIF+s8t1nYSL0re/1NpVHYmKQ0PtiE3LLVSVWLkO7mrAv4zsJapjQp0IdID3btQereuG0UDD5oNhRJp5Lan64jVXduQKranAxcBbw9fJLqCm7gtuHYjFRX/zCp4fcEJwobpp+RHtZ7G6mR+jlS9dMZwFedKLqH7yzMzKyU2yzMzKyUk4WZmZVysjAzs1JOFmZmVsrJwszMSv1/agnPytooFpwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.log10(song_read_counts.loc['GNAS_6', song_iPSC]+1), song_PSI.loc['GNAS_6', song_iPSC])\n",
    "plt.title('GNAS_6 in iPS cells', fontsize=20)\n",
    "plt.ylabel('$\\Psi$', fontsize=20)\n",
    "plt.xlabel('SJ reads coverage', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEkCAYAAADeqh2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debgcZZn+8e9NCHpwIQzgKIcto8giiCirqICILI7AIC4gjDgOzG8c3EVRUSKoIKgoI+hExTCCoiJmgiJhRsAFBYkGZAJGI6DkhBERAggBk/j8/nirodOp7uq9erk/19VXn1P1dtVb3efU0++uiMDMzKyRdcrOgJmZDT4HCzMzK+RgYWZmhRwszMyskIOFmZkVcrAwM7NCDhZmXSQpJF3Tw+Mfm53j2F6dY5DUu15Jd0i6o5xcjScHizEn6dmSPiXpF5LulbQye75e0ickvSDnNbOyf+CQ9OY6x638k3+kwbnXl7Q8S/fVJvL6aklXSLo7y+efJN0i6UJJb2jtynOP/2JJ35J0l6RHs+crJR3c6bF7TdI1VZ9J5fGgpJ9Ler+kiZzX7C/p25KWSfqLpPsk/VrSNyW9VZLKuBYbTOuWnQErR3Yj+FD2WAf4BfB14F7gKcBzgbcA75J0QkScW+dQsyRdGBEPtJGN1wIbAAEcLmmjiPhTnfzOBo4DVgDfBW4HngT8HfBKYB/ggjbyUDn+ycBpwD3Ad4C7gI2BnbNjX97kobYDHm43H034NnAdKX95LgDuAARsBhwOfBQ4VNKLImIlgKT3Z9tXAVcAi4HpwExgb+AI4Lxsv5mDxRj7EDALuBM4MiKurU0g6WnA20k39DxLgGcB7wU+0EYejgf+CnwSOBF4A/CpnHzsRQoUS4E9I2Jpzf7ppBt6WyS9mhQo/gc4PCIezDl+UyLiV+3mo8nj3w/c3yDJnIi4pvJLFgQXArsBRwEXSNoSOBV4AHhRRNxcfQBJ6wD7A6u7m3sbZq6GGkOS/g44GfgLcFBeoACIiLsj4v3AmXUO9e/AMuAdkjZrMQ87AHsA3wc+nuXluDrJ98qev1UbKLJ8royI/27l/FX5WCc7/8PAUbWBonL8Fo63VptFVbXdPpKOkPQzSQ9n1X0XS5ps4fgttVlExF3Apdmvu2XPuwPTgKtrA0X2mr9GxPxoYS4gSX8j6aOS/je7tvsl3STpDElPykl7uqRbJa3I0n5f0subPV+dPKyXVZ/9IqtSezhr2/gvSS/r5NjmYDGu3kgqVV4SEYuKEkdEvaqIh4EPAhOkKo1WHJ89z8mqnr4DbCvpxTlpK1VTz27xHM14Ianq5XLgPkmvkPReSW+TtGeXz/Vm4EJSNdG5wP+SquL+R9ITunyuapW2h8rNv/J+/p2kaR0fXJpJqsZ8P/AI8DngfFJJ8B3AJlVptwR+DpwE/BH4PKn6czvgCkn1vjA0Yw7wGVJ12n8C5wA/BHYEDuzguIarocZV5Zv6VV041hzgbcDRks6OiBuLXiDpicDRpOqUb1cd53BSEPlRzUuuyNIeJGkecDFwA7CklW+/deyaPf+BdMPbsSavPwSOiIg/dngeSDesXau/zWcN+0cChwLf6MI51iDpGaT3FeD67Pk64Heka71a0pxs368iop2qpwuBLYH3R8TpNeffGPhz1aYLsrRHRsTFVelmANcA50iaFxF/aCUDkjYAXkcKRLvXXoekjVo5nq3NJYvx9PTseap2h6StsmqT6sfb6x0oIv4KvIf0t/SJJs//GmBD4OsRsSLb9j3SDfsISRvWnGMK+Afgt6TG7IuAXwP3Z72jju7gG/LTsuf/RyohvYzUwL8DMB94CfDNNo9d65ycap8vZM+71SZu07HZZ/ZhSV8CbiFd489IQZaIeAg4BLgReDHwJVIp50FJP5D05mZLOkq95V6YHevjtfsj4p6IeCRLuxOp8fxb1YEiS7ccOAV4IvCq1i+bIJWgHiW1g9XmI7fjhDXPJYvxVFstUW0r0j9ttd8Bn653sIiYL+lK4OWSDo6Iop5DlaqGL1cdY5Wki4B3AseQqhCqz3G1pGeTSkV7k3op7QUckD3eIOnvI+LRgnPXqgQZkUoQN2W/L5L0D6SgtLekPSPipy0eu9aCnG13Zs8b5uxrR3UX4oeA3wDfAj5V3fYSEb8Edpa0C7Av8AJSG9JLssfxkvaNiPsKzrdH9jw/++LQSKVabwNJs3L2V6qrtis4zloi4gFJl5G+TNwo6VukEur1EdHL3mljw8FiPN0FbAus1bCa9aQRgKR1gWYbd99N+nZ5pqT59RJJ2g54EanK47qa3V8mBYvjqAkWWd7+SroB/Cg7lki9di4glQj+lQZBrY7KzfC2qkBROd+K7FreRPrm32mwWJ6zrdIe1HHbQWbf6t5QRSJiAVVBTNJupPdzJ9KXhrqlysyM7HmtUmqOSlXQ/tmjnic3caw8ryX1zDsK+HC27RFJlwDvbrVqy9bkaqjxVOn9tF+3DphVr1wAPAf4pwZJKw3b26pmEBlQqaLZQdILmzhnRMSVpJ5dAC9tI+uLs+e8Gzk8HkzWGtQ2iiLiZ8AJ2a/NvJ+V962ZHl2VLr9viwg1eLyx1XxDCu4RMSsing1sQWoX+3H2fEk7x7THOViMpzmkb7RHZN/0u+VkUg+pU0kD5taQ1YMfQ6pTPp9UV177qJRKWukVU+nu2s6I4x+S3outJa2Xs3+H7PmONo49rFp5PyulwwOybsjNpM3r8dZVEXFnRFxEqqL8DfAiN3J3xsFiDEXEb4GPAOsB32vwLX5Gne31jruMNMDu6eRXX7yKVBUxPyLeFBH/XPsAXk2qa39N1sMFSQdKOjxvcJykJ1ed64et5DfL8z2krpsbkAYqVh97f9LN5n5Sj6yRIGm3bLxG3hQg00lVOdDE+xkRPwd+Ajyv6nXVx9so6/1WqfL6EWm0fm7pU9KOSoNBWyJpE0m75+x6EqnDwirSWB5rk9ssxteppG+OHwSulfRzUo+Ze0lBYitSOwC0dhM+k1TV9KycfZUqqC/We3FEPCjpm8CxpOqDc0ntK2eTxkH8iPRNcRVpOotXZPm9HvhsC/ms9k7SQLUPSHoJ6X3YktQDazVwXNZbZ1RsSmof+qykH5N6TD0CPIPUvffppNH5pzZ5vKNJ3V4/JulV2c8CtgZeTvr87sjSHkXqsv0lSW8lfW7LSZ/lc0kluT2Bu1u8pkngOkm3krpA3wk8Ffj77HrOyRtwac1zsBhT2fiEWZK+Ruo2ui/pH/lJpGqI35IGV30lIn7RwnH/LOkU0mCrx0jamtSL6W7gsoLDfIEULI4jBYsLSVNT7E9qeH0JqRF0OalR/ZvAFyOirW+OEXF39q30ZFKA2IP0HnwXOD2nIX7YfZ/0Wb+c1AtqF1LAfQD4FWlg27nN3lwj4nZJzyd1oT6M1ObxCClAfJKqG39ELM26276FVNJ8Palx//9IQevfebztqhV3kBrk9yH9LW9M+uKzmDQA8OJ6L7TmqPMxTWZmNurcZmFmZoUcLMzMrJDbLGykSNqK1N7RjE+PWMO1Wc+4zcJGiqR9gKubTD4zIu7oXW7MRsdIBotsNLCZmbUoInIHY45sNdQoBkEzs15Sg2XX3cBtZmaFHCzMzKyQg4WZmRVysDAzs0IOFmZmVsjBwszMCjlYmJlZIQcLMzMrVOqgPEnnkxYnuTsidsjZL9Lc+geTlus8tpW1Fdo1d+EUZ81fzLLlK9h0xgQnHrANh+3czBLD3TkXUHj+fuax0fnayUerr2km/clzb+Zr19/J6gimSRy5++Z85LAdmzrmjPWnEwH3r1jJpjMm2HfbTbj6V3/s2ns7d+EUH75sEfc9vBKAGRPTmXXIc1o6ZqP3oHbfvttuwnduuovlK9L5Nlx/Oqe8svH5+v33NIqaeQ+H+X0udbqPbFWyPwP/WSdYHExaJOVg0kpmn4mIvKUTa18X7V7X3IVTvO/Sm1mxcvVj2yamT+P0w3fs+oead67p6wgEK1c/nv/a8/czj43O96oXTPKtn0+1lI9W895M+pPn3syF1/1+rdcevccWuQEj75iNdPLezl04xYmX3LTG5wnpcz7r1Ts1dcxG7wHQ1LVMnybOOiL/fL34e+rkpjiMN9Rm3sN+/9+2Q1Ld6T5KrYaKiB+SVrOq51BSIIlstbIZkp7RyzydNX/xWv94K1au5qz5i/tyrpV/jbVuLLXn72ceG53va9ff2XI+Ws17M+m/dv2dua+ttz3vmI108t6eNX/xWp8npM+52WM2eg+avZaVq+ufr9t/T5Wb4tTyFQQwtXwF77v0ZuYunOrpa8vUzHvY7//bbhv0NotJ0lq6FUuzbWuRdLykBZIWdHLCZctXtLS9F+cqStvPPDY67uo6pbdG+Wg1781sr5ePdvJXT7vvbTvvRbPpli1f0fbfUCfbi3RyUxzWG2oz72G//2+7bdCDRV5xKPcOEBGzI2KXiNilkxNuOmOipe29OFdR2n7msdFxp9WZdKxRPlrNezPb6+WjnfzV0+5728570Wy6TWdMtP031Mn2Ip3cFIf1htrMe9jv/9tuG/RgsRTYvOr3zYBlvTzhiQdsw8T0aWtsm5g+7bGG516fa/o6Yvq0NW9ytefvZx4bne/I3TdvOR+t5r2Z9EfuvnntyxpuzztmI528tycesM1anyekz7nZYzZ6D5q9lunT6p+v239PndwUN5iY3tL2QdHMe9jv/9tuG/QpyucBJ0i6mNTAfX9E3NXLE1YamvrRwFbvXEXn72cei863y5Z/01I+Ws17M+krjdjN9oaqPWYve0NVXtdJb6hm3oNOekN1++/pxAO2yW3IbeamWG+G7AYzZw+EZt7Dfv/fdlvZvaG+BuwDbAz8ATgFmA4QEZ/Pus5+FjiQ1HX2jRFR2CbRSW8oM+tcuz2aZp703dx6ZgG3n/GKrufT1tSoN9TIrpQ3itdlNur2OuMqpnLaJyZnTHDtSS8tIUfjZWC7zpqZVRv2ev1RNuhtFmY2Roa9Xn+UuRrKzMwAV0OZmVmHHCzMzKyQg4WZmRVysDAzs0IOFmZmVsjBwszMCjlYmJlZIQcLMzMr5GBhZmaFPN3HGBvGtY7NrBye7mNMDcPi8WbWX57uw9YyrGsdm1k5HCzG1LCudWxm5XCwGFPDvni8mfWXg8WY8iIzZtYK94YaU15kxsxa4d5QZjbU3AW8exr1hnKwMLOh5S7g3eWus2Y2kmbNW+Qu4H3iYGFmQ2nuwimWr1iZu89dwLvPwcLMhlKj0oO7gHefe0OZ2VBqVHoYly7g/Wzcd7AYEu7xYbamTWdMMJUTMDZcf/pY/G/UNu5PLV/B+y69GaAn1+9qqCFQ+aOYWr6C4PE/irkLp8rOmllp6g0sPeWVzykpR/3V7/ndHCyGgCf9M1vbYTtPcvrhOzI5YwIBkzMmxqrLbL/nd3M11BDwpH9m+Q7beXJsgkOtetVwvWrcd8liCHjSPzOr1e/53UoPFpIOlLRY0hJJJ+Xs30LS1ZIWSvqlpIPLyGeZPOmfmdXqdzVcqdN9SJoG/BrYH1gK3AAcGRG3VKWZDSyMiM9J2h64PCK2KjjuyE334d5QZtZrjab7KLvNYjdgSUTcBiDpYuBQ4JaqNAE8Nft5A2BZX3M4IMa5btbMyld2sJgE7qz6fSmwe02aWcCVkt4CPAl4Wd6BJB0PHN+DPNqAcOnKrDxlt1nkFXdq64+OBOZExGbAwcBXJK2V74iYHRG7RMQuPcinlcxjTczKVXawWApsXvX7ZqxdzfQm4BsAEfFT4InAxn3JnQ0MjzUxK1fZ1VA3AFtLmglMAa8DjqpJ83tgP2COpO1IweKPfc2lla6ZsSaupjLrnVJLFhGxCjgBmA/cCnwjIhZJOlXSIVmydwHHSboJ+Bpw7Mh1dbJCRWNNXE1l1lteKc+GQtGKaHudcVXuaNbJGRNce9JL+5nVoeLSmFUb5K6zZk2p3MDq3dg8JUrr+j1rqQ03BwsbGo3GmvR7npxR0KjTgIOF1Sq7N5RZV3hKlNa5NGatcLCwkTDu01W3wxNUWivcwG02poo6Ddj4cQO3WRX3AEqKOg2YVXPJwsaKv033RqMA7OA8PBqVLBwsbKx4PEb3NQrAgIPzEHE1lFnGPYC6/02/aN4ud88dDQ4WNlbGfTxGLwbitROAxyk4jwp3nbWxMu7jMXoxe2+jLrjunjs6HCxsrIz7eIxeVMM1CsDjHpxHiauhbOyM8xK1vaiGa6YLrntDDT/3hjIbI+46bI24N5SZAR6IZ+1zycKsBzwQzYaRB+WZ9ZGremxYNQoW7g1l1mW96J5qVjYHC7Mu8yhxG0UOFmZd5oFoNoocLGwgzF04xV5nXMXMk77LXmdcxdyFU11JW4a8gWgADz26auDyatYsN3Bb6VppEC5KOyi9kOYunOLDly3ivodXrrHdDd02yNzAbQOtlQbhRmkrgWRq+QqCxyfJK+Pb/GE7T7L+emsPY3JDtw0rBwsrXSsNwo3SDlovJDd02yhxsLDStdIg3CjtoN2c3dBto8TBwkrXysykjdIO2s15mGdcHfROBNZ/nhvKStfKfEVFafMav/t1c85rXD/98B0HosG9Fb1YIMmGn3tD2UgpqzfUKE3x4XXKx5fnhjLrobkLp3jXN25idc7f3DDeYGee9F3y/nsE3H7GK/qdHesjd50165FKiSIvUMBw9nwatLYfGwylBwtJB0paLGmJpJPqpHmNpFskLZL01X7n0ayevO661YbxBjvMDfPWO6U2cEuaBpwL7A8sBW6QNC8ibqlKszXwPmCviLhP0tPKya3Z2hqVHIb1BusFkixP2b2hdgOWRMRtAJIuBg4FbqlKcxxwbkTcBxARd/c9l2Z11FvTepo0lI3bFeO8TrnlK7saahK4s+r3pdm2as8Gni3pWknXSTow70CSjpe0QNKCHuXVbC31qmw++ZqdfLO1kVJ2ySKv1b22pXBdYGtgH2Az4EeSdoiI5Wu8KGI2MBtSb6juZ9VsbYNaZTMoEyra6Cg7WCwFNq/6fTNgWU6a6yJiJXC7pMWk4HFDf7Jo1tigVdl4UJ31QtnVUDcAW0uaKWk94HXAvJo0c4F9ASRtTKqWuq2vuTRrwqBMkTFoEyraaCi1ZBERqySdAMwHpgHnR8QiSacCCyJiXrbv5ZJuAVYDJ0bEn8rLtdnaevVtvp3qpGYmVHQ1lbXKI7jNuqAXU2S0O4VIUV5GaWqSdjhQ1ucR3GY91ovp0etVJ7396zc2rOYqGlQ3ztVUg7RA1rBxsDCro5U2iF5MkdEo0DS6yR228ySnH74jkzMmEKlEUV1qGLR1P/ppnANlp8ruDWU2kFptgzjxgG26Pj16vQF/FZWbXL2p3OtVrdQ77jBOTdKqcQ6UnXLJwixHq99Ai77NtyOvOqlWOze5cZ77yZMkts8lC7Mc7XwD7fZ4i+oBf/VKGO3c5AZ1IGE/9KIEOC4cLMxyDEpVTSUA1evB1O5NbtAGEvbLOAfKThV2nZX0lIh4sE/56Qp3nbVODWL3Unf5tF7raKU8ST8AHgD+SJq3aRnwnYi4vtsZ7RYHC+sG35xt3HQaLI4GLuDxSf8qL/hwRJzatVx2kYOFmVnrOl6DW9IuwCGkif42AnYHNiEtSHRdF/PaFQ4WNo5cErJOdVqyeHNEnFezbTrwaWB6RBxfs+9fI+JzHea5Iw4WNm4GsY3Fhk+n0328p3ZDNl34icCOOenf0Vr2zKxTHplsvdZM19ktJM0FFrL2wkRPk/Sh7GcBzwGe2cX8mVkTPDLZeq3ZcRavJLVZ1ApgVtXPYu2AYmY9NijjQmx0NRssvgcsYM1AIOAfST2lIK1HsQNwaNdyZ2ZN8chk67VmGrh/HxFb5GyfAK6OiD1qtv8mIrbubjZb4wZuG0fuDWWd6rQ31D9FxPk129YFPgU8Mac31Bsj4ssd5rkjDhZmZq3rxjiL55PaLTYHNgT2BP4WeHFE/KSLee0KBwszs9Y1ChaFbRaSXs+aI7grz6cNYqAwM7Pua3ZuqD8Dd5MauO8CLhvEkdsVLlmYmbWu0zaLp0bEAz3JWY84WJiZta7jNoth42BhZta6Tqf7MDOzMedgYWZmhbysqpk9xgP7rB63WZgZ4GnOzW0WZtYET3NujThYmBngac6tMQcLMwPqT2fuac4NHCzMLHPiAdswMX3aGts8zblVuDeUmQE81ojt3lCWp/TeUJIOBD5DWjzpixFxRp10RwDfBHaNiAUFx3RvKDOzFg1sbyhJ04BzgYOA7YEjJW2fk+4pwFuB6/ubQzMzg/LbLHYDlkTEbRHxF+Bi8pdlPQ04E3ikn5kzM7Ok7GAxCdxZ9fvSbNtjJO0MbB4R32l0IEnHS1ogqWEVlZmZta7sBu68urHHGhskrQOcDRxbdKCImA3Mzl7nBgtriqe3MGtO2cFiKWmp1orNgGVVvz8F2AG4RhLA04F5kg4pauQ2K1I7vcXU8hW879KbARwwzGqU2htK0rrAr4H9gCngBuCoiFhUJ/01wLvdG8q6Ya8zrmIqZ3Ty5IwJrj3ppS0fz6UUG3YdrcHdSxGxStIJwHxS19nzI2KRpFOBBRExr8z82Wjr5vQWLqXYqCu7GoqIuBy4vGbbh+qk3acfebLxsOmMidySRSvTW1RKE3nHqUzC52Bho6Ds3lBmpel0eotKaSIvUFR4Ej4bFaWXLMzK0un0FnlTetfyJHw2KhwsbKwdtvNk29VERaUGT8Jno8TVUGZtalRqmJwx4RXmbKSUPpFgL7jr7OAbhW6mXobURk2jrrMOFtZ3o3STHYWgZ1bhYGEDpduD4cysOwZ2inIbT17r2Wz4OFhY33mtZ7Ph42Bhfee1ns2Gj8dZWN95rWez4eMGbjMzA9zAbWZmHXKwMDOzQg4WZmZWyMHCzMwKOViYmVkhBwszMyvkYGFmZoUcLMzMrJCDhZmZFXKwMDOzQg4WZmZWyMHCzMwKOViYmVkhBwszMyvkYGFmZoUcLMzMrJCDhZmZFXKwMDOzQqUHC0kHSlosaYmkk3L2v1PSLZJ+Ken7krYsI59mZuOs1GAhaRpwLnAQsD1wpKTta5ItBHaJiOcClwBn9jeXZmZWdsliN2BJRNwWEX8BLgYOrU4QEVdHxMPZr9cBm/U5j2ZmY6/sYDEJ3Fn1+9JsWz1vAr6Xt0PS8ZIWSFrQxfyZmRmwbsnnV862yE0oHQ3sAuydtz8iZgOzs7S5xzAzs/aUHSyWAptX/b4ZsKw2kaSXAR8A9o6IR/uUNzMzy5RdDXUDsLWkmZLWA14HzKtOIGln4D+AQyLi7hLyaGY29koNFhGxCjgBmA/cCnwjIhZJOlXSIVmys4AnA9+UdKOkeXUOZ2ZmPaKI0avelxSjeF1mZr0kiYjIa0suvRrKzMyGgIOFmZkVcrAwM7NCDhZmZlbIwcLMzAo5WJiZWSEHCzMzK+RgYWZmhRwszMyskIOFmZkVcrAwM7NCDhZmZlbIwcLMzAo5WJiZWSEHCzMzK+RgYWZmhRwszMyskIOFmZkVcrAwM7NCDhZmZlZo3bIzYGbtm7twirPmL2bZ8hVsOmOCEw/YhsN2niw7WzaCFBFl56HrJMUoXpdZtbkLp3jfpTezYuXqx7ZNTJ/G6Yfv6IBhbZFERChvn6uhzIbUWfMXrxEoAFasXM1Z8xeXlCMbZQ4WZkNq2fIVLW0364SDhdmQ2nTGREvbzTrhYGE2pE48YBsmpk9bY9vE9GmceMA2JeXIRpl7Q9nYGLWeQ5W8j9I12eBybygbC+45ZFbMvaFs7LnnkFlnHCxsLLjnkFlnSg8Wkg6UtFjSEkkn5ex/gqSvZ/uvl7RV/3Npw849h8w6U2qwkDQNOBc4CNgeOFLS9jXJ3gTcFxHPAs4GPt7fXNoocM8hs86UXbLYDVgSEbdFxF+Ai4FDa9IcClyQ/XwJsJ+k3AYYs3oO23mS0w/fkckZEwiYnDHhxm2zFpTddXYSuLPq96XA7vXSRMQqSfcDGwH3VCeSdDxwfO+yasPusJ0nHRzM2lR2sMgrIdT2eW0mDRExG5gNqets51kzM7OKsquhlgKbV/2+GbCsXhpJ6wIbAPf2JXdmZgaUHyxuALaWNFPSesDrgHk1aeYBb8h+PgK4yiPuzMz6q9RqqKwN4gRgPjANOD8iFkk6FVgQEfOALwFfkbSEVKJ4XXk5NjMbT57uw8zMAE/3YWZmHXKwMDOzQg4WZmZWqOxxFj3jQd5mZt0zkg3cw07SgojYpex8dMMoXQv4egadr6d3XA1lZmaFHCzMzKyQg8Vgml12BrpolK4FfD2DztfTI26zMDOzQi5ZmJlZIQcLMzMr5GBREknnS7pb0v/W2S9J52Rrj/9S0vP7ncdWNHE9r8+u45eSfiJpp37nsRVF11OVbldJqyUd0a+8taqZa5G0j6QbJS2S9IN+5q9VTfytbSDpMkk3Zdfzxn7nsRWSNpd0taRbs/y+LSdN6fcDB4vyzAEObLD/IGDr7HE88Lk+5KkTc2h8PbcDe0fEc4HTGKCGuzrm0Ph6KmvIf5w0a/Igm0ODa5E0AzgPOCQingO8uk/5atccGn82/wbcEhE7AfsAn8yWQBhUq4B3RcR2wB7Av0naviZN6fcDB4uSRMQPabyI06HAf0ZyHTBD0jP6k7vWFV1PRPwkIu7Lfr2OtNDVwGri8wF4C/At4O7e56h9TVzLUcClEfH7LP2wX08AT1GaxuHJWdpV/chbOyLiroj4Rfbzg8CtpOWkq5V+P3CwGFx565OPygLSbwK+V3YmOiFpEvgH4PNl56ULng1sKOkaST+X9I9lZ6hDnwW2I626eTPwtoj4a7lZao6krYCdgetrdpV+PxjZuaFGQFNrjw8bSfuSgsWLys5Lhz4NvDciVo/APGTrAi8A9gMmgJ9Kui4ifl1uttp2AHAj8FLgmcB/S/pRRDxQbrYak/RkUkn17Tl5Lf1+4GAxuJpZn3yoSHou8EXgoIj4U9n56dAuwMVZoNgYOFjSqoiYW2622rIUuCciHgIekvRDYCdgWIPFG4EzshXQlki6HdgW+Fm52apP0nRSoLgoIi7NSVL6/cDVUINrHvCPWS+IPYD7I+KusjPVLklbAJcCxwzxN9bHRMTMiFotIx0AAAjISURBVNgqIrYCLgHePKSBAuC/gBdLWlfS+sDupHrzYfV7UikJSX8LbAPcVmqOGsjaVr4E3BoRn6qTrPT7gUsWJZH0NVJPjY0lLQVOAaYDRMTngcuBg4ElwMOkb0sDq4nr+RCwEXBe9m181aDMppmniesZGkXXEhG3SroC+CXwV+CLEdGwy3CZmvhsTgPmSLqZVH3z3oi4p6TsNmMv4BjgZkk3ZtveD2wBg3M/8HQfZmZWyNVQZmZWyMHCzMwKOViYmVkhBwszMyvkYGFmZoUcLMyGTDYtx0B3Y5S0laSQNKfsvFh3OFhYX0k6NruJVD8elfQ7SRfVm7pc0qyq9G8uOPZHGpx/fUnLs3Rf7dZ1mY06D8qzstwEVEY8P5U0MOko4FWS9ouIaxu8dpakC9uc6+e1wAakeXUOl7TRCEw9YtZzLllYWW6MiFnZ450RsTvwH8ATgLolA9II1k2A97Z53uNJo5Q/kZ3rDW0ex2ysOFhYQ9V1z5KeKekSSX+S9KCkKyXtkKXbRNJsSXdJekTSDdkMs634Uva8a4M0/06aQO0dklpaEyPL6x7A90mLFv0FOK7FPJK9FyFppqQTJN2SXfMdkt6fzfWDpFdL+pmkh5RWdvuspCfWOeZ+kq6QdG92rF9LOkPSBi3m7QBJl0u6J6ve+62ks5QWOMpLv5nSCmy/yc57b5bnD9akC0nXFLwfWzWZx/UlvU9pZb6HJP1Z0k8lHZmTVpLeoLS64h+zPN4pab6k1zZzPusOBwtr1lakOfb/lrRS2ZXAy4BrJG1NWtBoV+DrwDdIs5Z+T2kCwWZVpmFe2SDNw8AHSVNpf7SFY0MqVQDMyaqevgNsK+nFLR6n4hPAh0mzmX6eVGL5KHCKpLcCF5BKQp8H/o+0gtvZtQeR9C/Af5Oq4uaSpj+/l1R6+km9G33OcT4EXEGaCPC7wDnZ+d8NXCvpqTXpdyFVB76FFIA/A1wEPAjMau4taE12LT8GPgasBs4nvU+bAF/NaW/6KOnv7emkv6tPAf9DWsth0Ff0Gy0R4YcfdR+kIBHZ4wM1+z6Ybb+XdENcp2rfMdm+s2tec2y2fU7Oub6Q7bssZ9+sbN8/k77k3ES62Twv59gfyXn9E7N8Lgcmsm2vzNJ/pcX3ZE72ujuAyartM4B7gIeAPwLbVe17AnAL8CjwtKrtW2bbHgC2rTnPedl5Ztdsvyb9666xbd8s7U+AGXXe87Ortq1HWuo2gKNyrnHzmt8DuKbg/dgq5+9mTp2078n5fK4gBdzqz/RPpOm5188578Zl/3+M08MlC2vWHcAZNdsuyJ6fAJwYa65G9lXSUpbPq3O852U9nGZJ+pSkG0iBYBnwrkYZyc7zHlLQ+EST+X8NsCHw9YhYkW37HvAH4AhJGzZ5nGqnRcRUVb6Wk6aSXh/4XETcWrXvUVKpaz3SKm4VR2fbPhsRv6o5/gdI3/KPkfSEgry8NXs+LsvHYyJiDmkxoNdXbX4l6YY+LyLW6hUWEXfWbuuUpI1I17sgIs6sOd8jpJKUSB0dqq0kfTGozeMgzyQ7ctwbypp1Y0TU/sNWFl/5daS1gx8TaQW5P1B/re2dske13wMvjmwt6EYiYr6kK4GXSzo4Ii4veEmlbeLLVcdYJeki4J2kktA5ReetsSBnW+U9+XnOvkpgqX5Pnp89X1WbOCLuk7QQeAlp8Z6bGuRlT9JN9dWS8qpn1gM2qer9tUe2vZ/L2+4KTANC0qyc/dOz5+pgehGpmmyRpG8CPwB+GhH39zKjtjYHC2vWWv+c2c02d19mFY/fAGpdEBHHZo3BTyMttfoR4DJJe0bEw03k6d2kb8xnSppfL5Gk7UjLuP4q0mL31b5MChbH0XqwyLvuVU3sq35PKg3Y9RayqWwvarfYiPT/fEpBuieTqnYqx5tqkLbbNsqed6VxJ4YnV/38DuC3wD8BJ2WPVZIuB94VEUt6kVFbm6uhrFSR/CEiPgZ8EngujbvOVr/2ZlJV2HNIN5N6Kg3b26pmQCBwc7ZvB0kvbO8qOlIJKk+vs/8ZNekaHee+iFDB43dZ+kpV1WST+Qzqf7lsqgGex6/h7II8PtaLLiJWR8RnImInUueKVwHfBg4Brmiies66xMHCBsmppIbhEyTNbPI1J5N6SJ0KPKl2Z3YzOYbUcHo+qXtu7aNSKmm5G20XLMye96ndkfUceh7wCMXLnF4HbCjpOU2et1LCOqjJ9Pex5hrQlTxOo367VK2fkT6HtnqfRcTdEXFpRLyGVG33TGCHdo5lrXOwsIGRtXt8nFRNM6vJ1ywjlUieDrw9J8mrSNUf8yPiTRHxz7UPUhfMh4DXtDquoQsuJLU1vEXSs2r2nUYa3X5h1kDeSKVL7hckbVq7U9KTlNZurriM1GnhkDrjG2pLHD8DtpD08prtJ5N6dBWKiLtJbRC7SPqgpLVKKkpjeWZmPz8hG3+imjTTgb/Jfm2mutK6wG0WNmjOI/WGOlrSxyPiliZecyapqqn2ZguPV0F9sd6LI+LBrPH0WFJvnXNbynEHIuIOSW/PzvkLSd8gla72JjVa/4omRqtHxPclnQScDvwmq9O/nVT/v2V2vB8DB2bp/5I1hF9JGt/wL6TSxhNJDcz7seb94RPAAcB/Sfo6qRvyC4GZpK68+zR5yScAW5NKgsdI+jGpR9qm2Xl3BY7M8j5BGlNxh6Trgd9l+ds/SzuvuseZ9ZZLFjZQsm6tp5P+Nk9r8jV/JqdhNxssuDdwN+mbdCNfyJ77XhUVEeeRbsTXkUpC7yQ1+p8F7BkR9zZ5nI+Tek59lzTA7+2kUtMkMJtUCqhOv4BUhfQ5UkCp9AqbQc37GRHfBw4DFgGvI02TcgewG+km3uy1PkD6TN5CGpNSud59Sd2E30EaoAiptPdeUsB8IfA2UrfaB4B/xYPy+koRAz3TsZmZDQCXLMzMrJCDhZmZFXKwMDOzQg4WZmZWyMHCzMwKOViYmVkhBwszMyvkYGFmZoUcLMzMrND/BwUZ0oYHr+fMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.log10(mrna_per_event_song.loc['GNAS_6', song_iPSC]+1), song_PSI.loc['GNAS_6', song_iPSC])\n",
    "plt.title('GNAS_6 in iPS cells', fontsize=20)\n",
    "plt.ylabel('$\\Psi$', fontsize=20)\n",
    "plt.xlabel('mRNA molecules', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table of bimodality by given cell type\n",
    "\n",
    "Comparison of the intermediate exons that are bimodal by the quartile definition, before and after filtering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_dset(PSI_tab, subpop_list, mrna_counts, mrna_per_event, read_counts, coverage_tab, psi_int = 0.1, mrna_min = 10,\n",
    "             cell_min=0.5):\n",
    "    total = []\n",
    "    bimodal_total = []\n",
    "    filtered = []\n",
    "    bimodal_filtered = []\n",
    "    pvals = []\n",
    "    for subpop in subpop_list:\n",
    "        proport = get_proportions(PSI_tab, subpop, mrna_counts, mrna_per_event, read_counts, coverage_tab, \n",
    "                                  psi_int = psi_int, mrna_min = mrna_min, cell_min=cell_min)\n",
    "        \n",
    "        total.append(proport[0])\n",
    "        bimodal_total.append(proport[2])\n",
    "        filtered.append(proport[1])\n",
    "        bimodal_filtered.append(proport[4]) # filtered\n",
    "        #bimodal_filtered.append(proport[3]) # unfiltered\n",
    "        pvals.append(proport[7])\n",
    "        \n",
    "    return total, bimodal_total, filtered, bimodal_filtered, pvals\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total intermediate exons: 452\n",
      "Total intermediate exons that are bimodal by quartile definition: 108\n",
      "#################\n",
      "Selected exons: 97\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 6.425484125955408e-14\n",
      "Total intermediate exons: 381\n",
      "Total intermediate exons that are bimodal by quartile definition: 110\n",
      "#################\n",
      "Selected exons: 81\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 1.5871904004859756e-14\n",
      "Total intermediate exons: 476\n",
      "Total intermediate exons that are bimodal by quartile definition: 120\n",
      "#################\n",
      "Selected exons: 98\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 7.747197892632424e-15\n",
      "Total intermediate exons: 457\n",
      "Total intermediate exons that are bimodal by quartile definition: 173\n",
      "#################\n",
      "Selected exons: 43\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 3\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 1.7139526385648574e-06\n",
      "Total intermediate exons: 286\n",
      "Total intermediate exons that are bimodal by quartile definition: 77\n",
      "#################\n",
      "Selected exons: 33\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 1.5129644902405944e-05\n",
      "Total intermediate exons: 291\n",
      "Total intermediate exons that are bimodal by quartile definition: 78\n",
      "#################\n",
      "Selected exons: 36\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 5.495198262731851e-06\n",
      "Total intermediate exons: 342\n",
      "Total intermediate exons that are bimodal by quartile definition: 141\n",
      "#################\n",
      "Selected exons: 32\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 1.3522300405590075e-08\n",
      "Total intermediate exons: 254\n",
      "Total intermediate exons that are bimodal by quartile definition: 104\n",
      "#################\n",
      "Selected exons: 30\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 2\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 1.2457726158663918e-05\n",
      "Total intermediate exons: 213\n",
      "Total intermediate exons that are bimodal by quartile definition: 90\n",
      "#################\n",
      "Selected exons: 22\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 1\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 4.823278633114043e-05\n",
      "Total intermediate exons: 284\n",
      "Total intermediate exons that are bimodal by quartile definition: 102\n",
      "#################\n",
      "Selected exons: 37\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 1\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 4.328263870020988e-07\n",
      "Total intermediate exons: 616\n",
      "Total intermediate exons that are bimodal by quartile definition: 269\n",
      "#################\n",
      "Selected exons: 55\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 2.5476050557883045e-15\n",
      "Total intermediate exons: 220\n",
      "Total intermediate exons that are bimodal by quartile definition: 105\n",
      "#################\n",
      "Selected exons: 28\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 2\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 1.2111076400049347e-06\n",
      "Total intermediate exons: 188\n",
      "Total intermediate exons that are bimodal by quartile definition: 85\n",
      "#################\n",
      "Selected exons: 20\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 4\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 0.01345020242302973\n",
      "Total intermediate exons: 149\n",
      "Total intermediate exons that are bimodal by quartile definition: 51\n",
      "#################\n",
      "Selected exons: 27\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 1\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 5.78332849530742e-05\n",
      "Total intermediate exons: 7\n",
      "Total intermediate exons that are bimodal by quartile definition: 1\n",
      "#################\n",
      "Selected exons: 3\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 0.5714285714285712\n"
     ]
    }
   ],
   "source": [
    "total = [] \n",
    "bimodal_total = [] \n",
    "filtered = [] \n",
    "bimodal_filtered = [] \n",
    "pvals =[] \n",
    "\n",
    "\n",
    "cells_in_cluster = []\n",
    "\n",
    "cells_in_cluster.extend([len(x) for x in [chen_ES2i, chen_ES, chen_Epi, chen_MN]])\n",
    "cells_in_cluster.extend([len(x) for x in [lescroart_E6, lescroart_E7]])\n",
    "cells_in_cluster.extend([len(x) for x in [trapnell_M00, trapnell_M24, trapnell_M48, trapnell_M72]])\n",
    "cells_in_cluster.extend([len(x) for x in [song_iPSC, song_NPC, song_MN]])\n",
    "cells_in_cluster.append(len(shalek_PSI.columns))\n",
    "cells_in_cluster.append(len(das_PSI.columns))\n",
    "\n",
    "\n",
    "chen_processed = test_dset(chen_PSI, [chen_ES2i, chen_ES, chen_Epi, chen_MN], chen_mrna_counts, mrna_per_event_chen, \n",
    "          chen_read_counts, chen_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "lescroart_processed = test_dset(lescroart_PSI, [lescroart_E6, lescroart_E7], lescroart_mrna_counts, mrna_per_event_lescroart, \n",
    "          lescroart_read_counts, lescroart_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "trapnell_processed = test_dset(trapnell_PSI, [trapnell_M00, trapnell_M24, trapnell_M48, trapnell_M72], \n",
    "                               trapnell_mrna_counts, mrna_per_event_trapnell, trapnell_read_counts, \n",
    "                               trapnell_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "song_processed = test_dset(song_PSI, [song_iPSC, song_NPC, song_MN], song_mrna_counts, mrna_per_event_song, \n",
    "          song_read_counts, song_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "shalek_processed = test_dset(shalek_PSI, [shalek_PSI.columns], shalek_mrna_counts, mrna_per_event_shalek, \n",
    "          shalek_read_counts, shalek_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "das_processed = test_dset(das_PSI, [das_PSI.columns], das_mrna_counts, mrna_per_event_das, \n",
    "          das_read_counts, das_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "\n",
    "total.extend(chen_processed[0])\n",
    "total.extend(lescroart_processed[0])\n",
    "total.extend(trapnell_processed[0])\n",
    "total.extend(song_processed[0])\n",
    "total.extend(shalek_processed[0])\n",
    "total.extend(das_processed[0])\n",
    "\n",
    "bimodal_total.extend(chen_processed[1])\n",
    "bimodal_total.extend(lescroart_processed[1])\n",
    "bimodal_total.extend(trapnell_processed[1])\n",
    "bimodal_total.extend(song_processed[1])\n",
    "bimodal_total.extend(shalek_processed[1])\n",
    "bimodal_total.extend(das_processed[1])\n",
    "\n",
    "filtered.extend(chen_processed[2])\n",
    "filtered.extend(lescroart_processed[2])\n",
    "filtered.extend(trapnell_processed[2])\n",
    "filtered.extend(song_processed[2])\n",
    "filtered.extend(shalek_processed[2])\n",
    "filtered.extend(das_processed[2])\n",
    "\n",
    "bimodal_filtered.extend(chen_processed[3])\n",
    "bimodal_filtered.extend(lescroart_processed[3])\n",
    "bimodal_filtered.extend(trapnell_processed[3])\n",
    "bimodal_filtered.extend(song_processed[3])\n",
    "bimodal_filtered.extend(shalek_processed[3])\n",
    "bimodal_filtered.extend(das_processed[3])\n",
    "\n",
    "pvals.extend(chen_processed[4])\n",
    "pvals.extend(lescroart_processed[4])\n",
    "pvals.extend(trapnell_processed[4])\n",
    "pvals.extend(song_processed[4])\n",
    "pvals.extend(shalek_processed[4])\n",
    "pvals.extend(das_processed[4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "pval_adj = multipletests(pvals, method='fdr_bh')[1]\n",
    "\n",
    "cell_type = ['mES2i', 'mES', 'Epi', 'Motor neuron', 'Heart E6.75', 'Heart E7.25', \n",
    "             'Myoblast 00h', 'Myoblast 24h', 'Myoblast 48h', 'Myoblast 72h',\n",
    "             'iPSC', 'NPC', 'Motor neuron', 'BMDC', 'Olfactory neurons']\n",
    "dataset = ['Chen']*4 + ['Lescroart']*2 + ['Trapnell']*4 + ['Song']*3 + ['Shalek'] + ['Fletcher']\n",
    "organism = ['Mouse']*6+['Human']*7+['Mouse']*2\n",
    "\n",
    "bimodality_table = pd.DataFrame()\n",
    "bimodality_table['dataset'] = dataset\n",
    "bimodality_table['organism'] = organism\n",
    "bimodality_table['cell_type'] = cell_type\n",
    "bimodality_table['cells_in_cluster'] = cells_in_cluster\n",
    "bimodality_table['total_exons'] = total\n",
    "bimodality_table['bimodal_exons'] = bimodal_total\n",
    "bimodality_table['bimodal_percent'] = [str(round(x*100, 2))+'%' for x in np.array(bimodal_total)/np.array(total)]\n",
    "bimodality_table['selected_exons'] = filtered\n",
    "bimodality_table['selected_bimodal'] = bimodal_filtered\n",
    "bimodality_table['bimodal_percent_selected'] = [str(round(x*100, 2))+'%' for x in np.array(bimodal_filtered)/np.array(filtered)]\n",
    "bimodality_table['p-val'] = pvals\n",
    "bimodality_table['p-val (adj)'] = pval_adj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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>dataset</th>\n",
       "      <th>organism</th>\n",
       "      <th>cell_type</th>\n",
       "      <th>cells_in_cluster</th>\n",
       "      <th>total_exons</th>\n",
       "      <th>bimodal_exons</th>\n",
       "      <th>bimodal_percent</th>\n",
       "      <th>selected_exons</th>\n",
       "      <th>selected_bimodal</th>\n",
       "      <th>bimodal_percent_selected</th>\n",
       "      <th>p-val</th>\n",
       "      <th>p-val (adj)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Chen</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>mES2i</td>\n",
       "      <td>111</td>\n",
       "      <td>452</td>\n",
       "      <td>108</td>\n",
       "      <td>23.89%</td>\n",
       "      <td>97</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>6.425484e-14</td>\n",
       "      <td>2.409557e-13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Chen</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>mES</td>\n",
       "      <td>120</td>\n",
       "      <td>381</td>\n",
       "      <td>110</td>\n",
       "      <td>28.87%</td>\n",
       "      <td>81</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>1.587190e-14</td>\n",
       "      <td>7.935952e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Chen</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Epi</td>\n",
       "      <td>183</td>\n",
       "      <td>476</td>\n",
       "      <td>120</td>\n",
       "      <td>25.21%</td>\n",
       "      <td>98</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>7.747198e-15</td>\n",
       "      <td>5.810398e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Chen</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Motor neuron</td>\n",
       "      <td>74</td>\n",
       "      <td>457</td>\n",
       "      <td>173</td>\n",
       "      <td>37.86%</td>\n",
       "      <td>43</td>\n",
       "      <td>3</td>\n",
       "      <td>6.98%</td>\n",
       "      <td>1.713953e-06</td>\n",
       "      <td>3.213661e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Lescroart</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Heart E6.75</td>\n",
       "      <td>172</td>\n",
       "      <td>286</td>\n",
       "      <td>77</td>\n",
       "      <td>26.92%</td>\n",
       "      <td>33</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>1.512964e-05</td>\n",
       "      <td>2.063133e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Lescroart</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Heart E7.25</td>\n",
       "      <td>341</td>\n",
       "      <td>291</td>\n",
       "      <td>78</td>\n",
       "      <td>26.8%</td>\n",
       "      <td>36</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>5.495198e-06</td>\n",
       "      <td>9.158664e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Trapnell</td>\n",
       "      <td>Human</td>\n",
       "      <td>Myoblast 00h</td>\n",
       "      <td>58</td>\n",
       "      <td>342</td>\n",
       "      <td>141</td>\n",
       "      <td>41.23%</td>\n",
       "      <td>32</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>1.352230e-08</td>\n",
       "      <td>4.056690e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Trapnell</td>\n",
       "      <td>Human</td>\n",
       "      <td>Myoblast 24h</td>\n",
       "      <td>70</td>\n",
       "      <td>254</td>\n",
       "      <td>104</td>\n",
       "      <td>40.94%</td>\n",
       "      <td>30</td>\n",
       "      <td>2</td>\n",
       "      <td>6.67%</td>\n",
       "      <td>1.245773e-05</td>\n",
       "      <td>1.868659e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Trapnell</td>\n",
       "      <td>Human</td>\n",
       "      <td>Myoblast 48h</td>\n",
       "      <td>70</td>\n",
       "      <td>213</td>\n",
       "      <td>90</td>\n",
       "      <td>42.25%</td>\n",
       "      <td>22</td>\n",
       "      <td>1</td>\n",
       "      <td>4.55%</td>\n",
       "      <td>4.823279e-05</td>\n",
       "      <td>6.029098e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Trapnell</td>\n",
       "      <td>Human</td>\n",
       "      <td>Myoblast 72h</td>\n",
       "      <td>33</td>\n",
       "      <td>284</td>\n",
       "      <td>102</td>\n",
       "      <td>35.92%</td>\n",
       "      <td>37</td>\n",
       "      <td>1</td>\n",
       "      <td>2.7%</td>\n",
       "      <td>4.328264e-07</td>\n",
       "      <td>1.082066e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>Song</td>\n",
       "      <td>Human</td>\n",
       "      <td>iPSC</td>\n",
       "      <td>62</td>\n",
       "      <td>616</td>\n",
       "      <td>269</td>\n",
       "      <td>43.67%</td>\n",
       "      <td>55</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>2.547605e-15</td>\n",
       "      <td>3.821408e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>Song</td>\n",
       "      <td>Human</td>\n",
       "      <td>NPC</td>\n",
       "      <td>70</td>\n",
       "      <td>220</td>\n",
       "      <td>105</td>\n",
       "      <td>47.73%</td>\n",
       "      <td>28</td>\n",
       "      <td>2</td>\n",
       "      <td>7.14%</td>\n",
       "      <td>1.211108e-06</td>\n",
       "      <td>2.595231e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>Song</td>\n",
       "      <td>Human</td>\n",
       "      <td>Motor neuron</td>\n",
       "      <td>70</td>\n",
       "      <td>188</td>\n",
       "      <td>85</td>\n",
       "      <td>45.21%</td>\n",
       "      <td>20</td>\n",
       "      <td>4</td>\n",
       "      <td>20.0%</td>\n",
       "      <td>1.345020e-02</td>\n",
       "      <td>1.441093e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>Shalek</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>BMDC</td>\n",
       "      <td>13</td>\n",
       "      <td>149</td>\n",
       "      <td>51</td>\n",
       "      <td>34.23%</td>\n",
       "      <td>27</td>\n",
       "      <td>1</td>\n",
       "      <td>3.7%</td>\n",
       "      <td>5.783328e-05</td>\n",
       "      <td>6.673071e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>Fletcher</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Olfactory neurons</td>\n",
       "      <td>812</td>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "      <td>14.29%</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>5.714286e-01</td>\n",
       "      <td>5.714286e-01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      dataset organism          cell_type  cells_in_cluster  total_exons  \\\n",
       "0        Chen    Mouse              mES2i               111          452   \n",
       "1        Chen    Mouse                mES               120          381   \n",
       "2        Chen    Mouse                Epi               183          476   \n",
       "3        Chen    Mouse       Motor neuron                74          457   \n",
       "4   Lescroart    Mouse        Heart E6.75               172          286   \n",
       "5   Lescroart    Mouse        Heart E7.25               341          291   \n",
       "6    Trapnell    Human       Myoblast 00h                58          342   \n",
       "7    Trapnell    Human       Myoblast 24h                70          254   \n",
       "8    Trapnell    Human       Myoblast 48h                70          213   \n",
       "9    Trapnell    Human       Myoblast 72h                33          284   \n",
       "10       Song    Human               iPSC                62          616   \n",
       "11       Song    Human                NPC                70          220   \n",
       "12       Song    Human       Motor neuron                70          188   \n",
       "13     Shalek    Mouse               BMDC                13          149   \n",
       "14   Fletcher    Mouse  Olfactory neurons               812            7   \n",
       "\n",
       "    bimodal_exons bimodal_percent  selected_exons  selected_bimodal  \\\n",
       "0             108          23.89%              97                 0   \n",
       "1             110          28.87%              81                 0   \n",
       "2             120          25.21%              98                 0   \n",
       "3             173          37.86%              43                 3   \n",
       "4              77          26.92%              33                 0   \n",
       "5              78           26.8%              36                 0   \n",
       "6             141          41.23%              32                 0   \n",
       "7             104          40.94%              30                 2   \n",
       "8              90          42.25%              22                 1   \n",
       "9             102          35.92%              37                 1   \n",
       "10            269          43.67%              55                 0   \n",
       "11            105          47.73%              28                 2   \n",
       "12             85          45.21%              20                 4   \n",
       "13             51          34.23%              27                 1   \n",
       "14              1          14.29%               3                 0   \n",
       "\n",
       "   bimodal_percent_selected         p-val   p-val (adj)  \n",
       "0                      0.0%  6.425484e-14  2.409557e-13  \n",
       "1                      0.0%  1.587190e-14  7.935952e-14  \n",
       "2                      0.0%  7.747198e-15  5.810398e-14  \n",
       "3                     6.98%  1.713953e-06  3.213661e-06  \n",
       "4                      0.0%  1.512964e-05  2.063133e-05  \n",
       "5                      0.0%  5.495198e-06  9.158664e-06  \n",
       "6                      0.0%  1.352230e-08  4.056690e-08  \n",
       "7                     6.67%  1.245773e-05  1.868659e-05  \n",
       "8                     4.55%  4.823279e-05  6.029098e-05  \n",
       "9                      2.7%  4.328264e-07  1.082066e-06  \n",
       "10                     0.0%  2.547605e-15  3.821408e-14  \n",
       "11                    7.14%  1.211108e-06  2.595231e-06  \n",
       "12                    20.0%  1.345020e-02  1.441093e-02  \n",
       "13                     3.7%  5.783328e-05  6.673071e-05  \n",
       "14                     0.0%  5.714286e-01  5.714286e-01  "
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bimodality_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "bimodality_table.to_csv('selected_int_exons.csv', index=False, header=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Table of bimodality by given agglomerative clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total intermediate exons: 446\n",
      "Total intermediate exons that are bimodal by quartile definition: 118\n",
      "#################\n",
      "Selected exons: 94\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "0%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 4.089694836903887e-15\n",
      "Total intermediate exons: 402\n",
      "Total intermediate exons that are bimodal by quartile definition: 107\n",
      "#################\n",
      "Selected exons: 98\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 1\n",
      "1%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 1.7433419418207195e-14\n",
      "Total intermediate exons: 516\n",
      "Total intermediate exons that are bimodal by quartile definition: 136\n",
      "#################\n",
      "Selected exons: 76\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 1\n",
      "1%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 2.9693152867517824e-10\n",
      "Total intermediate exons: 364\n",
      "Total intermediate exons that are bimodal by quartile definition: 117\n",
      "#################\n",
      "Selected exons: 43\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "0%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 1.559634039870456e-08\n",
      "Total intermediate exons: 517\n",
      "Total intermediate exons that are bimodal by quartile definition: 146\n",
      "#################\n",
      "Selected exons: 61\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "0%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 3.472852039246133e-10\n",
      "Total intermediate exons: 286\n",
      "Total intermediate exons that are bimodal by quartile definition: 77\n",
      "#################\n",
      "Selected exons: 33\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "0%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 1.5129644902405944e-05\n",
      "Total intermediate exons: 291\n",
      "Total intermediate exons that are bimodal by quartile definition: 78\n",
      "#################\n",
      "Selected exons: 36\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "0%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 5.495198262731851e-06\n",
      "Total intermediate exons: 400\n",
      "Total intermediate exons that are bimodal by quartile definition: 142\n",
      "#################\n",
      "Selected exons: 41\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "0%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 4.519739931000002e-09\n",
      "Total intermediate exons: 251\n",
      "Total intermediate exons that are bimodal by quartile definition: 97\n",
      "#################\n",
      "Selected exons: 27\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "0%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 7.07430090611066e-07\n",
      "Total intermediate exons: 242\n",
      "Total intermediate exons that are bimodal by quartile definition: 97\n",
      "#################\n",
      "Selected exons: 32\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 1\n",
      "3%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 4.54160727484006e-07\n",
      "Total intermediate exons: 252\n",
      "Total intermediate exons that are bimodal by quartile definition: 101\n",
      "#################\n",
      "Selected exons: 24\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 1\n",
      "4%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 4.1129138247006845e-05\n",
      "Total intermediate exons: 616\n",
      "Total intermediate exons that are bimodal by quartile definition: 269\n",
      "#################\n",
      "Selected exons: 55\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "0%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 2.5476050557883045e-15\n",
      "Total intermediate exons: 212\n",
      "Total intermediate exons that are bimodal by quartile definition: 92\n",
      "#################\n",
      "Selected exons: 28\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 1\n",
      "4%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 7.265806571510512e-07\n",
      "Total intermediate exons: 168\n",
      "Total intermediate exons that are bimodal by quartile definition: 82\n",
      "#################\n",
      "Selected exons: 19\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 4\n",
      "21%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 0.008825301681599248\n",
      "Total intermediate exons: 149\n",
      "Total intermediate exons that are bimodal by quartile definition: 51\n",
      "#################\n",
      "Selected exons: 27\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 1\n",
      "4%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 5.78332849530742e-05\n",
      "Total intermediate exons: 7\n",
      "Total intermediate exons that are bimodal by quartile definition: 1\n",
      "#################\n",
      "Selected exons: 3\n",
      "Selected exons that are bimodal by Song et al. definition: 0\n",
      "0%\n",
      "Selected exons that are bimodal by quartile definition: 0\n",
      "0%\n",
      "#################\n",
      "Depletion of quartile bimodality in selected exon (p-value): 0.5714285714285712\n"
     ]
    }
   ],
   "source": [
    "total = [] \n",
    "bimodal_total = [] \n",
    "filtered = [] \n",
    "bimodal_filtered = [] \n",
    "pvals =[] \n",
    "\n",
    "\n",
    "cells_in_cluster = []\n",
    "\n",
    "cells_in_cluster.extend([len(x) for x in [chen_clust_filter[x][0].columns for x in range(len(chen_clust_filter))]])\n",
    "cells_in_cluster.extend([len(x) for x in [lescroart_E6, lescroart_E7]])\n",
    "cells_in_cluster.extend([len(x) for x in [trapnell_clust_filter[x][0].columns for x in range(len(trapnell_clust_filter))]])\n",
    "cells_in_cluster.extend([len(x) for x in [song_clust_filter[x][0].columns for x in range(len(song_clust_filter))]])\n",
    "cells_in_cluster.append(len(shalek_PSI.columns))\n",
    "cells_in_cluster.append(len(das_PSI.columns))\n",
    "\n",
    "\n",
    "chen_processed = test_dset(chen_PSI, [chen_clust_filter[x][0].columns for x in range(len(chen_clust_filter))], \n",
    "                           chen_mrna_counts, mrna_per_event_chen, \n",
    "          chen_read_counts, chen_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "lescroart_processed = test_dset(lescroart_PSI, [lescroart_E6, lescroart_E7], lescroart_mrna_counts, mrna_per_event_lescroart, \n",
    "          lescroart_read_counts, lescroart_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "trapnell_processed = test_dset(trapnell_PSI, [trapnell_clust_filter[x][0].columns for x in range(len(trapnell_clust_filter))], \n",
    "                               trapnell_mrna_counts, mrna_per_event_trapnell, trapnell_read_counts, \n",
    "                               trapnell_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "song_processed = test_dset(song_PSI, [song_clust_filter[x][0].columns for x in range(len(song_clust_filter))], \n",
    "                           song_mrna_counts, mrna_per_event_song, \n",
    "          song_read_counts, song_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "shalek_processed = test_dset(shalek_PSI, [shalek_PSI.columns], shalek_mrna_counts, mrna_per_event_shalek, \n",
    "          shalek_read_counts, shalek_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "das_processed = test_dset(das_PSI, [das_PSI.columns], das_mrna_counts, mrna_per_event_das, \n",
    "          das_read_counts, das_coverage_tab, psi_int = 0.2, mrna_min = 10)\n",
    "\n",
    "\n",
    "total.extend(chen_processed[0])\n",
    "total.extend(lescroart_processed[0])\n",
    "total.extend(trapnell_processed[0])\n",
    "total.extend(song_processed[0])\n",
    "total.extend(shalek_processed[0])\n",
    "total.extend(das_processed[0])\n",
    "\n",
    "bimodal_total.extend(chen_processed[1])\n",
    "bimodal_total.extend(lescroart_processed[1])\n",
    "bimodal_total.extend(trapnell_processed[1])\n",
    "bimodal_total.extend(song_processed[1])\n",
    "bimodal_total.extend(shalek_processed[1])\n",
    "bimodal_total.extend(das_processed[1])\n",
    "\n",
    "filtered.extend(chen_processed[2])\n",
    "filtered.extend(lescroart_processed[2])\n",
    "filtered.extend(trapnell_processed[2])\n",
    "filtered.extend(song_processed[2])\n",
    "filtered.extend(shalek_processed[2])\n",
    "filtered.extend(das_processed[2])\n",
    "\n",
    "bimodal_filtered.extend(chen_processed[3])\n",
    "bimodal_filtered.extend(lescroart_processed[3])\n",
    "bimodal_filtered.extend(trapnell_processed[3])\n",
    "bimodal_filtered.extend(song_processed[3])\n",
    "bimodal_filtered.extend(shalek_processed[3])\n",
    "bimodal_filtered.extend(das_processed[3])\n",
    "\n",
    "pvals.extend(chen_processed[4])\n",
    "pvals.extend(lescroart_processed[4])\n",
    "pvals.extend(trapnell_processed[4])\n",
    "pvals.extend(song_processed[4])\n",
    "pvals.extend(shalek_processed[4])\n",
    "pvals.extend(das_processed[4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "pval_adj = multipletests(pvals, method='fdr_bh')[1]\n",
    "\n",
    "cell_type = ['ES', 'Epi, early', 'Epi, late', 'Neuron, early', 'Neuron, late', 'Heart E6.75', 'Heart E7.25', \n",
    "             'Myoblast 00h', 'Myoblast 24h', 'Myoblast 48h', 'Myoblast 72h',\n",
    "             'iPSC', 'NPC', 'Motor neuron', 'BMDC', 'Olfactory neurons']\n",
    "dataset = ['Chen']*5 + ['Lescroart']*2 + ['Trapnell']*4 + ['Song']*3 + ['Shalek'] + ['Fletcher']\n",
    "organism = ['Mouse']*7+['Human']*7+['Mouse']*2\n",
    "\n",
    "bimodality_table = pd.DataFrame()\n",
    "bimodality_table['dataset'] = dataset\n",
    "bimodality_table['organism'] = organism\n",
    "bimodality_table['cell_type'] = cell_type\n",
    "bimodality_table['cells_in_cluster'] = cells_in_cluster\n",
    "bimodality_table['total_exons'] = total\n",
    "bimodality_table['bimodal_exons'] = bimodal_total\n",
    "bimodality_table['bimodal_percent'] = [str(round(x*100, 2))+'%' for x in np.array(bimodal_total)/np.array(total)]\n",
    "bimodality_table['selected_exons'] = filtered\n",
    "bimodality_table['selected_bimodal'] = bimodal_filtered\n",
    "bimodality_table['bimodal_percent_selected'] = [str(round(x*100, 2))+'%' for x in np.array(bimodal_filtered)/np.array(filtered)]\n",
    "bimodality_table['p-val'] = pvals\n",
    "bimodality_table['p-val (adj)'] = pval_adj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "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>dataset</th>\n",
       "      <th>organism</th>\n",
       "      <th>cell_type</th>\n",
       "      <th>cells_in_cluster</th>\n",
       "      <th>total_exons</th>\n",
       "      <th>bimodal_exons</th>\n",
       "      <th>bimodal_percent</th>\n",
       "      <th>selected_exons</th>\n",
       "      <th>selected_bimodal</th>\n",
       "      <th>bimodal_percent_selected</th>\n",
       "      <th>p-val</th>\n",
       "      <th>p-val (adj)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Chen</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>ES</td>\n",
       "      <td>217</td>\n",
       "      <td>446</td>\n",
       "      <td>118</td>\n",
       "      <td>26.46%</td>\n",
       "      <td>94</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>4.089695e-15</td>\n",
       "      <td>3.271756e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Chen</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Epi, early</td>\n",
       "      <td>98</td>\n",
       "      <td>402</td>\n",
       "      <td>107</td>\n",
       "      <td>26.62%</td>\n",
       "      <td>98</td>\n",
       "      <td>1</td>\n",
       "      <td>1.02%</td>\n",
       "      <td>1.743342e-14</td>\n",
       "      <td>9.297824e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Chen</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Epi, late</td>\n",
       "      <td>104</td>\n",
       "      <td>516</td>\n",
       "      <td>136</td>\n",
       "      <td>26.36%</td>\n",
       "      <td>76</td>\n",
       "      <td>1</td>\n",
       "      <td>1.32%</td>\n",
       "      <td>2.969315e-10</td>\n",
       "      <td>1.111313e-09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Chen</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Neuron, early</td>\n",
       "      <td>47</td>\n",
       "      <td>364</td>\n",
       "      <td>117</td>\n",
       "      <td>32.14%</td>\n",
       "      <td>43</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>1.559634e-08</td>\n",
       "      <td>3.564878e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Chen</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Neuron, late</td>\n",
       "      <td>22</td>\n",
       "      <td>517</td>\n",
       "      <td>146</td>\n",
       "      <td>28.24%</td>\n",
       "      <td>61</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>3.472852e-10</td>\n",
       "      <td>1.111313e-09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Lescroart</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Heart E6.75</td>\n",
       "      <td>172</td>\n",
       "      <td>286</td>\n",
       "      <td>77</td>\n",
       "      <td>26.92%</td>\n",
       "      <td>33</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>1.512964e-05</td>\n",
       "      <td>2.017286e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Lescroart</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Heart E7.25</td>\n",
       "      <td>341</td>\n",
       "      <td>291</td>\n",
       "      <td>78</td>\n",
       "      <td>26.8%</td>\n",
       "      <td>36</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>5.495198e-06</td>\n",
       "      <td>7.993016e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Trapnell</td>\n",
       "      <td>Human</td>\n",
       "      <td>Myoblast 00h</td>\n",
       "      <td>35</td>\n",
       "      <td>400</td>\n",
       "      <td>142</td>\n",
       "      <td>35.5%</td>\n",
       "      <td>41</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>4.519740e-09</td>\n",
       "      <td>1.205264e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Trapnell</td>\n",
       "      <td>Human</td>\n",
       "      <td>Myoblast 24h</td>\n",
       "      <td>89</td>\n",
       "      <td>251</td>\n",
       "      <td>97</td>\n",
       "      <td>38.65%</td>\n",
       "      <td>27</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>7.074301e-07</td>\n",
       "      <td>1.162529e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Trapnell</td>\n",
       "      <td>Human</td>\n",
       "      <td>Myoblast 48h</td>\n",
       "      <td>72</td>\n",
       "      <td>242</td>\n",
       "      <td>97</td>\n",
       "      <td>40.08%</td>\n",
       "      <td>32</td>\n",
       "      <td>1</td>\n",
       "      <td>3.12%</td>\n",
       "      <td>4.541607e-07</td>\n",
       "      <td>9.083215e-07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>Trapnell</td>\n",
       "      <td>Human</td>\n",
       "      <td>Myoblast 72h</td>\n",
       "      <td>35</td>\n",
       "      <td>252</td>\n",
       "      <td>101</td>\n",
       "      <td>40.08%</td>\n",
       "      <td>24</td>\n",
       "      <td>1</td>\n",
       "      <td>4.17%</td>\n",
       "      <td>4.112914e-05</td>\n",
       "      <td>5.062048e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>Song</td>\n",
       "      <td>Human</td>\n",
       "      <td>iPSC</td>\n",
       "      <td>62</td>\n",
       "      <td>616</td>\n",
       "      <td>269</td>\n",
       "      <td>43.67%</td>\n",
       "      <td>55</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>2.547605e-15</td>\n",
       "      <td>3.271756e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>Song</td>\n",
       "      <td>Human</td>\n",
       "      <td>NPC</td>\n",
       "      <td>73</td>\n",
       "      <td>212</td>\n",
       "      <td>92</td>\n",
       "      <td>43.4%</td>\n",
       "      <td>28</td>\n",
       "      <td>1</td>\n",
       "      <td>3.57%</td>\n",
       "      <td>7.265807e-07</td>\n",
       "      <td>1.162529e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>Song</td>\n",
       "      <td>Human</td>\n",
       "      <td>Motor neuron</td>\n",
       "      <td>67</td>\n",
       "      <td>168</td>\n",
       "      <td>82</td>\n",
       "      <td>48.81%</td>\n",
       "      <td>19</td>\n",
       "      <td>4</td>\n",
       "      <td>21.05%</td>\n",
       "      <td>8.825302e-03</td>\n",
       "      <td>9.413655e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>Shalek</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>BMDC</td>\n",
       "      <td>13</td>\n",
       "      <td>149</td>\n",
       "      <td>51</td>\n",
       "      <td>34.23%</td>\n",
       "      <td>27</td>\n",
       "      <td>1</td>\n",
       "      <td>3.7%</td>\n",
       "      <td>5.783328e-05</td>\n",
       "      <td>6.609518e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>Fletcher</td>\n",
       "      <td>Mouse</td>\n",
       "      <td>Olfactory neurons</td>\n",
       "      <td>812</td>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "      <td>14.29%</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>0.0%</td>\n",
       "      <td>5.714286e-01</td>\n",
       "      <td>5.714286e-01</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      dataset organism          cell_type  cells_in_cluster  total_exons  \\\n",
       "0        Chen    Mouse                 ES               217          446   \n",
       "1        Chen    Mouse         Epi, early                98          402   \n",
       "2        Chen    Mouse          Epi, late               104          516   \n",
       "3        Chen    Mouse      Neuron, early                47          364   \n",
       "4        Chen    Mouse       Neuron, late                22          517   \n",
       "5   Lescroart    Mouse        Heart E6.75               172          286   \n",
       "6   Lescroart    Mouse        Heart E7.25               341          291   \n",
       "7    Trapnell    Human       Myoblast 00h                35          400   \n",
       "8    Trapnell    Human       Myoblast 24h                89          251   \n",
       "9    Trapnell    Human       Myoblast 48h                72          242   \n",
       "10   Trapnell    Human       Myoblast 72h                35          252   \n",
       "11       Song    Human               iPSC                62          616   \n",
       "12       Song    Human                NPC                73          212   \n",
       "13       Song    Human       Motor neuron                67          168   \n",
       "14     Shalek    Mouse               BMDC                13          149   \n",
       "15   Fletcher    Mouse  Olfactory neurons               812            7   \n",
       "\n",
       "    bimodal_exons bimodal_percent  selected_exons  selected_bimodal  \\\n",
       "0             118          26.46%              94                 0   \n",
       "1             107          26.62%              98                 1   \n",
       "2             136          26.36%              76                 1   \n",
       "3             117          32.14%              43                 0   \n",
       "4             146          28.24%              61                 0   \n",
       "5              77          26.92%              33                 0   \n",
       "6              78           26.8%              36                 0   \n",
       "7             142           35.5%              41                 0   \n",
       "8              97          38.65%              27                 0   \n",
       "9              97          40.08%              32                 1   \n",
       "10            101          40.08%              24                 1   \n",
       "11            269          43.67%              55                 0   \n",
       "12             92           43.4%              28                 1   \n",
       "13             82          48.81%              19                 4   \n",
       "14             51          34.23%              27                 1   \n",
       "15              1          14.29%               3                 0   \n",
       "\n",
       "   bimodal_percent_selected         p-val   p-val (adj)  \n",
       "0                      0.0%  4.089695e-15  3.271756e-14  \n",
       "1                     1.02%  1.743342e-14  9.297824e-14  \n",
       "2                     1.32%  2.969315e-10  1.111313e-09  \n",
       "3                      0.0%  1.559634e-08  3.564878e-08  \n",
       "4                      0.0%  3.472852e-10  1.111313e-09  \n",
       "5                      0.0%  1.512964e-05  2.017286e-05  \n",
       "6                      0.0%  5.495198e-06  7.993016e-06  \n",
       "7                      0.0%  4.519740e-09  1.205264e-08  \n",
       "8                      0.0%  7.074301e-07  1.162529e-06  \n",
       "9                     3.12%  4.541607e-07  9.083215e-07  \n",
       "10                    4.17%  4.112914e-05  5.062048e-05  \n",
       "11                     0.0%  2.547605e-15  3.271756e-14  \n",
       "12                    3.57%  7.265807e-07  1.162529e-06  \n",
       "13                   21.05%  8.825302e-03  9.413655e-03  \n",
       "14                     3.7%  5.783328e-05  6.609518e-05  \n",
       "15                     0.0%  5.714286e-01  5.714286e-01  "
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bimodality_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "bimodality_table.to_csv('selected_int_exons_agg_clusters.csv', index=False, header=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}