R Chemistry in 3D

The three dimensional model of the acetylcholine molecule shown below is interactive, when your browser supports WebGL. You can use your mouse to rotate and zoom the model. In this post I will show you the main steps required to create such an interactive 3D model in R: easy peasy (more or less).

Your browser does not support the HTML5 canvas element.

Drag mouse to rotate model. Use mouse wheel or middle button to zoom it.

There are several R packages available for cheminformatics. However, as far as I could tell, there aren't any packages that are directly able to create 3D visualisations of molecules. Therefore, we need the capabilities of the rcdk package to load molecule files and 3D rendering capabilities of the rgl package.

Although it is possible to derive a 3D structure from a SMILES code (see e.g., the SMILES translator from the Computer-Aided Drug Design Group), I couldn't find such functionality in R packages. Therefore, we use an existing 3D model of the acetylcholine molecule from PubChem in this example. The first part of the script downloads the molecule and then loads it into the working memory.

	# required packages
	require(rcdk)
	require(rgl)
	# download 3D structure of acetylcholine from Pubchem:
	download.file("https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/cid/187/record/SDF/?record_type=3d&response_type=save&response_basename=Structure3D_CID_187",
	"acetylcholine.sdf", mode = "wb")
	# load the downloaded file:
	molfile <- load.molecules("acetylcholine.sdf")[[1]]

view raw molecule3d01.r hosted with ❤ by GitHub

Subsequently, we define some atom properties, which we will need for plotting later on: the atom symbol, the colour to be used in the model, and the Van der Waals radii of the atoms as listed at Wikipedia. Note that we only define these properties for the four atoms present in acetylcholine. If you want to apply the code to other molecules, you have to make sure all relevant atoms are specified. In our example we continue by getting the 3D atom coordinates from the molecule file, sticking it in a data.frame and merge it with the previously defined atom properties.

	# create a data.frame with atom properties
	atom.props <- data.frame(
	symbol = c("O", "N", "C", "H"),
	vanderwaals = c(1.52, 1.55, 1.7, 1.2), #Van der Waals radii of atoms
	colour = c("red", "blue", grey(0.1), "white")
	)
	# get atoms from the molecule
	atoms <- get.atoms(molfile)
	# get 3D locations of atoms as read from the file
	atom_dat <- data.frame(do.call('rbind', lapply(atoms, get.point3d)))
	names(atom_dat) <- c("x", "y", "z")
	# merge with atom info
	atom_dat$symbol <- unlist(lapply(atoms, get.symbol))
	atom_dat <- merge(atom_dat, atom.props, all.x = T, all.y = F)

view raw molecule3d02.r hosted with ❤ by GitHub

Now, all we have to do is to plot spheres at the positions of the atoms, with the right sizes and colours. For this we use the spheres3d function. We could also define our own (smoother) mesh of a sphere, but that will result in huge html files when saving the model with writeWebGL. We therefore stick with the lighter spheres3d function. Unfortunately, in Internet Explorer strange shadows appear (a bug?), but in other browsers it looks fine. Just see for yourself at the top of this page.

	# simply plot the atoms as coloured spheres.
	# multiply the Van der Waals radii with 0.75
	# to reduce overlap of atoms.
	mapply(function(x, y, z, vdw, col){
	spheres3d(x, y, z, + vdw*0.75, col = col, add = T)
	}, atom_dat$x, atom_dat$y, atom_dat$z,
	atom_dat$vanderwaals, atom_dat$colour)
	# Save as an interactive model on web page
	writeWebGL()
	# Open in default browser:
	browseURL(paste("file://", file.path(getwd(), "webGL/index.html"), sep = ""))
	# Tested in Chrome on Windows 7 machine with success
	# Tested in IE on Windows 7 machine with limited success
	# Tested in Chrome on Android device without success

view raw molecule3d03.r hosted with ❤ by GitHub

Although writeWebGL produces an html file, which is ready to go and can be uploaded to your website directly, it may need some editing when embedding it in a custom template. This requires some knowledge of the html and JavaScript languages. I won't be explaining the modifications I made here.

As a bonus, if your more a fan of the ball and stick model, below the code to add sticks to your molecule model.

	#################################
	# ball and stick model
	#################################
	# get bond info
	bonds <- get.bonds(molfile)
	bonds.atoms <- data.frame(do.call('rbind', lapply(bonds, function(x){
	atoms <- get.atoms(x)
	start <- get.point3d(atoms[[1]])
	end <- get.point3d(atoms[[2]])
	return (list(start.x = start[1], start.y = start[2], start.z = start[3],
	end.x = end[1], end.y = end[2], end.z = end[3]))
	})))
	# open rgl device
	open3d()
	# plot sticks
	mapply(function(x1, y1, z1, x2, y2, z2){
	shade3d(cylinder3d(data.frame(x = c(x1, x2), y = c(y1, y2), z = c(z1, z2)),
	radius = 0.1, sides = 16,
	e2=rbind(c(1,0,0),c(1,0,0))), col="grey")
	}, bonds.atoms$start.x, bonds.atoms$start.y, bonds.atoms$start.z,
	bonds.atoms$end.x, bonds.atoms$end.y, bonds.atoms$end.z)
	# plot balls
	mapply(function(x, y, z, vdw, col){
	spheres3d(x, y, z, vdw*0.25, col = col, add = T)
	}, atom_dat$x, atom_dat$y, atom_dat$z,
	atom_dat$vanderwaals, atom_dat$colour)

view raw molecule3d04.r hosted with ❤ by GitHub