train&predict.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Change pathogens into the set of taxa you want to predict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pathogens = set() # Taxa need to be in the form of URIs, i.e. <http://purl.obolibrary.org/obo/NCBITaxon_11676>\n",
    "# For example, we want to predict for Zika virus\n",
    "pathogens.add('<http://purl.obolibrary.org/obo/NCBITaxon_64320>')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading data the same way as in train.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "bacteria = set()\n",
    "virus = set()\n",
    "species_dict = {}\n",
    "with open('data/ncbitaxon/categories.dmp', 'r') as f:\n",
    "    for line in f:\n",
    "        items = line.strip().split('\\t')\n",
    "        species_dict['<http://purl.obolibrary.org/obo/NCBITaxon_'+items[2]+'>'] = '<http://purl.obolibrary.org/obo/NCBITaxon_' + items[1]+'>'\n",
    "        if items[0] == 'B':\n",
    "            bacteria.add('<http://purl.obolibrary.org/obo/NCBITaxon_'+items[2]+'>')\n",
    "        elif items[0] == 'V':\n",
    "            virus.add('<http://purl.obolibrary.org/obo/NCBITaxon_'+items[2]+'>')\n",
    "print(len(bacteria), len(virus), len(species_dict))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import pickle \n",
    "\n",
    "\n",
    "embedding_file = 'opa2vec/backup/embeddings_patho_ncbitaxon_short.out'\n",
    "data = pd.read_csv(embedding_file, header = None, sep = ' ', skiprows=0)\n",
    "embds_data = data.values\n",
    "patho_dict = dict(zip(embds_data[:,0],embds_data[:,1:]))\n",
    "print('finished reading embeddings 1')\n",
    "\n",
    "embedding_file = 'opa2vec2/backup/embeddings_HPiMPiGO_short.out'\n",
    "data = pd.read_csv(embedding_file, header = None, sep = ' ', skiprows=0)\n",
    "embds_data = data.values\n",
    "host_dict = dict(zip(embds_data[:,0],embds_data[:,1:]))\n",
    "print('finished reading embeddings 2')\n",
    "\n",
    "protein_set = set()\n",
    "with open('host_pheno_asso/human_pheno_asso_HPiMPiGO.txt','r') as f:\n",
    "    for line in f:\n",
    "        items = line.strip().split()\n",
    "        protein_set.add(items[0])\n",
    "print(len(protein_set))\n",
    "\n",
    "hosts = set()\n",
    "pathos = set()\n",
    "positives_set = set()\n",
    "count_dict = {}\n",
    "with open('data/hpidb2.score.txt', 'r') as f:\n",
    "    for line in f:\n",
    "        items = line.strip().split('\\t')\n",
    "        patho = '<http://purl.obolibrary.org/obo/NCBITaxon_' + items[1] + '>'\n",
    "        if ':' in items[2]:\n",
    "            if float(items[2].split('miscore:')[1]) >= 0.5:\n",
    "                hosts.add(items[0])\n",
    "                pathos.add(patho)\n",
    "                positives_set.add((items[0], patho))\n",
    "print(len(hosts), len(pathos), len(positives_set))\n",
    "\n",
    "positives = set()\n",
    "for pair in positives_set:\n",
    "    if (pair[0] in host_dict and  pair[0] in protein_set) and (pair[1] in patho_dict and pair[1] in virus):\n",
    "        positives.add(pair)\n",
    "print(len(positives))\n",
    "\n",
    "pathogens_counts = {}\n",
    "protein_counts = {}\n",
    "for positive in positives:\n",
    "    if positive[1] not in pathogens_counts:\n",
    "        pathogens_counts[positive[1]] = 0\n",
    "    pathogens_counts[positive[1]]+=1\n",
    "    if positive[0] not in protein_counts:\n",
    "        protein_counts[positive[0]] = 0\n",
    "    protein_counts[positive[0]]+=1\n",
    "print(len(pathogens_counts), len(protein_counts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "from keras.models import Sequential\n",
    "from keras.layers import Dense, Dropout, BatchNormalization, Activation\n",
    "from keras.utils import multi_gpu_model, Sequence, np_utils\n",
    "import math\n",
    "from keras.optimizers import SGD, Adam, RMSprop, Nadam\n",
    "from keras.callbacks import EarlyStopping, TensorBoard\n",
    "import scipy.stats as ss\n",
    "import matplotlib.pyplot as plt\n",
    "from keras.backend.tensorflow_backend import set_session\n",
    "from keras import backend as K\n",
    "from hyperopt import Trials, STATUS_OK, tpe\n",
    "from hyperas import optim\n",
    "from hyperas.distributions import choice, uniform\n",
    "\n",
    "class Generator(Sequence):\n",
    "    def __init__(self, x_set, y_set, batch_size):\n",
    "        self.x, self.y = x_set, y_set\n",
    "        self.batch_size = batch_size\n",
    "        self.nbatch = int(np.ceil(len(self.x) / float(self.batch_size)))\n",
    "        self.length = len(self.x)\n",
    "\n",
    "    def __len__(self):\n",
    "        return self.nbatch\n",
    "\n",
    "    def __getitem__(self, idx):\n",
    "        #print('idx: ',idx)\n",
    "        start = idx * self.batch_size\n",
    "        batch_len = min(self.batch_size, (self.length)-start)\n",
    "        X_batch_list = np.empty((batch_len, 600), dtype=np.float32)\n",
    "        y_batch_list = np.empty(batch_len, dtype=np.float32)\n",
    "       \n",
    "        for ids in range(start, min((idx + 1) * self.batch_size, self.length)):\n",
    "            array1 = host_dict[self.x[ids][0]]\n",
    "            array2 = patho_dict[self.x[ids][1]]\n",
    "            embds = np.concatenate([array1, array2])\n",
    "            X_batch_list[ids-start] = embds\n",
    "            y_batch_list[ids-start] = self.y[ids]\n",
    "           \n",
    "        return X_batch_list, y_batch_list\n",
    "\n",
    "\n",
    "# In[5]:\n",
    "\n",
    "\n",
    "config = tf.ConfigProto()\n",
    "config.gpu_options.allow_growth = True\n",
    "config.gpu_options.per_process_gpu_memory_fraction = 0.3\n",
    "sess = tf.Session(config=config)\n",
    "K.set_session(sess)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training and predicting, saving the results in predictions.tsv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rank_counts = {}\n",
    "write_counts = {}\n",
    "counter = 0\n",
    "output = open('predictions.tsv', 'w')\n",
    "output.close()\n",
    "for taxon in pathogens:\n",
    "    K.clear_session()\n",
    "    counter+=1\n",
    "    print('taxon ', counter)\n",
    "    val_pathos = set()\n",
    "    val_pathos.add(taxon)\n",
    "    train_pathos = set(list(pathogens_counts.keys())) - val_pathos\n",
    "    print('train pathos: ', len(train_pathos))\n",
    "    print(val_pathos)\n",
    "    \n",
    "    train_positives = []\n",
    "    val_positives = []\n",
    "    train_positives_set = set()\n",
    "    val_positives_set = set()\n",
    "    for items in positives:\n",
    "        if items[1] in train_pathos:\n",
    "            train_positives_set.add((items[0], items[1]))\n",
    "            train_positives.append((items[0], items[1], 1))\n",
    "        if items[1] in val_pathos:\n",
    "            val_positives_set.add((items[0], items[1]))\n",
    "            val_positives.append((items[0], items[1], 1))\n",
    "    print(len(train_positives), len(val_positives))\n",
    "\n",
    "    train_negatives = []\n",
    "    train_all_tuples = set()\n",
    "\n",
    "    for host in protein_set:\n",
    "        for patho in train_pathos:\n",
    "            if host in host_dict:\n",
    "                train_all_tuples.add((host, patho))\n",
    "    print(len(train_all_tuples))\n",
    "\n",
    "    for item in train_all_tuples:\n",
    "        if item not in train_positives_set:\n",
    "            train_negatives.append((item[0], item[1], 0))\n",
    "    print(len(train_negatives))\n",
    "\n",
    "    val_negatives = []\n",
    "    val_all_tuples = set()\n",
    "\n",
    "    for host in protein_set:\n",
    "        for patho in val_pathos:\n",
    "            if host in host_dict:\n",
    "                val_all_tuples.add((host, patho))\n",
    "    print(len(val_all_tuples))\n",
    "\n",
    "    for item in val_all_tuples:\n",
    "        if item not in val_positives_set:\n",
    "            val_negatives.append((item[0], item[1], 0))\n",
    "    print(len(val_negatives))\n",
    "\n",
    "    train_positives = np.repeat(np.array(list(train_positives)), len(train_negatives)//len(train_positives), axis = 0)\n",
    "    train_negatives = np.array(train_negatives)\n",
    "    triple_train = np.concatenate((train_positives, train_negatives), axis=0)\n",
    "    np.random.shuffle(triple_train)\n",
    "    triple_val = np.concatenate((val_positives, val_negatives), axis=0)\n",
    "\n",
    "    batch_size = 2**12\n",
    "    num_classes = 1\n",
    "\n",
    "    model = Sequential()\n",
    "    model.add(Dense(256, activation='relu', input_shape=(600,)))\n",
    "    model.add(BatchNormalization())\n",
    "    model.add(Dropout(0.515))\n",
    "    model.add(Dense(1, activation='sigmoid'))\n",
    "    #parallel_model = multi_gpu_model(model, 2)\n",
    "    parallel_model = model\n",
    "    parallel_model.compile(loss='binary_crossentropy',\n",
    "                  optimizer=Adam(),\n",
    "                  metrics=['accuracy'])\n",
    "\n",
    "    generator = Generator(triple_train[:len(triple_train),0:2], triple_train[:len(triple_train),2], batch_size)\n",
    "    val_generator = Generator(triple_val[:,0:2], triple_val[:,2], batch_size)\n",
    "\n",
    "    num_positive = len(val_positives)\n",
    "    history = parallel_model.fit_generator(generator=generator,\n",
    "                        epochs=30,\n",
    "                        steps_per_epoch = int(math.ceil(len(triple_train)/ batch_size)*0.04),\n",
    "                        verbose=2,\n",
    "                        validation_data=generator,\n",
    "                        validation_steps= 1, max_queue_size = 50,\n",
    "                        use_multiprocessing=True,\n",
    "                        workers = 5)\n",
    "    sim_list = parallel_model.predict_generator(generator=Generator(triple_val[:,0:2], triple_val[:,2], 1000), \n",
    "                                                verbose=2, steps=int(math.ceil(math.ceil(len(triple_val))/1000)), \n",
    "                                                max_queue_size = 20, use_multiprocessing=True, workers = 3)\n",
    "    \n",
    "    with open('predictions.tsv', 'a+') as f:\n",
    "        for i in range(triple_val.shape[0]):\n",
    "            f.write('%s\\t%s\\t%f\\n' % (triple_val[i][1],  triple_val[i][0], sim_list[i])) "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}