Torsion Angle vs. Time

[Michael Schwabe]

Hello, I am a huge fan of bio3d and have been using it as my predominant method for analyzing simulations. I am still new to scripting, but have managed to figure a lot of analysis steps by reading many of these issues. For example, the distance analysis provided by issue #446 has saved me a lot of time compared to similar analysis scripts used in VMD. I am currently struggling with the second part of issue #446 to determine the dihedral angle as a function of the simulation frame. This analysis can be done in VMD using the timeline feature, however, it is much slower and not as nice for analyzing the multitude of replicate simulations I have done.

In short, I would like to determine the phi and psi angles for a single residue over the course of the simulation. I have tried the solutions provided by #446, as well as #335, but both of these give me the angle for the first frame. I would like to write this output to a .csv file for use in a separate analysis. Here is what I currently have:

sele1 <- atom.select(pdb1, elety=c("C"), resno=69)
sele2 <- atom.select(pdb1, elety=c("N", "CA", "C"), resno=70)
inds.xyz = c(sele1$xyz, sele2$xyz)
phi <- torsion.xyz(pdb1$xyz[, inds.xyz])

sele3 <- atom.select(pdb1, elety=c("N", "CA", "C"), resno=70)
sele4 <- atom.select(pdb1, elety=c("N"), resno=71)
inds.xyz = c(sele3$xyz, sele4$xyz)
psi <- torsion.xyz(pdb1$xyz[, inds.xyz])

write.csv(phi, file = "R70_Phi.csv")
write.csv(psi, file = "R70_Psi.csv")

‌

[Xinqiu Yao]

Hi Michael,

‌

Great to hear you like Bio3D.

Regarding your question, it is most likely that the input pdb1 only contains one (the first) frame of coordinates. Can you print out the object to check? (For example, simply type pdb1 in the R session and press return. It will tell if your pdb is of single or multiple conformational models).

‌

‌

[Michael Schwabe]

Hi Xinqiu,

I have fit my starting structure (pdb1) to the simulation using the tutorial method:

setwd("/directory")
dcd1 <- read.dcd("200nsSimulation.dcd")
pdb1 <- read.pdb("ionized_protein.pdb")
ca.inds1 <- atom.select(pdb1, "protein")
xyz <- fit.xyz(fixed=pdb1$xyz, mobile=dcd1, fixed.inds=ca.inds1$xyz, mobile.inds=ca.inds1$xyz)

Do I need to write a new .pdb file with all the frames in the simulation, or can I run the calculation after the superimposition step?

[Xinqiu Yao]

Hi Michael,

‌

You don’t need to write a new .pdb file. When you call torsional calculations, use xyz rather than ‘pdb1’ to get the time series. For example,

phi <- torsion.xyz(xyz[, inds.xyz])

[Michael Schwabe]

Hi Xinqiu,

This was one of the first things I had tried, but I get the following error:

Warning message:
In dcd.header(trj, verbose) :
Check DCD header data is correct, particulary natom
Error in xyz[, inds.xyz] : invalid subscript type 'list'
Calls: torsion.xyz
Execution halted

Michael

[Xinqiu Yao]

It seems complaining the type of inds.xyz. That’s weird. Can you print the content of inds.xyz? It should be a plain numeric vector.

[Michael Schwabe]

Here is the script:

sele1 <- atom.select(pdb1, elety=c("C"), resno=60)
sele2 <- atom.select(pdb1, elety=c("N", "CA", "C"), resno=61)
inds.xyz <- c(sele1$xyz, sele2$xyz)
phi <- torsion.xyz(xyz[, inds.xyz])

sele3 <- atom.select(pdb1, elety=c("N", "CA", "C"), resno=61)
sele4 <- atom.select(pdb1, elety=c("N"), resno=62)
inds.xyz <- c(sele3$xyz, sele4$xyz)
psi <- torsion.xyz(xyz[, inds.xyz])

setwd("/directory")
write.csv(phi, file = "T61_Phi.csv")
write.csv(psi, file = "T61_Psi.csv")

I have tried changing sele$xyz to just sele as well and get the same error.

[Xinqiu Yao]

I recommend use combine.select() rather than concatenating. For example:

sele1 <- atom.select(pdb1, elety=c("C"), resno=60)
sele2 <- atom.select(pdb1, elety=c("N", "CA", "C"), resno=61)
inds <- combine.select(sele1, sele2, operator="+")
phi <- torsion.xyz(xyz[, inds$xyz])

‌

Note that in the last line, I use the dollar sign, inds$xyz, instead of period.

‌

Also, regarding your original question, I just realized torsion.xyz() is for a single snapshot, not a series, of structures.

You need to write a simple script to loop over simulation frames, each step calling torsion.xyz() for a single frame, and finally collect all angles:

phi <- apply(xyz[, inds$xyz], 1, torsion.xyz)

‌

[Michael Schwabe]

Thanks for all your help Xinqiu, this worked.