Motif_Finding_Contextual_Regression.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "XcXY390qJrXA"
      },
      "outputs": [],
      "source": [
        "# check GPU, colab sometimes doesn't give you any GPU...\n",
        "gpu_info = !nvidia-smi\n",
        "gpu_info = '\\n'.join(gpu_info)\n",
        "if gpu_info.find('failed') >= 0:\n",
        "  print('Not connected to a GPU')\n",
        "else:\n",
        "  print(gpu_info)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7ZPKaKvfCDjI"
      },
      "outputs": [],
      "source": [
        "# install h5py, since the original h5py on colab doesn't work\n",
        "pip install 'h5py==2.10.0' --force-reinstall"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9E78Qgodaek5"
      },
      "outputs": [],
      "source": [
        "# download data\n",
        "!wget -r -np -nH --reject \"index.html*\" --cut-dirs 6 \\\n",
        " https://krishna.gs.washington.edu/content/members/vagar/Xpresso/data/datasets/pM10Kb_1KTest/"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5C199MavOtpi"
      },
      "outputs": [],
      "source": [
        "# imports\n",
        "%tensorflow_version 1.x\n",
        "\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive')\n",
        "\n",
        "\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import tensorflow as tf\n",
        "\n",
        "import dill\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "from math import *\n",
        "import keras\n",
        "from keras import datasets, layers, models, Model,regularizers, initializers\n",
        "from keras.models import Model, load_model\n",
        "from keras.utils.vis_utils import plot_model\n",
        "from keras.optimizers import Adam, SGD\n",
        "from keras.constraints import maxnorm, nonneg\n",
        "from keras.activations import sigmoid\n",
        "from keras.layers import *\n",
        "from keras.metrics import *\n",
        "from keras.callbacks import Callback, ModelCheckpoint, EarlyStopping\n",
        "from hyperopt import hp, STATUS_OK"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "AF5Vk_RGtPhj"
      },
      "outputs": [],
      "source": [
        "# load data\n",
        "import h5py\n",
        "import os\n",
        "datadir=\"pM10Kb_1KTest\"\n",
        "trainfile = h5py.File(os.path.join(datadir, 'train.h5'), 'r')\n",
        "X_trainhalflife, X_trainpromoter, y_train, geneName_train = trainfile['data'], trainfile['promoter'], trainfile['label'], trainfile['geneName']\n",
        "validfile = h5py.File(os.path.join(datadir, 'valid.h5'), 'r')\n",
        "X_validhalflife, X_validpromoter, y_valid, geneName_valid = validfile['data'], validfile['promoter'], validfile['label'], validfile['geneName']\n",
        "testfile = h5py.File(os.path.join(datadir, 'test.h5'), 'r')\n",
        "X_testhalflife, X_testpromoter, y_test, geneName_test = testfile['data'], testfile['promoter'], testfile['label'], testfile['geneName']\n",
        "\n",
        "leftend = 3000\n",
        "rightend = 13500\n",
        "\n",
        "leftpos = 9800\n",
        "rightpos = 10200\n",
        "\n",
        "XTrain = X_trainpromoter[:,leftpos:rightpos,:]\n",
        "XVal = X_validpromoter[:,leftpos:rightpos,:]\n",
        "XTest = X_testpromoter[:,leftpos:rightpos,:]\n",
        "\n",
        "YTrain = y_train\n",
        "YVal = y_valid\n",
        "YTest = y_test"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nqeMeGtRifvF"
      },
      "outputs": [],
      "source": [
        "# generate data divisions\n",
        "import random\n",
        "import os\n",
        "import pickle\n",
        "\n",
        "num_of_data_division = 10\n",
        "\n",
        "train_num = len(XTrain)\n",
        "val_num = len(XVal)\n",
        "test_num = len(XTest)\n",
        "\n",
        "def genAndSaveDivisionLists(train_num, val_num, test_num):\n",
        "    total_num = train_num + val_num + test_num\n",
        "    total_set = range(0, total_num)\n",
        "    train_set = random.sample(total_set, train_num)\n",
        "    total_set = set(total_set) - set(train_set)\n",
        "    val_set = random.sample(total_set, val_num)\n",
        "    test_set = list(set(total_set) - set(val_set))\n",
        "    return {'train': train_set, 'valid': val_set, 'test': test_set}\n",
        "\n",
        "for i in range(0, num_of_data_division):\n",
        "    folder_name = \"/Your_Folder/run_\" + str(i) + \"/\"\n",
        "    os.mkdir(folder_name)\n",
        "    division_dict = genAndSaveDivisionLists(train_num, val_num, test_num)\n",
        "    with open(folder_name + \"division_dict\",\"wb\") as f:\n",
        "        pickle.dump(division_dict, f)\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mxDEkH15Ph85"
      },
      "outputs": [],
      "source": [
        "# load a specific data division\n",
        "import pickle\n",
        "\n",
        "total_x = np.concatenate([XTrain, XVal, XTest], axis=0)\n",
        "total_y = np.concatenate([YTrain, YVal, YTest], axis=0)\n",
        "total_gn = np.concatenate([geneName_train, geneName_valid, geneName_test], axis=0)\n",
        "\n",
        "division_set_num = 0\n",
        "\n",
        "folder_name = \"/Your_Folder/run_\" + str(division_set_num) + \"/\"\n",
        "\n",
        "with open(folder_name + \"division_dict_\" + str(division_set_num),\"rb\") as f:\n",
        "    division_dict = pickle.load(f)\n",
        "\n",
        "XTrain = total_x[division_dict['train'],:,:]\n",
        "XVal = total_x[division_dict['valid'],:,:]\n",
        "XTest = total_x[division_dict['test'],:,:]\n",
        "YTrain = total_y[division_dict['train']]\n",
        "YVal = total_y[division_dict['valid']]\n",
        "YTest = total_y[division_dict['test']]\n",
        "\n",
        "print(XTrain.shape)\n",
        "print(XVal.shape)\n",
        "print(XTest.shape)\n",
        "print(YTrain.shape)\n",
        "print(YVal.shape)\n",
        "print(YTest.shape)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6lSE-PgcP7FS"
      },
      "outputs": [],
      "source": [
        "# model\n",
        "def positiveAndSmallRegularizer(x):\n",
        "  return 1e-1 * (0 * keras.backend.mean(tf.square(x)) - keras.backend.mean(x))\n",
        "\n",
        "def strong_sigmoid(x):\n",
        "  return keras.backend.sigmoid(4 * x)\n",
        "\n",
        "class BiasLayer(Layer):\n",
        "  def __init__(self, *args, **kwargs):\n",
        "    super(BiasLayer, self).__init__(*args, **kwargs)\n",
        "\n",
        "  def build(self, input_shape):\n",
        "    self.bias = self.add_weight('bias',\n",
        "                                    shape=input_shape[1:],\n",
        "                                    initializer=initializers.Constant(0),\n",
        "                                    trainable=True)\n",
        "  def call(self, x):\n",
        "    return x + self.bias\n",
        "\n",
        "num_of_motifs = 64\n",
        "motif_by_num = 8\n",
        "\n",
        "layer_0 = Input(shape=XTrain.shape[1:],name='Input')\n",
        "\n",
        "layer_1_fun = Conv1D(num_of_motifs, motif_by_num, padding='valid', kernel_initializer='glorot_normal', kernel_constraint=nonneg(), activity_regularizer=regularizers.l1(0.00000002), name='Conv1D_1', bias=False)\n",
        "layer_1 = layer_1_fun(layer_0)\n",
        "\n",
        "layer_1_over_bp_num = Lambda(lambda x: x, name='layer_1_over_bp_num')(layer_1)\n",
        "\n",
        "motif_exists = Dense(num_of_motifs, kernel_initializer=initializers.Identity(), trainable=False, activity_regularizer=regularizers.l1(0), name='motif_exists')(layer_1)\n",
        "real_motif_counts = Lambda(lambda x: keras.backend.sum(x, axis=1), name='real_motif_counts')(motif_exists)\n",
        "motif_counts = Lambda(lambda x: keras.backend.mean(x, axis=1), name='motif_counts')(motif_exists)\n",
        "\n",
        "layer_2_pre = ReLU(name='ReLU_0')(layer_1)\n",
        "layer_2 = MaxPool1D(50,name='MaxPool_1')(layer_2_pre)\n",
        "\n",
        "layer_3 = Conv1D(16,7,padding='same', kernel_initializer='glorot_normal', kernel_regularizer=regularizers.l2(0), activation='relu',name='Conv1D_2')(layer_2)\n",
        "layer_4 = MaxPool1D(2,name='MaxPool_2')(layer_3)\n",
        "\n",
        "layer_5 = Conv1D(8,7,padding='same', kernel_initializer='glorot_normal', kernel_regularizer=regularizers.l2(0), activation='relu',name='Conv1D_3')(layer_4)\n",
        "layer_6 = MaxPool1D(2,name='MaxPool_3')(layer_5)\n",
        "\n",
        "layer_7 = Flatten(name='Flatten_1')(layer_6)\n",
        "\n",
        "layer_8 = Dense(64,name='Dense_1', kernel_regularizer=regularizers.l2(0.0001))(layer_7)\n",
        "layer_9 = ReLU(name='ReLU_1')(layer_8)\n",
        "layer_10 = Dropout(0.2,name='dropout_1')(layer_9)\n",
        "\n",
        "layer_11 = Dense(64,name='Dense_2', kernel_regularizer=regularizers.l2(0.0001))(layer_10)\n",
        "layer_12 = ReLU(name='ReLU_2')(layer_11)\n",
        "layer_13 = Dropout(0.2,name='dropout_2')(layer_12)\n",
        "\n",
        "layer_14 = Dense(num_of_motifs, name='Contextual_Weight', activity_regularizer=regularizers.l1(0.0001))(layer_13)\n",
        "\n",
        "layer_15 = Multiply(name='Multiply')([motif_counts,layer_14])\n",
        "\n",
        "layer_16 = Lambda(lambda x: keras.backend.sum(x, axis=1), name='Sum')(layer_15)\n",
        "layer_17 = BiasLayer(name='Bias')(layer_16)\n",
        "layer_18 = Lambda(lambda x: keras.backend.expand_dims(x, axis=1), name='Output')(layer_17)\n",
        "\n",
        "model = Model(inputs=layer_0,outputs=layer_18)\n",
        "model.summary()\n",
        "keras.utils.plot_model(model, \"model_with_shape_info.png\", show_shapes=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jGq7y0fXP-cO"
      },
      "outputs": [],
      "source": [
        "# training\n",
        "model.compile(SGD(lr=0.0005, momentum=0.9),'mean_squared_error', metrics=['mean_squared_error'])\n",
        "earlystop_cb = EarlyStopping(monitor='val_loss', patience=10, verbose=1, mode='min')\n",
        "check_cb = ModelCheckpoint(folder_name + 'bestparams_200_with_b_64_filters.h5', monitor='val_loss', verbose=1, save_best_only=True, mode='min')\n",
        "result = model.fit(XTrain, YTrain, batch_size=128, shuffle=\"batch\", epochs=100,\n",
        "        validation_data=[XVal, YVal], callbacks=[earlystop_cb, check_cb])\n",
        "mse_history = result.history['val_mean_squared_error']"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "cVd5s4AZQP4I"
      },
      "outputs": [],
      "source": [
        "# calculate accuracy\n",
        "best_file = folder_name + 'bestparams_200_with_b_64_filters.h5'\n",
        "model = models.load_model(best_file, custom_objects={'keras': keras, 'tf': tf, 'strong_sigmoid': strong_sigmoid, 'positiveAndSmallRegularizer': positiveAndSmallRegularizer,'motif_by_num':motif_by_num, 'BiasLayer':BiasLayer})\n",
        "YTrain_Pred = model.predict(XTrain)\n",
        "print('Correlation of training dataset: '+str(np.corrcoef(YTrain[:].transpose(),YTrain_Pred.transpose())[0,1]))\n",
        "YVal_Pred = model.predict(XVal)\n",
        "print('Correlation of validation dataset: '+str(np.corrcoef(YVal[:].transpose(),YVal_Pred.transpose())[0,1]))\n",
        "YTest_Pred = model.predict(XTest)\n",
        "print('Correlation of test dataset: '+str(np.corrcoef(YTest[:].transpose(),YTest_Pred.transpose())[0,1]))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "18ac_01PojbG"
      },
      "outputs": [],
      "source": [
        "# collect motifs from all runs\n",
        "import numpy as np\n",
        "\n",
        "total_x = np.concatenate([XTrain, XVal, XTest], axis=0)\n",
        "total_y = np.concatenate([YTrain, YVal, YTest], axis=0)\n",
        "total_gn = np.concatenate([geneName_train, geneName_valid, geneName_test], axis=0)\n",
        "\n",
        "all_motifs = []\n",
        "for division_set_num in range(0, 10):\n",
        "    folder_name = \"/Your_Folder/run_\" + str(division_set_num) + \"/\"\n",
        "    best_file = folder_name + 'bestparams_200_with_b_64_filters.h5'\n",
        "    model = models.load_model(best_file, custom_objects={'keras': keras, 'tf': tf, 'strong_sigmoid': strong_sigmoid, 'positiveAndSmallRegularizer': positiveAndSmallRegularizer,'motif_by_num':motif_by_num, 'BiasLayer':BiasLayer})\n",
        "    motif_out = keras.Model(inputs=model.input, outputs=model.get_layer('motif_counts').output)\n",
        "    NN_found_motifs = motif_out.predict(total_x)\n",
        "    all_motifs.append(NN_found_motifs)\n",
        "\n",
        "all_motifs = np.concatenate(all_motifs, axis=1)\n",
        "    "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DQU5fFv3044V"
      },
      "outputs": [],
      "source": [
        "# load JASPAR motifs\n",
        "from sklearn.linear_model import LinearRegression\n",
        "from sklearn.linear_model import Ridge\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "from sklearn.preprocessing import RobustScaler\n",
        "from scipy.stats import pearsonr\n",
        "import pickle\n",
        "import dill\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "\n",
        "CR_motif_name = []\n",
        "for i in range(0, 10):\n",
        "    for j in range(0, 64):\n",
        "        CR_motif_name.append(\"run_\" + str(i) + \"_motif_\" + str(j))\n",
        "\n",
        "JASPAR_fimo_table = pd.read_csv(\"/Your_Folder/all_JASPAR_fimo.csv\", index_col=0)\n",
        "JASPAR_fimo = JASPAR_fimo_table.values\n",
        "JASPAR_fimo_name = JASPAR_fimo_table.columns.values\n",
        "\n",
        "all_motifs_with_JASPAR = np.concatenate([all_motifs, JASPAR_fimo], axis=1)\n",
        "all_motif_name = []\n",
        "all_motif_name.extend(CR_motif_name)\n",
        "all_motif_name.extend(JASPAR_fimo_name)\n",
        "\n",
        "scaler = StandardScaler()\n",
        "scaler.fit(all_motifs_with_JASPAR)\n",
        "all_motifs_with_JASPAR_std = scaler.transform(all_motifs_with_JASPAR)\n",
        "\n",
        "print(all_motif_name)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "v5GI7xWTcG0g"
      },
      "outputs": [],
      "source": [
        "# feature selection with entropy\n",
        "from sklearn.feature_selection import RFE\n",
        "from sklearn.feature_selection import mutual_info_regression\n",
        "from sklearn.tree import DecisionTreeRegressor\n",
        "fs = RFE(estimator=DecisionTreeRegressor(), n_features_to_select=300, step=100, verbose=1)\n",
        "fs.fit(all_motifs_with_JASPAR_std, total_y)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "js16p3Z7SZTP"
      },
      "outputs": [],
      "source": [
        "# get selected features from feature selection\n",
        "information_selected_i = np.where(fs.support_)[0]\n",
        "information_score = fs.estimator_.feature_importances_"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "STa2vIYQtDo0"
      },
      "outputs": [],
      "source": [
        "# prediction accuracy check with selected\n",
        "selected_i = set(information_selected_i)\n",
        "selected_i = sorted(selected_i)\n",
        "print(selected_i)\n",
        "selected_motifs_with_JASPAR_std = all_motifs_with_JASPAR_std[:, selected_i]\n",
        "print(selected_motifs_with_JASPAR_std.shape)\n",
        "\n",
        "lr = LinearRegression()\n",
        "lr.fit(selected_motifs_with_JASPAR_std, total_y)\n",
        "print(lr.score(selected_motifs_with_JASPAR_std, total_y))"
      ]
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "background_execution": "on",
      "collapsed_sections": [],
      "name": "Motif_Finding_Contextual_Regression.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}