<jupyter_start><jupyter_text>Lab.08 / IBM3202 – Trajectory Analysis using MDanalysis Theoretical Aspects Now that you have already generated a molecular dynamics trajectory in the previous tutorial, it is crucial to obtain quantifiable insights about your molecular system. There are a handfull of metrics that can be employed to achieve this, here we are going to focus on two of the most popular ones. RMSD and RMSF and distances. Due to time constrains we are not going to cover more advanced metrics but they will be available as an appendix.  FIGURE 1. MD simulations of glycosylated SARS-CoV-2 spike protein attached to a membrane. Taken from the Amaro Lab at UCSD. ⚠️⚠️ The following section is an adapted excerpt from the introduction of the Cpptraj tutorial by <NAME>, available in the [this link](http://ambermd.org/tutorials/analysis/tutorial1/). ***Root Mean Square Deviation (RMSD)*** overview$RMSD$ measures the deviation of a target set of coordinates (i.e. a structure) to a reference set of coordinates, with $RMSD=0$ indicating a perfect overlap. Then it follows that if we have a MD trajectory one would expect that the lower the RMSD the less changes happen in the time scale studied. RMSD is defined as:$RMSD = \sqrt{\frac{\sum_{i = 0}^N m_i(X_i - Y_i)^2}{M}}$ Where **N** is the number of atoms,  $m_{i}$ is the mass of atom $i$, $X_i$ is the coordinate vector for target atom $i$, $Y_i$ is the coordinate vector for reference atom $i$, and $M$ is the total mass. If the $RMSD$ is not mass-weighted, for all $i$, $m_i = 1$, and $M = N$.When calculating $RMSD$ of a target to reference structure, there are two very important requirements as we will see soon in the practical part of this tutorial:1. The number of atoms in the target must match the number of atoms in the reference.2. The ordering of atoms in the target must match the ordering of atoms in the reference.  ***Root Mean Square Fluctuation (RMSF)*** overview As mentioned in the [MDanalysis user guide ](https://userguide.mdanalysis.org/stable/examples/analysis/alignment_and_rms/rmsf.html):> The root-mean-square-fluctuation ($RMSF$) of a structure is **the time average of the RMSD**. It is calculated according to the below equation, where $x_i$ is the coordinates of particle $i$ and $⟨x_i⟩$ is the ensemble average position of $i$: $ ρ_i=\sqrt{⟨(x_i−⟨x_i⟩)^2⟩}$> Where the $RMSD$ quantifies how much a structure diverges from a reference over time, the **$RSMF$ can reveal which areas of the system are the most mobile**. While $RMSD$ is frequently calculated to an initial state, the RMSF should be calculated to an average structure of the simulation. An area of the structure with high $RMSF$ values frequently diverges from the average, indicating high mobility. When $RMSF$ analysis is carried out on proteins, it is typically restricted to backbone or alpha-carbon atoms; these are more characteristic of conformational changes than the more flexible side-chains. FIGURE 2. RMSD and RMSF plots of the structural changes occuring due to antagonist binding in the ligand binding pocket of androgen receptor, elucidated through MD simulations.Sugunadevi S et al (2018) Front Pharmacology 9, 492   DistancesAs you may recall trajectory files stores the position of each individual atoms, the calculation of distances through an MD simulation is usually really straightforward.  MDanalysis package overview As defined in the [documentation](https://docs.mdanalysis.org/stable/documentation_pages/overview.html): > **MDAnalysis** is a Python package that provides classes to access data in molecular dynamics trajectories. It is object oriented so it treats atoms, groups of atoms, trajectories, etc as different objects. Each object has a number of operations defined on itself (also known as “methods”) and also contains values describing the object (“attributes”). For example, a **AtomGroup** object has a **center_of_mass()** method that returns the center of mass of the group of atoms. It also contains an attribute called residues that lists all the residues that belong to the group. Using methods such as **select_atoms()** (which uses CHARMM-style atom Selection commands) one can create new objects (in this case, another **AtomGroup**). Example of MDanalysis code A typical usage pattern is to iterate through a trajectory and analyze coordinates for every frame. In the following example the end-to-end distance of a protein and the radius of gyration of the backbone atoms are calculated: ```!pip3 install --upgrade MDAnalysis!pip install --upgrade MDAnalysisTestsimport MDAnalysisfrom MDAnalysis.tests.datafiles import PSF,DCD   test trajectoryimport numpy.linalgu = MDAnalysis.Universe(PSF,DCD)   always start with a Universenterm = u.select_atoms('segid 4AKE and name N')[0]   can access structure via segid (s4AKE) and atom namecterm = u.select_atoms('segid 4AKE and name C')[-1]   ... takes the last atom named 'C'bb = u.select_atoms('protein and backbone')   a selection (a AtomGroup)for ts in u.trajectory:   iterate through all frames    r = cterm.pos - nterm.pos   end-to-end vector from atom positions    d = numpy.linalg.norm(r)    end-to-end distance    rgyr = bb.radius_of_gyration()   method of a AtomGroup; updates with each frame    print "frame = %d: d = %f Angstroem, Rgyr = %f Angstroem" % (ts.frame, d, rgyr)```  Basic concepts of MD analysis1.   Universes and atom groups2.   Selections **Universe and AtomGroup**MDAnalysis is object oriented. Molecular systems consist of Atom objects (instances of the class MDAnalysis.core.groups.Atom), which are grouped in AtomGroup instances. You build the AtomGroup of your system by loading a topology (list of atoms and possibly their connectivity) together with a trajectory (coordinate information) into the central data structure, the Universe object:```u = MDAnalysis.Universe(PSF, DCD)print(u)``` **Selections**MDAnalysis comes with a fairly complete atom selection facility. Primarily, one uses the method select_atoms() of a Universe:```>>> CA = u.select_atoms("protein and name CA")>>> CA>>> ```but really any AtomGroup has a select_atoms() method:```>>> acidic = CA.select_atoms("resname ASP or resname GLU")>>> acidic>>> list(acidic.residues)[, , , ... ] ```See also All the selection keywords are described in the documentation.Numerical ranges can be written as first-last (or equivalently first:last 1), where the range is inclusive. For example, get residues with residue IDs 5 to 100:```>>> u.select_atoms("resid 5-100")>>> u.select_atoms("resid 5-100").n_residues96```Selections can be combined with boolean expressions. For example, to select the Cα atoms of all acidic residues [aspartic acid (“ASP”), glutamic acid (“GLU”), and histidines (named “HIS”, “HSD”, or “HSE”, depending on what force field is being used and what protonation state it is in)]:```>>> u.select_atoms("(resname ASP or resname GLU or resname HS*) and name CA")```We group with or separate selections by residue name (keyword resname). First either ASP, GLU, or any histidines are selected (we use “stemming” HS* to match any residue name that starts with “HS”). Then only those atoms whose name is “CA” are taken from the first set by an and selection. For convenience, the or in the first part of the selection can be taken implicitly with the shortcut syntax```>>> u.select_atoms("resname ASP GLU HS* and name CA")``` If you want to dig deeper into the selection syntaxis of MDanalysis you can read the full documentation [here](https://docs.mdanalysis.org/1.0.0/documentation_pages/selections.html) Experimental AspectsFor this tutorial we are going to use a MD trajectory from DNA binding domain of the HIV integrase 1, as you can see in the PDB [entry](https://www.rcsb.org/structure/1IHV) this structure was solved using NMR, and was found to form a dimer in solution.Here we are going to analyze 1000 frames of the integrase the monomeric and dimeric states and compare their RMSD, RMSF and measure distances. FIGURE 3. Cartoon representation of the structure of HIV integrase 1 (PDB 1IHV) Part 0 Downloading and Installing the required software  Installation we must install the softwares to perform this tutorial. Namely:- **MDAnalysis** for analyzing the data in molecular dynamics trajectories.- **py3Dmol** for visualization of the protein structure.<jupyter_code>!pip3 install --upgrade MDAnalysis
# Import MDanalysis
import MDAnalysis as mda
#from MDAnalysis.tests.datafiles import PSF, DCD, DCD2
from MDAnalysis.analysis import gnm
import matplotlib.pyplot as plt
%matplotlib inline
#Installing py3Dmol using pip
!pip install py3Dmol
#Importing py3Dmol for safety
import py3Dmol
!wget http://www.rcsb.org/pdb/files/1IHV.pdb.gz 
!gunzip 1IHV.pdb.gz
# We can visualize the dimeric conformation 
view=py3Dmol.view()
view.addModel(open('1IHV.pdb', 'r').read(),'pdb')
#Zooming into all visualized structures 
view.zoomTo()
#Here we set the background color as white
view.setBackgroundColor('white')
#Here we set the visualization style for chain B and C
view.setStyle({'cartoon': {'color':'purple'}})
#And we finally visualize the structures using the command below
view.show()<jupyter_output><empty_output><jupyter_text>Downloading MD trajectories<jupyter_code>#Here we copy to our Colab instance the trajectory files for the monomer and dimer
!wget https://github.com/pb3lab/ibm3202/raw/master/files/md_files/1ihv_dimer_protPBC.xtc
!wget https://github.com/pb3lab/ibm3202/raw/master/files/md_files/1ihv_mon_protPBC.xtc
!wget https://github.com/pb3lab/ibm3202/raw/master/files/md_files/1ihv_mon_protPBC.gro
!wget https://github.com/pb3lab/ibm3202/raw/master/files/md_files/1ihv_dimer_protPBC.gro<jupyter_output><empty_output><jupyter_text>Part I – Calculating $RMSD$ and $RMSF$  I.1 - RMSD !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Please note that the following cells are part of the tutorial of MDanalysis available [here](https://userguide.mdanalysis.org/stable/examples/analysis/alignment_and_rms/aligning_trajectory_to_frame.html) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!© Copyright 2019-2020, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, and <NAME>.  I.1A - Calculating $RMSD$ against a reference frame<jupyter_code>import MDAnalysis as mda
from MDAnalysis.analysis import align, rms<jupyter_output><empty_output><jupyter_text>First we need to load our trajectory files onto MDanalysis. This is done by creating an instance of an **Universe** object.<jupyter_code>#Here we create two Universes each containing the same trajectory of the monomeric trajectory, 
#one called mobile and the other ref which will be used as reference 
mobile = mda.Universe("/content/1ihv_mon_protPBC.gro", "/content/1ihv_mon_protPBC.xtc")
ref = mda.Universe("/content/1ihv_mon_protPBC.gro", "/content/1ihv_mon_protPBC.xtc")<jupyter_output><empty_output><jupyter_text>While `align.alignto` aligns single structures, or a frame of a trajectory, `align.AlignTraj` efficiently aligns an entire trajectory to a reference.We first check the $RMSD$ of our unaligned trajectory so we can compare results later. The code below sets the `mobile` trajectory to the last frame by indexing the last timestep, `ref` to the first frame by indexing the first timestep, and computes the root mean squared deviation between the $\alpha$-carbon positions.<jupyter_code>mobile.trajectory[-1]  # set mobile trajectory to last frame
ref.trajectory[0]  # set reference trajectory to first frame

mobile_ca = mobile.select_atoms('name CA')
ref_ca = ref.select_atoms('name CA')
rms.rmsd(mobile_ca.positions, ref_ca.positions, superposition=False)<jupyter_output><empty_output><jupyter_text>Now we can align the trajectory. We have already set ref to the first frame. In the cell below, we load the positions of the trajectory into memory so we can modify the trajectory in Python.<jupyter_code>aligner = align.AlignTraj(mobile, ref, select='name CA', in_memory=True).run()
mobile.trajectory[-1]  # set mobile trajectory to last frame
ref.trajectory[0]  # set reference trajectory to first frame

mobile_ca = mobile.select_atoms('name CA')
ref_ca = ref.select_atoms('name CA')
rms.rmsd(mobile_ca.positions, ref_ca.positions, superposition=False)<jupyter_output><empty_output><jupyter_text>**QUESTION❓:** How much does the RMSD before and after alignment compare?   I.1B - RMSD of a Universe with multiple selections over time It is more efficient to use the MDAnalysis.analysis.rms.RMSD class to calculate the $RMSD$ of an entire trajectory to a single reference point.The rms.RMSD class first performs a rotational and translational alignment of the target trajectory to the reference universe at ref_frame, using the atoms in select to determine the transformation. Then, without further alignment, the RMSD of each group in the `groupselections` argument is calculated.[Source](https://userguide.mdanalysis.org/stable/examples/analysis/alignment_and_rms/rmsd.html)<jupyter_code>#Here we create two Universes each containing the same trajectory of the monomeric trajectory, 
#one called mobile and the other ref which will be used as reference 
monomer_mobile = mda.Universe("/content/1ihv_mon_protPBC.gro", "/content/1ihv_mon_protPBC.xtc")
monomer_ref = mda.Universe("/content/1ihv_mon_protPBC.gro", "/content/1ihv_mon_protPBC.xtc")
rms.rmsd(monomer_mobile.select_atoms('backbone').positions,  # coordinates to align
         monomer_ref.select_atoms('backbone').positions,  # reference coordinates
         center=True,  # subtract the center of geometry
         superposition=True)  # superimpose coordinates
# Here we define two 
Loop1 = 'backbone and resid 227-240'
Loop2 = 'backbone and resid 252-257'
#Here we calculate the RMSD
R_rmsd = rms.RMSD(mobile,  # universe to align
             ref,  # reference universe or atomgroup
             select='backbone',  # group to superimpose and calculate RMSD
             groupselections=[Loop1, Loop2],  # groups for RMSD
             ref_frame=0)  # frame index of the reference
R_rmsd.run()<jupyter_output><empty_output><jupyter_text>The data is saved in R_rmsd.results.rmsd as an array. We can check the dimensions of the array using the *shape* attribute.<jupyter_code>R_rmsd.results.rmsd.shape<jupyter_output><empty_output><jupyter_text>The variable `R_rmsd.results.rmsd` has a row for each timestep. The first two columns of each row are the frame index of the time step, and the time (which is guessed in trajectory formats without timesteps). The third column is $RMSD$ of the `select` argument. The last few columns are the $RMSD$ of the groups in `groupselections`.  Plotting the data We can easily plot this data using the common data analysis package pandas. We turn the `R_rmsd.results.rmsd` array into a DataFrame and label each column below.<jupyter_code>import pandas as pd
df_rmsd_mono = pd.DataFrame(R_rmsd.results.rmsd,
                  columns=['Frame', 'Time (ns)','Backbone','Loop1','Loop2'])

df_rmsd_mono<jupyter_output><empty_output><jupyter_text>Here we use Plotly to easily create an interactive plot<jupyter_code>import plotly.graph_objects as go
import plotly.express as px
fig = px.line(df_rmsd_mono, x="Frame", y="Backbone",
        line_shape="spline", render_mode="svg", 
        labels={ "Backbone": "RMSD(Å)" })
fig.add_scatter(x=df_rmsd_mono["Frame"], y=df_rmsd_mono["Loop1"], name="Loop 1", showlegend=True )
fig.add_scatter(x=df_rmsd_mono["Frame"], y=df_rmsd_mono["Loop2"], name="Loop 2")
fig.add_scatter(x=df_rmsd_mono["Frame"], y=df_rmsd_mono["Backbone"], name="Backbone" )
fig.show()<jupyter_output><empty_output><jupyter_text>**QUESTION:** What is the range (in angstroms) of the RMSD fluctuations?  I.1B - Now is your turn to calculate the RMSD of the dimer<jupyter_code>#Here we create two Universes each containing the same trajectory of the monomeric trajectory, 
#one called mobile and the other ref which will be used as reference 
# dimer_mobile = mda.Universe("/content/1ihv_dimer_protPBC.gro", "/content/1ihv_dimer_protPBC.xtc")
# dimer_ref = mda.Universe("/content/1ihv_dimer_protPBC.gro", "/content/1ihv_dimer_protPBC.xtc")
# rms.rmsd(dimer_mobile.select_atoms('backbone').positions,  # coordinates to align
#          dimer_ref.select_atoms('backbone').positions,  # reference coordinates
#          center=True,  # subtract the center of geometry
#          superposition=True)  # superimpose coordinates
# Loop1A = 'backbone and resid 227-240'
# Loop2A = 'backbone and resid 252-257'
# R_rmsd_dimer = rms.RMSD(dimer_mobile,  # universe to align
#              dimer_ref,  # reference universe or atomgroup
#              select='backbone',  # group to superimpose and calculate RMSD
#              groupselections=[Loop1A, Loop2A],  # groups for RMSD
#              ref_frame=0)  # frame index of the reference
# R_rmsd_dimer.run()<jupyter_output><empty_output><jupyter_text>The data is saved in `R_rmsd.results.rmsd` as an array. We can check the dimensions of the array using the `shape` attribute.<jupyter_code>R_rmsd_dimer.results.rmsd.shape
import pandas as pd
#Here we create the pandas dataframe from the R_rmsd_dimer.rmsd object
df_rmsd_dimer = pd.DataFrame(R_rmsd_dimer.results.rmsd,
                  columns=['Frame', 'Time (ns)','Backbone','Loop1A','Loop2A'])

df_rmsd_dimer<jupyter_output><empty_output><jupyter_text>Lets plot the $RMSD$ over time for the Monomer and the Dimer<jupyter_code>import plotly.express as px

fig = px.line(df_rmsd_mono, x="Frame", y="Backbone",
        line_shape="spline", render_mode="svg", 
        labels={ "Backbone": "RMSD(Å)" })
fig.add_scatter(x=df_rmsd_mono["Frame"], y=df_rmsd_mono["Backbone"], name="Backbone Monomer" )
fig.add_scatter(x=df_rmsd_dimer["Frame"], y=df_rmsd_dimer["Backbone"], name="Backbone Dimer AVG" )
fig.show()<jupyter_output><empty_output><jupyter_text>Now is your turn to explore if there are any changes in the RMSD of the loops over time<jupyter_code>#Hints:
fig.add_scatter(x=df_rmsd_mono["Frame"], y=df_rmsd_mono["Loop1"], name="Loop 1", showlegend=True )
fig.add_scatter(x=df_rmsd_mono["Frame"], y=df_rmsd_mono["Loop2"], name="Loop 2")<jupyter_output><empty_output><jupyter_text>I.2 - RMSF  Now, we want to assess the average atomic fluctuations during the MD trajectories for both the monomeric and dimeric states of the integrase<jupyter_code>#First we need to make sure that our universes are properly aligned

aligner = align.AlignTraj(monomer_mobile, monomer_ref, select='name CA', in_memory=True).run()
aligner = align.AlignTraj(dimer_mobile, dimer_ref, select='name CA', in_memory=True).run()
#Here we create a selection of the previously aligned trajectory 
c_alphas_monomer = monomer_mobile.select_atoms('protein and name CA')
c_alphas_dimer = dimer_mobile.select_atoms('protein and name CA')
R_rmsf_mono = rms.RMSF(c_alphas_monomer).run()
R_rmsf_dimer = rms.RMSF(c_alphas_dimer).run()
rms.RMSF 
c_alphas_dimer.resids
import pandas as pd
#Here we create a pandas dataframe
df_rmsf_mono = pd.DataFrame(R_rmsf_mono.results.rmsf,
                  columns=['BackboneRMSF'])
df_rmsf_mono = df_rmsf_mono.assign(Residue = c_alphas_monomer.resids)

df_rmsf_dimer = pd.DataFrame(R_rmsf_dimer.results.rmsf,
                  columns=['BackboneRMSF'])
df_rmsf_dimer = df_rmsf_dimer.assign(Residue = c_alphas_dimer.resids)


df_rmsf_dimer_A = df_rmsf_dimer.head(52)
df_rmsf_dimer_B = df_rmsf_dimer.tail(52)


import plotly.express as px

fig = px.line(df_rmsf_mono, x="Residue", y="BackboneRMSF",
        line_shape="linear", render_mode="svg", 
        labels={ "BackboneRMSF": "RMSF(Å)" , "Residue":"Residue Number"}, color=None)
fig.add_scatter(x=df_rmsf_dimer_A["Residue"], y=df_rmsf_dimer_A["BackboneRMSF"], name="Dimer Chain A", line_shape="linear")
fig.add_scatter(x=df_rmsf_dimer_B["Residue"], y=df_rmsf_dimer_B["BackboneRMSF"], name="Dimer Chain B", line_shape="linear")
fig.add_scatter(x=df_rmsf_mono["Residue"], y=df_rmsf_mono["BackboneRMSF"], name="Monomer", line_shape="linear")

fig.show()<jupyter_output><empty_output><jupyter_text>**QUESTIONS❓**1. Is there any difference in local fluctuations between the monomeric and dimeric states?2. Which region exhibit a high atomic fluctuation through the trajectory?3. What are the structural features of these regions?<jupyter_code>#Visualize using Py3Dmol the monomeric state<jupyter_output><empty_output><jupyter_text>I.3 - Pairwise RMSD<jupyter_code>#Import modules
import MDAnalysis as mda
from MDAnalysis.analysis import diffusionmap, align<jupyter_output><empty_output><jupyter_text>Pairwise RMSDs are an effective way to quickly view similarities and differences in conformations (as measured by RMSD) across a entire trajectory and not only in comparison to just one reference frame.We are going to use the previously aligned trajectories **monomer_mobile** and **dimer_mobile** We can then calculate a pairwise $RMSD$ matrix with the `diffusionmap.DistanceMatrix` class, by using the default the rms.rmsd metric.<jupyter_code>#Monomer distance matrix calculation
matrix1 = diffusionmap.DistanceMatrix(monomer_mobile, select='name CA').run()
#Dimer distance matrix calculation
matrix2 = diffusionmap.DistanceMatrix(dimer_mobile, select='name CA').run()<jupyter_output><empty_output><jupyter_text>The results array is in `matrix.results.dist_matrix` as a square array with the shape (n_frames, n_frame).<jupyter_code>print(matrix1.results.dist_matrix.shape)
print(matrix2.results.dist_matrix.shape)<jupyter_output><empty_output><jupyter_text>We can use the common plotting package matplotlib to create a heatmap from this array.<jupyter_code>#Here we plot the Monomer RMSD matrix
plt.imshow(matrix1.results.dist_matrix, cmap='viridis')
plt.xlabel('Frame')
plt.ylabel('Frame')
plt.colorbar(label='RMSD (Angstrom)')

#Here we plot the Dimer RMSD matrix
plt.imshow(matrix2.results.dist_matrix, cmap='viridis')
plt.xlabel('Frame')
plt.ylabel('Frame')
plt.colorbar(label='RMSD (Angstrom)')<jupyter_output><empty_output><jupyter_text>Appendix A - Normal Mode and Principal Component Analysis    I - Normal Mode analysis long range contacts<jupyter_code>monomer_mobile = mda.Universe("/content/1ihv_mon_protPBC.gro", "/content/1ihv_mon_protPBC.xtc")
dimer_mobile = mda.Universe("/content/1ihv_dimer_protPBC.gro", "/content/1ihv_dimer_protPBC.xtc")
nma1 = gnm.GNMAnalysis(monomer_mobile, select='protein and name CA', cutoff=7.0)
nma1.run()
nma2 = gnm.GNMAnalysis(dimer_mobile,
                      select='protein and name CA',
                      cutoff=7.0)
nma2.run()

len(nma2.results)
%matplotlib inline
#sns.set_context('notebook')
%config InlineBackend.figure_format = 'retina'

## we plot the distribution of eigenvalues. The dominant conformation state is represented by the peak at:
eigenvalues1 = [res[1] for res in nma1.results]
eigenvalues2 = [res[1] for res in nma2.results]

histfig, histax = plt.subplots(nrows=2, sharex=True, sharey=True)
histax[0].hist(eigenvalues1)
histax[1].hist(eigenvalues2)

histax[1].set_xlabel('Eigenvalue')
histax[0].set_ylabel('Frequency (Monomer)')
histax[1].set_ylabel('Frequency (Dimer)');
plt.show()

import pandas as pd
import numpy as np

##Create Panda Dataframe Files
eu = pd.DataFrame({'LC': eigenvalues1,'NLC': eigenvalues2}) 

#Save Panda DataFrame
eu.to_csv('./DF_1.csv')


#inspect Dataframe
eu.head()
time1 = [res[0] for res in nma1.results]
time2 = [res[0] for res in nma2.results]
linefig, lineax = plt.subplots()
plt.plot(time1, eigenvalues1, label='Monomer')
plt.plot(time2, eigenvalues2, label='Dimer')
lineax.set_xlabel('Time (ps)')
lineax.set_ylabel('Eigenvalue')
plt.legend();
plt.show()<jupyter_output><empty_output><jupyter_text>II -PCA<jupyter_code>#Import things
import MDAnalysis as mda
from MDAnalysis.analysis import diffusionmap, align
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>**WARNING!!!**For best results, your trajectory should be aligned on your atom group selection before you run the analysis. Setting align=True will not give correct results in the PCA.<jupyter_code># Align the trajectory 
aligner1 = align.AlignTraj(monomer_mobile, monomer_mobile, select='backbone', in_memory=True).run()
aligner2 = align.AlignTraj(dimer_mobile, dimer_mobile, select='backbone', in_memory=True).run()<jupyter_output><empty_output><jupyter_text>Overview of the method**Principal component analysis (PCA)** is a statistical technique that decomposes a system of observations into linearly uncorrelated variables called principal components. These components are ordered so that the first principal component accounts for the largest variance in the data, and each following component accounts for lower and lower variance. PCA is often applied to molecular dynamics trajectories to **extract the large-scale conformational motions or “essential dynamics” of a protein**. The frame-by-frame conformational fluctuation can be considered a linear combination of the essential dynamics yielded by the PCA.In MDAnalysis, the method is as follows:> Optionally align each frame in your trajectory to the first frame.Construct a 3N x 3N covariance for the N atoms in your trajectory. Optionally, you can provide a mean; otherwise the covariance is to the averaged structure over the trajectory.Diagonalise the covariance matrix. The eigenvectors are the principal components, and their eigenvalues are the associated variance.Sort the eigenvalues so that the principal components are ordered by variance.<jupyter_code>import MDAnalysis as mda
import MDAnalysis.analysis.pca as pca
from MDAnalysis.coordinates.base import Timestep

import numpy as np
import os
import glob

import pandas as pd
import seaborn as sns

import matplotlib.pyplot as plt
import matplotlib.cm
import matplotlib.ticker as ticker
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Call the PCA functionYou can choose how many principal components to save from the analysis with n_components. The default value is None, which saves all of them. You can also pass a mean reference structure to be used in calculating the covariance matrix. With the default value of None, the covariance uses the mean coordinates of the trajectory.<jupyter_code>pcu1 = pca.PCA(monomer_mobile, select='protein and backbone',
             align=False, mean=None,
             n_components=None).run()
pcu2 = pca.PCA(dimer_mobile, select='protein and backbone',
             align=False, mean=None,
             n_components=None).run()<jupyter_output><empty_output><jupyter_text>The principal components are saved in pc.p_components. If you kept all the components, you should have an array of shape (natoms×3,natoms×3)<jupyter_code>backbone1 = monomer_mobile.select_atoms('protein and backbone')
n_bb1 = len(backbone1)
print('There are {} backbone atoms in the analysis'.format(n_bb1))
print(pcu1.p_components.shape)

backbone2 = dimer_mobile.select_atoms('protein and backbone')
n_bb2 = len(backbone2)
print('There are {} backbone atoms in the analysis'.format(n_bb2))
print(pcu2.p_components.shape)<jupyter_output><empty_output><jupyter_text>Get the variance of the first component<jupyter_code>pcu1.variance[0],pcu2.variance[0]<jupyter_output><empty_output><jupyter_text>This variance is somewhat meaningless by itself. It is much more intuitive to consider the variance of a principal component as a percentage of the total variance in the data. MDAnalysis also tracks the percentage cumulative variance in pc.cumulated_variance. As shown below, the first principal component contains 90.3% the total trajectory variance. The first three components combined account for 96.4% of the total variance.<jupyter_code>print(pcu1.cumulated_variance[0])
print(pcu1.cumulated_variance[2])
print(pcu2.cumulated_variance[0])
print(pcu2.cumulated_variance[2])
plt.plot(pcu1.cumulated_variance[:10])
plt.xlabel('Principal component')
plt.ylabel('Cumulative variance');

plt.plot(pcu2.cumulated_variance[:10])
plt.xlabel('Principal component')
plt.ylabel('Cumulative variance');<jupyter_output><empty_output><jupyter_text>Visualising projections into a reduced dimensional space The pc.transform() method transforms a given atom group into weights $w_i$ over each principal component $i$.$w_i(t)=(r(t)−r⎯⎯⎯)⋅u_i$$r(t)$ are the atom group coordinates at time t, r⎯⎯⎯ are the mean coordinates used in the PCA, and ui is the ith principal component eigenvector u.While the given atom group must have the same number of atoms that the principal components were calculated over, it does not have to be the same group.Again, passing n_components=None will tranform your atom group over every component. Below, we limit the output to projections over 5 principal components only.<jupyter_code>transformed1 = pcu1.transform(backbone1, n_components=5)


transformed2 = pcu2.transform(backbone2, n_components=5)
transformed1.shape, transformed2.shape<jupyter_output><empty_output><jupyter_text>The output has the shape (n_frames, n_components). For easier analysis and plotting we can turn the array into a DataFrame.<jupyter_code>df1 = pd.DataFrame(transformed1,
                  columns=['PC{}'.format(i+1) for i in range(5)])
df1['Time (ns)'] = df1.index * monomer_mobile.trajectory.dt
df1.head()

df2 = pd.DataFrame(transformed2,
                  columns=['PC{}'.format(i+1) for i in range(5)])
df2['Time (ns)'] = df2.index * dimer_mobile.trajectory.dt
df2.head()<jupyter_output><empty_output><jupyter_text>There are several ways we can visualise the data. Using the Seaborn’s PairGrid tool is the quickest and easiest way, if you have seaborn already installed.<jupyter_code>import seaborn as sns

g1 = sns.PairGrid(df1, hue='Time (ns)',
                 palette=sns.color_palette('Oranges_d',
                                           n_colors=len(df1)))
g1.map(plt.scatter, marker='.')
g2 = sns.PairGrid(df2, hue='Time (ns)',
                 palette=sns.color_palette('Oranges_d',
                                           n_colors=len(df2)))
g2.map(plt.scatter, marker='.')<jupyter_output><empty_output><jupyter_text>Another way to investigate the essential motions of the trajectory is to project the original trajectory onto each of the principal components, to visualise the motion of the principal component. The product of the weights $w_i(t)$ for principal component $_i$ with the eigenvector $u_i$ describes fluctuations around the mean on that axis, so the projected trajectory $r_i(t)$ is simply the fluctuations added onto the mean positions r⎯⎯⎯.$r_i(t)=w_i(t)×u_i+r$Below, we generate the projected coordinates of the first principal component. The mean positions are stored at pc.mean.<jupyter_code>pc1u1 = pcu1.p_components[:, 0]
trans1u1 = transformed1[:, 0]
projectedu1 = np.outer(trans1u1, pc1u1) + pcu1.mean
coordinatesu1 = projectedu1.reshape(len(trans1u1), -1, 3)

pc1u2 = pcu2.p_components[:, 0]
trans1u2 = transformed2[:, 0]
projectedu2 = np.outer(trans1u2, pc1u2) + pcu2.mean
coordinatesu2 = projectedu2.reshape(len(trans1u2), -1, 3)<jupyter_output><empty_output><jupyter_text>We can create a new universe from this to visualise the movement over the first principal component.<jupyter_code>!pip3 install nglview
import nglview as nv
proj1 = mda.Merge(backbone1)
proj1.load_new(coordinatesu1, order="fac")
view = nv.show_mdanalysis(proj1.atoms)
view<jupyter_output><empty_output><jupyter_text>If you have nglview installed, you can view the trajectory in the notebook. Otherwise, you can write the trajectory out to a file and use another program such as VMD. Below, we create a movie of the component.<jupyter_code>!pip install moviepy==0.2.2.11
!pip install imageio==1.6
from nglview.contrib.movie import MovieMaker
movie = MovieMaker(view, output='pc1u1.gif', in_memory=True)
movie.make()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.101871557336004, 'fraction_numerical': 0.01114171256885328, 'mean_word_length': 5.599132589838909, 'pattern_counts': {'":': 4, '<': 145, '<?xml version=': 0, '>': 179, 'https://': 10, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '34', 'dir': 'jupyter-structured-clean-dedup', 'id': '26830', 'license': 'MIT', 'n_tokens_mistral': 10319, 'n_tokens_neox': 9717, 'n_words': 3368, 'path': 'tutorials/lab08_MDanalysis.ipynb', 'repo_name': 'paulshamrat/biopy', 'size': '32013'}