PCA and free energy plots

[Former user]

Hi, I was wondering if the PCA function of Bio3D can be used to make a plot of the PC and the free energy of a system as sample by MD (specifically AMBER). In many publications my group has been reading, we often see 2D plots of PC's vs energy in kcal/mol. Can we use R and Bio3D to make similar plots with our amber generated trajectories? If so, can you point us in the right direction as to how to figure out how to do this?

Thanks, Jonathan

[Jonathan Gough]

Any help or insight or direction would be appreciated. Thanks.

[Lars Skjærven]

Hi Jonathan, Could you point to an example paper (perhaps with figure number or similar) ? I'm not quite sure I follow what you aim to do. Do you want the projection of a structure along a specific PC? In that case you should consult the PCA tutorial here, where you will find something like

plot(pc.xray$z[, 1], pc.xray$z[,2])

and replace "pc.xray$z[,2]" with whatever value you have from your simulation. right?

[Jonathan Gough]

Hi Lars,

Here are two examples. The second of which cites bio3D.

http://pubs.acs.org/doi/abs/10.1021/jp401587e Figure 6

or http://pubs.acs.org/doi/abs/10.1021/ct300284c figure 2

Thanks, Jonathan

[Lars Skjærven]

Thanks. Bio3d would help you with the PCA (of your simulation), and the plotting of the projections onto the PCs (see the link to the tutorial I mentioned above). You would then need to calculate the free energy - see the supplementary material for more details on this. There is no automated function in Bio3d to do this at the moment, but that would be a nice contribution ;)

[Xinqiu Yao]

Thanks Lar, for the response. It is a good idea to add function for free energy calculation (probably including both reweighting and plot). I would like to do it in near future, once I get more time from project at hand.

Hi Jonathan,

I think you can try the function discretize2d() in the "entropy" package, which basically builds a 2D histogram for your input (PC1-PC2). Then, normalize the histogram, calculate free energy with the formula DeltaG=-k_BTlog(P/P0), where P0 and P are probabilities of some reference bin and the bin you are calculating, respectively. Results can be plotted with e.g. heatmap().

You can also consider the package developed by the same group of the references you mentioned above. http://mccammon.ucsd.edu/computing/amdReweighting/ However, it is designed specifically for aMD reweighting and so I am not sure it is also suitable for your case.

[Lars Skjærven]

marking as resolved. we can add such a plot to the todo list

[Fabrício Bracht]

Hello Xin-Qiu Yao I've been trying to get the same plot that you suggested in this topic. I've calculated the PCA using a core fitted MD simulation. I've discretized using the discretize 2d with a number of bins calculated by the Scott's choice (https://en.wikipedia.org/wiki/Histogram#Scott.27s_normal_reference_rule). I've then normalized the resulting discretize2d result by divinding it by sum(discretize2d(MyPCA....)). I'm not sure, though, how I would calculate the free energy DeltaG you've mentioned, because I am not really sure which P and which P0 to use. If it is not asking too much, could you be a little more specific on that?

Thanks Fabrício

[Xinqiu Yao]

Hi Fabrício,

It doesn't matter which P0 you choose because only the relative free energy has a meaning, which is independent from P0. For example, you can choose the minimal non-zero probability as P0 (and so all free energy should be less than or equal to 0). Note that you should take the logarithm for non-zero 'P' only, and after the conversion, set DeltaG=0 at 'P=0' regions.

Hope it helps.