RMSD of two different trajectories in one graph

[Former user]

Hello,

I would like to graph two RMSDs on the same plot, but I have had an issue with this task. I have two different pdbs and two different dcds, but it only plots data from only one data set. Here's the two scripts that I have tried using.

library(bio3d) dcd <- read.dcd("outlook1all.dcd") pdb <- read.pdb("SRP4354all.pdb")

ca.inds <- atom.select(pdb, elety="CA") xyz <- fit.xyz(fixed=pdb$xyz, mobile=dcd, fixed.inds=ca.inds$xyz, mobile.inds=ca.inds$xyz) dim(xyz) == dim(dcd)

rd <- rmsd(xyz[1,ca.inds$xyz], xyz[,ca.inds$xyz])

plot(rd, typ="l", ylab="RMSD", xlab="Frame No.") #points(lowess(rd), typ="l", col="red", lty=2, lwd=2)

dcd2 <- read.dcd("outlookall.dcd") pdb2 <- read.pdb("SRP43.pdb")

ca2.inds <- atom.select(pdb2, elety="CA")

xyz2 <- fit.xyz(fixed=pdb$xyz, mobile=dcd2, fixed.inds=ca2.inds$xyz, mobile.inds=ca2.inds$xyz)

dim(xyz2) == dim(dcd2)

rd2 <- rmsd(xyz[1,ca2.inds$xyz], xyz[,ca2.inds$xyz])

points(rd2, typ="l", ylab="RMSD", xlab="Frame No.", col ="red") #points(lowess(rd2), typ="l", col="red", lty=2, lwd=2)

Here's the other one

library(bio3d) setwd("path") library(bio3d) files <- list.files("path", ".\.pdb$", full.names=TRUE) pdbs <- pdbaln("SRP4354all.pdb, SRP43.pdb") Note: Accessing online PDB files using 4 letter PDBID Error in pdbaln("SRP4354all.pdb, SRP43.pdb") : * Missing files: check filenames SRP4354all.pdb, SRP43.pdb pdbs <- pdbaln("SRP4354all.pdb", "SRP43.pdb") Reading PDB files: SRP4354all.pdb .

Extracting sequences

pdb/seq: 1 name: SRP4354all.pdb Error in if (fit) s$xyz <- pdbfit(s) : argument is not interpretable as logical In addition: Warning message: In seqaln(s, id = files, ...) : nothing to align

pdbs <- pdbaln("SRP4354all.pdb", "SRP43.pdb")

[Lars Skjærven]

So what's the trouble with the first section here?

Input to pdbaln should be vector of PDB file names. Try pdbaln(c("SRP4354all.pdb", "SRP43.pdb")).

[Mitchell Benton]

The first section kept loading the first trajectories RMSD into the graph instead of the second one, so I ended up graphing two copies of the same plot on the graph.

The pdbaln command line you gave me did work, but what would be the script for loading the two dcds so I can get the RMSDs?

[Lars Skjærven]

Hi Mitchell, What exactly do you want to calculate the RMSDs of? Your first script looks like a good start to calculate the rmsd along the two simulations, with reference to the starting structure (assuming the two pdb files contain the starting structures of the simulation). Or do you want to fit the two trajectories to each other and calculate the rmsd for a common atom subset?

[Mitchell Benton]

I would like to take each RMSD from the separate trajectory files and plot them on the same graph, so I can have a "side-by-side" analysis of the two RMSDs. I have both RMSDs by themselves, but not in a combined plot. The two different pdbs should be the appropriate starting structure.

[Lars Skjærven]

# rmsd of sim 1
dcd1 <- read.dcd("outlook1all.dcd") 
pdb1 <- read.pdb("SRP4354all.pdb")
ca.inds1 <- atom.select(pdb1, "calpha") 
xyz1 <- fit.xyz(fixed=pdb1$xyz, mobile=dcd1, fixed.inds=ca.inds1$xyz, mobile.inds=ca.inds1$xyz) 
rd1 <- rmsd(xyz1[1, ca.inds1$xyz], xyz1[, ca.inds1$xyz])

# rmsd of sim 2
dcd2 <- read.dcd("outlookall.dcd") 
pdb2 <- read.pdb("SRP43.pdb")
ca.inds2 <- atom.select(pdb2, "calpha") 
xyz2 <- fit.xyz(fixed=pdb2$xyz, mobile=dcd2, fixed.inds=ca.inds2$xyz, mobile.inds=ca.inds2$xyz) 
rd2 <- rmsd(xyz2[1, ca.inds2$xyz], xyz2[, ca.inds2$xyz])

ylim = c(0, max(c(rd1, rd2)))

plot(rd1, typ="l", ylab="RMSD", xlab="Frame No.", ylim=ylim)
lines(rd2, col=2)

Hope it helps. Note that this is not a bio3d specific issue, so you should consider to read up on plotting in R.