FALCOR_alg.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import colors\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "%matplotlib inline\n",
    "plt.rcParams['savefig.dpi'] = 300\n",
    "\n",
    "#stroke-miterlimit change this for svg file output\n",
    "%reload_ext autoreload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from fluctuation_code.MSS_MLE import *\n",
    "#imports the code set for completing the MLE optimization as a function of m, M\n",
    "#yay! completed writing this on\n",
    "\n",
    "#function set described below"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probability function of observing 'r' mutants in a given culture Given a rate m  $$p_0 = e^{-m}$$   $$p_r = \\frac{m}{r} \\sum_{i=0}^{r-1} \\frac{p_i}{(r-i+1)} $$   Maximum likelihood function (this is what needs to be maximized.   $$f(r|m) = \\prod_{i=1}^{c} f(r_i|m) = (p_0)^{c_0} (p_1)^{c_1} (p_2)^{c_2} (p_{r_{max}})^{c_{r_{max}}} $$    We also have the following MMS-MLE confidence intervals for the MSS-MLE method    $$CL_{+95\\%} = ln(m) + 1.96 \\sigma (e^{1.96 \\sigma})^{-0.315}$$  $$CL_{-95\\%} = ln(m) - 1.96 \\sigma (e^{1.96 \\sigma})^{+0.315}$$  with $$\\sigma_{ln(m)} \\approx \\frac{1.225m^{-0.315}}{\\sqrt{C}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code here is for developing the algorithm which executes automatic  MMS-MLE according to http://www.mitochondria.org/protocols/FALCOR.html "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "data as described in summary 2016_11_01 (good experimental quality, etc.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pr_ARRAY(m, r=0, P_R=[]) #takes trial value m mutations, with r observed number and computes probability\n",
    "Pr(m,r,P_R)              #non-array version of Pr_ARRAY\n",
    "\n",
    "mut_to_array(x,r_max)    #takes x colony values and converts to a counting array/histogram (mostly sparse usually)\n",
    "\n",
    "f_MLF(P_R, C)            #computes log of MLF function as np.log(P_R).dot(counts) \n",
    "MLF(m,Y,r_max=50)        #m is trial mutation rate (can be array), Y is mutation matrix, computes f_MLF\n",
    "                            #feature selects array or float depending on input\n",
    "                            #returns f_MLF value(s) for array of single m\n",
    "                            \n",
    "optimize_test            #returns x_1 and x_2 intervals that are 'best' for a given guess, used for finding optimum\n",
    "optimizer(Y,x_guess1, x_guess2, r_max, decay, draws, mult)\n",
    "                         #function takes many parameters and outputs limits, logm as values for the optimal 'm'\n",
    "                         \n",
    "confidence(m,cells,c=8)  #computes confidence intervals\n",
    "\n",
    "getM(x,cells,c=8)        #ties together the package. runs optimizer on a given set of X observations\n",
    "                         #returns confidence intervals and optimal 'm'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "def FALCOR_format(counts = [1,1,1,1,2,1,1,1], cellN =35000):\n",
    "    \n",
    "    conf, m = getM(counts, cellN, 16, r_max=170, decay=.5)\n",
    "\n",
    "    print('%.3f' % m)\n",
    "    print('%.3f' % conf[0])\n",
    "    print('%.3f' % conf[1])\n",
    "    \n",
    "    return\n",
    "\n",
    "#1686, 3912,266\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "110.859\n",
      "176.790\n",
      "57.155\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/aaroncp1an0/aaroncp1an0@gmail.com/lab_smolke/Afluctuationassay_tables/fluctuation_code/MSS_MLE.py:68: RuntimeWarning: divide by zero encountered in log\n",
      "  return np.log(P_R).dot(C)\n"
     ]
    }
   ],
   "source": [
    "FALCOR_format([0,2,4,3,6,1,3,0], cellN = 5*1E3*23)\n",
    "\n",
    "#15 (8.6,23)\n",
    "#13 (7, 20)\n",
    "#9 (4.5,15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def T_sig(m1,m2,culture1=8, culture2=8):\n",
    "    #degrees of freedom is C1+C2-2\n",
    "    sigma1 = np.e**(1.225*m1**-.315/np.sqrt(culture1))\n",
    "    sigma2 = np.e**(1.225*m2**-.315/np.sqrt(culture2))\n",
    "    return (np.log(m1) - np.log(m2)) / np.sqrt(sigma1**2/culture1 + sigma2**2/culture2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pval:   0.785175402994\n",
      "tval:   -0.277867561746\n"
     ]
    }
   ],
   "source": [
    "tv = T_sig(6.5,7.3, culture1=16, culture2=16)\n",
    "pval = stats.t.sf(np.abs(tv), 14)*2\n",
    "\n",
    "print('pval:   '+str(pval))\n",
    "print('tval:   '+str(tv))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.613\n",
      "9.805\n",
      "3.915\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/aaroncp1an0/aaroncp1an0@gmail.com/lab_smolke/Afluctuationassay_tables/fluctuation_code/MSS_MLE.py:68: RuntimeWarning: divide by zero encountered in log\n",
      "  return np.log(P_R).dot(C)\n"
     ]
    }
   ],
   "source": [
    "FALCOR_format([16,12,5,3,1,4,3,2,21,10,8,11,13,8,2,3,3,21,4,4,1,0,14,2,1\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    ],cellN=3660000\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[['1320000' '1500000' '1560000' '1560000' '1200000' '1320000']\n",
      " ['3,6,19,4,1,5,6,4' '0,2,2,0,0,16,0,6' '3,0,2,0,0,2,3,0' '1,1,1,2,0,2,2,0'\n",
      "  '2,5,6,7,12,15,4,6' '1,0,0,0,0,0,0,0']]\n",
      "OUTPUT\n",
      "[1320000, 1500000, 1560000, 1560000, 1200000, 1320000]\n",
      "[array([ 3,  6, 19,  4,  1,  5,  6,  4]), array([ 0,  2,  2,  0,  0, 16,  0,  6]), array([3, 0, 2, 0, 0, 2, 3, 0]), array([1, 1, 1, 2, 0, 2, 2, 0]), array([ 2,  5,  6,  7, 12, 15,  4,  6]), array([1, 0, 0, 0, 0, 0, 0, 0])]\n"
     ]
    }
   ],
   "source": [
    "filename='FALCOR_input.csv'\n",
    "filetext = np.genfromtxt(filename, dtype='U', delimiter=':')  #autostrip=True, )\n",
    "\n",
    "print(filetext) #.astype('str'))\n",
    "\n",
    "x = [int(i) for i in filetext[0]] #convert the first line to integer array\n",
    "y = [np.fromstring(j,dtype=int,sep=',') for j in filetext[1]]\n",
    "\n",
    "print('OUTPUT')\n",
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19.285\n",
      "28.439\n",
      "11.524\n",
      "\n",
      "\n",
      "5.442\n",
      "9.181\n",
      "2.487\n",
      "\n",
      "\n",
      "4.616\n",
      "7.918\n",
      "2.031\n",
      "\n",
      "\n",
      "5.777\n",
      "9.622\n",
      "2.718\n",
      "\n",
      "\n",
      "26.219\n",
      "37.832\n",
      "16.268\n",
      "\n",
      "\n",
      "0.947\n",
      "2.114\n",
      "0.180\n",
      "\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/aaroncp1an0/aaroncp1an0@gmail.com/lab_smolke/Afluctuationassay_tables/fluctuation_code/MSS_MLE.py:68: RuntimeWarning: divide by zero encountered in log\n",
      "  return np.log(P_R).dot(C)\n"
     ]
    }
   ],
   "source": [
    "for i in range(len(x)):\n",
    "    FALCOR_format(y[i],x[i])\n",
    "    print('\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "37817.8\t123\n"
     ]
    }
   ],
   "source": [
    "print('%.1f' % conf[0] + '\\t123')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'optimizer' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-12-1771bf71aa76>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmut_to_array\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcounts\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mlimits\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moptimizer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdecay\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m.1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdraws\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m70\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mmult\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlimits\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mlimits\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'g'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlimits\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mlimits\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'g'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'optimizer' is not defined"
     ]
    }
   ],
   "source": [
    "y=mut_to_array(counts)\n",
    "limits, m1 = optimizer(y,decay=.1,draws=70,mult=3)\n",
    "\n",
    "plt.plot(range(len(limits[:,0])),10**limits[:,0],color='g')\n",
    "plt.plot(range(len(limits[:,0])),10**limits[:,1],color='g')\n",
    "\n",
    "limits, m2 = optimizer(y,decay=.6,draws=70,mult=3)\n",
    "\n",
    "plt.plot(range(len(limits[:,0])),10**limits[:,0],color='b')\n",
    "plt.plot(range(len(limits[:,0])),10**limits[:,1],color='b')\n",
    "\n",
    "limits, m3 = optimizer(y,decay=.9,draws=70,mult=3)\n",
    "\n",
    "plt.plot(range(len(limits[:,0])),10**limits[:,0],color='r')\n",
    "plt.plot(range(len(limits[:,0])),10**limits[:,1],color='r')\n",
    "\n",
    "print(m1,m2,m3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}