Code/Machine_Learning/DL_Code.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "oqZII3qwkwKj"
   },
   "source": [
    "## A Mechanism-aware and Multi-omic Machine Learning Pipeline Characterises Yeast Cell Growth \n",
    "\n",
    "**Deep Learning Modelling Codebase**\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WGq1z00clQRV"
   },
   "source": [
    "We begin with importing the necessary files, we use Keras with Tensorflow 2.0 and the Shap Library for analysing feature importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 505
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 34723,
     "status": "ok",
     "timestamp": 1587994410029,
     "user": {
      "displayName": "Chris Culley",
      "photoUrl": "",
      "userId": "03467560960887582674"
     },
     "user_tz": -60
    },
    "id": "YxwmcuvpJ5eO",
    "outputId": "24aff8d1-f158-4ed9-879c-11fb9d989aeb"
   },
   "outputs": [],
   "source": [
    "import pyreadr\n",
    "import csv\n",
    "import tensorflow as tf\n",
    "from tensorflow.keras.callbacks import EarlyStopping\n",
    "import numpy as np\n",
    "import pandas as pa\n",
    "import os.path\n",
    "import matplotlib.pyplot as plts\n",
    "import seaborn as sns\n",
    "from tensorflow.keras import backend as F\n",
    "from tensorflow.keras.models import Sequential, Model, load_model\n",
    "from tensorflow.keras.layers import Dense, Dropout, Activation, Concatenate, Input\n",
    "from tensorflow.keras.optimizers import SGD\n",
    "from tensorflow.keras.constraints import max_norm\n",
    "from sklearn import preprocessing\n",
    "from tensorflow.keras.utils import plot_model\n",
    "import pandas as pd \n",
    "import shap\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "SEED = 120 # fixed for reproducability\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "kkznlSmgJ_t7"
   },
   "outputs": [],
   "source": [
    "#Checks that the reading of the data has worked\n",
    "def check_file_read_ok(data, name):\n",
    "    if data is None:\n",
    "        print(\"error with \" + name + \" data not read\")\n",
    "    else:\n",
    "        print(\"Success in loading \" + name)\n",
    "        print(data.shape)\n",
    "\n",
    "def remove_zero_entry_columns(data):\n",
    "    return data.loc[:, (data != 0).any(axis = 0)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "aQmWGt5eJ-BB"
   },
   "outputs": [],
   "source": [
    "\n",
    "#Name of the row we are trying to predict (growth rate)\n",
    "TARGET_NAME = 'log2relT'\n",
    "\n",
    "\n",
    "file_loc = ''\n",
    "\n",
    "with open(file_loc + 'testing_index.csv', 'r') as csvfile:\n",
    "    testing_index = []\n",
    "    for row in csv.reader(csvfile, delimiter=';'):\n",
    "        testing_index.append(row[0]) # careful here with [0]\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kRDDFfCvluCD"
   },
   "source": [
    "We next import all the datasets, these are a mix of: \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "*   **expression_data** - refers to the msb145172 dataset (Duibhir et al., 2014) containing transcription rates along with accompanying growth rates, 'log2relT' variable.  \n",
    "*   **metabolic_expression_data** - refers to a subset of the expression data which are used in a metabolic modelling stage\n",
    "* **reaction_data** - refers to the set of reaction flux rates after feeding in the expression data into a flux balance analyse optimised yeast metabolic model \n",
    "*  **iRF** - refers to the a subset of the features (both expression and reaction) that are selected for importance using the Iterative Random Forest algorithm \n",
    "* **SGL** - refers to the selected features selected by the sparse-group-lasso feature selection algorithm \n",
    "* **Gens** - refers to features selected by the NSGA-II algorithm \n",
    "* **independent_data** - refers to an experimentally independent test set captured in a environment that closely resembles the msb expression data and is put through metabolic modelling to create an independent set to test generalisability \n",
    "* **full_data** - contains the expression data and metabolic expression data combined\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 413
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 11302,
     "status": "ok",
     "timestamp": 1587995092107,
     "user": {
      "displayName": "Chris Culley",
      "photoUrl": "",
      "userId": "03467560960887582674"
     },
     "user_tz": -60
    },
    "id": "XFYKNk_WKC7e",
    "outputId": "4e096fe0-8728-4750-f40f-873070682fde"
   },
   "outputs": [],
   "source": [
    "#Load the data removing any\n",
    "full_data = pyreadr.read_r('data/completeDataset.RDS')[None]\n",
    "check_file_read_ok(full_data, \"full data\")\n",
    "\n",
    "# Here we extract the features highlighted by the genetic algorithm\n",
    "def get_the_genetic_algo_data():\n",
    "    gens = []\n",
    "    features = pa.read_csv(file_loc + \"data/genetic_feature_selection_features.csv\", header = None)\n",
    "    for x in range(features.shape[0]):\n",
    "        fet = features.iloc[x, :].values\n",
    "        gens.append(full_data[fet])\n",
    "\n",
    "    return gens\n",
    "\n",
    "gens = get_the_genetic_algo_data()\n",
    "\n",
    "full_data = remove_zero_entry_columns(full_data)\n",
    "expression_data = pyreadr.read_r('data/expressionOnly.RDS')[None]\n",
    "\n",
    "check_file_read_ok(expression_data, \"expression data\")\n",
    "expression_data = remove_zero_entry_columns(expression_data)\n",
    "\n",
    "\n",
    "metabolic_expression_data = pyreadr.read_r('data/metabolic_gene_data.RDS')[None]\n",
    "check_file_read_ok(metabolic_expression_data, \"met_expression data\")\n",
    "\n",
    "#Extract the target and drop target column from main data\n",
    "target_data = full_data[TARGET_NAME]\n",
    "full_data = full_data.drop(columns=TARGET_NAME)\n",
    "expression_data = expression_data.drop(columns=TARGET_NAME)\n",
    "\n",
    "\n",
    "iRF = pa.read_csv(file_loc + 'data/Features_Extracted_Using_iRF.csv', header = None)\n",
    "iRF.columns = ['Genes']\n",
    "check_file_read_ok(iRF, \"iRF data\")\n",
    "iRF = remove_zero_entry_columns(iRF)\n",
    "iRF = full_data[iRF.Genes]\n",
    "\n",
    "#SGL feature extracted data \n",
    "sgl = pa.read_csv(file_loc + 'data/Features_Extracted_Using_SGL.csv', header = None)\n",
    "sgl.columns = ['Genes']\n",
    "check_file_read_ok(sgl, \"sgl data\")\n",
    "sgl_list= []\n",
    "for x in sgl.Genes:\n",
    "    if x in full_data.columns:\n",
    "        sgl_list.append(x)   #We do this because some of the zero features have been removed from full_data\n",
    "\n",
    "sgl = full_data[sgl_list]\n",
    "knockouts = full_data['Row']\n",
    "\n",
    "full_data = full_data.drop(columns = 'Row')\n",
    "reaction_data = full_data.drop(columns = expression_data.columns.values)\n",
    "\n",
    "# Independent dataset\n",
    "independent_data = pd.read_csv(file_loc + 'data/completeDataset_ITS_gamma1.csv')\n",
    "independent_knockouts = independent_data['Row']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imputting missing values in the independent set \n",
    "\n",
    "\n",
    "We next look to impute any missing features in the experimental independent test set "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We need to impute the columns that are missing from the independent dataset \n",
    "to_be_predicted = list(set(full_data.columns) - set(independent_data.columns))\n",
    "# There are a couple in the independent that are not in our full dataset \n",
    "to_be_dropped_from_independent = list(set(independent_data.columns) - set(full_data.columns))\n",
    "\n",
    "\n",
    "independent_data = independent_data.drop(to_be_dropped_from_independent, axis = 1)\n",
    "targets_impute  = full_data[to_be_predicted]\n",
    "\n",
    "original_data_train = full_data.drop(to_be_predicted, axis = 1)\n",
    "data_to_be_predicted = independent_data.copy() \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We apply linear regression to predict the missing features using the original full (msb) dataset as training "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for x in range(targets_impute.shape[1]):\n",
    "    lm = LinearRegression()\n",
    "    lm.fit(original_data_train, targets_impute.iloc[:, x])\n",
    "    predictions = lm.predict(original_data_train)\n",
    "    predicted_values = lm.predict(data_to_be_predicted)\n",
    "    independent_data[targets_impute.columns[x]] = predicted_values "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We remove any duplicate gene knockouts that are already in the full dataset "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_rows = list(set(independent_knockouts) - set(knockouts))\n",
    "independent_lr = independent_data[independent_knockouts.isin(new_rows)]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "T7SjamvKoNYc"
   },
   "source": [
    "We later will use this experimentally independent data in varifying the generalisability of the models "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "CEdx7lbdbxoI"
   },
   "outputs": [],
   "source": [
    "#Independent data with linear regression imputting the missing values \n",
    "\n",
    "independent_target = pd.read_csv(file_loc + 'data/independent_target.csv')\n",
    "independent_genes = independent_knockouts\n",
    "\n",
    "def split_expression_and_fluxes(data):\n",
    "  return data[expression_data.columns], data[reaction_data.columns]\n",
    "\n",
    "independent_lr_expression, independent_lr_fluxes = split_expression_and_fluxes(independent_lr)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RTyZdJktoajA"
   },
   "source": [
    "We build a 3 layer feed-forward network using Keras for the single-view models and a multi-modal model which using a concatenation layer to join two transfer learned single-view models "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "v2LdwT-LLP1J"
   },
   "outputs": [],
   "source": [
    "\n",
    "def init_model_3layer(input_dim, learning_rate, epochs, momentum, neurons):\n",
    "    model = Sequential()\n",
    "    model.add(Dense(neurons, input_dim = input_dim, kernel_constraint=max_norm(3)))\n",
    "    model.add(Activation('sigmoid'))\n",
    "    model.add(Dropout(0.4))\n",
    "    model.add(Dense(neurons, input_dim = input_dim, kernel_constraint=max_norm(3)))\n",
    "    model.add(Activation('sigmoid'))\n",
    "    model.add(Dropout(0.4))\n",
    "    model.add(Dense(neurons, input_dim = input_dim, kernel_constraint=max_norm(3)))\n",
    "    model.add(Activation('sigmoid'))\n",
    "    model.add(Dropout(0.4))\n",
    "    model.add(Dense(1))\n",
    "    model.add(Activation('linear'))\n",
    "    rms = SGD(lr= learning_rate, decay= learning_rate / epochs, momentum=momentum)\n",
    "    model.compile(loss='mean_absolute_error', optimizer=rms, metrics = [\"mean_absolute_error\", 'mean_squared_error'])\n",
    "    return model\n",
    "\n",
    "def init_model(input_dim, learning_rate, epochs, momentum, neurons, trainable = True):\n",
    "    input = Input(shape = (input_dim,))\n",
    "    layer = Dense(neurons, activation='sigmoid', kernel_constraint=max_norm(3), name = \"expression_1\") (input)\n",
    "    layer = Dropout(rate=0.4) (layer)\n",
    "    layer = Dense(neurons, activation='sigmoid', kernel_constraint=max_norm(3), name = \"expression_2\") (layer)\n",
    "    layer = Dropout(rate=0.4) (layer)\n",
    "    predictions = Dense(1, activation='linear') (layer)\n",
    "    model = Model(inputs = input, outputs = predictions)\n",
    "    rms = SGD(lr= learning_rate, decay= learning_rate / epochs, momentum=momentum)\n",
    "    model.trainable = trainable\n",
    "    if (trainable) :\n",
    "        model.compile(loss='mean_squared_error', optimizer=rms, metrics = [\"mean_absolute_error\"])\n",
    "    return model\n",
    "\n",
    "def init_multi_model(input_dim,input_dim2, learning_rate, epochs, momentum, neurons, reaction_trained, expression_trained):\n",
    "    reaction_input = Input(shape = (input_dim,))\n",
    "    expression_input = Input(shape = (input_dim2,))\n",
    "    comb_layer = Concatenate()([reaction_trained(reaction_input), expression_trained(expression_input)])\n",
    "    comb_layer = Dense(neurons, activation='sigmoid', kernel_constraint=max_norm(3), name = \"last_hidden\") (comb_layer)\n",
    "    predictions = Dense(1, activation='linear') (comb_layer)\n",
    "    model = Model(inputs = [reaction_input,expression_input], outputs = predictions)\n",
    "    rms = SGD(lr= learning_rate, decay= learning_rate / epochs, momentum=momentum)\n",
    "    model.compile(loss='mean_squared_error', optimizer=rms, metrics = [\"mean_absolute_error\"])\n",
    "    return model\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Py4ifD7coui6"
   },
   "source": [
    "Parameter setting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "RhGGEP5FO5CX"
   },
   "outputs": [],
   "source": [
    "\n",
    "SEED = 120 # we fix the seed for reproducability \n",
    "number_of_instances = len(target_data) # how many datapoints do we have \n",
    "testing_index = list(map(int, testing_index)) # fixed testing indexes\n",
    "training_index = np.setxor1d(range(1,number_of_instances), testing_index) # not testing index \n",
    "epochs = 1000 # how many epochs we run the network for in training \n",
    "batches = 256 # the batch size in training \n",
    "momentum = 0.75 \n",
    "lrate = 0.005 \n",
    "neurons = 1000 # how many neurons in each hidden unit\n",
    "validation = 0.1 # we use a small validation set due to restricted dataset size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "fEoVLaAFpSGt"
   },
   "source": [
    "Normalising and general readying of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "FswezowaL04T"
   },
   "outputs": [],
   "source": [
    "\n",
    "target_data_train, target_data_test = target_data.drop(target_data.index[testing_index]), target_data.iloc[testing_index]\n",
    "target_data_test.to_csv(file_loc + 'data/test_growth_target.csv')\n",
    "\n",
    "\n",
    "#Split the data 80:20\n",
    "full_data_train, full_data_test = full_data.drop(full_data.index[testing_index]), full_data.iloc[testing_index, :]\n",
    "expression_data_train, expression_data_test = expression_data.drop(expression_data.index[testing_index]), expression_data.iloc[testing_index, :]\n",
    "metabolic_expression_data_train, metabolic_expression_data_test = metabolic_expression_data.drop(metabolic_expression_data.index[testing_index]), metabolic_expression_data.iloc[testing_index, :]\n",
    "reaction_data_train, reaction_data_test = reaction_data.drop(reaction_data.index[testing_index]), reaction_data.iloc[testing_index, :]\n",
    "target_data_train, target_data_test = target_data.drop(target_data.index[testing_index]), target_data.iloc[testing_index]\n",
    "iRF_train, iRF_test = iRF.drop(iRF.index[testing_index]), iRF.iloc[testing_index, :]\n",
    "sgl_train, sgl_test = sgl.drop(sgl.index[testing_index]), sgl.iloc[testing_index, :]\n",
    "\n",
    "\n",
    "#Preprocessing the data to mean of zero and unit variance - stopped using this for improved results\n",
    "full_scaler = preprocessing.StandardScaler().fit(full_data_train)\n",
    "full_data_scaled_train = full_scaler.transform(full_data_train).astype(np.float32)\n",
    "full_data_scaled_test = full_scaler.transform(full_data_test).astype(np.float32)\n",
    "\n",
    "expression_scaler = preprocessing.StandardScaler().fit(expression_data_train)\n",
    "expression_data_scaled_train = expression_scaler.transform(expression_data_train).astype(np.float32)\n",
    "expression_data_scaled_test = expression_scaler.transform(expression_data_test).astype(np.float32)\n",
    "independent_lr_expression_scaled = expression_scaler.transform(independent_lr_expression).astype(np.float32)\n",
    "np.savetxt(file_loc + 'data/independent_lr_expression_scaled.csv', independent_lr_expression_scaled, delimiter=\",\") # we save the independent test set to be used by the R code \n",
    "\n",
    "\n",
    "reaction_scaler = preprocessing.StandardScaler().fit(reaction_data_train)\n",
    "reaction_data_scaled_train = reaction_scaler.transform(reaction_data_train).astype(np.float32)\n",
    "reaction_data_scaled_test = reaction_scaler.transform(reaction_data_test).astype(np.float32)\n",
    "independent_lr_reaction_scaled = reaction_scaler.transform(independent_lr_fluxes).astype(np.float32)\n",
    "\n",
    "\n",
    "metabolic_expression_scaler = preprocessing.StandardScaler().fit(metabolic_expression_data_train)\n",
    "metabolic_expression_data_scaled_train = metabolic_expression_scaler.transform(metabolic_expression_data_train).astype(np.float32)\n",
    "metabolic_expression_data_scaled_test = metabolic_expression_scaler.transform(metabolic_expression_data_test).astype(np.float32)\n",
    "\n",
    "iRF_scaler = preprocessing.StandardScaler().fit(iRF_train)\n",
    "iRF_scaled_train = iRF_scaler.transform(iRF_train).astype(np.float32)\n",
    "iRF_scaled_test = iRF_scaler.transform(iRF_test).astype(np.float32)\n",
    "\n",
    "sgl_scale = preprocessing.StandardScaler().fit(sgl_train)\n",
    "sgl_scaled_train = sgl_scale.transform(sgl_train).astype(np.float32)\n",
    "sgl_scaled_test = sgl_scale.transform(sgl_test).astype(np.float32)\n",
    "\n",
    "target_train = target_data_train.astype(np.float32)\n",
    "target_test = target_data_test.astype(np.float32)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UFJmiRiNpwmI"
   },
   "source": [
    "Here we define the training experiments and the method to extract out the pre-trained model to be joined in a concatenation layer for the multi-modal experiments "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "3IEFerRHMMfp"
   },
   "outputs": [],
   "source": [
    "def run_experiment_full(name, index_in, data_train, data_test,requires_weight_save = False, requires_prediction_save = True, verbose = False):\n",
    "  model_full = init_model(data_train.shape[1], lrate,epochs, momentum, 1000)\n",
    "  if requires_weight_save:\n",
    "    if os.path.exists(file_loc + 'models/'+name+'_model' + str(index_in) + '.h5'): # should we load the weights? \n",
    "      model_full.load_weights(file_loc + 'models/'+name+'_model' + str(index_in) + '.h5')\n",
    "    else:\n",
    "      model_full.fit(x = data_train, y = target_train, epochs=epochs, batch_size=batches, validation_split=validation, verbose=verbose)\n",
    "      model_full.save_weights(file_loc + 'models/'+name+'_model' + str(index_in) + '.h5') # we save the weights\n",
    "  else: \n",
    "    model_full.fit(x = data_train, y = target_train, epochs=epochs, batch_size=batches, validation_split=validation, verbose=verbose)\n",
    "  score = model_full.evaluate(data_test, target_test, verbose=1)\n",
    "  if requires_prediction_save:\n",
    "  #  if not os.path.exists(file_loc +'predictions/'+name+'_Predictions_' + str(index_in) + '.csv'): # do we already have a prediction set \n",
    "       prediction = model_full.predict_on_batch(data_test)\n",
    "       np.savetxt(fname=file_loc +'predictions/'+name+'_Predictions_' + str(index_in) + \".csv\",  X=prediction, delimiter=',') # save the prediction set \n",
    "  return score\n",
    "\n",
    "# We remove the final layer from the pre-trained model and re-compile to allow it to be joined in a layer layer \n",
    "def load_older_model(data_train, fname, clipped = True):\n",
    "  model = init_model(data_train.shape[1], lrate,epochs, momentum, neurons)\n",
    "  model.load_weights(file_loc + 'models/'+fname+'_model.h5')\n",
    "  rms = SGD(lr= lrate, decay= lrate / epochs, momentum=momentum)\n",
    "  model.compile(loss='mean_squared_error', optimizer=rms, metrics = [\"mean_absolute_error\"])\n",
    "  model.summary()\n",
    "  if clipped:\n",
    "   model.trainable = True\n",
    "   model._layers.pop()\n",
    "   model._layers.pop()\n",
    "   model.outputs = [model.layers[-1].output]\n",
    "   model._layers[-1].output_nodes = []\n",
    "   model = Model(inputs=model.inputs, outputs=model.layers[-1].output)\n",
    "   model.compile(loss='mean_squared_error', optimizer=rms, metrics = [\"mean_absolute_error\"])\n",
    "   model.summary()\n",
    "  else:\n",
    "    return model\n",
    "  return model \n",
    "\n",
    "weights_model_run = []\n",
    "outputs = [] \n",
    "\n",
    "# The multi-modal model is used as a way of integrating two separate transfer learned view-point - fluxes and expressions \n",
    "def run_experiment_mm_full(name, index_in, data_train_1, data_train_2, data_test_1, data_test_2, name_1, name_2, requires_weight_save = False, requires_prediction_save = True, verbose = False):\n",
    "  model_1 = load_older_model(data_train_1, name_1)\n",
    "  model_2 = load_older_model(data_train_2, name_2)\n",
    "  model_full = init_multi_model(data_train_1.shape[1], data_train_2.shape[1], lrate, epochs, 0.75, 15, model_1, model_2)\n",
    "  if requires_weight_save:\n",
    "    if os.path.exists(file_loc + 'models/'+name+'_model' + str(index_in) + '.h5'):\n",
    "      model_full.load_weights(file_loc + 'models/'+name+'_model' + str(index_in) + '.h5')\n",
    "    else:\n",
    "      model_full.fit(x = [data_train_1, data_train_2], y = target_train, epochs=epochs, batch_size=batches, validation_split=validation,  verbose=verbose)\n",
    "      model_full.save_weights(file_loc + 'models/'+name+'_model' + str(index_in) + '.h5')\n",
    "  else: \n",
    "    model_full.fit(x = [data_train_1, data_train_2], y = target_train, epochs=epochs, batch_size=batches, validation_split=validation, verbose=verbose)\n",
    "  score = model_full.evaluate([data_test_1, data_test_2], target_test, verbose=1)\n",
    "  if requires_prediction_save:\n",
    "    if not os.path.exists(file_loc +'predictions/'+name+'_Predictions_' + str(index_in) + '.csv'):\n",
    "       prediction = model_full.predict_on_batch([data_test_1, data_test_2])\n",
    "       np.savetxt(fname=file_loc +'predictions/'+name+'_Predictions_' + str(index_in) + \".csv\",  X=prediction, delimiter=',')\n",
    "  if name_2 == 'expression':\n",
    "    weights_model_run.append(model_full.get_layer(name = \"last_hidden\").get_weights()[0])\n",
    "    outputs.append(model_full.layers[-1].get_weights())\n",
    "  return score\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "MPvD3VQ5qonS"
   },
   "source": [
    "We run each set of experiments for 100 runs so that we can gain confidence scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "svcWrlLaNNha"
   },
   "outputs": [],
   "source": [
    "def run_nsga(): \n",
    "  for i in range(0,10):\n",
    "    for x in range(0,9):\n",
    "          next = gens[5]\n",
    "          next_train, next_test = next.drop(next.index[testing_index]), next.iloc[testing_index, :]\n",
    "          scale = preprocessing.StandardScaler().fit(next_train)\n",
    "          next_scale_train = scale.transform(next_train).astype(np.float32)\n",
    "          next_scale_test = scale.transform(next_test).astype(np.float32)\n",
    "\n",
    "          score = run_experiment_full('NSGA-II', i*10+x, next_scale_train, next_scale_test)\n",
    "          print(score)\n",
    "\n",
    "def run_expression():\n",
    "  for i in range(0, 101):\n",
    "    score = run_experiment_full('expression', i, expression_data_scaled_train, expression_data_scaled_test)\n",
    "    print(score)\n",
    "\n",
    "#GEM\n",
    "def run_metabolic():\n",
    "  for i in range(0,101):\n",
    "    score = run_experiment_full('metabolic_expression', i, metabolic_expression_data_scaled_train, metabolic_expression_data_scaled_test)\n",
    "    print(score)\n",
    "\n",
    "def run_fluxes():\n",
    "  for i in range(0,101):\n",
    "    score = run_experiment_full('fluxes', i, reaction_data_scaled_train, reaction_data_scaled_test)\n",
    "    print(score)\n",
    "\n",
    "def run_concate_flu_ge():\n",
    "  for i in range(0,101):\n",
    "    score = run_experiment_full('concate_Flu_GE', i, full_data_scaled_train, full_data_scaled_test)\n",
    "    print(score)\n",
    "\n",
    "def run_irf():\n",
    "  for i in range(0,101):\n",
    "    score = run_experiment_full('iRF', i, iRF_scaled_train, iRF_scaled_test)\n",
    "    print(score)\n",
    "\n",
    "def run_sgl():\n",
    "  for i in range(0,101):\n",
    "    \n",
    "    score = run_experiment_full('SGL', i, sgl_scaled_train, sgl_scaled_test)\n",
    "    print(score)\n",
    "\n",
    "\n",
    "def run_mm_metabolic():\n",
    "  for i in range(0,101):\n",
    "    score = run_experiment_mm_full('multi_model_metabolic_expression', i, reaction_data_scaled_train, \n",
    "                       metabolic_expression_data_scaled_train, reaction_data_scaled_test, \n",
    "                       metabolic_expression_data_scaled_test, 'reaction', 'metabolic_expression')\n",
    "    print(score)\n",
    "\n",
    "def run_mm_full():\n",
    "  for i in range(0,101):\n",
    "    score = run_experiment_mm_full('multi_model_full_expression', i, reaction_data_scaled_train, \n",
    "                       expression_data_scaled_train, reaction_data_scaled_test, \n",
    "                       expression_data_scaled_test, 'reaction', 'expression', verbose = False)\n",
    "    print(score)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "B5pAXrTLrF3Y"
   },
   "source": [
    "# Training all of the models "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 807505,
     "status": "error",
     "timestamp": 1586441764589,
     "user": {
      "displayName": "Chris Culley",
      "photoUrl": "",
      "userId": "03467560960887582674"
     },
     "user_tz": -60
    },
    "id": "R83n6aNCAG2G",
    "outputId": "0e46b3bb-e9a1-4579-b285-24b0c88c0ae1"
   },
   "outputs": [],
   "source": [
    "run_expression()\n",
    "run_fluxes()\n",
    "run_mm_full()\n",
    "run_mm_metabolic()\n",
    "run_sgl()\n",
    "run_irf()\n",
    "run_concate_flu_ge()\n",
    "run_nsga()\n",
    "run_metabolic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "fGeaatPsrIoG"
   },
   "source": [
    "# Testing on the independent test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 5832,
     "status": "ok",
     "timestamp": 1587989568077,
     "user": {
      "displayName": "Chris Culley",
      "photoUrl": "",
      "userId": "03467560960887582674"
     },
     "user_tz": -60
    },
    "id": "_svIgLKK9PN0",
    "outputId": "9fdd33c8-f09c-410c-a7fd-7064ca2e594b"
   },
   "outputs": [],
   "source": [
    "# Here load in the multi-modal model \n",
    "\n",
    "expression_side = load_older_model(independent_lr_expression, 'expression')\n",
    "flux_side = load_older_model(independent_lr_fluxes, 'reaction')\n",
    "\n",
    "multi_modal_independent = init_multi_model(independent_lr_fluxes.shape[1], independent_lr_expression.shape[1], lrate, epochs, 0.75, 15, flux_side, expression_side)\n",
    "multi_modal_independent.load_weights(\"models/MM-Flu_GE.h5\")\n",
    "\n",
    "# Predict on the independent set multi-modal\n",
    "independent_predictions_lr_mm = multi_modal_independent.predict([ independent_lr_reaction_scaled, independent_lr_expression_scaled])\n",
    "\n",
    "# Predict on independent set expression only \n",
    "expression_model = load_older_model(independent_lr_expression, 'expression', clipped = False)\n",
    "independent_predictions_lr_e_dl = expression_model.predict(independent_lr_expression_scaled)\n",
    "\n",
    "checking = multi_modal_independent.predict([reaction_data_scaled_test, expression_data_scaled_test])\n",
    "np.savetxt('predictions/Independent/independent_lr_multi_modal_expression_fluxes.csv', independent_predictions_lr_mm, delimiter=\",\")\n",
    "np.savetxt('predictions/Independent/independent_lr_dl_expression_only.csv', independent_predictions_lr_e_dl, delimiter=\",\")\n",
    "\n",
    "double_gene_knockouts = independent_genes.apply(lambda x : '_' in x)\n",
    "np.savetxt('data/is_double_gene_knockout.csv', double_gene_knockouts, delimiter=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Om2-lgwqrOLz"
   },
   "source": [
    "# Calculating the feature importance using the SHAP method "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 53
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 11425,
     "status": "ok",
     "timestamp": 1586532869421,
     "user": {
      "displayName": "Chris Culley",
      "photoUrl": "",
      "userId": "03467560960887582674"
     },
     "user_tz": -60
    },
    "id": "WaFBN23Zr-oX",
    "outputId": "bd2c64cc-022c-4b18-bb61-8e0f0b20596e"
   },
   "outputs": [],
   "source": [
    "# Load the previously saved model \n",
    "multi_modal = load_model(\"models/MM-Flu_GE_n.h5\")\n",
    "\n",
    "# create the background set from which the mean can be generated\n",
    "background = [reaction_data_scaled_train, expression_data_scaled_train]\n",
    "\n",
    "# train the explainer \n",
    "explainer = shap.DeepExplainer(multi_modal, background)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DC_WLmB6r054"
   },
   "source": [
    "We calculate a normalised error for each data point so that we can determine which feature importance emphasis cause the issues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "W911f_Y3XwmQ"
   },
   "outputs": [],
   "source": [
    "predictions = multi_modal.predict([reaction_data_scaled_test, expression_data_scaled_test])\n",
    "errors = np.abs(predictions.reshape(228,) - target_test)\n",
    "norm_errors = (errors - np.mean(errors)) / np.std(errors)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "r8_XcgHus86R"
   },
   "source": [
    "Shap value calculation to gather feature importance for each of the inputs which export to csv format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "qpg7UYwd4LUD"
   },
   "outputs": [],
   "source": [
    "test_all = [reaction_data_scaled_test, expression_data_scaled_test]\n",
    "shap_values_all = explainer.shap_values(test_all)\n",
    "\n",
    "reaction_shap_values = pd.DataFrame(data = shap_values_all[0][0], columns = reaction_data.columns)\n",
    "reaction_shap_values['target'] = target_test.values\n",
    "reaction_shap_values['Normalised_Errors'] = norm_errors.values\n",
    "expression_shap_values = pd.DataFrame(data = shap_values_all[0][1], columns = expression_data.columns)\n",
    "expression_shap_values['Normalised_Errors'] = norm_errors.values\n",
    "expression_shap_values['target'] = target_test.values\n",
    "\n",
    "reaction_shap_values.to_csv('data/reaction_shap_values.csv')\n",
    "expression_shap_values.to_csv('data/expression_shap_values.csv')"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "authorship_tag": "ABX9TyN7xMNNGWNQ7DrdP0xLenI2",
   "collapsed_sections": [],
   "name": "yeast_growth_DL.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}