[gmx-users] ramachandran plot

Mon Nov 17 13:55:19 CET 2014

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.