case_studies_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "01075f5e-3b01-4859-872c-19267c936132",
   "metadata": {},
   "source": [
    "## Case Studies Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef82d3e6-d66f-45da-8a77-17256738d22a",
   "metadata": {},
   "outputs": [],
   "source": [
    "## import python package\n",
    "import os, sys\n",
    "path_list = os.getcwd().split('/')\n",
    "index = path_list.index('KGML-xDTD')\n",
    "script_path = '/'.join(path_list[:(index+1)] + ['model_evaluation','scripts'])\n",
    "sys.path.append(script_path)\n",
    "import pandas as pd\n",
    "import eval_utilities\n",
    "from KGML_xDTD import KGML_xDTD\n",
    "\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "from networkx.drawing.nx_agraph import graphviz_layout\n",
    "from matplotlib.backends.backend_pdf import PdfPages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff7513da-7015-406a-a9e0-a7b68ab153ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "## define some Classes and functions\n",
    "\n",
    "class Args:\n",
    "    def __init__(self, use_gpu: bool = True, gpu: int = 0):\n",
    "        self.use_gpu = use_gpu\n",
    "        self.gpu = gpu\n",
    "\n",
    "def my_draw_networkx_edge_labels(\n",
    "    G,\n",
    "    pos,\n",
    "    edge_labels=None,\n",
    "    label_pos=0.5,\n",
    "    font_size=10,\n",
    "    font_color=\"k\",\n",
    "    font_family=\"sans-serif\",\n",
    "    font_weight=\"normal\",\n",
    "    alpha=None,\n",
    "    bbox=None,\n",
    "    horizontalalignment=\"center\",\n",
    "    verticalalignment=\"center\",\n",
    "    ax=None,\n",
    "    rotate=True,\n",
    "    clip_on=True,\n",
    "    rad=0\n",
    "):\n",
    "    import matplotlib.pyplot as plt\n",
    "    import numpy as np\n",
    "\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    if edge_labels is None:\n",
    "        labels = {(u, v): d for u, v, d in G.edges(data=True)}\n",
    "    else:\n",
    "        labels = edge_labels\n",
    "    text_items = {}\n",
    "    for (n1, n2), label in labels.items():\n",
    "        (x1, y1) = pos[n1]\n",
    "        (x2, y2) = pos[n2]\n",
    "        (x, y) = (\n",
    "            x1 * label_pos + x2 * (1.0 - label_pos),\n",
    "            y1 * label_pos + y2 * (1.0 - label_pos),\n",
    "        )\n",
    "        pos_1 = ax.transData.transform(np.array(pos[n1]))\n",
    "        pos_2 = ax.transData.transform(np.array(pos[n2]))\n",
    "        linear_mid = 0.5*pos_1 + 0.5*pos_2\n",
    "        d_pos = pos_2 - pos_1\n",
    "        rotation_matrix = np.array([(0,1), (-1,0)])\n",
    "        ctrl_1 = linear_mid + rad*rotation_matrix@d_pos\n",
    "        ctrl_mid_1 = 0.5*pos_1 + 0.5*ctrl_1\n",
    "        ctrl_mid_2 = 0.5*pos_2 + 0.5*ctrl_1\n",
    "        bezier_mid = 0.5*ctrl_mid_1 + 0.5*ctrl_mid_2\n",
    "        (x, y) = ax.transData.inverted().transform(bezier_mid)\n",
    "\n",
    "        if rotate:\n",
    "            # in degrees\n",
    "            angle = np.arctan2(y2 - y1, x2 - x1) / (2.0 * np.pi) * 360\n",
    "            # make label orientation \"right-side-up\"\n",
    "            if angle > 90:\n",
    "                angle -= 180\n",
    "            if angle < -90:\n",
    "                angle += 180\n",
    "            # transform data coordinate angle to screen coordinate angle\n",
    "            xy = np.array((x, y))\n",
    "            trans_angle = ax.transData.transform_angles(\n",
    "                np.array((angle,)), xy.reshape((1, 2))\n",
    "            )[0]\n",
    "        else:\n",
    "            trans_angle = 0.0\n",
    "        # use default box of white with white border\n",
    "        if bbox is None:\n",
    "            bbox = dict(boxstyle=\"round\", ec=(1.0, 1.0, 1.0), fc=(1.0, 1.0, 1.0))\n",
    "        if not isinstance(label, str):\n",
    "            label = str(label)  # this makes \"1\" and 1 labeled the same\n",
    "\n",
    "        t = ax.text(\n",
    "            x,\n",
    "            y,\n",
    "            label,\n",
    "            size=font_size,\n",
    "            color=font_color,\n",
    "            family=font_family,\n",
    "            weight=font_weight,\n",
    "            alpha=alpha,\n",
    "            horizontalalignment=horizontalalignment,\n",
    "            verticalalignment=verticalalignment,\n",
    "            rotation=trans_angle,\n",
    "            transform=ax.transData,\n",
    "            bbox=bbox,\n",
    "            zorder=1,\n",
    "            clip_on=clip_on,\n",
    "        )\n",
    "        text_items[(n1, n2)] = t\n",
    "\n",
    "    ax.tick_params(\n",
    "        axis=\"both\",\n",
    "        which=\"both\",\n",
    "        bottom=False,\n",
    "        left=False,\n",
    "        labelbottom=False,\n",
    "        labelleft=False,\n",
    "    )\n",
    "\n",
    "    return text_items    \n",
    "\n",
    "def reorganize_df(edges):\n",
    "    temp = dict()\n",
    "    for index in range(len(edges)):\n",
    "        source, target = edges.loc[index,0],edges.loc[index,1]\n",
    "        if (source, target) not in temp:\n",
    "            temp[(source, target)] = dict()\n",
    "            temp[(source, target)]['edges'] = set([edges.loc[index,2].replace('biolink:','')])\n",
    "        else:\n",
    "            temp[(source, target)]['edges'].update(set([edges.loc[index,2].replace('biolink:','')]))\n",
    "\n",
    "    return pd.DataFrame([(key[0],key[1],'/'.join(list(value['edges']))) for key, value in temp.items()])\n",
    "\n",
    "def path_graph(path_res_dict, pair, title, tp_score, is_in_train_set, is_in_not_train_set, rad=0.07):\n",
    "    edges = reorganize_df(path_res_dict[pair])\n",
    "    g = nx.DiGraph()\n",
    "    g.add_edges_from([(x[0],x[1], {'relation': x[2].replace('biolink:','')}) for x in edges.to_numpy()])\n",
    "    labels= nx.get_edge_attributes(g,'relation')\n",
    "\n",
    "    f = plt.figure(figsize=(18, 12));\n",
    "    ax = f.add_subplot(111);\n",
    "    pos = graphviz_layout(g, prog='dot', args='-Grankdir=\"LR\"');\n",
    "    nx.draw_networkx_nodes(g, pos, node_size=10, linewidths=0.5, alpha=0.9)\n",
    "    nx.draw_networkx_edges(g, pos, connectionstyle=f'arc3,rad={rad}', width=0.5)\n",
    "    nx.draw_networkx_labels(g, pos, font_size=12, font_color='blue');\n",
    "    my_draw_networkx_edge_labels(g, pos, edge_labels=labels, rotate=True, font_size=6, font_weight='bold', rad = rad);\n",
    "    plt.text(0.005, 0.06, f\"Predicted Score: {tp_score}\", size=10, color='black', fontweight='bold', transform=ax.transAxes);\n",
    "    # plt.text(0.005, 0.03, f\"Is in train set: {is_in_train_set}\", size=10, color='black', fontweight='bold', transform=ax.transAxes);\n",
    "    # plt.text(0.005, 0, f\"Is in val/test set: {is_in_not_train_set}\", size=10, color='black', fontweight='bold', transform=ax.transAxes);\n",
    "    plt.text(0.005, 0.04, f\"Drug ID in KG2c: {pair[0]} ({title.split(' - ')[0]})\", size=10, color='black', fontweight='bold', transform=ax.transAxes);\n",
    "    plt.text(0.005, 0.02, f\"Disease ID in KG2c: {pair[1]} ({title.split(' - ')[1]})\", size=10, color='black', fontweight='bold', transform=ax.transAxes);\n",
    "    return [f,g]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "547504a8-8098-48d1-bc7c-43d3dbe7f25c",
   "metadata": {},
   "outputs": [],
   "source": [
    "## set up general parameters\n",
    "args = Args()\n",
    "\n",
    "## set up data path and model path\n",
    "data_path = '/'.join(path_list[:(index+1)] + ['model_evaluation','data'])\n",
    "model_path = '/'.join(path_list[:(index+1)] + ['model_evaluation','models'])\n",
    "\n",
    "## create a KGML-xDTD object (this step needs to take around 5 minutes because its need to load the required files (e.g. KG) and trained modules)\n",
    "xdtd = KGML_xDTD(args, data_path, model_path)\n",
    "\n",
    "## read training, validation and test datasets\n",
    "train_rf_3class = pd.read_csv(os.path.join(data_path, 'train_pairs.txt'), sep='\\t', header=0)\n",
    "val_rf_3class = pd.read_csv(os.path.join(data_path, 'val_pairs.txt'), sep='\\t', header=0)\n",
    "test_rf_3class = pd.read_csv(os.path.join(data_path, 'test_pairs.txt'), sep='\\t', header=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f029e566-4b4c-4a75-a994-290432295155",
   "metadata": {},
   "source": [
    "### Case Studies Analysis with Hemophilia B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2d90987-785d-436c-a77e-dcdc1358b876",
   "metadata": {},
   "outputs": [],
   "source": [
    "## predict top 100 potential drugs that could be used to treat a given disease\n",
    "res = xdtd.predict_top_N_drugs(disease_name='Hemophilia B', N=100)\n",
    "hemophilia_curie_id = xdtd._get_preferred_curie(name='Hemophilia B', curie_type='disease')\n",
    "## take the top 10 drug-Hemophilia B pairs that includes at most 5 results in training set\n",
    "tp_in_train_set = pd.DataFrame(res['drug_id'].isin(list(train_rf_3class.loc[(train_rf_3class['target'].isin([hemophilia_curie_id])) & (train_rf_3class['y']==1),'source'])))\n",
    "tp_not_in_train_set_list = list(val_rf_3class.loc[(val_rf_3class['target'].isin([hemophilia_curie_id])) & (val_rf_3class['y']==1),'source']) + list(test_rf_3class.loc[(test_rf_3class['target'].isin([hemophilia_curie_id])) & (test_rf_3class['y']==1),'source'])\n",
    "tp_not_in_train_set = pd.DataFrame(res['drug_id'].isin(tp_not_in_train_set_list))\n",
    "tn_in_train_set = pd.DataFrame(res['drug_id'].isin(list(train_rf_3class.loc[(train_rf_3class['target'].isin([hemophilia_curie_id])) & (train_rf_3class['y']==0),'source'])))\n",
    "tn_not_in_train_set_list = list(val_rf_3class.loc[(val_rf_3class['target'].isin([hemophilia_curie_id])) & (val_rf_3class['y']==0),'source']) + list(test_rf_3class.loc[(test_rf_3class['target'].isin([hemophilia_curie_id])) & (test_rf_3class['y']==0),'source'])\n",
    "tn_not_in_train_set = pd.DataFrame(res['drug_id'].isin(tn_not_in_train_set_list))\n",
    "random_pairs_in_train_set = pd.DataFrame(res['drug_id'].isin(list(train_rf_3class.loc[(train_rf_3class['target'].isin([hemophilia_curie_id])) & (train_rf_3class['y']==2),'source'])))\n",
    "random_pairs_not_in_train_set_list = list(val_rf_3class.loc[(val_rf_3class['target'].isin([hemophilia_curie_id])) & (val_rf_3class['y']==2),'source']) + list(test_rf_3class.loc[(test_rf_3class['target'].isin([hemophilia_curie_id])) & (test_rf_3class['y']==2),'source'])\n",
    "random_pairs_not_in_train_set = pd.DataFrame(res['drug_id'].isin(random_pairs_not_in_train_set_list))\n",
    "res = pd.concat([res,tp_in_train_set,tp_not_in_train_set,tn_in_train_set,tn_not_in_train_set,random_pairs_in_train_set,random_pairs_not_in_train_set], axis=1)\n",
    "res.columns = ['drug_id', 'drug_name', 'tn_score', 'tp_score', 'unknown_score', 'tp_in_train_set', 'tp_not_in_train_set', 'tn_in_train_set','tn_not_in_train_set','random_pairs_in_train_set','random_pairs_not_in_train_set']\n",
    "temp = pd.concat([res.loc[res['tp_in_train_set'],:][:5],res.loc[~res['tp_in_train_set'],:]]).reset_index(drop=True)\n",
    "hemophilia_b_top10 = temp.loc[:9,:].reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6472cf5-f134-40e4-a003-9a464ed4dccd",
   "metadata": {},
   "outputs": [],
   "source": [
    "## predict top 10 KG-based MOA paths for each of top 10 drug-Hemophilia B pairs that includes at most 5 results in training set\n",
    "path_res = dict()\n",
    "for row in hemophilia_b_top10.to_numpy():\n",
    "    drug_curie = row[0]\n",
    "    disease_curie = hemophilia_curie_id\n",
    "    temp = xdtd.predict_top_M_moa_paths(drug_curie=drug_curie, disease_curie=disease_curie, M=10)\n",
    "    path_res[(drug_curie, disease_curie)] = temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fb4a55f",
   "metadata": {},
   "outputs": [],
   "source": [
    "## integrate 10 paths into a subgraph for each drug-disease pair\n",
    "path_res_dict = dict()\n",
    "for pair in path_res:\n",
    "    path_segment_set = set()\n",
    "    if path_res[pair]:\n",
    "        for path in path_res[pair]:\n",
    "            path_segment = path[0].split('->')\n",
    "            temp = set([(path_segment[index],path_segment[index+2],path_segment[index+1]) for index in range(0,len(path_segment)-2,2)])\n",
    "            path_segment_set.update(temp)\n",
    "        path_res_dict[pair] = pd.DataFrame(path_segment_set)\n",
    "    else:\n",
    "        print(f'No KG-based paths with length up to 3 for pair {pair}')\n",
    "\n",
    "pdf = PdfPages(\"hemophilia_b_moa.pdf\")\n",
    "for pair in path_res_dict:\n",
    "    if path_res_dict[pair].shape[0] == 0:\n",
    "        continue\n",
    "    drug_name = eval_utilities.id_to_name(pair[0]).capitalize()\n",
    "    disease_name = eval_utilities.id_to_name(pair[1]).capitalize()\n",
    "    title = f\"{drug_name} - {disease_name}\"\n",
    "    tp_score = round(hemophilia_b_top10.loc[hemophilia_b_top10['drug_id']==pair[0],'tp_score'].item(),6)\n",
    "    is_in_train_set = str(hemophilia_b_top10.loc[hemophilia_b_top10['drug_id']==pair[0],'tp_in_train_set'].to_numpy().item())\n",
    "    is_in_not_train_set = str(hemophilia_b_top10.loc[hemophilia_b_top10['drug_id']==pair[0],'tp_not_in_train_set'].to_numpy().item())\n",
    "    fig, g = path_graph(path_res_dict, pair, title, tp_score, is_in_train_set, is_in_not_train_set)\n",
    "    pdf.savefig(fig, bbox_inches='tight',dpi=200)\n",
    "pdf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1731a1c-da9b-495a-99e8-5a6670b4db04",
   "metadata": {},
   "source": [
    "### Case Studies Analysis with Huntington disease"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e4388b8-3857-4df8-8caa-0b4ce8a8055d",
   "metadata": {},
   "outputs": [],
   "source": [
    "## predict top 100 potential drugs that could be used to treat a given disease\n",
    "res = xdtd.predict_top_N_drugs(disease_name='Huntington disease', N=100)\n",
    "huntington_curie_id = xdtd._get_preferred_curie(name='Huntington disease', curie_type='disease')\n",
    "## take the top 10 drug-Huntington disease pairs that includes at most 5 results in training set\n",
    "tp_in_train_set = pd.DataFrame(res['drug_id'].isin(list(train_rf_3class.loc[(train_rf_3class['target'].isin([huntington_curie_id])) & (train_rf_3class['y']==1),'source'])))\n",
    "tp_not_in_train_set_list = list(val_rf_3class.loc[(val_rf_3class['target'].isin([huntington_curie_id])) & (val_rf_3class['y']==1),'source']) + list(test_rf_3class.loc[(test_rf_3class['target'].isin([huntington_curie_id])) & (test_rf_3class['y']==1),'source'])\n",
    "tp_not_in_train_set = pd.DataFrame(res['drug_id'].isin(tp_not_in_train_set_list))\n",
    "tn_in_train_set = pd.DataFrame(res['drug_id'].isin(list(train_rf_3class.loc[(train_rf_3class['target'].isin([huntington_curie_id])) & (train_rf_3class['y']==0),'source'])))\n",
    "tn_not_in_train_set_list = list(val_rf_3class.loc[(val_rf_3class['target'].isin([huntington_curie_id])) & (val_rf_3class['y']==0),'source']) + list(test_rf_3class.loc[(test_rf_3class['target'].isin([huntington_curie_id])) & (test_rf_3class['y']==0),'source'])\n",
    "tn_not_in_train_set = pd.DataFrame(res['drug_id'].isin(tn_not_in_train_set_list))\n",
    "random_pairs_in_train_set = pd.DataFrame(res['drug_id'].isin(list(train_rf_3class.loc[(train_rf_3class['target'].isin([huntington_curie_id])) & (train_rf_3class['y']==2),'source'])))\n",
    "random_pairs_not_in_train_set_list = list(val_rf_3class.loc[(val_rf_3class['target'].isin([huntington_curie_id])) & (val_rf_3class['y']==2),'source']) + list(test_rf_3class.loc[(test_rf_3class['target'].isin([huntington_curie_id])) & (test_rf_3class['y']==2),'source'])\n",
    "random_pairs_not_in_train_set = pd.DataFrame(res['drug_id'].isin(random_pairs_not_in_train_set_list))\n",
    "res = pd.concat([res,tp_in_train_set,tp_not_in_train_set,tn_in_train_set,tn_not_in_train_set,random_pairs_in_train_set,random_pairs_not_in_train_set], axis=1)\n",
    "res.columns = ['drug_id', 'drug_name', 'tn_score', 'tp_score', 'unknown_score', 'tp_in_train_set', 'tp_not_in_train_set', 'tn_in_train_set','tn_not_in_train_set','random_pairs_in_train_set','random_pairs_not_in_train_set']\n",
    "temp = pd.concat([res.loc[res['tp_in_train_set'],:][:5],res.loc[~res['tp_in_train_set'],:]]).reset_index(drop=True)\n",
    "## remove chemotherapy drugs\n",
    "chemotherapy_medication = ['DOCETAXEL', 'CYCLOPHOSPHAMIDE', 'CISPLATIN', 'IMATINIB', 'METHOTREXATE']\n",
    "temp = temp.loc[~temp['drug_name'].isin(chemotherapy_medication),:].reset_index(drop=True)\n",
    "huntington_top10 = temp.loc[:9,:].reset_index(drop=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c70f8384-a454-4537-9b4e-76ec192270b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "## predict top 10 KG-based MOA paths for each of top 10 drug-Huntington disease pairs that includes at most 5 results in training set\n",
    "path_res = dict()\n",
    "for row in huntington_top10.to_numpy():\n",
    "    drug_curie = row[0]\n",
    "    disease_curie = huntington_curie_id\n",
    "    temp = xdtd.predict_top_M_moa_paths(drug_curie=drug_curie, disease_curie=disease_curie, M=10)\n",
    "    path_res[(drug_curie, disease_curie)] = temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6009a0b2-8e2f-4bb9-a786-4853dd40316f",
   "metadata": {},
   "outputs": [],
   "source": [
    "## integrate 10 paths into a subgraph for each drug-disease pair\n",
    "path_res_dict = dict()\n",
    "for pair in path_res:\n",
    "    path_segment_set = set()\n",
    "    if path_res[pair]:\n",
    "        for path in path_res[pair]:\n",
    "            path_segment = path[0].split('->')\n",
    "            temp = set([(path_segment[index],path_segment[index+2],path_segment[index+1]) for index in range(0,len(path_segment)-2,2)])\n",
    "            path_segment_set.update(temp)\n",
    "        path_res_dict[pair] = pd.DataFrame(path_segment_set)\n",
    "    else:\n",
    "        print(f'No KG-based paths with length up to 3 for pair {pair}')\n",
    "\n",
    "pdf = PdfPages(\"huntington_disease_moa.pdf\")\n",
    "for pair in path_res_dict:\n",
    "    if path_res_dict[pair].shape[0] == 0:\n",
    "        continue\n",
    "    drug_name = eval_utilities.id_to_name(pair[0]).capitalize()\n",
    "    disease_name = eval_utilities.id_to_name(pair[1]).capitalize()\n",
    "    title = f\"{drug_name} - {disease_name}\"\n",
    "    tp_score = round(hemophilia_b_top10.loc[hemophilia_b_top10['drug_id']==pair[0],'tp_score'].item(),6)\n",
    "    is_in_train_set = str(hemophilia_b_top10.loc[hemophilia_b_top10['drug_id']==pair[0],'tp_in_train_set'].to_numpy().item())\n",
    "    is_in_not_train_set = str(hemophilia_b_top10.loc[hemophilia_b_top10['drug_id']==pair[0],'tp_not_in_train_set'].to_numpy().item())\n",
    "    fig, g = path_graph(path_res_dict, pair, title, tp_score, is_in_train_set, is_in_not_train_set)\n",
    "    pdf.savefig(fig, bbox_inches='tight',dpi=200)\n",
    "pdf.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3.8.12 (mypersonal_env)",
   "language": "python",
   "name": "mypersonal_env"
  },
  "language_info": {
   "name": "python",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}