[R] "squashed" domains when doing PCA with bio3d

[Bianca Ihrig]

Hi everyone,

I have been doing PCA using bio3d. My protein of interest has two 
domains and there are two crystal structures, one in a closed (domains 
are close together) and one in an open state (domains are nearly 
opposite each other). I ran MD simulations on both states and did PCA. 
When doing PCA on either of the trajectories it is working fine but when 
I combine the two trajectories and all the frames to only one of the 
domains (as I would like to see the motion of the other domain for 
closed vs. open state).
But when doing so the second domain (the one I am not align to) is 
"squashed" when viewing the vectors in pymol. I dont understand why and 
I dont know how to fix it.
Many thanks in advance for your time and help!

This is the command I am using:

R
library(bio3d)

#read pdb file for open structure
pdbo <- read.pdb("ref-open(fr1800).pdb")
#to check: print(pdbo)
#read pdb file for closed structure
pdbc <- read.pdb("ref-closed(fr9990).pdb")

#read traj files
tc <- read.ncdf("A2-closed-1001fr-100ns.nc")
to <- read.ncdf("A2-open-100ns-1000fr.nc")

#combine the two traj (to and tc) into one called t2
t2 <- rbind(tc[2:1001,],to[2:1001,]) # 2000 frames 10413 coordinates

#Select residues of 1st domain (WHA) for alignment
whacac <- atom.select(pdbc, 
resno="9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,47,48,49,50,51,52,53,54,56,55,57,58,59,60,61,62,63,64,66,67,68,69,70,78,79,80,81,82", 
elety = "CA")
#Select CA of the two pdb files
cao <- atom.select(pdbo, elety = "CA") # 215 atoms
cac <- atom.select(pdbc, elety = "CA") # 215 atoms

# align the combined trajs to 1st domain of closed structure
xyz <- fit.xyz(fixed = pdbc$xyz, mobile = t2, fixed.inds = whacac$xyz, 
mobile.inds = whacac$xyz)
#to check use command: write.ncdf()

#####-- Do PCA over CA of the whole protein (not only domain 1)
pc <- pca.xyz(xyz[, cac$xyz])
png("pca_open-closed.png")
pc1 <- plot(pc)
dev.off()
#view pc vector in pymol
view.modes(pc, mode = 1, outprefix = "pc1", scale = 30, dual = FALSE, 
launch = TRUE)

Many thanks! Your help is much appreciated!