Mercurial > repos > chemteam > bio3d_dccm
view pca.R @ 2:c9ae0c833ded draft default tip
"planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/tree/master/tools/bio3d commit 3b99f08f22b9e0c16c0a0adc82f8c16c1a25cedf"
author | chemteam |
---|---|
date | Mon, 07 Oct 2019 12:44:45 -0400 |
parents | 386b4f4b6a85 |
children |
line wrap: on
line source
#!/usr/bin/env Rscript options(stringAsfactors = FALSE) args <- commandArgs(trailingOnly = TRUE) library(bio3d) dcdfile <- args[1] pdbfile <- args[2] dcd <- read.dcd(dcdfile) pdb <- read.pdb(pdbfile) method <- args[3] selection <- args[4] domain <- args[5] output <- args[6] pca_plot <- args[7] pca_cluster <- args[8] pc1_rmsf <- args[9] if (selection == "string") { inds <- atom.select(pdb, string = domain) } if (selection == "elety") { inds <- atom.select(pdb, elety = domain) } if (selection == "resid") { inds <- atom.select(pdb, resid = domain) } if (selection == "segid") { inds <- atom.select(pdb, segid = domain) } xyz <- fit.xyz(fixed=pdb$xyz, mobile=dcd, fixed.inds=inds$xyz, mobile.inds=inds$xyz) if (method == "FALSE") { pc <- pca.xyz(xyz[,inds$xyz], use.svd=FALSE) } if (method == "TRUE") { pc <- pca.xyz(xyz[,inds$xyz], use.svd=TRUE) } write.table(pc$au[,1:2:3], file = output, row.names = TRUE, col.names = FALSE, quote =FALSE, sep="\t") png(pca_plot) plot(pc, col=bwr.colors(nrow(xyz)) ) dev.off() png(pca_cluster) hc <- hclust(dist(pc$z[,1:2])) grps <- cutree(hc, k=2) plot(pc, col=grps) dev.off() png(pc1_rmsf) plot.bio3d(pc$au[,1], ylab="PC1 (A)", xlab="Residue Position", typ="l") points(pc$au[,2], typ="l", col="blue") dev.off()