cij <- dccm(trj) failed

[B Albert]

Hello

I am following the tutorial for network analysis in page:

http://thegrantlab.org/bio3d/tutorials/protein-structure-networks

with my own system. Here are my comamnds:

library(bio3d)
library(igraph)

pdb <- read.pdb("npt3d.pdb")
dcd <-read.dcd("top.dcd")

inds <- atom.select(pdb, resno = c(61:67, 153:165), elety = "CA")
trj <- fit.xyz(fixed = pdb$xyz, mobile = dcd,
               fixed.inds = inds$xyz, mobile.inds = inds$xyz)
cij <- dccm(trj)

but the job failed with following messages:

> inds <- atom.select(pdb, resno = c(61:67, 153:165), elety = "CA")
> trj <- fit.xyz(fixed = pdb$xyz, mobile = dcd,
+                fixed.inds = inds$xyz, mobile.inds = inds$xyz)
> cij <- dccm(trj)

 *** caught segfault ***
address 0x2ac299f669f8, cause 'memory not mapped'

Traceback:
 1: cov(dxyz)
 2: dccm.xyz(trj)
 3: dccm(trj)

Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
Selection: 

[Barry Grant]

Thanks for reporting this bug. It is one other R other users have reported with a number of other packages and indeed base functions e.g.: < http://stackoverflow.com/questions/25672974/caught-segfault-error-in-r >

The advise in a number of these posts is to uninstall all installed packages and start over. I am also wondering if you are running out of memory (the code should however not segfault but rather return you to the R prompt with a warning). Perhaps we could implement a catch of this with a solution like this one: < http://www.r-bloggers.com/bigcor-large-correlation-matrices-in-r/ >

How big is 'trj' object (i.e. type 'dim(trj)' or just 'trj' to report the number of rows and cols please).

[Xinqiu Yao]

Yes, it is possibly the 'trj' object is too big. Note that after calling fit.xyz() the returned value is the full atom 'xyz' not just for CA. Does your DCD file contain all heavy atoms? Are you going to build a network with CA atoms only? If both answers are yes, you might want to try

cij <- dccm(trj[,inds$xyz])

or

spdb <- trim(pdb, resno=pdb$atom[inds$atom, "resno"])
cij <- dccm(trj, grpby = spdb$atom[, "resno"])

The latter is to calculate correlations for all atoms, take maximum value for each residue pair, and return a reduced matrix to represent effective CA-CA couplings.

Let me know if this is relevant to your case.