计算化学公社

标题: rmsf和PCA的分析怎么做？ [打印本页]

作者
Author: 胡姐姐 时间: 2018-1-5 16:48
标题: rmsf和PCA的分析怎么做？
我做的是关于蛋白质的模拟，现在轨迹跑出来了。老师让我做一下rmsf和PCA分析，请问，大家都是怎么做出来的？纠结好久了。看网上有关于vmdICE可以进行rmsf的计算，但是window size 死活搞不懂

求大神指点啊

作者
Author: 霜晨月 时间: 2018-1-5 16:52
我一般用Gromacs，rmsf 可以用 gmx rmsf 命令，PCA的话，先用gmx covar生成本征向量和本征值，然后用gmx anaeig 进行分析。

作者
Author: 胡姐姐 时间: 2018-1-5 19:37

霜晨月发表于 2018-1-5 16:52
我一般用Gromacs，rmsf 可以用 gmx rmsf 命令，PCA的话，先用gmx covar生成本征向量和本征值，然后用gmx an ...

谢谢你啦

作者
Author: 胡姐姐 时间: 2018-1-5 21:39

霜晨月发表于 2018-1-5 16:52
我一般用Gromacs，rmsf 可以用 gmx rmsf 命令，PCA的话，先用gmx covar生成本征向量和本征值，然后用gmx an ...

请问，文件格式不对，怎么办？用什么转化呢？还是用gromacs重做一遍吗

作者
Author: 霜晨月 时间: 2018-1-5 22:29
用catdcd可以将dcd转换成gromacs所用的trr格式

作者
Author: 霜晨月 时间: 2018-1-5 22:31
我的方法可能太笨了，论坛上应该有高人知道如何直接分析NAMD的轨迹

作者
Author: 胡姐姐 时间: 2018-1-6 08:26

霜晨月发表于 2018-1-5 22:31
我的方法可能太笨了，论坛上应该有高人知道如何直接分析NAMD的轨迹

我觉得这个方法已经很好啦，还有更好的方法

作者
Author: 胡姐姐 时间: 2018-1-10 10:31

霜晨月发表于 2018-1-5 22:29
用catdcd可以将dcd转换成gromacs所用的trr格式

我现在已经将轨迹文件转换为trr形式，可是，如何将xst转化为gromacs的xtc？

作者
Author: fhh2626 时间: 2018-1-10 10:43
分子模拟软件不是黑箱，做科研也不是学会写个NAMD的输入文件就可以了。
你算一个量之前，至少要明白这个量是什么意思：
When a dynamical system fluctuates about some well-defined average position, the RMSD from the average over time can be referred to as the RMSF or root mean square fluctuation.
你如果知道RMSF就是一种RMSD，就知道用VMD自带的RMSD Visualizer点两下鼠标就可以算

PCA也一样，用NMWiz就可以算

作者
Author: 霜晨月 时间: 2018-1-10 11:58

胡姐姐发表于 2018-1-10 10:31
我现在已经将轨迹文件转换为trr形式，可是，如何将xst转化为gromacs的xtc？

Gromacs是用不到xst文件的，有trr和pdb就够了。
但是你将dcd转成trr之前，要先修正了PBC。否则，没有tpr文件，gromacs是无法修正PBC的。

作者
Author: 胡姐姐 时间: 2018-1-10 12:54

霜晨月发表于 2018-1-10 11:58
Gromacs是用不到xst文件的，有trr和pdb就够了。
但是你将dcd转成trr之前，要先修正了PBC。否则，没有tpr ...

请问，你有什么资料或文献推荐吗？大神，我还是没有听懂您说的

作者
Author: 霜晨月 时间: 2018-1-10 13:23
本帖最后由霜晨月于 2018-1-10 13:27 编辑

胡姐姐发表于 2018-1-10 12:54
请问，你有什么资料或文献推荐吗？大神，我还是没有听懂您说的

虽然我不是什么大神，但已经说得很明白了啊

用 gmx trjconv 可以把trr转成xtc。用不到xst文件的。详情可以看看gmx trjconv -h 和gmx covar、gmx anaeig的帮助内容。

不知道什么是PBC的话，可能就得详细看看NAMD或GMX的manual了。

作者
Author: 胡姐姐 时间: 2018-1-10 13:32

霜晨月发表于 2018-1-10 13:23
虽然我不是什么大神，但已经说得很明白了啊
用 gmx trjconv 可以把trr转成xtc。用不到xst文件的 ...

在我这里，您就是大神的！我知道大神说的很明白了，由于我平时知识积累的太少，所以没办法理解大神说的

，希望大神可以多多帮忙解惑，拜托啦

作者
Author: 胡姐姐 时间: 2018-1-10 13:36

fhh2626 发表于 2018-1-10 10:43
分子模拟软件不是黑箱，做科研也不是学会写个NAMD的输入文件就可以了。
你算一个量之前，至少要明白这个量 ...

那我是在算rmsf之前，先算一下整个轨迹的平均结构是吗？

作者
Author: 胡姐姐 时间: 2018-1-10 14:29

胡姐姐发表于 2018-1-10 13:32
在我这里，您就是大神的！我知道大神说的很明白了，由于我平时知识积累的太少，所以没办法理解大神说的{: ...

gmx grompp -f npt-nopr-md.mdp -c npt-pr.gro -p fws.top -o npt-nopr.tpr。  这是教程中指出的生成tpr文件。
1.我知道PBC是周期性边界条件，改的是边长和中心坐标吗？
2.npt.mdp是把续跑的这个文件
structure             all_wb_ion.psf
coordinates          all_wb_ion.pdb

set init_temp       100
set desired_temp    300
#temperature       $init_temp
set outputname       all-2nd
# Continuing a job from the restart files
binCoordinates    all-1st.restart.coor
binVelocities    all-1st.restart.vel  ;# remove the "temperature" entry if you use this!
extendedSystem    all-1st.xsc
firsttimestep    21520000
paraTypeCharmm       on
parameters          par_all36_prot.prm
parameters          toppar_water_ions.str
parameters          sub.final.str

exclude             scaled1-4
1-4scaling          1.0
cutoff                12.0
switching             on
switchdist          10.0
pairlistdist          14.0

timestep             2.0
rigidBonds          all
nonbondedFreq       1
fullElectFrequency    1
stepspercycle       20

wrapAll             on
wrapWater             on
wrapNearest          on

PME                   yes
PMEGridSizeX          72
PMEGridSizeY          72
PMEGridSizeZ          72

outputName          $outputname

restartfreq          2000    ;# 2000steps = every 1ps
dcdfreq             2000
xstFreq             2000
outputEnergies       500
outputPressure       500

langevin             on ;# do langevin dynamics
langevinDamping       5 ;# damping coefficient (gamma) of 5/ps
langevinTemp          $desired_temp
langevinHydrogen    on ;# don't couple langevin bath to hydrogens

useGroupPressure    yes ;# needed for rigidBonds
useFlexibleCell       no
useConstantArea       no
langevinPiston       on
langevinPistonTarget 1.01325 ;#  in bar -> 1 atm
langevinPistonPeriod 200
langevinPistonDecay 500
langevinPistonTemp    $desired_temp
run                   20000000

保存为mdp吗？
3.coord坐标文件改为gro？
4.拓扑文件是将 all_wb_ion.psf  改为top  对嘛？

作者
Author: 胡姐姐 时间: 2018-1-10 19:14

胡姐姐发表于 2018-1-10 14:29
gmx grompp -f npt-nopr-md.mdp -c npt-pr.gro -p fws.top -o npt-nopr.tpr。这是教程中指出的生成tpr ...

@霜晨月您看见我的问题了么

作者
Author: 霜晨月 时间: 2018-1-10 19:27

胡姐姐发表于 2018-1-10 14:29
gmx grompp -f npt-nopr-md.mdp -c npt-pr.gro -p fws.top -o npt-nopr.tpr。这是教程中指出的生成tpr ...

你要是用gmx进行RMSF或PCA分析的话，只要pdb和trr（或xtc）就行了，不用生成tpr。但pdb和trr的原子必须一一对应。
PBC问题，指的是如果你用NAMD跑出来的dcd轨迹里面，蛋白跨越了PBC边界（VMD里面看上去，蛋白是断掉的、非常古怪的结构），你得先修正了这个问题，把蛋白变成完整的、正常的结构，然后才能用catdcd转换成trr。

作者
Author: fhh2626 时间: 2018-1-10 23:40

胡姐姐发表于 2018-1-10 13:36
那我是在算rmsf之前，先算一下整个轨迹的平均结构是吗？

最简单的办法就是载入平均结构--载入轨迹（在同一个molecule下）--点align--点rmsd--点plot

作者
Author: 霜晨月 时间: 2018-1-11 09:34

fhh2626 发表于 2018-1-10 23:40
最简单的办法就是载入平均结构--载入轨迹（在同一个molecule下）--点align--点rmsd--点plot

这样算出来的应该是各帧相对于平均结构的RMSD，即time series of RMSD using 平均结构 as the reference，不是楼主要的RMSF。

作者
Author: fhh2626 时间: 2018-1-11 14:12
本帖最后由 fhh2626 于 2018-1-11 14:33 编辑

霜晨月发表于 2018-1-11 09:34
这样算出来的应该是各帧相对于平均结构的RMSD，即time series of RMSD using 平均结构 as the reference ...

嗯，你说得对，还要再平均一次，见笑了
RMSF一般是算每个residue的振动幅度的，不是算整体的

欢迎光临计算化学公社 (http://ccc.keinsci.com/)