Hi,
So if I understand you correctly, when using 2 topologies you see the
issue, but with 1 topology the issue goes away? If so then this is
probably a bug - I'll start looking into it immediately. Thanks for
the report. Let me know if I missed anything.
-Dan
On Wed, Nov 26, 2014 at 2:58 AM, Guillaume Roellinger <roellinger.bio.mx> wrote:
> Hi,
> Thanks for your quick answer and your help. I switched to covar and get also different eigenvalue results for both cases which seems to be consistent with the values obtained with mwcovar:
> - with the reference structures loaded as first:
> '792.76199', '45.47370', '16.99601', '9.82698', '8.86663
> - with the MD trajectories loaded as first:
> '0.83129', '0.00000', '0.00000', '0.00000', '0.00000'
>
> When doing a diff of both topologies, I get a slight difference in addition to the version date, which corresponds to the box dimensions:
> 1c1
> < %VERSION VERSION_STAMP = V0001.000 DATE = 11/26/14 08:52:57
> ---
>> %VERSION VERSION_STAMP = V0001.000 DATE = 11/26/14 08:54:07
> 23251c23251
> < 1.09471219E+02 8.23027770E+01 8.23027770E+01 8.23027770E+01
> ---
>> 1.09471219E+02 8.48211916E+01 8.48211916E+01 8.48211916E+01
>
> As specified in the current case, both structures have the same number of residues and the same sequence, so -as you said- just one topology is necessary for both conformations. However I am also working on other cases where some residues are mutated or the sequence of one conformation is longer than the other (and use masks to restrict the analysis to the sequence portion common to both conformations) and then it needs two different topologies.
>
> When loading just one topology (conf1.prmtop or conf2.prmtop) also with covar, I get the exact same result for the eigenvalue in both cases:
> '798.22756', '70.64851', '45.25603', '27.81210', '16.52141'
>
> I duplicated conf1.prmtop and renamed it in conf2.prmtop and once again I get two different results when the order of the trajin commands is changed (exactly the same results as described at the beginning of this message).
>
> Thanks for your help,
>
> G.R.
> ________________________________________
> Von: Daniel Roe <daniel.r.roe.gmail.com>
> Gesendet: Dienstag, 25. November 2014 19:08
> An: AMBER Mailing List
> Betreff: Re: [AMBER] Different PCA results when loading trajectories in a different order
>
> PS - Also, if this is the same protein but different conformations you
> should only need 1 topology (in fact if this is the case then
> conf1.prmtop and conf2.prmtop should only differ by their VERSION
> line). The only time you need a different topology file for two
> trajectories is if the atoms/residues do not match between the two
> (such as your PDB that only contains CA atoms).
>
> On Tue, Nov 25, 2014 at 11:03 AM, Daniel Roe <daniel.r.roe.gmail.com> wrote:
>> Hi,
>>
>> Out of curiosity do you see the same thing happen if you use
>> non-mass-weighted covariance matrices (i.e. covar instead of mwcovar)?
>> Note that typically in PCA just plain covariance matrices are used.
>> Mass-weighted covariance matrices are really only required if you want
>> to perform quasi-harmonic analysis.
>>
>> -Dan
>>
>> On Tue, Nov 25, 2014 at 10:49 AM, Guillaume Roellinger
>> <roellinger.bio.mx> wrote:
>>> Dear Amber users,
>>>
>>> Using cpptraj, I am trying to analyze two 10ns MDs in one principal component analysis (PCA). Both MDs are closely related to each other: it is the same protein but in two different conformations (conf1 and conf2) with same number of residues. I was interested to include the structures of each conformations just after equilibration as reference structures for the PCA. I autoimaged and removed water and ions from my trajectories and reference structures and generated the corresponding stripped topologies conf1.prmtop, conf2.prmtop for each conformation. Before calculating the mass-weighted matrix, I aligned all frames to an average structure generated from both trajectories (createAverageStructure_cpptraj file).
>>>
>>> The problem is that I get two very different results in my evecs file when switching the loading order between the reference structures and trajectories and I have no idea why.
>>> By loading the reference structure first (PCA_cpptraj_1), I get the following 5 first eigenvalues:
>>> 0.85663
>>> 3.57670
>>> 5.85045
>>> 7.69400
>>> 8.09997
>>> ...(and so on, increasing linearly)
>>>
>>> And by loading the MD trajectories first (PCA_cpptraj_2):
>>> 26.45372
>>> 2640022.22255
>>> 2642822.92226
>>> 2769903.48459
>>> 2844226.43589
>>> ...(and so on, continuing with values ~2e6, 3e6)
>>>
>>> If I try not to load any reference structures, I get almost the same results as when loading the reference structure first (PCA_cpptraj_1); so I checked the reference structures and they are both fine. I checked also the structures after the RMSD on the average structure with outtraj and in both cases the structures are correctly aligned.
>>> If necessary, I can provide a 12MB archive containing the input files and shortened MDs trajectories (20 first frames of each) which is sufficient to reproduce this behavior.
>>>
>>> I am using cpptraj V14.09 and calling cpptraj via $AMBERHOME/bin/cpptraj -i file_cpptraj
>>>
>>> Thanks for your help,
>>>
>>> Best regards,
>>>
>>> G.R.
>>>
>>> Input Files used:
>>>
>>> createAverageStructure_cpptraj:
>>> ------------------------------------------------------------
>>> parm conf1.prmtop
>>> parm conf2.prmtop?
>>> trajin conf1_MDtraj.mdcrd parmindex 0
>>> trajin conf2_MDtraj.mdcrd parmindex 1?
>>> rms RMSDfirst first out RMSDfirst.dat .CA
>>> average averageCA.pdb .CA
>>> ------------------------------------------------------------
>>>
>>>
>>> PCA_cpptraj_1 (loading reference structures first):
>>> ------------------------------------------------------------
>>> parm conf1.prmtop
>>> parm conf2.prmtop
>>> trajin conf1_afterEquil.rst7 parmindex 0
>>> trajin conf2_afterEquil.rst7 parmindex 1
>>> trajin conf1_MDtraj.mdcrd parmindex 0 # 500 frames
>>> trajin conf2_MDtraj.mdcrd parmindex 1 # 500 frames
>>> parm averageCA.prmtop
>>> reference averageCA.pdb parmindex 2
>>> rms RMSDaverage reference out RMSDaverage.dat .CA
>>> matrix mwcovar name mwcvmat out mwcvmat.dat .CA
>>> diagmatrix mwcvmat out evecs.dat vecs 200 name mymodes
>>> run
>>> readdata evecs.dat
>>> projection modes evecs.dat out project.dat beg 1 end 50 .CA
>>> ------------------------------------------------------------
>>>
>>> PCA_cpptraj_2 (loading MD first):
>>> ------------------------------------------------------------
>>> parm conf1.prmtop
>>> parm conf2.prmtop
>>> trajin conf1_MDtraj.mdcrd parmindex 0 # 500 frames
>>> trajin conf2_MDtraj.mdcrd parmindex 1 # 500 frames
>>> trajin conf1_afterEquil.rst7 parmindex 0
>>> trajin conf2_afterEquil.rst7 parmindex 1
>>> parm averageCA.prmtop
>>> reference averageCA.pdb parmindex 2
>>> rms RMSDaverage reference out RMSDaverage.dat .CA
>>> matrix mwcovar name mwcvmat out mwcvmat.dat .CA
>>> diagmatrix mwcvmat out evecs.dat vecs 200 name mymodes
>>> run
>>> readdata evecs.dat
>>> projection modes evecs.dat out project.dat beg 1 end 50 .CA
>>> ------------------------------------------------------------
>>> _______________________________________________
>>> 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 307
>> 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 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
>
> _______________________________________________
> 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 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 Nov 26 2014 - 08:00:03 PST
