[gmx-users] ramachandran plot
Urszula Uciechowska
urszula.uciechowska at biotech.ug.edu.pl
Tue Nov 18 09:30:09 CET 2014
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.uciechowska at 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.uciechowska at 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.uciechowska at 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-request at 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-request at 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-request at 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-request at 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-request at 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-request at 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/
More information about the gromacs.org_gmx-users
mailing list