Keras_MLP_PT.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from Bio import SeqIO\n",
    "from keras.callbacks import EarlyStopping, ModelCheckpoint\n",
    "from keras.layers import Dense, Dropout\n",
    "from keras.models import Sequential, load_model\n",
    "from keras.utils import to_categorical\n",
    "\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def containsOnly(seq, alphabet=\"ACDEFGHIKLMNPQRSTVWY\"):\n",
    "    for char in seq:\n",
    "        if char not in alphabet: \n",
    "            return False\n",
    "    return True\n",
    "\n",
    "\n",
    "def encode_labels(df):\n",
    "    idx, labels = pd.factorize(df, sort=True)\n",
    "    encoded = to_categorical(idx)\n",
    "    return encoded, labels\n",
    "\n",
    "\n",
    "def encode_kmers(kmers, alphabet=\"ACDEFGHIKLMNPQRSTVWY\", k=35):\n",
    "    encoded_kmers = []\n",
    "    for kmer in kmers:\n",
    "        kmer = kmer[:k]    # Make sure to take expected length\n",
    "        try:\n",
    "            idx = [alphabet.index(aa) for aa in kmer]\n",
    "            encoded_kmer = to_categorical(idx, len(alphabet))\n",
    "            encoded_kmer = np.array(encoded_kmer.flatten(), dtype=np.int)\n",
    "            encoded_kmers.append(encoded_kmer)\n",
    "        except ValueError:  # in case of non_amino_acid letters\n",
    "            print(kmer)\n",
    "    return np.array(encoded_kmers)\n",
    "\n",
    "\n",
    "def encode_datasets_PLS(file, k=35, sep=\"\\t\", header=0):\n",
    "    df = pd.read_csv(file, sep=sep, header=header)\n",
    "    x = encode_kmers(df.motif, k=k)\n",
    "    y, labels = encode_labels(df.name)\n",
    "    return x, y, labels\n",
    "\n",
    "\n",
    "def create_model(xshape, yshape):\n",
    "    model = Sequential()\n",
    "    model.add(Dense(256, input_dim=xshape, activation='relu'))\n",
    "    model.add(Dropout(0.5))\n",
    "    model.add(Dense(128, activation='relu'))\n",
    "    model.add(Dropout(0.5))\n",
    "    model.add(Dense(64, activation='relu'))\n",
    "    model.add(Dense(yshape, activation='softmax'))\n",
    "\n",
    "    model.compile(loss='categorical_crossentropy',\n",
    "                  optimizer='rmsprop',\n",
    "                  metrics=['accuracy'])\n",
    "\n",
    "    return model\n",
    "\n",
    "\n",
    "def PPR_motifs(accession, sequence, model, labels, bg=\"B\", k=35):\n",
    "    kmers = [sequence[i:i + k] for i in range(len(sequence) - k + 1)]\n",
    "    encoded_kmers = encode_kmers(kmers, k=k)\n",
    "    y_probs = model.predict(encoded_kmers)\n",
    "    y_classes = y_probs.argmax(axis=-1)\n",
    "\n",
    "    starts = np.where(labels[y_classes] != bg)\n",
    "    cls = y_classes[starts]\n",
    "    proba = [y_probs[i][y_classes[i]] for i in starts[0]]\n",
    "    motif = [sequence[s:s + k] for s in starts[0]]\n",
    "    d = {\"accession\": accession,\n",
    "         \"start\": starts[0],\n",
    "         \"end\": starts[0] + k,\n",
    "         \"name\": labels[cls],\n",
    "         \"score\": proba,\n",
    "         \"strand\": \"+\",\n",
    "         \"motif\": motif}\n",
    "\n",
    "    df = pd.DataFrame(d,\n",
    "                      columns=['accession', 'start', 'end', 'name', 'score', 'strand', 'motif'])\n",
    "\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline\n",
    "\n",
    "def plot_cm(model, x, y, lables, bb=None):\n",
    "    y_p = model.predict(x)\n",
    "    y_p_classes = y_p.argmax(axis=-1)\n",
    "    y_p_lables = np.array(labels[y_p_classes])\n",
    "    \n",
    "    y_classes = y.argmax(axis=-1)\n",
    "    y_lable = np.array(labels[y_classes])\n",
    "    cmatrix = pd.crosstab(y_lable, y_p_lables, \n",
    "                          rownames=[\"curatedLabel\"], \n",
    "                          colnames=[\"Prediction\"])\n",
    "    if bb:\n",
    "        cmatrix['B']['B'] = bb\n",
    "    sns.heatmap(cmatrix,\n",
    "                annot=True, fmt='',\n",
    "                cmap=plt.cm.Blues)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, labels = encode_datasets_PLS(\"dataset/new/Ath_167_19k.tsv\", k=35)\n",
    "\n",
    "x_train, x_valid, y_train, y_valid = train_test_split(x, y, test_size=0.25, random_state=42)\n",
    "\n",
    "assert x_train.shape[0] == y_train.shape[0]\n",
    "assert x_valid.shape[0] == y_valid.shape[0]\n",
    "\n",
    "xshape, yshape = x_train.shape[1], y_train.shape[1]\n",
    "\n",
    "print(\"x_train.shape: \\t{}\".format(x_train.shape))\n",
    "print(\"y_train.shape: \\t{}\".format(y_train.shape))\n",
    "print(\"x_valid.shape: \\t{}\".format(x_valid.shape))\n",
    "print(\"y_valid.shape: \\t{}\".format(y_valid.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = \"./Model/Model_Epoch_{epoch:02d}.h5\"\n",
    "\n",
    "if not os.path.exists(os.path.dirname(weights)):\n",
    "    os.mkdir(os.path.dirname(weights))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Train model\n",
    "model = create_model(xshape, yshape)\n",
    "\n",
    "check_pointer = ModelCheckpoint(filepath=weights, \n",
    "                                verbose=0, \n",
    "                                save_best_only=False)\n",
    "\n",
    "early_stopping = EarlyStopping(monitor='loss',\n",
    "                               min_delta=1e-5, \n",
    "                               patience=10, \n",
    "                               verbose=1, \n",
    "                               mode='auto')\n",
    "        \n",
    "model.fit(x_train, y_train, \n",
    "          batch_size=128,\n",
    "          epochs=50, \n",
    "          validation_data=(x_valid, y_valid),\n",
    "          callbacks=[early_stopping, check_pointer], \n",
    "          shuffle=True,\n",
    "          verbose=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Load Trained model\n",
    "model_bs0 = load_model(\"./Model/Model_Epoch_44.h5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "score, acc = model_bs0.evaluate(x, y, verbose=0)\n",
    "print('Test loss: {}'.format(score))\n",
    "print('Test accuracy: {}'.format(acc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_cm(model_bs0, x, y, labels, bb=3000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "\n",
    "def motifs_ath(model, prefix):\n",
    "    P = prefix + \"_P.tsv\"\n",
    "    L = prefix + \"_L.tsv\"\n",
    "    T = prefix + \"_T.tsv\"\n",
    "    U = prefix + \"_U.tsv\"\n",
    "    N = prefix + \"_N.tsv\"\n",
    "    fasta=\"dataset/ath_167/Athaliana_167_TAIR10.protein.fa\"\n",
    "    \n",
    "    with open(\"dataset/ath_167/p_ath_ppr.txt\") as f1:\n",
    "        ath_P = [line.rstrip() for line in f1]\n",
    "\n",
    "    with open(\"dataset/ath_167/pls_ath_ppr.txt\") as f2:\n",
    "        ath_L = [line.rstrip() for line in f2]\n",
    "\n",
    "    with open(\"dataset/ath_167/t_ath_tpr.txt\") as f3:\n",
    "        ath_T = [line.rstrip() for line in f3]\n",
    "\n",
    "    with open(\"dataset/ath_167/pt_ath.txt\") as f4:\n",
    "        ath_U = [line.rstrip() for line in f4]\n",
    "\n",
    "    for record in SeqIO.parse(fasta, \"fasta\"):\n",
    "        sequence = str(record.seq)[:-1] # remove \"*\" at the end of every sequence\n",
    "        accession = record.id\n",
    "    \n",
    "        if len(sequence) < 35 or not containsOnly(sequence):\n",
    "            continue\n",
    "        #ATG = r\"(AT.?G\\d{5})\"\n",
    "        ATG = r\"(AT.?G\\d{5}\\.1)\"\n",
    "        if re.search(ATG, accession):\n",
    "            match = re.search(ATG, accession)\n",
    "            acc = match.group()\n",
    "            try:\n",
    "                df = PPR_motifs(accession, sequence, model, labels, bg=\"B\", k=35)\n",
    "            except ValueError:\n",
    "                continue\n",
    "        else:\n",
    "            continue\n",
    "        if acc[:9] in ath_P:\n",
    "            df.to_csv(P, mode='a', sep=\"\\t\", header=False, index=False)\n",
    "        elif acc[:9] in ath_L:\n",
    "            df.to_csv(L, mode='a', sep=\"\\t\", header=False, index=False)\n",
    "        elif acc[:9] in ath_T:\n",
    "            df.to_csv(T, mode='a', sep=\"\\t\", header=False, index=False)\n",
    "        elif acc[:9] in ath_U:\n",
    "            df.to_csv(U, mode='a', sep=\"\\t\", header=False, index=False)\n",
    "        else:\n",
    "            df.to_csv(N, mode='a', sep=\"\\t\", header=False, index=False)\n",
    "    return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "motifs_ath(model_bs0, prefix=\"dataset/ath_167/000\")"
   ]
  }
 ],
 "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
}