Torsion PCA issue

[Karan Kapoor]

I am trying to do PCA of a trajectory in the torsion (phi-psi) space and have some related question.

First of all I am checking if I can calculate these angles for my starting PDB file. But instead of getting 194 angles for the 194 selected residues, I am only getting 145 angles. Why is this happening? I have attached the starting PDB for your reference.

pdbfile1 = "protein.pdb"
pdb1 = read.pdb(pdbfile1)

sel1 = atom.select(pdb1, "protein", elety=c("C"), resno=1)
sel2 = atom.select(pdb1, "protein", elety=c("N", "CA", "C"), resno=2:195)
selphi <- combine.select(sel1, sel2, operator="+")

sel3 = atom.select(pdb1, "protein", elety=c("N", "CA", "C"), resno=2:195)
sel4 = atom.select(pdb1, "protein", elety=c("N"), resno=196)
selpsi <- combine.select(sel3, sel4, operator="+")

phipdb <- torsion.xyz(pdb1$xyz[, selphi$xyz])
psipdb <- torsion.xyz(pdb1$xyz[, selpsi$xyz])

I am then using the following to calculate these angles for my trajectory containing 1000 frames.

phidcd <- apply(xyz[, selphi$xyz], 1, torsion.xyz)
psidcd <- apply(xyz[, selpsi$xyz], 1, torsion.xyz)

I get a vector of 145x1000 dimensions for each of the above. How can I now combine these matrices to do the PCA with pca.tor? I am trying the following but not sure if it is the correct way.

#Taking matrix transpose
phidcd_t = t(phidcd)
psidcd_t = t(psidcd)

#Combining the columns
tor = cbind(phidcd_t, psidcd_t)
pc.tor <- pca.tor(tor)

Would it be possible to get the residue-wise contributions (loadings) for this PCA similar to the cartesian PCA (pc.car$au)? Right now, plot.pca.loadings(pc.tor) gives index-wise contributions which is not what I want and pc.tor$au gives NULL.

Additionally, can I also output the trajectories representing each PC? I get the following error when using mktrj.

mktrj(pc.tor, pc=1, file="pc1_CA.pdb") 

Error in write.pdb(xyz = coor, file = file, pdb = pdb, ...) :
write.pdb: 'length(xyz)' must be divisable by 3.
In addition: Warning message:
In as.xyz(t(coor)) : number of cartesian coordinates not a multiple of 3

Thank you for your time!

[Xinqiu Yao]

Hi,

‌

To calculate a single PDB, torsion.pdb() is recommended. For trajectory, we have a version under “new_funs/” in the bio3d repository, called xyz2torsion.R, which is not included in the main package yet but should work well.

‌

torsion.xyz() by default increases atom numbers by 4. You can change it to a smaller number, but be careful with the result as the atom order may not be the correct one. torsion.pdb() and xyz2torsion() would be better choices.

‌

For torsion PCA, you must be careful about periodicity. Check out this paper for a more detailed discussion: https://aip.scitation.org/doi/10.1063/1.4998259