###
### Example of basic molecular dynamics trajectory analysis
###
### Authors Xin-Qiu Yao
### Lars Skjaerven
### Barry J Grant
###
require(bio3d); require(graphics);
pause <- function() {
cat("Press ENTER/RETURN/NEWLINE to continue.")
readLines(n=1)
invisible()
}
#############################################
## #
## Basic analysis of HIVpr trajectory data #
## #
#############################################
pause()
Press ENTER/RETURN/NEWLINE to continue.
# Read example trajectory file
trtfile <- system.file("examples/hivp.dcd", package="bio3d")
trj <- read.dcd(trtfile)
NATOM = 198
NFRAME= 117
ISTART= 0
last = 117
nstep = 117
nfile = 117
NSAVE = 1
NDEGF = 0
version 24
Reading (x100)
|======================================================================| 100%
# Read the starting PDB file to determine atom correspondence
pdbfile <- system.file("examples/hivp.pdb", package="bio3d")
pdb <- read.pdb(pdbfile)
# Whats in the new pdb object
print(pdb)
Call: read.pdb(file = pdbfile)
Total Models#: 1
Total Atoms#: 198, XYZs#: 594 Chains#: 2 (values: A B)
Protein Atoms#: 198 (residues/Calpha atoms#: 198)
Nucleic acid Atoms#: 0 (residues/phosphate atoms#: 0)
Non-protein/nucleic Atoms#: 0 (residues: 0)
Non-protein/nucleic resid values: [ none ]
Protein sequence:
PQITLWQRPLVTIKIGGQLKEALLDTGADDTVLEEMSLPGRWKPKMIGGIGGFIKVRQYD
QILIEICGHKAIGTVLVGPTPVNIIGRNLLTQIGCTLNFPQITLWQRPLVTIKIGGQLKE
ALLDTGADDTVLEEMSLPGRWKPKMIGGIGGFIKVRQYDQILIEICGHKAIGTVLVGPTP
VNIIGRNLLTQIGCTLNF
+ attr: atom, xyz, helix, sheet, calpha,
call
pause()
Press ENTER/RETURN/NEWLINE to continue.
# How many rows (frames) and columns (coords) present in trj
dim(trj)
[1] 117 594
ncol(trj) == length(pdb$xyz)
[1] TRUE
pause()
Press ENTER/RETURN/NEWLINE to continue.
# Trajectory Frame Superposition on Calpha atoms
ca.inds <- atom.select(pdb, elety = "CA")
xyz <- fit.xyz(fixed = pdb$xyz, mobile = trj,
fixed.inds = ca.inds$xyz,
mobile.inds = ca.inds$xyz)
# Root Mean Square Deviation (RMSD)
rd <- rmsd(xyz[1, ca.inds$xyz], xyz[, ca.inds$xyz])
plot(rd, typ = "l", ylab = "RMSD", xlab = "Frame No.")
points(lowess(rd), typ = "l", col = "red", lty = 2, lwd = 2)
summary(rd)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 1.834 1.967 2.151 2.167 4.112
pause()
Press ENTER/RETURN/NEWLINE to continue.
# Root Mean Squared Fluctuations (RMSF)
rf <- rmsf(xyz[, ca.inds$xyz])
plot(rf, ylab = "RMSF", xlab = "Residue Position", typ="l")
pause()
Press ENTER/RETURN/NEWLINE to continue.
# Principal Component Analysis
pc <- pca.xyz(xyz[, ca.inds$xyz])
plot(pc, col = bwr.colors(nrow(xyz)))
pause()
Press ENTER/RETURN/NEWLINE to continue.
# Cluster in PC space
hc <- hclust(dist(pc$z[, 1:2]))
grps <- cutree(hc, k = 2)
plot(pc, col = grps)
pause()
Press ENTER/RETURN/NEWLINE to continue.
# Cross-Correlation Analysis
cij <- dccm(xyz[, ca.inds$xyz])
plot(cij)
## view.dccm(cij, pdb, launch = TRUE)