Activity3B_protein-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Molecular Dynamics Simulation of a Protein\n",
    "\n",
    "In this activity, we will replicate the first published simulation of a protein from 1977.\n",
    "\n",
    "*McCammon, J. A.; Gelin, B. R.; Karplus, M. Nature 1977, 267 (5612), 585–590.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mdtraj as md\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import nglview as ngl\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify the names of the parameter/topology file and trajectory file\n",
    "topology_file = \n",
    "trajectory_file = \n",
    "\n",
    "# Load trajectorry using mdtraj\n",
    "traj = md.load()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## View Trajectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "visualize = ngl.show_mdtraj(traj)\n",
    "visualize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can get the atom and bond information using mdtraj. This loads the bond and atom information\n",
    "# into a dataframe (from the Python library pandas). MolSSI has a pandas tutorial if you want to learn\n",
    "# more - https://education.molssi.org/python-data-analysis/02-pandas/index.html\n",
    "atoms, bonds = traj.topology.to_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can slice this the same way we sliced numpy arrays, but you have to add .iloc to the variable name.\n",
    "\n",
    "# Display information for the first 10 atoms\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We are going to analyze an N-C alpha bond in a phenylalanine \n",
    "# residue (Phe22) - http://www.bmrb.wisc.edu/referenc/commonaa.php?phe\n",
    "\n",
    "print(traj.topology.select('resid 21 and name N'))\n",
    "print(traj.topology.select('resid 21 and name CA'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "atoms.iloc[[349, 351]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting C alpha N bond"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bond_indices = [349, 351]\n",
    "NClength = md.compute_distances(traj, [bond_indices])\n",
    "\n",
    "t_start = 5000 # start at step 5000\n",
    "t_step = 500 # plot over 0.5 ps\n",
    "\n",
    "t_end = t_start + t_step\n",
    "\n",
    "plt.plot(NClength[t_start:t_end], color='Tomato')\n",
    "plt.title(r'N-C$\\alpha$ bond length')\n",
    "plt.xlabel('Time (fs)')\n",
    "plt.ylabel(r'Bond length (nm)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert the plot above to angstroms to match the plot from 1977"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make the same type of plot for the N-H bond for Phe22.\n",
    "# First you will need to identify the atom number for the H atom.\n",
    "# What is different about this plot?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting C alpha N C angle\n",
    "\n",
    "angle_indices = [349, 351, 367]\n",
    "bondangle = md.compute_angles(traj, [angle_indices])\n",
    "\n",
    "# add code to convert from radians to degrees\n",
    "\n",
    "plt.plot(bondangle[t_start:t_end], color='Goldenrod')\n",
    "plt.title(r'N-C$\\alpha$-C bond angle')\n",
    "plt.xlabel('Time (fs)')\n",
    "plt.ylabel('Bond angle (radians)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert the phi, psi, and omega torsions from radians to degrees and then plot them vs. time.\n",
    "# Fix the omega torsion plot to look more human-readable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}