Network visualization – part 4: 3D networks

In the fourth and final part of my graph visualization series, I’ll show how to create 3D network plots. 3D plots are more than just pretty plots – they allow you to rotate, scale, and zoom in and out of the network. These options may help identify interesting interaction patterns. Unfortunately, there are not many 3D network visualization tools (especially not free ones). I’m also not aware of any tools that have a R library.

So how are we going to create a 3D plot from R?

Well, we need to be clever: we will pretend that our graph represents a chemical structure and use Jmol, an open-source 3D viewer for chemical structures, to visualize it. As there is no direct link between R and Jmol, the only way to visualize our network is to create a corresponding Jmol file (.mol2 file) in R and open it in Jmol. This is similar to what we have done for Gephi.

As before, we will use the weighted network of characters’ coappearances in Victor Hugo’s novel “Les Miserables” (LesMiserables.txt). Unfortunately, the number of molecule (graph) properties that we can load in Jmol directly from .mol2 file is limited, so we will only use one property – node degree – and we will use it to define a node size. Note that various molecule (graph) properties in Jmol can be set through Jmol scripts (see below), but we won’t focus on those in this blog.

First, let’s load our network and calculate degrees of all nodes:

library("igraph") rm(list = ls()) ########################################################################## # Read data dataSet <- read.table("lesmis.txt", header = FALSE, sep = "t") # Create a graph. Use simplify to ensure that there are no duplicated edges or self loops gD <- simplify(graph.data.frame(dataSet, directed=FALSE)) # Calculate degree for all nodes degAll <- degree(gD, v = V(gD), mode = "all")

Next, we'll calculate the node coordinates using layout function from igraph library:

coord3D < - layout.fruchterman.reingold(gD, dim = 3)

Now, we are ready to make the .mol2 file. We'll call it 3D_lesmis.mol2. We'll use the igraph library functions vcount() and ecount(gD) to compute the number of nodes and edges. For more details about .mol2 file format, see: mol2.pdf

nType <- 1 write.table("@MOLECULE", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE) write.table("Les Miserables", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE) write.table(data.frame(vcount(gD), ecount(gD), nType), file = "3D_lesmis.mol2", sep = "t", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE) write.table("SMALL", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE) write.table("NO_CHARGESn", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE) write.table("@ATOM", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE)

Next, we will add a list of nodes into our 3D_lesmis.mol2 file. Each node is represented by an atom, and the atom type will define the color of the node. We will define nodes as follows: nodes with degree one will be represented as helium atoms (light blue); degree two as sodium atoms (purple); degree 3-5 as oxygen atoms (red); degree 6-10 as gold (yellow); degree 11-15 as phosphorus (orange); degree larger than 15 as chlorine atoms (green). For a full color scheme, see: jscolors. In this step, we will also assign node coordinates to each atom:

isType <- 1 for (i in 1:vcount(gD)) { if (degAll[i] == 1) isIn <- "He" if (degAll[i] == 2) isIn <- "Na" if ((degAll[i] > 2) & (degAll[i] 5) & (degAll[i] 10) & (degAll[i] 15) isIn <- "Cl" hlpL <- data.frame(i, V(gD)[i]$name, coord3D[i, 1], coord3D[i, 2], coord3D[i, 3], isIn, isType, 0.0 ) write.table(hlpL, file = "3D_lesmis.mol2", sep = "t", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE) }

Finally, we will add a list of edges in the 3D_lesmis.mol2 file:

write.table("@BOND", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE) isType <- 1 for (i in 1:ecount(gD)) { hlpL <- data.frame(i, get.edge(gD,i)[1], get.edge(gD,i)[2], isType) write.table(hlpL, file = "3D_lesmis.mol2", sep = "t", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE) }

Here is the resulting file: 3D_lesmis.mol2. The resulting network is shown in Figure 1A:

Figure 1

By default, the color of the edges are inferred from the color of the nodes. We can change that easily as follows: do a right click anywhere in the network window to open a pop-up menu. Then click on "Color," then "Bonds," and then select the color of choice (I used white). Now all edges are the same color (white), as shown in Figure 1B. I find the edges too thick compared to nodes sizes. To change edge thickness, go to the "Display" menu, select "Bond," and then the type you want (I chose "Wireframe"). As we can see in Figure 1C, our network looks nice. Unsurprisingly (as we used the same layout), it looks similar to the networks we created in Cytoscape and Gephi. However, here we can rotate the network and see the interactions from various angles and detail levels without a need to create multiple network plots. There is one more thing we are missing - node labels. They are a part of the .mol2 file, but are not visible by default. To show labels, go to the "Display" menu, select "Label," and then select "Name." Super easy! The resulting network is shown in the Figure 2A.

Jmol allows us to do additional customizations with Jmol scripts. The scripts are very intuitive and they can be used to define different node (atom) and edge (bond) sizes, edge colors, as well as the node label size. To use Jmol scripts, we need Jmol console. You can find it at "File" menu -> "Console."

When using Jmol scripts, the first thing to do is to select a node (or edge) on which we want to work. When a node/edge is selected, we can define its attributes. For example, to change the size of helium atoms, we use a command select: select _He" (there must be an underscore before the atom symbol). Then, we can define its size with a command spacefill. We can also define the edge thickness between any two helium atoms using the wireframe command, as well as the color of this edge with the color bonds command. Finally, we can define the size of the label accompanying the atom/node with the label command.

So let's how it works.

select _He; spacefill 0.2; wireframe 10; color bonds blue; label 10; select _Na; spacefill 0.4; wireframe 12; color bonds purple; label 12; select _O; spacefill 0.6; wireframe 14; color bonds red; label 14; select _Au; spacefill 0.8; wireframe 16; color bonds yellow; label 16; select _P; spacefill 1.0; wireframe 18; color bonds orange; label 18; select _Cl; spacefill 1.5; wireframe 20; color bonds green; label 20;

Figure 2B shows the results of the script above (with zoomed in). In case that we are not interested in all labels, e.g., for nodes that do not have many interacting partners, we can remove some of the labels using the label hide command (Figure 2C and 2D):

select _He; label hide; select _Na; label hide; select _O; label hide

Figure 2

#############################
I've noticed that some parts of the code are not displaying property, so here is a full version of the code (available at gist)

library("igraph")
rm(list = ls())
##########################################################################
# Read data
dataSet <- read.table("lesmis.txt", header = FALSE, sep = "t")

# Create a graph. Use simplify to ensure that there are no duplicated edges or self loops
gD <- simplify(graph.data.frame(dataSet, directed=FALSE))

# Calculate degree for all nodes
degAll <- degree(gD, v = V(gD), mode = "all")

coord3D <- layout.fruchterman.reingold(gD, dim = 3) 
nType <- 1

write.table("@<TRIPOS>MOLECULE", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE)
write.table("Les Miserables", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE)
write.table(data.frame(vcount(gD), ecount(gD), nType), file = "3D_lesmis.mol2", sep = "t", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE)
write.table("SMALL", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE)
write.table("NO_CHARGESn", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE)
write.table("@<TRIPOS>ATOM", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE)

isType <- 1
for (i in 1:vcount(gD))
{
  if (degAll[i] == 1)
    isIn <- "He"
  
  if (degAll[i] == 2)
    isIn <- "Na"
  
  if ((degAll[i] > 2) & (degAll[i] <= 5))
    isIn <- "O"
  
  if ((degAll[i] > 5) & (degAll[i] <= 10))
    isIn <- "Au"
  
  if ((degAll[i] > 10) & (degAll[i] <= 15))
    isIn <- "P"
  
  if (degAll[i] > 15)
    isIn <- "Cl"
  
  hlpL <- data.frame(i, V(gD)[i]$name, coord3D[i, 1], coord3D[i, 2], coord3D[i, 3], isIn, isType, 0.0 )
  
  write.table(hlpL, file = "3D_lesmis.mol2", sep = "t", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE)
}

write.table("@<TRIPOS>BOND", file = "3D_lesmis.mol2", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE)
isType <- 1

for (i in 1:ecount(gD))
{
  hlpL <- data.frame(i, get.edge(gD,i)[1], get.edge(gD,i)[2], isType)
  
  write.table(hlpL, file = "3D_lesmis.mol2", sep = "t", row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE)
}


################################
# Visualization adjustments using jmol scripts
#
# _He; spacefill 0.2; wireframe 10; color bonds blue;  label 10;
# select _Na; spacefill 0.4; wireframe 12; color bonds purple;  label 12;
# select _O; spacefill 0.6; wireframe 14; color bonds red;  label 14;
# select _Au; spacefill 0.8; wireframe 16; color bonds yellow;  label 16;
# select _P; spacefill 1.0; wireframe 18; color bonds orange;  label 18;
# select _Cl; spacefill 1.5; wireframe 20; color bonds green;  label 20;
#
# select _He; label hide;
# select _Na; label hide;
# select _O; label hide
################################