run_algs_real_data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import math_utils, sys_utils\n",
    "import numpy as np\n",
    "import os\n",
    "import dci\n",
    "import yaml\n",
    "import pandas as pd\n",
    "\n",
    "# additional imports (for plotting the difference graph in this notebook, not required for running DCI)\n",
    "import networkx as nx\n",
    "import graphviz\n",
    "\n",
    "SAVE_FOLDER = 'data/real_data/'# folder where results will be saved\n",
    "\n",
    "# the inputs below must be present in the correct folder\n",
    "# data corresponding to condition 1 (in this case gene expression of activated CD8+ T-cells)\n",
    "filename1 = 'data/real_data/input/gene_exp_sarkar_act.csv'\n",
    "# data corresponding to condition 2 (in this case gene expression of naive CD8+ T-cells)\n",
    "filename2 = 'data/real_data/input/gene_exp_sarkar_naive.csv'\n",
    "# undirected Markov difference network (estimated via KLIEP https://arxiv.org/pdf/1304.6803.pdf)\n",
    "difference_undirected_filename = 'data/real_data/input/diff_graph_sarkar_center_filt0.005_full.csv'\n",
    "\n",
    "# this notebook will estimate the difference DAG based on the two tuning parameters below for demonstration purposes\n",
    "# in practice the DCI algorithm should be run across a combination of tuning parameters (which was done in the paper)\n",
    "# therefore results based on running this notebook are different from what is shown in the paper\n",
    "DDAG_SKEL_ALPHA = 0.1\n",
    "DDAG_ALPHA = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def determine_changed_nodes(mn_diff):\n",
    "    # determine which nodes have any change from state 1 to state 2\n",
    "    # a node i is included in the changed node set if there exists a node j such that \n",
    "    # precision(i,j) in state 1 is not equal to precision(i,j) in state 2\n",
    "    # markov (undirected) difference graph gives how precision matrix changed across the two states\n",
    "    # i.e a zero for entry (i,j) means that there is no change in precision(i,j) across the two states\n",
    "    return set(np.where(np.sum(mn_diff, axis=0) !=0)[0])\n",
    "    \n",
    "def read_data(center=True):\n",
    "    # read in the data and center the data around 0 if center is True (default)\n",
    "    X1 = pd.read_csv(filename1, delimiter=',', index_col=0).T\n",
    "    X2 = pd.read_csv(filename2, delimiter=',', index_col=0).T\n",
    "    gene_names = X1.columns.values\n",
    "    \n",
    "    if center:\n",
    "        X1 = X1 - X1.mean(axis=0)\n",
    "        X2 = X2 - X2.mean(axis=0)\n",
    "    \n",
    "    # read in markov difference network\n",
    "    mn_diff = np.loadtxt(difference_undirected_filename, delimiter = ',')\n",
    "    # determine which nodes changed\n",
    "    changed_nodes = determine_changed_nodes(mn_diff)\n",
    "    # get all edges with nonzero precision\n",
    "    est_dug = set(math_utils.upper_tri_ixs_nonzero(mn_diff))\n",
    "    return X1.values, X2.values, est_dug, changed_nodes, gene_names\n",
    "\n",
    "def make_graph(skeleton, oriented_edges, gene_names, known_edges=set()):\n",
    "    # create a graph for plotting\n",
    "    unoriented_edges = skeleton - oriented_edges - {(j, i) for i, j in oriented_edges}\n",
    "    # make a directed graph\n",
    "    g = nx.DiGraph()\n",
    "    for i, j in oriented_edges:\n",
    "        color = 'black'\n",
    "        g.add_edge(gene_names[i], gene_names[j], color=color, penwidth=3)\n",
    "\n",
    "    for i, j in unoriented_edges:\n",
    "        color = 'black'\n",
    "        g.add_edge(gene_names[i], gene_names[j], arrowhead='none', color=color, penwidth=3)\n",
    "    \n",
    "    for i, j in known_edges - oriented_edges - unoriented_edges:\n",
    "        if (i, j) not in skeleton and (j, i) not in skeleton:\n",
    "            g.add_edge(gene_names[i], gene_names[j], arrowhead='none', color='gray')\n",
    "    return g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# read in data for DCI algorithm\n",
    "X1, X2, est_dug, changed_nodes, gene_names = read_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimating skeleton of the difference DAG (running Algorithm 2 of DCI)...\n",
      "This might take a bit of time...\n"
     ]
    }
   ],
   "source": [
    "# run DCI algorithm 1 (to limit the number of hypotheses tested limit max_set_size to 3)\n",
    "print('Estimating skeleton of the difference DAG (running Algorithm 2 of DCI)...')\n",
    "print('This might take a bit of time...')\n",
    "retained_edges, _, _, _ = dci.estimate_ddag_skeleton(X1, X2, est_dug, changed_nodes, DDAG_SKEL_ALPHA, max_set_size = 3, verbose=False)\n",
    "# save results\n",
    "RES_FOLDER = SAVE_FOLDER + 'dci_results/'\n",
    "sys_utils.ensure_dirs([RES_FOLDER])\n",
    "yaml.dump(retained_edges, open(RES_FOLDER + 'estimated_ddag_skeleton.yaml', 'w'), indent=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Assigning edge directions (running Algorithm 3 of DCI)...\n",
      "This might take a bit of time...\n"
     ]
    }
   ],
   "source": [
    "# run DCI algorithm 2\n",
    "print('Assigning edge directions (running Algorithm 3 of DCI)...')\n",
    "print('This might take a bit of time...')\n",
    "est_ddag = dci.estimate_ddag(X1, X2, retained_edges, changed_nodes, DDAG_ALPHA, max_set_size = 3, verbose=False)\n",
    "# save results\n",
    "yaml.dump(est_ddag, open(RES_FOLDER + 'estimated_ddag.yaml', 'w'), indent=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plot the graph if desired (requires graphviz, pydot and networkx packages in python)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'data/real_data/dci_results/graph_thres.gv.png'"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('Plot the graph if desired (requires graphviz, pydot and networkx packages in python)')\n",
    "g = make_graph(retained_edges, est_ddag, gene_names)\n",
    "fn = RES_FOLDER + 'graph_thres.gv'\n",
    "nx.nx_pydot.write_dot(g, fn)\n",
    "graphviz.render('dot', 'png', fn)"
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}