Many projects have collected samples from different regions and from different depths of the ocean. Some, such as the pioneering study by Craig Venter [@venter_environmental_2004] have pioneered metagenomic sequencing and others, such as Tara Oceans Expedition have collected large amounts of data with global ecological questions in mind. On the Tara Oceans (8th and 9th expedition for this vessel ) researchers used a small sailboat outfitted with a lab and filtration supplies to collect samples from many different size fractions of microorganisms in the oceans over three years. They collected these samples to look at the different kinds of microorganisms present in different parts of the oceans to look at their composition and to observe their spatial patterns and distribution.
The scientists collected the samples and then used either targeted sequencing (amplicon approach using primers for specific targets such as the ribosomal genes and then amplifying these targets using PCR) or using metagenomic sequencing (where all the genetic material in a sample is sequenced) of each of the size fractions.
After the sequencing and quality checking of the samples was done, the sequences were taxonomically classified (different approaches for the different targets, see here for the details in @brum_patterns_2015 and @sunagawa_structure_2015). After that the data could be made into a species occurrence table where the rows are different sites and the columns are the observations of the different organisms at each site [@lima-mendez_determinants_2015].
Many of these microbial species in these types of studies have not yet been characterized in the lab. Thus, to know more about the organisms and their interactions, we can observe which ones occur at the same sites or under the same kinds of environmental conditions. One way to do that is by using co-occurrence networks where you examine which organisms occur together at which sites. The more frequently that organisms co-occur at the same site, the stronger the interaction predicted among these organisms. For a review of some of the different kinds of techniques and software for creating interaction networks see: @weiss_correlation_2016.
These kinds of analyses can be useful for studies where the organisms have not yet been characterized in the lab because these analyses can provide insights about the communities and how the organisms within them are interacting. These analyses can be exploratory, so that we can see which organisms warrant further insights and perhaps experimental work . We can also learn about how the overall community is organized (community structure) by looking at some of the network properties (that is the overall way that the organisms are co-occurring and the properties of the network seen this way).
In this analysis we are using a Tara Ocean data and we have data from the bacterial dataset [@sunagawa_structure_2015] and also from the viral dataset [@brum_patterns_2015]. They have been examined in @lima-mendez_determinants_2015 and we have used the original relative abundances to visualize the data. Data were retrieved from: http://www.raeslab.org/companion/ocean-interactome.html
We will run this example using RCy3 [@shannon_rcytoscape:_2013] to drive the visualization of these networks in Cytoscape [@shannon_cytoscape:_2003] using CyREST [@ono_cyrest:_2015].
library(RCy3) library(igraph) library(RColorBrewer)
To run this example Cytoscape software must be running. In Cytoscape we will also need Allegro-plugin for this example.
To begin we create a connection in R that we can use to manipulate the networks and then we will delete any windows that were already in Cytoscape so that we don't use up all of our memory.
cy <- CytoscapeConnection() deleteAllWindows(cy)
We will read in a species co-occurrence matrix that was calculated using Spearman Rank coefficient. (If interested in seeing how this was done please see scripts and the raw data in inst/data-raw)
## scripts for processing located in "inst/data-raw/" prok_vir_cor <- read.delim("./data/virus_prok_cor_abundant.tsv")
There are many different ways to work with graphs in R. We will use both the igraph [@csardi_igraph_2006] and the graph [@gentleman_graph:_2016] package to work with our network with Cytoscape.
The igraph package is used to convert the co-occurrence dataframe into a network that we can send to Cytoscape. In this case our graph is undirected (so "directed = FALSE") since we do not have any information about the direction of the interactions.
graph_vir_prok <- igraph::simplify(igraph::graph.data.frame(prok_vir_cor, directed = FALSE))
Since these are data from small, microscopic organisms that were sequenced using shotgun sequencing, we rely on the classification of the sequences to know what kind of organisms are in the samples. In this case the bacterial viruses (bacteriophage), were classified by Basic Local Alignment Search Tool (BLAST http://blast.ncbi.nlm.nih.gov/Blast.cgi) by searching for their closest sequence in the RefSeq database (see methods in @brum_patterns_2015). The prokaryotic taxonomic classifications were determined using the SILVA database.
phage_id_affiliation <- read.delim("./data/phage_ids_with_affiliation.tsv") bac_id_affi <- read.delim("./data/prok_tax_from_silva.tsv")
In preparation for sending the networks to Cytoscape we will add in the taxonomic data. Some of the organisms do not have taxonomic classifications associated with them so we have described them as "not_class" for not classified. We do that because we have had problems sending "NA"s to Cytoscape from RCy3.
genenet.nodes <- as.data.frame(igraph::vertex.attributes(graph_vir_prok)) ## not all have classification, so create empty columns genenet.nodes$phage_aff <- rep("not_class", nrow(genenet.nodes)) genenet.nodes$Tax_order <- rep("not_class", nrow(genenet.nodes)) genenet.nodes$Tax_subfamily <- rep("not_class", nrow(genenet.nodes)) for (row in seq_along(1:nrow(genenet.nodes))){ if (genenet.nodes$name[row] %in% phage_id_affiliation$first_sheet.Phage_id_network){ id_name <- as.character(genenet.nodes$name[row]) aff_to_add <- unique(subset(phage_id_affiliation, first_sheet.Phage_id_network == id_name, select = c(phage_affiliation, Tax_order, Tax_subfamily))) genenet.nodes$phage_aff[row] <- as.character(aff_to_add$phage_affiliation) genenet.nodes$Tax_order[row] <- as.character(aff_to_add$Tax_order) genenet.nodes$Tax_subfamily[row] <- as.character(aff_to_add$Tax_subfamily) } }
## do the same for proks genenet.nodes$prok_king <- rep("not_class", nrow(genenet.nodes)) genenet.nodes$prok_tax_phylum <- rep("not_class", nrow(genenet.nodes)) genenet.nodes$prok_tax_class <- rep("not_class", nrow(genenet.nodes)) for (row in seq_along(1:nrow(genenet.nodes))){ if (genenet.nodes$name[row] %in% bac_id_affi$Accession_ID){ aff_to_add <- unique(subset(bac_id_affi, Accession_ID == as.character(genenet.nodes$name[row]), select = c(Kingdom, Phylum, Class))) genenet.nodes$prok_king[row] <- as.character(aff_to_add$Kingdom) genenet.nodes$prok_tax_phylum[row] <- as.character(aff_to_add$Phylum) genenet.nodes$prok_tax_class[row] <- as.character(aff_to_add$Class) } }
Add to the network the data related to the connections between the organisms, the edge data, and then prepare to send the nodes and edges to Cytoscape using the function cyPlot()
.
genenet.edges <- data.frame(igraph::as_edgelist(graph_vir_prok)) names(genenet.edges) <- c("name.1", "name.2") genenet.edges$Weight <- igraph::edge_attr(graph_vir_prok)[[1]] genenet.edges$name.1 <- as.character(genenet.edges$name.1) genenet.edges$name.2 <- as.character(genenet.edges$name.2) genenet.nodes$name <- as.character(genenet.nodes$name) ug <- cyPlot(genenet.nodes,genenet.edges)
Now we will send the network from R to Cytoscape.
cw <- CytoscapeWindow("Tara oceans", graph = ug, overwriteWindow = TRUE)
displayGraph(cw) layoutNetwork(cw) fitContent(cw)
saveImage(cw, "co-occur0", "png", h = 1200) knitr::include_graphics("./co-occur0.png")
We would like to get an overview of the different phylum of bacteria that are in the network. One way is to colour the different nodes based on their phylum classification. The package Rcolorbrewer will be used to generate a set of good colours for the nodes.
families_to_colour <- unique(genenet.nodes$prok_tax_phylum) families_to_colour <- families_to_colour[!families_to_colour %in% "not_class"] node.colour <- RColorBrewer::brewer.pal(length(families_to_colour), "Set3")
Use the colours from Rcolorbrewer to colour the nodes in Cytoscape.
setNodeColorRule(cw, "prok_tax_phylum", families_to_colour, node.colour, "lookup", default.color = "#ffffff")
displayGraph(cw) layoutNetwork(cw) fitContent(cw)
saveImage(cw, "co-occur0_1", "png", h = 1200) knitr::include_graphics("./co-occur0_1.png")
Next we would like to change the shape of the node to reflect whether the nodes are viral or prokaryotic in origin. In this dataset all of the viral node names start with "ph_", thus we can set the viral nodes to be diamond-shaped by looking for all the nodes that start with "ph" in the network.
shapes_for_nodes <- c("DIAMOND") phage_names <- grep("ph_", genenet.nodes$name, value = TRUE) setNodeShapeRule(cw, "label", phage_names, shapes_for_nodes)
displayGraph(cw) fitContent(cw)
saveImage(cw, "co-occur1", "png", h = 1200) knitr::include_graphics("./co-occur1.png")
The classification of the viral data was done in a very conservative manner so not many of the viral nodes were identified. However, if we do want to add some of this information to our visualization we can colour the edges of the viral nodes by family. The main families that were identified in this dataset are the Podoviridae, the Siphoviridae and the Myoviridae (for more info see NCBI Podoviridae, NCBI Myoviridae, and NCBI Siphoviridae)
setDefaultNodeBorderWidth(cw, 5) families_to_colour <- c(" Podoviridae", " Siphoviridae", " Myoviridae") node.colour <- RColorBrewer::brewer.pal(length(families_to_colour), "Dark2") setNodeBorderColorRule(cw, "Tax_subfamily", families_to_colour, node.colour, "lookup", default.color = "#000000")
displayGraph(cw) fitContent(cw)
saveImage(cw, "co-occur2", "png", h = 1200) knitr::include_graphics("./co-occur2.png")
After doing all of this colouring to the network we would like to layout the network in a way that allows us to more easily see which nodes are connected without overlap. To do the layout we will use the Cytoscape plugin Allegro.
When using RCy3 to drive Cytoscape, if we are not sure what the current values are for a layout or we are not sure what kinds of values are accepted for the different parameters of our layout, we can investigate using the RCy3 functions getLayoutPropertyNames()
and then getLayoutPropertyValue()
.
getLayoutNames(cw) getLayoutPropertyNames(cw, layout.name = "allegro-spring-electric") getLayoutPropertyValue(cw, "allegro-spring-electric","gravity") getLayoutPropertyValue(cw, "allegro-spring-electric","maxIterations") getLayoutPropertyValue(cw, "allegro-spring-electric","noOverlapIterations")
Once we decide on the properties we want, we can go ahead and set them like this:
setLayoutProperties(cw, layout.name = "allegro-spring-electric", list(gravity = 100, scale = 6)) layoutNetwork(cw, layout.name = "allegro-spring-electric") fitContent(cw)
saveImage(cw, "co-occur3", "png", h = 1200) knitr::include_graphics("./co-occur3.png")
One thing that might be interesting to visualize is nodes that are connected to many different nodes and nodes that are connected to few other nodes. The number of other nodes to which one node is connected is called degree. We can use a gradient of size to quickly visualize nodes that have high degree.
## initiate a new node attribute ug2 <- initNodeAttribute(graph = ug, "degree", "numeric", 0.0) ## degree from graph package for undirected graphs not working well, ## so instead using igraph to calculate this from the original graph nodeData(ug2, nodes(ug2), "degree") <- igraph::degree(graph_vir_prok) cw2 <- CytoscapeWindow("Tara oceans with degree", graph = ug2, overwriteWindow = TRUE)
displayGraph(cw2) layoutNetwork(cw2)
degree_control_points <- c(min(igraph::degree(graph_vir_prok)), mean(igraph::degree(graph_vir_prok)), max(igraph::degree(graph_vir_prok))) node_sizes <- c(20, 20, 80, 100, 110) # number of control points in interpolation mode, # the first and the last are for sizes "below" and "above" the attribute seen. setNodeSizeRule(cw2, "degree", degree_control_points, node_sizes, mode = "interpolate") layoutNetwork(cw2, "force-directed")
fitContent(cw2) Sys.sleep(10) # to make sure content is fit before taking an image fitContent(cw2) saveImage(cw2, "co-occur_degree", "png", h=1200)
knitr::include_graphics("./co-occur_degree.png")
The visualization displays several different areas where there are highly connected nodes that are in the same bacterial phylum. We will select one of these nodes, all of the nodes connected to this node, its first neighbours, and then the nodes connected to the first neighbours. One node that is in a group of highly connected nodes is the cyanobacterial node "GQ377772". We will select it and its first and second neighbours and then make a new network from these nodes and their connections.
selectNodes(cw2, "GQ377772") # selects specific nodes getSelectedNodes(cw2) selectFirstNeighborsOfSelectedNodes(cw2) getSelectedNodes(cw2)
Now select the second neighbours of node "GQ377772".
selectFirstNeighborsOfSelectedNodes(cw2) getSelectedNodes(cw2)
Create subnetwork from these nodes and their edges.
newnet <- createWindowFromSelection(cw2, "subnet", "TRUE") layoutNetwork(newnet, "force-directed")
fitContent(newnet) Sys.sleep(10) saveImage(newnet, "co-occur_subnet", "png", h = 1200)
knitr::include_graphics("./co-occur_subnet.png")
This has been a very basic introduction to exploring co-occurence networks in Cytoscape using RCy3.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.