read.dcd error for MD trajectory

[Mamta Mohan]

Hello,

I am following tutorial : http://thegrantlab.org/bio3d/html/pca.tor.html

Example code run fine on trajectory provided as an example.

I decided to use code on my dcd file and I am having issues. The issue I am having is:

pc.tordcd <- pca.tor(tor.dcd[,gaps.dcd$f.inds]) Warning message: In as.xyz(xyz) : number of cartesian coordinates not a multiple of 3

To understand why I am having issue, I tried to back track the code and this is what I found:

dim(dcd) [1] 5000 4017

4017/3 = 1339 = NATOM

dcd <- read.dcd(dcdfile) NATOM = 1339 NFRAME= 5000 ISTART= 0 last = 5000 nstep = 5000 nfile = 5000 NSAVE = 1 NDEGF = 0 version 24 |================================================| 100% Warning message: In dcd.header(trj, verbose) : Check DCD header data is correct, particulary natom

for example trajectory these values are:

dim(trj) [1] 117 594 594/3 = 198 = NATOM trj <- read.dcd(system.file("examples/hivp.dcd", package="bio3d")) NATOM = 198 NFRAME= 117 ISTART= 0 last = 117 nstep = 117 nfile = 117 NSAVE = 1 NDEGF = 0 version 24 |==================================================| 100%

From the read.dcd it is apparent that dcd file header is not read correctly. Because protein is labelled with two labels, I think read.dcd skips both labelled residues and than gap.inspect, torsion.xyz and pca.tor all are fine in there own way but fail because of issue propagating from read.dcd.

Is there a way to address this issue so I can read dcd file correctly for analysis.

I would appreciate your response.

[Lars Skjærven]

How do you calculate tor.dcd? And how does gaps.dcd relate to tor.dcd?

[Mamta Mohan]

Dear Lars,

I followed example listed in the link: http://thegrantlab.org/bio3d/html/pca.tor.html for the error I posted before.

I tried it again today just to eliminate if I made any mistake.

read.dcd does put out a warning as I listed before. Number of atoms are 4017 and 4017/3 = 1339 Ignoring this warning is causing issue in pc.tor which show that it has 1336 atoms

pc.tor is getting calculated using gaps index applied to tor.dcd file. gaps.dcd also has 1336 in f.inds, which is the index of non-gap value.

This index is getting used with tor.dcd which has 1339, which is a mis-match.

So there are issues. One which is in read.dcd and another seems to be in gaps call.

I can provide more information if you want.

I am listing below the output of commands:

library("bio3d", lib.loc="~/R/x86_64-pc-linux-gnu-library/3.3") setwd("/home/mamta/SysLyzMtsl/NEUCluster/DataAnalysis/Bio3D") dcdfile <- "/home/mamta/SysLyzMtsl/NEUCluster/DataAnalysis/Bio3D/stripped.dcd" dcd <- read.dcd(dcdfile) NATOM = 1339 NFRAME= 5000 ISTART= 0 last = 5000 nstep = 5000 nfile = 5000 NSAVE = 1 NDEGF = 0 version 24 |=================================================| 100% Warning message: In dcd.header(trj, verbose) : Check DCD header data is correct, particulary natom tor.dcd <- t(apply(dcd, 1, torsion.xyz, atm.inc=1)) gaps <- gap.inspect(tor.dcd) pc.tor <- pca.tor(tor.dcd[,gaps$f.inds]) Warning message: In as.xyz(xyz) : number of cartesian coordinates not a multiple of 3

[Lars Skjærven]

Hi Mamta, This is just a warning and not an error. Your script seems fine, but you should take note of the atm.inc=1 argument. Also, make sure you have the desired atom selection. Perhaps only backbone atoms is what you want to include in the torsion.xyz() call? Please see the documentation for torsion.xyz: http://thegrantlab.org/bio3d/html/torsion.xyz.html L

[Mamta Mohan]

Thank you Lars. I will.