05_datafilterABA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nre-do of ABA ISH filtering to match ST protocol and test if it makes a difference from before\\n\\nShaina Lu\\nZador & Gillis Labs\\nApril 2020\\n'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "re-do of ABA ISH filtering to match ST protocol and test if it makes a difference from before\n",
    "\n",
    "Shaina Lu\n",
    "Zador & Gillis Labs\n",
    "April 2020\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import h5py\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in raw ABA ISH Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "ALLEN_PATH = \"/data/slu/allen_adult_mouse_ISH/allen_adultmouse_manthan_v3.h5\"\n",
    "ONTOLOGY_PATH = \"/data/slu/allen_adult_mouse_ISH/ontologyABA.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['gene.id.info',\n",
       " 'gene.id.info.cols',\n",
       " 'vox.df',\n",
       " 'vox.df.cols',\n",
       " 'vox.meta',\n",
       " 'vox.meta.cols']"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#open gene by voxel hdf5 file\n",
    "f = h5py.File(ALLEN_PATH, 'r')\n",
    "list(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#unpack each matrix within hdf5 file\n",
    "metacols = list(map(lambda x: x.decode('UTF-8'), list(f['vox.meta.cols'])))\n",
    "dfcols = list(map(lambda x: x.decode('UTF-8'),list(f['vox.df.cols'])))\n",
    "\n",
    "meta = pd.DataFrame(np.char.decode(np.array(f['vox.meta']).T), columns=metacols)\n",
    "voxdf = pd.DataFrame(np.array(f['vox.df']).T, columns=dfcols)\n",
    "geneIDName = pd.DataFrame(np.char.decode(np.array(f['gene.id.info']).T), columns=[\"expt_id\",\"gene_symbol\",\"entrez_id\",\"specimen_id\",\"plane\"])\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(26079, 5)\n",
      "(159326, 26078)\n",
      "(159326, 4)\n"
     ]
    }
   ],
   "source": [
    "#note one more gene in geneIDName than in the cols of voxdf, 1 non-uniq in saggittal\n",
    "print(geneIDName.shape)\n",
    "print(voxdf.shape)\n",
    "print(meta.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#read in ontology file\n",
    "ontology = pd.read_csv(ONTOLOGY_PATH)\n",
    "ontology = ontology.drop([ontology.columns[5], ontology.columns[6]], axis=1)\n",
    "ontology = ontology.fillna(-1)  #make root's parent -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__what's happening with the duplicate in geneIDName?__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "expt_id        26078\n",
       "gene_symbol    19942\n",
       "entrez_id      19428\n",
       "specimen_id     5440\n",
       "plane              2\n",
       "dtype: int64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geneIDName.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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>expt_id</th>\n",
       "      <th>gene_symbol</th>\n",
       "      <th>entrez_id</th>\n",
       "      <th>specimen_id</th>\n",
       "      <th>plane</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2744</th>\n",
       "      <td>78153149</td>\n",
       "      <td>Avp</td>\n",
       "      <td>11998</td>\n",
       "      <td>77632911</td>\n",
       "      <td>sagittal</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2745</th>\n",
       "      <td>78153149</td>\n",
       "      <td>Avp</td>\n",
       "      <td>11998</td>\n",
       "      <td>77632911</td>\n",
       "      <td>sagittal</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       expt_id gene_symbol entrez_id specimen_id     plane\n",
       "2744  78153149         Avp     11998    77632911  sagittal\n",
       "2745  78153149         Avp     11998    77632911  sagittal"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#there are two rows from the sagittal plane that are duplicated across all entries\n",
    "dupfilt = geneIDName.duplicated(keep=False)\n",
    "geneIDName.loc[dupfilt, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#remove one of these rows\n",
    "geneIDName = geneIDName.drop(labels=2745, axis='index')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(26078, 5)\n",
      "expt_id        26078\n",
      "gene_symbol    19942\n",
      "entrez_id      19428\n",
      "specimen_id     5440\n",
      "plane              2\n",
      "dtype: int64\n"
     ]
    }
   ],
   "source": [
    "print(geneIDName.shape)\n",
    "print(geneIDName.nunique())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filter I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__remove spots that do not map to the brain first here since there are many this will speed up subsequent propogation of ontology__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#filter out voxels that don't map to the brain\n",
    "voxbrain = voxdf.loc[meta.ids != 'NA',:]\n",
    "metabrain = meta.loc[meta.ids != 'NA',:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(62529, 26078)\n",
      "(62529, 4)\n"
     ]
    }
   ],
   "source": [
    "print(voxbrain.shape)\n",
    "print(metabrain.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# propagate ontology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/slu/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:15: FutureWarning: `item` has been deprecated and will be removed in a future version\n",
      "  from ipykernel import kernelapp as app\n"
     ]
    }
   ],
   "source": [
    "#do this in a sloppy way that i know works\n",
    "#propogate first then look-up\n",
    "#goes up to the point where all labels are root\n",
    "proponto = pd.DataFrame()\n",
    "proponto[['id','1']] = ontology[['id', 'parent']]\n",
    "\n",
    "curr = 1\n",
    "while True:\n",
    "    #get series representing next level up\n",
    "    newcol = []\n",
    "    for i in range(proponto.shape[0]):\n",
    "        if proponto.iloc[i,curr] == -1:   #if already at parent of root, just give root again (-1)\n",
    "            newval = -1\n",
    "        else:\n",
    "            newval = ontology.loc[ontology[\"id\"] == int(proponto.iloc[i,curr]), \"parent\"].item()\n",
    "        newcol.append(newval)\n",
    "\n",
    "    if len(proponto[str(curr)].value_counts()) < 3:  #all are root after this point\n",
    "        break\n",
    "    proponto[str(curr+1)] = pd.Series(newcol)     #add series to df\n",
    "    curr += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.0      1281\n",
       " 997.0       6\n",
       "Name: 10, dtype: int64"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#at 10 levels up, all labels are root\n",
    "#proponto.head()\n",
    "proponto.loc[:,'10'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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>id</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>10</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>577</td>\n",
       "      <td>369</td>\n",
       "      <td>322</td>\n",
       "      <td>453</td>\n",
       "      <td>315</td>\n",
       "      <td>695</td>\n",
       "      <td>688</td>\n",
       "      <td>567</td>\n",
       "      <td>8</td>\n",
       "      <td>997</td>\n",
       "      <td>997</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>657</td>\n",
       "      <td>345</td>\n",
       "      <td>322</td>\n",
       "      <td>453</td>\n",
       "      <td>315</td>\n",
       "      <td>695</td>\n",
       "      <td>688</td>\n",
       "      <td>567</td>\n",
       "      <td>8</td>\n",
       "      <td>997</td>\n",
       "      <td>997</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1114</td>\n",
       "      <td>402</td>\n",
       "      <td>669</td>\n",
       "      <td>315</td>\n",
       "      <td>695</td>\n",
       "      <td>688</td>\n",
       "      <td>567</td>\n",
       "      <td>8</td>\n",
       "      <td>997</td>\n",
       "      <td>997</td>\n",
       "      <td>997</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>606</td>\n",
       "      <td>886</td>\n",
       "      <td>254</td>\n",
       "      <td>315</td>\n",
       "      <td>695</td>\n",
       "      <td>688</td>\n",
       "      <td>567</td>\n",
       "      <td>8</td>\n",
       "      <td>997</td>\n",
       "      <td>997</td>\n",
       "      <td>997</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>472</td>\n",
       "      <td>426</td>\n",
       "      <td>403</td>\n",
       "      <td>278</td>\n",
       "      <td>477</td>\n",
       "      <td>623</td>\n",
       "      <td>567</td>\n",
       "      <td>8</td>\n",
       "      <td>997</td>\n",
       "      <td>997</td>\n",
       "      <td>997</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     id    1    2    3    4    5    6    7    8    9   10\n",
       "0   577  369  322  453  315  695  688  567    8  997  997\n",
       "1   657  345  322  453  315  695  688  567    8  997  997\n",
       "2  1114  402  669  315  695  688  567    8  997  997  997\n",
       "3   606  886  254  315  695  688  567    8  997  997  997\n",
       "4   472  426  403  278  477  623  567    8  997  997  997"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "proponto = proponto.astype(int)  #typecast as integer\n",
    "proponto = proponto.replace(to_replace=-1, value=997)  #replace root's parent (-1) with root (997)\n",
    "proponto.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__now propagate full binary ontology for voxels__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "fullontST = pd.DataFrame(index=metabrain.index, columns=ontology.id, dtype=np.int64)\n",
    "fullontST.columns = fullontST.columns.astype(str)\n",
    "fullontST = fullontST.fillna(0) #make all entries 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(metabrain.ids.unique())): #for each unique brain area\n",
    "    tempareas = proponto.loc[proponto.id == int(metabrain.ids.unique()[i]), :] #get row representing that area propagated\n",
    "    #tempareas = tempareas.values.flatten()\n",
    "    #print(tempareas)\n",
    "    for val in tempareas.values.flatten():\n",
    "        #set row,col entries to 1 \n",
    "        fullontST.loc[metabrain.ids == metabrain.ids.unique()[i], str(val)] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "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>id</th>\n",
       "      <th>577</th>\n",
       "      <th>657</th>\n",
       "      <th>1114</th>\n",
       "      <th>606</th>\n",
       "      <th>472</th>\n",
       "      <th>1117</th>\n",
       "      <th>10714</th>\n",
       "      <th>632</th>\n",
       "      <th>484682479</th>\n",
       "      <th>312782586</th>\n",
       "      <th>...</th>\n",
       "      <th>356</th>\n",
       "      <th>153</th>\n",
       "      <th>684</th>\n",
       "      <th>312782554</th>\n",
       "      <th>484682520</th>\n",
       "      <th>1076</th>\n",
       "      <th>234</th>\n",
       "      <th>476</th>\n",
       "      <th>1051</th>\n",
       "      <th>619</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>9284</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9285</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9286</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9287</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9349</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 1287 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "id    577  657  1114  606  472  1117  10714  632  484682479  312782586  ...  \\\n",
       "9284  0.0  0.0   0.0  0.0  0.0   0.0    0.0  0.0        0.0        0.0  ...   \n",
       "9285  0.0  0.0   0.0  0.0  0.0   0.0    0.0  0.0        0.0        0.0  ...   \n",
       "9286  0.0  0.0   0.0  0.0  0.0   0.0    0.0  0.0        0.0        0.0  ...   \n",
       "9287  0.0  0.0   0.0  0.0  0.0   0.0    0.0  0.0        0.0        0.0  ...   \n",
       "9349  0.0  0.0   0.0  0.0  0.0   0.0    0.0  0.0        0.0        0.0  ...   \n",
       "\n",
       "id    356  153  684  312782554  484682520  1076  234  476  1051  619  \n",
       "9284  0.0  0.0  0.0        0.0        0.0   0.0  0.0  0.0   0.0  0.0  \n",
       "9285  0.0  0.0  0.0        0.0        0.0   0.0  0.0  0.0   0.0  0.0  \n",
       "9286  0.0  0.0  0.0        0.0        0.0   0.0  0.0  0.0   0.0  0.0  \n",
       "9287  0.0  0.0  0.0        0.0        0.0   0.0  0.0  0.0   0.0  0.0  \n",
       "9349  0.0  0.0  0.0        0.0        0.0   0.0  0.0  0.0   0.0  0.0  \n",
       "\n",
       "[5 rows x 1287 columns]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fullontST.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(62529, 1287)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fullontST.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__how many of the voxels map to ventricular systems or fiber tracts?__\n",
    "these areas were removed in ST data so remove them here too"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "#ventricular systesm - 2 voxels map here\n",
    "print(fullontST.loc[:,\"73\"].sum())\n",
    "#fiber tracts - 0 map here\n",
    "print(fullontST.loc[:,\"1009\"].sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "79019     3.0\n",
       "120053    4.0\n",
       "dtype: float64"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#which are to two voxels that map to ventricular systems and where else do they map?\n",
    "vsvox = fullontST.loc[fullontST[\"73\"]==1,\"73\"].index\n",
    "fullontST.loc[vsvox,:].sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#find the non-zero cols\n",
    "index1 = fullontST.loc[vsvox[0],:].to_numpy().nonzero()\n",
    "index2 = fullontST.loc[vsvox[1],:].to_numpy().nonzero()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Index(['73', '997', '129'], dtype='object', name='id')\n",
      "Index(['145', '73', '997', '153'], dtype='object', name='id')\n"
     ]
    }
   ],
   "source": [
    "#these are all ventricular areas\n",
    "print(fullontST.columns[index1])\n",
    "print(fullontST.columns[index2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/slu/miniconda3/lib/python3.7/site-packages/pandas/core/frame.py:4097: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  errors=errors,\n"
     ]
    }
   ],
   "source": [
    "#remove the two ventricular voxels\n",
    "fullontST.drop(labels=vsvox, axis=0, inplace=True)\n",
    "voxbrain.drop(labels=vsvox, axis=0, inplace=True)\n",
    "metabrain.drop(labels=vsvox, axis=0, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(62527, 1287)\n",
      "(62527, 26078)\n",
      "(62527, 4)\n"
     ]
    }
   ],
   "source": [
    "print(fullontST.shape)\n",
    "print(voxbrain.shape)\n",
    "print(metabrain.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filter II"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__filter genes, replace NA's (here -1) with 0 and remove genes that don't express in X% voxels__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/slu/miniconda3/lib/python3.7/site-packages/pandas/core/frame.py:4258: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  method=method,\n"
     ]
    }
   ],
   "source": [
    "#treats 0's as not expressing as well\n",
    "voxbrain.replace(to_replace=-1, value=0, inplace=True)            #make missing values 0\n",
    "nonzerocount = voxbrain.astype(bool).sum(axis=0)    #gets number of spots for each gene that have non-zero expression\n",
    "\n",
    "voxbrain = voxbrain.loc[:, nonzerocount > (.001*voxbrain.shape[0])]  #if gene has expression in at least 0.1% of samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(62527, 26008)\n"
     ]
    }
   ],
   "source": [
    "voxbrain.head()\n",
    "print(voxbrain.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: this is only two less genes than before in allenadultmouseISH/filter_one so shouldn't dramatically change anything"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write to new hdf5 part I\n",
    "__write filtered voxels, meta data, and propagated ontology to new hdf5 file__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "OUT_PATH = \"ABAISH_filt_v6.h5\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "voxbrain.to_hdf(OUT_PATH, key=\"voxbrain\", mode='w')\n",
    "metabrain.to_hdf(OUT_PATH, key=\"metabrain\")\n",
    "fullontST.to_hdf(OUT_PATH, key=\"propontology\")\n",
    "geneIDName.to_hdf(OUT_PATH, key=\"geneIDName\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Average duplicated genes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "#I want to average all genes that are repeated regardless of plane. \n",
    "#initialize dataframe with only non-repeated series (cols)\n",
    "nodupfilt = geneIDName.duplicated(subset=\"gene_symbol\", keep=False) #marks all duplicates as True\n",
    "nodupgenes = geneIDName.loc[~nodupfilt,:] #14970\n",
    "indexers = nodupgenes.loc[nodupgenes.expt_id.isin(list(voxbrain)),\"expt_id\"] #solves NAN cols\n",
    "avgvoxbrain = voxbrain.loc[:,indexers]\n",
    "#rename the columns using gene_symbols\n",
    "temp = dict(zip(nodupgenes.expt_id.tolist(),nodupgenes.gene_symbol.tolist()))\n",
    "avgvoxbrain.rename(columns=temp, inplace=True)\n",
    "\n",
    "#get the remaining gene_symbols that are duplicated\n",
    "dupgenes = geneIDName.loc[nodupfilt,:]\n",
    "#get unique gene_symbols from duplicated\n",
    "uniqfilt = dupgenes.duplicated(subset=\"gene_symbol\", keep='first')\n",
    "uniqgenes = dupgenes.loc[~uniqfilt,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4972"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#how many genes are duplicated?\n",
    "dupgenes.gene_symbol.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(uniqgenes.shape[0]):  #loop through unique duplicated genes\n",
    "    #get all series ids for current uniq gene symbol\n",
    "    currseries = dupgenes.loc[dupgenes.gene_symbol == uniqgenes.iloc[i,1], \"expt_id\"].values\n",
    "    curravg = voxbrain.loc[:,currseries].mean(axis=1) #get voxels and average them\n",
    "    curravg = curravg.rename(uniqgenes.iloc[i,1])     #rename and add to df\n",
    "    avgvoxbrain = pd.concat([avgvoxbrain, curravg], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(62527, 19934)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "avgvoxbrain.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write to new hdf5 part II"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__write averaged vox brain to hdf5__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "OUT_PATH2 = \"ABAISH_filt_v6_avgdup.h5\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "#write to second out file\n",
    "avgvoxbrain.to_hdf(OUT_PATH2, key=\"avgvoxbrain\", mode='w')\n",
    "metabrain.to_hdf(OUT_PATH2, key=\"metabrain\")\n",
    "fullontST.to_hdf(OUT_PATH2, key=\"propontology\")\n",
    "geneIDName.to_hdf(OUT_PATH2, key=\"geneIDName\")"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}