updated analysis of four trajectories

codes/four-trajectories-230611.ipynb

@@ -0,0 +1,297 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Installation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# NLP\n",
+    "! python -m pip install nltk==3.5\n",
+    "! python -m pip install numpy matplotlib\n",
+    "\n",
+    "# MDAnalysis\n",
+    "! pip install --upgrade MDAnalysis\n",
+    "! pip install --upgrade MDAnalysisTests\n",
+    "! pip install --upgrade MDAnalysisData\n",
+    "\n",
+    "# mdtraj, nglview, cython, pytraj, tsplot\n",
+    "# gnuplot, prody\n",
+    "! pip install mdtraj\n",
+    "! pip install nglview\n",
+    "! pip install cython --upgrade\n",
+    "! pip install pytraj\n",
+    "! pip install tsplot\n",
+    "\n",
+    "#!sudo apt-get install gnuplot-x11\n",
+    "#!pip install prody"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Mount Google Drive"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from google.colab import drive\n",
+    "drive.mount('/content/gdrive')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Let's make a folder first. We need to import the os and path library\n",
+    "import os\n",
+    "from pathlib import Path \n",
+    "\n",
+    "#Then, we define the path of the folder we want to create.\n",
+    "#Notice that the HOME folder for a hosted runtime in colab is /content/\n",
+    "mdpath = Path(\"/content/gdrive/MyDrive/works/psmb8/3unf\")\n",
+    "#mdpath = Path(\"/content\")\n",
+    "\n",
+    "\n",
+    "#Now, we create the folder using the os.mkdir() command\n",
+    "#The if conditional is just to check whether the folder already exists\n",
+    "#In which case, python returns an error\n",
+    "if os.path.exists(mdpath):\n",
+    "  print(\"path already exists\")\n",
+    "if not os.path.exists(mdpath):\n",
+    "  os.mkdir(mdpath)\n",
+    "  print(\"path was succesfully created\")\n",
+    "\n",
+    "# Change path\n",
+    "#First, we will change to the new folder. We will use python now :)\n",
+    "os.chdir(mdpath)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## RMSD"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import MDAnalysis as mda\n",
+    "from MDAnalysis.analysis import rms\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "# Load the trajectories and topologies\n",
+    "u1 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-wt-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-wt-20ns/md_1_all.xtc')\n",
+    "u2 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-q8c-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-q8c-20ns/md_1_all.xtc')\n",
+    "u3 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-r67c-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-r67c-20ns/md_1_all.xtc')\n",
+    "u4 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-n84c-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-n84c-20ns/md_1_all.xtc')\n",
+    "\n",
+    "# Select the protein atoms for RMSD calculation\n",
+    "ref1 = u1.select_atoms('protein')\n",
+    "ref2 = u2.select_atoms('protein')\n",
+    "ref3 = u3.select_atoms('protein')\n",
+    "ref4 = u4.select_atoms('protein')\n",
+    "\n",
+    "# Calculate the RMSD\n",
+    "R1 = rms.RMSD(u1, ref1, select='backbone', groupselections=['protein'])\n",
+    "R1.run()\n",
+    "R2 = rms.RMSD(u2, ref2, select='backbone', groupselections=['protein'])\n",
+    "R2.run()\n",
+    "R3 = rms.RMSD(u3, ref3, select='backbone', groupselections=['protein'])\n",
+    "R3.run()\n",
+    "R4 = rms.RMSD(u4, ref4, select='backbone', groupselections=['protein'])\n",
+    "R4.run()\n",
+    "\n",
+    "# Plot the RMSD\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(R1.rmsd[:, 1], R1.rmsd[:, 2], label='wild_type')\n",
+    "ax.plot(R2.rmsd[:, 1], R2.rmsd[:, 2], label='mutant1')\n",
+    "ax.plot(R3.rmsd[:, 1], R3.rmsd[:, 2], label='mutant2')\n",
+    "ax.plot(R4.rmsd[:, 1], R4.rmsd[:, 2], label='mutant3')\n",
+    "ax.legend()\n",
+    "ax.set_xlabel('Time (ps)')\n",
+    "ax.set_ylabel('RMSD (Å)')\n",
+    "\n",
+    "# Save the plot as a PNG file\n",
+    "plt.savefig('rmsd_plot.png', dpi=300)\n",
+    "\n",
+    "# Show the plot\n",
+    "plt.show()\n"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## RMSF CA"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import MDAnalysis as mda\n",
+    "from MDAnalysis.analysis import rms\n",
+    "from MDAnalysis.analysis.rms import RMSF\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "# Load the trajectories and topologies\n",
+    "u1 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-wt-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-wt-20ns/md_1_all.xtc')\n",
+    "u2 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-q8c-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-q8c-20ns/md_1_all.xtc')\n",
+    "u3 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-r67c-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-r67c-20ns/md_1_all.xtc')\n",
+    "u4 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-n84c-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-n84c-20ns/md_1_all.xtc')\n",
+    "\n",
+    "# Select the C-alpha atoms\n",
+    "calpha1 = u1.select_atoms('protein and name CA')\n",
+    "calpha2 = u2.select_atoms('protein and name CA')\n",
+    "calpha3 = u3.select_atoms('protein and name CA')\n",
+    "calpha4 = u4.select_atoms('protein and name CA')\n",
+    "\n",
+    "# Align the protein to the reference structure\n",
+    "ref1 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-wt-20ns/6wzu_solv_ions.pdb')\n",
+    "ref2 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-q8c-20ns/6wzu_solv_ions.pdb')\n",
+    "ref3 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-r67c-20ns/6wzu_solv_ions.pdb')\n",
+    "ref4 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-n84c-20ns/6wzu_solv_ions.pdb')\n",
+    "R1 = rms.RMSD(u1, ref1, select='protein and name CA', center=True, superposition=True)\n",
+    "R1.run()\n",
+    "R2 = rms.RMSD(u2, ref2, select='protein and name CA', center=True, superposition=True)\n",
+    "R2.run()\n",
+    "R3 = rms.RMSD(u3, ref3, select='protein and name CA', center=True, superposition=True)\n",
+    "R3.run()\n",
+    "R4 = rms.RMSD(u4, ref4, select='protein and name CA', center=True, superposition=True)\n",
+    "R4.run()\n",
+    "\n",
+    "# Calculate RMSF for each trajectory\n",
+    "RMSF1 = RMSF(calpha1).run()\n",
+    "RMSF2 = RMSF(calpha2).run()\n",
+    "RMSF3 = RMSF(calpha3).run()\n",
+    "RMSF4 = RMSF(calpha4).run()\n",
+    "\n",
+    "# Plot RMSF values of all trajectories on the same plot\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(RMSF1.rmsf, label='wild type')\n",
+    "ax.plot(RMSF2.rmsf, label='mutant1')\n",
+    "ax.plot(RMSF3.rmsf, label='mutant2')\n",
+    "ax.plot(RMSF4.rmsf, label='mutant3')\n",
+    "ax.set_xlabel('Residue')\n",
+    "ax.set_ylabel('RMSF (nm)')\n",
+    "ax.legend()\n",
+    "\n",
+    "# Set the y-axis limits based on the range of RMSF values\n",
+    "ymin = 0\n",
+    "ymax = max(RMSF1.rmsf.max(), RMSF2.rmsf.max(), RMSF3.rmsf.max(), RMSF4.rmsf.max()) + 0.1\n",
+    "ax.set_ylim(ymin, ymax)\n",
+    "\n",
+    "# Save the plot as a high-resolution PNG image\n",
+    "fig.savefig('rmsf_ca.png', dpi=300)\n",
+    "\n",
+    "plt.show()\n"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## RG"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import MDAnalysis as mda\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "# Load the trajectory and topology files for all systems\n",
+    "u1 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-wt-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-wt-20ns/md_1_all.xtc')\n",
+    "u2 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-q8c-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-q8c-20ns/md_1_all.xtc')\n",
+    "u3 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-r67c-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-r67c-20ns/md_1_all.xtc')\n",
+    "u4 = mda.Universe('/content/gdrive/MyDrive/works/bcl2/5n20-n84c-20ns/6wzu_solv_ions.pdb', '/content/gdrive/MyDrive/works/bcl2/5n20-n84c-20ns/md_1_all.xtc')\n",
+    "\n",
+    "# Select only protein atoms\n",
+    "protein_sel1 = u1.select_atoms('protein')\n",
+    "protein_sel2 = u2.select_atoms('protein')\n",
+    "protein_sel3 = u3.select_atoms('protein')\n",
+    "protein_sel4 = u4.select_atoms('protein')\n",
+    "\n",
+    "# Initialize arrays to store Rg values and time\n",
+    "Rg1 = np.zeros(len(u1.trajectory))\n",
+    "Rg2 = np.zeros(len(u2.trajectory))\n",
+    "Rg3 = np.zeros(len(u3.trajectory))\n",
+    "Rg4 = np.zeros(len(u4.trajectory))\n",
+    "time1 = np.zeros(len(u1.trajectory))\n",
+    "time2 = np.zeros(len(u2.trajectory))\n",
+    "time3 = np.zeros(len(u3.trajectory))\n",
+    "time4 = np.zeros(len(u4.trajectory))\n",
+    "\n",
+    "# Loop over all frames in trajectory and calculate Rg\n",
+    "for ts in u1.trajectory:\n",
+    "    Rg1[ts.frame] = protein_sel1.radius_of_gyration()\n",
+    "    time1[ts.frame] = u1.trajectory.time\n",
+    "\n",
+    "for ts in u2.trajectory:\n",
+    "    Rg2[ts.frame] = protein_sel2.radius_of_gyration()\n",
+    "    time2[ts.frame] = u2.trajectory.time\n",
+    "\n",
+    "for ts in u3.trajectory:\n",
+    "    Rg3[ts.frame] = protein_sel3.radius_of_gyration()\n",
+    "    time3[ts.frame] = u3.trajectory.time\n",
+    "\n",
+    "for ts in u4.trajectory:\n",
+    "    Rg4[ts.frame] = protein_sel4.radius_of_gyration()\n",
+    "    time4[ts.frame] = u4.trajectory.time\n",
+    "\n",
+    "# Plot Rg values of all systems on the same plot\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.plot(time1, Rg1, label='wild type')\n",
+    "ax.plot(time2, Rg2, label='mutant1')\n",
+    "ax.plot(time3, Rg3, label='mutant2')\n",
+    "ax.plot(time4, Rg4, label='mutant3')\n",
+    "ax.set_xlabel('Time (ps)')\n",
+    "ax.set_ylabel('Radius of gyration (nm)')\n",
+    "ax.legend()\n",
+    "\n",
+    "# Save the plot as a high-resolution PNG image\n",
+    "fig.savefig('rg.png', dpi=300)\n",
+    "\n",
+    "plt.show()\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "name": "python"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}