Mercurial > repos > chemteam > bio3d_pca_visualize

comparison pca.R @ 0:65bfd1b90b96 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 cd0830e5e3502721fa355cc8e3fedc331201a6e4
author chemteam
date Tue, 26 Feb 2019 08:29:24 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:65bfd1b90b96
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
11 dcd <- read.dcd(dcdfile)
12 pdb <- read.pdb(pdbfile)
13
14 method <- args[3]
15 selection <- args[4]
16 domain <- args[5]
17
18 output <- args[6]
19 pca_plot <- args[7]
20 pca_cluster <- args[8]
21 pc1_rmsf <- args[9]
22
23
24 if (selection == "string") {
25 inds <- atom.select(pdb, string = domain)
26 }
27 if (selection == "elety") {
28 inds <- atom.select(pdb, elety = domain)
29 }
30 if (selection == "resid") {
31 inds <- atom.select(pdb, resid = domain)
32 }
33 if (selection == "segid") {
34 inds <- atom.select(pdb, segid = domain)
35 }
36 xyz <- fit.xyz(fixed=pdb$xyz, mobile=dcd, fixed.inds=inds$xyz, mobile.inds=inds$xyz)
37
38 if (method == "FALSE") {
39 pc <- pca.xyz(xyz[,inds$xyz], use.svd=FALSE)
40 }
41 if (method == "TRUE") {
42 pc <- pca.xyz(xyz[,inds$xyz], use.svd=TRUE)
43 }
44
45 write.table(pc$au[,1:2:3], file = output, row.names = TRUE, col.names = FALSE, quote =FALSE, sep="\t")
46
47 png(pca_plot)
48 plot(pc, col=bwr.colors(nrow(xyz)) )
49 dev.off()
50
51 png(pca_cluster)
52 hc <- hclust(dist(pc$z[,1:2]))
53 grps <- cutree(hc, k=2)
54 plot(pc, col=grps)
55 dev.off()
56
57 png(pc1_rmsf)
58 plot.bio3d(pc$au[,1], ylab="PC1 (A)", xlab="Residue Position", typ="l")
59 points(pc$au[,2], typ="l", col="blue")
60 dev.off()
61