DSP_genomics.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "DSP genomics.ipynb",
      "provenance": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/matthewberry/uiuc_com_dsp/blob/master/DSP_genomics.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8ZlQjjAmJuSY",
        "colab_type": "text"
      },
      "source": [
        "## Using This Notebook\n",
        "\n",
        "This notebook is an interactive environment that combines explanatory text with executable code. It provides computational tools useful in genomics, and you will familiarize yourself with them by stepping through a series of analyses that examine related data from multiple angles. This suite of tools and example can then serve as a starting point for a project of your own design.\n",
        "\n",
        "If you are new to notebooks, you might find this introduction helpful: [Overview of Colaboratory Features](https://colab.research.google.com/notebooks/basic_features_overview.ipynb) You might also want to refer to the [Python 3 documentation](https://docs.python.org/3/).\n",
        "\n",
        "**Important Note**: After a period of inactivity (Google does not specify exactly how long), Google will disconnect your notebook from the virtual machine that had been running it. When you return, Google will connect to a new virtual machine. Any data files you saved to your Google Drive will remain, but any variables or methods defined in your previous virtual machine will have to be reloaded. (You will know when this happens, because cells that previously ran without error will suddenly stop working, and the notebook will lose its connection to your Google Drive.) To reload the variable and method definitions, and to restore the connection to your Google Drive, you can simply re-run the cells that perform those tasks. This notebook will explain which cells might need to be re-run.\n",
        "\n",
        "### Before You Proceed\n",
        "\n",
        "In order to run the notebook, you will need your own copy of it, along with your own copy of the data. Here are the steps to follow:\n",
        "\n",
        "1. If you have not already enabled Google Apps @ Illinois, which allows you to use Google Drive, Google Docs, and so on with your illinois.edu account, [enable Google Apps @ Illinois](https://answers.uillinois.edu/illinois/page.php?id=55049).\n",
        "2. Check the currently active Google account on this notebook. Look near the top-right corner of the screen for either a `Sign In` button or a round icon containing either a letter or your profile photo. If you see a `Sign In` button, click it and follow the prompts to sign in with your illinois.edu account. Otherwise, click the icon to open a popup that identifies the currenly active Google account. If the account is not your illinois.edu account, switch to your illinois.edu account (you might have to click `Add Account` if it is not already an option in the list).\n",
        "3. In the `File` menu above, select `Save a copy in Drive...`. This step will create a new browser tab containing your copy of the notebook. At this point, you can close the old browser tab that contained the original copy of the notebook.\n",
        "4. Open the `File` menu and select `Locate in Drive`, which will show you where you can find the notebook if you need to open it later.\n",
        "5. [Click here](https://drive.google.com/drive/folders/1vsIBpVR0xi56u0WtihppTXrWscDBO70U?usp=sharing) to open the master copy of the example data in a new tab. In the new tab, again make sure the active Google account is your illinois.edu account. Click on the small triangle that appears after the folder name near the top of the screen. From the menu that appears, select `Add to My Drive`.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RczciUOYuOzp",
        "colab_type": "text"
      },
      "source": [
        "## Installation\n",
        "\n",
        "The cell below installs software required to perform the analyses. Run the cell and wait for it to complete, which will take about two minutes. You'll see lots of text output as the cell runs, but there's no need to read it unless the following cell fails.\n",
        "\n",
        "You **will** need to re-run this cell whenever your notebook connects to a new virtual machine."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HgrIxPJT-UiT",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!pip3 install -I pyyaml==5.1.2\n",
        "!pip3 install xmlrunner==1.7.7 redis==3.3.8 lifelines==0.22.8\n",
        "!pip3 install git+https://github.com/KnowEnG/KnowEnG_Pipelines_Library.git@mjberry/update_dependencies\n",
        "!pip3 install git+https://github.com/KnowEnG/Data_Cleanup_Pipeline.git@mjberry/create_package\n",
        "!pip3 install git+https://github.com/KnowEnG/General_Clustering_Pipeline.git@mjberry/create_package\n",
        "!pip3 install git+https://github.com/KnowEnG/Samples_Clustering_Pipeline.git@mjberry/create_package\n",
        "!pip3 install git+https://github.com/KnowEnG/Feature_Prioritization_Pipeline.git@mjberry/create_package\n",
        "!pip3 install git+https://github.com/KnowEnG/Gene_Prioritization_Pipeline.git@mjberry/create_package\n",
        "!pip3 install git+https://github.com/KnowEnG/Geneset_Characterization_Pipeline.git@mjberry/create_package\n",
        "!pip3 install git+https://github.com/KnowEnG/Spreadsheets_Transformation.git@mjberry/create_package"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ftgYqoTDuzOu",
        "colab_type": "text"
      },
      "source": [
        "## Environment Setup\n",
        "\n",
        "The cell below sets up the environment for running the analyses. Run the cell and wait for it to complete, which will only take a few seconds. You won't see any text output this time.\n",
        "\n",
        "You probably will not need to call any of the methods defined in this method, and you probably will not need to edit anything in the cell.\n",
        "\n",
        "You **will** need to re-run the cell whenever your notebook connects to a new virtual machine.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_ifxh2nzm173",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import csv\n",
        "import os\n",
        "import shutil\n",
        "from tempfile import mkdtemp\n",
        "import urllib.request\n",
        "\n",
        "from IPython.display import HTML\n",
        "from lifelines import KaplanMeierFitter\n",
        "from lifelines.statistics import multivariate_logrank_test\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd\n",
        "\n",
        "from kndatacleanup import data_cleanup\n",
        "from knfeatureprioritization import feature_prioritization\n",
        "from kngeneprioritization import gene_prioritization\n",
        "from kngenesetcharacterization import geneset_characterization\n",
        "from knsamplesclustering import samples_clustering\n",
        "from kngeneralclustering import general_clustering\n",
        "from knspreadsheetstransformation.spreadsheets_transformation_toolbox import \\\n",
        "    get_cluster_binary_dataframe\n",
        "\n",
        "%matplotlib inline\n",
        "\n",
        "NETWORK_DIR_PATH = '/network/'\n",
        "\n",
        "REDIS_PARAMS = {\n",
        "    'host': 'knowredis.knoweng.org',\n",
        "    'password': 'KnowEnG',\n",
        "    'port': 6379\n",
        "}\n",
        "\n",
        "NUM_CPUS = 2\n",
        "\n",
        "def fetch_network(edge_file_path):\n",
        "    \"\"\"Given the local path to an edge file, ensures that the edge file exists\n",
        "    on disk, downloading it from AWS if necessary.\n",
        "\n",
        "    Arguments:\n",
        "        edge_file_path (str): The local path to an edge file as found in the\n",
        "            data returned by `get_interaction_networks` and\n",
        "            `get_gene_property_networks`.\n",
        "\n",
        "    Returns:\n",
        "        None\n",
        "\n",
        "    \"\"\"\n",
        "    if not os.path.isfile(edge_file_path):\n",
        "        url = \"https://s3.amazonaws.com/KnowNets/KN-20rep-1706/\" + \\\n",
        "            \"userKN-20rep-1706/\" + edge_file_path[len(NETWORK_DIR_PATH):]\n",
        "        os.makedirs(os.path.dirname(edge_file_path), exist_ok=True)\n",
        "        with urllib.request.urlopen(url) as response:\n",
        "            with open(edge_file_path, 'wb') as out_file:\n",
        "                shutil.copyfileobj(response, out_file)\n",
        "\n",
        "def fetch_network_metadata():\n",
        "    \"\"\"Downloads from AWS the network overview metadata files required to\n",
        "    implement the network utility methods.\n",
        "\n",
        "    Arguments:\n",
        "        None\n",
        "\n",
        "    Returns:\n",
        "        None\n",
        "\n",
        "    \"\"\"\n",
        "    filenames = ['db_contents.txt', 'species_desc.txt', 'edge_type.txt']\n",
        "    for filename in filenames:\n",
        "        out_file_path = os.path.join(NETWORK_DIR_PATH, filename)\n",
        "        if not os.path.isfile(out_file_path):\n",
        "            url = \"https://s3.amazonaws.com/KnowNets/KN-20rep-1706/\" + \\\n",
        "                \"userKN-20rep-1706/\" + filename\n",
        "            with urllib.request.urlopen(url) as response:\n",
        "                with open(out_file_path, 'wb') as out_file:\n",
        "                    shutil.copyfileobj(response, out_file)\n",
        "\n",
        "def get_path_to_newest_file_having_prefix(search_dir_path, prefix):\n",
        "    \"\"\"Finds all files in `search_dir_path` whose name begins with `prefix`.\n",
        "    Returns the newest of these matching files, or returns None if no matching\n",
        "    file exists.\n",
        "\n",
        "    Arguments:\n",
        "        search_dir_path (str): The local path to the directory to search.\n",
        "        prefix (str): The string used to filter the files in `search_dir_path`.\n",
        "\n",
        "    Returns:\n",
        "        str: The path to the newest matching file, or None if no matching files\n",
        "            exist.\n",
        "\n",
        "    \"\"\"\n",
        "    matches = [os.path.join(search_dir_path, name) \\\n",
        "        for name in os.listdir(search_dir_path) \\\n",
        "        if name.startswith(prefix)]\n",
        "    # ensure they're all files\n",
        "    matches = [m for m in matches if os.path.isfile(m)]\n",
        "    return_val = None\n",
        "    if matches:\n",
        "        return_val = sorted(matches, \\\n",
        "            key=lambda path: os.path.getctime(path), reverse=True)[0]\n",
        "    return return_val\n",
        "\n",
        "def get_cleaned_file_path(original_file_path, results_dir_path):\n",
        "    \"\"\"Given the name of a file passed to `kndatacleanup.data_cleanup`,\n",
        "    along with the `results_dir_path` passed to `kndatacleanup.data_cleanup`,\n",
        "    returns the path at which the cleaned version of the input can be found.\n",
        "\n",
        "    Arguments:\n",
        "        original_file_path (str): The path to the file that was passed to\n",
        "            `kndatacleanup.data_cleanup`.\n",
        "        results_dir_path (str): The path to the results directory that was\n",
        "            passed to `kndatacleanup.data_cleanup`.\n",
        "\n",
        "    Returns:\n",
        "        str: The path to the cleaned version of `original_file_path` that was\n",
        "            or would be produced by `kndatacleanup.data_cleanup`.\n",
        "\n",
        "    \"\"\"\n",
        "    original_name = os.path.basename(original_file_path)\n",
        "    original_name_root = os.path.splitext(original_name)[0]\n",
        "    return os.path.join(results_dir_path, original_name_root + \"_ETL.tsv\")\n",
        "\n",
        "def get_gene_map_file_path(original_file_path, results_dir_path):\n",
        "    \"\"\"Given the name of an omics file passed to `kndatacleanup.data_cleanup`,\n",
        "    along with the `results_dir_path` passed to `kndatacleanup.data_cleanup`,\n",
        "    returns the path at which the mapping of gene names to gene identifiers\n",
        "    can be found.\n",
        "\n",
        "    Arguments:\n",
        "        original_file_path (str): The path to the omics file that was passed to\n",
        "            `kndatacleanup.data_cleanup`.\n",
        "        results_dir_path (str): The path to the results directory that was\n",
        "            passed to `kndatacleanup.data_cleanup`.\n",
        "\n",
        "    Returns:\n",
        "        str: The path to the gene-name mapping file that was or would be\n",
        "            produced by `kndatacleanup.data_cleanup`.\n",
        "\n",
        "    \"\"\"\n",
        "    original_name = os.path.basename(original_file_path)\n",
        "    original_name_root = os.path.splitext(original_name)[0]\n",
        "    return os.path.join(results_dir_path, original_name_root + \"_MAP.tsv\")\n",
        "\n",
        "os.makedirs(NETWORK_DIR_PATH, exist_ok=True)\n",
        "fetch_network_metadata()\n",
        "\n",
        "!rm -rf /content/sample_data"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "j4DCoKLVwBvT",
        "colab_type": "text"
      },
      "source": [
        "## Knowledge Network Utility Methods\n",
        "\n",
        "The cell below defines several utility methods for working with the knowledge network. These methods are used in the example analyses and might be useful to you in your project. Run the cell and wait for it to complete, which will only take a second or so. It won't produce any text output.\n",
        "\n",
        "You probably will not need to edit anything within the cell.\n",
        "\n",
        "A later cell shows how to use the knowledge network utility methods.\n",
        "\n",
        "You **will** need to re-run the cell whenever your notebook connects to a new virtual machine.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JKn3Uxx6rLce",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def get_network_species():\n",
        "    \"\"\"Returns information about the species found in the knowledge network.\n",
        "\n",
        "    Arguments:\n",
        "        None\n",
        "\n",
        "    Returns:\n",
        "        list: A list in which each element is a dictionary. Each dictionary has\n",
        "            keys 'id' (which can be passed to other methods that require a\n",
        "            `species_id`), 'short_latin_name', 'latin_name', 'familiar_name',\n",
        "            and 'group_name'.\n",
        "\n",
        "    \"\"\"\n",
        "    return_val = []\n",
        "    species_file_path = os.path.join(NETWORK_DIR_PATH, 'species_desc.txt')\n",
        "    with open(species_file_path) as csvfile:\n",
        "        for row in csv.reader(csvfile, delimiter='\\t'):\n",
        "            return_val.append({\n",
        "                'id': row[0],\n",
        "                'short_latin_name': row[1],\n",
        "                'latin_name': row[2],\n",
        "                'familiar_name': row[3],\n",
        "                'group_name': row[5]\n",
        "            })\n",
        "    return return_val\n",
        "\n",
        "def display_network_species():\n",
        "    \"\"\"Displays a table of the species found in the knowledge network.\n",
        "\n",
        "    Arguments:\n",
        "        None\n",
        "\n",
        "    Returns:\n",
        "        None\n",
        "\n",
        "    \"\"\"\n",
        "    html_string = \"<table><tr><th>Familiar Name (Latin Name)</th><th>Species Id</th></tr>\"\n",
        "    for species in get_network_species():\n",
        "        html_string += \"<tr><td>\" + species['familiar_name'] + \" (\" + \\\n",
        "            species['latin_name'] + \")</td><td>\" + species['id'] + \"</td></tr>\"\n",
        "    html_string += \"</table>\"\n",
        "    return HTML(html_string)\n",
        "\n",
        "def get_edge_type_name_to_pretty_name():\n",
        "    \"\"\"Returns a dictionary in which the keys are edge type names and the values\n",
        "    are pretty network names.\n",
        "\n",
        "    Arguments:\n",
        "        None\n",
        "\n",
        "    Returns:\n",
        "        dict: A dictionary in which the keys are edge type names and the values\n",
        "            are pretty network names.\n",
        "\n",
        "    \"\"\"\n",
        "    return_val = {}\n",
        "    file_path = os.path.join(NETWORK_DIR_PATH, 'edge_type.txt')\n",
        "    with open(file_path) as csvfile:\n",
        "        for row in csv.DictReader(csvfile, delimiter='\\t'):\n",
        "            return_val[row['et_name']] = row['pretty_name']\n",
        "    return return_val\n",
        "\n",
        "def get_interaction_networks(species_id):\n",
        "    \"\"\"Given a `species_id`, returns information about the interaction networks\n",
        "    available in the knowledge network.\n",
        "\n",
        "    Arguments:\n",
        "        species_id (int or str): The id for the species of interest, as returned\n",
        "            by `get_network_species` or displayed by `display_network_species`.\n",
        "\n",
        "    Returns:\n",
        "        list: A list in which each element is a dictionary. Each dictionary has\n",
        "            two keys, 'name' and 'edge_file_path'.\n",
        "\n",
        "    \"\"\"\n",
        "    species_id = str(species_id) # user-friendliness\n",
        "    return_val = []\n",
        "    contents_file_path = os.path.join(NETWORK_DIR_PATH, 'db_contents.txt')\n",
        "    with open(contents_file_path) as csvfile:\n",
        "        edge_type_name_to_pretty_name = get_edge_type_name_to_pretty_name()\n",
        "        for row in csv.DictReader(csvfile, delimiter='\\t'):\n",
        "            if row['n1_type'] == 'Gene' and row['taxon'] == species_id:\n",
        "                return_val.append({\n",
        "                    'name': edge_type_name_to_pretty_name[row['et_name']],\n",
        "                    'edge_file_path': os.path.join(\\\n",
        "                        NETWORK_DIR_PATH, 'Gene', species_id, row['et_name'], \\\n",
        "                        species_id + '.' + row['et_name'] + '.edge')\n",
        "                })\n",
        "    return return_val\n",
        "\n",
        "def display_interaction_networks(species_id):\n",
        "    \"\"\"Given a `species_id`, displays information about the interaction\n",
        "    networks available in the knowledge network.\n",
        "\n",
        "    Arguments:\n",
        "        species_id (int or str): The id for the species of interest, as returned\n",
        "            by `get_network_species` or displayed by `display_network_species`.\n",
        "\n",
        "    Returns:\n",
        "        None\n",
        "\n",
        "    \"\"\"\n",
        "    html_string = \"<table><tr><th>Interaction Network Name</th><th>Edge File Path</th></tr>\"\n",
        "    for network in get_interaction_networks(species_id):\n",
        "        html_string += \"<tr><td>\" + network['name'] + \"</td><td>\" + \\\n",
        "            network['edge_file_path'] + \"</td></tr>\"\n",
        "    html_string += \"</table>\"\n",
        "    return HTML(html_string)\n",
        "\n",
        "def get_gene_property_networks(species_id):\n",
        "    \"\"\"Given a `species_id`, returns information about the gene-property\n",
        "    networks available in the knowledge network.\n",
        "\n",
        "    Arguments:\n",
        "        species_id (int or str): The id for the species of interest, as returned\n",
        "            by `get_network_species` or displayed by `display_network_species`.\n",
        "\n",
        "    Returns:\n",
        "        list: A list in which each element is a dictionary. Each dictionary has\n",
        "            two keys, 'name' and 'edge_file_path'.\n",
        "\n",
        "    \"\"\"\n",
        "    species_id = str(species_id) # user-friendliness\n",
        "    return_val = []\n",
        "    contents_file_path = os.path.join(NETWORK_DIR_PATH, 'db_contents.txt')\n",
        "    with open(contents_file_path) as csvfile:\n",
        "        edge_type_name_to_pretty_name = get_edge_type_name_to_pretty_name()\n",
        "        for row in csv.DictReader(csvfile, delimiter='\\t'):\n",
        "            if row['n1_type'] == 'Property' and row['taxon'] == species_id:\n",
        "                return_val.append({\n",
        "                    'name': edge_type_name_to_pretty_name[row['et_name']],\n",
        "                    'edge_file_path': os.path.join(\\\n",
        "                        NETWORK_DIR_PATH, 'Property', species_id, row['et_name'], \\\n",
        "                        species_id + '.' + row['et_name'] + '.edge')\n",
        "                })\n",
        "    return return_val\n",
        "\n",
        "def display_gene_property_networks(species_id):\n",
        "    \"\"\"Given a `species_id`, displays information about the gene-property\n",
        "    networks available in the knowledge network.\n",
        "\n",
        "    Arguments:\n",
        "        species_id (int or str): The id for the species of interest, as returned\n",
        "            by `get_network_species` or displayed by `display_network_species`.\n",
        "\n",
        "    Returns:\n",
        "        None\n",
        "\n",
        "    \"\"\"\n",
        "    html_string = \"<table><tr><th>Interaction Network Name</th><th>Edge File Path</th></tr>\"\n",
        "    for network in get_gene_property_networks(species_id):\n",
        "        html_string += \"<tr><td>\" + network['name'] + \"</td><td>\" + \\\n",
        "            network['edge_file_path'] + \"</td></tr>\"\n",
        "    html_string += \"</table>\"\n",
        "    return HTML(html_string)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aVvYGbnRxE-9",
        "colab_type": "text"
      },
      "source": [
        "### Using the Knowledge Network Utility Methods\n",
        "\n",
        "The three cells below show how `display_network_species`, `display_interaction_networks`, and `display_gene_property_networks` can be called to view information about the knowledge network. This information can be useful in configuring analyses, as you'll see later.\n",
        "\n",
        "These methods are based on three other methods defined in the cell above, `get_network_species`, `get_interaction_networks`, and `get_gene_property_networks`. The \"get\" versions return the same information as the \"display\" versions, but the \"get\" versions return it in a format convenient for use in code instead of a format that's easy to read.\n",
        "\n",
        "You **will not** need to re-run these three cells whenever your notebook connects to a new virtual machine."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dSMHjVi6tnNc",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# display all species in the knowledge network\n",
        "display_network_species()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uep6yYfe83ly",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# display interaction networks for rat (species id 10116)\n",
        "display_interaction_networks('10116')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cNl8XB-T84x5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# display gene property networks for roundworm (species id 6239)\n",
        "display_gene_property_networks('6239')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BmTUV9Lp_rK0",
        "colab_type": "text"
      },
      "source": [
        "## Analytics Methods\n",
        "\n",
        "The cell below defines methods for running clustering, prioritization, and gene-set characterization. Run the cell and wait for it to complete, which will take a second or so. It won't produce any output.\n",
        "\n",
        "You **will** need to re-run the cell whenever your notebook connects to a new virtual machine.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0M6MGR3K_pd3",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "def do_clustering(\\\n",
        "    omics_file_path, phenotype_file_path, results_dir_path, num_clusters, \\\n",
        "    species_id, interaction_network_edge_file_path, network_influence, \\\n",
        "    num_bootstraps, bootstrap_sample_fraction):\n",
        "    \"\"\"Performs a clustering upon the samples found in `omics_file_path`.\n",
        "\n",
        "    Arguments:\n",
        "        omics_file_path (str): The path to the omics file.\n",
        "        phenotype_file_path (str): The path to a file containing phenotype data\n",
        "            on the same samples as found in `omics_file_path`, or None if no\n",
        "            phenotype data are to be analyzed. If analyzed, each phenotype will\n",
        "            scored for statistically significant differences between the\n",
        "            clusters.\n",
        "        results_dir_path (str): The path to a directory where results files\n",
        "            should be stored.\n",
        "        num_clusters (int): The number of clusters to create.\n",
        "        species_id (int or str): The id for the species of interest, as returned\n",
        "            by `get_network_species` or displayed by `display_network_species`,\n",
        "            or None not using an `interaction_network_edge_file_path`.\n",
        "        interaction_network_edge_file_path (str): The path to an interaction\n",
        "            network edge file, to use a knowledge-guided approach to clustering,\n",
        "            or else None.\n",
        "        network_influence (float): A number between 0 and 1 that specifies the\n",
        "            amount to which network data should influence the results, or None\n",
        "            if not using an `interaction_network_edge_file_path`.\n",
        "        num_bootstraps (int): A number of bootstrap iterations to run. Use 0 for\n",
        "            no bootstrapping.\n",
        "        bootstrap_sample_fraction (float): A number between 0 and 1 that\n",
        "            specifies what fraction of the data should be used in each bootstrap\n",
        "            iteration, or None if not using bootstrapping.\n",
        "\n",
        "    Returns:\n",
        "        None\n",
        "\n",
        "    \"\"\"\n",
        "    try:\n",
        "        species_id = str(species_id) # user-friendliness\n",
        "        os.makedirs(results_dir_path, exist_ok=True)\n",
        "\n",
        "        if interaction_network_edge_file_path is None:\n",
        "            pipeline_type = 'general_clustering_pipeline'\n",
        "        else:\n",
        "            fetch_network(interaction_network_edge_file_path)\n",
        "            pipeline_type = 'samples_clustering_pipeline'\n",
        "\n",
        "        cleanup_parameters = {\n",
        "            'spreadsheet_name_full_path': omics_file_path,\n",
        "            'pipeline_type': pipeline_type,\n",
        "        'results_directory': results_dir_path\n",
        "        }\n",
        "        if phenotype_file_path is not None:\n",
        "            cleanup_parameters['phenotype_name_full_path'] = phenotype_file_path\n",
        "        if interaction_network_edge_file_path is not None:\n",
        "            cleanup_parameters.update({\n",
        "                'gg_network_name_full_path': interaction_network_edge_file_path,\n",
        "                'taxonid': species_id,\n",
        "                'source_hint': '',\n",
        "                'redis_credential': {\n",
        "                    'host': REDIS_PARAMS['host'],\n",
        "                    'port': REDIS_PARAMS['port'],\n",
        "                    'password': REDIS_PARAMS['password']\n",
        "                }\n",
        "            })\n",
        "        data_cleanup.run_pipelines(cleanup_parameters, \\\n",
        "            data_cleanup.SELECT[pipeline_type])\n",
        "\n",
        "        clustering_parameters = {\n",
        "            'spreadsheet_name_full_path': get_cleaned_file_path(\\\n",
        "                omics_file_path, results_dir_path),\n",
        "            'results_directory': results_dir_path,\n",
        "            'processing_method': 'parallel',\n",
        "            'parallelism': NUM_CPUS,\n",
        "            'number_of_clusters': num_clusters,\n",
        "            'run_directory': results_dir_path,\n",
        "            'tmp_directory': './tmp'\n",
        "        }\n",
        "        if phenotype_file_path is not None:\n",
        "            clustering_parameters.update({\n",
        "                'phenotype_name_full_path': get_cleaned_file_path(\\\n",
        "                    phenotype_file_path, results_dir_path),\n",
        "                'threshold': 15\n",
        "            })\n",
        "\n",
        "        method_prefix = ''\n",
        "        if num_bootstraps > 0:\n",
        "            clustering_parameters.update({\n",
        "                'number_of_bootstraps': num_bootstraps,\n",
        "                'rows_sampling_fraction': 1.0,\n",
        "                'cols_sampling_fraction': bootstrap_sample_fraction\n",
        "            })\n",
        "            method_prefix = 'cc_'\n",
        "\n",
        "        if interaction_network_edge_file_path is not None:\n",
        "            clustering_parameters.update({\n",
        "                'gg_network_name_full_path': interaction_network_edge_file_path,\n",
        "                'rwr_max_iterations': 100,\n",
        "                'rwr_convergence_tolerence': 1.0e-4,\n",
        "                'rwr_restart_probability': network_influence,\n",
        "                'top_number_of_genes': 100,\n",
        "                'nmf_conv_check_freq': 50,\n",
        "                'nmf_max_invariance': 200,\n",
        "                'nmf_max_iterations': 10000,\n",
        "                'nmf_penalty_parameter': 1400,\n",
        "                'method': method_prefix + 'net_nmf'\n",
        "            })\n",
        "            samples_clustering.SELECT[clustering_parameters['method']](\\\n",
        "                clustering_parameters)\n",
        "        else:\n",
        "            clustering_parameters.update({\n",
        "                'top_number_of_rows': 100,\n",
        "                'affinity_metric': 'euclidean',\n",
        "                'linkage_criterion': 'ward',\n",
        "                'method': method_prefix + 'hclust'\n",
        "            })\n",
        "            general_clustering.SELECT[clustering_parameters['method']](\\\n",
        "                clustering_parameters)\n",
        "    except:\n",
        "        print(\"Something went wrong! Check the debugging information below, \" + \\\n",
        "            \"and look for log output in \" + results_dir_path)\n",
        "        raise\n",
        "    else:\n",
        "        print(\"Find results in \" + results_dir_path)\n",
        "\n",
        "def do_prioritization(\\\n",
        "    omics_file_path, phenotype_file_path, results_dir_path, \\\n",
        "    correlation_measure, missing_value_strategy, num_exported_features, \\\n",
        "    num_response_correlated_features, species_id, \\\n",
        "    interaction_network_edge_file_path, network_influence):\n",
        "    \"\"\"Prioritizes the features (genes or otherwise) found in `omics_file_path`\n",
        "    for each phenotype found in `phenotype_file_path`.\n",
        "\n",
        "    Arguments:\n",
        "        omics_file_path (str): The path to the omics file.\n",
        "        phenotype_file_path (str): The path to a file containing phenotype data\n",
        "            on the same samples as found in `omics_file_path`.\n",
        "        results_dir_path (str): The path to a directory where results files\n",
        "            should be stored.\n",
        "        correlation_measure (str): Either 't_test' for binary or categorical\n",
        "            phenotypes or 'pearson' for numeric phenotypes.\n",
        "        missing_value_strategy (str): Governs how to handle missing values in\n",
        "            `omics_file_path`. Options are 'average' to use the average value\n",
        "            for the feature among the other samples, 'remove' to drop any\n",
        "            samples with missing values, or 'reject' to fail if any missing\n",
        "            values are found (perhaps as a sanity check if you believe missing\n",
        "            values were prevented upstream).\n",
        "        num_exported_features (int): The number of top features per phenotype to\n",
        "            include in the matrix that can be passed to `do_characterization`.\n",
        "        num_response_correlated_features (int): The number of top features to\n",
        "            retain from the first stage of the analysis if using an\n",
        "            `interaction_network_edge_file_path`.\n",
        "        species_id (int or str): The id for the species of interest, as returned\n",
        "            by `get_network_species` or displayed by `display_network_species`,\n",
        "            or none if not using an `interactive_network_edge_file_path`.\n",
        "        interaction_network_edge_file_path (str): The path to an interaction\n",
        "            network edge file, to use a knowledge-guided approach to\n",
        "            prioritization, or else None.\n",
        "        network_influence (float): A number between 0 and 1 that specifies the\n",
        "            amount to which network data should influence the results, or None\n",
        "            if not using an `interaction_network_edge_file_path`.\n",
        "\n",
        "    Returns:\n",
        "        None\n",
        "\n",
        "    \"\"\"\n",
        "    try:\n",
        "        species_id = str(species_id) # user-friendliness\n",
        "        os.makedirs(results_dir_path, exist_ok=True)\n",
        "\n",
        "        if interaction_network_edge_file_path is None:\n",
        "            pipeline_type = 'feature_prioritization_pipeline'\n",
        "        else:\n",
        "            fetch_network(interaction_network_edge_file_path)\n",
        "            pipeline_type = 'gene_prioritization_pipeline'\n",
        "\n",
        "        cleanup_parameters = {\n",
        "            'spreadsheet_name_full_path': omics_file_path,\n",
        "            'phenotype_name_full_path': phenotype_file_path,\n",
        "            'pipeline_type': pipeline_type,\n",
        "            'correlation_measure': correlation_measure, # t_test, pearson, edgeR\n",
        "            'impute': missing_value_strategy, # average, remove, reject\n",
        "            'results_directory': results_dir_path\n",
        "        }\n",
        "        if interaction_network_edge_file_path is not None:\n",
        "            cleanup_parameters.update({\n",
        "                'taxonid': species_id,\n",
        "                'source_hint': '',\n",
        "                'redis_credential': {\n",
        "                    'host': REDIS_PARAMS['host'],\n",
        "                    'port': REDIS_PARAMS['port'],\n",
        "                    'password': REDIS_PARAMS['password']\n",
        "                }\n",
        "            })\n",
        "        data_cleanup.run_pipelines(cleanup_parameters, \\\n",
        "            data_cleanup.SELECT[pipeline_type])\n",
        "\n",
        "        prioritization_parameters = {\n",
        "            'correlation_measure': correlation_measure,\n",
        "            'spreadsheet_name_full_path': get_cleaned_file_path(\\\n",
        "                omics_file_path, results_dir_path),\n",
        "            'phenotype_name_full_path': get_cleaned_file_path(\\\n",
        "                phenotype_file_path, results_dir_path),\n",
        "            'results_directory': results_dir_path,\n",
        "            'top_gamma_of_sort': num_exported_features,\n",
        "            'max_cpu': NUM_CPUS\n",
        "        }\n",
        "        if interaction_network_edge_file_path is not None:\n",
        "            prioritization_parameters.update({\n",
        "                'gg_network_name_full_path': interaction_network_edge_file_path,\n",
        "                'rwr_max_iterations': 100,\n",
        "                'rwr_convergence_tolerence': 1.0e-4,\n",
        "                'rwr_restart_probability': network_influence,\n",
        "                'top_beta_of_sort': num_response_correlated_features,\n",
        "                'method': 'net_correlation'\n",
        "            })\n",
        "            gene_prioritization.net_correlation(prioritization_parameters)\n",
        "        else:\n",
        "            prioritization_parameters.update({\n",
        "                'top_beta_of_sort': num_exported_features,\n",
        "                'method': 'correlation',\n",
        "            })\n",
        "            feature_prioritization.correlation(prioritization_parameters)\n",
        "    except:\n",
        "        print(\"Something went wrong! Check the debugging information below, \" + \\\n",
        "            \"and look for log output in \" + results_dir_path)\n",
        "        raise\n",
        "    else:\n",
        "        print(\"Find results in \" + results_dir_path)\n",
        "\n",
        "def do_characterization(\\\n",
        "    gene_matrix_file_path, results_dir_path, species_id, \\\n",
        "    gene_property_edge_file_path, interaction_network_edge_file_path, \\\n",
        "    network_influence):\n",
        "    \"\"\"Compares user-submitted gene sets to those found in a gene-property\n",
        "    network from the knowledge network.\n",
        "\n",
        "    Arguments:\n",
        "        gene_matrix_file_path (str): The path to the gene matrix that defines\n",
        "            one or more gene sets.\n",
        "        results_dir_path (str): The path to a directory where results files\n",
        "            should be stored.\n",
        "        species_id (int or str): The id for the species of interest, as returned\n",
        "            by `get_network_species` or displayed by `display_network_species`.\n",
        "        gene_property_edge_file_path: The path to a gene-property network edge\n",
        "            file.\n",
        "        interaction_network_edge_file_path (str): The path to an interaction\n",
        "            network edge file, to use a knowledge-guided approach to\n",
        "            characterization, or else None.\n",
        "        network_influence (float): A number between 0 and 1 that specifies the\n",
        "            amount to which network data should influence the results, or None\n",
        "            if not using an `interaction_network_edge_file_path`.\n",
        "\n",
        "    Returns:\n",
        "        None\n",
        "\n",
        "    \"\"\"\n",
        "    try:\n",
        "        species_id = str(species_id) # user-friendliness\n",
        "        os.makedirs(results_dir_path, exist_ok=True)\n",
        "\n",
        "        fetch_network(gene_property_edge_file_path)\n",
        "\n",
        "        cleanup_parameters = {\n",
        "            'spreadsheet_name_full_path': gene_matrix_file_path,\n",
        "            'pipeline_type': 'geneset_characterization_pipeline',\n",
        "            'results_directory': results_dir_path,\n",
        "            'taxonid': species_id,\n",
        "            'source_hint': '',\n",
        "            'redis_credential': {\n",
        "                'host': REDIS_PARAMS['host'],\n",
        "                'port': REDIS_PARAMS['port'],\n",
        "                'password': REDIS_PARAMS['password']\n",
        "            }\n",
        "        }\n",
        "        data_cleanup.run_pipelines(cleanup_parameters, \\\n",
        "            data_cleanup.SELECT['geneset_characterization_pipeline'])\n",
        "\n",
        "        characterization_parameters = {\n",
        "            'spreadsheet_name_full_path': get_cleaned_file_path(\\\n",
        "                gene_matrix_file_path, results_dir_path),\n",
        "            'gene_names_map': get_gene_map_file_path(\\\n",
        "                gene_matrix_file_path, results_dir_path),\n",
        "            'results_directory': results_dir_path,\n",
        "            'pg_network_name_full_path': gene_property_edge_file_path,\n",
        "            'max_cpu': NUM_CPUS\n",
        "        }\n",
        "        if interaction_network_edge_file_path is None:\n",
        "            characterization_parameters.update({\n",
        "                'method': 'fisher'\n",
        "            })\n",
        "            geneset_characterization.fisher(characterization_parameters)\n",
        "        else:\n",
        "            fetch_network(interaction_network_edge_file_path)\n",
        "            characterization_parameters.update({\n",
        "                'method': 'DRaWR',\n",
        "                'rwr_max_iterations': 500,\n",
        "                'rwr_convergence_tolerence': 1.0e-4,\n",
        "                'rwr_restart_probability': network_influence,\n",
        "                'gg_network_name_full_path': interaction_network_edge_file_path\n",
        "            })\n",
        "            geneset_characterization.DRaWR(characterization_parameters)\n",
        "    except:\n",
        "        print(\"Something went wrong! Check the debugging information below, \" + \\\n",
        "            \"and look for log output in \" + results_dir_path)\n",
        "        raise\n",
        "    else:\n",
        "        print(\"Find results in \" + results_dir_path)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qRxFpjpTM0b8",
        "colab_type": "text"
      },
      "source": [
        "## Connect to Google Drive\n",
        "\n",
        "The cell below enables this notebook to use your Google Drive for file storage. Subsequent cells will use this access to load the example files you copied earlier and to save results of the example analyses. You might also find this helpful in running your own analyses.\n",
        "\n",
        "Run the cell and click on the link that appears in the output. On the linked page, select your illinois.edu account and grant the requested permissions. The page will then display a code. Copy the code and paste it in the box that appears in the output below. Then press Enter.\n",
        "\n",
        "You **will** need to re-run the cell whenever your notebook connects to a new virtual machine.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gTqz9ZNXvlG_",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from google.colab import drive\n",
        "GDRIVE_MOUNT_PATH = '/content/gdrive'\n",
        "drive.mount(GDRIVE_MOUNT_PATH)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EC9dFiPKnLmv",
        "colab_type": "text"
      },
      "source": [
        "## Setting File Locations\n",
        "\n",
        "In the cell below, we will tell the notebook where the example files can be found and where the results should be saved.\n",
        "\n",
        "To confirm the locations, find the file-folder icon in the narrow bar that runs along the left side of this screen. Click it to reveal the `Files` panel.\n",
        "\n",
        "In the `Files` panel, you should see one folder named `gdrive`. Click the arrow next to the `gdrive` folder to expand it, and continue navigating through the folders until you find the `Genomics Data Science Project - example analyses inputs` folder copied previously. Right-click on `Genomics Data Science Project - example analyses inputs` and select `Copy path`. Paste the value into the cell below, and compare it to the value assigned to `INPUT_DATA_DIR_PATH`. If the values are different, replace the pre-coded value with the one you pasted.\n",
        "\n",
        "This notebook is configured to store results in a folder named `Genomics Data Science Project - example analyses outputs` alongside the folder of input data. The folder will be created if it does not already exist. If you would rather store the results elsewhere, you can change the value of `OUTPUT_DATA_DIR_PATH` below.  \n",
        "\n",
        "Once you have made any changes to `INPUT_DATA_DIR_PATH` and `OUTPUT_DATA_DIR_PATH`, run the cell.\n",
        "\n",
        "If at any point you open the `Files` panel or click its `REFRESH` button and do not see `gdrive`, you might need to re-run the previous cell.\n",
        "\n",
        "Note this cell also specifies the output directories that will be used for the different analyses in the example. They are defined here, in the last quick-running cell before the analyses below, because the variables will need to be refreshed if your notebook connects to a new virtual machine.\n",
        "\n",
        "You **will** need to re-run the cell whenever your notebook connects to a new virtual machine, but the values assigned to `INPUT_DATA_DIR_PATH` and `OUTPUT_DATA_DIR_PATH` will not change unless you move the folders within your Google Drive.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HnBOrd8rMm2A",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "INPUT_DATA_DIR_PATH = '/content/gdrive/My Drive/Genomics Data Science Project - example analyses inputs'\n",
        "OUTPUT_DATA_DIR_PATH = os.path.join(\\\n",
        "    os.path.dirname(INPUT_DATA_DIR_PATH),\n",
        "    'Genomics Data Science Project - example analyses outputs')\n",
        "os.makedirs(OUTPUT_DATA_DIR_PATH, exist_ok=True)\n",
        "\n",
        "CLUSTERING1_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'clustering1')\n",
        "CLUSTERING2_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'clustering2')\n",
        "CLUSTERING3_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'clustering3')\n",
        "CLUSTERING4_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'clustering4')\n",
        "CLUSTERING5_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'clustering5')\n",
        "CLUSTERING6_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'clustering6')\n",
        "\n",
        "PRIORITIZATION1_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'prioritization1')\n",
        "PRIORITIZATION2_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'prioritization2')\n",
        "PRIORITIZATION3_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'prioritization3')\n",
        "\n",
        "CHARACTERIZATION1_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'characterization1')\n",
        "CHARACTERIZATION2_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'characterization2')\n",
        "CHARACTERIZATION3_DIR_PATH = os.path.join(OUTPUT_DATA_DIR_PATH, 'characterization3')\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zWsE083wAk-2",
        "colab_type": "text"
      },
      "source": [
        "## Loading Pre-Computed Results (optional)\n",
        "\n",
        "The following sections will guide you through a series of analyses. Some of the steps are computationally intensive and require more than a few minutes to run. To avoid waiting, you may load pre-computed results into your `OUTPUT_DATA_DIR_PATH` at this point, by running the cell below and waiting until it finishes, which will take about 7 minutes. You will then be able to inspect the data, and you will still be able to run any of the remaining cells if you wish.\n",
        "\n",
        "You **will not** need to re-run this cell whenever your notebook connects to a new virtual machine."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fVsPsLnlCDJk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "src_top_dir = os.path.join(INPUT_DATA_DIR_PATH, 'precomputed results')\n",
        "dst_top_dir = OUTPUT_DATA_DIR_PATH\n",
        "\n",
        "# shutil.copytree's dirs_exist_ok not introduced until py3.8, but colab uses 3.6\n",
        "for root, dirs, files in os.walk(src_top_dir):\n",
        "    for src_dir in dirs:\n",
        "        src_dir_path = os.path.join(root, src_dir)\n",
        "        dst_dir_path = src_dir_path.replace(src_top_dir, dst_top_dir, 1)\n",
        "        os.makedirs(dst_dir_path, exist_ok=True)\n",
        "    for src_file in files:\n",
        "        src_file_path = os.path.join(root, src_file)\n",
        "        dst_file_path = src_file_path.replace(src_top_dir, dst_top_dir, 1)\n",
        "        shutil.copy2(src_file_path, dst_file_path)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3t-CbUY6sP2K",
        "colab_type": "text"
      },
      "source": [
        "## Clustering\n",
        "\n",
        "The following four cells will use standard clustering techniques to group samples according to different omics data. Run each cell; note each cell contains a comment with an estimated running time. You'll see some output describing the inputs and results.\n",
        "\n",
        "As each of these cells finishes, it will store the results to your Google Drive. For that reason, you **will not** need to re-run these cells whenever your notebook connects to a new virtual machine."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ub32Jwwksg31",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes about 9 minutes\n",
        "do_clustering(\\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering1_genecopynumber.tsv'), \\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering_clinical_data.tsv'), \\\n",
        "    CLUSTERING1_DIR_PATH, 8, None, None, None, 0, None)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZKX0UCvmskVz",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes about 3 minutes\n",
        "do_clustering(\\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering2_exp_HiSeqV2.tsv'), \\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering_clinical_data.tsv'), \\\n",
        "    CLUSTERING2_DIR_PATH, 13, None, None, None, 0, None)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bV7S9H7FsngW",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes less than 1 minute\n",
        "do_clustering(\\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering3_hMethyl.tsv'), \\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering_clinical_data.tsv'), \\\n",
        "    CLUSTERING3_DIR_PATH, 19, None, None, None, 0, None)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BFJcuLclspo5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes less than 1 minute\n",
        "do_clustering(\\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering4_RPPA_RBN.tsv'), \\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering_clinical_data.tsv'), \\\n",
        "    CLUSTERING4_DIR_PATH, 8, None, None, None, 0, None)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NWMbtPJZssUT",
        "colab_type": "text"
      },
      "source": [
        "### Network-Based Clustering\n",
        "\n",
        "This fifth clustering analysis incorporates the knowledge network in order to improve results over sparse omics data. As with the above clustering analyses, run the cell and wait until it completes.\n",
        "\n",
        "As with all of the analysis cells, it will store the results to your Google Drive. For that reason, you **will not** need to re-run this cell whenever your notebook connects to a new virtual machine."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3mnPXgrTss9S",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes almost 2 hours\n",
        "do_clustering(\\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering5_mutation.tsv'), \\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering_clinical_data.tsv'), \\\n",
        "    CLUSTERING5_DIR_PATH, 14, '9606', \\\n",
        "    '/network/Gene/9606/hn_IntNet/9606.hn_IntNet.edge', 0.5, 0, None)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VwpR_2p-bwjo",
        "colab_type": "text"
      },
      "source": [
        "### Cluster-of-Clusters Analysis (COCA)\n",
        "\n",
        "This sixth clustering analysis operates upon the cluster assignments generated by the previous five clustering analyses, along with the results of an additional clustering based on miRNA data. The additional clustering is not part of this notebook, but its results are provided in the folder of inputs so that this cell can load them. Again, run the cell and wait until it finishes.\n",
        "\n",
        "As with all of the analysis cells, it will store the results to your Google Drive. For that reason, you **will not** need to re-run this cell whenever your notebook connects to a new virtual machine."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "A8Kdpf2McHdS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes about 8 minutes\n",
        "\n",
        "# gather the outputs from the five previous clustering analyses, along with the\n",
        "# miRNA clusters\n",
        "raw_coca_inputs = [\n",
        "    get_path_to_newest_file_having_prefix(CLUSTERING1_DIR_PATH, 'samples_label_by_cluster'),\n",
        "    get_path_to_newest_file_having_prefix(CLUSTERING2_DIR_PATH, 'samples_label_by_cluster'),\n",
        "    get_path_to_newest_file_having_prefix(CLUSTERING3_DIR_PATH, 'samples_label_by_cluster'),\n",
        "    get_path_to_newest_file_having_prefix(CLUSTERING4_DIR_PATH, 'samples_label_by_cluster'),\n",
        "    get_path_to_newest_file_having_prefix(CLUSTERING5_DIR_PATH, 'samples_label_by_cluster'),\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering6_miRNA.tsv')\n",
        "]    \n",
        "\n",
        "# assemble the raw inputs into a single file formatted like an omics file\n",
        "os.makedirs(CLUSTERING6_DIR_PATH, exist_ok=True)\n",
        "coca_input_file_path = os.path.join(CLUSTERING6_DIR_PATH, 'input.tsv')\n",
        "temp_dir_path = mkdtemp()\n",
        "try:\n",
        "    for input in raw_coca_inputs:\n",
        "        shutil.copy(input, temp_dir_path)\n",
        "    coca_input_df = get_cluster_binary_dataframe(\\\n",
        "        [os.path.basename(input) for input in raw_coca_inputs], temp_dir_path).T\n",
        "    coca_input_df.to_csv(coca_input_file_path, sep='\\t')\n",
        "finally:\n",
        "    shutil.rmtree(temp_dir_path)\n",
        "\n",
        "do_clustering(\\\n",
        "    coca_input_file_path, \\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering_clinical_data.tsv'), \\\n",
        "    CLUSTERING6_DIR_PATH, 13, None, None, None, 200, 0.8)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "icBaJd0bG5eU",
        "colab_type": "text"
      },
      "source": [
        "We can also perform a survival analysis on the identified clusters."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iU2YPY6TG4xJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes about 1 second\n",
        "\n",
        "# load the file containing the clinical data\n",
        "# survival data are found in columns _OS_IND (boolean representing event) and\n",
        "# _OS (float indicating time)\n",
        "phenotype_df = pd.read_csv(\\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_clustering_clinical_data.tsv'), \\\n",
        "    sep='\\t', index_col=0, header=0)\n",
        "\n",
        "# load the cluster assignments\n",
        "cluster_labels_file_path = get_path_to_newest_file_having_prefix(\\\n",
        "    CLUSTERING6_DIR_PATH, 'samples_label_by_cluster')\n",
        "cluster_labels_df = pd.read_csv(cluster_labels_file_path, \\\n",
        "    sep='\\t', index_col=0, header=None, names=['cluster'])\n",
        "# reorder cluster_labels_df to match the sample order in phenotype_df\n",
        "\n",
        "combined_df = pd.concat([phenotype_df['_OS_IND'], phenotype_df['_OS'], \\\n",
        "        cluster_labels_df], axis=1, sort=True)\n",
        "\n",
        "# retain only the samples that have a time and a cluster\n",
        "combined_df.dropna(subset=['_OS', 'cluster'], inplace=True)\n",
        "\n",
        "# fill missing values in event (0, for censored)\n",
        "combined_df['_OS_IND'].fillna(value=0, inplace=True)\n",
        "\n",
        "# calculate p-value\n",
        "test_stats = multivariate_logrank_test(combined_df['_OS'].values, \\\n",
        "    combined_df['cluster'].values, combined_df['_OS_IND'].values)\n",
        "\n",
        "# draw plot\n",
        "fig = plt.figure()\n",
        "ax = fig.gca()\n",
        "\n",
        "kmf = KaplanMeierFitter()\n",
        "\n",
        "for name, grouped_df in combined_df.groupby('cluster'):\n",
        "    kmf.fit(grouped_df[\"_OS\"], grouped_df[\"_OS_IND\"], \\\n",
        "        label='Cluster ' + str(int(name)))\n",
        "    kmf.plot(ax=ax, show_censors=True)\n",
        "\n",
        "plt.title ('P-value = %s' %(test_stats.p_value))\n",
        "plt.xlabel('Time (days)');"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6mRscDmTd9UF",
        "colab_type": "text"
      },
      "source": [
        "## Gene Prioritization\n",
        "\n",
        "The following three cells will analyze gene expression data to determine the genes most associated with phenotypes of interest.\n",
        "\n",
        "In the first of the three prioritization cells, the phenotypes are PANCAN disease types, and the method is a standard prioritization technique. Run the cell and wait until it finishes.\n",
        "\n",
        "As with all of the analysis cells, it will store the results to your Google Drive. For that reason, you **will not** need to re-run this cell whenever your notebook connects to a new virtual machine."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oTXqLBcSfOgj",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes about 6 minutes\n",
        "do_prioritization(\\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_prioritization_expr.tsv'), \\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'pancan_disease_types.gXc'), \\\n",
        "    PRIORITIZATION1_DIR_PATH, 't_test', 'average', 100, 50, None, None, None)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UNkmEKAbfq7l",
        "colab_type": "text"
      },
      "source": [
        "In the second of the three prioritization cells, the phenotypes are again the PANCAN disease types, but the method incorporates the knowledge network. Again, run the cell and wait until it finishes.\n",
        "\n",
        "As with all of the analysis cells, it will store the results to your Google Drive. For that reason, you **will not** need to re-run this cell whenever your notebook connects to a new virtual machine."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dugvbIPxgBkq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes about 14 minutes\n",
        "do_prioritization(\\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_prioritization_expr.tsv'), \\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'pancan_disease_types.gXc'), \\\n",
        "    PRIORITIZATION2_DIR_PATH, 't_test', 'average', 100, 50, '9606', \\\n",
        "    '/network/Gene/9606/hn_IntNet/9606.hn_IntNet.edge', 0.5)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Mn1YlqtogXoX",
        "colab_type": "text"
      },
      "source": [
        "In the third of the three prioritization cells, the phenotypes are the COCA cluster assignments. The standard method is used.\n",
        "\n",
        "Run the cell and wait until it finishes.\n",
        "\n",
        "As with all of the analysis cells, it will store the results to your Google Drive. For that reason, you **will not** need to re-run this cell whenever your notebook connects to a new virtual machine."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Y9E9efqLgnUb",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes about 6 minutes\n",
        "do_prioritization(\\\n",
        "    os.path.join(INPUT_DATA_DIR_PATH, 'in_prioritization_expr.tsv'), \\\n",
        "    get_path_to_newest_file_having_prefix(CLUSTERING6_DIR_PATH, 'samples_label_by_cluster'), \\\n",
        "    PRIORITIZATION3_DIR_PATH, 't_test', 'average', 100, 50, None, None, None)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KWH5pdx8lVhh",
        "colab_type": "text"
      },
      "source": [
        "## Gene-Set Characterization\n",
        "\n",
        "The final three cells compare the top genes found by the gene prioritization analyses with gene sets from the Gene Ontology database and a set of known cancer drivers.\n",
        "\n",
        "Each cell in this sequence corresponds to one of the three gene prioritization analyses above. Run each cell and wait for it to finish.\n",
        "\n",
        "As with all of the analysis cells, these will store the results to your Google Drive. For that reason, you **will not** need to re-run these cells whenever your notebook connects to a new virtual machine."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sqoPge_4mcHp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes less than 1 minute\n",
        "for property_edge_file_path in [\\\n",
        "        '/network/Property/9606/gene_ontology/9606.gene_ontology.edge', \\\n",
        "        '/network/Property/9606/cancer_driver_genes/9606.cancer_driver_genes.edge']:\n",
        "    do_characterization(\\\n",
        "        get_path_to_newest_file_having_prefix(PRIORITIZATION1_DIR_PATH, 'top_'), \\\n",
        "        os.path.join(CHARACTERIZATION1_DIR_PATH, os.path.basename(property_edge_file_path)), \\\n",
        "        '9606', property_edge_file_path, None, None)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uEyycIpIFa9N",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes less than 1 minute\n",
        "for property_edge_file_path in [\\\n",
        "        '/network/Property/9606/gene_ontology/9606.gene_ontology.edge', \\\n",
        "        '/network/Property/9606/cancer_driver_genes/9606.cancer_driver_genes.edge']:\n",
        "    do_characterization(\\\n",
        "        get_path_to_newest_file_having_prefix(PRIORITIZATION2_DIR_PATH, 'top_'), \\\n",
        "        os.path.join(CHARACTERIZATION2_DIR_PATH, os.path.basename(property_edge_file_path)), \\\n",
        "        '9606', property_edge_file_path, None, None)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lGDg9MiMFepB",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# takes less than 1 minute\n",
        "for property_edge_file_path in [\\\n",
        "        '/network/Property/9606/gene_ontology/9606.gene_ontology.edge', \\\n",
        "        '/network/Property/9606/cancer_driver_genes/9606.cancer_driver_genes.edge']:\n",
        "    do_characterization(\\\n",
        "        get_path_to_newest_file_having_prefix(PRIORITIZATION3_DIR_PATH, 'top_'), \\\n",
        "        os.path.join(CHARACTERIZATION3_DIR_PATH, os.path.basename(property_edge_file_path)), \\\n",
        "        '9606', property_edge_file_path, None, None)"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}