1 - Connectomics Getting Started.ipynb

{
 "metadata": {
  "name": "1 - Connectomics Getting Started"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "\n# Connectomics Challenge\n\nThis challenge requires to discover directed neural connections from time series of neural imagery depicting neuron activation patterns.\n\nMore at http://www.kaggle.com/c/connectomics\n\n## Links and Resources\n\n### Connectomics Challenge Resources\n\n * [Description of the connectomics challenge data set](http://www.kaggle.com/c/connectomics/data)\n * [Connectomics Knowledgebase](http://connectomics.chalearn.org/help/data)\n * [Simple tool helping to select an approproate Network Reconstruction Algorithm](http://www.network-reconstruction.org/)\n\n### Papers\n\n * [Transfer Entropy reconstruction and labeling of neuronal connections from simulated calcium imaging](http://arxiv.org/abs/1309.4287)\n * [Model-Free Reconstruction of Excitatory Neuronal Connectivity from Calcium Imaging Signals](http://www.ploscompbiol.org/article/info%3Adoi%2F10.1371%2Fjournal.pcbi.1002653)\n * [Granger causality and transfer entropy are equivalent for\nGaussian variables](http://www.sussex.ac.uk/Users/lionelb/downloads/NCOMP/publications/tegc_PRL_pre.pdf)\n * [Introduction to Neuroscience](http://neuroscience.uth.tmc.edu/s1/introduction.html)\n * [Imaging Calcium in Neurons](http://www.cell.com/neuron/abstract/S0896-6273%2812%2900172-9)\n \n\nWe start by a listing of the installed python packages using pip freeze, so all of this can be reproduced. We're also setting matplotlib to inline output, and load the cythonmagic extension for IPython (see http://ipython.org/ipython-doc/dev/config/extensions/cythonmagic.html )\n"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "!pip freeze",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "Cython==0.20.1\r\nJinja2==2.7.2\r\nMarkupSafe==0.18\r\nPygments==1.6\r\nSphinx==1.2.1\r\nTheano==0.6.0\r\nargparse==1.2.1\r\nbackports.ssl-match-hostname==3.4.0.2\r\ndistribute==0.7.3\r\ndocutils==0.11\r\nh5py==2.2.1\r\nipython==1.1.0\r\nmatplotlib==1.3.1\r\nnetworkx==1.8.1\r\nnose==1.3.0\r\nnumexpr==2.3.1\r\nnumpy==1.8.0\r\npandas==0.13.0\r\npatsy==0.2.1\r\npebl==1.01\r\npydot==1.0.28\r\npymc==2.3\r\npyparsing==2.0.1\r\npython-dateutil==2.2\r\npytz==2013.9\r\npyzmq==14.0.1\r\nscikit-learn==0.14.1\r\nscipy==0.13.3\r\nsix==1.5.2\r\nstatsmodels==0.5.0\r\ntables==3.1.0\r\ntornado==3.2\r\nwsgiref==0.1.2\r\n"
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%matplotlib inline\n%load_ext cythonmagic",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "\nNext, download and extract the data packages linked at http://www.kaggle.com/c/connectomics/data to a\nsuitable data directory"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# Some core imports\nimport os\nimport sys\nfrom subprocess import call\nimport time\nimport pandas\nimport numpy as np\nimport h5py\n# These are IPython specific\nfrom IPython.display import display, clear_output\n",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# Core configuration\n\ndatadir = \"/data/quick/connectomics\"\nmirror = \"http://www.causality.inf.ethz.ch/data/connectomics\"\n\nnbdir = os.getcwd()\n\nsources = [\n           \n           \"http://www.kaggle.com/c/connectomics/download/normal-1.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/normal-2.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/normal-3.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/normal-4.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/normal-4-lownoise.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/normal-4-lownoise.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/small.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/highcc.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/lowcc.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/lowcon.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/highcon.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/normal-3-highrate.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/test.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/valid.tgz\",\n           \"http://www.kaggle.com/c/connectomics/download/sampleSubmission.csv.zip\"\n          ]\n\n",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "print \"Creating data directory\"\ncall(['mkdir', '-p', datadir + \"/downloads\"])\ncall(['mkdir', '-p', datadir + \"/tmp\"])\nos.chdir(datadir + \"/tmp\")",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "Creating data directory\n"
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "\ncmd = ['wget', '-O', 'target', 'src']\n\nfor src in sources:\n    clear_output()\n    print \"Downloading %s\" % (src)\n    sys.stdout.flush()\n    parts = src.rsplit('/', 1)\n    filename = parts[-1]\n    target = \"%s/tmp/%s\" % (datadir, filename)\n    mirror_src = \"%s/%s\" % (mirror, filename)\n    try:\n        os.unlink(target)\n    except:\n        pass\n    cmd[2] = target\n    cmd[3] = mirror_src\n    #print ' '.join(cmd)\n    call(cmd)\n    tgz_cmd = [\"tar\", \"-xf\", filename]\n    print \"Extracting %s\" % (filename)\n    sys.stdout.flush()\n    call(tgz_cmd)\n    print \"Removing %s\" % (filename)\n    sys.stdout.flush()\n    os.unlink(target)\nclear_output()\nprint \"Downloads finished\"",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "Downloads finished\n"
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "## Making the data accessible\n\nWe are going to make the data a bit more accessible by loading it into a hierarchical HDF5 data store using the h5py library ( see http://docs.h5py.org/en/latest/ ). While pandas also supports HDF5, it does not natively support multidimensional tensors/matrices which we will use later on in conjunction with Theano.\n"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "def import_network(store, name, f_file, p_file, n_file):\n    print \"Importing %s\" % (name)\n    fdata = pandas.read_csv(f_file, sep=',', skipinitialspace=True, header=None, prefix='N', dtype=np.float32)\n    concatenated = pandas.concat([fdata[c] for c in fdata.columns])\n    timecount = len(fdata)\n    neuroncount = len(fdata.columns)\n    fdata = None\n    tdata = concatenated.values.reshape((neuroncount, timecount))\n    g = store.create_group(name)\n    dset = g.create_dataset(\"fluorescence\", data=tdata, compression='lzf')\n    \n    g['fluorescence'].dims[0].label = 'neuron'\n    g['fluorescence'].dims[1].label = 'time'\n    pdata = pandas.read_csv(p_file,  sep=',', skipinitialspace=True, header=None, names=['x','y'], dtype=np.float64)\n    concatenated = pandas.concat([pdata[c] for c in pdata.columns])\n    tensor = concatenated.values.reshape((2, len(pdata))).T\n    pset = g.create_dataset(\"networkPositions\", data=tensor)\n    if ((n_file is not None) and os.path.isfile(n_file)):\n        ndata = pandas.read_csv(n_file,  sep=',', skipinitialspace=True, header=None, names=['i','j', 'weight'], dtype={ 'i' : np.int32, 'j' : np.int32, 'weight' : np.float32})\n        connection_matrix = np.zeros((neuroncount, neuroncount), dtype=np.float32)\n        for idx in range(len(ndata)):\n            i = int(ndata.i[idx])\n            j = int(ndata.j[idx])\n            weight = float(ndata.weight[idx])\n            if (weight==-1.0):\n                weight = 0.0\n            connection_matrix[i-1,j-1] = weight\n        g.create_dataset(\"connectionMatrix\", data=connection_matrix, compression='lzf')\n    return g",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "store = h5py.File(datadir + '/store.h5')\n",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "\nfor i in range(1,7):\n    import_network(store, 'small-%d' % (i),\n                   'fluorescence_iNet1_Size100_CC0%dinh.txt' % (i), \n                   'networkPositions_iNet1_Size100_CC0%dinh.txt'  % (i),\n                   'network_iNet1_Size100_CC0%dinh.txt'  % (i))\n    \nfor nname in ['normal-1', 'normal-2', 'normal-3', 'normal-4', 'normal-4-lownoise', 'lowcc', 'highcc', 'highcon', 'test', 'valid']:\n        import_network(store, nname, \n                       \"%s/fluorescence_%s.txt\" % (nname, nname),\n                       \"%s/networkPositions_%s.txt\" % (nname, nname),\n                       \"%s/network_%s.txt\" % (nname, nname))\nprint \"Done\"\nprint store",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "Importing small-1\nImporting small-2"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting small-3"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting small-4"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting small-5"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting small-6"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting normal-1"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting normal-2"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting normal-3"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting normal-4"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting normal-4-lownoise"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting lowcc"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting highcc"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting highcon"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting test"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nImporting valid"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\nDone"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\n<HDF5 file \"store.h5\" (mode r+)>\n"
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "That's it, we have the raw data in a nice accessible format in a HDF5 store and can access the individual datasets like this:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "store['normal-1']['connectionMatrix']\n",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 29,
       "text": "<HDF5 dataset \"connectionMatrix\": shape (1000, 1000), type \"<f4\">"
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "## Network Connectivity Postprocessing\n\nNext, we are going to calculate a neuron distance matrix and \"causal distance\" matrices which determine how close a neuron is from another one, both in spatial distance and how many \"hops\" would be neccessary for a spike to traverse from one neuron to another."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "\nfrom math import sqrt\n\ndef create_spatial_distance_matrix(net):\n    \n    positions = net['networkPositions'].value\n    \n    def distance(n1, n2):\n        x1 = positions[n1,0]\n        y1 = positions[n1,1]\n        x2 = positions[n2,0]\n        y2 = positions[n2,1]\n        return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))\n    nneurons = positions.shape[0]\n    distanceMatrix = np.zeros((nneurons, nneurons), dtype=np.float32)\n    for n1 in range(nneurons):\n        for n2 in range(positions.shape[0]):\n            distanceMatrix[n1,n2] = distance(n1,n2)\n    try:\n        del net['spatialDistanceMatrix']\n    except:\n        pass\n    net.create_dataset(\"spatialDistanceMatrix\", data=distanceMatrix)\n    \n",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 63
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "\n\ndef update_signal_distance_matrix(net):\n    \n    if (not 'connectionMatrix' in net.keys()):\n        return False\n    connections = net['connectionMatrix'].value\n    try:\n        del net['signalDistanceMatrix']\n    except:\n        pass\n    inf = float('inf')\n    \n    signalMatrix = np.zeros_like(connections, dtype=np.float32)\n    signalMatrix[:,:] = connections[:,:]\n    signalMatrix[signalMatrix==0.0] = inf\n    for i in range(signalMatrix.shape[0]):\n        signalMatrix[i,i] = 0.0\n    \n    def signal(src, via):\n        depth = signalMatrix[src,via]\n        signal = signalMatrix[via,:] + float(depth)\n        signalMatrix[src,:] = np.minimum(signalMatrix[src,:], signal)\n    \n    \n    for i in range(1000):\n        orig = signalMatrix.copy()\n        for src in range(signalMatrix.shape[0]):\n            for via in range(signalMatrix.shape[0]):\n                if (via==target):\n                    continue\n                if (signalMatrix[src,via]>=inf):\n                    continue\n                signal(src, via)\n        if (np.all(signalMatrix==orig)):\n            break\n    \n    net.create_dataset(\"signalDistanceMatrix\", data=signalMatrix)\n    print \"Done updating\"\n",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 105
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "from sys import stdout\nfor nwname in store.keys():\n    net = store[nwname]\n    print \"Network: %s\" % (nwname)\n    if (('networkPositions' in net.keys())):\n        print \"Calculating spatial distance matrix\"\n        stdout.flush()\n        create_spatial_distance_matrix(net)\n        print \"Calculating signal distance matrix\"\n        stdout.flush()\n        update_signal_distance_matrix(net)",
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Let's verify that the code is correct, by counter-checking with the shortest paths given by the NetworkX Library:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import networkx as nx\n\ndef create_nx_graph(net):\n    ''' Load a neural network with known connectivity as a NetworkX directed graph'''\n    if (not 'connectionMatrix' in net):\n        return None\n    G = nx.DiGraph() \n    conn = net['connectionMatrix']\n    for i in range(conn.shape[0]):\n        G.add_node(i)\n    for i in range(conn.shape[0]):\n        for j in range(conn.shape[0]):\n            if (conn[i,j]>0.0):\n                G.add_edge(i,j)\n    return G\n\n",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 89
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "\n# Simple test case: Check all connections of the small-1 network\nG = create_nx_graph(store['small-1'])\nfor i in range(100):\n    for j in range(100):\n        try:\n            plen = len(nx.shortest_path(G, i, j))-1 \n            assert plen == store['small-1']['signalDistanceMatrix'][i,j]\n        except:\n            assert store['small-1']['signalDistanceMatrix'][i,j]==float('inf')\nprint \"Test passed\"",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "Test passed\n"
      }
     ],
     "prompt_number": 104
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "store.close()",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    }
   ],
   "metadata": {}
  }
 ]
}