PCA and DCCM

Protocol PCA
1.Selection of the binding site for docking:
 a. Identify the binding site for search of the protein where the ligand is going to be docked.
 b. Define a search box around the selected site and run the docking software selected for the search.
2. Preparation of the protein-ligand complex:
 a. Extract the protein and ligand coordinates from the best docking pose file.
 b. Prepare the protein and ligand structures for molecular dynamics simulations by adding missing atoms, assigning protonation states,   and neutralizing the system with counterions if necessary.
 c. Assign force field parameters to the protein and ligand molecules.
 d. Add water molecules to solvate the protein-ligand complex and create a periodic simulation box. The size of the box should be large enough to avoid artificial interactions between periodic images of the system.
3. Preparation of the simulation system:
 a. Minimize the energy of the system using an energy minimization algorithm to remove any steric clashes or bad contacts.
 b. Equilibrate the system using a suitable equilibration protocol (e.g. NVT or NPT) to allow the water and ions to equilibrate around the protein-ligand complex.
3. Molecular dynamics simulations:
 a. Run several production molecular dynamics simulations of the protein-ligand complex in explicit solvent in NPT, with different initial velocities or random seed values.
 b. For each simulation, save the coordinates and velocities of the system at regular intervals for later analysis.
4. PCA analysis using Rstudio
Use the following commands in Rstudio. You require the bio3d library.
#### Commands for RStudio
# load library
library(bio3d)
#Trajectories
mydcdfile <- "m3.dcd"
#Reference pdb 
mypdbfile <- "m3.pdb"
#Load trajectories and reference
dcd <- read.dcd(mydcdfile)
pdb <- read.pdb(mypdbfile)
#Select CA atoms
ca.inds <- atom.select(pdb, elety="CA")
#Calculate PCA
xyz <- fit.xyz(fixed=pdb$xyz, mobile=dcd,fixed.inds=ca.inds$xyz,mobile.inds=ca.inds$xyz)
pc <- pca.xyz(xyz[,ca.inds$xyz])
#Save PCA
p1 <- mktrj.pca(pc, pc=1, b=pc$au[,1], file="pc1.pdb")
p2 <- mktrj.pca(pc, pc=2,b=pc$au[,2], file="pc2.pdb")
p3 <- mktrj.pca(pc, pc=3,b=pc$au[,3], file="pc3.pdb")
###################################################
Protocol DDCM
1. System Preparation.
 a. Follow steps 1-3 in previous section
2. DDCM Analysis
Use the cpptraj input below:

#CPPTRAJ SCRIPT
parm XX.top
trajin XX.crd 1 last 1
#select the residues of interest
strip !(:1-548) 
autoimage
rms first
#calculate DDCM on CA atoms
matrix correl :2-273,275-547@CA :2-273,275-547@CA out correl.gnu
go
quit

Plot Results using gnuplot
#GNUPLOT SCRIPT:

set terminal png
set border lw 3
set output 'XXX.png'
set view map
set title 'XXX' font ',18'
set ylabel 'Residue number' font ',16'
set xlabel 'Residue number' font ',16'
set cbrange [-1:1]
set xrange [1:546]
set yrange [1:546]
set palette defined ( -1 "blue", 0 "white", 1 "red" )
#in case of calculate the DCCM difference between two system, combine the data for the .gnu output to generate the difference file
plot 'correl.gnu' using 1:2:3 with points pointtype 6 palette linewidth 1 notitle