PSA_PCA_ANOVA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "outputs": [],
   "source": [
    "import os\n",
    "import MDAnalysis as mda\n",
    "from MDAnalysis.analysis import psa\n",
    "import plotly.graph_objects as go\n",
    "import mdtraj as md\n",
    "from tqdm import tqdm"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-08-14T01:11:25.426433800Z",
     "start_time": "2023-08-14T01:11:24.805349800Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [],
   "source": [
    "trajectory_file = 'I:/Nick/200ns/100ns/2/md_0_1.xtc'\n",
    "topology_file = 'I:/Nick/200ns/100ns/2/EKB1m.pdb'\n",
    "output_dir = 'I:/Nick/200ns/100ns/2'\n",
    "chunk_size_ns = 10"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-08-14T01:26:40.373173900Z",
     "start_time": "2023-08-14T01:26:40.366663400Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [
    {
     "ename": "OSError",
     "evalue": "No such file: I:/Nick/200ns/100ns/2/md_0_1.xtc",
     "output_type": "error",
     "traceback": [
      "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m",
      "\u001B[1;31mOSError\u001B[0m                                   Traceback (most recent call last)",
      "Cell \u001B[1;32mIn[9], line 1\u001B[0m\n\u001B[1;32m----> 1\u001B[0m traj \u001B[38;5;241m=\u001B[39m \u001B[43mmd\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mload\u001B[49m\u001B[43m(\u001B[49m\u001B[43mtrajectory_file\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mtop\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mtopology_file\u001B[49m\u001B[43m)\u001B[49m\n\u001B[0;32m      2\u001B[0m total_frames \u001B[38;5;241m=\u001B[39m traj\u001B[38;5;241m.\u001B[39mn_frames\n\u001B[0;32m      3\u001B[0m frame_rate \u001B[38;5;241m=\u001B[39m traj\u001B[38;5;241m.\u001B[39mtimestep\n",
      "File \u001B[1;32m~\\anaconda3\\envs\\MD-things\\lib\\site-packages\\mdtraj\\core\\trajectory.py:409\u001B[0m, in \u001B[0;36mload\u001B[1;34m(filename_or_filenames, discard_overlapping_frames, **kwargs)\u001B[0m\n\u001B[0;32m    403\u001B[0m     \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mIOError\u001B[39;00m(\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mSorry, no loader for filename=\u001B[39m\u001B[38;5;132;01m%s\u001B[39;00m\u001B[38;5;124m (extension=\u001B[39m\u001B[38;5;132;01m%s\u001B[39;00m\u001B[38;5;124m) \u001B[39m\u001B[38;5;124m'\u001B[39m\n\u001B[0;32m    404\u001B[0m                 \u001B[38;5;124m'\u001B[39m\u001B[38;5;124mwas found. I can only load files \u001B[39m\u001B[38;5;124m'\u001B[39m\n\u001B[0;32m    405\u001B[0m                 \u001B[38;5;124m'\u001B[39m\u001B[38;5;124mwith extensions in \u001B[39m\u001B[38;5;132;01m%s\u001B[39;00m\u001B[38;5;124m'\u001B[39m \u001B[38;5;241m%\u001B[39m (\n\u001B[0;32m    406\u001B[0m                     filename_or_filenames[\u001B[38;5;241m0\u001B[39m], extension, FormatRegistry\u001B[38;5;241m.\u001B[39mloaders\u001B[38;5;241m.\u001B[39mkeys()))\n\u001B[0;32m    408\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m loader\u001B[38;5;241m.\u001B[39m\u001B[38;5;18m__name__\u001B[39m \u001B[38;5;129;01mnot\u001B[39;00m \u001B[38;5;129;01min\u001B[39;00m [\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mload_dtr\u001B[39m\u001B[38;5;124m'\u001B[39m]:\n\u001B[1;32m--> 409\u001B[0m         \u001B[43m_assert_files_exist\u001B[49m\u001B[43m(\u001B[49m\u001B[43mfilename_or_filenames\u001B[49m\u001B[43m)\u001B[49m\n\u001B[0;32m    410\u001B[0m \u001B[38;5;28;01melse\u001B[39;00m:\n\u001B[0;32m    411\u001B[0m     _assert_files_or_dirs_exist(filename_or_filenames)\n",
      "File \u001B[1;32m~\\anaconda3\\envs\\MD-things\\lib\\site-packages\\mdtraj\\core\\trajectory.py:100\u001B[0m, in \u001B[0;36m_assert_files_exist\u001B[1;34m(filenames)\u001B[0m\n\u001B[0;32m     98\u001B[0m \u001B[38;5;28;01mfor\u001B[39;00m fn \u001B[38;5;129;01min\u001B[39;00m filenames:\n\u001B[0;32m     99\u001B[0m     \u001B[38;5;28;01mif\u001B[39;00m \u001B[38;5;129;01mnot\u001B[39;00m (os\u001B[38;5;241m.\u001B[39mpath\u001B[38;5;241m.\u001B[39mexists(fn) \u001B[38;5;129;01mand\u001B[39;00m os\u001B[38;5;241m.\u001B[39mpath\u001B[38;5;241m.\u001B[39misfile(fn)):\n\u001B[1;32m--> 100\u001B[0m         \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mIOError\u001B[39;00m(\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mNo such file: \u001B[39m\u001B[38;5;132;01m%s\u001B[39;00m\u001B[38;5;124m'\u001B[39m \u001B[38;5;241m%\u001B[39m fn)\n",
      "\u001B[1;31mOSError\u001B[0m: No such file: I:/Nick/200ns/100ns/2/md_0_1.xtc"
     ]
    }
   ],
   "source": [
    "traj = md.load(trajectory_file, top=topology_file)\n",
    "total_frames = traj.n_frames\n",
    "frame_rate = traj.timestep"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-08-14T01:26:42.041510200Z",
     "start_time": "2023-08-14T01:26:41.959939100Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "chunk_size_frames = int(chunk_size_ns / frame_rate)\n",
    "\n",
    "for chunk_start in range(0, total_frames, chunk_size_frames):\n",
    "    chunk_end = min(chunk_start + chunk_size_frames, total_frames)\n",
    "    chunk = traj[chunk_start:chunk_end]\n",
    "    output_file = f'{output_dir}/chunk_{chunk_start}-{chunk_end}.xtc'\n",
    "    chunk.save(output_file)\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import plotly.subplots as sp\n",
    "# ANOVA HEATMAP\n",
    "fig = sp.make_subplots(rows=5, cols=4, subplot_titles=[f\"df{i}\" for i in range(1,21)],\n",
    "                       shared_xaxes=True, shared_yaxes=True, vertical_spacing=0.03, horizontal_spacing=0.03,\n",
    "                       column_width=[0.8]*4, row_heights=[0.2, 0.2, 0.2, 0.2, 0.2], specs=[[{'type': 'heatmap'}]*4]*5,\n",
    "                       print_grid=False, )\n",
    "fig.update_layout(height=1000, width=1200, title='Significance of differences in means', font=dict(size=10))\n",
    "\n",
    "for i in range(1, 21):\n",
    "    results = []\n",
    "    for j in range(1, 21):\n",
    "        row = []\n",
    "        for c in range(0, 15):\n",
    "            f_value, p_value = stats.f_oneway(dfs[j-1].iloc[:, c], dfs[i-1].iloc[:, c])\n",
    "            if p_value < 0.05:\n",
    "                row.append(1)\n",
    "            else:\n",
    "                row.append(0)\n",
    "        results.append(row)\n",
    "    results = np.array(results)\n",
    "    percentage = 100 * np.count_nonzero(results) / results.size\n",
    "    print(f\"Percentage of significant's over the total for df{i}: {percentage:.2f}%\")\n",
    "\n",
    "    heatmap = go.Heatmap(z=results, x=[i for i in range(1,16)], y=[(j)*10 for j in range(1,21)], colorscale='Viridis')\n",
    "    row = (i-1) // 4 + 1\n",
    "    col = (i-1) % 4 + 1\n",
    "    fig.add_trace(heatmap, row=row, col=col)\n",
    "\n",
    "fig.update_xaxes(title_text='Variant', row=5, col=3, title_font=dict(size=20), title_standoff=2.5)\n",
    "fig.update_yaxes(title_text='Duration (ns)', row=3, col=1, title_font=dict(size=20))\n",
    "\n",
    "fig.show()\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "outputs": [],
   "source": [
    "def run_psa(folder_path):\n",
    "    # Get a list of all .xtc files in the folder\n",
    "    traj_files = [f for f in os.listdir(folder_path) if f.endswith('.xtc')]\n",
    "\n",
    "    # Find the PDB file in the same directory as the trajectory files\n",
    "    pdb_file = None\n",
    "    for f in os.listdir(folder_path):\n",
    "        if f.endswith('.pdb'):\n",
    "            pdb_file = os.path.join(folder_path, f)\n",
    "            break\n",
    "\n",
    "    if pdb_file is None:\n",
    "        print(f'No PDB file found in {folder_path}')\n",
    "        return\n",
    "\n",
    "    # Create a list of universe objects for each trajectory file\n",
    "    trajectories = []\n",
    "    labels = []\n",
    "    for traj_file in traj_files:\n",
    "        traj_path = os.path.join(folder_path, traj_file)\n",
    "        label = os.path.splitext(traj_file)[0].split('_')[-1]\n",
    "        labels.append(label)\n",
    "        u = mda.Universe(pdb_file, traj_path)\n",
    "        trajectories.append(u)\n",
    "\n",
    "    # Create a reference universe object from the PDB file\n",
    "    ref = mda.Universe(pdb_file)\n",
    "\n",
    "    # Initialize a PSA object with the trajectories, labels, and reference\n",
    "    ps = psa.PSAnalysis(trajectories,\n",
    "                        labels=labels,\n",
    "                        reference=ref,\n",
    "                        select='name CA',\n",
    "                        path_select='name CA')\n",
    "\n",
    "    # Run the Hausdorff distance metric\n",
    "    ps.run(metric='hausdorff')\n",
    "\n",
    "    # Get the array map from the PSA object\n",
    "    array_map = ps.get_output()\n",
    "\n",
    "    # Create the figure and add the heatmap trace\n",
    "    fig = go.Figure(data=go.Heatmap(z=np.array(array_map), colorscale='RdYlBu', showscale=True, reversescale=False))\n",
    "\n",
    "    # Set the x-axis and y-axis labels\n",
    "    fig.update_layout(xaxis_title='Simulation Chunk', yaxis_title='Simulation Chunk')\n",
    "\n",
    "    # Save the plot with the relevant name\n",
    "    plot_name = f'Hausdorff_{os.path.basename(folder_path)}'\n",
    "    fig.write_html(f'{plot_name}.html')\n",
    "\n",
    "def run_psa_on_all_directories(root_folder):\n",
    "    for dirpath, dirnames, filenames in os.walk(root_folder):\n",
    "        for dirname in dirnames:\n",
    "            folder_path = os.path.join(dirpath, dirname)\n",
    "            run_psa(folder_path)\n"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-08-14T00:45:36.527499600Z",
     "start_time": "2023-08-14T00:45:36.525999100Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "outputs": [],
   "source": [
    "root_folder = '/mnt/f/200ns/100ns.2'\n",
    "run_psa_on_all_directories(root_folder)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-08-14T00:45:38.090113300Z",
     "start_time": "2023-08-14T00:45:38.075610500Z"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}