Re: [AMBER] Principal Component Analysis

From: Daniel Roe <daniel.r.roe.gmail.com>
Date: Wed, 2 Dec 2015 11:00:49 -0700

A couple of comments on your input script:

On Wed, Dec 2, 2015 at 3:23 AM, Beale, John <John.Beale.stlcop.edu> wrote:
> trajin nacamd.nc 1 last 1
> reference nac.inpcrd
> center origin :1-35
> image origin center

The center and image commands are not relevant to the analysis you are
performing here and will only slow things down.

> strip :WAT,Na+
> rms reference mass :1-35.CA,C,N out RMSD.out
> matrix covar name matrixdat .CA out covmat-ca.dat
> analyze matrix matrixdat out evecs-ca.dat vecs 5
> analyze matrix matrixdat name evecs-ca vecs 5

As Jason mentioned, there is no reason to do this twice. If you want
to both write out the eigenvalues/eigenvectors and save them
internally as a data set just put all the keywords with a single
command (also I suggest that you use 'diagmatrix' instead of 'analyze
matrix'), e.g.

diagmatrix matrixdat out evecs-ca.dat name evecs-ca vecs 5

> analyze modes fluct out analyzemodesflu.dat name evecs-ca beg 1 end 5
> analyze modes displ out analyzemodesdis.dat name evecs-ca beg 1 end 5

FYI, you don't need the 'analyze' prefix for the 'modes' analysis command.

> projection modes evecs-ca.dat out pca12-ca beg 1 end 3 .CA

I suspect this caused you problems the first time you ran it. The
issue is that by default cpptraj (like ptraj before it) queues up all
action and analysis commands until trajectory processing (i.e. a
'Run') is initiated, either by the end of the script or by 'go'/'run'.
The actions are executed during a Run, and analyses are executed after
the Run is complete. This means that with this input your 'projection'
command is actually being executed just after the 'matrix' command.
Since you need to first calculate your eigenvectors before can project
your coordinates along them, this 'project' command will fail the
first time because 'evecs-ca.dat' hasn't been created yet. If you want
to make sure things happen in the order you want, you will have to put
a 'run' command before 'projection' to make sure the eigenvectors are
calculated before the 'projection' command.

If I were going to write the script I'd do it like so:

trajin nacamd.nc
reference nac.inpcrd
strip :WAT,Na+
rms reference mass :1-25.CA,C,N out RMSD.out
matrix covar name matrixdat .CA out covmat-ca.dat
# Save the RMS-fit coordinates as a data set for later projection
createcrd MyCrd
# Calculate eigenvectors from covariance matrix
diagmatrix matrixdat out evecs-ca.dat name evecs-ca vecs 5
# Run, first pass
run
# Project coordinates from MyCrd on eigenvector data set
crdaction MyCrd projection MyProjection evecs evecs-ca beg 1 end 3 .CA
# Histogram first 2 PC projections
runanalysis hist MyProjection:1 MyProjection:2 out PC_1_2.gnu bins 200

For more examples on using cpptraj to perform PCA see scripts in the
supporting info of http://pubs.acs.org/doi/abs/10.1021/ct400862k and
http://pubs.acs.org/doi/abs/10.1021/jp4125099.

Hope this helps,

-Dan

-- 
-------------------------
Daniel R. Roe, PhD
Department of Medicinal Chemistry
University of Utah
30 South 2000 East, Room 307
Salt Lake City, UT 84112-5820
http://home.chpc.utah.edu/~cheatham/
(801) 587-9652
(801) 585-6208 (Fax)
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Dec 02 2015 - 10:30:03 PST
Custom Search