RE: [AMBER] How to obtain dialsplot

From: Niel Henriksen <niel.henriksen.utah.edu>
Date: Wed, 12 May 2010 10:11:44 -0600

I use the program "R". It can be a somewhat difficult process to get the installation working nicely.
There is a bit of a learning curve and easier programs probably exist, but I haven't found one with the flexibility of R.

If you're still interested:

1) download and install R
2) Add the plotrix, circular, and cairo packages
3) Create output files of dihedral vs. time (I usually do all the dihedrals for nucleic acids)
3) Example run command: /path/to/Rbin/R CMD BATCH --no-save --no-restore "--args ./ radials.png" input-script.R

Below is an input script for a 34 residue RNA sequence. The output is a large .png with 34 X 9 dials (alpha,beta,gamma,delta,epsilon,zeta,chi,nu1,nu2). The input consists of 34 X 9 files named a1.dat, b1.dat ... a34.dat, b34.dat, etc. You must have the file present with numbers in it, even if the dihedral doesn't exist (ie, put "0 0" in a1.dat).
There are all kinds of modifications you can make to the script to suit your needs. You'll have to read the (extensive) documentation.
Good luck
--Niel

(kudos to Jordan Monnet for figuring this out)
input-script.R:
############# BEGINNING OF THE SCRIPT #####################

args = (commandArgs(TRUE))

##### LOADING PACKAGES ######
require(plotrix)
require(circular)
require(Cairo)
#############################

######### OUTPUT ############
CairoPNG(args[2],width=1800,height=6800)
par(mfrow=c(34,9))
# I assume we need 200*200px for each plot
#############################

for(j in 1:34){ # loop on bases
  print(j)
  for(i in c("a","b","g","d","e","z","c","nu1-","nu2-")){ # loop on parameters
    ########## DATA #############
    dihe = paste(i,j,sep="")
    file = paste(args[1],dihe,sep="")
    data = read.table(paste(file,"dat",sep="."),colClasses="numeric",h=F)

    colnames(data) = c("Time", "Angle")
    #############################

    ######## PLOT ALL ###########
    par(cex.lab=0.4, cex.axis=0.8, cex.main=1.4, cex.sub=0.4,new=F)
    radial.plot(data$Time,
        rad(data$Angle), # convert degree into radian
        labels = c("0",90, "-/+ 180", -90),
        label.pos = c(0,pi/2,pi,3*pi/2),
        start = pi/2,
        clockwise = TRUE,
        rp.type = "p",
        lwd = 0.005,
        show.centroid = FALSE,
        show.radial.grid = FALSE,
        show.grid.label = FALSE,
        show.grid = FALSE,
        main = dihe)
    #############################
  }
}

#############################
dev.off()
#############################

###########################################################
# ..Jordan Monnet.. }><(({*> #
# __/ #
###########################################################


________________________________________
From: amber-bounces.ambermd.org [amber-bounces.ambermd.org] On Behalf Of Homa Azizian [homa.azizian.anu.edu.au]
Sent: Tuesday, May 11, 2010 8:10 PM
To: AMBER Mailing List
Subject: [AMBER] How to obtain dialsplot

Hi amber users,
I want to creat dialsplot of chi1 and chi2 angle. I could not find how to produce it after geting the out put file from ptraj.
Thank's for replying in advance
Homa
-----
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber

_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed May 12 2010 - 09:30:05 PDT
Custom Search