{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/anndata/_core/anndata.py:21: FutureWarning: pandas.core.index is deprecated and will be removed in a future version.  The public classes are available in the top-level namespace.\n",
      "  from pandas.core.index import RangeIndex\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as pl\n",
    "from matplotlib import rcParams\n",
    "import scanpy as sc\n",
    "import seaborn as sns\n",
    "from scipy.sparse import csr_matrix\n",
    "from scipy.stats import rankdata\n",
    "import pickle\n",
    "import statistics as stat\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load each lineage - either low or high clustering resolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "endo = sc.read_h5ad(\"../data/endo.dpt.h5ad\")\n",
    "hep = sc.read_h5ad(\"../data/hep.dpt.h5ad\")\n",
    "neur = sc.read_h5ad(\"../data/neur.dpt.h5ad\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "endo = sc.read_h5ad(\"../data/endo.dpt.hires.h5ad\")\n",
    "hep = sc.read_h5ad(\"../data/hep.dpt.hires.h5ad\")\n",
    "neur = sc.read_h5ad(\"../data/neur.dpt.hires.h5ad\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Useful functions for pseudobulk aggregation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def aggregate(sc_sub):\n",
    "    \"\"\"\n",
    "    input: sc_sub, a scanpy h5ad object, [cells x genes]\n",
    "    output: mat, a csr matrix, [1 x genes], summed over cells\n",
    "    \"\"\"\n",
    "    mat = csr_matrix.sum(sc_sub.X, 0)\n",
    "    return(mat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_normalize(mat):\n",
    "    \"\"\"\n",
    "    input: mat, a matrix, [1 x genes], containing raw counts\n",
    "    output: logcp10k, a csr matrix, [1 x genes], = log10(counts per 10k UMI + 1)\n",
    "    \"\"\"\n",
    "    cp10k = mat * 1e4 / np.sum(mat, 1)\n",
    "    logcp10k = np.log10(cp10k + 1)\n",
    "    return(logcp10k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scale(mat, center=True):\n",
    "    \"\"\"\n",
    "    input: mat, a [levels x gene] csr matrix\n",
    "           center, boolean indicating whether or not you want to center the data\n",
    "    output: scaled_mat, a [levels x gene] csr matrix where each gene has zero mean, unit variance across all levels\n",
    "    \"\"\"\n",
    "    if center:\n",
    "        mat = mat - np.mean(mat, 0)\n",
    "    stdevs = np.std(mat, 0)\n",
    "    stdevs[stdevs == 0] = 1\n",
    "    scaled_mat = mat / stdevs\n",
    "    return(scaled_mat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Master function for pseudobulk aggregation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_pseudobulk(sc_obj, stratifier=\"individual\", pt_label=\"dpt_pseudotime\", nbins=5):\n",
    "    \"\"\"\n",
    "    input: sc_obj, a full h5ad object\n",
    "           stratifier, the string identifier for the category over which pseudobulk will be aggregated (along with pseudotime)\n",
    "           pt_label, the string identifier for pseudotime, ie 'dpt_pseudotime'\n",
    "           nbins, the number of pseudotime bins to compute\n",
    "    output: Y, an [NGT, 1] array where N=number of levels to stratifier, G=num genes, T=nbins\n",
    "            X, an [NGT, 3] matrix (N,G,T as above), where the first column is time, second column is level, third is gene\n",
    "            leveldict, a dictionary mapping ints in second column of X to unique values in sc_obj.obs[stratifier]\n",
    "            genedict, a dictionary mapping ints in third column of X to gene names in sc_obj\n",
    "            cellcounts, \n",
    "    \"\"\"\n",
    "    levels = np.unique(sc_obj.obs[stratifier])\n",
    "    bin_width = 1/nbins\n",
    "    nL = len(levels)\n",
    "    nG = sc_obj.X.shape[1]\n",
    "    pseudobulk = np.zeros([nbins*nL, nG])\n",
    "    X = np.zeros([nbins*nL*nG, 3])\n",
    "    cellcounts = np.zeros(nbins*nL, dtype=int)\n",
    "    empty = []\n",
    "    for t in range(nbins):\n",
    "        if t < nbins-1:\n",
    "            t_subset = (t*bin_width <= sc_obj.obs[pt_label]) & (sc_obj.obs[pt_label] < (t+1)*bin_width)\n",
    "        else:\n",
    "            t_subset = (t*bin_width <= sc_obj.obs[pt_label]) & (sc_obj.obs[pt_label] <= (t+1)*bin_width)\n",
    "        for l in range(nL):\n",
    "            l_subset = sc_obj.obs[stratifier] == levels[l]\n",
    "            sub = sc_obj[l_subset & t_subset]\n",
    "            # get cell counts for this subset\n",
    "            cellcounts[t*nL + l] = sub.X.shape[0]\n",
    "            # get metadata, store in X\n",
    "            t_sub = np.median(sub.obs[pt_label])\n",
    "            X[(t*nL*nG+l*nG):(t*nL*nG+(l+1)*nG),0] = t_sub\n",
    "            X[(t*nL*nG+l*nG):(t*nL*nG+(l+1)*nG),1] = l\n",
    "            X[(t*nL*nG+l*nG):(t*nL*nG+(l+1)*nG),2] = range(nG)\n",
    "            # get expression data, store in pseudobulk\n",
    "            sub_sum = aggregate(sub)\n",
    "            sub_cpm = log_normalize(sub_sum)\n",
    "            pseudobulk[t*nL + l, :] = sub_cpm\n",
    "            # if there is a sample with zero cells, we'll remember that and delete it \n",
    "            if cellcounts[t*nL + l] == 0:\n",
    "                empty.append(t*nL + l)\n",
    "    if sum(cellcounts==0) > 0:\n",
    "        print(\"Removing \" + str(sum(cellcounts==0)) + \" empty bins\")\n",
    "    X = X[~np.isnan(X[:,0])]\n",
    "    pseudobulk = np.delete(pseudobulk, empty, axis=0)\n",
    "    pseudobulk = scale(pseudobulk)\n",
    "    Y = pseudobulk.flatten()\n",
    "    leveldict = {i:levels[i] for i in range(nL)}\n",
    "    genedict = {sc_obj.var_names[i]:i for i in range(nG)}\n",
    "    \n",
    "    return({\"Y\":Y,\"X\":X,\"level_dict\":leveldict, \"gene_dict\":genedict, \"cell_counts\":cellcounts})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Endothelial Lineage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratify by individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "endo_ind = get_pseudobulk(endo, stratifier=\"individual\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/endo.dpt.ind.X.txt\", endo_ind['X'])\n",
    "np.savetxt(\"../data/endo.dpt.ind.Y.txt\", endo_ind['Y'])\n",
    "with open(\"../data/endo.dpt.ind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(endo_ind['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/endo.dpt.hires.ind.X.txt\", endo_ind['X'])\n",
    "np.savetxt(\"../data/endo.dpt.hires.ind.Y.txt\", endo_ind['Y'])\n",
    "with open(\"../data/endo.dpt.hires.ind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(endo_ind['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratify by batch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "endo_batch = get_pseudobulk(endo, stratifier=\"Batch\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/endo.dpt.batch.X.txt\", endo_batch['X'])\n",
    "np.savetxt(\"../data/endo.dpt.batch.Y.txt\", endo_batch['Y'])\n",
    "with open(\"../data/endo.dpt.batch.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(endo_batch['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/endo.dpt.hires.batch.X.txt\", endo_batch['X'])\n",
    "np.savetxt(\"../data/endo.dpt.hires.batch.Y.txt\", endo_batch['Y'])\n",
    "with open(\"../data/endo.dpt.hires.batch.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(endo_batch['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratify by both batch and individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "endo.obs['batchind'] = [endo.obs['individual'][i] + '--' + endo.obs['Batch'][i] for i in range(endo.X.shape[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3335: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/work-zfs/abattle4/josh/ebpilot/myenv/lib/python3.6/site-packages/ipykernel/__main__.py:6: RuntimeWarning: invalid value encountered in true_divide\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3335: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/work-zfs/abattle4/josh/ebpilot/myenv/lib/python3.6/site-packages/ipykernel/__main__.py:6: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Removing 2 empty bins\n"
     ]
    }
   ],
   "source": [
    "endo_batchind = get_pseudobulk(endo, stratifier=\"batchind\", nbins=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/endo.dpt.batchind.X.txt\", endo_batchind['X'])\n",
    "np.savetxt(\"../data/endo.dpt.batchind.Y.txt\", endo_batchind['Y'])\n",
    "with open(\"../data/endo.dpt.batchind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(endo_batchind['gene_dict'], f)\n",
    "    f.close()\n",
    "with open(\"../data/endo.dpt.batchind.linedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(endo_batchind['level_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/endo.dpt.hires.batchind.X.txt\", endo_batchind['X'])\n",
    "np.savetxt(\"../data/endo.dpt.hires.batchind.Y.txt\", endo_batchind['Y'])\n",
    "with open(\"../data/endo.dpt.hires.batchind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(endo_batchind['gene_dict'], f)\n",
    "    f.close()\n",
    "with open(\"../data/endo.dpt.hires.batchind.linedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(endo_batchind['level_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hepatocyte lineage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratified by individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Removing 1 empty bins\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3335: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/work-zfs/abattle4/josh/ebpilot/myenv/lib/python3.6/site-packages/ipykernel/__main__.py:6: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    }
   ],
   "source": [
    "hep_ind = get_pseudobulk(hep, stratifier=\"individual\", nbins=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/hep.dpt.ind.X.txt\", hep_ind['X'])\n",
    "np.savetxt(\"../data/hep.dpt.ind.Y.txt\", hep_ind['Y'])\n",
    "with open(\"../data/hep.dpt.ind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(hep_ind['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/hep.dpt.hires.ind.X.txt\", hep_ind['X'])\n",
    "np.savetxt(\"../data/hep.dpt.hires.ind.Y.txt\", hep_ind['Y'])\n",
    "with open(\"../data/hep.dpt.hires.ind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(hep_ind['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratified by batch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "hep_batch = get_pseudobulk(hep, stratifier=\"Batch\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/hep.dpt.batch.X.txt\", hep_batch['X'])\n",
    "np.savetxt(\"../data/hep.dpt.batch.Y.txt\", hep_batch['Y'])\n",
    "with open(\"../data/hep.dpt.batch.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(hep_batch['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/hep.dpt.hires.batch.X.txt\", hep_batch['X'])\n",
    "np.savetxt(\"../data/hep.dpt.hires.batch.Y.txt\", hep_batch['Y'])\n",
    "with open(\"../data/hep.dpt.hires.batch.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(hep_batch['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratified by batch and individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "hep.obs['batchind'] = [hep.obs['individual'][i] + '--' + hep.obs['Batch'][i] for i in range(hep.X.shape[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3335: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/work-zfs/abattle4/josh/ebpilot/myenv/lib/python3.6/site-packages/ipykernel/__main__.py:6: RuntimeWarning: invalid value encountered in true_divide\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3335: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/work-zfs/abattle4/josh/ebpilot/myenv/lib/python3.6/site-packages/ipykernel/__main__.py:6: RuntimeWarning: invalid value encountered in true_divide\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3335: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/work-zfs/abattle4/josh/ebpilot/myenv/lib/python3.6/site-packages/ipykernel/__main__.py:6: RuntimeWarning: invalid value encountered in true_divide\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3335: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/work-zfs/abattle4/josh/ebpilot/myenv/lib/python3.6/site-packages/ipykernel/__main__.py:6: RuntimeWarning: invalid value encountered in true_divide\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3335: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/home-2/jpopp4@jhu.edu/.local/lib/python3.6/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/work-zfs/abattle4/josh/ebpilot/myenv/lib/python3.6/site-packages/ipykernel/__main__.py:6: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Removing 5 empty bins\n"
     ]
    }
   ],
   "source": [
    "hep_batchind = get_pseudobulk(hep, stratifier=\"batchind\", nbins=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/hep.dpt.batchind.X.txt\", hep_batchind['X'])\n",
    "np.savetxt(\"../data/hep.dpt.batchind.Y.txt\", hep_batchind['Y'])\n",
    "with open(\"../data/hep.dpt.batchind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(hep_batchind['gene_dict'], f)\n",
    "    f.close()\n",
    "with open(\"../data/hep.dpt.batchind.linedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(hep_batchind['level_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/hep.dpt.hires.batchind.X.txt\", hep_batchind['X'])\n",
    "np.savetxt(\"../data/hep.dpt.hires.batchind.Y.txt\", hep_batchind['Y'])\n",
    "with open(\"../data/hep.dpt.hires.batchind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(hep_batchind['gene_dict'], f)\n",
    "    f.close()\n",
    "with open(\"../data/hep.dpt.hires.batchind.linedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(hep_batchind['level_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neuronal lineage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratified by individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "neur_ind = get_pseudobulk(neur, stratifier=\"individual\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/neur.dpt.ind.X.txt\", neur_ind['X'])\n",
    "np.savetxt(\"../data/neur.dpt.ind.Y.txt\", neur_ind['Y'])\n",
    "with open(\"../data/neur.dpt.ind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(neur_ind['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/neur.dpt.hires.ind.X.txt\", neur_ind['X'])\n",
    "np.savetxt(\"../data/neur.dpt.hires.ind.Y.txt\", neur_ind['Y'])\n",
    "with open(\"../data/neur.dpt.hires.ind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(neur_ind['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratified by batch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "neur_batch = get_pseudobulk(neur, stratifier=\"Batch\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/neur.dpt.batch.X.txt\", neur_batch['X'])\n",
    "np.savetxt(\"../data/neur.dpt.batch.Y.txt\", neur_batch['Y'])\n",
    "with open(\"../data/neur.dpt.batch.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(neur_batch['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/neur.dpt.hires.batch.X.txt\", neur_batch['X'])\n",
    "np.savetxt(\"../data/neur.dpt.hires.batch.Y.txt\", neur_batch['Y'])\n",
    "with open(\"../data/neur.dpt.hires.batch.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(neur_batch['gene_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stratified by batch and individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "neur.obs['batchind'] = [neur.obs['individual'][i] + '--' + neur.obs['Batch'][i] for i in range(neur.X.shape[0])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "neur_batchind = get_pseudobulk(neur, stratifier=\"batchind\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/neur.dpt.batchind.X.txt\", neur_batchind['X'])\n",
    "np.savetxt(\"../data/neur.dpt.batchind.Y.txt\", neur_batchind['Y'])\n",
    "with open(\"../data/neur.dpt.batchind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(neur_batchind['gene_dict'], f)\n",
    "    f.close()\n",
    "with open(\"../data/neur.dpt.batchind.linedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(neur_batchind['level_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt( \"../data/neur.dpt.hires.batchind.X.txt\", neur_batchind['X'])\n",
    "np.savetxt(\"../data/neur.dpt.hires.batchind.Y.txt\", neur_batchind['Y'])\n",
    "with open(\"../data/neur.dpt.hires.batchind.genedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(neur_batchind['gene_dict'], f)\n",
    "    f.close()\n",
    "with open(\"../data/neur.dpt.hires.batchind.linedict.pickle\", 'wb') as f:\n",
    "    pickle.dump(neur_batchind['level_dict'], f)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:myenv]",
   "language": "python",
   "name": "conda-env-myenv-py"
  },
  "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}