Mercurial > repos > chemteam > bio3d_pca

comparison rmsd.R @ 0:f56fb2762abb draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
-1:000000000000 0:f56fb2762abb
1 #!/usr/bin/env Rscript
2
3 options(stringAsfactors = FALSE)
4 args <- commandArgs(trailingOnly = TRUE)
5
6 library(bio3d)
7
8 dcdfile <- args[1]
9 pdbfile <- args[2]
10 selection <- args[3]
11
12 dcd <- read.dcd(dcdfile)
13 pdb <- read.pdb(pdbfile)
14
15
16 if (selection == "string") {
17 domain <- args[4]
18 output <- args[5]
19 rmsd_plot <- args[6]
20 rmsd_hist <- args[7]
21 inds <- atom.select(pdb, string = domain)
22 }
23 if (selection == "resno") {
24 res1 <- args[4]
25 res2 <- args[5]
26 output <- args[6]
27 rmsd_plot <- args[7]
28 rmsd_hist <- args[8]
29 inds <- atom.select(pdb, resno=res1:res2)
30 }
31 if (selection == "elety") {
32 domain <- args[4]
33 output <- args[5]
34 rmsd_plot <- args[6]
35 rmsd_hist <- args[7]
36 inds <- atom.select(pdb, elety = domain)
37 }
38 if (selection == "resid") {
39 domain <- args[4]
40 output <- args[5]
41 rmsd_plot <- args[6]
42 rmsd_hist <- args[7]
43 inds <- atom.select(pdb, resid = domain)
44 }
45 if (selection == "segid") {
46 domain <- args[4]
47 output <- args[5]
48 rmsd_plot <- args[6]
49 rmsd_hist <- args[7]
50 inds <- atom.select(pdb, segid = domain)
51 }
52
53 xyz <- fit.xyz(fixed=pdb$xyz, mobile=dcd, fixed.inds=inds$xyz, mobile.inds=inds$xyz)
54
55 rd <- rmsd(xyz[1,inds$xyz], xyz[,inds$xyz])
56
57 write.table(rd, file = output, row.names = TRUE, col.names = FALSE, quote =FALSE, sep="\t")
58
59 png(rmsd_plot)
60 plot(rd, typ="l", ylab="RMSD (Å)", xlab="Frame No.")
61 points(lowess(rd), typ="l", col="red", lty=2, lwd=2)
62 dev.off()
63
64 png(rmsd_hist)
65 hist(rd, breaks=40, freq=FALSE, main="RMSD Histogram", xlab="RMSD")
66 lines(density(rd), typ="l", col="red", lty=2, lwd=2)
67 dev.off()