Re: [AMBER] Running PCA

From: Nicole Ippolito <ippolito.nicole89.gmail.com>
Date: Tue, 29 Jul 2014 10:07:55 -0500

Thank you for the detailed explanation and the reference. I believe it is
working.

Nicole


On Tue, Jul 29, 2014 at 9:31 AM, Daniel Roe <daniel.r.roe.gmail.com> wrote:

> PS - For a concrete example of performing PCA with cpptraj see
> http://pubs.acs.org/doi/abs/10.1021/jp4125099 and supporting material
> therein.
>
> -Dan
>
>
> On Tue, Jul 29, 2014 at 8:16 AM, Daniel Roe <daniel.r.roe.gmail.com>
> wrote:
>
> > Hi,
> >
> > You're missing some elements from your script, and have a few things out
> > of order. By default in cpptraj when you specify an action or an analysis
> > it goes into either the action list or analysis list respectively; it is
> > not immediately executed. I will go through your input and try to explain
> > what's happening:
> >
> > On Fri, Jul 25, 2014 at 11:56 AM, Nicole Ippolito <
> > ippolito.nicole89.gmail.com> wrote:
> >
> >> trajin ~/BRAF/MD/V600E505H/sor/V600E505H_sor_aMD_0_300ns.crd.gz
> >> createcrd crd1
> >> run
> >>
> >
> > Here you are essentially loading a trajectory into the coords data set
> > named 'crd1'. The 'run' statement means "process all currently loaded
> > trajectories with any initialized actions", and 'createcrd' is the only
> > action set up here.
> >
> >
> >> rms first .CA
> >> reference V600E505H_sor_aMD_avg_pdb.pdb
> >>
> >
> > These statements seem out of order. They are placed after the run
> > statement, and so they are not processed until the first 'run' is
> > completed. Presumably your intent was to remove overall rotation and
> > translation from your coordinates with the 'rms' statement. If that is
> the
> > case then it should come before the createcrd statement. Also, since you
> > also have a later 'rms' command that fits to a reference, this makes me
> > think this first 'rms' was for creating an average structure to use as a
> > reference, and the 'average' command is missing.
> >
> > The 'reference' command is actually what kills the run: if you check your
> > output you will see the message "Error: reference: Could not set up
> > trajectory". This indicates that "V600E505H_sor_aMD_avg_pdb.pdb" does not
> > exist; perhaps you meant to create an average?
> >
> >
> >> crdaction crd1 rms reference .CA
> >> crdaction crd1 matrix out matrix.dat name mymatrix covar .CA
> >>
> > runanalysis diagmatrix mymatrix out evecs.dat vecs 20 name mymodes
> >>
> >
> > These statements seem fine as long as a reference is successfully loaded
> > and your 'crd1' coordinates is properly loaded.
> >
> >
> >> crdaction crd1 projection modes evecs.dat start 1 stop 3 .CA myproj
> >>
> >
> > This will not work, since you named your modes internally as 'mymodes'
> but
> > you're trying to use the file name 'evecs.dat'. Change 'evecs.dat' to
> > 'mymodes'.
> >
> >
> >> hist myproj:1 bins 200 out PC.hist.agr
> >> hist myproj:2 bins 200 out PC.hist.agr
> >>
> >
> > This seems like it should work as well.
> >
> > So to summarize, your final script should look something like this:
> >
> > trajin ~/BRAF/MD/V600E505H/sor/V600E505H_sor_aMD_0_300ns.crd.gz
> > # Fit to first frame for averaging
> > rms first .CA
> > # Average the trajectory
> > average V600E505H_sor_aMD_avg_pdb.pdb pdb
> > # Save the rms-fit coords
> > createcrd crd1
> > # First run
> > run
> > # Load the generated average as reference
> > reference V600E505H_sor_aMD_avg_pdb.pdb
> > # Fit saved coords to the reference
> > crdaction crd1 rms reference .CA
> > # Calculate covariance matrix from fit coords
> > crdaction crd1 matrix out matrix.dat name mymatrix covar .CA
> > # Diagonalize matrix to get eigenvectors/eigenvalues
> > runanalysis diagmatrix mymatrix out evecs.dat vecs 20 name mymodes
> > # Project fit coords along eigenvectors 1-3
> > crdaction crd1 projection modes evecs.dat start 1 stop 3 .CA myproj
> > # Set up histograms for PC projections
> > hist myproj:1 bins 200 out PC.hist.agr
> > hist myproj:2 bins 200 out PC.hist.agr
> > hist myproj:3 bins 200 out PC.hist.agr
> > # Run analyses
> > runanalysis
> >
> > Note that the final 'runanalysis' is not strictly necessary as cpptraj
> > will try to run any actions or analyses that are set up before quitting
> > when reading a script, but I include it for completeness.
> >
> > Hope this helps,
> >
> > -Dan
> >
> >
> >
> >
> >>
> >>
> >> CPPTRAJ: Trajectory Analysis. V14.01
> >> ___ ___ ___ ___
> >> | \/ | \/ | \/ |
> >> _|_/\_|_/\_|_/\_|_
> >> Reading 'parm' as Amber Topology
> >> INPUT: Reading Input from STDIN
> >> [trajin ~/BRAF/MD/V600E505H/sor/V600E505H_sor_aMD_0_300ns.crd.gz ]
> >> Reading
> >> '/home/aubndi/BRAF/MD/V600E505H/sor/V600E505H_sor_aMD_0_300ns.crd.gz' as
> >> Amber Trajectory
> >> [createcrd crd1 ]
> >> CREATECRD: Saving coordinates from Top parm to "crd1"
> >> [run]
> >> ---------- RUN BEGIN -------------------------------------------------
> >>
> >> PARAMETER FILES:
> >> 0: 'parm', 36610 atoms, 11002 res, box: Orthogonal, 10731 mol, 10722
> >> solvent, 60000 frames
> >>
> >> INPUT TRAJECTORIES:
> >> 0: 'V600E505H_sor_aMD_0_300ns.crd.gz' is an AMBER trajectory, Parm parm
> >> (Orthogonal box) (reading 60000 of 60000)
> >> Coordinate processing will occur on 60000 frames.
> >> TIME: Run Initialization took 0.0000 seconds.
> >>
> >> BEGIN TRAJECTORY PROCESSING:
> >> .....................................................
> >> ACTION SETUP FOR PARM 'parm' (1 actions):
> >> 0: [createcrd crd1 ]
> >> ----- V600E505H_sor_aMD_0_300ns.crd.gz (1-60000, 1) -----
> >> 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% Error: reference: Could not set
> up
> >> trajectory.
> >> 1 errors encountered reading input.
> >> 100% Complete.
> >>
> >> Read 60000 frames and processed 60000 frames.
> >> TIME: Trajectory processing: 1743.3814 s
> >> TIME: Avg. throughput= 34.4159 frames / second.
> >>
> >> ACTION OUTPUT:
> >>
> >> DATASETS:
> >> 1 data set:
> >> crd1 "crd1" (coordinates), size is 60000 (25138 MB) Box Coords,
> >> 36610 atoms
> >> ---------- RUN END ---------------------------------------------------
> >> [rms first .CA]
> >> RMSD: (.CA), reference is first frame (.CA), with fitting.
> >> [reference V600E505H_sor_aMD_avg_pdb.pdb]
> >> TIME: Total execution time: 1743.8619 seconds.
> >> _______________________________________________
> >> AMBER mailing list
> >> AMBER.ambermd.org
> >> http://lists.ambermd.org/mailman/listinfo/amber
> >>
> >
> >
> >
> > --
> > -------------------------
> > Daniel R. Roe, PhD
> > Department of Medicinal Chemistry
> > University of Utah
> > 30 South 2000 East, Room 201
> > Salt Lake City, UT 84112-5820
> > http://home.chpc.utah.edu/~cheatham/
> > (801) 587-9652
> > (801) 585-6208 (Fax)
> >
>
>
>
> --
> -------------------------
> Daniel R. Roe, PhD
> Department of Medicinal Chemistry
> University of Utah
> 30 South 2000 East, Room 201
> 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
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Jul 29 2014 - 08:30:02 PDT
Custom Search