4_Project_CHL.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kLse5cbrtAn0",
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "**Computational Health Laboratory Project, A.Y. 2021/2022**\n",
    "\n",
    "**Authors:** Niko Dalla Noce, Alessandro Ristori, Andrea Zuppolini"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EnunYcRW8FBh",
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# **CHL Project, Pathway Analysis**\n",
    "Starting fron one or more genes, extract from interaction databases the genes they interact with. Using the expanded gene set, perform pathway analysis and obtain all disease pathways in which the genes appear. Merge the pathways to obtain a larger graph. Perform further network analysis to extract central biomarkers and communities beyond pathways. Compute a distance between the initial gene set and the various pathways (diseases)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VBEF2viCEfdx",
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## **Colab setup**\n",
    "Takes care of the project setup on Colab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "a5orYga8hr_f",
    "outputId": "af8b10d2-1f8b-4ed7-c3df-db1356806b15",
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "if 'google.colab' in str(get_ipython()):\n",
    "    import subprocess\n",
    "    from google.colab import drive\n",
    "    out_clone = subprocess.run([\"git\", \"clone\", \"https://github.com/nikodallanoce/ComputationalHealthLaboratory\"], text=True, capture_output=True)\n",
    "    print(\"{0}{1}\".format(out_clone.stdout, out_clone.stderr))\n",
    "    %pip install -U PyYAML\n",
    "    %pip install gseapy\n",
    "    drive.mount(\"/content/drive/\")\n",
    "    %cp \"/content/drive/Shareddrives/CHL/config.yml\" \"/content/ComputationalHealthLaboratory\"\n",
    "    %cd ComputationalHealthLaboratory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zExTy0EvuEaV",
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## **Obtain all the genes that interacts with the starting one**\n",
    "Starting from a gene, retrieve its neighbours and the interactions between them.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [],
   "source": [
    "import requests\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from config import ACCESS_KEY, BASE_URL"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Some utilities methods needed to create and clean the interactions dataset."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [],
   "source": [
    "def build_starting_gene_interactions(interactions_dataset: str) -> pd.DataFrame:\n",
    "    # Build the interactions dataframe, choose and rename the needed columns\n",
    "    gene_starting_interactions = pd.read_table(interactions_dataset)\n",
    "    gene_starting_interactions.rename(columns={\"Official Symbol Interactor A\": \"InteractorA\",\n",
    "                                      \"Official Symbol Interactor B\": \"InteractorB\"}, inplace=True)\n",
    "    gene_starting_interactions = gene_starting_interactions[[\"InteractorA\", \"InteractorB\"]]\n",
    "\n",
    "    # Put uppercase all the genes inside each interaction\n",
    "    gene_starting_interactions[\"InteractorA\"] = gene_starting_interactions[\"InteractorA\"].str.upper()\n",
    "    gene_starting_interactions[\"InteractorB\"] = gene_starting_interactions[\"InteractorB\"].str.upper()\n",
    "    return gene_starting_interactions\n",
    "\n",
    "\n",
    "def remove_duplicated_interactions(interactions_dataset: pd.DataFrame, verbose: bool = True) -> pd.DataFrame:\n",
    "    # Look for duplicated interactions\n",
    "    duplicated_interactions = pd.DataFrame(np.sort(\n",
    "        interactions_dataset[[\"InteractorA\", \"InteractorB\"]].values, 1)).duplicated()\n",
    "    if verbose:\n",
    "        print(\"Duplicated interactions:\\n{0}\".format(duplicated_interactions.value_counts()))\n",
    "\n",
    "    # Delete such interactions from the dataset\n",
    "    cleaned_interactions_dataset = interactions_dataset[~duplicated_interactions.values]\n",
    "    return cleaned_interactions_dataset\n",
    "\n",
    "\n",
    "def remove_self_loop_interactions(interactions_dataset: pd.DataFrame, verbose: bool = True) -> pd.DataFrame:\n",
    "    # Look for interactions where both proteins are the same\n",
    "    same_proteins_interactions = pd.Series(interactions_dataset[[\"InteractorA\", \"InteractorB\"]].nunique(axis=1) == 1)\n",
    "    if verbose:\n",
    "        print(\"Useless interactions:\\n{0}\".format(same_proteins_interactions.value_counts()))\n",
    "\n",
    "    # Delete such interactions from the dataset\n",
    "    cleaned_interactions_dataset = interactions_dataset[~same_proteins_interactions.values]\n",
    "    return cleaned_interactions_dataset"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "We expect that you have retrieved the gene interactions from BioGRID in the tab3 format, if so, then build the pandas dataframe with the starting gene’s interactions."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "data": {
      "text/plain": "    InteractorA InteractorB\n149        NSP8         SON\n150        NSP9         SON\n151        ORF6         SON\n152        ORF8         SON\n153        CCNF         SON",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>InteractorA</th>\n      <th>InteractorB</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>149</th>\n      <td>NSP8</td>\n      <td>SON</td>\n    </tr>\n    <tr>\n      <th>150</th>\n      <td>NSP9</td>\n      <td>SON</td>\n    </tr>\n    <tr>\n      <th>151</th>\n      <td>ORF6</td>\n      <td>SON</td>\n    </tr>\n    <tr>\n      <th>152</th>\n      <td>ORF8</td>\n      <td>SON</td>\n    </tr>\n    <tr>\n      <th>153</th>\n      <td>CCNF</td>\n      <td>SON</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "starting_interactions = build_starting_gene_interactions(\"datasets/BIOGRID.tab3.txt\")\n",
    "starting_interactions.tail()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "We need to remove duplicated and self-loops interactions from the dataframe."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "starting_interactions = remove_duplicated_interactions(starting_interactions)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "starting_interactions = remove_self_loop_interactions(starting_interactions)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Save the cleaned dataframe into a csv file."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "starting_interactions.to_csv(\"datasets/geneset.csv\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Obtain all the proteins that interact with our starting gene."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "proteins_list = set(starting_interactions[\"InteractorA\"])\n",
    "proteins_list.update(starting_interactions[\"InteractorB\"])\n",
    "proteins_list.remove(\"SON\")  # Do not consider the starting gene\n",
    "proteins_list = list(proteins_list)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## **Expand the interactions dataset**\n",
    "Expand the dataset using the proteins obtained from the previous step."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [],
   "source": [
    "request_url = BASE_URL + \"/interactions\"\n",
    "data = {}\n",
    "\n",
    "step = 5\n",
    "for i in range(0, len(proteins_list), step):\n",
    "    end = i+step\n",
    "    if end >= len(proteins_list):\n",
    "        end = len(proteins_list)\n",
    "    \n",
    "    # List of genes to search for\n",
    "    gene_list = proteins_list[i:end]\n",
    "\n",
    "    params = {\n",
    "        \"accesskey\": ACCESS_KEY,\n",
    "        \"format\": \"json\",  # Return results in TAB2 format\n",
    "        \"geneList\": \"|\".join(gene_list),  # Must be | separated\n",
    "        \"searchNames\": \"true\",  # Search against official names\n",
    "        \"includeInteractors\": \"true\",  # Set to true to get any interaction involving EITHER gene, set to false to get interactions between genes\n",
    "        \"includeInteractorInteractions\": \"false\",  # Set to true to get interactions between the geneList’s first order interactors\n",
    "        \"includeEvidence\": \"false\",  # If false \"evidenceList\" is evidence to exclude, if true \"evidenceList\" is evidence to show\n",
    "        \"selfInteractionsExcluded\": \"true\", # If true no self-interactions will be included\n",
    "    }\n",
    "\n",
    "    r = requests.get(request_url, params=params)\n",
    "    interactions = r.json()\n",
    "    \n",
    "    # Check if the interactions are more than the allowed number\n",
    "    if len(interactions)==10000:\n",
    "      assert False\n",
    "\n",
    "    # Create a hash of results by interaction identifier\n",
    "    for interaction_id, interaction in interactions.items():\n",
    "        data[interaction_id] = interaction"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "data": {
      "text/plain": "        InteractorA InteractorB\n3305885        CCNF       ZBTB1\n3305886        CCNF       ZGPAT\n3305887        CCNF      ZNF638\n3305888        CCNF      ZNF687\n3305889        CCNF       ZWINT",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>InteractorA</th>\n      <th>InteractorB</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>3305885</th>\n      <td>CCNF</td>\n      <td>ZBTB1</td>\n    </tr>\n    <tr>\n      <th>3305886</th>\n      <td>CCNF</td>\n      <td>ZGPAT</td>\n    </tr>\n    <tr>\n      <th>3305887</th>\n      <td>CCNF</td>\n      <td>ZNF638</td>\n    </tr>\n    <tr>\n      <th>3305888</th>\n      <td>CCNF</td>\n      <td>ZNF687</td>\n    </tr>\n    <tr>\n      <th>3305889</th>\n      <td>CCNF</td>\n      <td>ZWINT</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load the data into a pandas dataframe\n",
    "dataset = pd.DataFrame.from_dict(data, orient=\"index\")\n",
    "\n",
    "# Re-order the columns and select only the columns we want to see\n",
    "columns = [\"OFFICIAL_SYMBOL_A\", \"OFFICIAL_SYMBOL_B\"]\n",
    "dataset = dataset[columns]\n",
    "\n",
    "# Rename the columns and make all the values uppercase\n",
    "dataset = dataset.rename(columns={\"OFFICIAL_SYMBOL_A\": \"InteractorA\", \"OFFICIAL_SYMBOL_B\": \"InteractorB\"})\n",
    "dataset[\"InteractorA\"] = dataset[\"InteractorA\"].str.upper()\n",
    "dataset[\"InteractorB\"] = dataset[\"InteractorB\"].str.upper()\n",
    "\n",
    "# Print the dataframe\n",
    "dataset.tail()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Drop duplicated interactions, they're not interesting from out point of view."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Duplicated interactions:\n",
      "False    79296\n",
      "True     25281\n",
      "dtype: int64\n"
     ]
    }
   ],
   "source": [
    "dataset = remove_duplicated_interactions(dataset)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Drop self-loops since they're useless for our analysis."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Useless interactions:\n",
      "False    79283\n",
      "True        13\n",
      "dtype: int64\n"
     ]
    }
   ],
   "source": [
    "dataset = remove_self_loop_interactions(dataset)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Unify the interactions from the starting gene with the ones obtained by the requests to the BioGrid dataset."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [],
   "source": [
    "dataset = pd.concat([dataset, starting_interactions])"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of nodes: 13010\n"
     ]
    }
   ],
   "source": [
    "nodes = pd.concat([dataset[\"InteractorA\"], dataset[\"InteractorB\"]]).unique()\n",
    "print(\"Number of nodes: {0}\".format(len(nodes)))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "At last, save the interactions and nodes into csv file for pathway enrichment."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [],
   "source": [
    "# Save interactions and nodes dataset to csv\n",
    "dataset.to_csv(\"datasets/interactions.csv\")\n",
    "pd.DataFrame(nodes).to_csv(\"datasets/genes.csv\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## **Pathway enrichment**\n",
    "Find all the diseases linked to the nodes retrieved by the previous step."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "outputs": [],
   "source": [
    "import gseapy as gp"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Obtain all the pathways connected to our nodes, for our case we are going to use the DisGeNET dataset."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "if os.path.exists(\"datasets/diseases_pathways.csv\"):\n",
    "    df_diseases = pd.read_csv(\"datasets/diseases_pathways.csv\", sep=\",\", index_col=0)\n",
    "else:\n",
    "    enr = gp.enrichr(gene_list=pd.DataFrame(nodes),\n",
    "                      gene_sets=['DisGeNET'],  # Datasets from the gp.get_library_name() method\n",
    "                      organism='Human',\n",
    "                      description='DEGs_up_1d',\n",
    "                      outdir='test'\n",
    "                  )\n",
    "\n",
    "    # Keep those pathways with an adjusted p-value < 0.1\n",
    "    df_diseases = enr.results[enr.results[\"Adjusted P-value\"] < 0.1][[\"Term\", \"Overlap\", \"P-value\", \"Adjusted P-value\", \"Genes\"]]"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "outputs": [
    {
     "data": {
      "text/plain": "                                      Term Overlap   P-value  \\\n584                   Chronic otitis media   55/69  0.005896   \n585  Inadequate arch length for tooth size   47/58  0.005953   \n586                         Tooth Crowding   47/58  0.005953   \n587       Tooth mass arch size discrepancy   47/58  0.005953   \n588                 Tooth size discrepancy   47/58  0.005953   \n\n     Adjusted P-value                                              Genes  \n584          0.098950  IGHM;CD81;WIPF1;FMR1;DOCK8;CHD7;JMJD1C;COMT;GT...  \n585          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...  \n586          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...  \n587          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...  \n588          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...  ",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>Term</th>\n      <th>Overlap</th>\n      <th>P-value</th>\n      <th>Adjusted P-value</th>\n      <th>Genes</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>584</th>\n      <td>Chronic otitis media</td>\n      <td>55/69</td>\n      <td>0.005896</td>\n      <td>0.098950</td>\n      <td>IGHM;CD81;WIPF1;FMR1;DOCK8;CHD7;JMJD1C;COMT;GT...</td>\n    </tr>\n    <tr>\n      <th>585</th>\n      <td>Inadequate arch length for tooth size</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n    </tr>\n    <tr>\n      <th>586</th>\n      <td>Tooth Crowding</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n    </tr>\n    <tr>\n      <th>587</th>\n      <td>Tooth mass arch size discrepancy</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n    </tr>\n    <tr>\n      <th>588</th>\n      <td>Tooth size discrepancy</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_diseases.tail()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Save the pathways in a csv file just like for the interactions and nodes."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "outputs": [],
   "source": [
    "df_diseases.to_csv(\"datasets/diseases_pathways.csv\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Build a dict with all the diseases, this will be helpful when we will need to work on the graph."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "outputs": [],
   "source": [
    "diseases = dict()\n",
    "\n",
    "for i, disease in df_diseases.iterrows():\n",
    "    disease_genes = disease['Genes'].split(\";\")\n",
    "    term = disease['Term']\n",
    "    diseases[i] = {\"name\": term, \"genes\": disease_genes}"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## **Protein-Protein network**\n",
    "Build the protein-to-protein network and link each node to its diseases.\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "outputs": [],
   "source": [
    "import networkx as nx"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Build the graph and fill it with its nodes (the proteins coming from the dataset)."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "outputs": [],
   "source": [
    "# Build the graph\n",
    "protein_graph = nx.Graph(name='Protein Interactions Graph')\n",
    "\n",
    "# Build the nodes\n",
    "for node in nodes:\n",
    "    protein_graph.add_node(node, diseases=[])  # Each node will have a list with the disease pathways it belongs to"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Insert into the nodes their respective diseases."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "outputs": [],
   "source": [
    "for i, disease in diseases.items():\n",
    "    disease_genes = disease['genes']\n",
    "    for gene in disease_genes:\n",
    "        protein_graph.nodes[gene][\"diseases\"].append(i)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "There could be nodes without any diseases, they still need to be kept into the network."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "outputs": [],
   "source": [
    "nodes_no_disease = list()\n",
    "for node in protein_graph.nodes:\n",
    "    if len(protein_graph.nodes[node][\"diseases\"])==0:\n",
    "        nodes_no_disease.append(str(node))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Nodes without diseases: 5101\n"
     ]
    }
   ],
   "source": [
    "print(\"Nodes without diseases: {0}\".format(len(nodes_no_disease)))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Then, build the edges, it's straightforward as the nodes are known, the edges' weights will be the number of diseases in common between the nodes at each end."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "outputs": [],
   "source": [
    "def intersection(lst1, lst2):\n",
    "    inters = list()\n",
    "    if not (len(lst1) == 0 or len(lst2) == 0):\n",
    "        set1 = set(lst1)\n",
    "        inters = [elem for elem in lst2 if elem in set1]\n",
    "    return inters"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "outputs": [],
   "source": [
    "for _, interaction in dataset.iterrows():\n",
    "    first_protein, second_protein = interaction[0], interaction[1]  # Proteins involved in the interaction\n",
    "\n",
    "    # Retrieve the proteins' diseases\n",
    "    prot1_dis = protein_graph.nodes()[interaction[0]]['diseases']\n",
    "    prot2_dis = protein_graph.nodes()[interaction[1]]['diseases']\n",
    "    \n",
    "    # Build the edge\n",
    "    protein_graph.add_edge(first_protein, second_protein, weight=len(intersection(prot1_dis, prot2_dis)))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "There could be edges without any disease, as we did for the nodes, they still need to be kept."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "outputs": [],
   "source": [
    "edges_no_disease = list()\n",
    "for edge in protein_graph.edges:\n",
    "    if protein_graph.edges[edge][\"weight\"]==0:\n",
    "        edges_no_disease.append(str(edge))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Edges without diseases: 45404\n"
     ]
    }
   ],
   "source": [
    "print(\"Edges without diseases: {0}\".format(len(edges_no_disease)))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "At last, save the graph."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "outputs": [],
   "source": [
    "nx.write_gpickle(protein_graph, 'datasets/protein_graph.gpickle')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## **Disease Metrics**\n",
    "Metrics needed to compare the various diseases."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Load the graph if already built previously."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "outputs": [],
   "source": [
    "import os\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "if os.path.exists(\"datasets/protein_graph.gpickle\"):\n",
    "    protein_graph = nx.read_gpickle(\"datasets/protein_graph.gpickle\")\n",
    "elif not \"protein_graph\" in locals():\n",
    "    raise ValueError(\"It was not possible to find the graph, build it from the previous steps\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Size of the largest pathway component:** Fraction of disease proteins that lie in the disease's largest pathway component (i.e., the relative size of the largest connected component (LCC) of the disease)."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "outputs": [],
   "source": [
    "def largest_conn_comp(diseases_dict: dict) -> list:\n",
    "    lcc_score = list()\n",
    "    for _, disease in tqdm(diseases_dict.items()):\n",
    "        sub_graph = protein_graph.subgraph(disease['genes'])  # Subgraph of the current disease\n",
    "        largest_cc = max(nx.connected_components(sub_graph), key=len)\n",
    "        lcc_score.append(len(largest_cc) / len(sub_graph.nodes()))\n",
    "    \n",
    "    return lcc_score"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "outputs": [
    {
     "data": {
      "text/plain": "  0%|          | 0/589 [00:00<?, ?it/s]",
      "application/vnd.jupyter.widget-view+json": {
       "version_major": 2,
       "version_minor": 0,
       "model_id": "896c78225425472bac372ca7b84a8023"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "if \"LCC Score\" in df_diseases.columns:\n",
    "    df_diseases[\"LCC Score\"] = largest_conn_comp(diseases)\n",
    "else:\n",
    "    df_diseases.insert(len(df_diseases.columns), \"LCC Score\", largest_conn_comp(diseases), True)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "outputs": [
    {
     "data": {
      "text/plain": "                                      Term Overlap   P-value  \\\n584                   Chronic otitis media   55/69  0.005896   \n585  Inadequate arch length for tooth size   47/58  0.005953   \n586                         Tooth Crowding   47/58  0.005953   \n587       Tooth mass arch size discrepancy   47/58  0.005953   \n588                 Tooth size discrepancy   47/58  0.005953   \n\n     Adjusted P-value                                              Genes  \\\n584          0.098950  IGHM;CD81;WIPF1;FMR1;DOCK8;CHD7;JMJD1C;COMT;GT...   \n585          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...   \n586          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...   \n587          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...   \n588          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...   \n\n     LCC Score  \n584   0.018182  \n585   0.234043  \n586   0.234043  \n587   0.234043  \n588   0.234043  ",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>Term</th>\n      <th>Overlap</th>\n      <th>P-value</th>\n      <th>Adjusted P-value</th>\n      <th>Genes</th>\n      <th>LCC Score</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>584</th>\n      <td>Chronic otitis media</td>\n      <td>55/69</td>\n      <td>0.005896</td>\n      <td>0.098950</td>\n      <td>IGHM;CD81;WIPF1;FMR1;DOCK8;CHD7;JMJD1C;COMT;GT...</td>\n      <td>0.018182</td>\n    </tr>\n    <tr>\n      <th>585</th>\n      <td>Inadequate arch length for tooth size</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n      <td>0.234043</td>\n    </tr>\n    <tr>\n      <th>586</th>\n      <td>Tooth Crowding</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n      <td>0.234043</td>\n    </tr>\n    <tr>\n      <th>587</th>\n      <td>Tooth mass arch size discrepancy</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n      <td>0.234043</td>\n    </tr>\n    <tr>\n      <th>588</th>\n      <td>Tooth size discrepancy</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n      <td>0.234043</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_diseases.tail()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Distance of pathway components:** For each pair of pathway components, we calculate the average shortest path length between each set of proteins, and then, the average of this is taken over all pairs of the components."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "outputs": [],
   "source": [
    "from numpy.ma.core import mean\n",
    "\n",
    "def distance_pathway_comps(diseases_dict: dict) -> list:\n",
    "    dpc_score = list()\n",
    "    for _, disease in tqdm(diseases_dict.items()):\n",
    "        sub_graph = protein_graph.subgraph(disease['genes'])\n",
    "        conn_comps = list(nx.connected_components(sub_graph))\n",
    "        distances = list()\n",
    "        for i, comp in enumerate(conn_comps):\n",
    "            for j in range(i+1, len(conn_comps)):\n",
    "                dist = 0\n",
    "                for first_comp_protein in comp:\n",
    "                    for second_comp_protein in conn_comps[j]:\n",
    "                        dist += nx.shortest_path_length(protein_graph, source=first_comp_protein, target=second_comp_protein)\n",
    "                \n",
    "                distances.append(dist / (len(comp) * len(conn_comps[j])))\n",
    "\n",
    "        dpc_score.append(mean(distances))\n",
    "    \n",
    "    return dpc_score"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "if os.path.exists(\"datasets/mean_distances.csv\"):\n",
    "    df_mean_distances = pd.read_csv(\"datasets/mean_distances.csv\", sep=\",\", index_col=0)\n",
    "elif not \"df_mean_distances\" in locals():\n",
    "    df_mean_distances = pd.DataFrame(distance_pathway_comps(diseases))\n",
    "    df_mean_distances.to_csv('datasets/mean_distances.csv')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "outputs": [
    {
     "data": {
      "text/plain": "                                      Term Overlap   P-value  \\\n584                   Chronic otitis media   55/69  0.005896   \n585  Inadequate arch length for tooth size   47/58  0.005953   \n586                         Tooth Crowding   47/58  0.005953   \n587       Tooth mass arch size discrepancy   47/58  0.005953   \n588                 Tooth size discrepancy   47/58  0.005953   \n\n     Adjusted P-value                                              Genes  \\\n584          0.098950  IGHM;CD81;WIPF1;FMR1;DOCK8;CHD7;JMJD1C;COMT;GT...   \n585          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...   \n586          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...   \n587          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...   \n588          0.099228  AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...   \n\n     LCC Score  DPC Score  \n584   0.018182   2.678114  \n585   0.234043   2.693694  \n586   0.234043   2.693694  \n587   0.234043   2.693694  \n588   0.234043   2.693694  ",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>Term</th>\n      <th>Overlap</th>\n      <th>P-value</th>\n      <th>Adjusted P-value</th>\n      <th>Genes</th>\n      <th>LCC Score</th>\n      <th>DPC Score</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>584</th>\n      <td>Chronic otitis media</td>\n      <td>55/69</td>\n      <td>0.005896</td>\n      <td>0.098950</td>\n      <td>IGHM;CD81;WIPF1;FMR1;DOCK8;CHD7;JMJD1C;COMT;GT...</td>\n      <td>0.018182</td>\n      <td>2.678114</td>\n    </tr>\n    <tr>\n      <th>585</th>\n      <td>Inadequate arch length for tooth size</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n      <td>0.234043</td>\n      <td>2.693694</td>\n    </tr>\n    <tr>\n      <th>586</th>\n      <td>Tooth Crowding</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n      <td>0.234043</td>\n      <td>2.693694</td>\n    </tr>\n    <tr>\n      <th>587</th>\n      <td>Tooth mass arch size discrepancy</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n      <td>0.234043</td>\n      <td>2.693694</td>\n    </tr>\n    <tr>\n      <th>588</th>\n      <td>Tooth size discrepancy</td>\n      <td>47/58</td>\n      <td>0.005953</td>\n      <td>0.099228</td>\n      <td>AMER1;SETD5;NOTCH3;TRIO;RPL10;SATB2;GNAI3;PLOD...</td>\n      <td>0.234043</td>\n      <td>2.693694</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_diseases[\"DPC Score\"] = df_mean_distances\n",
    "df_diseases.tail()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Save the updated pathways with their scores."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "outputs": [],
   "source": [
    "df_diseases.to_csv(\"datasets/diseases_scores.csv\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## **Network biomarkers**\n",
    "Compute the central nodes in the graph by taking into account their normalized degree."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "outputs": [],
   "source": [
    "import networkx.algorithms.centrality as nx_c"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Compute the degree of each node and choose those that satisfy a treshold, the starting gene must be considered too."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "outputs": [
    {
     "data": {
      "text/plain": "          centrality\nKIAA1429    0.223922\nESR2        0.175801\nESR1        0.174879\nFANCD2      0.162733\nMYC         0.156815\nKIF14       0.138596\nHIST1H4A    0.127220\nBRD4        0.117534\nEED         0.112537\nCIT         0.111308\nEFTUD2      0.107233\nNR2C2       0.106772\nEP300       0.103928\nNXF1        0.102237\nHSPA8       0.100930\nTRIM28      0.097778\nVCP         0.085018\nPHB         0.084172\nKIF20A      0.080483\nPRC1        0.079176\nKIF23       0.076255\nNAA40       0.075025\nRBM39       0.071489\nSNRNP70     0.070874\nORF8        0.069644\nECT2        0.069644\nOBSL1       0.068875\nZC3H18      0.065032\nRNF2        0.064955\nSON         0.011223",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>centrality</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>KIAA1429</th>\n      <td>0.223922</td>\n    </tr>\n    <tr>\n      <th>ESR2</th>\n      <td>0.175801</td>\n    </tr>\n    <tr>\n      <th>ESR1</th>\n      <td>0.174879</td>\n    </tr>\n    <tr>\n      <th>FANCD2</th>\n      <td>0.162733</td>\n    </tr>\n    <tr>\n      <th>MYC</th>\n      <td>0.156815</td>\n    </tr>\n    <tr>\n      <th>KIF14</th>\n      <td>0.138596</td>\n    </tr>\n    <tr>\n      <th>HIST1H4A</th>\n      <td>0.127220</td>\n    </tr>\n    <tr>\n      <th>BRD4</th>\n      <td>0.117534</td>\n    </tr>\n    <tr>\n      <th>EED</th>\n      <td>0.112537</td>\n    </tr>\n    <tr>\n      <th>CIT</th>\n      <td>0.111308</td>\n    </tr>\n    <tr>\n      <th>EFTUD2</th>\n      <td>0.107233</td>\n    </tr>\n    <tr>\n      <th>NR2C2</th>\n      <td>0.106772</td>\n    </tr>\n    <tr>\n      <th>EP300</th>\n      <td>0.103928</td>\n    </tr>\n    <tr>\n      <th>NXF1</th>\n      <td>0.102237</td>\n    </tr>\n    <tr>\n      <th>HSPA8</th>\n      <td>0.100930</td>\n    </tr>\n    <tr>\n      <th>TRIM28</th>\n      <td>0.097778</td>\n    </tr>\n    <tr>\n      <th>VCP</th>\n      <td>0.085018</td>\n    </tr>\n    <tr>\n      <th>PHB</th>\n      <td>0.084172</td>\n    </tr>\n    <tr>\n      <th>KIF20A</th>\n      <td>0.080483</td>\n    </tr>\n    <tr>\n      <th>PRC1</th>\n      <td>0.079176</td>\n    </tr>\n    <tr>\n      <th>KIF23</th>\n      <td>0.076255</td>\n    </tr>\n    <tr>\n      <th>NAA40</th>\n      <td>0.075025</td>\n    </tr>\n    <tr>\n      <th>RBM39</th>\n      <td>0.071489</td>\n    </tr>\n    <tr>\n      <th>SNRNP70</th>\n      <td>0.070874</td>\n    </tr>\n    <tr>\n      <th>ORF8</th>\n      <td>0.069644</td>\n    </tr>\n    <tr>\n      <th>ECT2</th>\n      <td>0.069644</td>\n    </tr>\n    <tr>\n      <th>OBSL1</th>\n      <td>0.068875</td>\n    </tr>\n    <tr>\n      <th>ZC3H18</th>\n      <td>0.065032</td>\n    </tr>\n    <tr>\n      <th>RNF2</th>\n      <td>0.064955</td>\n    </tr>\n    <tr>\n      <th>SON</th>\n      <td>0.011223</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nodes_degree = pd.DataFrame.from_dict(nx_c.degree_centrality(protein_graph), orient=\"index\", columns=[\"centrality\"])\n",
    "nodes_degree = nodes_degree.sort_values(by='centrality', ascending=False)\n",
    "biomarkers = pd.concat([nodes_degree.iloc[:29], nodes_degree[nodes_degree.index==\"SON\"]]) # Insert the starting gene\n",
    "biomarkers"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Save the biomarkers into a csv file."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "biomarkers.to_csv(\"datasets/biomarkers.csv\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## **Community analysis**\n",
    "Find the communities inside the graph and look how the one in which the starting gene fell is connected to the various diseases."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Find all the communities with the Louvain method (beware that is not deterministic) and discard those with only one node inside."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "outputs": [],
   "source": [
    "louvain_communities = list(nx.algorithms.community.louvain_communities(protein_graph))\n",
    "communities = [community for community in louvain_communities if len(community) > 1]  # Discard those communities with one node"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of communities: 9\n"
     ]
    }
   ],
   "source": [
    "print(\"Number of communities: {0}\".format(len(communities)))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Compute the mean number of nodes of those communities kept after the pruning on their size."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "outputs": [],
   "source": [
    "def mean_size_communities(communities: list) -> float:\n",
    "    mean_size = 0\n",
    "    for community in communities:\n",
    "        mean_size += len(community)\n",
    "\n",
    "    mean_size /= len(communities)\n",
    "    return mean_size"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean size of communities: 681.1111111111111\n"
     ]
    }
   ],
   "source": [
    "print(\"Mean size of communities: {0}\".format(str(mean_size_communities(communities))))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Compute the number of diseases in a community."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "outputs": [],
   "source": [
    "def diseases_in_community(protein_graph: nx.Graph, community: set) -> set:\n",
    "    diseases_community = set()\n",
    "    for protein in list(community):\n",
    "        diseases_protein = protein_graph.nodes[protein][\"diseases\"]\n",
    "        diseases_community.update(diseases_protein)\n",
    "\n",
    "    return diseases_community"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of diseases in community 5: 585\n"
     ]
    }
   ],
   "source": [
    "diseases_comm = diseases_in_community(protein_graph, communities[5])\n",
    "print(\"Number of diseases in community {0}: {1}\".format(5, len(diseases_comm)))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's also see how many disease, on average, where linked to those communities with only one node."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "outputs": [],
   "source": [
    "def mean_diseases_communities_size_n(communities: list, protein_graph: nx.Graph, n: int = 1) -> float:\n",
    "    mean_diseases = 0\n",
    "    n_size_commmunities = 0\n",
    "    for community in communities:\n",
    "        community = list(community)\n",
    "        if len(community) == n:\n",
    "            protein = community[0]\n",
    "            protein_diseases = protein_graph.nodes[protein]['diseases']\n",
    "            n_size_commmunities += 1\n",
    "            mean_diseases += len(protein_diseases)\n",
    "\n",
    "    if n_size_commmunities == 0:\n",
    "        print(\"There are no communities with {0} nodes\".format(n))\n",
    "    else:\n",
    "        mean_diseases /= n_size_commmunities\n",
    "\n",
    "    return mean_diseases"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean diseases for those communities with one node: 2.2324127906976745\n"
     ]
    }
   ],
   "source": [
    "mean_diseases_one_node = mean_diseases_communities_size_n(louvain_communities, protein_graph)\n",
    "print(\"Mean diseases for those communities with one node: {0}\".format(str(mean_diseases_one_node)))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "At last compute the mean number of nodes for the diseases."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "outputs": [],
   "source": [
    "def mean_genes_diseases(diseases: dict) -> float:\n",
    "    mean_size = 0\n",
    "    for _, disease in diseases.items():\n",
    "        mean_size += len(disease['genes'])\n",
    "\n",
    "    mean_size /= len(diseases.keys())\n",
    "    return mean_size"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean number of genes for disease: 246.55857385398983\n"
     ]
    }
   ],
   "source": [
    "print(\"Mean number of genes for disease: {0}\".format(mean_genes_diseases(diseases)))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's see if the communities do not share any nodes."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "outputs": [],
   "source": [
    "def are_communities_distinct(communities: list) -> bool:\n",
    "    for i, first_community in enumerate(communities):\n",
    "        for j in range(i+1, len(communities)):\n",
    "            second_community = communities[j]\n",
    "            if len(intersection(first_community, second_community))>0:\n",
    "                return False\n",
    "\n",
    "    return True"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "outputs": [
    {
     "data": {
      "text/plain": "True"
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "are_communities_distinct(communities)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Compute how many proteins are shared between each community and disease."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "outputs": [],
   "source": [
    "def communities_metrics(communities: list, diseases: dict) -> pd.DataFrame:\n",
    "    df_ranks = list()\n",
    "    for i, community in enumerate(communities):\n",
    "        tot_genes = dict()\n",
    "        shared_genes = dict()\n",
    "        for k, disease in diseases.items():\n",
    "            genes = disease['genes']\n",
    "            shared_genes_community = intersection(genes, community)\n",
    "            tot_genes[k] = len(genes)\n",
    "            shared_genes[k] = len(shared_genes_community)\n",
    "\n",
    "        for j in range(len(tot_genes)):\n",
    "            n_genes, n_shared_genes = tot_genes[j], shared_genes[j]\n",
    "            if n_shared_genes > 1:\n",
    "                df_ranks.append({\"Community\": i, \"Disease\": diseases[j]['name'], \"Shared genes\": n_shared_genes,\n",
    "                                 \"Disease genes\": n_genes, \"Community size\": len(community)})\n",
    "\n",
    "    df_ranks = pd.DataFrame(df_ranks)\n",
    "\n",
    "    # Ratio of the shared genes (between community and disease pathway) and the number of genes in the disease\n",
    "    df_ranks[\"Ratio disease\"] = df_ranks['Shared genes'] / df_ranks['Disease genes']\n",
    "\n",
    "    # Ratio of the shared genes (between community and disease pathway) and the size of the community\n",
    "    df_ranks['Ratio community'] = df_ranks[\"Shared genes\"] / df_ranks[\"Community size\"]\n",
    "\n",
    "    # Relevance of the  based on the previous computed metrics\n",
    "    df_ranks[\"Relevance\"] = df_ranks[\"Ratio disease\"] * df_ranks[\"Ratio community\"]\n",
    "    return df_ranks"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "outputs": [
    {
     "data": {
      "text/plain": "      Community                 Disease  Shared genes  Disease genes  \\\n1118          1  Tooth size discrepancy             7             47   \n1668          2  Tooth size discrepancy             3             47   \n2212          3  Tooth size discrepancy             4             47   \n3315          5  Tooth size discrepancy             6             47   \n3854          6  Tooth size discrepancy             3             47   \n4398          7  Tooth size discrepancy             3             47   \n4962          8  Tooth size discrepancy            12             47   \n\n      Community size  Ratio disease  Ratio community  Relevance  \n1118             796       0.148936         0.008794   0.001310  \n1668             886       0.063830         0.003386   0.000216  \n2212             716       0.085106         0.005587   0.000475  \n3315            1132       0.127660         0.005300   0.000677  \n3854             620       0.063830         0.004839   0.000309  \n4398             450       0.063830         0.006667   0.000426  \n4962             463       0.255319         0.025918   0.006617  ",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>Community</th>\n      <th>Disease</th>\n      <th>Shared genes</th>\n      <th>Disease genes</th>\n      <th>Community size</th>\n      <th>Ratio disease</th>\n      <th>Ratio community</th>\n      <th>Relevance</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>1118</th>\n      <td>1</td>\n      <td>Tooth size discrepancy</td>\n      <td>7</td>\n      <td>47</td>\n      <td>796</td>\n      <td>0.148936</td>\n      <td>0.008794</td>\n      <td>0.001310</td>\n    </tr>\n    <tr>\n      <th>1668</th>\n      <td>2</td>\n      <td>Tooth size discrepancy</td>\n      <td>3</td>\n      <td>47</td>\n      <td>886</td>\n      <td>0.063830</td>\n      <td>0.003386</td>\n      <td>0.000216</td>\n    </tr>\n    <tr>\n      <th>2212</th>\n      <td>3</td>\n      <td>Tooth size discrepancy</td>\n      <td>4</td>\n      <td>47</td>\n      <td>716</td>\n      <td>0.085106</td>\n      <td>0.005587</td>\n      <td>0.000475</td>\n    </tr>\n    <tr>\n      <th>3315</th>\n      <td>5</td>\n      <td>Tooth size discrepancy</td>\n      <td>6</td>\n      <td>47</td>\n      <td>1132</td>\n      <td>0.127660</td>\n      <td>0.005300</td>\n      <td>0.000677</td>\n    </tr>\n    <tr>\n      <th>3854</th>\n      <td>6</td>\n      <td>Tooth size discrepancy</td>\n      <td>3</td>\n      <td>47</td>\n      <td>620</td>\n      <td>0.063830</td>\n      <td>0.004839</td>\n      <td>0.000309</td>\n    </tr>\n    <tr>\n      <th>4398</th>\n      <td>7</td>\n      <td>Tooth size discrepancy</td>\n      <td>3</td>\n      <td>47</td>\n      <td>450</td>\n      <td>0.063830</td>\n      <td>0.006667</td>\n      <td>0.000426</td>\n    </tr>\n    <tr>\n      <th>4962</th>\n      <td>8</td>\n      <td>Tooth size discrepancy</td>\n      <td>12</td>\n      <td>47</td>\n      <td>463</td>\n      <td>0.255319</td>\n      <td>0.025918</td>\n      <td>0.006617</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "communities_rank = communities_metrics(communities, diseases)\n",
    "communities_rank[communities_rank[\"Disease\"]==\"Tooth size discrepancy\"]"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "communities_rank.to_csv(\"datasets/communities_metrics.csv\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "We need to know in wich community our starting gene is."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "outputs": []
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "outputs": [],
   "source": [
    "def look_for_gene_community(protein: str, communities: list) -> int:\n",
    "    for i, community in enumerate(communities):\n",
    "        if protein in community:\n",
    "            return i\n",
    "\n",
    "    return -1"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The gene SON is in community 1\n"
     ]
    }
   ],
   "source": [
    "protein_community = \"SON\"\n",
    "gene_community = look_for_gene_community(protein_community, communities)\n",
    "if gene_community == -1:\n",
    "    raise Exception(\"The gene {0} is not in one of the communities\".format(protein_community))\n",
    "else:\n",
    "    print(\"The gene {0} is in community {1}\".format(protein_community, gene_community))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Knowing in which community our starting gene is, we can retrieve all the diseases inisde such community and rank them by their relevance."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "outputs": [],
   "source": [
    "# Retrieve the diseases inside the community in which the gene is in and keep those that share at least ten genes with the community\n",
    "disease_rank = communities_rank[(communities_rank[\"Community\"]==gene_community) &\n",
    "                             (communities_rank[\"Shared genes\"] > 10)].sort_values(by=\"Relevance\", ascending=False).drop([\"Community\", \"Shared genes\",\n",
    "                                                                                                                         \"Disease genes\", \"Community size\"], axis=1)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "outputs": [
    {
     "data": {
      "text/plain": "                           Disease  Ratio disease  Ratio community  Relevance\n1020                 Fibroid Tumor       0.053333         0.015075   0.000804\n579   Invasive carcinoma of breast       0.045455         0.017588   0.000799\n715     Adenomatous Polyposis Coli       0.044068         0.016332   0.000720\n959         Malignant mesothelioma       0.045082         0.013819   0.000623\n979                   Mesothelioma       0.039855         0.013819   0.000551",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>Disease</th>\n      <th>Ratio disease</th>\n      <th>Ratio community</th>\n      <th>Relevance</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>1020</th>\n      <td>Fibroid Tumor</td>\n      <td>0.053333</td>\n      <td>0.015075</td>\n      <td>0.000804</td>\n    </tr>\n    <tr>\n      <th>579</th>\n      <td>Invasive carcinoma of breast</td>\n      <td>0.045455</td>\n      <td>0.017588</td>\n      <td>0.000799</td>\n    </tr>\n    <tr>\n      <th>715</th>\n      <td>Adenomatous Polyposis Coli</td>\n      <td>0.044068</td>\n      <td>0.016332</td>\n      <td>0.000720</td>\n    </tr>\n    <tr>\n      <th>959</th>\n      <td>Malignant mesothelioma</td>\n      <td>0.045082</td>\n      <td>0.013819</td>\n      <td>0.000623</td>\n    </tr>\n    <tr>\n      <th>979</th>\n      <td>Mesothelioma</td>\n      <td>0.039855</td>\n      <td>0.013819</td>\n      <td>0.000551</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "disease_rank.tail()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "outputs": [],
   "source": [
    "disease_rank.to_csv(\"datasets/community_gene_metrics.csv\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Compute the distances from a communitiy to the other ones."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def communities_distance(communities_rank: pd.DataFrame,\n",
    "                         first_index: int,\n",
    "                         second_index: int,\n",
    "                         metric: str = \"Relevance\") -> float:\n",
    "    # Retrieve all the diseases linked to the communities\n",
    "    tmp_first = communities_rank[communities_rank[\"Community\"] == first_index][[\"Disease\", metric]]\n",
    "    tmp_second = communities_rank[communities_rank[\"Community\"] == second_index][[\"Disease\", metric]]\n",
    "\n",
    "    # Shared diseases between the two communities\n",
    "    shared_diseases = tmp_first.merge(tmp_second, left_on=\"Disease\", right_on=\"Disease\")\n",
    "\n",
    "    # Diseases that are not shared between the two communities\n",
    "    tmp_first = tmp_first[~tmp_first['Disease'].isin(shared_diseases[\"Disease\"])]\n",
    "    tmp_second = tmp_second[~tmp_second['Disease'].isin(shared_diseases[\"Disease\"])]\n",
    "\n",
    "    distance = 0\n",
    "    for _, disease in shared_diseases.iterrows():\n",
    "        distance += np.power(disease[metric + \"_x\"] - disease[metric + \"_y\"], 2)\n",
    "\n",
    "    for _, first_disease in tmp_first.iterrows():\n",
    "        distance += np.power(first_disease[metric], 2)\n",
    "\n",
    "    for _, second_disease in tmp_second.iterrows():\n",
    "        distance += np.power(second_disease[metric], 2)\n",
    "\n",
    "    return distance"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "son_comm = look_for_gene_community(\"SON\", communities)\n",
    "for i in range(len(communities)):\n",
    "    print(\"Distance to community {0}: {1}\".format(i, communities_distance(communities_rank, son_comm, i)))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Retrieve all the diseases associated with the starting gene."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Retrieve all the disease associated with the starting node\n",
    "son_index_diseases = protein_graph.nodes()['SON']['diseases']\n",
    "son_diseases = df_diseases.loc[son_index_diseases, 'Term']"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Look for those disease pathways where the staring gene is not present."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Retrieve their metrics relative to the community in which the starting gene is located\n",
    "community_diseases_no_son = disease_rank[~disease_rank['Disease'].isin(son_diseases)].sort_values(by='Relevance', ascending=False)\n",
    "community_diseases_no_son"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now look for those disease pathways, and their metrics, related to the gene community where such gene belongs to those pathways."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "community_diseases_son = disease_rank[disease_rank['Disease'].isin(son_diseases)]\n",
    "community_diseases_son"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Save the communities into a csv file."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "df_communities = pd.DataFrame(communities)\n",
    "df_communities.to_csv(\"datasets/communities.csv\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## **Draw the protein, disease and community graphs**\n",
    "Plot the graphs to show the results of our work."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "outputs": [],
   "source": [
    "from pyvis.network import Network"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Methods used to build the pyvis network."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "outputs": [],
   "source": [
    "def __build_network_disease__(protein_graph: nx.Graph, genes_disease: list = None, protein: str = None, plot_other_edges: bool = True) -> Network:\n",
    "    net = Network(width=1080, height=720)\n",
    "    node_index = dict()\n",
    "    for i, node in enumerate(protein_graph.nodes()):\n",
    "        shape = \"dot\"\n",
    "        size = 8\n",
    "        color = \"grey\"\n",
    "        node_index[node] = i\n",
    "        if not genes_disease is None and node in genes_disease:\n",
    "            color = \"orange\"\n",
    "            size = 16\n",
    "            shape = \"diamond\"\n",
    "\n",
    "        if node == protein:\n",
    "            color = \"slateblue\"\n",
    "            size = 32\n",
    "\n",
    "        net.add_node(i, label=node, size=size, color=color, shape=shape)\n",
    "\n",
    "    for edge_from, edge_to in protein_graph.edges():\n",
    "        if not genes_disease is None and ((edge_from == protein and edge_to in genes_disease) or (edge_from in genes_disease and edge_to == protein)) and protein not in genes_disease:\n",
    "            net.add_edge(node_index[edge_from], node_index[edge_to], color=\"blue\", value=1)\n",
    "        elif not genes_disease is None and edge_from in genes_disease and edge_to in genes_disease:\n",
    "            net.add_edge(node_index[edge_from], node_index[edge_to], color=\"orangered\", value=1)\n",
    "        elif plot_other_edges:\n",
    "            net.add_edge(node_index[edge_from], node_index[edge_to], color=\"grey\")\n",
    "\n",
    "    return net\n",
    "\n",
    "\n",
    "def __build_network_protein__(protein_graph: nx.Graph, protein: str = None, plot_other_edges: bool = True) -> Network:\n",
    "    net = Network(width=1080, height=720)\n",
    "    node_index = dict()\n",
    "    for i, node in enumerate(protein_graph.nodes()):\n",
    "        node_index[node] = i\n",
    "        if node != protein:\n",
    "            net.add_node(i, label=node, size=8)\n",
    "        else:\n",
    "            net.add_node(i, label=node, size=16, color=\"slateblue\")\n",
    "\n",
    "    for edge_from, edge_to in protein_graph.edges():\n",
    "        if edge_from == protein or edge_to == protein:\n",
    "            net.add_edge(node_index[edge_from], node_index[edge_to], color=\"blue\", value=1)\n",
    "        elif plot_other_edges:\n",
    "            net.add_edge(node_index[edge_from], node_index[edge_to], color=\"grey\")\n",
    "\n",
    "    return net"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Plot the protein graph, you can also plot a disease's genes in the same graph. If a node is a diamond, then it belongs to the disease."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "outputs": [],
   "source": [
    "def plot_protein_network(protein_graph: nx.Graph, genes_disease: list = None, biomarkers: list = None, protein: str = None, plot_other_edges: bool = True) -> None:\n",
    "    if biomarkers is not None:\n",
    "        plot_graph = protein_graph.subgraph(biomarkers)\n",
    "    else:\n",
    "        plot_graph = protein_graph.copy()\n",
    "\n",
    "    net = __build_network_disease__(plot_graph, genes_disease, protein, plot_other_edges)\n",
    "    net.toggle_drag_nodes(False)\n",
    "    net.show_buttons(['physics'])\n",
    "    net.force_atlas_2based(spring_strength=0.02)\n",
    "    net.show(\"protein_graph.html\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "outputs": [],
   "source": [
    "plot_protein_network(protein_graph, biomarkers=biomarkers.index, protein=\"SON\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Plot the disease graph, you can also plot a disease's genes on the same plot, the latters will have a diamond shape."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "outputs": [],
   "source": [
    "def plot_disease_network(protein_graph: nx.Graph, genes_disease: list, protein: str = None) -> None:\n",
    "    sub_graph = protein_graph.subgraph(genes_disease)\n",
    "    net = __build_network_protein__(sub_graph, protein)\n",
    "    net.toggle_drag_nodes(False)\n",
    "    net.show_buttons(['physics'])\n",
    "    net.force_atlas_2based(spring_strength=0.02)\n",
    "    net.show(\"disease_graph.html\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "outputs": [],
   "source": [
    "plot_disease_network(protein_graph, diseases[5][\"genes\"], \"SON\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Plot a community graph and color the chosen protein, if any is passed, and the edges linked to it."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "outputs": [],
   "source": [
    "def plot_community_protein(protein_graph: nx.Graph, communities: list, protein: str = None, plot_other_edges: bool = True) -> None:\n",
    "    if protein is not None:\n",
    "        community = communities[look_for_gene_community(protein, communities)]\n",
    "    else:\n",
    "        community = np.random.randint(0, len(communities))\n",
    "\n",
    "    sub_graph = protein_graph.subgraph(community)\n",
    "    net = __build_network_protein__(sub_graph, protein, plot_other_edges)\n",
    "    net.toggle_drag_nodes(False)\n",
    "    net.show_buttons(['physics'])\n",
    "    net.force_atlas_2based(spring_strength=0.02)\n",
    "    net.show(\"community_protein_graph.html\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "outputs": [],
   "source": [
    "plot_community_protein(protein_graph, communities, \"SON\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Plot a community and color the disease's genes passed to the method."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "outputs": [],
   "source": [
    "def plot_community_disease(protein_graph: nx.Graph, genes_disease: list, community: set, protein: str = None, plot_other_edges: bool = True) -> None:\n",
    "    sub_graph = protein_graph.subgraph(community)\n",
    "    net = __build_network_disease__(sub_graph, genes_disease, protein, plot_other_edges)\n",
    "    net.toggle_drag_nodes(False)\n",
    "    net.show_buttons(['physics'])\n",
    "    net.force_atlas_2based(spring_strength=0.02)\n",
    "    net.show(\"community_disease_graph.html\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "outputs": [],
   "source": [
    "plot_community_disease(protein_graph, diseases[228][\"genes\"], communities[look_for_gene_community(\"SON\", communities)], \"SON\", True)\n",
    "\n",
    "# Chosen disease pathway to plot\n",
    "# SON -> Intellectual disability: 5, Undergrowth: 102, Strabismus: 16\n",
    "# No SON -> Small head: 1, Epilepsy: 228, Hyperreflexia: 118\""
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Plot all the communities."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "outputs": [],
   "source": [
    "import seaborn as sn\n",
    "\n",
    "def plot_communities(protein_graph: nx.Graph, communities: list, protein: str = None) -> None:\n",
    "    palette = sn.color_palette(\"tab10\", len(communities)).as_hex()\n",
    "    net = Network(width=1080, height=720)\n",
    "    node_index = dict()\n",
    "    index = 0\n",
    "    gene_community = look_for_gene_community(protein, communities)\n",
    "    if gene_community == -1:\n",
    "        raise Warning(\"Protein {0} not found, no community will be highlighted with gold\".format(protein))\n",
    "\n",
    "    for i in range(len(communities)):\n",
    "        community = communities[i]\n",
    "        color = palette[i]\n",
    "        sub_graph = protein_graph.subgraph(community)\n",
    "        for _, node in enumerate(sub_graph.nodes()):\n",
    "            node_index[node] = index\n",
    "            if node == protein:\n",
    "                net.add_node(index, label=node, size=32)\n",
    "            else:\n",
    "                net.add_node(index, label=node, size=8)\n",
    "\n",
    "            index += 1\n",
    "\n",
    "        for edge_from, edge_to in sub_graph.edges():\n",
    "            net.add_edge(node_index[edge_from], node_index[edge_to], color=color)\n",
    "\n",
    "    net.show_buttons(['physics'])\n",
    "    net.force_atlas_2based(spring_strength=0.02)\n",
    "    net.show(\"communities.html\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "outputs": [],
   "source": [
    "plot_communities(protein_graph, communities, \"SON\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Retrieve all the genes, inside a community and belonging to a disease pathway, that interact with the starting gene."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def disease_genes_community_interacting_with(community_graph: nx.Graph, genes_disease: list = None, protein: str = None) -> list:\n",
    "    \"\"\"\n",
    "    Retrieve all the genes that belongs to a disease pathway that interact with the passed protein\n",
    "    :param community_graph: subgraph of the protein-to-protein network, represents the community in which the protein was found, we advise to use the look_for_gene_community() method to find such community\n",
    "    :param genes_disease: genes that belong to a disease pathway\n",
    "    :param protein: name of the protein\n",
    "    :return: list of genes, which lie inside the community, that interact with the protein and belong to the disease pathway\n",
    "    \"\"\"\n",
    "    genes_interacting = list()\n",
    "    for edge_from, edge_to in community_graph.edges():\n",
    "        if (edge_from in genes_disease or edge_to in genes_disease) and (edge_from == protein or edge_to == protein):\n",
    "            if not edge_from == protein:\n",
    "                genes_interacting.append(edge_from)\n",
    "            else:\n",
    "                genes_interacting.append(edge_to)\n",
    "\n",
    "    return genes_interacting"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "community = communities[look_for_gene_community(\"SON\", communities)]\n",
    "gene_interactions_disease = disease_genes_community_interacting_with(protein_graph.subgraph(community), diseases[228][\"genes\"], \"SON\")\n",
    "gene_interactions_disease"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Project_CHL.ipynb",
   "provenance": []
  },
  "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}