In spholmes/GraphsTutorial: Practical Tutorial on Graph Manipulation

library(learnr)
knitr::opts_chunk$set(echo = TRUE)

Co-occurrence graphs

We are going to look at microbial abundances of different strains/species of bacteria in biological samples.

These samples come from a longitudinal study of mice from several litters.

Bacterial community example using a `phyloseq` object

We load the phyloseq package first and read the phyloseq object that integrates all the data.

library("phyloseq")
#Read an external file
micephyloseq  = readRDS(url("https://web.stanford.edu/class/bios221/data/ps1.rds"))
#Internal file from the package is also available 
# data(micephyloseq)
micephyloseq

We could like to connect the samples that have a lot of bacteria in common. The best way to measure similarity between species abundances is the "Jaccard distance".

We start by computing the Jaccard distance between samples and set a threshold for the distance to correspond to an edge in the "co-occurrence" graph.

Exercise: Add a line of code that makes the shape give the litter information:

data(micephyloseq)
micephyloseq
ig<-make_network(micephyloseq, max.dist=0.3)
plot_network(ig,micephyloseq,color="host_subject_id")

ig<-make_network(micephyloseq, max.dist=0.3)
plot_network(ig,micephyloseq,color="host_subject_id",shape="family_relationship")

We can replicate what the function is doing by hand, using the tidygraph package, you have to add a few lines to remove the isolates.

We start by loading the libraries

library("tidygraph")
library("igraph")
library("ggraph")

d1 <- as.matrix(phyloseq::distance(micephyloseq, method="jaccard"))
d1 <- d1 +diag(ncol(d1))
threshM <- (d1<0.35)
gr <- graph_from_adjacency_matrix(threshM,  mode = "lower")
sampledata = data.frame( sample_data(micephyloseq))
V(gr)$id = sampledata[names(V(gr)), "host_subject_id"]
V(gr)$litter = sampledata[names(V(gr)), "family_relationship"]
tigr<- gr%>% as_tbl_graph()
tigr
ggraph(tigr,layout="fr")+ geom_edge_link() +
  geom_node_point(aes(color=id,shape=litter))+
  theme_graph()

tigr %>%
  activate(nodes) %>% 
  filter(!node_is_isolated()) %>%
  ggraph(,layout="fr") + geom_edge_link() +
  geom_node_point(aes(color=id,shape=litter))+
  theme_graph()