Re: [AMBER] PCA eigenvalues

From: Nisha Amarnath Jonniya <phd1601271002.iiti.ac.in>
Date: Fri, 24 Nov 2017 17:09:01 +0530

Dear all,

    As I was also doing PCA analysis. This mail helped me too but I got
error in execution of hist: command .

These are my script.

parm complex_solv.prmtop
trajin trajout10-51.mdcrd
rms first :1-273&!.H=
average crdset ex_average
createcrd ex_trajectories
run
crdaction ex_trajectories rms ref ex_average :1-273&!.H=
crdaction ex_trajectories matrix covar name ex_covar :1-273&!.H=
runanalysis diagmatrix ex_covar out evececs 20 name exEvecs:1-273&!.H=
readdata evececs name Evecs
runanalysis modes eigenval name Evecs out evalfrac.dat

crdaction ex_trajectories projection ex modes Evecs beg 1 end 3
:1-273&!.H= crdframes 84000

hist ex:1 bins 100 out 5drb_hist.gnu norm name ex_1
hist ex:2 bins 100 out 5drb_hist.gnu norm name ex_2
hist ex:3 bins 100 out 5drb_hist.gnu norm name ex_3
run
readdata evececs name Evecs
parm complex_solv.prmtop
parmstrip !(:1-273&!.H=)
parmwrite out ex_modes.prmtop
runanalysis modes name Evecs trajout ex_mode1.nc pcmin -100 pcmax 100 tmode
1 trajoutmask :1-273&!.H= trajoutfmt netcdf
runanalysis modes name Evecs trajout ex_mode2.nc pcmin -100 pcmax 100 tmode
2 trajoutmask :1-273&!.H= trajoutfmt netcdf
runanalysis modes name Evecs trajout ex_mode3.nc pcmin -100 pcmax 100 tmode
3 trajoutmask :1-273&!.H= trajoutfmt netcdf



I am getting error in hist command with following error


 Data set 'ex:1' not found.
Error: Dataset ex:1 not found.
Error: Could not setup analysis [hist]
    1 errors encountered reading input.
Warning: File 'ex_hist.gnu' has no sets containing data.




Can you please guide me further how to rectify this problem.






On Fri, Nov 24, 2017 at 12:26 AM, Lizelle Lubbe <LBBLIZ002.myuct.ac.za>
wrote:

> Hi Dan,
>
> Thanks so much for the great advice you sent earlier.
> It worked out well but I noticed the following warning in the output:
> Warning: Frame 25 Coordinates out of bounds (100). What does it mean and
> should I be worried about it?
>
> The projection of PC1 vs PC2 was written as a gnu file like you suggested
> and displays a scatter of frames. However, the plot has a black background
> and purple dots for each frame which is really not very clear.
> Is there any way that I can change the colour or write it in another
> format? I can't seem to edit the plot at all in gnuplot.
>
> Regards
>
> Lizelle
>
> ________________________________________
> From: Daniel Roe <daniel.r.roe.gmail.com>
> Sent: 20 November 2017 05:53:21 PM
> To: AMBER Mailing List
> Subject: Re: [AMBER] PCA eigenvalues
>
> s.dat vHi,
>
> On Sat, Nov 18, 2017 at 10:05 AM, Lizelle Lubbe <LBBLIZ002.myuct.ac.za>
> wrote:
> > Hi Amber users,
> >
> > I have performed PCA analysis on my MD trajectories in cpptraj and would
> like to plot the fraction of total variance vs eigenvalue as shown in
> figure s9 of Roe et al 2014 (j phys chem b, 118, 3543-3552).
>
> To do this you'll want to calculate all principal components (i.e.
> eigenvectors) from your covariance matrix, so 3 * the number of
> selected atoms. Then you run the 'modes' analysis command with the
> 'eigenval' keyword on the resulting modes data, e.g.
>
> runanalysis diagmatrix myCovar out evececs 246 name myEvecs :1-4&!.H=
> runanalysis modes eigenval name myEvecs out evalfrac.dat
>
> > I also want to have a scatter plot of pc1 vs pc2 and pc2 vs pc3 to
> visualize clusters in this space.
> > Could someone please guide me in how this can be accomplished?
>
> For these, just give desired the principal component projection data
> set names to the 'hist' command, e.g. for the projection of PC 1 vs 2:
>
> hist Ndom:1 Ndom:2 bins 100 out Ndom_hist.gnu norm name Ndom_1_2
>
> > Also, how should the pca histogram analysis be interpreted?
>
> This is a lot trickier. As you probably know principal components
> represent axes that in this case describe variance in the selected
> atoms in your system. So PCs can tell you about what the dominant
> modes of motion in your system look like, but generally may not tell
> you how your system actually moves. Personally I like to visualize the
> modes of motion with the 'trajout' keyword of 'modes' (which it
> appears you are doing). It can also help to look at the so-called
> 'porcupine plots' - I use the nmwiz plugin for VMD (it appears you are
> doing that as well).
>
> Hope this helps,
>
> -Dan
>
> >
> > This is my pca script:
> >
> > trajin traj.nc
> >
> > rms first :1-720&!.H=
> >
> > average crdset Ndom_average
> >
> > createcrd Ndom_trajectories
> >
> > run
> >
> > crdaction Ndom_trajectories rms ref Ndom_average :1-720&!.H=
> >
> > crdaction Ndom_trajectories matrix covar name Ndom_covar :1-720&!.H=
> >
> > runanalysis diagmatrix Ndom_covar out Ndom_evecs.dat vecs 20 name
> NdomEvecs nmwiz nmwizvecs 3 nmwizfile Ndom_pca.nmd nmwizmask :1-720&!.H=
> >
> > crdaction Ndom_trajectories projection Ndom modes NdomEvecs beg 1 end 3
> :1-720&!.H= crdframes 1,3000
> >
> > hist Ndom:1 bins 100 out Ndom_hist.agr norm name Ndom_1
> >
> > hist Ndom:2 bins 100 out Ndom_hist.agr norm name Ndom_2
> >
> > hist Ndom:3 bins 100 out Ndom_hist.agr norm name Ndom_3
> >
> > run
> >
> > clear all
> >
> > readdata Ndom_evecs.dat name Evecs
> >
> > parm NDOM_CPLX_rnb_stripped_box.prmtop
> >
> > parmstrip !(:1-720&!.H=)
> >
> > parmwrite out Ndom_modes.prmtop
> >
> > runanalysis modes name Evecs trajout Ndom_mode1.nc pcmin -100 pcmax 100
> tmode 1 trajoutmask :1-720&!.H= trajoutfmt netcdf
> >
> > runanalysis modes name Evecs trajout Ndom_mode2.nc pcmin -100 pcmax 100
> tmode 2 trajoutmask :1-720&!.H= trajoutfmt netcdf
> >
> > runanalysis modes name Evecs trajout Ndom_mode3.nc pcmin -100 pcmax 100
> tmode 3 trajoutmask :1-720&!.H= trajoutfmt netcdf
> >
> > Kind regards
> >
> > Lizelle Lubbe
> > PhD Chemical Biology candidate
> > University of Cape Town
> > Disclaimer - University of Cape Town This e-mail is subject to UCT
> policies and e-mail disclaimer published on our website at
> http://www.uct.ac.za/about/policies/emaildisclaimer/ or obtainable from
> +27 21 650 9111. If this e-mail is not related to the business of UCT, it
> is sent by the sender in an individual capacity. Please report security
> incidents or abuse via csirt.uct.ac.za
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
>
>
>
> --
> -------------------------
> Daniel R. Roe
> Laboratory of Computational Biology
> National Institutes of Health, NHLBI
> 5635 Fishers Ln, Rm T900
> Rockville MD, 20852
> https://www.lobos.nih.gov/lcb
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
> Disclaimer - University of Cape Town This e-mail is subject to UCT
> policies and e-mail disclaimer published on our website at
> http://www.uct.ac.za/about/policies/emaildisclaimer/ or obtainable from
> +27 21 650 9111. If this e-mail is not related to the business of UCT, it
> is sent by the sender in an individual capacity. Please report security
> incidents or abuse via csirt.uct.ac.za
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>



-- 
Nisha Amarnath Jonniya
PhD Research Scholar
Biosciences and Biomedical Engineering
Indian Institute of Technology, Indore
India
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Nov 24 2017 - 04:00:03 PST
Custom Search