data_process.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "# Python modules used for API access...\n",
    "from chembl_webresource_client.new_client import new_client\n",
    "from tqdm import tqdm\n",
    "import pickle\n",
    "import numpy as np\n",
    "from functions import PROCESS\n",
    "import matplotlib.pyplot as plt\n",
    "import collections\n",
    "\n",
    "# Initialize PROCESS object\n",
    "drugs = pd.read_csv('data/primary_data/KIBA_drugs.csv')\n",
    "obj = PROCESS(drugs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Module I: Data retrieval and processing\n",
    "Use this block to retrieve/process your data by providing a drug file\n",
    "\"\"\"\n",
    "# Step 1: Extaract terget data for 52k drugs from CHEMBL\n",
    "main_dict = obj.extractDrugTargets(drugs['drugs'].tolist())\n",
    "obj.save_obj(main_dict, 'data/main_dict.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2: Extract Drug-Target pairs for Homo-Sapiens with IC50 values\n",
    "main_dict = obj.load_obj('data/drug_target.pkl')\n",
    "DTI = obj.getDTI(main_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": "100%|██████████| 961/961 [37:54<00:00,  2.37s/it]\n"
    }
   ],
   "source": [
    "# Step 3: Get mapping from CHEMBL to uniprot and save at data/out.txt\n",
    "obj.chembl2uniprot('data/DTI_index.csv', 'data/chembl2uniprot961.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Step 4: Screen Data based on units, sequence-availiability\n",
    "DTI = pd.read_csv('data/primary_data/DTI2.csv')\n",
    "#clean_DTI = obj.screenDTI(DTI, 'data/chembl2uniprot.txt', 'data/clean_DTI.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": "100%|██████████| 76872/76872 [00:01<00:00, 45112.92it/s]\n"
    }
   ],
   "source": [
    "# Step 5: Transform all activities into one unit (uM)\n",
    "final_clean_DTI = pd.read_csv('data/clean_DTI.csv')[['target','drug','IC50','unit']]\n",
    "for i in tqdm(range(len(final_clean_DTI))):\n",
    "    if final_clean_DTI.at[i,'unit']=='pM':\n",
    "        final_clean_DTI.at[i,'IC50']=final_clean_DTI.at[i,'IC50']*1000000   # pM to uM\n",
    "        final_clean_DTI.at[i,'unit']='uM'\n",
    "    elif final_clean_DTI.at[i,'unit']=='nM':\n",
    "        final_clean_DTI.at[i,'IC50']=final_clean_DTI.at[i,'IC50']*1000  # nM to uM\n",
    "        final_clean_DTI.at[i,'unit']='uM'\n",
    "    elif final_clean_DTI.at[i,'unit']=='mM':\n",
    "        final_clean_DTI.at[i,'IC50']=final_clean_DTI.at[i,'IC50']/1000  # mM to uM\n",
    "        final_clean_DTI.at[i,'unit']='uM'\n",
    "    elif final_clean_DTI.at[i,'unit']=='uM':\n",
    "        pass\n",
    "final_clean_DTI.to_csv('data/final_clean_DTI.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 5: Download fasta files for all proteins\n",
    "mapping = pd.read_csv('data/chembl2uniprot.txt', header=None, sep='\\t')\n",
    "mapping = mapping[mapping[0].isin(DTI_index['target'].unique())]\n",
    "for uniprot_ID in tqdm(mapping[1].tolist()):\n",
    "    wget.download('https://www.uniprot.org/uniprot/'+uniprot_ID+'.fasta', out='data/fasta_968')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Counter({False: 61624, True: 15146})"
     },
     "metadata": {},
     "execution_count": 57
    }
   ],
   "source": [
    "collections.Counter(DTI_index.duplicated(['target','drug']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Counter({'intermediate': 25589, 'inactive': 28748, 'active': 7287})"
     },
     "metadata": {},
     "execution_count": 67
    }
   ],
   "source": [
    "collections.Counter(DTI_index['activity'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": "100%|██████████| 61624/61624 [01:03<00:00, 976.61it/s]\n100%|██████████| 968/968 [00:03<00:00, 250.94it/s]\n"
    }
   ],
   "source": [
    "# Step 6\n",
    "# Import and process level 2 data\n",
    "drug_smiles = pd.read_csv('data/primary_data/smiles.csv')[['drug','smile']]  # Drug index dataframe\n",
    "drug_smiles['seq_char_count'] = drug_smiles['smile'].apply(lambda x: len(str(x)))   # Add column with character count\n",
    "drug_smiles = drug_smiles[drug_smiles['smile'].notnull()]   # Remove drugs with no smile strings\n",
    "\n",
    "DTI_index = pd.read_csv('data/primary_data/final_clean_DTI.csv')[['target','drug','IC50','unit']]\n",
    "DTI_index = DTI_index[DTI_index['drug'].isin(drug_smiles['drug'].tolist())]     # Exclude drugs for which smiles are not available\n",
    "DTI_index = DTI_index[~DTI_index.duplicated(['target','drug'])]     # Remove duplicate datapoints\n",
    "DTI_index = DTI_index.reset_index(drop=True)\n",
    "\n",
    "mapping = pd.read_csv('data/primary_data/chembl2uniprot.txt', header=None, sep='\\t')\n",
    "mapping = mapping[mapping[0].isin(DTI_index['target'].unique())]    # Select targets that are present in data (DTI_index)\n",
    "targets = os.listdir('data/fasta_968')  # List all target fasta files\n",
    "\n",
    "# Categorize drug-target pairs by IC50 values (Make labels)\n",
    "act = []\n",
    "bct = []\n",
    "for i in tqdm(range(DTI_index.shape[0])):\n",
    "    if DTI_index['IC50'][i]<=0.1:\n",
    "        act.append('active')\n",
    "    elif DTI_index['IC50'][i]>0.1 and DTI_index['IC50'][i]<=30:\n",
    "        act.append('intermediate')\n",
    "    elif DTI_index['IC50'][i]>30:\n",
    "        act.append('inactive')\n",
    "    bct.append(mapping[mapping[0]==DTI_index['target'][i]][1].values[0])\n",
    "    \n",
    "DTI_index['activity'] = act\n",
    "DTI_index['target_uniprot'] = bct\n",
    "\n",
    "# Step 7\n",
    "# Fetch fasta sequences form files and create target index dataframe\n",
    "def fetchFasta(targets):\n",
    "    target_seq = pd.DataFrame(columns=['target_uniprot','target_chembl','seq'])\n",
    "    for fasta in tqdm(targets):\n",
    "        if fasta.split('.')[0] in mapping[1].tolist():\n",
    "            f = open('data/fasta_968/'+fasta,'r')\n",
    "            lines = \"\".join(line.strip() for line in f.readlines()[1:])\n",
    "            dict = {'target_uniprot':fasta.split('.')[0],'target_chembl':mapping[(mapping[1]==fasta.split('.')[0])][0].values[0], 'seq':lines}\n",
    "            target_seq = target_seq.append(dict, True)\n",
    "            f.close()\n",
    "    return target_seq\n",
    "\n",
    "target_seq = fetchFasta(targets)\n",
    "# Length of sequence in train data.\n",
    "#target_seq['seq_char_count'] = target_seq['seq'].apply(lambda : len(x))\n",
    "\n",
    "DTI_index.to_csv('data/DTI_index.csv')\n",
    "target_seq.to_csv('data/target_seq.csv')\n",
    "drug_smiles.to_csv('data/drug_smiles.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Module II: Feature retrieval and processing\n",
    "Use this block to retrieve/process your data by providing a drug file\n",
    "\"\"\"\n",
    "# Extract drug descriptors\n",
    "attribs = [\n",
    "      'MaxAbsPartialCharge', 'MinPartialCharge', 'MinAbsPartialCharge', 'HeavyAtomMolWt', 'MaxAbsEStateIndex', 'NumRadicalElectrons',\n",
    "      'NumValenceElectrons', 'MinAbsEStateIndex', 'MaxEStateIndex', 'MaxPartialCharge', 'MinEStateIndex', 'ExactMolWt', 'MolWt', 'BalabanJ',\n",
    "      'BertzCT', 'Chi0', 'Chi0n', 'Chi0v', 'Chi1', 'Chi1n', 'Chi1v', 'Chi2n', 'Chi2v', 'Chi3n', 'Chi3v', 'Chi4n', 'Chi4v',\n",
    "      'HallKierAlpha', 'Ipc', 'Kappa1', 'Kappa2', 'Kappa3', 'LabuteASA', 'PEOE_VSA1', 'PEOE_VSA10', 'PEOE_VSA11', 'PEOE_VSA12', 'PEOE_VSA13', \n",
    "      'PEOE_VSA14', 'PEOE_VSA2', 'PEOE_VSA3', 'PEOE_VSA4', 'PEOE_VSA5', 'PEOE_VSA6', 'PEOE_VSA7', 'PEOE_VSA8', 'PEOE_VSA9', 'SMR_VSA1',\n",
    "      'SMR_VSA10', 'SMR_VSA2', 'SMR_VSA3', 'SMR_VSA4', 'SMR_VSA5', 'SMR_VSA6', 'SMR_VSA7', 'SMR_VSA8', 'SMR_VSA9', 'SlogP_VSA1', 'SlogP_VSA10', 'SlogP_VSA11', 'SlogP_VSA12',      'SlogP_VSA2','SlogP_VSA3', 'SlogP_VSA4', 'SlogP_VSA5', 'SlogP_VSA6', 'SlogP_VSA7', 'SlogP_VSA8', 'SlogP_VSA9', 'TPSA', 'EState_VSA1', 'EState_VSA10', \n",
    "      'EState_VSA11', 'EState_VSA2',      'EState_VSA3', 'EState_VSA4', 'EState_VSA5', 'EState_VSA6', 'EState_VSA7', 'EState_VSA8', 'EState_VSA9',\n",
    "      'VSA_EState1', 'VSA_EState10', 'VSA_EState2', 'VSA_EState3', 'VSA_EState4', 'VSA_EState5', 'VSA_EState6', 'VSA_EState7', 'VSA_EState8',\n",
    "      'VSA_EState9', 'FractionCSP3', 'HeavyAtomCount', 'NHOHCount', 'NOCount', 'NumAliphaticCarbocycles', 'NumAliphaticHeterocycles',\n",
    "      'NumAliphaticRings', 'NumAromaticCarbocycles', 'NumAromaticHeterocycles', 'NumAromaticRings', 'NumHAcceptors', 'NumHDonors', 'NumHeteroatoms',\n",
    "      'NumRotatableBonds', 'NumSaturatedCarbocycles','NumSaturatedHeterocycles', 'NumSaturatedRings', 'RingCount', 'MolLogP', 'MolMR'\n",
    "    ]\n",
    "\n",
    "fset_mol=pd.DataFrame()\n",
    "for chembl_id in tqdm(drugs['drugs']):\n",
    "    descriptors, status = obj.mol2des(chembl_id, attribs)\n",
    "    if status:\n",
    "        fset_mol = fset_mol.append(descriptors)\n",
    "fset_mol.to_csv('fset_drugs.csv')\n",
    "\n",
    "# Extract protein descriptors\n",
    "#des = obj.extractProteinDes('AAC', 'data/chembl2uniprot.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": "100%|██████████| 14/14 [53:05<00:00, 227.52s/it]\n"
    }
   ],
   "source": [
    "# Extract protein descriptors\n",
    "chembl2uniprot = pd.read_csv('data/chembl2uniprot961.txt', sep='\\t', header=None)\n",
    "for desType in tqdm(['AAC', 'CKSAAP', 'CTDC', 'CTDD', 'CTDT', 'DDE', 'DPC', 'GAAC', 'GDPC', 'GTPC', 'TPC', 'QSOrder', 'KSCTriad', 'CTriad']):\n",
    "    des = obj.extractProteinDes(desType, 'data/chembl2uniprot961.txt', 'data/fset_protein')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract smile strings for every drug\n",
    "smile = []\n",
    "for drug in tqdm(final_clean_DTI['drug'].unique()):\n",
    "    smile.append(new_client.activity.filter(molecule_chembl_id=drug).only(['canonical_smiles'])[0]['canonical_smiles'])\n",
    "smile_df = pd.DataFrame(final_DTI['drug'].unique(), columns=['drug'])\n",
    "smile_df['smile'] = smile\n",
    "smile_df.to_csv('data/smiles.csv')"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9-final"
  },
  "orig_nbformat": 2,
  "kernelspec": {
   "name": "python36964bitdrugnn2venv371fef40068b4ff6bd9b65af96a197ea",
   "display_name": "Python 3.6.9 64-bit ('drugnn2': venv)"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}