Mercurial > repos > chemteam > bio3d_pca
diff rmsd.R @ 0:f56fb2762abb draft
planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/tree/master/tools/bio3d commit 580d80c3fd856bb3ae18ef7b5458eb3b5e2c7374
author | chemteam |
---|---|
date | Mon, 08 Oct 2018 12:49:36 -0400 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rmsd.R Mon Oct 08 12:49:36 2018 -0400 @@ -0,0 +1,67 @@ +#!/usr/bin/env Rscript + +options(stringAsfactors = FALSE) +args <- commandArgs(trailingOnly = TRUE) + +library(bio3d) + +dcdfile <- args[1] +pdbfile <- args[2] +selection <- args[3] + +dcd <- read.dcd(dcdfile) +pdb <- read.pdb(pdbfile) + + +if (selection == "string") { + domain <- args[4] + output <- args[5] + rmsd_plot <- args[6] + rmsd_hist <- args[7] + inds <- atom.select(pdb, string = domain) +} +if (selection == "resno") { + res1 <- args[4] + res2 <- args[5] + output <- args[6] + rmsd_plot <- args[7] + rmsd_hist <- args[8] + inds <- atom.select(pdb, resno=res1:res2) +} +if (selection == "elety") { + domain <- args[4] + output <- args[5] + rmsd_plot <- args[6] + rmsd_hist <- args[7] + inds <- atom.select(pdb, elety = domain) +} +if (selection == "resid") { + domain <- args[4] + output <- args[5] + rmsd_plot <- args[6] + rmsd_hist <- args[7] + inds <- atom.select(pdb, resid = domain) +} +if (selection == "segid") { + domain <- args[4] + output <- args[5] + rmsd_plot <- args[6] + rmsd_hist <- args[7] + inds <- atom.select(pdb, segid = domain) +} + +xyz <- fit.xyz(fixed=pdb$xyz, mobile=dcd, fixed.inds=inds$xyz, mobile.inds=inds$xyz) + +rd <- rmsd(xyz[1,inds$xyz], xyz[,inds$xyz]) + +write.table(rd, file = output, row.names = TRUE, col.names = FALSE, quote =FALSE, sep="\t") + +png(rmsd_plot) +plot(rd, typ="l", ylab="RMSD (Å)", xlab="Frame No.") +points(lowess(rd), typ="l", col="red", lty=2, lwd=2) +dev.off() + +png(rmsd_hist) +hist(rd, breaks=40, freq=FALSE, main="RMSD Histogram", xlab="RMSD") +lines(density(rd), typ="l", col="red", lty=2, lwd=2) +dev.off()