Principal Component Analysis of the Trp-Cage MD Simulation

library(bio3d)
traj <- read.ncdf("simulation_325/trajectory.netcdf")

## [1] "Reading file simulation_325/trajectory.netcdf"
## [1] "Produced by program: cpptraj"
## [1] "File conventions AMBER version 1.0"
## [1] "Frames: 250000"
## [1] "Atoms: 304"

reference_pdb <- read.pdb("../1L2Y_Conf_01.pdb") # Experimental structure

##   HEADER    DE NOVO PROTEIN                         25-FEB-02   1L2Y

traj_pdb <- read.pdb("traj.pdb") # Frame from trajectory, required because atom indices differ from the experimental

print(reference_pdb$xyz)

## 
##    Total Frames#: 1
##    Total XYZs#:   912,  (Atoms#:  304)
## 
##     [1]  -8.901  4.127  -0.555  <...>  -1.341  10.903  1.473  [912] 
## 
## + attr: Matrix DIM = 1 x 912

print(traj)

## 
##    Total Frames#: 250000
##    Total XYZs#:   912,  (Atoms#:  304)
## 
##     [1]  2.378  0.638  0.629  <...>  -93.672  27.127  -140.77  [228000000] 
## 
## + attr: Matrix DIM = 250000 x 912

fixed.inds <- atom.select(reference_pdb, "backbone", elety="CA")
mobile.inds <- atom.select(traj_pdb, "backbone", elety="CA")

xyz <- fit.xyz(
  fixed=reference_pdb$xyz, mobile=traj,
  fixed.inds=fixed.inds$xyz,
  mobile.inds=mobile.inds$xyz
)
pc <- pca.xyz(xyz)
plot(pc, col=bwr.colors(nrow(xyz)))

p1 <- mktrj.pca(pc, pc=1, b=pc$au[,1], file="pc1.pdb", resid=traj_pdb$atom$resid, resno=traj_pdb$atom$resno, eleno=traj_pdb$atom$eleno, elety=traj_pdb$atom$elety)
p2 <- mktrj.pca(pc, pc=2,b=pc$au[,2], file="pc2.pdb", resid=traj_pdb$atom$resid, resno=traj_pdb$atom$resno, eleno=traj_pdb$atom$eleno, elety=traj_pdb$atom$elety)
p3 <- mktrj.pca(pc, pc=3,b=pc$au[,2], file="pc3.pdb", resid=traj_pdb$atom$resid, resno=traj_pdb$atom$resno, eleno=traj_pdb$atom$eleno, elety=traj_pdb$atom$elety)
p4 <- mktrj.pca(pc, pc=4,b=pc$au[,2], file="pc4.pdb", resid=traj_pdb$atom$resid, resno=traj_pdb$atom$resno, eleno=traj_pdb$atom$eleno, elety=traj_pdb$atom$elety)