Zhao et al. Figure 9 analysis

ZhaoFig9

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def doEM(exonCounts, exonLens, numIt=5000):\n",
    "    # the first and last exons are compatible with \n",
    "    # both isoforms\n",
    "    compatAB = exonCounts[0] + exonCounts[2]\n",
    "    # the middle exon is compatible only with B\n",
    "    compatB = exonCounts[1] \n",
    "    la = exonLengths[0] + exonLengths[2]\n",
    "    lb = exonLengths[0] + exonLengths[1] + exonLengths[2]\n",
    "    \n",
    "    alphas = [0.5, 0.5]\n",
    "    for i in xrange(numIt):\n",
    "        s = (alphas[0] / la) + (alphas[1] / lb)\n",
    "        sa = compatAB * ((alphas[0] / la) / (s))\n",
    "        sb = compatAB * ((alphas[1] / lb) / (s))\n",
    "        sb += compatB\n",
    "        alphas = [sa, sb]\n",
    "        \n",
    "    return alphas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# They specify only that exons 1 and 2 are twice as long as exon 3\n",
    "exonLengths = [100, 100, 50]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# reads assigned to A = 180.0, # reads assigned to B = 8.89318162514e-322\n"
     ]
    }
   ],
   "source": [
    "# Subpanel B (top)\n",
    "exonCountsBT = [120, 0, 60]\n",
    "resBT = doEM(exonCountsBT, exonLengths)\n",
    "print(\"# reads assigned to A = {}, # reads assigned to B = {}\".format(resBT[0], resBT[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# reads assigned to A = 18.0, # reads assigned to B = 7.11454530011e-322\n"
     ]
    }
   ],
   "source": [
    "# Subpanel B (bottom)\n",
    "exonCountsBB = [12, 0, 6]\n",
    "resBB = doEM(exonCountsBB, exonLengths)\n",
    "print(\"# reads assigned to A = {}, # reads assigned to B = {}\".format(resBB[0], resBB[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# reads assigned to A = 0.0899460323806, # reads assigned to B = 299.910053968\n"
     ]
    }
   ],
   "source": [
    "# Subpanel C (top)\n",
    "exonCountsCT = [120, 120, 60]\n",
    "resCT = doEM(exonCountsCT, exonLengths)\n",
    "print(\"# reads assigned to A = {}, # reads assigned to B = {}\".format(resCT[0], resCT[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# reads assigned to A = 0.00899460323806, # reads assigned to B = 29.9910053968\n"
     ]
    }
   ],
   "source": [
    "# Subpanel C (bottom)\n",
    "exonCountsCB = [12, 12, 6]\n",
    "resCB = doEM(exonCountsCB, exonLengths)\n",
    "print(\"# reads assigned to A = {}, # reads assigned to B = {}\".format(resCB[0], resCB[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}