diff --git a/Analyzer_peptides.ipynb b/Analyzer_peptides.ipynb
new file mode 100644
index 0000000..9a93547
--- /dev/null
+++ b/Analyzer_peptides.ipynb
@@ -0,0 +1,786 @@
+{
+ "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": [
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " Time | \n",
+ " Type | \n",
+ " Protein | \n",
+ " Peptide | \n",
+ " AccType | \n",
+ " DonType | \n",
+ " WaterIdx | \n",
+ " DistanceAWat | \n",
+ " DistanceDWat | \n",
+ " AngleDon | \n",
+ " ... | \n",
+ " PosAtoms | \n",
+ " PosAtomsIdx | \n",
+ " ProtIsPos | \n",
+ " RecRingType | \n",
+ " LigRingType | \n",
+ " RecRingAtoms | \n",
+ " RecAtomsIdx | \n",
+ " LigRingAtoms | \n",
+ " LigRingAtomsIdx | \n",
+ " Offset | \n",
+ "
\n",
+ " \n",
+ " Frame | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " 0 | \n",
+ " 0 | \n",
+ " waterbridge | \n",
+ " TYR473_R | \n",
+ " GLU3_P | \n",
+ " O3 | \n",
+ " O.co2 | \n",
+ " 844.0 | \n",
+ " 4.062278 | \n",
+ " 2.598884 | \n",
+ " 100.127071 | \n",
+ " ... | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ "
\n",
+ " \n",
+ " 0 | \n",
+ " 0 | \n",
+ " hbond | \n",
+ " ASN487_R | \n",
+ " ASN4_P | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " ... | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ "
\n",
+ " \n",
+ " 0 | \n",
+ " 0 | \n",
+ " waterbridge | \n",
+ " LYS417_R | \n",
+ " ARG6_P | \n",
+ " Ng+ | \n",
+ " N3 | \n",
+ " 801.0 | \n",
+ " 2.696090 | \n",
+ " 3.243224 | \n",
+ " 137.087178 | \n",
+ " ... | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ "
\n",
+ " \n",
+ " 0 | \n",
+ " 0 | \n",
+ " waterbridge | \n",
+ " LYS417_R | \n",
+ " ARG6_P | \n",
+ " N3 | \n",
+ " Ng+ | \n",
+ " 808.0 | \n",
+ " 2.895963 | \n",
+ " 3.821806 | \n",
+ " 121.030379 | \n",
+ " ... | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ "
\n",
+ " \n",
+ " 0 | \n",
+ " 0 | \n",
+ " waterbridge | \n",
+ " TYR421_R | \n",
+ " ARG6_P | \n",
+ " O2 | \n",
+ " Ng+ | \n",
+ " 808.0 | \n",
+ " 2.845505 | \n",
+ " 3.821806 | \n",
+ " 121.030379 | \n",
+ " ... | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ " NaN | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
5 rows × 37 columns
\n",
+ "
"
+ ],
+ "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": [
+ ""
+ ]
+ },
+ "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
+}
diff --git a/plipMD_peptide_V1.1.py b/plipMD_peptide_V1.1.py
new file mode 100644
index 0000000..8984781
--- /dev/null
+++ b/plipMD_peptide_V1.1.py
@@ -0,0 +1,194 @@
+import MDAnalysis as md
+
+import pandas as pd
+
+from prody import confProDy,parsePDB,writePDB
+silence_prody=confProDy(verbosity='none')
+
+from plip.structure.preparation import PDBComplex
+
+import progressbar
+
+import os
+import sys
+import warnings
+warnings.filterwarnings("ignore")
+
+
+def plipmd(topol=None,traj=None):
+
+ traj=list(traj.strip('[]').split(','))
+
+ u = md.Universe(topol,traj)
+
+ print ('\nINFO: your system contains {} segments with labels {} \n'.format(len(u.segments),list(u.segments.segids)))
+
+ if len (u.segments.segids) ==1:
+ print ('''
+ WARNING: Only one segment was identified. Proceeding to manual segment definition\n\n
+
+ You must define the peptide residues using MDAnalysis format.
+
+ For instance: resid 1:20
+
+ This means peptide is from residue 1 to 20. Only peptide residues are required.
+
+ Peptide residues will be added to "P" chain.
+
+ Extra residues will be consider "RECEPTOR" (R). \n\n
+ ''')
+
+ pep_residues=input('Type the peptide residues (example - resid 1:20 -):')
+ pep_segment = u.add_Segment(segid='P')
+ pep = u.select_atoms(pep_residues)
+ pep.residues.segments=pep_segment
+
+ rec_segment = u.add_Segment(segid='R')
+ rec=u.select_atoms('protein and (not segid P)')
+ rec.residues.segments=rec_segment
+
+ print ('\nINFO: your system contains {} segments with labels {} \n'.format(len(u.segments),list(u.segments.segids)))
+
+
+ peptide_segment=input ('\n\n1) Type the segment of your peptide:\n>')
+
+ receptor_segment=input ('\n\n1) Type the segment of your receptor:\n>')
+
+ sol_name=input ('\n2) Type the ResName of your water model. Must be 3-4 letter code (example - WAT or SOL or TIP3 -):\n>')
+
+ for res in u.residues:
+ if res.resname==sol_name:
+ res.resname='HOH'
+ if 'HI' in res.resname or 'HSD' in res.resname:
+ res.resname='HIS'
+ if 'CY' in res.resname:
+ res.resname='CYS'
+ for atom in u.atoms:
+ if atom.name=='OH2':
+ atom.name='OW'
+
+ System=u.select_atoms('protein or resname HOH',updating=True)
+ System=System.select_atoms('segid {} or (around 8 segid {})'.format(peptide_segment,peptide_segment),updating=True)
+
+ table=pd.DataFrame()
+ index=0
+ print ('\nINFO: Your trajectory lenght is:{} steps\n'.format(range(len(u.trajectory))))
+ start=int(input('4) Type the starting STEP to analyze:\n>'))
+ finish=int(input('\n5) Type the ending STEP to analyze:\n>'))
+ bar=progressbar.ProgressBar(max_value=finish)
+ print ('\n\n----- ----- ----- RUNNING THE ANALYSIS ----- ----- -----\n\n')
+
+ for i in range(start,finish):
+ try:
+ name='frame_tmp.pdb'
+ PDB= md.Writer(name, multiframe=False)
+ for ts in u.trajectory[i:i+1]:
+ PDB.write(System)
+
+ #Prody for converting peptide ATOM to HETATM
+ pdb_file = parsePDB(name)
+ labels=[True if i.getSegname()==peptide_segment else False for i in pdb_file]
+ pdb_file.setFlags('hetatm',labels)
+ writePDB(name,pdb_file)
+
+ plip_job = PDBComplex()
+ plip_job.load_pdb(name)
+ plip_job.analyze()
+
+ for key in plip_job.interaction_sets:
+ if key.split(':')[1]==peptide_segment:
+ for interaction in plip_job.interaction_sets.get(key).all_itypes:
+ if interaction.reschain != peptide_segment:
+ table.loc[index,'Frame']=ts.frame
+ table.loc[index,'Time']=ts.time
+ interaction_type=str(type(interaction)).split('.')[-1].replace("'>","")
+ table.loc[index,'Type']=interaction_type
+ table.loc[index,'Protein']=interaction.restype+str(interaction.resnr)+'_'+interaction.reschain
+ table.loc[index,'Peptide']=interaction.restype_l+str(interaction.resnr_l)+'_'+interaction.reschain_l
+
+ if interaction_type == 'hbond':
+ table.loc[index,'Acceptor']=interaction.atype
+ table.loc[index,'AcceptorIdx']=interaction.a.idx
+ table.loc[index,'Donor']=interaction.dtype
+ table.loc[index,'DonorIdx']=interaction.d.idx
+ table.loc[index,'DistanceAD']=interaction.distance_ad
+ table.loc[index,'DistanceAH']=interaction.distance_ah
+ table.loc[index,'Angle']=interaction.angle
+ table.loc[index,'Force']=interaction.type
+ table.loc[index,'ProtIsDon']=interaction.protisdon
+
+ elif interaction_type == 'pication':
+ table.loc[index,'Charge']=interaction.charge.type
+ table.loc[index,'ChargedAtoms']=",".join([i.type for i in interaction.charge.atoms])
+ table.loc[index,'Force']=interaction.type
+ table.loc[index,'RingType']=interaction.ring.type
+ table.loc[index,'RingAtoms']=",".join([i.type for i in interaction.ring.atoms])
+ table.loc[index,'RingAtomsIdx']=",".join([str(i.idx) for i in interaction.ring.atoms])
+
+ elif interaction_type=='saltbridge':
+ table.loc[index,'NegAtoms']=",".join([i.type for i in interaction.negative.atoms])
+ table.loc[index,'NegAtomsIdx']=",".join([str(i.idx) for i in interaction.negative.atoms])
+ table.loc[index,'PosAtoms']=",".join([i.type for i in interaction.positive.atoms])
+ table.loc[index,'PosAtomsIdx']=",".join([str(i.idx) for i in interaction.positive.atoms])
+ table.loc[index,'Distance']=interaction.distance
+ table.loc[index,'ProtIsPos']=interaction.protispos
+
+ elif interaction_type == 'hydroph_interaction':
+ table.loc[index,'RecAtom']=interaction.bsatom.type
+ table.loc[index,'RecAtomIdx']=interaction.bsatom.idx
+ table.loc[index,'LigAtom']=interaction.ligatom.type
+ table.loc[index,'LigAtomIdx']=interaction.ligatom.idx
+ table.loc[index,'Distance']=interaction.distance
+
+ elif interaction_type == 'pistack':
+ table.loc[index,'RecRingType']=interaction.proteinring.type
+ table.loc[index,'LigRingType']=interaction.ligandring.type
+ table.loc[index,'RecRingAtoms']=",".join([i.type for i in interaction.proteinring.atoms])
+ table.loc[index,'RecAtomsIdx']=",".join([str(i.idx) for i in interaction.proteinring.atoms])
+ table.loc[index,'LigRingAtoms']=",".join([i.type for i in interaction.ligandring.atoms])
+ table.loc[index,'LigRingAtomsIdx']=",".join([str(i.idx) for i in interaction.ligandring.atoms])
+ table.loc[index,'Distance']=interaction.distance
+ table.loc[index,'Angle']=interaction.angle
+ table.loc[index,'Offset']=interaction.offset
+
+ elif interaction_type == 'waterbridge':
+ table.loc[index,'AccType']=interaction.atype
+ table.loc[index,'DonType']=interaction.dtype
+ table.loc[index,'WaterIdx']=interaction.water_orig_idx
+ table.loc[index,'DistanceAWat']=interaction.distance_aw
+ table.loc[index,'DistanceDWat']=interaction.distance_dw
+ table.loc[index,'AngleDon']=interaction.d_angle
+ table.loc[index,'AngleWat']=interaction.w_angle
+ table.loc[index,'ProtIsDon']=interaction.protisdon
+
+ elif interaction_type == 'halogenbond':
+ table.loc[index,'Acceptor']=interaction.acctype
+ table.loc[index,'Donor']=interaction.acctype
+ table.loc[index,'Distance']=interaction.distance
+ table.loc[index,'DonAngle']=interaction.don_angle
+ table.loc[index,'AccAngle']=interaction.acc_angle
+
+ elif interaction_type=='metal_complex':
+ table.loc[index,'MetalType']=interaction.metal.type
+ table.loc[index,'Idx']=interaction.metal.idx
+ table.loc[index,'TargetType']=interaction.target_type
+ table.loc[index,'FunctGroup']=interaction.target.fgroup
+ table.loc[index,'Geometry']=interaction.geometry
+ table.loc[index,'Distance']=interaction.distance
+ table.loc[index,'Location']=interaction.location
+
+ index=index+1
+ bar.update(i+1)
+ os.remove(name)
+
+ except Exception:
+ continue
+
+ print ('\n\n----- ----- ----- SAVING THE RESULTS, PLEASE WAIT ----- ----- -----\n\n')
+ table.set_index(['Frame','Time'], inplace=True)
+ table.sort_index(inplace=True)
+ table.to_excel('Interactions_Table.xlsx')
+ print ('\n\n***** ***** ***** ALL DONE, DATA SAVED ON: Interactions_Table.xlsx ***** ***** *****\n\n')
+if __name__ == "__main__":
+
+ plipmd(sys.argv[1],sys.argv[2])