{ "cells": [ { "cell_type": "markdown", "id": "aada14c4-3ed3-4e1b-b641-273b38965440", "metadata": {}, "source": [ "# Sachs et al. 2005" ] }, { "cell_type": "code", "execution_count": 1, "id": "cb30bfe1-f462-4aad-9660-c120d5904da5", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/usr/local/Caskroom/miniconda/base/envs/cstrees/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n", " from .autonotebook import tqdm as notebook_tqdm\n" ] } ], "source": [ "import random\n", "import time\n", "\n", "from causallearn.search.ConstraintBased.PC import pc\n", "from causallearn.utils.cit import chisq\n", "\n", "import pandas as pd\n", "import cslearn.scoring as sc\n", "import cslearn.learning as ctl\n", "import numpy as np\n", "\n", "%load_ext autoreload\n", "%autoreload 2" ] }, { "cell_type": "markdown", "id": "0c255490", "metadata": {}, "source": [ "## Read data" ] }, { "cell_type": "code", "execution_count": 2, "id": "3baad6d3-2e1c-48ec-9c94-77815df91db3", "metadata": {}, "outputs": [], "source": [ "sachs = pd.read_csv('../data/sachs_obs_data.csv')" ] }, { "cell_type": "markdown", "id": "395b3bc0", "metadata": {}, "source": [ "As a pre-processing step, we make all variables binary." ] }, { "cell_type": "code", "execution_count": 3, "id": "e99fc198-0ccc-4aaf-8f0d-673283d5775d", "metadata": {}, "outputs": [], "source": [ "# Binarize data and add a row at the top of the dataframe containing variable cardinalities\n", "sachsnp = sachs.to_numpy()\n", "\n", "sachs2 = np.zeros([len(sachs), len(list(sachs.columns))], int)\n", "\n", "for i in range(len(list(sachs.columns))):\n", " sachs2[:,i] = pd.cut(sachsnp[:,i], 2, labels = False)\n", " \n", "# add row with state space cardinalities\n", "sachs2states = np.zeros([len(sachs) + 1, len(list(sachs.columns))], int)\n", "sachs2states[0,:] = [2 for i in range(len(list(sachs.columns)))]\n", "for i in range(len(sachs)):\n", " sachs2states[i+1,:] = sachs2[i,:]\n", " \n", "sachsdf = pd.DataFrame(sachs2states, columns = list(sachs.columns))" ] }, { "cell_type": "markdown", "id": "74bbbcc9", "metadata": {}, "source": [ "## MCMC sampling" ] }, { "cell_type": "markdown", "id": "f0c1f2ad", "metadata": {}, "source": [ "We run without restricting the parent sets according to a CPDAG (poss_cvars=None).\n", "To run with this restriction uncomment the commented lines of code and let the poss_cvars=poss_cvars in sc.order_score_tables." ] }, { "cell_type": "code", "execution_count": 4, "id": "45e65886-a4fe-456b-ac11-2bc036ce17da", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Context score tables: 100%|██████████| 11/11 [00:00<00:00, 17.21it/s]\n", "Creating #stagings tables: 100%|██████████| 11/11 [00:00<00:00, 197.15it/s]\n", "Order score tables: 100%|██████████| 11/11 [01:12<00:00, 6.58s/it]\n", "Gibbs order sampler: 100%|██████████| 5000/5000 [00:00<00:00, 7991.33it/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Computation time in seconds: 73.76801800727844\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "# pcgraph = pc(sachsdf[1:].values, 0.9, \"chisq\", node_names=sachsdf.columns)\n", "# poss_cvars = ctl.causallearn_graph_to_posscvars(pcgraph, labels=sachsdf.columns)\n", "# print(\"Possible context variables per node:\", poss_cvars)\n", "np.random.seed(1)\n", "random.seed(1)\n", "start = time.time()\n", "score_table, context_scores, context_counts = sc.order_score_tables(sachsdf,\n", " max_cvars=2,\n", " alpha_tot=1.0,\n", " method=\"BDeu\",\n", " poss_cvars=None)\n", "\n", "orders, scores = ctl.gibbs_order_sampler(5000, score_table)\n", "end = time.time()\n", "print('Computation time in seconds:', end - start)" ] }, { "cell_type": "markdown", "id": "4ed66998", "metadata": {}, "source": [ "Optimal variable ordering" ] }, { "cell_type": "code", "execution_count": 5, "id": "4a50b757-8366-4566-944c-bac9210e0262", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['Raf', 'PIP3', 'PLCg', 'PIP2', 'PKA', 'JNK', 'PKC', 'Akt', 'Erk', 'Mek', 'p38']\n" ] } ], "source": [ "sachsmaporder = orders[scores.index(max(scores))]\n", "print(sachsmaporder)" ] }, { "cell_type": "markdown", "id": "2d1612d2", "metadata": {}, "source": [ "Get optimal tree for ordering" ] }, { "cell_type": "code", "execution_count": 6, "id": "c318a7d9-7926-4dc6-a8a7-d41d295f01d3", "metadata": {}, "outputs": [], "source": [ "sachsopttree = ctl._optimal_cstree_given_order(sachsmaporder, context_scores)" ] }, { "cell_type": "markdown", "id": "c57dc9ac", "metadata": {}, "source": [ "Dataframe representation" ] }, { "cell_type": "code", "execution_count": 7, "id": "08d28dcc", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
RafPIP3PLCgPIP2PKAJNKPKCAktErkMekp38
022222222222
1*----------
2**---------
3*0*--------
4*10--------
5*11--------
6***0-------
70**1-------
81**1-------
9**0**------
10**1**------
11******-----
12****0**----
13****1**----
14*******0---
15*******1---
160********--
171********--
18******0***-
19******1*0*-
20******1*1*-
21-----------
\n", "
" ], "text/plain": [ " Raf PIP3 PLCg PIP2 PKA JNK PKC Akt Erk Mek p38\n", "0 2 2 2 2 2 2 2 2 2 2 2\n", "1 * - - - - - - - - - -\n", "2 * * - - - - - - - - -\n", "3 * 0 * - - - - - - - -\n", "4 * 1 0 - - - - - - - -\n", "5 * 1 1 - - - - - - - -\n", "6 * * * 0 - - - - - - -\n", "7 0 * * 1 - - - - - - -\n", "8 1 * * 1 - - - - - - -\n", "9 * * 0 * * - - - - - -\n", "10 * * 1 * * - - - - - -\n", "11 * * * * * * - - - - -\n", "12 * * * * 0 * * - - - -\n", "13 * * * * 1 * * - - - -\n", "14 * * * * * * * 0 - - -\n", "15 * * * * * * * 1 - - -\n", "16 0 * * * * * * * * - -\n", "17 1 * * * * * * * * - -\n", "18 * * * * * * 0 * * * -\n", "19 * * * * * * 1 * 0 * -\n", "20 * * * * * * 1 * 1 * -\n", "21 - - - - - - - - - - -" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sachsopttree.to_df()" ] }, { "cell_type": "markdown", "id": "b13dd013", "metadata": {}, "source": [ "## LDAG representation" ] }, { "cell_type": "markdown", "id": "5b40c1cd", "metadata": {}, "source": [ "Print the LDAG representation of the learned CStree model" ] }, { "cell_type": "code", "execution_count": 8, "id": "8ecd8ffb-d221-4a1b-b1c0-958799a7b1cd", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Raf\n", "\n", "Raf\n", "\n", "\n", "\n", "PKA\n", "\n", "PKA\n", "\n", "\n", "\n", "Raf->PKA\n", "\n", "\n", "[[0]]\n", "\n", "\n", "\n", "Mek\n", "\n", "Mek\n", "\n", "\n", "\n", "Raf->Mek\n", "\n", "\n", "\n", "\n", "\n", "Akt\n", "\n", "Akt\n", "\n", "\n", "\n", "PKA->Akt\n", "\n", "\n", "\n", "\n", "\n", "Erk\n", "\n", "Erk\n", "\n", "\n", "\n", "Akt->Erk\n", "\n", "\n", "\n", "\n", "\n", "p38\n", "\n", "p38\n", "\n", "\n", "\n", "Erk->p38\n", "\n", "\n", "[[0]]\n", "\n", "\n", "\n", "PIP3\n", "\n", "PIP3\n", "\n", "\n", "\n", "PIP2\n", "\n", "PIP2\n", "\n", "\n", "\n", "PIP3->PIP2\n", "\n", "\n", "\n", "\n", "\n", "PIP2->PKA\n", "\n", "\n", "\n", "\n", "\n", "PLCg\n", "\n", "PLCg\n", "\n", "\n", "\n", "PLCg->PIP2\n", "\n", "\n", "[[0]]\n", "\n", "\n", "\n", "JNK\n", "\n", "JNK\n", "\n", "\n", "\n", "PLCg->JNK\n", "\n", "\n", "\n", "\n", "\n", "PKC\n", "\n", "PKC\n", "\n", "\n", "\n", "PKC->p38\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ ">" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "LDAG = sachsopttree.to_LDAG() \n", "agraph = LDAG.plot_graphviz(args=\"-Goneblock=True\",\n", " prog=\"circo\")\n", "agraph \n", "#agraph.draw('sachs_CStree_LDAG.png')\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.3" } }, "nbformat": 4, "nbformat_minor": 5 }