Hi Tsjerk, I was wondering how did you get the input.dat file? did you got it from grace/xmgrace? I was able to save the plot in .ps or .png format file. However whenever I tried to run kde2d.r I was getting an error: Please see below: ./kde2d.r ramachandran.png ramachandran-2.png ./kde2d.r: line 11: syntax error near unexpected token `(' ./kde2d.r: line 11: `args <- commandArgs()'
what could be wrong? best regards Urszula > Hi Urszula, > > Save the following as kde2d.r and make it executable. Then you can run > > ./kde2d.r rama.dat rama.png > > Cheers, > > Tsjerk > > ####kde2d.r#### > > #!/usr/bin/env Rscript > > # Script for drawing 2D combined linear/circular KDE plots. > # > # (c)2014 Tsjerk A. Wassenaar > # Friedrich-Alexander University of Erlangen-Nuremberg > > > args <- commandArgs() > args <- args[-(1:match("--args", args))] > > > # MASS library is needed for bandwidths > require(MASS) > > > # This 2D circular KDE function was developed for correct handling > # of bivariate angular data as part of the orientational analysis > # used with DAFT (Manuscript submitted). > kde2d <- function (x, y, h, n = 25, xlim, ylim, circular=TRUE, phase=0) > { > if (length(y) != length(x)) > stop("data vectors must be the same length") > > n <- rep(n, length.out = 2L) > phase <- rep(phase, length.out = 2L) > circular <- rep(circular, length.out = 2L) > > s <- which(!is.na(x) | !is.na(y)) > x <- x[s] > y <- y[s] > > nx <- length(x) > > if (any(!is.finite(x)) || any(!is.finite(y))) > stop("missing or infinite values in the data are not allowed") > > > if (circular[1]) > x <- (x+180)%%360-180 > > if (circular[2]) > y <- (y+180)%%360-180 > > > if (missing(xlim)) > xlim <- if (circular[1]) c(-180,180) else range(x) > > if (missing(ylim)) > ylim <- if (circular[2]) c(-180,180) else range(y) > > if (any(!is.finite(c(xlim,ylim)))) > stop("only finite values are allowed in 'xlim' and 'ylim'") > > > h <- if (missing(h)) > c(bandwidth.nrd(x), bandwidth.nrd(y)) > else rep(h, length.out = 2L) > if (any(h <= 0)) > stop("bandwidths must be strictly positive") > h <- h/4 > > > if (circular[1]) > { > gx <- seq.int(-180, 180, length.out = n[1L]) > ax <- ((outer(gx, x-phase[1], "-")+180)%%360-180)/h[1L] > gx <- gx + phase[1] > } > else > { > gx <- seq.int(xlim[1], xlim[2], length.out = n[1L]) > ax <- outer(gx, x, "-")/h[1L] > } > > > if (circular[2]) > { > gy <- seq.int(-180, 180, length.out = n[2L]) > ay <- ((outer(gy, y-phase[2], "-")+180)%%360-180)/h[2L] > gy <- gy + phase[2] > } > else > { > gy <- seq.int(ylim[1], ylim[2], length.out = n[2L]) > ay <- outer(gy, y, "-")/h[2L] > } > > > z <- tcrossprod(matrix(dnorm(ax), , nx), matrix(dnorm(ay), , nx))/(nx > * > h[1L] * h[2L]) > > > list(x = gx, y = gy, z = z) > } > > > # Color gradients for coloring density plots > > rng <- seq(0,1,length.out=100) > > white2red <- rgb(1,rev(rng),rev(rng)) > blue2white <- rgb(rng,rng,1) > > # Rainbow > magenta2red <- rgb(1,0,rev(rng)) > red2yellow <- rgb(1,rng,0) > yellow2green <- rgb(rev(rng),1,0) > green2cyan <- rgb(0,1,rng) > cyan2blue <- rgb(0,rev(rng),1) > > red2orange <- rgb(1,0.5*rng,0) > orange2yellow <- rgb(1,0.5+0.5*rng,0) > > > > data <- read.table(args[1]) > d <- kde2d(data[,1],data[,2]) > > > ## Plotting density images > > # Here the density images are prepared and written to file. The first one > (WT) is explained in detail: > > #---Start of image--- > > # Wildtype > > # Start a 1200x1200 png file, with a font size (for labels) of 40 points > # Changing the resolution requires changing the font size to keep the > relative size. > png( > args[2], > width=1200, > height=1200, > pointsize=40 > ) > > # Write an image of the density for WT > image( > d$x, d$y, d$z, # Image data X/Y and intensity (Z) > zlim=c(0,max(d$z)), # Limits for coloring > xlab=expression(phi), # X-label text. expression(beta) gives the > greek > character > ylab=expression(psi), # Y-label text. > main="", # Main title. Can be set to "" to suppress > title. > bty="n", # Surrounding box type. To draw a thicker box, > the box is here suppressed, using "n" for None. > col=c(white2red,red2orange,orange2yellow) > ) > > box(lwd=3) # Add a box with a thicker line. Change this to > match the resolution of the imge. > # Add points > points(data[,1], data[,2], pch=".", cex=1) > > # Add contour lines > contour( > d$x, d$y, d$z, > zlim=c(0,max(d$z)), > add=TRUE, # If add=FALSE (default), then 'contour' starts > a new plot. > labels="", # Suppress labels on contour lines > labcex=0.001, # To avoid gaps in the controu lines (because > of > the empty label), set the label size very small > nlevels=20, # Number of lines to draw, between range of > 'zlim' > lwd=3 # Width of contour lines. Set to 3 to match png > resolution. Change when changing the resolution. > ) > > # Close the image file, saving it. > dev.off() > > #---End of image--- > > > > > On Mon, Nov 17, 2014 at 1:44 PM, Urszula Uciechowska < > urszula.uciechow...@biotech.ug.edu.pl> wrote: > >> Dear Tsjerk, >> >> Could you send me the script again? as I did not find the attachment in >> your previous email. >> >> best regards >> Urszula >> >> >> > Hi Urszula, >> > >> > Attached is an R script for drawing 2D circular or combined >> > circular/linear >> > KDEs. It works as a command line script, taking as arguments a file >> with >> > 2-column data, and an output image file name (./kde2d.r input.dat >> > output.png). It's pretty easy to adapt to suit your needs, even if you >> > don't know R specifically. I wrote the routine for a manuscript we >> just >> > submitted, where the method is explained in some detail, and if you >> like >> > the images and care to use it for publications, I would be obliged if >> you >> > could reference that paper ("High-Throughput Simulations of Dimer and >> > Trimer Assembly of Membrane Proteins. The DAFT approach."). If you run >> > into >> > problems or have further questions, please let me know. >> > >> > Cheers, >> > >> > Tsjerk >> > >> > On Fri, Nov 14, 2014 at 4:31 PM, Urszula Uciechowska < >> > urszula.uciechow...@biotech.ug.edu.pl> wrote: >> > >> >> Could you please send me the code? It can be on my private email. >> >> >> >> Thanks a lot >> >> >> >> Ursuzla >> >> >> >> > Hi Urszula, >> >> > >> >> > It's a while ago that I made that one, and I don't have the code at >> >> hand. >> >> > But it's a combination of a density plot (kde2d) with the points >> laid >> >> > over, >> >> > and a polygon to highlight the forbidden region. These days I'm >> doing >> >> > circular 2D KDEs, whic is more correct. I can send the R code for >> thay >> >> if >> >> > you want. >> >> > >> >> > Best, >> >> > >> >> > Tsjerk >> >> > On Nov 14, 2014 12:34 PM, "Urszula Uciechowska" < >> >> > urszula.uciechow...@biotech.ug.edu.pl> wrote: >> >> > >> >> >> Dear Gromacs users, >> >> >> >> >> >> I generated the ramachandran plot and would like to colour it in >> >> xmgrace >> >> >> as it was shown in tutorial: >> >> >> http://nmr.chem.uu.nl/~adrien/course/molmod/analysis2.html >> >> >> >> >> >> Does anyone know how to do it? >> >> >> >> >> >> Best regards >> >> >> Urszula Uciechowska >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> ----------------------------------------- >> >> >> Ta wiadomość została wysłana z serwera Uniwersytetu Gdańskiego >> >> >> http://www.ug.edu.pl/ >> >> >> >> >> >> -- >> >> >> Gromacs Users mailing list >> >> >> >> >> >> * Please search the archive at >> >> >> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before >> >> >> posting! >> >> >> >> >> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >> >> >> >> >> >> * For (un)subscribe requests visit >> >> >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users >> or >> >> >> send a mail to gmx-users-requ...@gromacs.org. >> >> >> >> >> > -- >> >> > Gromacs Users mailing list >> >> > >> >> > * Please search the archive at >> >> > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before >> >> > posting! >> >> > >> >> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >> >> > >> >> > * For (un)subscribe requests visit >> >> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users >> or >> >> send >> >> > a mail to gmx-users-requ...@gromacs.org. >> >> > >> >> >> >> >> >> University of Gdansk and Medical Univesity of Gdansk >> >> Department of Molecular and Cellular Biology >> >> ul. Kladki 24 >> >> 80-822 Gdansk >> >> Poland >> >> >> >> >> >> ----------------------------------------- >> >> Ta wiadomość została wysłana z serwera Uniwersytetu Gdańskiego >> >> http://www.ug.edu.pl/ >> >> >> >> -- >> >> Gromacs Users mailing list >> >> >> >> * Please search the archive at >> >> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before >> >> posting! >> >> >> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >> >> >> >> * For (un)subscribe requests visit >> >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or >> >> send a mail to gmx-users-requ...@gromacs.org. >> >> >> > >> > >> > >> > -- >> > Tsjerk A. Wassenaar, Ph.D. >> > -- >> > Gromacs Users mailing list >> > >> > * Please search the archive at >> > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before >> > posting! >> > >> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >> > >> > * For (un)subscribe requests visit >> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or >> send >> > a mail to gmx-users-requ...@gromacs.org. >> > >> >> >> University of Gdansk and Medical Univesity of Gdansk >> Department of Molecular and Cellular Biology >> ul. Kladki 24 >> 80-822 Gdansk >> Poland >> >> >> ----------------------------------------- >> Ta wiadomość została wysłana z serwera Uniwersytetu Gdańskiego >> http://www.ug.edu.pl/ >> >> -- >> Gromacs Users mailing list >> >> * Please search the archive at >> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before >> posting! >> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >> >> * For (un)subscribe requests visit >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or >> send a mail to gmx-users-requ...@gromacs.org. >> > > > > -- > Tsjerk A. Wassenaar, Ph.D. > -- > Gromacs Users mailing list > > * Please search the archive at > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before > posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send > a mail to gmx-users-requ...@gromacs.org. > University of Gdansk and Medical Univesity of Gdansk Department of Molecular and Cellular Biology ul. Kladki 24 80-822 Gdansk Poland ----------------------------------------- Ta wiadomość została wysłana z serwera Uniwersytetu Gdańskiego http://www.ug.edu.pl/ -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.