Plotting in R
This page will walk you through loading CHAP output data into R and show you how to visualise the radius profile of the permeation pathway. It also introduces a set of turnkey scripts bundled with CHAP that you can use to quickly visualise CHAP data. These scripts are a good starting point for getting an overview of the data CHAP generates and can be adapted to suit your needs. It is assumed that you have some basic understanding of the R language and the ggplot2 plotting system.
Example: Plotting Time-averaged Pathway Profiles
As CHAP writes all relevant data into a single JSON file, you need to load a JSON parser first. There are several available in R; here we will use
jsonlite, which can be installed by typing
output.json file from CHAP can then be loaded via
# load JSON parser: library(jsonlite) # load CHAP data: dat <- fromJSON(readLines("output.json", n = 1), flatten = FALSE)
flatten is set to
FALSE to preserve the structure of the JSON data. You can now use
names(dat) to get an overview of the data contained in
The data on the time-averaged radius profile can be found in
dat$pathwayProfile and can be plotted using
ggplot by typing (make sure to load
plt <- ggplot(data = as.data.frame(dat$pathwayProfile), aes(x = s, y = radiusMean)) + geom_line() + xlab("s (nm)") + ylab("R (nm)")
This will create a plot object named
plt, which only contains the mean radius profile averaged over all trajectory frames, but it is very easy to add a confidence interval:
plt + geom_ribbon(aes(x = s, ymin = radiusMean - radiusSd, ymax = radiusMean + radiusSd), alpha = 0.2)
The above will create a grey area around the mean line, which indicates the one standard deviation confidence interval.
In addition to the radius profile, it is interesting to also look at the residues lining the permeation pathway. A data frame for the residue data in the CHAP output can be created by typing
res.data <- as.data.frame(dat$residueSummary)
and the residue positions can be added to the radius profile via:
plt + geom_point(data = res.data[res.data$poreFacing.mean > 0.5,], aes(x = s.mean, y = rho.mean, colour = hydrophobicity)) + scale_colour_distiller(palette = "BrBG", name = "H (a.u.)", limits = c(-max(abs(data$hydrophobicity)), max(abs(data$hydrophobicity))))
Note that this will only plot residues that are facing the pore more than 50% of the time and that residues are coloured by their respective hydrophobicity.
The above commands will generate a plot that (together with some further visual tweaking) will look like the following:
A similar procedure can be used to create plots of other pathway properties such as its hydrophobicity or the local solvent density. Take a look at
chap/scripts/plotting/R/chap_plot_pathway_profile.R for an overview of what other pathway properties can be plotted.
Turnkey Plotting Scripts
To get a quick overview of the results of a CHAP run you can use the plotting scripts available under
chap/scripts/plotting/R/. There are three scripts,
chap_plot_profile_time_series.R, each of which generates a set of figures that display time-averaged pathway profiles, scalar-valued time series, and vector-valued time series respectively.
The scripts require R to be installed and the
Rscript command to be available. They also need the
ggplot2 libraries and will install them automatically from CRAN if they are not found. The easiest way to use the scripts is to copy them into your working directory, where you can run them by simply typing
./chap_plot_pathway_profile.R -filename output.json -dpi 300 ./chap_plot_scalar_time_series.R -filename output.json -dpi 300 ./chap_plot_profile_time_series.R -filename output.json -dpi 300
This will read data from
output.json and create PNG figures with a resolution of 300 dots per inch. Both arguments are optional and default to
300 respectively. To avoid copying the scripts to your working directoy each time you generate new data, you may consider adding the script location to your PATH. Of course, you can also open and run the scripts from an R IDE such as RStudio.