CytoscapeReactomeFI.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cytoscape Reactome FI Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Overview\n",
    "--------\n",
    "This example demonstrates a workflow using ReactomeFIViz CyREST API to perform a comparison pathway and network analysis for TCGA BRCA (breast invasive carcinoma) and OV (ovarian serous cystadenocarcinoma) mutation data.\n",
    "\n",
    "The steps executed are as follows:\n",
    "\n",
    "1. Download each cohort mutation annotation file (MAF) from the Broad Firehose WS server.\n",
    "\n",
    "2. Build the two networks using the ReactomeFIViz CyREST API buildFISubNetwork resource.\n",
    "\n",
    "3. Cluster the network into network modules using the cluster resource.\n",
    "\n",
    "4. Perform pair-wise module overlapping analysis.\n",
    "\n",
    "5. Load the Reactome pathways.\n",
    "\n",
    "6. Partition the network modules generated from two networks into the following sets:\n",
    "\n",
    "  * Modules with a significant proportion of genes in common\n",
    "\n",
    "  * Unshared modules for each of the two cancer types\n",
    "\n",
    "7. Perform pathway enrichment analysis for each unshared gene module.\n",
    "\n",
    "8. Export the diagram for an unshared pathway for each cancer type.\n",
    "\n",
    "Prerequisites\n",
    "------------\n",
    "* Python version 2.6 or later.\n",
    "\n",
    "* Python pandas package version 0.22.0 or later.\n",
    "\n",
    "* The latest version of Cytoscape with the latest version of the\n",
    "  ReactomeFIViz app.\n",
    "\n",
    "Installation\n",
    "-----------\n",
    "* Install the required packages in a Python environment, e.g. in the\n",
    "  [Anaconda](https://docs.anaconda.com/anaconda/) virtual environment:\n",
    "\n",
    "      conda create -n cytoscape pandas scipy statsmodels\n",
    "      source activate cytoscape\n",
    "\n",
    "* Clone the example GitHub repository, e.g.:\n",
    "  \n",
    "      mkdir ~/reactome-fi\n",
    "      cd ~/reactome-fi\n",
    "      git clone https://github.com/reactome-fi/workflows.git\n",
    "      cd workflows\n",
    "\n",
    "* In Cytoscape, select `Help > Check for Updates` to ensure that\n",
    "  ReactomeFIViz is current.\n",
    "\n",
    "Usage\n",
    "-----\n",
    "* Start Cytoscape.\n",
    "\n",
    "* In a console, start Jupyter on the notebook directory, e.g.:\n",
    "\n",
    "      ipython notebook --notebook-dir=~/reactome-fi/workflows\n",
    "\n",
    "* Run all notebook cells in the web browser Jupyter window.\n",
    "\n",
    "*Note*: Building the networks clears an existing Cytoscape session,\n",
    "if it exists."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "\n",
    "**Utility functions.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Python utilities.\n",
    "import sys\n",
    "import os\n",
    "import re\n",
    "from functools import reduce\n",
    "import requests\n",
    "import itertools\n",
    "import pandas as pd\n",
    "from scipy.stats import hypergeom\n",
    "from statsmodels.sandbox.stats.multicomp import multipletests\n",
    "from IPython.display import (display, HTML)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Constants.\n",
    "\n",
    "# The number of background genes N is the Reactome FI\n",
    "# network gene size.\n",
    "BACKGROUND_GENE_CNT = 12283\n",
    "\n",
    "# The REST service url.\n",
    "SERVICE_URL = 'http://localhost:1234'\n",
    "\n",
    "# The default minimum cluster module size.\n",
    "DEF_MIN_MODULE_SIZE = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_url(*params, **kwargs):\n",
    "    \"\"\"\n",
    "    :param params: the URL path component strings\n",
    "    :option app: the CyREST application name\n",
    "    :return: the REST app URL\n",
    "    \"\"\"\n",
    "    path = [SERVICE_URL]\n",
    "    app = kwargs.get('app')\n",
    "    if app is not None:\n",
    "        path.append(app)\n",
    "    path.append('v1')\n",
    "    path.extend(params)\n",
    "    return '/'.join(path)\n",
    "\n",
    "# Returns the CyREST Reactome FI url.\n",
    "get_fi_url = lambda *params: get_url(*params, app='reactomefiviz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def parse_fi_table_response(resp, parsers, index=None):\n",
    "    \"\"\"\n",
    "    Returns the data frame for the given CyREST Reactome\n",
    "    FI response. The response `data` property object must\n",
    "    be a JSON object with properties *tableHeaders* and\n",
    "    *tableContent*. If the response `data` is empty, then\n",
    "    this function returns an empty data frame with columns\n",
    "    given by the *parsers* keys.\n",
    "    \n",
    "    The required *parsers* dictionary argument associates\n",
    "    a parser with each column.\n",
    "    \n",
    "    :param resp: the CyREST response\n",
    "    :param parsers: the column parser dictionary\n",
    "    :option index: the index column name\n",
    "    :return: the parsed data frame\n",
    "    \"\"\"\n",
    "    # The response JSON data object.\n",
    "    data = resp.json()['data']\n",
    "    # The data columns.\n",
    "    columns = data.get('tableHeaders') if data else parsers.keys()\n",
    "    parsers_list = [parsers[col] for col in columns]\n",
    "    # Parses a content row.\n",
    "    parse_row = lambda row: tuple(parsers_list[i](value)\n",
    "                                  for i, value in enumerate(row))\n",
    "    # The parsed content list.\n",
    "    content = map(parse_row, data['tableContent']) if data else []\n",
    "    # Return the data frames.\n",
    "    return pd.DataFrame.from_records(content, index=index, columns=columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def download_maf(cohort, out=None):\n",
    "    \"\"\"\n",
    "    Downloads the MAF file for the given cancer type cohort.\n",
    "    The default output file is `<cohort>.maf` in the\n",
    "    sandbox directory.\n",
    "    \n",
    "    *Note*: As of April 2018, the Firebrowse server is unstable.\n",
    "    Calling this method is not recommended until the server\n",
    "    stabilizes.\n",
    "    \n",
    "    :param cohort: the cancer type cohort name\n",
    "    :option out: the optional output file path\n",
    "    :return: the output file path\n",
    "    \"\"\"\n",
    "    # Download the MAF.\n",
    "    url = 'http://firebrowse.org/api/v1/Analyses/Mutation/MAF'\n",
    "    maf_file = out if out else os.path.join(sandbox, \"%s.maf\" % cohort) \n",
    "    print(\"Downloading the %s MAF file to %s...\" % (cohort, maf_file))\n",
    "    params = dict(format='tsv', cohort=cohort, page_size=200)\n",
    "    eof = False\n",
    "    page = 1\n",
    "    with open(maf_file, 'w') as f:\n",
    "        while not eof:\n",
    "            params['page'] = page\n",
    "            resp = requests.get(url, params=params)\n",
    "            if not resp.ok:\n",
    "                print(\"Error encountered downloading the %s MAF file. Please retry.\" %\n",
    "                      cohort)\n",
    "                resp.raise_for_status()\n",
    "            text = resp.text\n",
    "            if text:\n",
    "                f.write(text)\n",
    "                page = params['page'] = page + 1\n",
    "            else:\n",
    "                eof = True\n",
    "            print('+', end='')\n",
    "    print('')\n",
    "    print(\"MAF file downloaded.\")\n",
    "    return maf_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_network(maf_file, cutoff=.01):\n",
    "    \"\"\"\n",
    "    Builds the network from the given MAF file.\n",
    "    Only gene modules whose sample size exceeds a\n",
    "    threshold are included. The threshold is the\n",
    "    greater of 2 and sample count * cutoff, where\n",
    "    sample count is the total number of MAF samples.\n",
    "    \n",
    "    :param maf_file: the downloaded MAF file\n",
    "    :option cutoff: the minimum proportion of samples\n",
    "    \"\"\"\n",
    "    # Clear the current Cytoscape session, if any.\n",
    "    requests.delete(get_url('session')).ok\n",
    "    # Read the MAF into a data frame.\n",
    "    maf_df = pd.read_csv(maf_file, sep='\\t',\n",
    "                         usecols=['Tumor_Sample_Barcode'])\n",
    "    # The number of samples.\n",
    "    sample_cnt = len({tsb for tsb in maf_df.Tumor_Sample_Barcode})\n",
    "    print(\"Sample Count: %d\" % sample_cnt)\n",
    "    sample_cutoff = max(2, int(0.01*sample_cnt))\n",
    "    print(\"Building the FI network with sample cut-off %d...\" %\n",
    "          sample_cutoff)\n",
    "    # Build the FI network with 1% sample cut-off.\n",
    "    body = dict(fiVersion='2016', format='MAF', file=maf_file,\n",
    "                enteredGenes='', chooseHomoGenes=False,\n",
    "                userLinkers=False, showUnLinked=False,\n",
    "                fetchFIAnnotations=True,\n",
    "                sampleCutoffValue=sample_cutoff)\n",
    "    requests.post(get_fi_url('buildFISubNetwork'), json=body)\n",
    "    print(\"The FI network is loaded to Cytoscape.\")\n",
    "    return maf_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The cluster content value parsers.\n",
    "cluster_parsers = {\n",
    "    'Module': int, 'Nodes in Module': int, 'Node Percentage': float,\n",
    "    'Samples in Module': int, 'Sample Percentage': float,\n",
    "    'Node List': lambda s: s.split(',')\n",
    "}\n",
    "\n",
    "def cluster():\n",
    "    \"\"\"\n",
    "    Cluster the network currently displayed in Cytoscape into\n",
    "    gene modules.\n",
    "    \n",
    "    :return: a data frame with index *Module* and columns\n",
    "      *module_size* and *genes*\n",
    "    \"\"\"\n",
    "    # Cluster the currently displayed network.\n",
    "    print(\"Clustering the FI network...\")\n",
    "    resp = requests.get(get_fi_url('cluster'))\n",
    "    # Parse the response JSON into a data frame.\n",
    "    parsed = parse_fi_table_response(resp, parsers=cluster_parsers,\n",
    "                                     index='Module')\n",
    "    # Rename the columns.\n",
    "    rename_opts = {'Node List': 'Genes', 'Nodes in Module': 'Module Size'}\n",
    "    renamed = parsed.rename(columns=rename_opts)\n",
    "    # Retain only the genes and size columns.\n",
    "    sliced = renamed.loc[:, ['Module Size', 'Genes']]\n",
    "    print(\"The cluster table is available in Cytoscape.\")\n",
    "    return sliced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prepare(cohort):\n",
    "    \"\"\"\n",
    "    Prepares the given cancer type cohort for analysis.\n",
    "    \n",
    "    :param cohort: the cancer type cohort name\n",
    "    :return: the clustered data frame\n",
    "    \"\"\"\n",
    "    # Clear the current Cytoscape session, if any.\n",
    "    try:\n",
    "        requests.delete(get_url('session')).ok\n",
    "    except requests.exceptions.ConnectionError:\n",
    "        print(\"Error: Connection refused: is Cytoscape started?\")\n",
    "        raise StandardError()\n",
    "    maf_file = os.path.join(sandbox, \"%s.maf\" % cohort)\n",
    "    # Download the MAF, if necessary.\n",
    "    if (not os.path.exists(maf_file)):\n",
    "        # Download disabled due to Firebrowse server instability.\n",
    "        #download_maf(cohort, out=maf_file)\n",
    "        raise IOError(\"The %s MAF file was not found: %s\" %\n",
    "                      (cohort, maf_file)) \n",
    "    # Build the network.\n",
    "    build_network(maf_file)\n",
    "    # Cluster into gene modules.\n",
    "    return cluster()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def overlap_pvalue(N, overlap, gs1, gs2):\n",
    "    \"\"\"\n",
    "    :param N: the number of background genes\n",
    "    :param overlap: the genes in common\n",
    "    :param gs1: the first gene set\n",
    "    :param gs2: the second gene set\n",
    "    :return: the probability of gene set overlap from\n",
    "      a background population of *N* genes\n",
    "    \"\"\"\n",
    "    # The hypergeometric distribution parameters.\n",
    "    k = len(overlap) - 1\n",
    "    n1 = len(gs1)\n",
    "    n2 = len(gs2)\n",
    "    # Return the overlap probability.\n",
    "    return hypergeom.sf(k, N, n1, n2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def filter_on_module_size(inputs, **kwargs):\n",
    "    \"\"\"\n",
    "    Returns subsets of the given data frames with module\n",
    "    count no greater than a threshold value. The threshold\n",
    "    is given by the *max_module_count* option, with default\n",
    "    the minimum module count of the input data frames.\n",
    "    The module selection criterion is a gene module size\n",
    "    which ensures that no more than *max_module_count*\n",
    "    modules are selected from each data frame.\n",
    "    \n",
    "    :param inputs: the input {name: data frame} dictionary\n",
    "    :param kwargs: the following options:\n",
    "    :option min_module_size: the minimum module size (default 3)\n",
    "    :option max_module_count: the maximum number of modules\n",
    "    :return: the corresponding filtered data frames\n",
    "    \"\"\"\n",
    "    dfs = inputs.values()\n",
    "    # The size cutoff.\n",
    "    max_module_cnt_opt = kwargs.get('max_module_count')\n",
    "    min_df_size = min(df.index.size for df in dfs)\n",
    "    if max_module_cnt_opt:\n",
    "        n = min(max_module_cnt_opt, min_df_size)\n",
    "    else:\n",
    "        n = min_df_size\n",
    "    min_module_size_opt = kwargs.get('min_module_size')\n",
    "    if min_module_size_opt:\n",
    "        min_module_size = min_module_size_opt\n",
    "    else:\n",
    "        min_module_size = DEF_MIN_MODULE_SIZE\n",
    "\n",
    "    # The n largest data frames.\n",
    "    largest = [df['Module Size'].nlargest(n) for df in dfs]\n",
    "    sizes = [ds.iloc[n - 1] for ds in largest]\n",
    "    cutoff = max(min_module_size, *sizes)\n",
    "    return {name: df.loc[df['Module Size'].ge(cutoff), :]\n",
    "            for name, df in inputs.items()}\n",
    "\n",
    "def slice_on_node_list(name, df):\n",
    "    \"\"\"\n",
    "    Slices the given data frame vertically on the *genes* column.\n",
    "    Renames the *Module* index to the given name and the *genes*\n",
    "    column to '<name> Genes', e.g. 'BRCA Genes' for *name* 'BRCA'.\n",
    "    Adds a column '<name> Size' with the module size.\n",
    "    Enriches the genes and adds a column'\n",
    "    \n",
    "    :param name: the cancer type cohort name\n",
    "    :param df: the gene modules data frame\n",
    "    :return: the data frame with renamed index and column\n",
    "    \"\"\"\n",
    "    rename_dict = {'Genes': \"%s Genes\" % name}\n",
    "    sliced_df = df.loc[:, ['Genes']].rename(columns=rename_dict)\n",
    "    sliced_df.index.names = [name]\n",
    "    size_col = \"%s Size\" % name\n",
    "    size = df.Genes.apply(len)\n",
    "    return sliced_df.assign(**{size_col: size})\n",
    "\n",
    "def cartesian(df1, df2):\n",
    "    \"\"\"\n",
    "    Takes the cartesian product of the input data frames.\n",
    "    The resulting data frame has a mult-index with levels\n",
    "    given by the respective input data frame indexes.\n",
    "    \n",
    "    :param df1: the first data frame\n",
    "    :param df2: the second data frame\n",
    "    :return: the cartesian product of the input data frames\n",
    "    \"\"\"\n",
    "    rows = itertools.product(df1.iterrows(), df2.iterrows())\n",
    "    values = (left.append(right) for (_, left), (_, right) in rows)\n",
    "    indexes = [df1.index, df2.index]\n",
    "    index_names = [index.name for index in indexes]\n",
    "    index_values = [index.values for index in indexes]\n",
    "    multi_index = pd.MultiIndex.from_product(index_values, names=index_names)\n",
    "    return pd.DataFrame(values, index=multi_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def append_overlap(n, cross_df):\n",
    "    \"\"\"\n",
    "    Calculates the pair-wise overlap p-value for the given\n",
    "    gene module cross-product.\n",
    "    The return value is the cross-product augmented with\n",
    "    two columns:\n",
    "    * _p-value_ - the hypergeometric test p value\n",
    "    * _FDR_ - the FDR corrected for multiple tests\n",
    "    \n",
    "    :param n: the background gene count\n",
    "    :param cross_df: the cross-product data frame\n",
    "    :return: the augmented data frame\n",
    "    \"\"\"\n",
    "    # Determines the shared genes.\n",
    "    intersect = lambda gs1, gs2: set(gs1).intersection(set(gs2))\n",
    "    # Select only the gene columns.\n",
    "    gene_cols = [col for col in cross_df.columns if col.endswith(\"Genes\")]\n",
    "\n",
    "    # The gene set intersections for each row.\n",
    "    genesets = lambda row: [row[col] for col in gene_cols]\n",
    "    shared = [intersect(*genesets(row))\n",
    "              for _, row in cross_df.iterrows()]\n",
    "    shared_sizes = [len(gs) for gs in shared]\n",
    "    # Augment the cross-product with the shared genes and gene count.\n",
    "    kwargs = {'Shared Genes': shared, 'Shared Size': shared_sizes}\n",
    "    shared_df = cross_df.assign(**kwargs)\n",
    "\n",
    "    # Augment the shared data frame with the pvalues.\n",
    "    pvalue = lambda row: overlap_pvalue(BACKGROUND_GENE_CNT,\n",
    "                                        row['Shared Genes'],\n",
    "                                        *genesets(row))\n",
    "    pvals = shared_df.agg(pvalue, axis=1)\n",
    "    overlap_df = shared_df.assign(**{'p-value': pvals})\n",
    "\n",
    "    # Correct the p-values for multiple comparison hypothesis\n",
    "    # testing by applying the Benjamini–Hochberg FDR procedure.\n",
    "    _, corrected, _, _ = multipletests(pvals.values, method='fdr_bh')\n",
    "    overlap_df['FDR'] = corrected\n",
    "    return overlap_df\n",
    "\n",
    "def partition_shared(overlap_df, inputs, cutoff=.01):\n",
    "    \"\"\"\n",
    "    Partitions the gene modules into three data series:\n",
    "    * A shared data series for gene modules with significant\n",
    "      overlap\n",
    "    * Unshared data series for each cancer type\n",
    "    The shared data series has a multi-index of the cancer types.\n",
    "    The unshared data series has index `Module`.\n",
    "    Each data series values are gene lists.\n",
    "    The gene modules are filtered by p-value (see *cutoff* option below).\n",
    "    \n",
    "    :param overlap_df: the input data frame is the BRCA x OV\n",
    "      module cartesian cross-product augmented with a pair-wise\n",
    "      overlap *FDR*\n",
    "    :param inputs: the {name: <data frame>} inputs\n",
    "    :option cutoff: the corrected p-value cutoff value which\n",
    "      determines a shared module (default .01)\n",
    "    :return the {'shared': <series>, 'unshared': {name: <series>}}\n",
    "      data series dictionary for each *name* in inputs\n",
    "    \"\"\"\n",
    "    # The shared gene modules have FDR <= cutoff\n",
    "    criterion = overlap_df['FDR'].le(cutoff)\n",
    "    shared_df = overlap_df.loc[criterion]\n",
    "\n",
    "    all_modules = [set(values)\n",
    "                   for values in zip(*overlap_df.index.values)]\n",
    "    shared_modules = [set(values)\n",
    "                      for values in zip(*shared_df.index.values)]\n",
    "    unshared_modules = [all_modules[i].difference(shared_modules[i])\n",
    "                        for i in range(len(all_modules))]\n",
    "    unshared_idx_dict = {name: unshared_modules[i]\n",
    "                         for i, name in enumerate(overlap_df.index.names)}\n",
    "    unshared = {name: df.loc[unshared_idx_dict[name], :]\n",
    "                for name, df in inputs.items()}\n",
    "    for name, df in unshared.items():\n",
    "        print(\"Unshared %s module count: %d\" % (name, df.index.size))\n",
    "\n",
    "    return dict(shared=shared_df, unshared=unshared)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def discover_overlap(n, inputs, **kwargs):\n",
    "    \"\"\"\n",
    "    Performs pair-wise overlap analysis on the given cluster\n",
    "    data frames.\n",
    "    \n",
    "    :param n: the background gene count\n",
    "    :param inputs: the {name: data frame} inputs dictionary\n",
    "    :param kwargs: the filter_on_module_size options \n",
    "    :return the {'shared': <series>, 'unshared': {name: <series>}}\n",
    "      data series dictionary for each *name* in inputs\n",
    "    \"\"\"\n",
    "    # Select the largest modules.\n",
    "    filtered = filter_on_module_size(inputs, **kwargs)\n",
    "    # Prep pair-wise analysis by taking the cartesian product\n",
    "    # of the two filtered data frames renamed with unique columns.\n",
    "    sliced = {name: slice_on_node_list(name, df)\n",
    "              for name, df in filtered.items()}\n",
    "    cross_df = cartesian(*sliced.values())\n",
    "    overlap_df = append_overlap(n, cross_df)\n",
    "    return overlap_df\n",
    "\n",
    "def print_overlap(overlap_df, cutoff=None, cancer=None, module=None, show_genes=False, format='html'):\n",
    "    \"\"\"\n",
    "    Represents the given overlap data frame as a data frame\n",
    "    suitable for printing. If both cancer and module are specificed \n",
    "    only the required rows will be output.\n",
    "    \n",
    "    :param overlap_df: the raw overlap data frame\n",
    "    :option cutoff: the FDR cut-off value\n",
    "    :option cancer: one of BRCA or OV\n",
    "    :option module: module index for the specified cancer\n",
    "    :option show_genes: control if genes should be output\n",
    "    :option format: 'html' or 'text' (default is html)\n",
    "    :return: a data frame suitable for printing\n",
    "    \"\"\"\n",
    "    # Filter the overlap FDR. Set this to a lower value to\n",
    "    # Restrict the table size.\n",
    "    if cutoff:\n",
    "        filtered_df = overlap_df[overlap_df.FDR.le(cutoff)]\n",
    "    else:\n",
    "        filtered_df = overlap_df\n",
    "    rename_opts = {name: \"%s Module\" % name\n",
    "                   for name in filtered_df.index.names}\n",
    "    flat_df = filtered_df.reset_index().rename(rename_opts, axis=1)\n",
    "\n",
    "    if show_genes :\n",
    "        columns = [col for col in flat_df.columns]\n",
    "    else :\n",
    "        columns = [col for col in flat_df.columns if not col.endswith(\"Genes\")]\n",
    "    cohort_col_grps = [[col for col in columns if col.startswith(name)]\n",
    "                       for name in filtered_df.index.names]\n",
    "    cohort_cols = reduce(lambda x,y: x + y, cohort_col_grps)\n",
    "    non_cohort_cols = [col for col in columns if col not in cohort_cols]\n",
    "    ordered_cols = cohort_cols + non_cohort_cols\n",
    "    printable_df = flat_df.loc[:, ordered_cols]\n",
    "    \n",
    "    # Perform further filter for output\n",
    "    filter = None\n",
    "    if (cancer and module) :\n",
    "        if (cancer == 'BRCA') :\n",
    "            filter = printable_df['BRCA Module'].eq(module)\n",
    "        elif (cancer == 'OV') :\n",
    "            filter = printable_df['OV Module'].eq(module)\n",
    "    if filter is not None :\n",
    "        printable_df = printable_df[filter]\n",
    "    \n",
    "    if format == 'text':\n",
    "        print(printable_df.to_string(index=False))\n",
    "    elif format == 'html':\n",
    "        start = '<h4>Module Overlap'\n",
    "        if cutoff:\n",
    "            style = 'font-size:normal;font-weight:normal;'\n",
    "            middle = \"<span style=%s> (FDR <= %s)</span>\" % (style, cutoff)\n",
    "        else:\n",
    "            middle = ''\n",
    "        end = '</h4>'\n",
    "        heading = start + middle + end\n",
    "        display(HTML(heading))\n",
    "        display(HTML(printable_df.to_html(index=False)))\n",
    "    else:\n",
    "        raise ValueError(\"Unrecognized overlap print format: %s\" % format)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "enrichment_parsers = {\n",
    "    'ReactomePathway': lambda p: p,\n",
    "    'RatioOfProteinInPathway': float,\n",
    "    'NumberOfProteinInPathway': int,\n",
    "    'ProteinFromGeneSet': int,\n",
    "    'P-value': float,\n",
    "    'FDR': float,\n",
    "    'HitGenes': lambda s: s.split(',')\n",
    "}\n",
    "\n",
    "def enrich_genes(genes, cutoff=.001):\n",
    "    \"\"\"\n",
    "    Perform pathway enrichment analysis on the given genes.\n",
    "    The resulting data frame has index *Pathway* and\n",
    "    columns *p-value* and *FDR*\n",
    "    \n",
    "    :param genes: the gene list or set to enrich\n",
    "    :option cutoff: the FDR cut-off (default .001)\n",
    "    :return: the result data frame, or an empty data frame\n",
    "      if the genes could not be enriched\n",
    "    \"\"\"\n",
    "    data = ','.join(genes)\n",
    "    # Perform the Reactome enrichment analysis.\n",
    "    resp = requests.post(get_fi_url('ReactomePathwayEnrichment'),\n",
    "                         data=data)\n",
    "    resp.raise_for_status()\n",
    "    parsed_df = parse_fi_table_response(resp, enrichment_parsers,\n",
    "                                        index='ReactomePathway')\n",
    "    # Rename the index and P-value column for consistency.\n",
    "    parsed_df.index.rename('Pathway', inplace=True)\n",
    "    renamed_df = parsed_df.rename({'P-value': 'p-value'}, axis=1)\n",
    "    # We only want the p-value and FDR.\n",
    "    sliced_df = renamed_df.loc[:, ['p-value', 'FDR']]\n",
    "    # Sort by FDR.\n",
    "    sorted_df = sliced_df.sort_values(by='FDR')\n",
    "    # Apply the cutoff, if it exists.\n",
    "    return sorted_df[sorted_df.FDR.le(cutoff)] if cutoff else sorted_df\n",
    "\n",
    "def enrich_modules(name, genes_ds, **kwargs):\n",
    "    \"\"\"\n",
    "    Perform pathway enrichment on the given gene modules.\n",
    "    \n",
    "    :param name: the cohort name\n",
    "    :param genes_ds: the data series with module number index\n",
    "      and gene list value\n",
    "    :param kwargs: the `enrich_genes` options\n",
    "    :return: a data frame with Module index and columns\n",
    "      *Genes* and *Pathways*, where each *Pathways* value\n",
    "      is the enrichment result data frame\n",
    "    \"\"\"\n",
    "    print(\"Enriching %d %s modules...\" % (genes_ds.index.size, name))\n",
    "    # Gets the enriched pathways.\n",
    "    enrich = lambda genes: enrich_genes(genes, **kwargs)\n",
    "    pathways = genes_ds.apply(enrich)\n",
    "    print(\"%s enrichment completed.\" % name)\n",
    "    return genes_ds.to_frame().assign(Pathways=pathways)\n",
    "\n",
    "def _accum_pathways(accum, index, row):\n",
    "    pathways_df = row.Pathways\n",
    "    pathways = pathways_df.index.get_values()\n",
    "    dup_module = [index] * pathways_df.index.size\n",
    "    multi_tuples = list(zip(pathways, dup_module))\n",
    "    multi_names = ('Pathway', 'Module')\n",
    "    multi_index = pd.MultiIndex.from_tuples(multi_tuples, names=multi_names)\n",
    "    multi_values = pathways_df.values #list(zip(*pathways_df.values))\n",
    "    multi_cols = pathways_df.columns.get_values()\n",
    "    multi_df = pd.DataFrame(data=multi_values, columns=multi_cols,\n",
    "                            index=multi_index)\n",
    "    return multi_df if accum.empty else accum.append(multi_df)\n",
    "\n",
    "def transpose_pathways(enrichment_df):\n",
    "    accumulate = lambda accum, iter_tuple: _accum_pathways(accum, *iter_tuple)\n",
    "    transposed_df = reduce(accumulate, enrichment_df.iterrows(),\n",
    "                           pd.DataFrame({'Dummy' : []}))\n",
    "    return transposed_df\n",
    "\n",
    "def distinct_pathways(inputs):\n",
    "    \"\"\"\n",
    "    Returns the distinct pathways from the given input data sets.\n",
    "    \n",
    "    :param inputs: the input {name, partition data frame} dictionary\n",
    "    :return: the distinct {name: pathways data series} dictionary\n",
    "    \"\"\"\n",
    "    extract_pathways = lambda table: set(table.index.get_values())\n",
    "    pathway_ds_dict = {name: df.Pathways.apply(extract_pathways)\n",
    "                       for name, df in inputs.items()}\n",
    "    pathways_dict = {name: reduce(lambda s1, s2: s1.union(s2), p)\n",
    "                     for name, p in pathway_ds_dict.items()}\n",
    "    for name, pathways in pathways_dict.items():\n",
    "        print(\"%s has %d pathways.\" % (name, len(pathways)))\n",
    "\n",
    "    intersect = lambda s1, s2: s1.intersection(s2)\n",
    "    common = reduce(intersect, pathways_dict.values())\n",
    "    print(\"There are %d pathways in common between %s and %s.\" %\n",
    "          (len(common), *inputs.keys()))\n",
    "\n",
    "    distinct_dict = {name: sorted(pathways.difference(common))\n",
    "                for name, pathways in pathways_dict.items()}\n",
    "    for name, pathways in distinct_dict.items():\n",
    "        print(\"%s has %d distinct unshared pathways.\" %\n",
    "              (name, len(pathways)))\n",
    "    return distinct_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_hierarchy_node(pathway, hierarchy):\n",
    "    \"\"\"\n",
    "    Returns the Reactome hierarchy node for the given pathway.\n",
    "    \n",
    "    :param pathway: the pathway to check\n",
    "    :param hierarchy: the Reactome pathway hierarchy to check\n",
    "    :return: the hierarchy node\n",
    "    \"\"\"\n",
    "    if hierarchy['name'] == pathway:\n",
    "        return hierarchy\n",
    "    children = hierarchy.get('children')\n",
    "    if children:\n",
    "        for subtree in children:\n",
    "            target = get_hierarchy_node(pathway, subtree)\n",
    "            if target:\n",
    "                return target\n",
    "\n",
    "def export_diagram(db_id, pathway, genes, out_dir=None):\n",
    "    \"\"\"\n",
    "    Exports a diagram PDF for the given pathway. The PDF\n",
    "    is placed in the target output directory. The file name\n",
    "    capitalizes the pathway name and removes spaces and\n",
    "    punctuation, e.g. `DEKBindsTFAP2Homodimers`.\n",
    "    \n",
    "    :param db_id: the Reactome pathway db id\n",
    "    :param pathway: the pathway name to export\n",
    "    :param genes: the hit genes for the pathway\n",
    "    :option out_dir: the target directory (default current directory)\n",
    "    :return: the exported PDF file name\n",
    "    \"\"\"\n",
    "    # Re-enrich the genes in order to get the proper diagram\n",
    "    # highlighting.\n",
    "    enrich_genes(genes)\n",
    "    if not out_dir:\n",
    "        out_dir = os.getcwd()\n",
    "    separator = re.compile('[^\\w]+')\n",
    "    capitalized = [word[0].upper() + word[1:]\n",
    "                   for word in separator.split(pathway) if word]\n",
    "    base_name = ''.join(capitalized)\n",
    "    file_name = os.path.join(out_dir, \"%s.pdf\" % base_name)\n",
    "    body = dict(dbId=db_id, pathwayName=pathway, fileName=file_name)\n",
    "    requests.post(get_fi_url('exportPathwayDiagram'), json=body)\n",
    "    print(\"Exported pathway '%s' to %s.\" % (pathway, file_name))\n",
    "    return file_name\n",
    "\n",
    "def export_first_pathway_diagram(pathways, df, hierarchy):\n",
    "    \"\"\"\n",
    "    Exports the diagram for the first pathway which has one.\n",
    "    \n",
    "    :param pathways: the pathways to choose from\n",
    "    :param df: the enrichment data frame\n",
    "    :param hierarchy: the Reactome pathway hierarchy\n",
    "    :return: the exported diagram file name\n",
    "    \"\"\"\n",
    "    for pathway in pathways:\n",
    "        node = get_hierarchy_node(pathway, hierarchy)\n",
    "        if node and node['hasDiagram']:\n",
    "            db_id = node['dbId']\n",
    "            genes = next(\n",
    "                row.Genes for _, row in df.iterrows()\n",
    "                if pathway in row.Pathways.index\n",
    "            )\n",
    "            return export_diagram(db_id, pathway, genes, sandbox)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "\n",
    "**Prepare the network.**\n",
    "\n",
    "* Download the cohort mutation annotation file (MAF) from the Broad Firehose WS server.\n",
    "\n",
    "* Build the network using the ReactomeFIViz CyREST API buildFISubNetwork resource.\n",
    "\n",
    "* Cluster the network into network modules using the cluster resource.\n",
    "\n",
    "(cf. User Guide\n",
    "[Mutation Analysis](http://wiki.reactome.org/index.php?title=ReactomeFIViz&redirect=no#Gene_Set.2FMutation_Analysis)\n",
    "step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The sandbox directory contains the MAF files and\n",
    "# generated diagrams. The default is the working\n",
    "# directory. Change this to your workflows git\n",
    "# repository location if you don't start the\n",
    "# notebook in that directory.\n",
    "sandbox = os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare the BRCA cohort for analysis.\n",
    "# This takes a while.\n",
    "brca_df = prepare('BRCA')\n",
    "print(\"BRCA module count: %d\" % brca_df.index.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Expected Result*: Cytoscape displays the BRCA FI network\n",
    "and clustered gene modules table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare the OV cohort for analysis.\n",
    "# This takes a while.\n",
    "ov_df = prepare('OV')\n",
    "print(\"OV module count: %d\" % ov_df.index.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Expected Result*: Cytoscape displays the OV FI network\n",
    "and clustered gene modules table."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "\n",
    "**Perform pair-wise module overlapping analysis.**\n",
    "\n",
    "Determine the significance of pair-wise module overlap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# The analysis input.\n",
    "inputs = {'BRCA': brca_df, 'OV': ov_df}\n",
    "\n",
    "# Discover cluster modules with significant overlap.\n",
    "# The default is to select no more than 20 modules \n",
    "# per cancer type and with at least 3 genes per module. \n",
    "# Reset the options below to change these limits.\n",
    "# Note that we already limited the module selection\n",
    "# by setting a sample size threshold when we built\n",
    "# the network.\n",
    "module_opts = dict(max_module_count=None, min_module_size=None)\n",
    "overlap_df = discover_overlap(BACKGROUND_GENE_CNT, inputs,\n",
    "                              **module_opts)\n",
    "\n",
    "# Make a data frame suitable for printing.\n",
    "# Reset the FDR *cutoff* to limit the number of modules\n",
    "# printed, or set to None to print all modules.\n",
    "print_overlap(overlap_df, cutoff=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "\n",
    "**Check overlap for a specific module in a cancer type.**\n",
    "\n",
    "From the above, Two BRCA modules, 13 and 17, are not listed with FDR < 0.05. Check whether or not these two modules have any overlap with OV modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Check BRCA Module 13 overlap with OV modules.\n",
    "# By preliminary analysis, if shared size = 0, FDR should be 1.0. Therefore we\n",
    "# choose FDR cutoff = 0.99 to get a simple output\n",
    "print_overlap(overlap_df, cancer = 'BRCA', module = 13, cutoff = 0.99, show_genes = True)\n",
    "\n",
    "# Check BRCA Module 17 overlap with OV modules.\n",
    "print_overlap(overlap_df, cancer = 'BRCA', module = 17, cutoff = 0.99, show_genes = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "\n",
    "**Partition the modules.**\n",
    "\n",
    "Partition the network modules into the following sets:\n",
    "\n",
    "* Modules with a significant proportion of genes in common\n",
    "\n",
    "* Unshared modules for each of the two cancer types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Partition into shared and unshared.\n",
    "partition = partition_shared(overlap_df, inputs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "\n",
    "**Load the Reactome pathways.**\n",
    "\n",
    "(cf. the User Guide\n",
    "[Explore Pathways](http://wiki.reactome.org/index.php?title=ReactomeFIViz&redirect=no#Explore_Reactome_Pathways)\n",
    "step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resp = requests.get(get_fi_url('pathwayTree'))\n",
    "reactome_tree = resp.json()['data']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Expected Result*: Reactome heirarchy displayed in the Cytoscape Control Panel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "\n",
    "**Perform pathway enrichment on the unshared gene modules.**\n",
    "\n",
    "(cf. User Guide\n",
    "[Enrichment Analysis](http://wiki.reactome.org/index.php?title=ReactomeFIViz&redirect=no#Pathway_Enrichment_Analysis)\n",
    "step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Perform pathway enrichment on the gene modules.\n",
    "rename_opts = {'Shared Genes': 'Genes'}\n",
    "columns = ['Genes', 'p-value', 'FDR']\n",
    "shared_df = partition['shared'].rename(columns=rename_opts).loc[:, columns]\n",
    "groups = dict(Shared=shared_df)\n",
    "groups.update(partition['unshared'])\n",
    "enriched = {name: enrich_modules(name, df.Genes)\n",
    "            for name, df in groups.items()}\n",
    "enriched_unshared = {name: enriched[name] for name in inputs.keys()}\n",
    "distinct = distinct_pathways(enriched_unshared)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the unshared pathways. \n",
    "display(HTML('<h4>Unshared Module Pathways:</h4>'))\n",
    "for name, pathways in distinct.items():\n",
    "    display(HTML(\"<h5>%s:</h5>\" % name))\n",
    "    transposed_df = transpose_pathways(enriched[name])\n",
    "    grouped_df = transposed_df.groupby(level='Pathway', sort=False).min()\n",
    "    grouped_df.sort_values(by='FDR', inplace=True)\n",
    "    criterion = grouped_df.index.isin(pathways)\n",
    "    distinct_df = grouped_df[criterion]\n",
    "    display(HTML(distinct_df.reset_index().to_html(index=False)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Expected Result*: Enrichment results are successively displayed in the Cytoscape Table Panel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The cancer type.\n",
    "cancer = 'BRCA'\n",
    "\n",
    "# The enriched BRCA module numbers.\n",
    "module_numbers = enriched[cancer].index.get_values()\n",
    "display(HTML('<h4>Enriched %s module numbers:</h4>'% cancer))\n",
    "print(module_numbers)\n",
    "\n",
    "# Show the enriched pathways for a BRCA module.\n",
    "# Change the module number to see a different module.\n",
    "module_number = module_numbers[0]\n",
    "df = enriched[cancer].loc[module_number]\n",
    "display(HTML(\"<h4>Module %d Genes:</h4>\" % module_number))\n",
    "print(','.join(df.Genes))\n",
    "print()\n",
    "display(HTML(\"<h4>Module %d Pathways:</h4>\" % module_number))\n",
    "display(HTML(df.Pathways.loc[:, ['p-value', 'FDR']].reset_index().to_html(index=False)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___\n",
    "\n",
    "**Export an unshared pathway diagram for each cancer type.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export one diagram per cancer type.\n",
    "for name, pathways in distinct.items():\n",
    "    print(\"Exporting a %s pathway...\" % name)\n",
    "    df = enriched[name]\n",
    "    export_first_pathway_diagram(pathways, df, reactome_tree)"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}