Analysis of AMBER molecular dynamics simulation results (conformational analysis)–HIV protease-inhibitor complex (6)
RMSD RMSF b-facto calculation
RMSD
RMSD measures the deviation of a target set of coordinates (i.e. a structure) to a reference set of coordinates, with
R
m
S
D.
=
0.0
\mathrm{RMSD}=0.0
RMSD=0.0 indicating a perfect overlap. RMSD is defined as:
R
m
S
D.
=
∑
i
=
0
N
[
m
i
?
(
x
i
?
Y
i
)
2
]
m
R M S D=\sqrt{\frac{\sum_{i=0}^N\left[m_i *\left(X_i-Y_i\right)^2\right]}{M}}
RMSD=M∑i=0N?[mi(XiYi?)2]?
B-Factoryl temperature factor | Debye-Waller factor
Function: The B factor reflects the “diffusion” of the electron density of atoms in the crystal. This “diffusion” actually reflects the conformational state of the protein molecule in the crystal. The higher the B factor, the greater the “diffusion”, and the corresponding The conformation of the site is more unstable. In crystallographic data, the B factor is generally given in atomic units, and we can convert it to the B factor of the corresponding residue to analyze the conformational stability of the residue.
RMSF
RMSF and B factor conversion formula:
R
m
S
f
2
=
3
B
8
π
2
R M S F^2=\frac{3 B}{8 \pi^2}
RMSF2=8π23B?
tleap.in
parm ../top/com.top trajin ../md/md1.crd 1 last 1 trajin ../md/md2.crd 1 last 1 reference top/com.pdb center :1-198 mass image center familiar rms reference @CA ##rmsd rms reference mass out rms_pro.dat :1-198@CA,C,N time 0.002 rms reference mass out rms_lig.dat :199 time 0.002 ##b-factor atomicfluct out bfactor_residue.dat :1-198@CA,C,N byres bfactor atomicfluct out bfactor_allatom.dat :1-198@CA,C,N byatom bfactor ## rmsf atomicfluct out rmsf_residue.dat :1-198@CA,C,N byres atomicfluct out rmsf_allatom.dat :1-198@CA,C,N byatom
implement
cpptraj -i tleap.in
visualization
xmgrace rms_pro.dat
Calculate rms2d
tleap.in
parm ../top/com.top trajin ../md/md2.crd rms2d mass :1-198@CA,C,N out 2drms_pro.gnu time 0.0002
The output image in gnuplot format can be opened with gnuplot
Distance angle calculation
tleap.in
parm ../top/com.top # trajin ../md/md1.crd trajin ../md/md2.crd ## MOL_199@O5 ILE_149@H ILE_149@N distance dist1 :149@N :199@O5 out dis.dat time 0.0002 ## MOL_199@O5 ILE_149@H ILE_149@N angle angle1 :199@O5 :149@H :149@N out angle.dat time 0.0002
implement
cpptraj -i tleap.in
SASA calculation
Accessible surface area (ASA, accessible surface area) or solvent-accessible surface area (SASA, solvent-accessible surface area) is the surface area of biomolecules accessible to solvents. Measurements of ASA are usually described in square angstroms (the standard unit of measurement in molecular biology). ASA was first described by Lee & Richards in 1971 and is sometimes referred to as the Lee-Richards molecular surface. ASA is usually calculated using the “rolling ball” algorithm developed by Shrake & Rupley in 1973. The algorithm uses a sphere (of the solvent) of a specific radius to “probe” the surface of the molecule.
The solvent-accessible surface is shown on the left as compared to the van der Waals surface. The van der Waals surface given by the atomic radius is shown in red. The accessible surface is drawn with dashed lines and was created by rolling the center of the sphere (blue) traced along the van der Waals surface of the sphere. Note that the scale bar for the probe radius depicted here is smaller than the typical 1.4A. The center of a solvent molecule as it rolls on the van der Waals surface of the protein, as shown in the right figure below (Gromiha and Ahmad, 2005). In general, water spheres are assumed to be solvent molecules with a radius of 1.4A. The solute molecules are represented by a set of interlocking spheres with appropriate van der Waals radii assigned to each atom, and the solvent molecules roll along the outer shell of the van der Waals surface in a convenient sectioning plane. Thus, the ASA of an atom of radius r is the region on the surface of a sphere of radius R = r + rsolv at which, at every point, the center of a solvent molecule can be brought into contact with that atom without penetrating any Other atomic solute molecules.
parm ../top/com.top # trajin ../md/md1.crd trajin ../md/md2.crd surf out surf.dat surf :1-10 out surf_10.dat surf :11-20 out surf_20.dat surf :1-20 out surf_all.dat
implement
cpptraj -i tleap.in
Hydrogen bond calculation
# Hydrogen bond analysis with cpptraj # Load topology and trajectory #trajin ../md1.crd #trajin ../md2.crd #trajin ../md3.crd #trajin ../md4.crd trajin ../md/md3.crd # hbond of protein and ligand hbond angle 120.000 dist 3.500 donormask: 1-306 acceptormask: 307 out nhb.dat avgout avghb.dat hbond angle 120.000 dist 3.500 donor mask: 224 acceptor mask: 1-223 out nhb.dat avgout avghb.dat #hbond of protein and bridging water #hbond angle 120.000 dist 3.500 donormask:1-223 acceptormask:232 out nhb.dat avgout avghb.dat #hbond angle 120.000 dist 3.500 donormask:232 acceptormask:1-223 out nhb.dat avgout avghb.dat #hbond of liqand and bridging water #hbond angle 120.000 dist 3.500 donormask:224 acceptormask:232 out nhb.dat avgout avghb.dat #hbond angle 120.000 dist 3.500 donormask:232 acceptormask:224 out nhb.dat avgout avghb.dat
cluster calculation
The word cluster can be translated into clusters or clustering. Although it is essentially the same in mathematics, it is still somewhat different when used specifically. Appropriate distinctions can be made clearer and avoid some ambiguities.
In the context of MD, we say cluster analysis, which generally means that, given a configuration, which contains multiple atoms, we assign these atoms to different clusters according to certain rules (usually distance) In this way, the system can be divided into some clusters, and each cluster can be analyzed to obtain some information.
The process of clustering, taking molecular docking as an example. After docking, each docking configuration is sorted by scoring from high to low.
- The first configuration is the first kind, and the number of the first kind of configuration is one.
- Compare the rmsd of the second configuration with the highest-scoring configuration in the first category. If rmsd is less than 2, the second configuration belongs to the first category, and the number of configurations in this category is increased by one; if rmsd is greater than 2, a new category is created, and the second configuration belongs to the newly created category, and the configuration of this category The quantity is one.
- Compare the rmsd of the nth configuration with the configuration with the highest score in each of the above classes in turn. If the rmsd of the mth class is less than 2, it belongs to the mth class. rmsd; if the rmsd with all classes is greater than 2, create a new class, the nth configuration belongs to the newly created class, and the number of configurations in this class is one.
tleap.in
parm ../top/com.top # trajin ../md/md1.crd trajin ../md/md2.crd cluster hieragglo epsilon 1.0 clusters 5 averagelinkage rms mass @CA,C,N out frame2cluster.dat summary summary.dat repout Cluster repfmt pdb
Based on quantum computing, semi-empirical, qm-mmpbsa
tleap.in
& general startframe=1, endframe=1000, interval=5, verbose=2, keep_files=2, / &gb igb=1, ifqnt=1, qmcharge_com=1, qm_residues='27,28,30,32,47,81,84,124,129,148,149,180,183,199', qm_theory='PM6-D', qmcharge_lig=0, qmcharge_rec=1, /
MMPBSA.py -O -i tleap.in -o qm_mmgbsa.dat -cp ../top/com.top -rp ../top/pro.top -lp ../top/lig.top -y ../md/md2.crd > qm_mmgbsa.log
Find residue id Find residues at 2.5 Angstroms
parm ../top/com.top trajin ../md/md2.crd 1 last 1 reference ../top/com.pdb center :1-198 mass image center familiar rms reference @CA ### rmsf atomicfluct out rmsf_residue.dat :199<@2.5 byres
visualization
xmgrace *.dat