Help with loading of multiple replica simulations

[Mohit Mazumder]

Hi,

I am relatively new to the module. I have a very basic question !! I have two different (random seed) simulations data of 150 ns each on a single protein. How do i load multiple replicates (in my case two) for the system and do consensus correlation matrix and consensus network. I have tried the tutorial and read the suggested papers, still not able to do it.

[Xinqiu Yao]

Since you have only two replicates, I wouldn't suggest do consensus network analysis as we did in the paper, which usually requires 5 or more replicates. The simplest way is just to put trajectories together and calculate one correlation matrix. Calculate a contact map from the trajectories also. And then, you can use the cna() function to build network in a similar way to the original 'Dynamical Network Analysis' approach proposed by Luthey-Schulten's lab (PNAS 2009). For example, cna(cij, cutoff.cij=0.4, cm=cm), where cij is the correlation matrix and cm is the contact map.

Alternatively, you may think of splitting your trajectories into several simulation chunks (There are multiple ways to do that, e.g. VMD, Bio3D, AmberTools, etc.). Each chunk has smaller time scale (e.g. 50ns) but then you have multiple replicates and can use the consensus method. To do that, you need to:

1. Calculate cij for each chunk.
2. Calculate contact map from the cumulative trajectories.
3. Call filter.dccm() function to get the consensus matrix.
4. Call cna() function to build the network.

For example, suppose you have all chunks stored in 'trajs',

files <- list.files('trajs')

# this is your topology pdb file.
pdb <- read.pdb('top.pdb') 

chunks <- lapply(files, read.ncdf, at.sel=atom.select(pdb, 'calpha'))
traj <- do.call(rbind, chunks)
cijs <- lapply(chunks, dccm)
cm <- cmap.xyz(traj, dcut=10, scut=1, pcut=0.75)

# calculate consensus cij
cij <- filter.dccm(cijs, cutoff.cij=0.4, cmap=cm)

# build the network (Note the cutoff of cij here is 0 because we have already removed weak couplings from above calling of the **filter.dccm()**
net <- cna(cij, cutoff.cij=0)

Let me know if you have more questions.

[Mohit Mazumder]

Dear Mr Yao,

Thank you so much for the detailed help !!! Much appreciated.

I have been able to load the file but got this error while calculating consensus cij

cij <- filter.dccm(cijs, cutoff.cij=0.4, cmap=cm)

Error in if (cmap) { : argument is not interpretable as logical In addition: Warning message: In if (cmap) { : the condition has length > 1 and only the first element will be used

[Xinqiu Yao]

Hi,

You need to install the latest DEVELOPMENT version bio3d to use filter.dccm() with externally calculated cmap. See here for how to install.

For the released version, proper way is filter.dccm(cijs, cutoff.cij=0.4, cmap=TRUE, xyz=traj, dcut=10, scut=1, pcut=0.75), which will calculate cmap internally.

[Mohit Mazumder]

Hi,

I tried and installed the required version . Now it show :

cm <- cmap.xyz(traj, dcut=10, scut=1, pcut=0.75) Error in if (!(dims[2L]%%3 == 0)) warning("number of cartesian coordinates not a multiple of 3") : argument is of length zero

[Xinqiu Yao]

Check by dim(traj). Is it correct? If not, check the files. Does it contain file names you want? (You may need to slightly modify the command shown above, i.e. use files <- list.files('trajs', full.names=TRUE) and try again).

[Barry Grant]

We need a posting guide that tells folks posting issues to explicitly include enough detail so we can help them - give session.info() specs, detail steps and code to reproduce, details of actual results and expected results if appropriate.

Remember, no one can help if we aren't clear what the issue is and how it arose.

[Xinqiu Yao]

Cool idea! We can put the message to the template that will show up whenever a new issue is created (I think there is a way to do it in the settings). Will try soon.

[Xinqiu Yao]

Have added a temp user guide for issue creation (Try to create a new issue to see the effect). You can also modify the message via the Settings.

[Barry Grant]

Thanks Xinqiu! I have added a little more details to it. Hopefully this will make things easer in the future.