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', '') |
