In spholmes/GraphsTutorial: Practical Tutorial on Graph Manipulation

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

Testing a graph for association with a factor variable

For univariate data in two groups, we can test the differences by looking at the number of runs within each group.

Seeing the number of runs in a one-dimensional, two-sample, nonparametric Wald-Wolfowitz test can indicate whether the two groups have the same distributions.

library(ggplot2)
dfbr=data.frame(measure=c(rnorm(15,0.9),rnorm(15,1.8)),
  group=as.factor(c(rep("men",15),rep("women",15))))
ggplot(dfbr,aes(x=measure,group=group,y=0)) + ylim(-0.25,+0.25) +
  geom_point(aes(col=group,x=measure,y=0,shape=group),size=5,alpha=0.6)+
  scale_color_manual(values=c("blue","red"))+
  theme_bw() + geom_hline(yintercept = 0) +
  theme(panel.border = element_blank(),
  axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        axis.title.x=element_blank(),
        panel.grid.major.y = element_blank() ,
        panel.grid.minor.y = element_blank() )+ coord_fixed()

Testing the relation between a graph and vertex factors

The samples have factor covariates

The nodes in the graph correspond to the samples and are associated to factors that specify the litter and individual ids of the mice from which the samples were collected.

library("phyloseq")
library("phyloseqGraphTest")
library("igraph")
data(micephyloseq)
#micephyloseq  = readRDS(url("https://web.stanford.edu/class/bios221/data/ps1.rds"))
micephyloseq
sampledata = data.frame( sample_data(micephyloseq))
d1 = as.matrix(phyloseq::distance(micephyloseq, method="jaccard"))
gr = graph.adjacency(d1,  mode = "undirected", weighted = TRUE)
net = igraph::mst(gr)
V(net)$id = sampledata[names(V(net)), "host_subject_id"]
V(net)$litter = sampledata[names(V(net)), "family_relationship"]

The minimum spanning tree based on Jaccard dissimilarity and annotated with the mice ID and litter factors.

library(ggnetwork)
gnet=ggnetwork(net)
ggplot(gnet, aes(x = x, y = y, xend = xend, yend = yend))+
  geom_edges(color = "darkgray") +
  geom_nodes(aes(color = id, shape = litter)) + theme_blank()+
  theme(legend.position="bottom")

library("phyloseqGraphTest")
gt = graph_perm_test(micephyloseq, "host_subject_id", distance="jaccard",
                     type="mst",  nperm=1000)
plot_test_network(gt)

gt
plot_permutations(gt)
gt$pval

The permutation histogram of the number of pure edges in the network obtained from the minimal spanning tree with Jaccard similarity.

net = make_network(micephyloseq, max.dist = 0.35)
sampledata = data.frame(sample_data(micephyloseq))
V(net)$id = sampledata[names(V(net)), "host_subject_id"]
V(net)$litter = sampledata[names(V(net)), "family_relationship"]
netg = ggnetwork(net)

ggplot(netg, aes(x = x, y = y, xend = xend, yend = yend)) +
  geom_edges(color = "darkgray") +
  geom_nodes(aes(color = id, shape = litter)) + theme_blank()+
    theme(plot.margin = unit(c(0, 5, 2, 0), "cm"))+
    theme(legend.position = c(1.4, 0.3),legend.background = element_blank(),
          legend.margin=margin(0, 3, 0, 0, "cm"))+
         guides(color=guide_legend(ncol=2))

gt = graph_perm_test(micephyloseq, "family_relationship",
        grouping = "host_subject_id",
        distance = "jaccard", type = "mst", nperm= 1000)
gt$pval

plot_permutations(gt)

gtnn1 = graph_perm_test(micephyloseq, "family_relationship",
                      grouping = "host_subject_id",
                      distance = "jaccard", type = "knn", knn = 1)
gtnn1$pval