leave_one_out_val.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "af7f5051",
   "metadata": {},
   "source": [
    "Leave one out validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "883ab838",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import tensorflow\n",
    "from tensorflow import keras\n",
    "import matplotlib.pyplot as plt\n",
    "from Bio import SeqIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0441370f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def from_files_to_vectors(fasta_path, positive=True):\n",
    "    if positive:\n",
    "        proteins = list(SeqIO.parse(fasta_path+'pos_25.fasta', \"fasta\"))\n",
    "    else:\n",
    "        proteins = list(SeqIO.parse(fasta_path+'neg_25.fasta', \"fasta\"))\n",
    "    extension = \".out\"\n",
    "    files = [\"aac\", \"dpc\", \"ctdc\", \"ctdt\", \"ctdd\"]\n",
    "    if positive:\n",
    "        names = \"_pos\"\n",
    "    else:\n",
    "        names = \"_neg\"\n",
    "    for i in range(len(files)):\n",
    "        files[i] += names\n",
    "    datasets = [[] for el in files]\n",
    "    for i in range(len(files)):\n",
    "        with open(fasta_path+files[i]+extension) as f:\n",
    "            lines = f.readlines()[1:]\n",
    "            check_prot = 0\n",
    "            for line in lines:\n",
    "                information = line.split('\\t')\n",
    "                if not information[0] == proteins[check_prot].id:\n",
    "                    print(\"Error in protein order! Return\")\n",
    "                    return datasets\n",
    "                datasets[i].append(numpy.array([float(el) for el in information[1:]]))\n",
    "                check_prot += 1\n",
    "        datasets[i] = numpy.array(datasets[i])\n",
    "    return datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "271908b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "pos_datasets = from_files_to_vectors(\"./Adhesin_data/25_similarity/pos/\", positive=True)\n",
    "neg_datasets = from_files_to_vectors(\"./Adhesin_data/25_similarity/neg/\", positive=False)\n",
    "\n",
    "y_pos = numpy.ones(pos_datasets[0].shape[0])\n",
    "y_neg = numpy.zeros(neg_datasets[0].shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8eea4c1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# attach datasets in order to obtain a matrix of (n, 20+400+39+39+195) features\n",
    "\n",
    "# keep in mind the number of virulent factors and the number of not virulent factors\n",
    "rows = 0\n",
    "n_pos = y_pos.shape[0]\n",
    "n_neg = y_neg.shape[0]\n",
    "rows = n_pos + n_neg\n",
    "print('pos:', n_pos)\n",
    "print('neg:', n_neg)\n",
    "\n",
    "# feature vectors dimensions\n",
    "columns = 0\n",
    "for i in range(len(pos_datasets)):\n",
    "    tmp_dim = pos_datasets[i].shape[1]\n",
    "    print(i+1, '-th feature dim:', tmp_dim)\n",
    "    columns += tmp_dim\n",
    "\n",
    "# data matrix to process\n",
    "X = numpy.zeros((rows, columns))\n",
    "print('Data matrix dimension:', X.shape)\n",
    "for i in range(n_pos):\n",
    "    X[i] = numpy.concatenate([pos_datasets[j][i] for j in range(5)])\n",
    "for i in range(n_neg):\n",
    "    X[n_pos+i] = numpy.concatenate([neg_datasets[j][i] for j in range(5)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c57fbcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# permutation\n",
    "numpy.random.seed(990)\n",
    "\n",
    "y = numpy.concatenate((y_pos, y_neg), axis=0)\n",
    "c = numpy.random.permutation(numpy.arange(y.shape[0]))\n",
    "y = y[c]\n",
    "X = X[c] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0eb6f200",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import accuracy_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfc508c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remember: standardization or scaling AFTER train-test split, otherwise there's data leakage!!!\n",
    "X_train = X[:int(X.shape[0]*.5)]\n",
    "X_val = X[int(X.shape[0]*.5):int(X.shape[0]*.75)]\n",
    "X_test = X[int(X.shape[0]*.75):]\n",
    "\n",
    "y_train = y[:int(y.shape[0]*.5)]\n",
    "y_val = y[int(y.shape[0]*.5):int(y.shape[0]*.75)]\n",
    "y_test = y[int(y.shape[0]*.75):]\n",
    "\n",
    "print('Training data shape:', X_train.shape, y_train.shape)\n",
    "print('Validation data shape:', X_val.shape, y_val.shape)\n",
    "print('Test data shape:', X_test.shape, y_test.shape)\n",
    "\n",
    "print('\\nTraining virulent factors and not-virulent factors:', int(sum(y_train)), int(y_train.shape[0]-sum(y_train)))\n",
    "print('Validation virulent factors and not-virulent factors:', int(sum(y_val)), int(y_val.shape[0]-sum(y_val)))\n",
    "print('Test virulent factors and not-virulent factors:', int(sum(y_test)), int(y_test.shape[0]-sum(y_test)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02431db8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# PCA\n",
    "\n",
    "stdX = numpy.zeros(X_train.shape)\n",
    "stdX_val = numpy.zeros(X_val.shape)\n",
    "stdX_test = numpy.zeros(X_test.shape)\n",
    "\n",
    "means = numpy.zeros(X_train.shape[1])\n",
    "std_devs = numpy.zeros(X_train.shape[1])\n",
    "\n",
    "for j in range(X_train.shape[1]):\n",
    "    column = X_train[:,j]\n",
    "    means[j] = numpy.mean(column)\n",
    "    std_devs[j] = numpy.std(column)\n",
    "    stdX[:,j] = (column - means[j]) / std_devs[j]\n",
    "    stdX_val[:,j] = (X_val[:,j] - means[j]) / std_devs[j]\n",
    "    stdX_test[:,j] = (X_test[:,j] - means[j]) / std_devs[j]\n",
    "\n",
    "numpy.save('means', means)\n",
    "numpy.save('std_devs', std_devs)\n",
    "\n",
    "covariance_matrix = numpy.cov(stdX.T)\n",
    "\n",
    "eigen_values, eigen_vectors = numpy.linalg.eig(covariance_matrix)\n",
    "eigen_values = numpy.real(eigen_values)\n",
    "eigen_vectors = numpy.real(eigen_vectors)\n",
    "\n",
    "variance_explained = []\n",
    "for i in eigen_values:\n",
    "    variance_explained.append((i/sum(eigen_values))*100)\n",
    "\n",
    "cumulative_variance_explained = numpy.cumsum(variance_explained)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96d41cf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.title(\"Explained variance vs Number of components\")\n",
    "\n",
    "plt.xlabel(\"Number of components\")\n",
    "plt.ylabel(\"Cumulative explained variance\")\n",
    "plt.grid(color='gray', linewidth=.4)\n",
    "\n",
    "plt.plot(range(len(cumulative_variance_explained)), cumulative_variance_explained)\n",
    "\n",
    "plt.show()\n",
    "#plt.savefig('Explained variance.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f532275",
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 400 # more or less like adhesin\n",
    "print('Principal components:', K)\n",
    "print('Discarded components:', columns-K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ab24221",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.title(\"Explained variance vs first \" + str(K) + \" components\")\n",
    "\n",
    "plt.xlabel(\"Number of components\")\n",
    "plt.ylabel(\"Cumulative explained variance\")\n",
    "plt.grid(color='gray', linewidth=.4)\n",
    "\n",
    "plt.plot(range(K), cumulative_variance_explained[:K])\n",
    "\n",
    "plt.show()\n",
    "#plt.savefig('Explained variance.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "351a44e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Project using first K components\n",
    "\n",
    "projection_matrix = numpy.real((eigen_vectors.T[:][:K]).T)\n",
    "print(projection_matrix.shape)\n",
    "\n",
    "numpy.save('projection_matrix', projection_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c7d4fe5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# project matrices\n",
    "X_train = stdX.dot(projection_matrix)\n",
    "X_val = stdX_val.dot(projection_matrix)\n",
    "X_test = stdX_test.dot(projection_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "106a765d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorflow.keras import regularizers\n",
    "\n",
    "class neural_network:\n",
    "    def __init__(self):\n",
    "        input = tensorflow.keras.Input(shape=(K,))\n",
    "\n",
    "        dense = tensorflow.keras.layers.Dense(units=10, activation='sigmoid')(input) # it shows that a linear activation is sufficiet to get the same accuracy of the sigmoid (so it's like a matrix multiplication with a non-linearity at the end...)\n",
    "        #batch = tensorflow.keras.layers.BatchNormalization()(dense)\n",
    "        #activation = tensorflow.keras.activations.sigmoid(batch)\n",
    "        #drop = tensorflow.keras.layers.Dropout(.1)(activation)\n",
    "\n",
    "        output = tensorflow.keras.layers.Dense(1, activation='sigmoid')(dense)\n",
    "\n",
    "        model = tensorflow.keras.models.Model(inputs=input, outputs=output)\n",
    "        model.compile(optimizer=tensorflow.keras.optimizers.Adam(learning_rate=0.0001), loss='binary_crossentropy', metrics=['accuracy', 'mse', 'recall'])\n",
    "        self.model = model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ac169e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "nn = neural_network()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e89ce31",
   "metadata": {},
   "outputs": [],
   "source": [
    "history = nn.model.fit(\n",
    "    x=X_train, \n",
    "    y=y_train,\n",
    "    batch_size=64,\n",
    "    epochs=1000,\n",
    "    verbose=1,\n",
    "    validation_data=(X_val, y_val),\n",
    "    shuffle=True,\n",
    "    callbacks=[tensorflow.keras.callbacks.EarlyStopping(\n",
    "    restore_best_weights=True,\n",
    "    patience=20\n",
    "        )]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19bcf9f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.title('Accuracy during training')\n",
    "plt.plot(range(len(history.history['loss'])), history.history['accuracy'], label='training accuracy', color='green')\n",
    "plt.plot(range(len(history.history['loss'])), history.history['val_accuracy'], label='validation accuracy', color='red')\n",
    "plt.grid(color='gray', linewidth=.4)\n",
    "\n",
    "plt.legend(loc=\"lower right\")\n",
    "plt.xlabel(\"epochs\")\n",
    "plt.ylabel(\"accuracy\")\n",
    "plt.xlim(0, len(history.history['loss'])+1)\n",
    "plt.ylim(.45, 1.)\n",
    "plt.savefig('acc.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0dce97c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.title('Loss during training')\n",
    "plt.plot(range(len(history.history['loss'])), history.history['loss'], label='training loss', color='green')\n",
    "plt.plot(range(len(history.history['loss'])), history.history['val_loss'], label='validation loss', color='red')\n",
    "plt.grid(color='gray', linewidth=.4)\n",
    "\n",
    "plt.legend(loc=\"upper right\")\n",
    "plt.xlabel(\"epochs\")\n",
    "plt.ylabel(\"loss\")\n",
    "plt.savefig('loss.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f347a42d",
   "metadata": {},
   "outputs": [],
   "source": [
    "nn.model.evaluate(x=X_test, y=y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10166181",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate prediction distributions\n",
    "predictions = nn.model.predict(X_test)\n",
    "ground_truth = y_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d4c2cc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get true and false negative and positives\n",
    "results = numpy.c_[predictions, ground_truth]\n",
    "TP = results[numpy.where(((results[:,1] == 1) & (results[:,0] > .50)))]\n",
    "TN = results[numpy.where((results[:,1] == 0) & (results[:,0] < .50))]\n",
    "FP = results[numpy.where(((results[:,1] == 0) & (results[:,0] > .50)))]\n",
    "FN = results[numpy.where((results[:,1] == 1) & (results[:,0] < .50))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d94a12a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=[9,5])\n",
    "plt.hist(numpy.append(TP[:,0], TN[:,0]), bins=20, alpha=0.5, label = 'True positives and True negatives')\n",
    "plt.hist(numpy.append(FP[:,0], FN[:,0]), bins=20, alpha=0.5, label = 'False positives and False negatives')\n",
    "plt.xticks(numpy.arange(0, 1.05, step=0.05))\n",
    "plt.legend()\n",
    "plt.title('Virulent predictions on test set')\n",
    "plt.xlabel('Predicted value')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "603e77af",
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracy = (len(TP) + len(TN)) / len(results)\n",
    "sensitivity = len(TP) / (len(TP) + len(FN))\n",
    "precision = len(TP) / (len(TP) + len(FP))\n",
    "f1_score = 2 * (precision * sensitivity) / (precision + sensitivity)\n",
    "FPR = len(FP) / (len(FP) + len(TN))\n",
    "FNR = len(FN)/ (len(TP)+len(FN))\n",
    "TNR = len(TN)/(len(TN)+len(FP))\n",
    "FDR = len(FP)/(len(FP)+len(TP))\n",
    "TPR = len(TP)/(len(TP)+len(FN))\n",
    "NPV = len(TN)/(len(TN)+len(FN))\n",
    "\n",
    "print(\"Accuracy:\", accuracy)\n",
    "print(\"Sensitivity (Recall):\", sensitivity)\n",
    "print(\"Precision:\", precision)\n",
    "print(\"F1-score:\", f1_score)\n",
    "print(\"False Positive Rate (FPR):\", FPR)\n",
    "print(\"False Negative Rate (FNR):\", FNR)\n",
    "print(\"True Negative Rate (TNR):\", TNR)\n",
    "print(\"False Discovery Rate (FDR):\", FDR)\n",
    "print(\"True Positive Rate (FDR):\", TPR)\n",
    "print(\"Negative Predictive Value (NPV):\", NPV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9157e775",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "import numpy as np\n",
    "from sklearn import metrics\n",
    "from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score\n",
    "from sklearn.metrics import classification_report\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "y_pred_class = predictions >= 0.51\n",
    "cm = confusion_matrix(ground_truth, y_pred_class)\n",
    "TN, FP, FN, TP = cm.ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05f27fc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "cm.ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe065d0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "confusion_matrix = metrics.confusion_matrix(ground_truth, y_pred_class)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8526b42",
   "metadata": {},
   "outputs": [],
   "source": [
    "cm_display = metrics.ConfusionMatrixDisplay(confusion_matrix = confusion_matrix, display_labels = [\"Not Adhesin\", \"Adhesin\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "047d4148",
   "metadata": {},
   "outputs": [],
   "source": [
    "cm_display.plot(cmap='Blues')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0a0847a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Supponendo che cm_display.plot(cmap='Blues') sia già stato eseguito\n",
    "cm_display.plot(cmap='Blues')\n",
    "\n",
    "# Salva il grafico in un file PNG a 300 DPI\n",
    "plt.savefig('conf_matrix_adh.png', dpi=330)\n",
    "\n",
    "# Mostra il grafico\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}