Mercurial > repos > chemteam > bio3d_rmsf

comparison rmsf.R @ 0:e838317708a6 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:50:23 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:e838317708a6
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 rmsf_plot <- args[6]
20 inds <- atom.select(pdb, string = domain)
21 }
22 if (selection == "resno") {
23 res1 <- args[4]
24 res2 <- args[5]
25 output <- args[6]
26 rmsf_plot <- args[7]
27 inds <- atom.select(pdb, resno=res1:res2)
28 }
29 if (selection == "elety") {
30 domain <- args[4]
31 output <- args[5]
32 rmsf_plot <- args[6]
33 inds <- atom.select(pdb, elety = domain)
34 }
35 if (selection == "resid") {
36 domain <- args[4]
37 output <- args[5]
38 rmsf_plot <- args[6]
39 inds <- atom.select(pdb, resid = domain)
40 }
41 if (selection == "segid") {
42 domain <- args[4]
43 output <- args[5]
44 rmsf_plot <- args[6]
45 inds <- atom.select(pdb, segid = domain)
46 }
47
48 xyz <- fit.xyz(fixed=pdb$xyz, mobile=dcd, fixed.inds=inds$xyz, mobile.inds=inds$xyz)
49
50 rf <- rmsf(xyz[,inds$xyz])
51
52 write.table(rf, file = output, row.names = TRUE, col.names = FALSE, quote =FALSE, sep="\t")
53
54 png(rmsf_plot)
55 plot(rf, ylab="RMSF", xlab="Residue Position", typ="l")
56 dev.off()
57