Analyzing molecular trajectories with MDAnalysis in SAMSON

Facebooktwittergoogle_plusredditpinterestlinkedinmailFacebooktwittergoogle_plusredditpinterestlinkedinmail

MDAnalysis is a well-known Python library for analyzing molecular dynamics trajectories. In this tutorial, we demonstrate how it can be used to analyze trajectories saved in Path nodes in SAMSON and superimpose structures.

Analyzing trajectories

Let us, for example, compute the radius of gyration and the RMSD for a protein with trajectories saved in Path in SAMSON. To use the functionality of MDAnalysis, we need to export trajectories we would like to analyze from SAMSON into a format supported by MDAnalysis (e.g., XYZ format). These trajectory files are then provided to MDAnalysis alongside with a topology file.

Let us first import MDAnalysis modules that we are going to use in this tutorial:

import MDAnalysis
import MDAnalysis.analysis.rms

We will export trajectories from SAMSON (Path in terms of SAMSON) each in a separate xyz-file. For that, we define the following function, which has as parameters a Path and a prefix filename for xyz-files:

def exportTrajectories(sbpath, filename):
    '''
    Exports trajectories from a Path node 'sbpath' to xyz files with a name starting with 'filename'
    '''
    trajectory_files = []                                               # a list of trajectory files
    sbpath.currentStep = 0                                              # set currentStep for the Path to 0
 
    for step in range(sbpath.numberOfSteps):                            # loop over steps in the Path
        sbpath.currentStep = step                                       # increment currentStep
 
        fn = filename + str(step) + '.xyz'                              # a name of a file 
        trajectory_files.append(fn)                                     # append list of trajectory files
 
        indexer = SAMSON.getNodes('n.t a')                              # get a node indexer for all atoms
 
        SAMSON.exportToFile(indexer, fn, '')                            # export current trajectory into a file 'fn'
 
    return trajectory_files                                             # return list of trajectory files

To compute the radius of gyration and RMSD we define the following functions, which use the functionality of MDAnalysis and have as an input parameter MDAnalysis.Universe. Note, that the basic unit of length in MDAnalysis is the Angstrom.

def computeRgyr(u):
    '''
    Compute the radius of gyration thanks to MDAnalysis
    '''
    bb = u.select_atoms('protein and backbone')                         # a selection (AtomGroup)
 
    for ts in u.trajectory:                                             # loop over trajectories
        rgyr = bb.radius_of_gyration()                                  # get the radius of gyration
        print("Rgyr [frame: {0}] = {1} A".format(ts.frame, rgyr))
 
def computeRMSD(u):
    '''
    Compute the RMSD between backbone atoms
    '''
    bb = u.select_atoms('backbone')                                     # a selection (AtomGroup)
    u.trajectory[0]                                                     # choose the first frame
    frame_0 = bb.positions                                              # coordinates of first frame
 
    for i in range(1, u.trajectory.n_frames):                           # loop over trajectories
        u.trajectory[i]                                                 # forward to the next frame
        frame_i = bb.positions                                          # coordinates of the current frame
        rmsd = MDAnalysis.analysis.rms.rmsd(frame_0, frame_i)           # get the radius of gyration
        print("RMSD [frame: {0}] = {1} A".format(i, rmsd))

Now, we can use these function to compute the desired parameters for a chosen Path. We consider that a Path for which we want to compute these parameters is the first one.

sbpaths = SAMSON.getNodes('n.t path')                                   # get all Path nodes
 
if sbpaths.size != 0:
    sbpath = sbpaths[0]                                                 # get the first Path
 
    trajectory_files = exportTrajectories(sbpath, '/path/to/trajectories/traj-')
 
    topology_file = '/path/to/topology/protein.pdb'                     # a topology file for this molecule
 
    u = MDAnalysis.Universe(topology_file, trajectory_files)            # read from a list of trajectories
 
    computeRgyr(u)
    computeRMSD(u)

Superposition of structure

Let us consider a case, when we have two structures (e.g., two structural models in terms of SAMSON) in the active Document in SAMSON and we want to superimpose them in a way that minimizes the RMSD. For that, we can use functions from the MDAnalysis.analysis.align module (see MDAnalysis: Superposition of structure). In this example, we consider the first molecule to be the reference one, and the second one to be the one which should be superimposed on the first one.

Let us import MDAnalysis modules that we are going to use in this tutorial:

import MDAnalysis
from MDAnalysis.analysis import align, rms

Next, we need to export structures from SAMSON into files in a format supported by MDAnalysis (e.g., pdb-format). Note, that in this example each structure is in a separate structural model in SAMSON. In this case, we need to get an indexer of all structural models and export an indexer of each structural model into a separate file.

indexer = SAMSON.getNodes('n.t sm')                                     # get an indexer of structural models            
 
ref = indexer[0]                                                        # get the first structural model    
mob = indexer[1]                                                        # get the second structural model    
ref_indexer = ref.getNodes()                                            # get an indexer of nodes for the first structural model
mob_indexer = mob.getNodes()
 
filename_ref = '/path/to/files/ref.pdb'
filename_mob = '/path/to/files/mob.pdb'
 
SAMSON.exportToFile(ref_indexer, filename_ref, '')                      # export the reference structure
                                                                        # the third input parameter is for options used for importing: '' is for default
SAMSON.exportToFile(mob_indexer, filename_mob, '')                      # export the structure which should be rotated

Now, we can open and process these structures in MDAnalysis:

u_ref = MDAnalysis.Universe(filename_ref)
u_mob = MDAnalysis.Universe(filename_mob)

For example, we can check the C-alpha RMSD between these structures:

rms.rmsd(u_mob.atoms.CA.positions, u_ref.atoms.CA.positions)

We need to take care of translations. For that, we shift molecules by their centers of mass:

u_ref0 = u_ref.atoms.CA.positions - u_ref.atoms.CA.center_of_mass()     # get a structure shifted by its center of mass
u_mob0 = u_mob.atoms.CA.positions - u_mob.atoms.CA.center_of_mass()     # get a structure shifted by its center of mass

Next, we can compute the rotation matrix that superposes mob on ref while minimizing the CA-RMSD:

R, rmsd = align.rotation_matrix(u_mob0, u_ref0)
print(rmsd)
print(R)

Now, we can superimpose mob on ref inside SAMSON using the rotation matrix and centers of masses for these structures. For that, we need to convert the rotation matrix and centers of masses into SAMSON Types. Note, that the basic unit of length in MDAnalysis is the Angstrom.

u_ref_com = u_ref.atoms.CA.center_of_mass()
u_mob_com = u_mob.atoms.CA.center_of_mass()
 
ref_center_of_mass = Type.vector3(Quantity.angstrom(u_ref_com[0]), Quantity.angstrom(u_ref_com[1]), Quantity.angstrom(u_ref_com[2]))
mob_center_of_mass = Type.vector3(Quantity.angstrom(u_mob_com[0]), Quantity.angstrom(u_mob_com[1]), Quantity.angstrom(u_mob_com[2]))
 
rotation_matrix = Type.matrix33(R.tolist())                             # create a rotation matrix in SAMSON
                                                                        # R is an ndarray and should be converted into a list

To translate and rotate the mob structure, we need to apply translations and the rotation to each atom of the structure. For that, we need to get an indexer of all atoms in the mob structure.

mob_indexer_a = mob.getNodes('n.t a')                                   # get an indexer of all atoms in the mob
 
for a in mob_indexer_a:                                                 # loop over all atoms
    a.setPosition( a.getPosition() - mob_center_of_mass)                # translate by mob's center of mass
    a.setPosition( rotation_matrix * a.getPosition() )                  # rotate
    a.setPosition( a.getPosition() + ref_center_of_mass)                # translate by ref's center of mass

Alternatively, one can superimpose mob on ref using MDAnalysis, write the resulting structure into a pdb-file and import it into SAMSON:

mob.atoms.translate(-mob.atoms.CA.center_of_mass())
mob.atoms.rotate(R)
mob.atoms.translate(ref.atoms.CA.center_of_mass())
mob.atoms.write('/path/to/files/mob_on_ref.pdb')
 
SAMSON.importFromFile('/path/to/files/mob_on_ref.pdb', '')
Facebooktwittergoogle_plusredditpinterestlinkedinmailFacebooktwittergoogle_plusredditpinterestlinkedinmail

Leave a Reply

Your email address will not be published.

Time limit is exhausted. Please reload CAPTCHA.

  • Subscribe

    Subscribe to make sure you receive all our updates, tips and tutorials