Analyzer_peptides.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib as mt\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import re\n",
    "from itertools import groupby\n",
    "\n",
    "import MDAnalysis as mda\n",
    "from MDAnalysis.analysis import align, rms\n",
    "from MDAnalysis.tests.datafiles import PSF, DCD, PDB_small\n",
    "\n",
    "from statsmodels.nonparametric.smoothers_lowess import lowess\n",
    "import progressbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Time</th>\n",
       "      <th>Type</th>\n",
       "      <th>Protein</th>\n",
       "      <th>Peptide</th>\n",
       "      <th>AccType</th>\n",
       "      <th>DonType</th>\n",
       "      <th>WaterIdx</th>\n",
       "      <th>DistanceAWat</th>\n",
       "      <th>DistanceDWat</th>\n",
       "      <th>AngleDon</th>\n",
       "      <th>...</th>\n",
       "      <th>PosAtoms</th>\n",
       "      <th>PosAtomsIdx</th>\n",
       "      <th>ProtIsPos</th>\n",
       "      <th>RecRingType</th>\n",
       "      <th>LigRingType</th>\n",
       "      <th>RecRingAtoms</th>\n",
       "      <th>RecAtomsIdx</th>\n",
       "      <th>LigRingAtoms</th>\n",
       "      <th>LigRingAtomsIdx</th>\n",
       "      <th>Offset</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Frame</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>waterbridge</td>\n",
       "      <td>TYR473_R</td>\n",
       "      <td>GLU3_P</td>\n",
       "      <td>O3</td>\n",
       "      <td>O.co2</td>\n",
       "      <td>844.0</td>\n",
       "      <td>4.062278</td>\n",
       "      <td>2.598884</td>\n",
       "      <td>100.127071</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>hbond</td>\n",
       "      <td>ASN487_R</td>\n",
       "      <td>ASN4_P</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>waterbridge</td>\n",
       "      <td>LYS417_R</td>\n",
       "      <td>ARG6_P</td>\n",
       "      <td>Ng+</td>\n",
       "      <td>N3</td>\n",
       "      <td>801.0</td>\n",
       "      <td>2.696090</td>\n",
       "      <td>3.243224</td>\n",
       "      <td>137.087178</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>waterbridge</td>\n",
       "      <td>LYS417_R</td>\n",
       "      <td>ARG6_P</td>\n",
       "      <td>N3</td>\n",
       "      <td>Ng+</td>\n",
       "      <td>808.0</td>\n",
       "      <td>2.895963</td>\n",
       "      <td>3.821806</td>\n",
       "      <td>121.030379</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>waterbridge</td>\n",
       "      <td>TYR421_R</td>\n",
       "      <td>ARG6_P</td>\n",
       "      <td>O2</td>\n",
       "      <td>Ng+</td>\n",
       "      <td>808.0</td>\n",
       "      <td>2.845505</td>\n",
       "      <td>3.821806</td>\n",
       "      <td>121.030379</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 37 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       Time         Type   Protein Peptide AccType DonType  WaterIdx  \\\n",
       "Frame                                                                  \n",
       "0         0  waterbridge  TYR473_R  GLU3_P      O3   O.co2     844.0   \n",
       "0         0        hbond  ASN487_R  ASN4_P     NaN     NaN       NaN   \n",
       "0         0  waterbridge  LYS417_R  ARG6_P     Ng+      N3     801.0   \n",
       "0         0  waterbridge  LYS417_R  ARG6_P      N3     Ng+     808.0   \n",
       "0         0  waterbridge  TYR421_R  ARG6_P      O2     Ng+     808.0   \n",
       "\n",
       "       DistanceAWat  DistanceDWat    AngleDon  ...  PosAtoms  PosAtomsIdx  \\\n",
       "Frame                                          ...                          \n",
       "0          4.062278      2.598884  100.127071  ...       NaN          NaN   \n",
       "0               NaN           NaN         NaN  ...       NaN          NaN   \n",
       "0          2.696090      3.243224  137.087178  ...       NaN          NaN   \n",
       "0          2.895963      3.821806  121.030379  ...       NaN          NaN   \n",
       "0          2.845505      3.821806  121.030379  ...       NaN          NaN   \n",
       "\n",
       "      ProtIsPos  RecRingType LigRingType  RecRingAtoms  RecAtomsIdx  \\\n",
       "Frame                                                                 \n",
       "0           NaN          NaN         NaN           NaN          NaN   \n",
       "0           NaN          NaN         NaN           NaN          NaN   \n",
       "0           NaN          NaN         NaN           NaN          NaN   \n",
       "0           NaN          NaN         NaN           NaN          NaN   \n",
       "0           NaN          NaN         NaN           NaN          NaN   \n",
       "\n",
       "       LigRingAtoms  LigRingAtomsIdx Offset  \n",
       "Frame                                        \n",
       "0               NaN              NaN    NaN  \n",
       "0               NaN              NaN    NaN  \n",
       "0               NaN              NaN    NaN  \n",
       "0               NaN              NaN    NaN  \n",
       "0               NaN              NaN    NaN  \n",
       "\n",
       "[5 rows x 37 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table=pd.read_excel('Pep4/Results/Interactions_Table.xlsx',index_col=[0])\n",
    "table.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['saltbridge', 'hbond', 'hydroph_interaction', 'pistack', 'waterbridge']\n"
     ]
    }
   ],
   "source": [
    "interaction_types=list(set(table['Type']))\n",
    "print (interaction_types)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "time=table['Time'].drop_duplicates()\n",
    "index=table.index.drop_duplicates()\n",
    "hbond=[]\n",
    "saltbridge=[]\n",
    "pication=[]\n",
    "hydroph=[]\n",
    "pistack=[]\n",
    "wb=[]\n",
    "for x in index: \n",
    "    hbond.append(len([i for i in table.loc[x,'Type'] if i=='hbond']))\n",
    "    wb.append(len([i for i in table.loc[x,'Type'] if i=='waterbridge']))\n",
    "    saltbridge.append(len([i for i in table.loc[x,'Type'] if i=='saltbridge']))\n",
    "    pication.append(len([i for i in table.loc[x,'Type'] if i=='pication']))\n",
    "    hydroph.append(len([i for i in table.loc[x,'Type'] if i=='hydroph_interaction']))\n",
    "    pistack.append(len([i for i in table.loc[x,'Type'] if i=='pistack']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "plt.rcParams['axes.linewidth'] = 1.5\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,5))\n",
    "\n",
    "ax.scatter([i/100 for i in index],[float('nan') if x==0 else x for x in hbond],c='k',s=0.5, marker='.')\n",
    "#ax.scatter([i/1000 for i in time],wb,c='C4',s=1, marker='x',label='Water Bridge')\n",
    "\n",
    "plt.title ('Pep1',fontsize=24,fontweight='bold',color='k')\n",
    "plt.xlabel ('Time (ns)',fontsize=20,fontweight='bold')\n",
    "plt.ylabel ('Frequency\\n(π-Cation)',fontsize=20,fontweight='bold')\n",
    "\n",
    "plt.tick_params ('both',width=2,labelsize=14)\n",
    "\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "#plt.tight_layout()\n",
    "\n",
    "plt.savefig('Pep1_Hbond.png',dpi=300,format='png',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "table=table[table['Type']!='hydroph_interaction']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interaction_types=list(set(table['Type']))\n",
    "heatmap=pd.DataFrame(index=sorted(list(set(table['Protein']))), columns=sorted(set(table['Peptide'])))\n",
    "\n",
    "for i in list(set(table['Protein'])):\n",
    "    residue=sorted(list(table[table['Protein']==i]['Peptide']))\n",
    "    groups=groupby(residue)\n",
    "    for x in groups:\n",
    "        heatmap.loc[i,x[0]]=len(list(x[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap.columns=[i.split('_')[0] for i in heatmap.columns]\n",
    "heatmap.index=[i.split('_')[0] for i in heatmap.index]\n",
    "\n",
    "cols=list(heatmap.columns)\n",
    "cols.sort(key=lambda res: int(re.split('(\\d+)',res)[1]))\n",
    "ndx=list(heatmap.index)\n",
    "ndx.sort(key=lambda res: int(re.split('(\\d+)',res)[1]))\n",
    "heatmap = heatmap[cols]\n",
    "heatmap=heatmap.reindex(ndx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap=heatmap.transpose().fillna(0)\n",
    "heatmap = heatmap[(heatmap.T != 0).any()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,3))\n",
    "\n",
    "sns.heatmap (heatmap,yticklabels=heatmap.index, xticklabels=heatmap.columns,vmax=6500,cmap='Reds',cbar_kws=dict(label='Frequency',shrink=1,orientation='vertical',spacing='uniform',pad=0.02))\n",
    "\n",
    "plt.title('Peptide 6',size='18',weight='bold',color='hotpink')\n",
    "plt.ylabel('Peptide residue',fontsize=14,fontweight='bold')\n",
    "plt.xlabel('RBD residue',fontsize=14,fontweight='bold')\n",
    "plt.xticks (rotation=90,fontsize=5)\n",
    "plt.yticks (fontsize=8)\n",
    "\n",
    "#ax.xaxis.tick_top()\n",
    "plt.tick_params ('both',width=1.5)\n",
    "plt.savefig('Pep6/Results/Interactions_HM.png',dpi=300,format='png',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "inter=[]\n",
    "for i in table.index:\n",
    "    try:\n",
    "        inter.append([i,len(table.loc[i,'Time'])])\n",
    "    except:\n",
    "        pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "inter.sort(key=lambda x: x[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2, 30]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inter[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Universe with 49135 atoms>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u=mda.Universe ('Pep4/Results/equilibration.gro','Pep4/Results/production_fit.xtc')\n",
    "u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "for res in u.residues:\n",
    "    if res.resname=='TIP3':\n",
    "        res.resname='HOH'\n",
    "    if 'HI' in res.resname or 'HSD' in res.resname:\n",
    "        res.resname='HIS'\n",
    "    if 'CY' in res.resname:\n",
    "        res.resname='CYS'\n",
    "for atom in u.atoms:\n",
    "    if atom.name=='OH2':\n",
    "        atom.name='OW'\n",
    "\n",
    "pep_segment = u.add_Segment(segid='P')\n",
    "pep = u.select_atoms('resid 1:23')\n",
    "pep.residues.segments=pep_segment\n",
    "\n",
    "rec_segment = u.add_Segment(segid='R')\n",
    "rec=u.select_atoms('protein and (not segid P)')\n",
    "rec.residues.segments=rec_segment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/mnt/c/Users/angel/Linux_programs/miniconda3/envs/AnalysisMD/lib/python3.7/site-packages/MDAnalysis/coordinates/PDB.py:916: UserWarning: Found no information for attr: 'altLocs' Using default value of ' '\n",
      "  \"\".format(attrname, default))\n",
      "/mnt/c/Users/angel/Linux_programs/miniconda3/envs/AnalysisMD/lib/python3.7/site-packages/MDAnalysis/coordinates/PDB.py:916: UserWarning: Found no information for attr: 'icodes' Using default value of ' '\n",
      "  \"\".format(attrname, default))\n",
      "/mnt/c/Users/angel/Linux_programs/miniconda3/envs/AnalysisMD/lib/python3.7/site-packages/MDAnalysis/coordinates/PDB.py:916: UserWarning: Found no information for attr: 'occupancies' Using default value of '1.0'\n",
      "  \"\".format(attrname, default))\n",
      "/mnt/c/Users/angel/Linux_programs/miniconda3/envs/AnalysisMD/lib/python3.7/site-packages/MDAnalysis/coordinates/PDB.py:916: UserWarning: Found no information for attr: 'tempfactors' Using default value of '0.0'\n",
      "  \"\".format(attrname, default))\n"
     ]
    }
   ],
   "source": [
    "with mda.Writer(\"Pep4/Results/Frame.pdb\", u) as PDB:\n",
    "    for ts in u.trajectory[2:3]:\n",
    "        PDB.write(u.atoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MD analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMSD heatmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u=mda.Universe ('Spike_RBD/Results/equilibration.gro','Spike_RBD/Results/production_fit.xtc')\n",
    "u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pep_segment = u.add_Segment(segid='P')\n",
    "pep = u.select_atoms('resid 1:23')\n",
    "pep.residues.segments=pep_segment\n",
    "\n",
    "rec_segment = u.add_Segment(segid='R')\n",
    "rec=u.select_atoms('protein and (not segid P)')\n",
    "rec.residues.segments=rec_segment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Receptor=u.select_atoms('segid R and (resid 333:526)')\n",
    "print (Receptor)\n",
    "\n",
    "Peptide=u.select_atoms('segid P')\n",
    "print (Peptide)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bar=progressbar.ProgressBar(max_value=len(u.trajectory[:5000:100]))\n",
    "RMSD_hmap=pd.DataFrame()\n",
    "for i in range(len(u.trajectory[:5000:100])):\n",
    "    for j in range(len(u.trajectory[:5000:100])):\n",
    "        bb = Receptor.select_atoms('backbone')\n",
    "        u.trajectory[i]\n",
    "        A = bb.positions.copy() # coordinates of first frame\n",
    "        u.trajectory[j]         # forward to last frame\n",
    "        B = bb.positions.copy()  # coordinates of last frame\n",
    "        rmsd=rms.rmsd(A, B, center=True)\n",
    "        RMSD_hmap.loc[i,j]=rmsd\n",
    "    bar.update(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RMSD_hmap.to_excel('Spike_RBD/Results/RMSD_Hmap.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,8))\n",
    "\n",
    "n = 10\n",
    "while len(RMSD_hmap)/n > 6:\n",
    "    n += 10\n",
    "    \n",
    "ax=sns.heatmap(RMSD_hmap,square=True,xticklabels=n,yticklabels=n,vmin=0, vmax=2.2,cmap='RdBu_r',cbar_kws=dict(label='RMSD (Å)',shrink=1,orientation='vertical',spacing='uniform',pad=0.02))\n",
    "\n",
    "plt.title('apo-RBD',size=26,weight='bold',color='k')\n",
    "plt.ylabel('Time (ns)',fontsize=22,fontweight='bold', rotation=90)\n",
    "plt.xlabel('Time (ns)',fontsize=22,fontweight='bold', rotation=0)\n",
    "\n",
    "#ax.xaxis.tick_top()\n",
    "plt.tick_params ('both',width=2,labelsize=18)\n",
    "cax = plt.gcf().axes[-1]\n",
    "cax.tick_params(labelsize=16)\n",
    "ax.figure.axes[-1].yaxis.label.set_size(22)\n",
    "plt.tight_layout()\n",
    "plt.savefig('Spike_RBD/Results/RMSD_HM_RBD.png',quality=95,dpi=300,format='png',bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMSF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system=u.select_atoms('segid R and (resid 333:526)')\n",
    "calphas = system.select_atoms(\"name CA\")\n",
    "rmsfer = rms.RMSF(calphas, verbose=True).run(stop=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['axes.linewidth'] = 1.5\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.plot(calphas.resnums, rmsfer.rmsf,linewidth=1.5,color='k')\n",
    "plt.xlabel (' Residue αC',fontsize=16,fontweight='bold')\n",
    "plt.ylabel ('RMSF (Å)',fontsize=16,fontweight='bold')\n",
    "plt.tick_params ('both',width=2,labelsize=12)\n",
    "plt.grid (axis='y')\n",
    "plt.tight_layout()\n",
    "plt.savefig('Spike_RBD/Results/RBD_RMSF.png',dpi=600,format='png',transparent=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save=pd.DataFrame(rmsfer.rmsf,index=calphas.resnums)\n",
    "save.to_csv('Spike_RBD/Results/RBD_RMSF.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.add_TopologyAttr('tempfactors')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmsf=[]\n",
    "for atom in system.atoms:\n",
    "    rmsf.append(save.loc[atom.resid,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with mda.Writer(\"Spike_RBD/Results/RBD.pdb\", system.n_atoms) as PDB:\n",
    "    for ts in u.trajectory[0:1]:\n",
    "        system.atoms.tempfactors = rmsf\n",
    "        PDB.write(system.atoms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "array=[[min(save[0]), max(save[0])],[min(save[0]) , max(save[0])]]\n",
    "plt.imshow(array,cmap='jet')\n",
    "m=plt.colorbar(orientation='vertical',aspect=10)\n",
    "m.set_label ('RMSF (Å)',fontsize=20,fontweight='bold')\n",
    "\n",
    "m.ax.tick_params(labelsize=16) \n",
    "\n",
    "plt.savefig('Spike_RBD/Results/Bfactor_bar.png',quality=95,dpi=600,format='png',transparent=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system=u.select_atoms('protein')\n",
    "with mda.Writer(\"Spike_RBD/Results/RBD_snapshots.pdb\", system.n_atoms) as PDB:\n",
    "    for ts in u.trajectory[0:5000:1000]:\n",
    "        PDB.write(system.atoms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "AnalysisMD",
   "language": "python",
   "name": "analysismd"
  },
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}