How can I filter residue cross correlation plots in bio3d?

[Former user]

I am using bio3d package to generate residue cross correlation plots. I am interested in cross correlations that are not caused by the residues touching each other. Can anyone tell me ‘filter’ these plots by blanking out :

1. all identities (the diagonals)

2. those residues that touch each other

[Barry Grant]

How about calculating a contact map, with whatever distance criteria you want, and then invert it (make the 1s into 0s and vice-versa). Then you can multiply your correlation plot by your inverted contact map to effectively mask out (set to zero) these values.

Or more simply using the 1s in your original contact map to set the elements of your correlation matrix to zero.

[Xinqiu Yao]

Hi,

I think Barry's method is the easiest one. Here I give some more details:

# Suppose your cross correlation matrix is named 'cij', 
# MD trajectory 'xyz' and MD topology 'pdb'.

# Calculate contact map with the MD trajectory
noh <- atom.select(pdb, 'noh') # this is to exclude hydrogen atoms
cm <-   cmap(xyz[, noh$xyz], grpby=pdb$atom[noh$atom, 'resno'], scut=0, pcut=0.7, mask.lower=FALSE)

# Alternatively, calculate contact map with a PDB
#cm <- cmap(pdb, scut=0, pcut=1, mask.lower=FALSE)

# Note that you can always try other parameters in above function. Type 'help(cmap)' for more details.
diag(cm) <- 1

# Filter cij with inverse of cmap
new.cij <- cij
new.cij[cm==1] <- 0

# Plot cij
plot(new.cij, sse=pdb)
# Customize your plot by checking 'help(plot.dccm)'

Let me know if you have more questions!

[Khurram Mushtaq]

Thank you very much Barry Grant and Xin-Qiu Yao!

I am actually new to this package. So, I generated residue cross correlation plots as described on [http://thegrantlab.org/bio3d/tutorials/normal-mode-analysis]. What I understand that it calculated cross correlation based on modes of your protein not the MD trajectory. So, can you please tell me how to preform this analysis based on modes.

Also, I have MD trajectories but that is in .xtc/.trr (gromacs) format. Is there any way to read those files into 'xyz'. Sorry, for the naive questions, I am actually new to this.

[Xinqiu Yao]

Hi Khurram,

You can calculate contact map directly from pdb also. For example, cm <- cmap(pdb). Again, check help(cmap) to learn the meaning of arguments and adjust them if necessary.

bio3d doesn't read .xtc/.trr format currently. Convert them to DCD (Charmm format) or NetCDF (Amber format) files first. There are plenty of software to do the conversion I believe. For example, read you trajectories with VMD in one format and write them out in another (check VMD manual for more details). Also, you may try openbabel, which is more straightforward.

[Khurram Mushtaq]

Dear Xin-Qiu Yao,

Thank you very much for your inputs. I gave the following code:

pdb <- read.pdb("ct.pdb")

modes <- nma(pdb)

cij <- dccm(modes)

cm <- cmap(pdb)

diag(cm) <- 1

new.cij <- cij

new.cij[cm==1] <- 0

plot(new.cij, sse=pdb)

This is what my output look like. I hope I am successful in blanking the residues that touch each other, what do you say? Do you any other comments regarding it.

[Xinqiu Yao]

Hi Khurram,

This looks good. A couple things to mention:

First, you may still have couplings shown above for residues that are close to each other in sequence. Note that by default 'cmap()' doesn't calculate contact map for neighboring residues. To mask off these residue pairs, you may replace the corresponding command to:

cm <- cmap(pdb, scut=0, mask.lower=FALSE)

Note that even if you do not want to exclude residue pairs that are neighbors in sequence, the mask.lower=FALSE is still needed in the call of 'cmap()'. Otherwise, your 'new.cij' will be asymmetric.

Second, it seems your pdb doesn't contain secondary structure information. Try sse = dssp(pdb) and then plot(new.cij, sse=sse). This will add SSE annotation on the top and right to facilitate your analysis.

Good luck!

[Khurram Mushtaq]

Dear Xin-Qiu Yao,

Thank you very much! Really appreciate your help. It made my analysis easy and interesting.