issue with comparing MD PCs and X-ray PCs

[Former user]

Following instructions in NMA DHFR part 1 tutorial, I repeatedly get stuck at the following point with fitting for my system. I was able to replicate the DHFR tutorial without any issues. What could be the reason for this error? I cannot figure it out.

trj <- trj[, atom2xyz(md.inds)] trj <- fit.xyz(fixed=pdbs$xyz[1, ], mobile=trj, + fixed.inds=core$c0.5A.xyz, mobile.inds=core$c0.5A.xyz) Error in mobile[, mobile.inds] : subscript out of bounds

Also in md.inds <- pdb2aln.ind(aln=pdbs, pdb=pdb, id="md", inds=gaps.res$f.inds) how can I specify alignment to specific structure in the pdbs?

Is this the correct use?

md.inds <- pdb2aln.ind(aln=pdbs, pdb=pdb, id="md", inds=gaps.res$f.inds,aln.id=41) where 41 is the structure with same sequence as MD simulation pdb structure.

Thank you.

[Xinqiu Yao]

Hi,

About pdb2aln.ind, to specify a particular structure in the pdbs, use the ID of that structure NOT the number. In this case, proper command should be

md.inds <- pdb2aln.ind(aln=pdbs, pdb=pdb, id="md", inds=gaps.res$f.inds, aln.id=pdbs$id[41])

I am not sure about the reasons causing the errors due to the lack of MD data. Maybe Lars can help with this?

[Lars Skjærven]

Hey, I double checked the DHFR example. It works on my side, but perhaps you are missing the trajectory file and the PDB? Did you evaluate the following lines in the tutorial:

pdb <- read.pdb("md-traj/1rx2-CA.pdb")
trj <- read.ncdf("md-traj/1rx2_5ns.nc")

You probably need to download those from the repository: pdb, traj

Note that the code in the DHFR example might not work in the for another protein system. Below I draft an alternative approach which shoud work better.

Xin-qiu: I'm not sure I understood the correct usage of pdb2aln.ind. I found it difficult to use, and thus drafted an alternative, which is basically a wrapper for your function. See usage below. What do you think?

pdbsinds2 <- function(pdbs, pdb, inds=NULL) {
  if(is.null(inds))
    f.inds <- gap.inspect(pdbs$resno)$f.inds
  else {
    f.inds <- inds
  }

  inds.b <- pdb2aln.ind(aln=pdbs, pdb=pdb, id="new", inds=f.inds)
  inds.a <- inds.b

  inds.a <- f.inds[!is.na(inds.b)]
  inds.b <- inds.b[!is.na(inds.b)]

  out <- list(a=list(atom=inds.a, xyz=atom2xyz(inds.a)),
              b=list(atom=inds.b, xyz=atom2xyz(inds.b)))
  return(out)
}

pdb <- read.pdb("md-traj/1rx2-CA.pdb")
trj <- read.ncdf("md-traj/1rx2_5ns.nc")

# find indices to new pdb to core
inds <- pdbsinds2(pdbs, pdb, core$atom)
trj <- fit.xyz(pdbs$xyz[1,], trj, inds$a$xyz, inds$b$xyz)

# find indices for projection
inds <- pdbsinds2(pdbs, pdb)
proj <- project.pca(trj[,inds$b$xyz], pc.xray)
cols <- densCols(proj[,1:2])

plot(proj[,1:2], col=cols, pch=16,
     ylab="Prinipcal Component 2", xlab="Principal Component 1",
     xlim=range(pc.xray$z[,1]), ylim=range(pc.xray$z[,2]))
points(pc.xray$z[,1:2], col=1, pch=1, cex=1.1)
points(pc.xray$z[,1:2], col=grps.rd, pch=16)

[Xinqiu Yao]

Hi Lars,

Good idea! So, the return value contains both input and output indices, is it correct? We can integrate it with the original pdb2aln.ind function. Thanks for the update!

[Xinqiu Yao]

Hi,

I've updated pdb2aln.ind() (see commit) and now it can be directly used the same way as in above:

inds <- pdb2aln.ind(pdbs, pdb, core$atom)
trj <- fit.xyz(pdbs$xyz[1,], trj, inds$a$xyz, inds$b$xyz)

[Lars Skjærven]

Seems to work!