plotting issue related with cmap library

[Former user]

I have used the tutorial to get contact map over trajectory. Here is the R script:

library(bio3d)
pdb <- read.pdb("step1_pdbreader.pdb")
inds <- atom.select(pdb, "calpha")
ref.cont <- cmap( pdb$xyz[inds$xyz], dcut=8, scut=0, mask.lower = FALSE )
plot.cmap(ref.cont)
dev.copy(png,'ref_contact_map.png')
dev.off()

trj <- read.dcd("tenframes.dcd")
sum.cont <- NULL
for(i in 1:nrow(trj)) {

        #contact map for frame i
        cont <- cmap(trj[i,inds$xyz], dcut=8, scut=0, mask.lower = FALSE)
        #product with reference
        prod.cont <- ref.cont * cont
        sum.cont <- c(sum.cont, sum(prod.cont, na.rm=TRUE))
}

plot(sum.cont, type="l")
dev.copy(png,'sum.cont_over_trajectories_trj.png')
dev.off()

the structure I am calculating contact map for has four helices composed of two different type of chains. Meaning that structure is a tetramer like ABAB. I want to plot cmap in a way that as opposed to residue index, i want to have residue numbers in x and y axes and i want to have only A and B chains (so i want to combine the contact information of two identical A chains and two identical B chains). But i don't know how.. Thanks in advance

[Barry Grant]

How about supplying a RCSB PDB id so we can reproduce your example code? Thanks!

[Former user]

It is a modified pdb file. I attached both pdb and dcd files. Thanks a lot

[Xinqiu Yao]

Hi,

The plot.cmap() function doesn't support 'resno' option yet. But you can customized axis tick labels by following scripts:

plot.cmap(ref.cont, axes=FALSE)
box()
at <- axTicks(1); at[1] = 1
labels <- pdb$atom[pdb$calpha, 'resno'][at]
axis(1, at=at, labels=labels)
axis(2, at=at, labels=labels)

To plot for only two chains or any part of pdb, use atom.select() to get subset of atoms. For example,

inds <- atom.select(pdb, "calpha", chain=c('A', 'B'))
ref.cont <- cmap( pdb$xyz[inds$xyz], dcut=8, scut=0, mask.lower = FALSE )
plot.cmap(ref.cont)

However, it seems that your pdb has a unique chain ID 'P' for all residues. So, you may have to fix this first. Try to find the boundary of chains and label each part with distinct chain ID (e.g. 'A', 'B', 'C', or 'D'). This can be easily done with many editors such as vim.

Let me know if this solve your problem or not.

[Former user]

instead of chain names, i have used segname as PROA and PROB. It worked! But my question was how to plot PROA and PROB separately in the same graph. so x and y axes will be like from residue number 133 ... 341 133 ... 341
Huge thanks for the x-y labeling!

[Xinqiu Yao]

What about this:

pdbA <- atom.select(pdb, 'calpha', segid='PROA', value=TRUE)
pdbB <- atom.select(pdb, 'calpha', segid='PROB', value=TRUE)
ref.contA <- cmap( pdbA$xyz, dcut=8, scut=0, mask.lower = TRUE )
ref.contB <- cmap( pdbB$xyz, dcut=8, scut=0, mask.lower = FALSE )
ref.contB[upper.tri(ref.contB)] <- NA

plot.cmap(ref.contA, axes=FALSE)
box()
at <- axTicks(1); at[1] = 1
labels <- pdbA$atom[pdbA$calpha, 'resno'][at]
axis(1, at=at, labels=labels)
at <- axTicks(2); at[1] = 1
labels <- pdbB$atom[pdbB$calpha, 'resno'][at]
axis(2, at=at, labels=labels)

plot.cmap(ref.contB, add=TRUE, col='blue')

[Former user]

that helps me see the contacts in two different chains in the same plot. i guess i was trying to do the same by continuing the data for chain B at the end of chain A's. Thanks a lot Xin-Qiu.