coverage_factors.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculation of coverage factors of bipolar cell types in inner and outer plexiform layer<br/> based on convex hulls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import loadmat\n",
    "from scipy.spatial import ConvexHull\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from shapely.geometry import Polygon\n",
    "import shapely.ops\n",
    "%matplotlib inline\n",
    "matplotlib.rcParams.update({'font.size': 14})\n",
    "sns.set_style(\"whitegrid\")\n",
    "matplotlib.rcParams['font.family'] = 'sans-serif'\n",
    "matplotlib.rcParams['font.sans-serif'] = ['Arial']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "skeleton_data=loadmat('data/skeletons_OPL_final.mat')\n",
    "skeleton_ids=skeleton_data['kn_allSkeletons_e2006_IDs'].flatten()\n",
    "skeletons=skeleton_data['kn_allSkeletons_e2006'].flatten()\n",
    "M=np.loadtxt('data/coordinate_rotation')\n",
    "cell_ids=np.loadtxt('data/BC_IDs_new')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "soma_data=loadmat('data/soma_positions.mat')\n",
    "soma_pos=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_allSomata']\n",
    "soma_internal_ids=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_cellIDs']\n",
    "soma_line_ids=soma_data['kn_e2006_ALLSKELETONS_FINAL2012_cellIDs_pure_forSomata']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "BC_excluded=np.array([691,709,827,836])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def coverage(selection,layer):\n",
    "    area=[]\n",
    "    total_nodes=np.empty((0,2))\n",
    "    polygons=[]\n",
    "    for cell in cell_ids[selection,0]:\n",
    "        if not cell in BC_excluded:\n",
    "            nodes_cell=np.empty((0,3))\n",
    "            for skel in np.where(skeleton_ids==cell)[0]:\n",
    "                nodes=skeletons[skel].item()[list(skeletons[skel].dtype.names).index('nodes')][:,:3]\n",
    "                edges=skeletons[skel].item()[list(skeletons[skel].dtype.names).index('edges')]\n",
    "                nodes=nodes[np.unique(edges)-1,:]\n",
    "                if nodes.shape[0]<2:\n",
    "                    continue\n",
    "                nodes=nodes[:,:3]/[16.5,16.5,25]\n",
    "                nodes=np.dot(M,nodes.T).T\n",
    "                soma=np.dot(M,soma_pos[soma_line_ids[0,np.where(soma_internal_ids==cell_ids[cell_ids[:,0]==390,1])[1][0]]-1,:3].T).T\n",
    "                if layer=='OPL':\n",
    "                    nodes_cell=np.concatenate((nodes_cell,nodes[nodes[:,0]<soma[0],0:3]),axis=0)\n",
    "                elif layer=='IPL':\n",
    "                    nodes_cell=np.concatenate((nodes_cell,nodes[nodes[:,0]>soma[0],0:3]),axis=0) \n",
    "                else:\n",
    "                    print(\"Layer has to be 'IPL' or 'OPL'\")\n",
    "                    return(np.nan,np.nan,np.nan)\n",
    "            nodes_cell=nodes_cell[:,1:]/[2000/33,40]\n",
    "            if nodes_cell.shape[0]>2:\n",
    "                total_nodes=np.concatenate((total_nodes,nodes_cell),axis=0)\n",
    "                hull=ConvexHull(nodes_cell)\n",
    "                polygons.append(Polygon(nodes_cell[hull.vertices,:]))    \n",
    "    for i in range(len(polygons)):\n",
    "        area.append(polygons[i].area)\n",
    "    area=np.array(area)\n",
    "    hull=ConvexHull(total_nodes)\n",
    "    overlaps=[]\n",
    "    for i in range(len(polygons)):\n",
    "        for j in range(i+1,len(polygons)):\n",
    "            overlaps.append([cell_ids[selection,0][i],cell_ids[selection,0][j],polygons[i].intersection(polygons[j]).area])\n",
    "    overlaps=np.array(overlaps)\n",
    "    return(np.mean(area),np.std(area)/np.sqrt(len(area)-1),np.sum(area)/shapely.ops.cascaded_union(polygons).area)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "BC_data=np.zeros((14,8))\n",
    "BC_data[:,0]=np.arange(58,72)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for i in range(len(BC_data)):\n",
    "    BC_data[i,1]=np.sum((cell_ids[:,4]==BC_data[i,0])&(np.in1d(cell_ids[:,0],BC_excluded,invert=True)))\n",
    "    BC_data[i,2:5]=coverage(cell_ids[:,4]==BC_data[i,0],'OPL')\n",
    "    BC_data[i,5:8]=coverage(cell_ids[:,4]==BC_data[i,0],'IPL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "BC_data=pd.DataFrame(BC_data,columns=['type','n','OPL_hull_mean','OPL_hull_SEM','OPL_coverage','IPL_hull_mean','IPL_hull_SEM','IPL_coverage'])"
   ]
  },
  {
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}