Amber Workshop - tutorial 8

(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2005

Case Study: All Atom Structure Prediction and Folding Simulations
of a Stable Protein (Folding Trp-Cage Peptide)

We will start our analysis of the results by plotting temperature, total energy, kinetic energy and potential energy from the output file. Examining these properties will allow us to check that nothing strange happened during the MD. We can check that the temperature rose smoothly and equilibrated at 325 K. Similarly changes in the kinetic energy and potential energy should be smooth and should level off at equilibrated values. Any sudden spikes in the temperature or energy plots are indicative of a problem with our simulation protocol, bad starting structure, too long a time step, inappropriate parameters etc.

In order to extract the data from the output files we can use the following perl script that automates the process (process_mdout.perl). We can then process the output files in one go by listing them all on the command line as follows:

mkdir analysis
cd analysis
../process_mdout.perl ../heat1.out ../heat2.out ../heat3.out ../heat4.out ../heat5.out ../heat6.out ../heat7.out ../equil1.out ../equil2.out ../equil3.out ../equil4.out ../equil5.out ../equil6.out ../equil7.out ../equil8.out ../equil9.out ../equil10.out

Here is a tar file containing all of the files produced by process_mdout.perl - analysis.tar.gz (2.3 Mb)

We can now use a graphing program, such as xmgr which we used in tutorial 1, to plot the temperature and energy data.

This looks fine, our temperature has equilibrated around 325 K as we hoped. The heating phase also shows a nice gradual increase in temperature, with no sudden jumps so it looks like things are okay here.

Next we look at the energy plots. We are most interested in the potential energy plot since we need to locate the structure with the lowest energy so that we can reproduce the RMSd plots they show in the paper. The structure with the lowest energy we will be labelling our folded state.

This looks good, the kinetic energy is stable as we'd expect since the temperature plot, which is directly proportional to the kinetic energy, was stable. The total and potential energy plots show an initial decrease and then stabilise implying that our system has folded up to a state more stable than the starting linear structure. By plotting just the potential energy we can see the decrease more clearly. We can also plot a running average over 10ps as they did in the paper which makes things even clearer. Note, in xmgr you can add a running average by going to Data->Transformations->Running averages.

This is what they show as Figure 1(A) in the paper. Note our initial energy is higher than that shown in the paper. This is most likely due to differences in the generalised Born term since they used AMBER 6 while we used AMBER 8. AMBER 7 and 8 have slightly different Born radii to AMBER 6. It is also possible that we have not folded to the same structure as they do in the paper. Further analysis will be able to tell us more.

The energy initially increases very sharply, this is due to us heating the system but then it starts to decrease again as our system folds up. Overall we see a decrease of around 40 kcal/mol over the 50 ns. This is in line with what they state in paragraph 4 of the paper.

The next step is to locate the structure with the lowest energy, this we will call our folded state and we will use it to reproduce Figure 1 B from the paper. We can locate the lowest energy by looking at the summary.EPTOT file that process_mdout.perl produced. But first we need to strip off the first 50ps of the data as this represents the heating phase. Here is the file with this section removed. (summary_EPTOT_stripped.dat.gz). The following awk script will quickly find the lowest value in the file for us. It will print a value every time it finds a new minimum value, the last value printed will be the lowest energy found.

>gunzip summary_EPTOT_stripped.dat.gz
>cat summary_EPTOT_stripped.dat | awk '{if($2<min) {min=$2;print $1" "min}}'

NOTE: your answers may be different here if you did your own simulation as variations in machine precision mean you will have explored a slightly different region of phase space.

So, the lowest energy was -547.4913 kcal/mol and it occurred 42.505 ns into our simulation. We can find out where this is in our output files by using grep:

So it occurred on step 1227500 of our 9th equilibration step. We wrote to the mdcrd file every 500 steps (ntwx=500) during the equilibration so this should represent frame 2455 (1227500/500) of equil9.mdcrd. Lets extract that structure using ptraj and take a look at it. Here is the ptraj script (extract_frame9_2455.trajin):

> $AMBERHOME/exe/ptraj TC5b.prmtop < extract_frame9_2455.trajin > extract_frame9_2455.out

This looks pretty good, the system is well folded with the tryptophan in a cage structure vaguely similar to that shown in Figure 2 of the paper. At this point we could compare this structure to the crystal structure but we have to remember that when the paper was published the experimental structure had not been determined so we shall leave this step until last. First off we will try to generate figure 1B which shows the backbone RMSd between the trajectory structures and this lowest energy structure. Here is ptraj script for this:

Check the output from ptraj for errors. The only errors you should see are ones telling you that frame 5001 of each of the trajectories is corrupt. This is not really an error since there were only 5000 frames in each mdcrd file anyway so we can safely ignore it.

Our rmsd plot shows significantly more fluctuation than that shown in the paper although we do see the formation of a plateau around 16 to 38 ns. However, this does not show the stability over 50 ns that they suggest in paragraph 5 of the paper. This might suggest many things. Chances are our trajectory has not folded up to the native structure but has instead folded to an alternative structure. This structure is fairly stable but obviously must unfold again after 38 ns. Indeed, if you watch equil8.mdcrd, equil9.mdcrd and equil10.mdcrd in vmd you will see the non-alpha helix part of the structure keep unfolding and refolding as the system tries to locate a more stable structure. The difference we see from the paper may be due to the implicit solvent model we chose or due to some other minor difference in the parameters. What this data tells us is that most likely we haven't run our simulation for long enough since the end point of the trajectory shows that the system wants to unfold and fold up differently. The starting structure we have used may be slightly different than that used in the paper and so while their's folded quickly ours gets itself trapped in a high energy minimum.

It is also possible that the slightly different Born radii we have used (they used the AMBER 6 Born radii remember) may have delayed the folding process. There are a number of ways we could extend this simulation further and try and fold up correctly. We could simply run a longer trajectory and wait for the system to leave the high energy basin it is trapped in. This may take a long time. Alternatively we could try a simulated annealing approach where we heat the system up and then let it cool again. We could also try a replica exchange simulation where we have many replicas and exchange properties between them. This avoids the problems of getting trapped in high energy minima. For the record I have obtained from the authors the structure that was published in the paper and can confirm that it gives a lower energy than our lowest energy structure. As such we have just been unlucky that our simulation wasn't long enough for the system to fold completely. Feel free to try extending it yourself if you wish or try some of the alternative approaches I mention above. There is a lesson to be learnt here. You can never really run a simulation for "too" long. Just because the system hasn't folded in 50 ns doesn't mean it won't fold in the next 50 ns so be aware of this when you are considering what your results imply.

Even though our system does not seem to have reached a completely stable structure it has at least folded up to a structure that shows stability over more than 20ns. Since we don't have the time to keep trying different simulation protocols or extend our simulation we shall go ahead with what we have and try several different types of analysis that we can do on our 50 ns of trajectory data.

Lets start by comparing our lowest energy structure to the NMR models that are now available (1L2Y.pdb). Like all NMR structures this pdb consists of an ensemble of possible structures, in this case 38. Ideally we should average these structures and compare the average structure to our theoretical structure. However, in the interests of speed I will skip this step and we can just compare the first structure in the file to our lowest energy structure. We can align the two structures and measure the RMSD using VMD. This is acceptable because, as the image below shows, residues 3 to 19 of the NMR structures are well conserved and it is just the ends of the chain that differ largely between the 38 structures.

Note: I am using VMD v1.8.3 here - If you are using a different version your menus may be slightly different. Note you need v1.8.2 or later for the RMSD fitting to work.

At the moment the structures won't line up correctly. This is because the coordinate frames are different.

In order to compare the structures we need to RMSD fit and align them. Fortunately VMD contains a tool that makes this very easy. You can access it by clicking Extensions->Analysis->RMSD Calculator:

We can now do the alignment. We will align the backbone's only so make sure the Backbone only box is highlighted. Also select residues 3 to 19 by typing residue 2 to 18 in the text box (REMEMBER VMD NUMBERS FROM ZERO). Then click Align to align the structures and then click RMSD to measure the rmsd between the two structures.

So we have a backbone RMSD of 3.72 angstroms. This is a far cry from the 1.4 angstroms they quote in the paper which means we have almost certainly not folded to the same structure. Here is the alignment as shown in VMD (after I have changed the graphical representation options to make it easier to compare. The NMR structure is shown in blue):

Indeed, the comparison looks very bad. The Tryptophan looks to be in a very different position at almost 90 degrees to the experimental position. However, this isn't the whole story. Here we fitted to residues 3 to 19 (2 to 18 in VMD speak.) What happens if we repeat the process and fit just the residues involved in the alpha helix, 3 to 9?

Open the RMSD fitting tool again and this time enter "residue 2 to 8" in the text box, then hit align and then RMSD:

We get a backbone RMSD for the alpha helix of just 0.7 angstroms. Indeed it fits very well. Things start to go wrong around residue 9 where the backbone twists the wrong way. After this the structures are very different. There is one thing you should note, however, take a close look at the ends of the chain in each of the two structures:

Clearly there is a difference in the hydrogen bonding at the ends of the protein chain. Our structure has at least two, maybe 3 hydrogen bonds that are acting to zip the ends of the structure closed. In the NMR structure on the other hand the hydrogen bonding is very different. So this may be one explanation as to why the structures are so different, perhaps these hydrogen bonds are what are keeping us trapped in the high energy minimum. Lets take a look then at how we can use AMBER's ptraj to look at hydrogen bonds over the course of our trajectory. This may also shed some light on why we get an initial structure formed early on in our trajectory that stays stable for almost 15 ns and yet has no real secondary structure (no alpha helix). If you watch the trajectory you should be able to see this structure and also observe how around 15 to 16 ns it suddenly snaps open and allows the system to refold with the alpha helix in place. We will also look at how we can use the MMTSB toolset to do some cluster analysis on our trajectory. This will allow us to assign the structures over time to various families so that we can comment on their lifetimes.

6.1.1) Converting MDCRD's to Binpos
The first thing we will do is convert our compressed text format mdcrd files into binpos format. This is a binary version of a trajectory file. We could do the subsequent steps easily enough with mdcrd files but using binpos is significantly faster. At the same time we will concatenate our 10 mdcrd files into a single binpos file. Here is the ptraj script to do it:

mdcrd_to_binpos.ptraj

trajin equil1.mdcrd.gz
trajin equil2.mdcrd.gz
trajin equil3.mdcrd.gz
trajin equil4.mdcrd.gz
trajin equil5.mdcrd.gz
trajin equil6.mdcrd.gz
trajin equil7.mdcrd.gz
trajin equil8.mdcrd.gz
trajin equil9.mdcrd.gz
trajin equil10.mdcrd.gz

trajout equil_1-10.binpos binpos

If all worked you should now have a 174 mb equil_1-10.binpos file in your working directory.

6.1.2) Calculating an Average Structure
Next we shall calculate an average structure over our entire trajectory (excluding the heating part).

Remember an average is just that, an average of the coordinates visited during the trajectory. The system is rotating and translating so you cannot hope to learn anything from this average. It is not a physically meaningful structure. Remember this when you are averaging coordinates.

However, if we are cleverer about calculating our average we can obtain more meaningful results. First of all we can do a mass weighted backbone RMS fit of every frame to the first frame in the trajectory, this will remove the rotation and translation aspect. Secondly we can look at a specific part of a trajectory where the average structure may have more meaning. Lets try this by looking at the section between 28 and 35 ns. This is frames 28000 to 35000 of the trajectory file.

This is much better and much more revealing. Try opening it yourself in VMD and exploring the structure. Load the TC5b.prmtop file first before loading the pdb file in this way the bonding will be sensible. Notice how some of the residues are very small, notably some of the hydrogen bonds lengths are tiny. This suggests that these residues are very dynamic and rotate a lot during this section of the trajectory. If you rotate something around a bond your average coordinate will be a point on that centre of rotation. However, note that the Tryptophan looks almost perfect. This must hardly move over this section of the trajectory. Similarly the backbone is well formed here which indicates that the folded part of the structure stays well defined between 28 and 35 ns. This is in line with what the RMSD plot suggested.

6.1.3) Hydrogen Bonding Analysis
Next we shall look at how we can use ptraj to monitor the presence of hydrogen bonds during our trajectory.

analyse_hbond.ptraj

trajin equil_1-10.binpos

#-- Donors from standard amino acids
donor mask :GLN@OE1
donor mask :GLN@NE2
donor mask :ASN@OD1
donor mask :ASN@ND2
donor mask :TYR@OH
donor mask :ASP@OD1
donor mask :ASP@OD2
donor mask :GLU@OE1
donor mask :GLU@OE2
donor mask :SER@OG
donor mask :THR@OG1
donor mask :HIS@ND1
donor mask :HIE@ND1
donor mask :HID@NE2

#-- Acceptors from standard amino acids
acceptor mask  :ASN@ND2 :ASN@HD21
acceptor mask  :ASN@ND2 :ASN@HD22
acceptor mask  :TYR@OH  :TYR@HH
acceptor mask  :GLN@NE2 :GLN@HE21
acceptor mask  :GLN@NE2 :GLN@HE22
acceptor mask  :TRP@NE1 :TRP@HE1
acceptor mask  :LYS@NZ  :LYS@HZ1
acceptor mask  :LYS@NZ  :LYS@HZ2
acceptor mask  :LYS@NZ  :LYS@HZ3
acceptor mask  :SER@OG  :SER@HG
acceptor mask  :THR@OG1 :THR@HG1
acceptor mask  :ARG@NH2 :ARG@HH21
acceptor mask  :ARG@NH2 :ARG@HH22
acceptor mask  :ARG@NH1 :ARG@HH11
acceptor mask  :ARG@NH1 :ARG@HH12
acceptor mask  :ARG@NE  :ARG@HE
acceptor mask  :HIS@NE2 :HIS@HE2
acceptor mask  :HIE@NE2 :HIE@HE2
acceptor mask  :HID@ND1 :HID@HD1
acceptor mask  :HIP@ND1,NE2 :HIP@HE2,HD1

#-- Backbone donors and acceptors for this particular molecule
#   N-H for prolines do not exist so are not in the mask
#  
donor mask @O
acceptor mask  :2-11,13-16,20@N :2-20@H
#Terminal residues have different atom names
donor mask @OXT
acceptor mask :1@N :1@H1
acceptor mask :1@N :1@H2
acceptor mask :1@N :1@H3

#-- series hbt is just a placeholder to ensure we get the full analysis. If you don't
#have the word series you don't get a full analysis.
hbond print .05 series hbt

This is quite a complex ptraj script so some explanation is required. The first section lists the standard amino acid side chain H bond donors and acceptors. This should be included in any ptraj H bond script if you are looking for all the hydrogen bonds. Note, only H bonds in residues defined here will be found. Hence if you have a co-enzyme in your system you will need to manually add its donors and acceptors to this list.

The next section (#-- Backbone donors...) lists which backbone atoms can act as donors and acceptors. This is in addition to the specific amino acid donors and acceptors listed before this section. Backbone atoms are tricky for a number of reasons. Firstly prolines do not have protonated backbone nitrogens:

and so must be excluded from the acceptor mask. Secondly the N and C termini have different atom naming conventions and also extra atoms, more protons and an OXT atom at the C terminus. These considerations have to be made for each specific system. I.e. the line

print .05 means "print only hydrogen bonds that are resident for more than 5% of the simulation."

series hbt "tells ptraj to print out information related to the time series. In other words we get a nice graphical plot of where in the trajectory specific H-bonds exist. hbt is just a name for this series, in this case just a temporary container."

WARNING: H-bond analysis can be very memory intensive. For example this analysis on just a 20 amino acid peptide chain with no water requires 622 MB of memory. On a P4 3.0GHz machine with 1 GB of ram it takes around 1 minute to run. If you have less than 1 GB of ram it may take considerably longer.

Before you look at the actual result always check to make sure your selections actually worked and gave you the number of atoms you were expecting. Pay special attention to any selections that are ignored. E.g.

In this case the warning is benign since we have no histidine residues present in our TRPcage structure.

If we are certain things ran as we expected we can take a look at the results. These are printed at the very end of the output file. In our case it reads:

This table should be reasonably self explanatory. An interesting section is that on the far right of the table which shows the percentage occupancy during each 10th (5 ns) of the trajectory file. You should immediately see that there are almost no hydrogen bonds during the first 15 ns (30%) of the trajectory. Those hydrogen bonds that do exist in this first section are typically not present in the remainder of the trajectory. Take for example the two hydrogen bond pairs that are only present during the first 15 ns of the trajectory. These are atoms 226 - 94(GLY15:O-TRP6:H) and 176 - 104(GLY10:O-TRP6:HE). You can visualise these in VMD. Fire up VMD and load the TC5b.prmtop and equil_1-10.binpos. It will load quicker if you turn off the display of the molecule (double click the D by it's filename in "VMD Main") and also minimise VMD.

Change to the first structure in the trajectory. Then go to the OpenGL display, press "2" to change to "bond label" mode. Then click the two pairs of H-O atoms so you get their distances labelled.

Now play through the trajectory and watch what happens. You may want to change the step size to 10 so that it steps through the trajectory at a reasonable speed. Here is an MPEG movie of the trajectory (TRPcage_h_bonds.mpg 14.5 MB). Notice how these two hydrogen bonds quickly form and remain formed for the first 15 ns of the simulation. These two bonds act to stop the alpha helix forming. Eventually they break and the alpha helix forms very rapidly. These two hydrogen bonds never reform after this point. The first major secondary structure motif to form is the alpha helix, this is in line with what they describe in the paper only they see helix formation much earlier than we do. In our simulation the formation of the two hydrogen bonds discussed above is sufficient to keep us trapped in a high energy minimum for many nano seconds, so delaying the overall folding.

Lets take a look at some of the other hydrogen bonds during the last 30 ns of the trajectory. They mention two unusual ones in the paper:

| 250 :16@O | 104 :6@HE1 103 :6@NE1 | 19.92 2.850 ( 0.09) 22.71 (11.50) 2.9 ( 2.4) 19 | -o-x- |
| 116 :6@O | 178 :11@H 177 :11@N | 14.42 2.885 ( 0.08) 28.88 (14.13) 1.5 ( 0.9) 11 | .-..--.|

However they are only present for about 20 % and 14 % of the trajectory respectively and appear to be transient. The first one disappearing at the end and the second one while present is obviously not there all the time. This suggests that these hydrogen bonds are not as strong or as well formed as in the structure the system folds to in the paper where they quote populations of 92 and 72 % respectively. The difference in percentages here may be due to us using a tighter criterion in the ptraj analysis than they used in the paper but more than likely it is due to our structure not properly folding.

Note we also see another H bond formed with the same donor as the second H bond they refer to in the paper. However, our acceptor is the neighbouring glycine residue (res 10). This is more or less uniformly present during the last 30 ns of the trajectory and may indicate that the Tryptophan carbonyl is between the two glycine residues rather than locked primarily to the second one. You can check this using VMD if you want. This may be a symptom of our Tryptophan being in a different orientation due to the misfolding.

| 116 :6@O | 171 :10@H 170 :10@N | 21.21 2.880 ( 0.08) 27.43 (13.53) 1.7 ( 1.4) 28 | .-------|

6.1.4) Dihedral Analysis using Ptraj
If you compare our theoretical folded structure with the NMR one an obvious discrepancy is in the position of the Tryptophan. In particular in the position of the tryptophan ring with respect to the protein backbone. Although we created our starting structure in what I believe is an identical method to how they did in the paper it is always possible that we did not create the same starting structure. While the final folded structure should be independent of the starting structure the timescale on which the folding occurs may not be independent. Unfortunately we don't have time in this tutorial to try running several simulations from different starting structures to test this. However, one thing we can do is take a look at the various backbone related dihedral angles involving the Tryptophan. From this information we may be able to make some predictions about whether the folding pathway is likely to be sensitive to the initial positioning of the Tryptophan.

measure_trp_angles.ptraj

trajin equil_1-10.binpos
dihedral trpphi :5@C :6@N :6@CA :6@C out trp_phi.dat
dihedral trppsi :6@N :6@CA :6@C :7@N out trp_psi.dat
dihedral trpchi1 :6@N :6@CA :6@CB :6@CG out trp_chi1.dat
dihedral trpchi2 :6@CA :6@CB :6@CG :6@CD2 out trp_chi2.dat

You may need to go through the files and manually adjust any sections where the angle has jumped from -180 to +180 so that it looks clean when plotted, I won't go into details on how to do this here but suffice to say you can do it in xmgrace with some code like:

This plot is very informative. Notice how the phi angle stays more or less constant during the entire simulation. Secondly notice how psi starts out around 160 degrees and then drops sharply to around -50 degress after about 12.5 ns. This corresponds to the point where the alpha helix starts to form. Hence we may get very different behaviour if we start our simulation with a structure where the values of the phi and psi angles are reversed. Perhaps this is what they had in the paper? If you have the inclination, and the computer time available to you, you could try this and see what results you obtain.

The chi1 and chi2 angles are also very informative. Here is the RMSD plot again, compare it to the chi1 and chi2 plots above:

Notice how the chi1 and chi2 angle changes appear to be correlated with the changes in RMSD. This suggests that changes in the Tryptophan structure are central to the stability of each of the structures sampled by the TRPcage peptide. This correlation could be investigated further by calculating a cross correlation between the RMSD and chi1 and chi2 angles. Perhaps you could try this out if you are interested. Numerical recipes contains example code for calculating correlation functions.

The final type of analysis we will do in this tutorial is called cluster analysis. Here we will use an algorithm that goes through the trajectory and locates clusters of structures that are the same. From this it creates centroids describing each cluster and then gives an RMSD for each structure in the trajectory with respect to each identified cluster.

In order to complete this section of the tutorial you will require the MMTSB toolset. This is not part of AMBER but can be downloaded free of charge from http://mmtsb.scripps.edu/software/mmtsbtoolset.html. If you are running this tutorial at an AMBER workshop then the MMTSB toolset should already be installed on your machine. You can check by seeing if the environment variable MMTSBDIR is set. E.g.

If it isn't set then you probably need to install the MMTSB toolset. Refer to the documentation on the website. Essentially you should untar the toolset to somewhere like /usr/local/mmtsb, set MMTSBDIR to this directory and then add $MMTSBDIR/bin and $MMTSBDIR/perl to your path. Once that is done we can continue.

WARNING: Do not try this clustering over a remote filesystem or you will be waiting a very long time!!! If your current working directory is on an NFS share I highly suggest that you copy all of the files you'll need to a local scratch directory and then copy the results back when you are done. Also, avoid doing an 'ls' in the directory containing the pdb files or you'll end up waiting forever for the list to scroll past.

The clustering process creates a large number of files, namely a pdb file for every frame in the trajectory. As such we don't want to 'pollute' our working directory with lots of files so we shall create some sub directories to do the processing in:

This will give you 50,000 pdb files in the PDBfit directory. Each one representing a single snapshot from the trajectory that has been RMSD fitted to residues 2 to 20 of the first structure of the NMR data.

This step is not strictly necessary but it will make it easier to locate specific pdb files later on.

Ptraj does not write leading zeros on the pdb names which means when you sort the directory it is listed as:

Rather than numerically. The following c shell script will rename the files with leading zeros:

fix_numbering_pdb.csh

#! /bin/csh
set ff = *.pdb.1
set fnam = $ff:r
set numfil = `ls -1 *.pdb.???? | wc -l`
if( $numfil != 0)then
  foreach fnam (*.pdb.????)
     set fr=$fnam:r
     set fnum=$fnam:e
     mv $fnam $fr.0$fnum
     echo $fnam $fr.0$fnum
  end
endif
set numfil = `ls -1 *.pdb.??? | wc -l`
if( $numfil != 0)then
  foreach fnam (*.pdb.???)
     set fr=$fnam:r
     set fnum=$fnam:e
     mv $fnam $fr.00$fnum
     echo $fnam $fr.00$fnum
  end
endif
set numfil = `ls -1 *.pdb.?? | wc -l`
if( $numfil != 0)then
  foreach fnam (*.pdb.??)
     set fr=$fnam:r
     set fnum=$fnam:e
     mv $fnam $fr.000$fnum
     echo $fnam $fr.000$fnum
  end
endif
set numfil = `ls -1 *.pdb.? | wc -l`
if( $numfil != 0)then
  foreach fnam (*.pdb.?)
     set fr=$fnam:r
     set fnum=$fnam:e
     mv $fnam $fr.0000$fnum
     echo $fnam $fr.0000$fnum
  end
endif

To do the clustering we will use the mmtsb tool kclust. This tool requires a file with the list of structures to cluster. In this case we will call the file clustfils. We will fill it with all of the structures from our trajectory. If you only wanted to do the cluster analysis on a portion of the trajectory you could of course only include the structures you are interested in within this file. Lets create the file:

Remove first line since kclust can only handle < 50000 structures. We will fit structures 2 to 50,000. We need to remember that we did this as it means our pdbs actually start counting from two. This means that the ID's in the centroid files will be off by one. You will see what I mean later. Here is the clustfils file.

The clustering command should take around 3 to 5 minutes to run. However, it is very very memory intensive and may take considerably longer on your machine if you do not have much memory available. In fact, if you run out of memory and do not have enough swap space, the analysis will not work at all. It seems to take a bit less than 1 Gb. On a P4/3GHz/512 Gb RAM, the machine pages nearly unceasingly, the entire run taking 10 minutes. If your machine is paging (watch a monitor like xosview), be sure to do NOTHING else on the machine while the job runs, even moving the mouse around. This will save the machine from having to page your application in and out of swap at the same time as running and paging your clustering job.

The main output from this analysis so far is the file Centroid_6 (Centroid_6.gz 478kB, don't forget to unzip this file if you download it.) which contains:
- the centroids in pdb format
- the members of each centroid
- the rmsd distance of each member from its centroid.

Take a quick look at the Centroid_6 file in a text editor so you can see what the format looks like. Note the offset in the ID and pdb numbering:
97 trp.pdb.00098 5.221610
ID 97 is actually pdb file number 98. Be careful of this.

The centroids themselves are not physically meaningful since they are effectively a mathematical construct based on the member of a cluster. However, the structure that has the lowest RMSD to each centroid is meaningful and much easier to understand. In this way we can look at the structure that is closest to the centroid that represents each of the clusters found.

extract_centroids.awk

BEGIN{b0=2;}

{centind=index($1,"#Centroid");
   c=$2;
   getline;centind=index($0,"#Centroid");
   FIL0 = sprintf("centroid%2.2d.member.dat",c)
   while(centind != 1){
     print $1,$3 > FIL0 ;
     getline;centind=index($0,"#Centroid");
   }

   numcent=0;    print $2,$1,NR-b0; 
   c=$2;
   getline; endrec = index("End",$0);
   while( endrec != 1 ){
     FIL = sprintf("centroid%2.2d.pdb",c);
     print > FIL;
     getline; endrec = index($0,"#End");
   }
   b0=NR+2;
 }

Note: The tee command allows us to pipe the output to a file Centroid_6.stats and to standard out (the screen) at the same time.

This command produces a pdb file for each centroid, a member file for each centroid and a stats file with the number of members in each centroid. You can download all these files as a tar archive here: centroid_data.tar.gz (326 kB)

1 #Centroid 11732
2 #Centroid 6179
3 #Centroid 26772
4 #Centroid 2735
5 #Centroid 98
6 #Centroid 2483

Two of them (1 & 3) are highly populated, three are marginally populated (2,4 & 6) and one is probably insignificant (5).

Individual pdb (centroid??.pdb) and member files (centroid??.member.dat) for each centroid are produced. View the pdb files in VMD but keep in mind that these are much like an average structure and do not represent a physically meaningful configuration, they are only the centroid of the cluster.

You can find the structure that is nearest the cluster centre (centroid), so you can see a more physically relevant structure of the centroid. Sort the member file numerically on the second column (rms distance from the centroid).

So structure ID 21340 has the lowest RMSD to centroid3. So we will copy this to make our bestmember for this centroid. Remember though that our ID and pdb file are offset by one since we deleted the first pdb file from the data set. Hence structure ID 21340 refers to pdb file trp.pdb.21341:

Now view this best member pdb analyzing for features of interest such as the results of the H-bond analysis.

Repeat above for the other clusters. Here is a tar file containing the 6 best members: best_members.tar.gz

The member files for each of the centroids can be read into xmgrace or xmgr and plotted. Some work is required to cleanup the plot to make it easy to read. This can be done by loading each member data set in the following order:

(the order of most populated):
centroid03.member.dat
centroid01.member.dat
centroid02.member.dat
centroid04.member.dat
centroid06.member.dat
centroid05.member.dat

Bring up Set Appearance window from xmgrace Plot menu
Select all sets
set Line properties Type to None
set Symbol properties to Type circle and adjust Size to 16
click Apply (at bottom)

Select each set, and change the Symbol colors as follows (press Apply
after selecting each color):

G0.S0 black (leave it alone, it is already black)
G0.S1 red
G0.S2 blue
G0.S3 green
G0.S4 yellow
G0.S5 orange

Some significant patterns appear in this clustering analysis. Note that cluster 6 is only at the very start of the simulation and can be ignored from now on. Note the progression of one cluster to another, where blue is the only transition structure to black except a short transition via green. The initial meta-stable cluster is the red cluster and the transition to the dominant cluster is through the clusters green, yellow, and blue.

From here one can do many other clusterings by varying the clustering parameters and the range of structures over which the clustering is done. Once a region of the trajectory is selected as interesting, clustering can be done on just that region by creating a list of pdb files (clustfils) contained in that region.

Secondly, one can change the radius for clustering, smaller radii produce more (smaller) clusters, often unmanageable and insignificantly different from each other. However, without some trial and error, one does not know what a productive and informative value of the radius should be.

Thirdly, one can explore the clustering of pdb files that have been aligned in different ways. Our alignment was to residues 2 to 20 of the first structure. One could try aligning to another part of the structure, or another structure such as the lowest energy or even one of the closest members of the first clustering analysis.

One more analysis of the clusters produced above is to process the clusters through ptraj to find out other information about the clusters. The simplest analysis is rmsd to some structure such as the centroid of the cluster, but we probably will not gain any new information since the clustering produced that information anyway.

The best approach would be to create trajectories containing only the members of each cluster and doing analysis on each trajectory such as rms, and hbond. I leave this step to the more enthusiastic tutee.