Description of vignette

Tara Oceans

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].

How to examine organisms that occur together: co-occurence networks

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.

What can we find out by creating co-occurence networks?

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).

What kind of data are used in this vignette?

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

Set up Cytoscape and R connection

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].

Requirements

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)

Read in data

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))

Read in taxonomic classification

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")

Add the taxonomic classifications to the network and then send network to Cytoscape

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)

Send network to Cytoscape using RCy3

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")

Colour network by prokaryotic phylum

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")

Set node shape to reflect virus or prokaryote

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")

Colour edges of phage nodes

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")

Do layout to minimize overlap of nodes.

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")

Look at network properties

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)

Size by degree

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")

Select an interesting node and make a subnetwork from it

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")

Conclusion

This has been a very basic introduction to exploring co-occurence networks in Cytoscape using RCy3.

References



tmuetze/Bioconductor_RCy3_the_new_RCytoscape documentation built on May 31, 2019, 4:39 p.m.