From dafa895eefb26575aa6a43062868e42be02232fe Mon Sep 17 00:00:00 2001 From: paulshamrat Date: Sun, 11 Jun 2023 11:54:46 -0400 Subject: [PATCH] updated analysis of four trajectories --- codes/four-trajectories-230611.ipynb | 297 +++++++++++++++++++++++++++ 1 file changed, 297 insertions(+) create mode 100644 codes/four-trajectories-230611.ipynb diff --git a/codes/four-trajectories-230611.ipynb b/codes/four-trajectories-230611.ipynb new file mode 100644 index 0000000..fead393 --- /dev/null +++ b/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 +}