CANDO_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CANDO tutorial\n",
    "*Updated January 2025*\n",
    "\n",
    "This notebook will walk you through the basics of using the Computational Analysis of Novel Drug Opportunities platform, or CANDO. This includes preparing data, setting up a CANDO object, making therapeutic predictions, and benchmarking the platform.\n",
    "\n",
    "## Table of contents\n",
    "\n",
    "- [Introduction](#Introduction)\n",
    "- [Basic tutorial](#Basic-tutorial)\n",
    "  - [Getting started](#Getting-started)\n",
    "  - [Accessing data](#Accessing-data)\n",
    "  - [Drug prediction](#Drug-prediction)\n",
    "  - [Benchmarking](#Benchmarking)\n",
    "- [Advanced topics](#Advanced-topics)\n",
    "  - [Customizing dataset](#Customizing-dataset)\n",
    "  - [Advanced drug prediction](#Advanced-drug-prediction)\n",
    "\n",
    "We recommend starting with the Getting started section, but if you wish to skip directly to a later section, or if you re-start the notebook and must re-create your CANDO object from scratch, please run the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cando as cnd # shortened to cnd for ease of typing/clarity\n",
    "import os\n",
    "\n",
    "while os.getcwd().split(os.sep)[-1] == 'tutorial': \n",
    "    os.chdir('..')\n",
    "\n",
    "if not os.path.exists('tutorial'):\n",
    "    cnd.get_tutorial()\n",
    "    cnd.get_data(v='test.0',org='tutorial')\n",
    "\n",
    "os.chdir(\"tutorial\")\n",
    "\n",
    "# Define variables with filepaths for our data\n",
    "cmpd_map='cmpds-v2.2.tsv' # Drug identifiers\n",
    "ind_map='cmpds2inds-v2.2.tsv' # Drug-indication associations\n",
    "matrix_file='tutorial_matrix-approved.tsv' # Drug-protein interaction matrix\n",
    "\n",
    "# Define additional variables used later in the tutorial; will be discussed as they appear\n",
    "dist_metric='cosine'\n",
    "ncpus=1\n",
    "\n",
    "# Create CANDO object\n",
    "cando = cnd.CANDO(cmpd_map, ind_map, matrix=matrix_file, compound_set='approved', compute_distance=True, \n",
    "                  dist_metric=dist_metric, ncpus=ncpus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This introduction will overview the fundamental ideas of the **C**omputational **A**nalysis of **N**ovel **D**rug **O**pportunities (CANDO) platform. If you want to continue directly to the practical portion of the tutorial, please skip to the [Getting started](#Getting-started) section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CANDO is a similarity-based drug discovery platform, meaning it predicts new drugs for given indications based on their similarity to existing drugs within that indication. It calculates similarity by comparing the multitarget interactions of a given compound to those of the drugs of interest. Any quantifiable interaction between a compound and biological system can be used for this purpose, and multiple types of interactions can even be combined. For this tutorial, we will be using a series of compound-protein binding scores as our interaction signature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drug-protein binding scores can be generated through a variety of methods. CANDO has multiple in-house pipelines to create such scores based on small molecule drug chemical structure and protein binding pocket data (generated by [the Zhang group's protein binding site prediction algorithm, COACH](https://zhanggroup.org/COACH/)). These pipelines, CANDOCK and BANDOCK, can be used to create a protein library or add new proteins you think might be relevant to an existing protein library. Once these binding scores are generated,  they are arranged into a vector or list corresponding to each compound: the *proteomic signature* of that compound. These proteomic signatures are arranged into a drug-protein interaction matrix. Once this matrix is created, alongside a corresponding lists of compounds and any known compound-indication associations, CANDO can be used to make novel drug predictions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To prepare for novel drug or indication prediction, CANDO first creates a compound-compound similarity matrix. Similarity scores are calculated based on the compound interaction signatures, with more similar interaction signatures having a higher similarity score. Finally, for each drug with a known effect, a sorted list of the most similar compounds to that drug is created. To predict a new drug for a given indication, the most similar compounds to all drugs already in that indication are compared. If a non-indicated drug is very similar to multiple inidcated drugs, that is, it appears in the top ranks of multiple similar lists, it will be predicted as a possible new drug in that indication. Likewise, to predict a new indication for a given compound, the indications of the drugs most similar to it are found. If multiple drugs similar to the compound have the same inidcation, the indication is returned as a potential use for that compound."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Benchmarking is useful for determining how trustworthy these predictions are. In most cases, the overall performance of any major release of CANDO will already be benchmarked. However, if you customize the set of proteins/interaction scores, indications, and/or compounds available to CANDO, or if you want to know CANDO's performance on a subset of indications or compounds of interest, re-benchmarking will allow you to estimate the trustworthiness of your results more specifically. CANDO has multiple built-in benchmarking processes to allow rapid and easy benchmarking assessments to be carried out as necessary.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic tutorial\n",
    "\n",
    "This section will go through the basics of using CANDO for simple drug discovery purposes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting started\n",
    "\n",
    "The first step of using CANDO is importing the package. We will also be importing the os module for use in the next step. \n",
    "\n",
    "If this code throws an ImportError, check to see that you have cando.py downloaded and all relevant modules are installed. CANDO can be installed using `pip install cando-py` or using an Anaconda environment. See the [CANDO Github page](https://github.com/ram-compbio/CANDO) for more installation information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cando as cnd # shortened to cnd for ease of typing/clarity\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, you must obtain all input data required to run CANDO. Most of the basic functions of CANDO require only three files:\n",
    "\n",
    "- **c_map** (compound map) - str, the filename of a file containing information on all drugs to be examined\n",
    "  - Each compound has two identifiers (a CANDO identifier and a DrugBank identifier), its name, and its approval status (e.g. approved or investigational) listed\n",
    "- **i_map** (indication map) - str, the filename of a file containing information on all indications to be examined\n",
    "  - Each line maps one CANDO compound identifier to an indication, represented by its name and two identifiers (a MeSH identifier and a CANDO identifier)\n",
    "- **matrix** - str, the filename of a file containing drug-protein binding scores (i.e. the protein interaction signature)\n",
    "  - Must contain the same number of compounds as the compound map file\n",
    "    \n",
    "In this case, we will be using real compound and indication maps, but we will use a **tutorial matrix** with a limited number of drug-protein interactions so that we can generate predictions rapidly. Other datasets can be downloaded using the command `cnd.get_data()` or [via this website](http://protinfo.compbio.buffalo.edu/cando/data/v2.2+/).\n",
    "\n",
    "If you have previously downloaded an older version of this tutorial, run `cnd.clear_cache()` before the other commands below to clear out the old tutorial data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In case this box is run 2+ times, navigates up to the first non-tutorial folder\n",
    "while os.getcwd().split(os.sep)[-1] == 'tutorial': \n",
    "    os.chdir('..')\n",
    "\n",
    "# Pull data for this tutorial\n",
    "cnd.get_tutorial()\n",
    "cnd.get_data(v='test.0',org='tutorial')\n",
    "os.chdir(\"tutorial\") # navigates to the tutorial folder\n",
    "\n",
    "# Define variables with filepaths for our data\n",
    "cmpd_map='cmpds-v2.2.tsv' # Drug identifiers\n",
    "ind_map='cmpds2inds-v2.2.tsv' # Drug-indication associations\n",
    "matrix_file='tutorial_matrix-approved.tsv' # Drug-protein interaction matrix\n",
    "\n",
    "# Define additional variables used later in the tutorial; will be discussed as they appear\n",
    "dist_metric='cosine'\n",
    "ncpus=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CANDO is an object-oriented program, and the majority of its methods require the creation of a CANDO object before they can be used. The CANDO object organizes and processes all drug, indication, and matrix data. To create a CANDO object, you will use the command `cnd.CANDO()` and pass the three files discussed above as arguments. We will be using a couple of optional arguments in addition to ensure our CANDO object is set up how we want it:\n",
    "\n",
    "1. **compound_set** - string, indicates which compounds to include in the CANDO object. This is \"all\" by default, which allows for novel drug prediction. Using \"approved\" instead allows us to focus on drug repurposing (finding new uses for existing drugs) instsead.\n",
    "2. **compute_distance** - boolean, determines whether CANDO will compute the distances between the compounds based on the similarity of their interaction signatures. Unless you have a saved distance file (we do not), you have to use compute_distance For larger protein sets. Note, either `compute_distance` or `read_dists` must be used.\n",
    "3. **save_dists** - string, takes a filename to save the distances calculated (e.g. \"big_protein_dists.tsv\"). Must be used alongside `compute_distance`. This is useful for larger protein sets (e.g. >100,000 proteins) where it might be costly to repeatedly compute these distances. *Since we are using a small protein set, we will not use this.*\n",
    "4. **read_dists** - string, takes the filename of a saved distance file (e.g. \"big_protein_dists.tsv\"). This can only be used if you have previously computed and saved distances. *Since we have not done that, we will not use this.*\n",
    "5. **dist_metric** - string, distance metric used for compute_distance. This is \"rmsd\" (root mean squared distance) by default, but \"cosine\" can also be used. This is only used if compute_distance is True.\n",
    "6. **ncpus** - integer, parallelizes certain processes across the specified number of CPUs, making them faster. You can increase or decrease this number based on the number of CPUs available to you.\n",
    "\n",
    "There are many other optional argments available for the CANDO object instantiation. If you want to know more, please consult [the CANDO documentation](https://github.com/ram-compbio/CANDO/blob/master/docs/CANDO-v2.2.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando = cnd.CANDO(cmpd_map, ind_map, matrix=matrix_file, compound_set='approved', compute_distance=True, \n",
    "                  dist_metric=dist_metric, ncpus=ncpus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating or instantiating the CANDO object, as above, makes all of the functions of CANDO available. However, before continuing onto the actual usage of CANDO, we will first cover how to access information about the major objects that make up CANDO. This will allow you to ensure that all data was integrated into the CANDO object as expected and that you are creating the predictions you intend to in later steps.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing data\n",
    "\n",
    "CANDO contains a large quantity of data, smooth access to which is essential to usage of its prediction functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CANDO consists of multiple objects, each working in tandem to create a drug prediction pipeline. The major objects are:\n",
    "\n",
    "- **CANDO object** - object, contains and organizes most variables and methods necessary to run CANDO\n",
    "- **Compound** - object, represents a single small molecule compound/drug\n",
    "- **Indication** - object, represents a single disease or other indication\n",
    "- **Protein** - object, represents a protein in the proteomic signature\n",
    "\n",
    "Other objects also exist, but are less commonly used, including **Compound_pair**, **Pathway**, and **ADR**. You will not use these objects in most applications\n",
    "\n",
    "We will start by going over the CANDO object, which we just created/instantiated. The CANDO object contains lists of all objects in your model. The CANDO object also has many other properties that might be useful; consult [the CANDO documentation](https://github.com/ram-compbio/CANDO/blob/master/docs/CANDO-v2.2.pdf) for more information.\n",
    "\n",
    "- **cando.compounds** - list, all Compound objects (taken from the compound map, `cmpd_map`)\n",
    "  - Only compounds in this list can be associated with indications or predicted as novel therapies\n",
    "- **cando.indications** - list, all Indication objects (taken from the indication map, `ind_map`)\n",
    "  - Only indications in this list can be associated with drugs or have novel therapies predicted\n",
    "- **cando.proteins** - list, all Protein objects (taken from the drug-protein interaction matrix, `matrix_file`)\n",
    "\n",
    "To start with, we will check how many objects are in each group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(cando.compounds), 'compounds')\n",
    "print(len(cando.indications), 'indications')\n",
    "print(len(cando.proteins), 'proteins')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this tutorial, our input files have 2449 compounds, 2214 indications, and 64 proteins. As a reminder, **most real matrices contain significantly more than 64 proteins**, but the protein count has been greatly reduced for the purposes of this tutorial.\n",
    "\n",
    "Let's look into the Compound objects in a little more depth. We can find a Compound if we know its CANDO ID, but we are unlikely to know that information offhand. Luckily, we can look up our compound by name using the `search_compound` function. This function will search for drugs with a similar name to your search term and return the top results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Searching for \"siprofloxasin\"...')\n",
    "cando.search_compound('siprofloxasin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Despite some mispellings, we found \"ciprofloxacin,\" the antibiotic we were looking for! Once we have our compound ID, we can use the `get_compound` function to get the corresponding Compound object without going through the cando.compounds list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cipro = cando.get_compound(421)\n",
    "print('The variable cipro contains the compound:', cipro.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Compound contains the information CANDO \"knows\" about the drug, metabolite, or small molecule it represents. This includes the following:\n",
    "\n",
    "- **compound.id_** - integer, the CANDO ID of the compound (i.e. 421)\n",
    "- **compound.name** - string, the name of the compound (i.e. ciprofloxacin)\n",
    "- **compound.status** - the clinical trial status of the compound (from the compound map)\n",
    "  - All included compounds will be \"approved\" if `compound_set = approved`; otherwise, can be \"approved\" or \"other\"\n",
    "- **compound.sig** - list, the (proteomic) signature of the drug, meaning all protein binding scores in order (from the matrix file)\n",
    "- **compound.indications** - list of string, the MeSH IDs of the indications associated with the drug (from the indication map)\n",
    "  - The drug in question is considered effective against all indications in this list\n",
    " \n",
    "Other variables may also be associated with the Compound object. Consult [the CANDO documentation](https://github.com/ram-compbio/CANDO/blob/master/docs/CANDO-v2.2.pdf) for more information.\n",
    "\n",
    "To explore ciprofloxacin, we will look into some of the properties listed above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(cipro.name, cipro.id_)\n",
    "print(cipro.name, 'is', cipro.status)\n",
    "print('It is associated with', len(cipro.indications), 'indications')\n",
    "print('Its signature contains', len(cipro.sig), 'proteins')\n",
    "print('Signature:', cipro.sig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ciprofloxacin has existed since the 1980s, and it has accumulated a large number of indication associations in that time. Its signature represents its interactions with the 64 proteins we included in our matrix. These scores are from 0.0 to 1.0; we can see that ciprofloxacin does not have a high predicted binding affinity with any of these proteins."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's look at Indication objects. CANDO uses [**Me**dical **S**ubject **H**eadings (MeSH)](https://meshb.nlm.nih.gov/) IDs to internally catalogue indications. In order to predict new drugs for a given indication, you must first choose the most appropriate MeSH term. While you could search MeSH itself for the right term, not all terms have an Indication associated with them in CANDO. Searching CANDO for a MeSH ID ensures that your term is already present in your model.\n",
    "\n",
    "Finding an Indication is similar to finding a Compound. You can start by using `cando.search_indication` to find the closest text match to the disease you are interested in. If we wanted to find pneumonia, we could start with the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.search_indication('Pneumonia')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many diseases have multiple comorbidities, etiologies, or otherwise related indications, so it is very common to get multiple results for an indication search. In this case, it seems like viral, aspiration, and bacterial pneumonia should have different causes and possibly different treatments, so it makes sense to focus on a specific type rather than the general term. If we want to study bacterial pneumonia, we would look at \"Pneumonia, Bacterial\" with a MeSH ID of [MESH:D018410](https://meshb.nlm.nih.gov/record/ui?ui=D018410). \n",
    "\n",
    "We can use this identifier to access the Indication itself using `get_indication`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pneu_bact = cando.get_indication(\"MESH:D018410\")\n",
    "print('The variable pneu_bact contains the indication:', pneu_bact.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like a Compound, an Indication contains various information about the indication as it exists in your CANDO object. This includes:\n",
    "- **indication.id_** - string, the MeSH identifier of the indication (i.e. MESH:D018410)\n",
    "- **indication.name** - string, the English name of the indication (i.e. Pneumonia, Bacterial)\n",
    "- **indication.compounds** - list of integer, the CANDO ID of all compounds associated with the indication\n",
    "  - If a compound is listed here, it is considered effective against our indication\n",
    "\n",
    "Other properties may also be associated with the Indication object. Consult [the CANDO documentation](https://github.com/ram-compbio/CANDO/blob/master/docs/CANDO-v2.2.pdf) for more information.\n",
    "\n",
    "We can find out more about how many compounds known to treat bacterial pneumonia using these properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pneu_bact.name, pneu_bact.id_)\n",
    "print(len(pneu_bact.compounds), 'compounds are associated with bacterial pneumonia')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By combining Compound and Indication properties, we can get even more information about what drugs are associated with bacterial pneumonia in our CANDO object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The following compounds are associated with bacterial pneumonia:')\n",
    "print(*[cando.get_compound(cmpd).name for cmpd in pneu_bact.compounds], sep=', ')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you check the list above, you will see a familiar compound. Bacterial pneumonia is one of the 77 indications associated with ciprofloxacin.\n",
    "\n",
    "Now that we can access drug and indication information, we can continue on to completing practical predictions and assessments in the [Drug prediction](#Drug-prediction) and [Benchmarking](#Benchmarking) sections.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Drug prediction\n",
    "\n",
    "One of the most common ways to use CANDO is to predict new drugs for an indication of interest. Let's say we want to predict new drugs for human immunodeficiency virus (HIV). We can start by searching for \"HIV\" to find the relevant MeSH term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.search_indication('HIV')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HIV is associated with the development of multiple other conditions, so multiple results are returned. However, if we are specifically looking for drugs effective against HIV as a *virus*, not treatments for the resulting conditions, \"HIV Infections\" is the result we are looking for. This is associated with [MESH:D015658](https://meshb.nlm.nih.gov/record/ui?ui=D015658).\n",
    "\n",
    "Because CANDO is a similarity-based drug discovery method, it is generally more effective at predicting new drugs for indications that already have multiple drugs associated with them. In addition, if there are no drugs associated with the indication we want to predict effective drugs for, we would have to use a different approach to predict novel drugs (explored in the [Advanced topics](#Advanced-topics)).\n",
    "\n",
    "To ensure we have enough drugs to create effective predictions against HIV, we should find the Indication associated with HIV infections and check how many drugs are already associated with it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hiv = cando.get_indication(\"MESH:D015658\")\n",
    "\n",
    "print(len(hiv.compounds), 'compounds are associated with', hiv.name)\n",
    "print(*[cando.get_compound(cmpd).name for cmpd in hiv.compounds], sep=', ')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 32 drugs associated with HIV Infections in our data, which makes sense as it is a well-studied disease. Most seem to be protease inhibitors (saquinavir, fosamprenavir) or nucleoside analogs (zalcitabine, didanosine). \n",
    "\n",
    "We can now make predictions for other drugs/compounds that may be useful against HIV using the `canpredict_compounds` function. Given an indication, canpredict_compounds takes the lists of most similar compounds to each drug already associated with that indication. It then checks to see which compounds are similar to multiple indicated drugs; compounds that are similar to more indicated drugs are ranked higher, and ties are broken based on average rank within the similar lists. Possible parameters include:\n",
    "\n",
    "- **ind_id** - string, the indication for which new compounds should be predicted (i.e. MESH:D015658)\n",
    "- **n** - integer, the number of similar compounds to be considered per indicated compoud (10 by default)\n",
    "- **topX** - integer, the maximum number of predictions to be outputted (10 by default)\n",
    "- **consensus** - boolean, outputs only compounds that appear in multiple similar lists if True (True by default)\n",
    "- **keep_associated** - boolean, includes indicated drugs in the output list if True (False by default)\n",
    "- **cmpd_set** - \"all\" \"approved\" or \"other\", determines which compound set to use (\"all\" by default)\n",
    "  - Note, if the CANDO object was created with only approved compounds, \"all\" will still only include aproved compounds\n",
    "- **save** - string, name of an output file to save results (tsv format, nothing/not saved by default)\n",
    "\n",
    "Because we have 32 drugs associated with HIV, the full output of canpredict_compounds would be very long. Therefore, we will only look at the top 10 results (`topX = 10`). There are plenty of opportunities for a compound to appear in at least 2 of the 32 similar lists, so looking at only the top 10 most similar compounds (`n = 10`) and excluding any compounds that appear in only one similar list (`consensus = True`) are also appropriate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.canpredict_compounds(\"MESH:D015658\", n=10, topX=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These results are a little complex at first glance, so let's break down some key columns:\n",
    "\n",
    "- **score1** - integer, the number of similar lists in which the compound appears at or above rank n (here, rank 10)\n",
    "  - This is the primary ranking factor\n",
    "- **score2** - float, the average rank of each compound in the \"score1\" similar lists it appears in\n",
    "  - When score1 is tied for two compounds, this breaks the tie\n",
    "- **id** - integer, the CANDO ID of the compound; can be used to directly access the compound using `cando.get_compound()`\n",
    "- **approved** - boolean, true if the compound is listed as approved in the input data; false otherwise\n",
    "- **name** - string, the name of the compound\n",
    "\n",
    "The results can change greatly if we alter the input parameters. Instead of considering only the top 10 most similar compounds to each indicated compounds, let's change it to the top 25 (`n = 25`) and see what changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.canpredict_compounds(\"MESH:D015658\", n=25, topX=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the first compound is the same (darunavir), the second and third, paclitaxel and decitabine, are new; previously, these two did not even appear in the top 10 predictions. Meanwhile, the previous third place compound, entecavir, is not in the top 10 anymore. To find out where it ended up, we can see more results for the same assessment using `topX = 25`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.canpredict_compounds(\"MESH:D015658\", n=25, topX=25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns out entecavir only appeared in one more similar list when we considered the top 25 instead of the top 10 most similar compounds, pushing it down to rank 12.\n",
    "\n",
    "If we want to, we can also see where compounds already associated with HIV appear relative to newly predicted compounds using the same assessment using `keep_associated = True`. Technically, compounds associated with HIV are at a slight disadvantage, since they cannot appear in their own similar lists, but, because they are effective against HIV, we would still expect them to be pretty highly ranked overall if CANDO is functioning well. Let's run this assessment with `n = 10` and `topX = 10`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.canpredict_compounds(\"MESH:D015658\", n=10, topX=10, \n",
    "                           keep_associated=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If an asterisk * appears next to the label in the \"approved\" column, that means that drug is already associated with HIV. We can see there are three such compounds that appear in this list: amprenavir, saquinavir, and nelfinavir. As it turns out, CANDO ranks darunavir higher than every compound approved to treat HIV in this prediction.\n",
    "\n",
    "Before moving on, let's repeat this assessment on an indication with many fewer approved drugs. When we searched \"HIV\" earlier, one of our results was \"HIV Wasting Syndrome,\" MESH:D019247. Let's look at that indication now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hiv_ws = cando.get_indication('MESH:D019247')\n",
    "print(hiv_ws.name, hiv_ws.id_)\n",
    "print(len(hiv_ws.compounds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are only 4 drugs associated with HIV wasting syndrome in our dataset. This will likely reduce the number of results available through `canpredict_compounds`. Let's try running it with the same parameters as we did previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.canpredict_compounds('MESH:D019247', n=10, topX=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, though we still requested the top 10 predictions, only five appeared because only five compounds appear in multiple similar lists of the four drugs associated with HIV wasting syndrome. Stanozolol, the top ranked prediction, appears in only two similar lists (score1 is 2), but it appears in rank 2 on average (score2 is 1.0; rank 1 is 0.0).\n",
    "\n",
    "If we want to increase the number of predictions we receive, we can try a couple things. First, as we have done previously, we can increase `n` to 25, which is likely to increase the number of compounds appearing in multiple similar lists. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.canpredict_compounds('MESH:D019247', n=25, topX=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a full 10 results, but the ranks of most of our previous results (aside from testosterone propionate) have fallen. If we want to extend the previous list instead of changing the results, we can keep `n = 10` but set `consensus = False`. This will include in the ranks compounds that only appear in one similar list, but are highly ranked in that list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.canpredict_compounds('MESH:D019247', n=10, topX=10, consensus=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, this output has the same top ranks as our first results, but the table has been extended to include additional compounds.\n",
    "\n",
    "While CANDO is more often used to predict novel drugs for an indication, it can also be used for the reverse, predicting new indications for an existing compound. This could be useful, for example, if one discovers a new small molecule in nature and wonders if it could have any pharmaceutical applications. `canpredict_indications` finds the most similar drugs to the compound in question, determines what those drugs are indicated to treat, and then returns any indications associated with multiple most similar drugs as potential uses. \n",
    "\n",
    "Key parameters of `canpredict_indications` are largely the same as those of `canpredict_compounds`, with three exceptions: the first argument is the Compound object, as opposed to the indication id; there is no `keep_associated` parameter; and there is an additional `sorting` variable that takes a string to determine whether to rank the outputtted indications by probability (`\"prob\"`) or score (`\"score\"`).\n",
    "\n",
    "Let's look at the antimicrobial paromomycin as an example. First, we have to find its id:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.search_compound('paromomycin')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we get the Compound object and check which indications paromomycin is already associated with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "paro = cando.get_compound(1225)\n",
    "for indic in paro.indications:\n",
    "    print(cando.get_indication(indic).name, cando.get_indication(indic).id_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like paromomycin is associated with six infection-related indications. We will use the top 10 most similar compounds to paromomycin (`n = 10`) and print the top 10 results (`topX = 10`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cando.canpredict_indications(paro, n=10, topX=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of paromomycin's existing indications are infections, so it is unsurprising that all of the predictions are also infections. As with the output of `canpredict_compounds`, the score indicates how many times an indication appears associated with the top 10 (or top n) most similar compounds to paromomycin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section covered the most basic ways to predict a new drug for an indication and an new indication for a drug. To learn about additional drug discovery methods within CANDO, you may continue to the [Advanced topics](#Advanced-topics) section or consult [the CANDO documentation](https://github.com/ram-compbio/CANDO/blob/master/docs/CANDO-v2.2.pdf).\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Benchmarking\n",
    "\n",
    "Benchmarking is essential for assessing how well CANDO is working, particularly when using a new dataset or prediction function, but also when looking at predictions for a specific indication in greater depth. Benchmarking ensures that our platform has predictive value and that we know how much to trust the predictions we are working with.\n",
    "\n",
    "CANDO has a primary built-in benchmarking protocol, as well as multiple older benchmarking functions. We will focus on the primary benchmarking protocol, `canbenchmark_new`, here. This method uses a leave-one-out benchmarking protocol on every indication with two or more associated drugs. Every drug is withheld from its indication, and CANDO predicts novel drugs for that indication based on the remaining drugs. If the left out drug appears in these predictions, that is considered a success, as CANDO is able to predict the efficacy of that drug.\n",
    "\n",
    "The primary benchmarking function calculates four metrics:\n",
    "\n",
    "- **New average indication accuracy (nAIA)** - float, average % of indicated drugs that can be predicted within a given cutoff\n",
    "  - The \"new\" distinguishes these metrics from old average indication accuracy, as calculated by our older benchmarking functions\n",
    "  - A **Control nAIA** is also provided; this is the nAIA we would anticipate if CANDO worked as well as random chance\n",
    "- **Pairwise accuracy (PA)** - float, % of drugs that appear within a given cutoff in the similarity lists of drugs with the same indication\n",
    "  - PA assesses similarity list quality, whereas nAIA assesses prediction quality\n",
    "- **Indication coverage (IC)** - integer, the number of indications with a non-zero indication accuracy\n",
    "- **Normalized discounted cumulative gain (nNDCG)** - float, scores # indicated drugs predicted within a given cutoff with prioritization of early recall\n",
    "\n",
    "Higher scores in these metrics should lead to higher confidence in our predictions. Additional metrics, such as precision at a given rank or area under the receiver operating characteristic (AUROC/AUC), may also be calculated. See the [the CANDO documentation](https://github.com/ram-compbio/CANDO/blob/master/docs/CANDO-v2.2.pdf) or [this paper on evaluating drug repurposing technologies](https://www.biorxiv.org/content/10.1101/2020.12.03.410274v1) for more information about metrics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`canbenchmark_new` has a couple key arguments:\n",
    "- **file_name** - string, a descriptive name to be assigned to all benchmarking files created\n",
    "- **n** - integer, the number of similar compounds to be considered per indicated compoud (10 by default)\n",
    "  - Use the same value for `n` as you used in `canpredict_compounds` to ensure your benchmarking results are applicable\n",
    "- **indications** - list of string, a list of the subset of indications to be examined\n",
    "  - If this argument is not used, all indications with at least two associated drugs will be assessed\n",
    "- **associated** - boolean, if True (default), excludes all drugs that are not associated with an indication from consideration\n",
    " \n",
    "As with other CANDO functions, additional arguments are available. `canbenchmark_new` has not yet been added to the documentation, but some arguments are explained in the function header within cando.py.\n",
    "\n",
    "We can choose to assess CANDO in two ways: we can assess the overall performance of CANDO, or we can assess a subset of indications of interest. The former is especially useful when we are using an unusual dataset that might affect the performance of CANDO. In this tutorial, we are using a matrix with only 64 proteins, which might alter performance, so it makes sense to run an overall assessment. We will give this assessment the name \"tutorial_assessment\" and use n=10, as we did when running `canpredict_compounds`. **Note that this function will take a long time to run.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = cando.canbenchmark_new('tutorial_assessment', n=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The primary results of `canbenchmark_new` are printed out directly for you to read. Based on the nAIA results, across all indications, CANDO is able to recall 8.6% of indicated drugs in the top 10 compounds (out of 1464), 13.3% in the top 25, and 23.8% in the top 100. CANDO also outperforms random chance (control-nAIA) at every cutoff. In addition, the indication coverage (IC) tells us that 523 indications (out of 1595 assessed; see top100% column) have non-zero performance at the top 10 threshold, 721 at the top 25 threshold, and 982 at the top 100 threshold. Note that these results correspond to the 64-protein matrix used in this tutorial (plus the compound and indication mapping); as a general rule, performance is better when a larger protein set is considered.\n",
    "\n",
    "In addition, multiple additional files that may be of interest are generated:\n",
    "\n",
    "- **Summary file** - TSV, contains overall metrics; top-level folder with a name starting \"summary\"\n",
    "- **Raw results** - CSV, contains the rank at which each indicated drug was predicted; raw_results folder\n",
    "- **Pairwise results** - CSV, contains the ranks of each pair of drugs with the same indications in each others' similarity lists; pairwise_results folder\n",
    "- **Results analysed** - TSV, contains indication accuracy or NDCG scores for each individual indication; results_analysed_named folder\n",
    "\n",
    "These files can be viewed in the tutorial folder, and they can be opened in a text editor like Notepad or in certain spreadsheet editors, like Excel. For now, we will just look at the results of the summary file, which should reflet the results printed out when running `canbenchmark_new` above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('summary-tutorial_assessment-10-associated.tsv', 'r') as f:\n",
    "    print(f.read()) # Note: do NOT do this with the longer, non-summary files; it may crash or hang"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The other reason one might want to complete benchmarking is to determine how well CANDO performs on a given indication. For example, earlier we were looking into HIV infections, MESH:D015658. Let's replicate those results here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.canpredict_compounds(\"MESH:D015658\", n=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have the same 10 results as before, but now we want to estimate how trustworthy these results are. We can use `canbenchmark_new` with the same `n` value as we used in our prediction to get such an estimate. Since we *only* want to assess how well CANDO performs on HIV here, we can pass a list with only the MeSH ID we are interested in to the `indications` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = cando.canbenchmark_new(\"HIV_assessment\", n=10, indications=[\"MESH:D015658\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We found that 9.375% of drugs associated with HIV infections, or 3 of the 32 drugs, were ranked within the top 10 compounds predicted for this indication, as compared to 2449 compounds ranked overall. If you recall, this is also what we found when we ran `canpredict_compounds` with `keep_associated=True`. When we look at the top 25 compounds, this increases to 21.875%, or 7 out of 32 drugs.\n",
    "\n",
    "Based on these results, we can expect that not all compounds effective against HIV infections are going to be in our top 10 predictions, which is unsurprising. However, by comparing our nAIA to the control nAIA, we can see that CANDO is performing far better than random chance on predicting effective drugs for HIV. CANDO also performs better on HIV infections than on indications in general. Therefore, it is reasonable to use CANDO to predict novel drugs for HIV infections."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Historically, we have assessed the performance of CANDO based on calculating accuracies for the similar lists of every compound within an indication. This is the idea behind the original CANDO benchmarking algorithm, `canbenchmark`. Although the newer `canbenchmark_new` is more relevant to most applications, the original algorithm is still preserved for posterity and potential in the development of CANDO.\n",
    "\n",
    "As an example of the original `canbenchmark`, let's look at hyperparathyroidism, MESH:D006961. There are four drugs associated with this indication: cinacalcet (884), paricalcitol (786), dihydrotachysterol (938), and calcitriol (32). As it turns out, cinacalcet does not have any of the other three drugs in its top 10 most similar compounds. However, paricalcitol, dihydrotachysterol, and calcitriol are all in each other's top 10 most similar compounds lists. Since three out of the four compounds' similar lists have at least one other indicated compound in the top 10, we calculate an indication accuracy of 75% at the top 10 threshold for hyperparathyroidism. This can be repeated for multiple rank thresholds (top 25, 100, etc) and for every other indication with at least two associated compounds. From this, the unweighted AIA, weighted APA, and IC can be calculated for these thresholds. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Doing this work by hand is tedious and repetitive, which is why we created `canbenchmark` to rapidly calculate these metrics for every indication. `canbenchmark` only requires one argument, a string that determines what the output files it creates will be named. The protein (or other interaction) matrix, compound list, and compond-indication associations given to the CANDO object as input may all affect the bnechmarking performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.canbenchmark('tutorial')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Besides the abbreviated results it prints out, `canbenchmark` creates three additional files: a raw results file that contains information on every compound-indication pair (found at raw_results/raw_results-tutorial.csv), an analyzed results file that contains every individual indication accuracy (found at results_analysed_named/results_analysed_named-tutorial.csv), and a summary file containing overall metrics (found at summary-tutorial.tsv). You can open these results from your file navigator into spreadsheet software like Excel or text editors like Notepad, and these files can also be directly opened and examined via your code, as below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('summary-tutorial.tsv', 'r') as f:\n",
    "    print(f.read()) # Note: do NOT do this with the longer, non-summary files; it may crash or hang"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, we have an AIA of 20% at the top 10 threshold, meaning that, across all indications, about 20% of compounds had at least one other compound with the same indication in their top 10 most similar lists. The APA is higher because the chance of recovering an indicated compound in the top ranks is higher when there are more indicated compounds, so more heavily weighting larger indications leads to a better result. Indication coverage is 819, meaning that 819 out of 1595 indications with at least two associated compounds had an indication accuracy greater than 0%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced topics\n",
    "\n",
    "This section will cover more advanced and customizable ways of working with CANDO. Note that this section will not go into the same depth as the basic tutorial, as it is assumed you already understand the basics of CANDO."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Customizing dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interaction matrix\n",
    "\n",
    "The example interaction matrix, matrix_file, has already been downloaded via get_tutorial(). This function may take anywhere from ~1-10 mins on 3 cores, depending on your computer's processor.\n",
    "\n",
    "In this step we will generate a matrix of 2,449 approved compounds by 64 proteins, populated with the corresponding interaction score betwen each drug and protein. The final matrix will have drugs as the columns (indexed according to the compound mapping file), and proteins as the indices (indexed by PDB and chain ID).\n",
    "\n",
    "The function generate_matrix() first creates a dataframe for all chemical fingerprint similarity scores comparing each compound in the specified version library v (or only approved drugs with approved_only=True) to every potential binding site ligand from the PDB using RDKit to compute the chemical fingerprints in which fp denotes the type/radius of fingerprint (rd_ecfp4, rd_ecfp8, etc). The vector type, vect, can be binary (\"1024_bit\" or \"2048_bit\") or integer (\"int\") and denote the presence or absence or count of molecular substructures in the molecule, respectively. If using integer vectors, dist should be set to \"dice\", but binary vectors can be set to Tanimoto (\"tani\"). Next, it creates a dataframe of all potential binding sites for each protein in the specified protein library, protlib, with their corresponding binding site scores from the specified binding site prediction method, bs (\"coach\", \"cof\", \"ssite\", \"tms\"). The function then iterates over all drugs and protein binding sites to populate the matrix with the best score based on the following input parameters: 1) i_score - the scoring protocol of choice, which can be 'C', 'dC', 'P', 'CxP', 'dCxP', where C is the fingerprint similarity score (either Tanimoto or Sorenson-Dice coefficient), dC is the percentile of the C score for the compound compared to all ligands in the library, and P is binding site score associated with the ligand predicted by the COACH algorithm, 2) c_cutoff and p_cutoff - set cutoffs which ignore any C or P scores below each threshold, respectively, and 3) percentile_cutoff - similar to c_cutoff but with the dC score (overrides c_cutoff if not None). These scores serve as a proxy for binding strength/probability of the drug and protein target. We then output the matrix to a tsv file, out_file, which in this case is \"tutorial_matrix-all.tsv\" (out_path can be set to write the file to a specific directory).\n",
    "\n",
    "We will create a matrix using the 'CxP' protocol, which just chooses the top Dice score to any ligand predicted to bind to a given protein, without any cutoffs. This function is parallelized, so setting the variable ncpus will change the number of processors that are used for this function. NOTE: percentile cutoff protocols ('dC', 'dCxP') take much longer to compute than 'C' and 'P' protocols. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate example cando interaction matrix (2,449 compounds x 64 proteins)\n",
    "cnd.generate_matrix(v=\"test.0\", fp=\"rd_ecfp4\", vect=\"int\",\n",
    "                    dist=\"dice\", org=\"tutorial\", bs=\"coach\",\n",
    "                    c_cutoff=0.0, p_cutoff=0.0, percentile_cutoff=0.0,\n",
    "                    i_score=\"dCxP\", out_file='', out_path=\".\",\n",
    "                    nr_ligs=True, approved_only=True,\n",
    "                    lib_path='', prot_path='', lig_name=False, ncpus=ncpus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Novel compound\n",
    "\n",
    "The CANDO platform contains an extensive library of approved drugs and other compounds from DrugBank. However, if you wish to predict indications or similar drugs for a compound that is not present in our library, we make it possible with the `generate_signature()` function.\n",
    "\n",
    "First, you must have the compound properly formatted in mol file format. There are many programs that provide conversion between many chemical file formats, such as OpenBabel.\n",
    "\n",
    "Next, run `generate_signature()`. This will populate a tsv file with Tanimoto/Sorenson-Dice similarity scores of the provided compound to all binding site ligands in our database. These values will be used for the generation of the drug-proteome signature. The input parameters for this function are very similar to those for `generate_matrix()`, however the first argument is the path to the compound structure file in mol format. The output signature file will be saved in tsv format with the name you provide for \"out_file\" (with the appended path from \"out_path\").\n",
    "\n",
    "**NOTE: the interaction score protocol must match the input matrix protocol, which was 'CxP' as above. If we used a different protocol, these scores would not be directly comparable.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmpd_file = \"lmk235.mol\"\n",
    "signature_file = \"lmk235_signature.tsv\"\n",
    "\n",
    "cnd.generate_signature(cmpd_file, fp=\"rd_ecfp4\", vect=\"int\", dist=\"dice\", \n",
    "                      org=\"tutorial\", bs=\"coach\", c_cutoff=0.0, p_cutoff=0.0, \n",
    "                      percentile_cutoff=0.0, i_score=\"CxP\", out_file=signature_file, \n",
    "                      out_path=\".\", nr_ligs=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then must add the compound to the existing `CANDO` object with the `add_cmpd()` function. The inputs are the newly generated signature file and the desired name (which below is \"lmk-235\"). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.add_cmpd(signature_file, new_name='lmk-235')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that our new compound is added to the platform, we can see what other compounds to which it is similar. We can use the `similar_compounds()` function from before to print those compounds. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lmk235 = cando.get_compound(10674)\n",
    "cando.similar_compounds(lmk235, n=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can predict potential indications for which our new compound may be useful. We can use the `canpredict_indications()` as before to print those results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.canpredict_indications(lmk235, n=25, topX=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interestingly, both Hypertension (MESH:D006973) and Pain (MESH:D010146) are top predictions - the analgesic and hypotensive properties of lmk-235 are both supported by in vivo studies in the literature. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compound library\n",
    "\n",
    "In order to create a new compound set, you must first have a TSV (tab separated values) file that contains one of two chemical file types:\n",
    "\n",
    "1. SMILES - `file_type='smi'`. The file must have the SMILES string as the first column and the corresponding compound name in the second column, e.g.\n",
    "    `C1CNCCN(C1)S(=O)(=O)C2=CC=CC3=C2C=CN=C3 fasudil`\n",
    "    \n",
    "2. Mol - `file_type='mol'`. The file must have the name of the file, without the file extension. In addition, the path to the files must be given in the argument `cmpd_dir`.\n",
    "\n",
    "In our example, we will create a new compound library of select tyrosine kinase inhibitors (TKIs) that we will then use to create a new matrix with just these drugs. First, we generate the library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnd.add_cmpds(\"tki_set-test.smi\", file_type='smi', fp=\"rd_ecfp4\", vect=\"int\", cmpd_dir=\".\", v='tki')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now generated a new compound library, located at `./tki/cmpds`. This location is partially hardcoded, so it will also create a directory in your current working directory, `./`, named after the `v` you choose, `./tki`. This makes it easier for the end-user.\n",
    "\n",
    "Notice we chose to use the rdkit fingerprint ecfp4, `fp=\"rd_ecfp4\"`, and the vector type int, `vect=\"int\"`. We need to keep track of this for our matrix generation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we will generate a matrix using this **tki** compound library. The generate matrix function can take in a customized matrix, as opposed to the pregenerated/downlaoded matrices based upon our predefined `v`, e.g. v2.2, v2.3, etc. By setting `v` to the same name you used in the previous `add_cmpds` function, you can load those cmpds to create a new matrix.\n",
    "\n",
    "We will set `v='tki'` and `lib_path='.'`. This means we will use the **tki** library, which is located in the current working dir. In addition, sicne we created our compound library using fingerprint type ecfp4, we need to make sure we use it here, as well. The remaining arguemnts has already been discussed in the previous [Interaction matrix](#Interaction-matrix) section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnd.generate_matrix(v=\"tki\", lib_path='.',\n",
    "                    fp=\"rd_ecfp4\", vect=\"int\",\n",
    "                    dist=\"dice\", org=\"tutorial\", bs=\"coach\",\n",
    "                    c_cutoff=0.0, p_cutoff=0.0, percentile_cutoff=0.0,\n",
    "                    i_score=\"CxP\", out_file='', out_path=\"\",\n",
    "                    nr_ligs=True, approved_only=False,\n",
    "                    lig_name=False, ncpus=ncpus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a brand new matrix with all of our new TKIs scored against our set of 64 test proteins. This file is located at `./tki/matrices/`. Again, these paths are mostly hardcoded. You can control where the matrix is written and what the name is using the `out_path` and `out_file` arguments. When these arugments are set to '' the hardcoded paths are used, as shown in this example.\n",
    "\n",
    "Now we will load the newly created matrix into a CANDO object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set compound mapping variable to the new TKI compound mapping file\n",
    "tki_map = 'tki/mappings/cmpds-tki.tsv'\n",
    "tki_matrix = 'tki/matrices/rd_ecfp4-int-dice-tutorial-coach-c0.0-p0.0-CxP.tsv'\n",
    "# Create CANDO object using the new compound mapping and TKI matrix\n",
    "tki_cando = cnd.CANDO(tki_map, ind_map, matrix=tki_matrix, compound_set='all', compute_distance=True, \n",
    "                  dist_metric=dist_metric, ncpus=ncpus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the new CANDO object data and the most similar compounds to the first TKI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print cando object stats\n",
    "print('compounds', len(tki_cando.compounds))\n",
    "print('indications', len(tki_cando.indications))\n",
    "print('proteins', len(tki_cando.proteins))\n",
    "print('')\n",
    "\n",
    "# print first TKI name and signature\n",
    "c = tki_cando.compounds[0]\n",
    "print(c.name, len(c.sig))\n",
    "print(c.sig)\n",
    "print('')\n",
    "\n",
    "# top5 most similar compounds to first TKI\n",
    "for s in c.similar[0:5]:\n",
    "    print(s[0].name, round(s[1], 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced drug prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### De novo prediction\n",
    "\n",
    "Sometimes there are no drugs associated with an indication (whether it be in reality or in your input indication mapping); in those cases we can use the `canpredict_denovo()` function to suggest putative candidates. Basically, this function counts/sums the number of protein interaction scores above a set threshold for each compound and ranks them according to frequency and strength. This is particularly useful for pathogen proteomes (like SARS-CoV-2 or bacterial proteomes) or for finding which compounds target a subset of proteins of interest (e.g. kinases). In our case, we have a diverse set of proteins, which would render the results meaningless. Instead, we can input a list of protein IDs in which we are interested. We have precompiled a list of bacterial proteins already present in the sample matrix for convenience. Note: if we had read in an indication-genes mapping file, we could use the `ind_id=` parameter to automatically select all proteins associated with said indication. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bacterial_proteins = [\"3mk7C\", \"3eziA\", \"1u2mC\", \"3atsA\", \"2x4mA\", \n",
    "                      \"4wliA\", \"1t4aA\", \"4zxkA\", \"1zhhA\", \"1eb0A\", \n",
    "                      \"4nqwB\", \"2gqrA\", \"2qlcA\", \"3rf1A\", \"2xpwA\"]\n",
    "\n",
    "cando.canpredict_denovo(method='sum', threshold=0.6, topX=30, proteins=bacterial_proteins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Machine learning\n",
    "\n",
    "The \"proteomic vectors\" within CANDO lend themselves well to machine learning to perhaps learn more complex relationships between the proteins within the vector and their impacts on the treatment of diseases. CANDO has built-in ML algorithms that allow for two main functionalities:\n",
    "1. Benchmark the platform using a hold-one-out protocol very similar to canbenchmark\n",
    "2. Make predictions for novel or non-associated compounds that may be therapeutic for a given disease\n",
    "\n",
    "The ML module currently supports 2 algorithms: random forests and logistic regression. The models are trained on drugs approved for the disease (positive classes) and an equal number of randomly selected \"neutral samples\", which are drugs/compounds not approved for the disease (negative samples). Random seeds may be set to ensure the same compounds are used in training. \n",
    "\n",
    "We have the option to benchmark the platform with an ML algorithm - this module outputs files very similar to canbenchmark. For this tutorial, we will skip this function as it requires a great deal of time to complete (training a separate model for EVERY drug-disease association, basically). The command to do so with a logistic regression classifier is below, feel free to run it! The `'out='` flag defines the name of the output files. Again, only diseases with 2+ compounds associated are benchmarked. \n",
    "\n",
    "`cnd.ml(method='rf', benchmark=True, seed=50, out='test_rf')`\n",
    "\n",
    "We can also use this module to predict if a certain compound may be therapeutic for a given disease. We can use the \n",
    "`'predict='` flag to specify a list of compounds that we wish to predict with the classifier. Let's use three drugs, imatinib, buprenorphine, and lisdexamfetamine, and see if they are predicted to have antibiotic activity using a random forest classifier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "baci = cando.get_indication('MESH:D001424')\n",
    "\n",
    "imat = cando.get_compound(503)\n",
    "bup = cando.get_compound(797)\n",
    "lamf = cando.get_compound(1121)\n",
    "\n",
    "cando.ml(method='rf', effect=baci, benchmark=False, seed=50, predict=[imat, bup, lamf])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The accuracy is quite poor, so fine-tuning and hyperparameter optimization is necessary to enhance performance. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Custom protein sets\n",
    "\n",
    "It may be useful for some users to probe compound-protein interaction similarity, but only in the context of a few particular proteins (e.g. set of kinase inhibitors). Instead of generating a new matrix with all of these proteins and their corresponding interaction values, which can begin to take up a lot of storage if done multiple times, the `'protein_set='` flag can be specified during the instantiation of the CANDO object. This flag contains the path to the protein subset the user wishes to use, which is simply a list of UniProt protein IDs. The CANDO object will automatically check for each ID if it either simply matches any UniProt IDs within the matrix or if that UniProt ID is associated with any PDB chains within the matrix (based on a mapping from the SIFTs project). If there are matches, the CANDO object will now contain Compound objects with only those protein interaction values in their signatures. Below is how to use this functionality - this time we can search for the bacterial proteins from the previous `canpredict_denovo()` section, but filter the proteins from the start using a list of UniProt IDs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando_subset = cnd.CANDO(cmpd_map, ind_map, matrix=matrix_file, compound_set='approved',\n",
    "                         compute_distance=True, protein_set=protein_set,\n",
    "                         dist_metric=dist_metric, ncpus=ncpus)\n",
    "\n",
    "print(\"Number of proteins in new signature =\", len(cando_subset.proteins))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The signature was successfully edited to 15 proteins. Note: this does not nececessarily mean the each UniProt ID had a corresponding PDB match -- multiple PDB chains can be associated with a given UniProt ID. If the matrix contains UniProt IDs in the first column instead of PDB IDs, as in this case, the \"Direct UniProt matches\" value would be incremented if they match.\n",
    "\n",
    "We can also repeat all benchmarks and predictive algorithms with the new signatures. Below is the default benchmarking results with the new signatures. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando_subset.canbenchmark('test_subset')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's repeat the ML code from above, but this time with the reduced protein subset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "baci = cando_subset.get_indication('MESH:D001424')\n",
    "\n",
    "imat = cando_subset.get_compound(503)\n",
    "bup = cando_subset.get_compound(797)\n",
    "lamf = cando_subset.get_compound(1121)\n",
    "\n",
    "cando_subset.ml(method='rf', effect=baci, benchmark=False, \n",
    "                seed=50, predict=[imat, bup, lamf])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the output probability significantly change for one of the compounds, lisdexamfetamine, from 0.340 to 0.600, illustrating the impact of protein subset composition on model behavior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Virtual screening\n",
    "\n",
    "Though the CANDO platform is mainly intended for multitarget drug discovery and repurposing, there still exists the option to check the top compound hits for a given protein. This is accomplished via the `virtual_screen()` function. Let's test the top hits for two of the known bacterial proteins from above, namely \"1u2mC\", \"3atsA\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cando.virtual_screen(\"1u2mC\")\n",
    "cando.virtual_screen(\"3atsA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Streptozocin is the top hit for 1u2mC, which suggests it shares significant structural similarity to a ligand known/predicted to bind to this protein. Similarly, kanamycin is the 5th hit for 3atsA. These scores/ranks can change significantly based on scoring protocol used for the matrix - changing i_score in `generate_matrix()` to \"dC\", for example, would greatly affect the output. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Credits: CANDO tutorial created by Zackary Falls; revised by Melissa Van Norden*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}