-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathTrajectory Analysis using Bio3D.txt
64 lines (37 loc) · 1.3 KB
/
Trajectory Analysis using Bio3D.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#Bio3D installation
> install.packages("bio3d", dependencies=TRUE)
> library(bio3d)
*#Trajectory Preparation
> pdb<-read.pdb("md.pdb")
> dcd<-read.dcd("md.dcd", cell=FALSE)
#Select CA atoms of the protein
> ca.inds<-atom.select(pdb,"protein",elety="CA")
#Fit the trajectory to the selection of protein atoms
> trj<-fit.xyz(pdb$xyz, dcd, fixed.inds=ca.inds$xyz, mobile.inds=ca.inds$xyz)
#Trim the PDB and DCD to include only the selected protein CA atoms
> protpdb<-trim.pdb(pdb,ca.inds)
> protdcd<-trim(trj, col.inds=ca.inds$xyz)
# Print the trimmed PDB coordinates
> print(protpdb$xyz)
#Print the trimmed trajectory coordinates
> print(protdcd)
##Trajectory Analysis
#RMSD
> rd<-rmsd(protlig_pdb, protlig_dcd, fit=TRUE, ncore=1) #change the "ncore" if you have multiple cores available
> plot(rd, typ='l',ylab="RMSD (Å)",xlab="Frame No.")
> points(lowess(rd), typ='l', col="red", lty=1, lwd=3)
> dev.copy(jpeg,filename="rmsd.jpg")
> dev.off()
#RMSF
> rf<-rmsf(protlig_dcd)
> plot(rf, ylab="RMSF (Å)",xlab="Residue")
> dev.copy(jpeg,filename="rmsf.jpg")
> dev.off()
#Principal component Analysis (PCA)
> pc<-pca.xyz(protlig_dcd)
> plot(pc,col=bwr.colors(nrow(dcd)))
> dev.copy(jpeg, filename="pca,jpg")
> dev.off()
#Dynamic Cross-Correlation Matrix (DCCM)
> dc<-dccm.xyz(protlig_dcd)
> plot.dccm(dc)