+
+
+
+
+
+
+Deliverable 1
+Nicola Kriefall
+2/26/2020
+ +
+
+
+
+
+
+Deliverable 2 & 3
+
+
+Overview
+For my project, I am analyzing bacterial microbiomes hosted by corals in the Florida Keys. Tissue samples of corals were collected at 6 sites throughout the Florida Keys in 2015, 2017 (just after Hurricane Irma), and 2018. I will be determining whether we can see an impact of Hurricane Irma in disrupted microbial networks at the Irma timepoint. Sequencing data collected from these samples were analyzed following Callahan et al., 2016 to identify which types of bacteria are present in the coral samples. The scripts I use to create the input data (.csv files) is available here
+Between deliverable 1 and now, I expanded my investigation from only 1 of 6 sites to all of the sites combined, and added the 2018 timepoint. I also added a few network metrics.
+
+
+Workflow
+Begin by installing necessary packages.
+# Install Required packages
+#install.packages("igraph")
+library("igraph")##
+## Attaching package: 'igraph'## The following objects are masked from 'package:stats':
+##
+## decompose, spectrum## The following object is masked from 'package:base':
+##
+## union#install.packages("qgraph")
+library("qgraph")## Registered S3 methods overwritten by 'huge':
+## method from
+## plot.sim BDgraph
+## print.sim BDgraph#install.packages("MCL")
+library("MCL")
+#install.packages("vegan")
+library("vegan")## Loading required package: permute##
+## Attaching package: 'permute'## The following object is masked from 'package:igraph':
+##
+## permute## Loading required package: lattice## This is vegan 2.5-6##
+## Attaching package: 'vegan'## The following object is masked from 'package:igraph':
+##
+## diversity# Install SpiecEasi package
+#install.packages("devtools")
+library("devtools") ## Loading required package: usethis##
+## Attaching package: 'devtools'## The following object is masked from 'package:permute':
+##
+## check#install_github("zdk123/SpiecEasi")
+library("SpiecEasi")##
+## Attaching package: 'SpiecEasi'## The following object is masked from 'package:igraph':
+##
+## make_graphlibrary("ggplot2")
+library("ggpubr")## Loading required package: magrittr
+
+Reading in & renaming data
+Data is organized with coral samples in rows, and the types of bacteria the corals host in columns. The data are counts of the number of sequence reads that matched each type of bacteria.
+data_1516 <- read.csv("seq.trim.sid.1516.no0 copy.csv",row.names=1) #from 2015&2016
+data_17 <- read.csv("seq.trim.sid.17.no0 copy.csv",row.names=1) #from 2017
+data_18 <- read.csv("seq.trim.sid.18.no0 copy.csv",row.names=1) #from 2018
+
+#removing samples with less than an average of 2 reads per OTU - recommended by authors
+data_1516_2 <- data_1516[,-which(colMeans(data_1516) <= 2)]
+data_17_2 <- data_17[,-which(colMeans(data_17) <= 2)]
+data_18_2 <- data_18[,-which(colMeans(data_18) <= 2)]
+
+Building the network
+This is a network of correlations of the abundances of the types of bacteria across samples. Following instructions here.
+# SparCC network
+sparcc.matrix_1516 <- sparcc(data_1516_2)
+sparcc.matrix_17 <- sparcc(data_17_2)
+sparcc.matrix_18 <- sparcc(data_18_2)
+
+sparcc.cutoff <- 0.3
+sparcc.adj_1516 <- ifelse(abs(sparcc.matrix_1516$Cor) >= sparcc.cutoff, 1, 0)
+sparcc.adj_17 <- ifelse(abs(sparcc.matrix_17$Cor) >= sparcc.cutoff, 1, 0)
+sparcc.adj_18 <- ifelse(abs(sparcc.matrix_18$Cor) >= sparcc.cutoff, 1, 0)
+
+# Add OTU names to rows and columns
+rownames(sparcc.adj_1516) <- colnames(data_1516_2)
+colnames(sparcc.adj_1516) <- colnames(data_1516_2)
+
+rownames(sparcc.adj_17) <- colnames(data_17_2)
+colnames(sparcc.adj_17) <- colnames(data_17_2)
+
+rownames(sparcc.adj_18) <- colnames(data_18_2)
+colnames(sparcc.adj_18) <- colnames(data_18_2)
+
+# Build network from adjacency
+sparcc.net_1516 <- graph.adjacency(sparcc.adj_1516,
+ mode = "undirected",
+ diag = FALSE)
+sparcc.net_17 <- graph.adjacency(sparcc.adj_17,
+ mode = "undirected",
+ diag = FALSE)
+sparcc.net_18 <- graph.adjacency(sparcc.adj_18,
+ mode = "undirected",
+ diag = FALSE)
+
+# Use sparcc.net for the rest of the method
+net_1516 <- sparcc.net_1516
+net_17 <- sparcc.net_17
+net_18 <- sparcc.net_18
+# Hub detection
+net_1516.cn <- closeness(net_1516)## Warning in closeness(net_1516): At centrality.c:2784 :closeness centrality is
+## not well-defined for disconnected graphsnet_1516.bn <- betweenness(net_1516)
+net_1516.pr <- page_rank(net_1516)$vector
+net_1516.hs <- hub_score(net_1516)$vector
+net_1516.hs.sort <- sort(net_1516.hs, decreasing = TRUE)
+net_1516.hs.top5 <- head(net_1516.hs.sort, n = 5)
+
+net_17.cn <- closeness(net_17)## Warning in closeness(net_17): At centrality.c:2784 :closeness centrality is not
+## well-defined for disconnected graphsnet_17.bn <- betweenness(net_17)
+net_17.pr <- page_rank(net_17)$vector
+net_17.hs <- hub_score(net_17)$vector
+net_17.hs.sort <- sort(net_17.hs, decreasing = TRUE)
+net_17.hs.top5 <- head(net_17.hs.sort, n = 5)
+
+net_18.cn <- closeness(net_18)## Warning in closeness(net_18): At centrality.c:2784 :closeness centrality is not
+## well-defined for disconnected graphsnet_18.bn <- betweenness(net_18)
+net_18.pr <- page_rank(net_18)$vector
+net_18.hs <- hub_score(net_18)$vector
+net_18.hs.sort <- sort(net_18.hs, decreasing = TRUE)
+net_18.hs.top5 <- head(net_18.hs.sort, n = 5)
+
+Plotting the network
+# Function 2: Plot network with node size scaled to hubbiness
+plot.net <- function(net, scores, outfile, title) {
+# Convert node label from names to numerical IDs.
+ features <- V(net)$name
+ col_ids <- seq(1, length(features))
+ V(net)$name <- col_ids
+ node.names <- features[V(net)] # Nodes’ color.
+ V(net)$color <- "white"
+ # Define output image file.
+ outfile <- paste(outfile, "jpg", sep=".") # Image properties.
+ jpeg(outfile, width = 4800, height = 9200, res = 300, quality = 100)
+ par(oma = c(4, 1, 1, 1))
+ # Main plot function.
+ plot(net, vertex.size = (scores*5)+4, vertex.label.cex = 1)
+ title(title, cex.main = 4)
+ # Plot legend containing OTU names.
+ #labels = paste(as.character(V(net)), node.names, sep = ") ")
+ #legend("bottom", legend = labels, xpd = TRUE, ncol = 5, cex = 1.2)
+ dev.off()
+}
+
+plot.net(net_1516, net_1516.hs, outfile = "network_1516", title = "Network - 2015 & 2016")## quartz_off_screen
+## 2plot.net(net_17, net_17.hs, outfile = "network_17", title = "Network - 2017")## quartz_off_screen
+## 2plot.net(net_18, net_18.hs, outfile = "network_18", title = "Network - 2018")## quartz_off_screen
+## 2
+
+Looking at results
+Opening up the .jpg files created in the step above (network_15.jpg & network_17.jpg), it looks like there are a lot less connections in the microbiome in 2017 (Hurricane Irma timepoint), which is exactly what I was hypothesizing. I’m going to need to figure out how to quantify that though.
+
+
+Extracting some network metrics
+
+
+Eigenvector centrality
+eigen.1516 <- eigen_centrality(net_1516)
+eigen.1516.tosort <- data.frame(eigen.1516[["vector"]])
+eigen.1516.tosort$vector <- eigen.1516.tosort$eigen.1516...vector...
+eigen.1516.sorted <- eigen.1516.tosort[order(-eigen.1516.tosort$vector),]
+eigen.1516.sorted$sqs <- rownames(eigen.1516.sorted)
+eigen.1516.top10 <- eigen.1516.sorted[1:10,]
+sqs_order <- eigen.1516.top10$sqs
+eigen.1516.top10$sqs <- factor(eigen.1516.top10$sqs,levels=sqs_order)
+gg.eig.1516 <- ggplot(eigen.1516.top10,aes(x=sqs,y=vector))+
+ geom_bar(stat="identity")+
+ theme(axis.text.x=element_text(angle=90,hjust=1))+
+ ylab("Eigenvector centrality")+
+ xlab("Bacterial sequence")+
+ ggtitle("2015 & 2016 (pre-Irma)")
+
+eigen.17 <- eigen_centrality(net_17)
+eigen.17.tosort <- data.frame(eigen.17[["vector"]])
+eigen.17.tosort$vector <- eigen.17.tosort$eigen.17...vector...
+eigen.17.sorted <- eigen.17.tosort[order(-eigen.17.tosort$vector),]
+eigen.17.sorted$sqs <- rownames(eigen.17.sorted)
+eigen.17.top10 <- eigen.17.sorted[1:10,]
+sqs_order <- eigen.17.top10$sqs
+eigen.17.top10$sqs <- factor(eigen.17.top10$sqs,levels=sqs_order)
+gg.eig.17 <- ggplot(eigen.17.top10,aes(x=sqs,y=vector))+
+ geom_bar(stat="identity")+
+ theme(axis.text.x=element_text(angle=90,hjust=1))+
+ ylab("Eigenvector centrality")+
+ xlab("Bacterial sequence")+
+ ggtitle("2017 (Irma)")
+
+eigen.18 <- eigen_centrality(net_18)
+eigen.18.tosort <- data.frame(eigen.18[["vector"]])
+eigen.18.tosort$vector <- eigen.18.tosort$eigen.18...vector...
+eigen.18.sorted <- eigen.18.tosort[order(-eigen.18.tosort$vector),]
+eigen.18.sorted$sqs <- rownames(eigen.18.sorted)
+eigen.18.top10 <- eigen.18.sorted[1:10,]
+sqs_order <- eigen.18.top10$sqs
+eigen.18.top10$sqs <- factor(eigen.18.top10$sqs,levels=sqs_order)
+gg.eig.18 <- ggplot(eigen.18.top10,aes(x=sqs,y=vector))+
+ geom_bar(stat="identity")+
+ theme(axis.text.x=element_text(angle=90,hjust=1))+
+ ylab("Eigenvector centrality")+
+ xlab("Bacterial sequence")+
+ ggtitle("2018 (post-Irma)")
+
+ggarrange(gg.eig.1516,gg.eig.17,gg.eig.18)
Two important takeaways: 1. It looks like no two bacteria (‘sqs’ on the x-axis) are the most central nodes across all three timepoints 2. Steepest drop in eigen centrality during the 2017 timepoint - during Hurricane Irma!
+
+
+Eigenvector centrality
+btwn.1516.tosort <- data.frame(betweenness(net_1516))
+btwn.1516.tosort$vector <- btwn.1516.tosort$betweenness.net_1516.
+btwn.1516.sorted <- btwn.1516.tosort[order(-btwn.1516.tosort$vector),]
+btwn.1516.sorted$sqs <- rownames(btwn.1516.sorted)
+btwn.1516.top10 <- btwn.1516.sorted[1:10,]
+sqs_order <- btwn.1516.top10$sqs
+btwn.1516.top10$sqs <- factor(btwn.1516.top10$sqs,levels=sqs_order)
+gg.btwn.1516 <- ggplot(btwn.1516.top10,aes(x=sqs,y=vector))+
+ geom_bar(stat="identity")+
+ theme(axis.text.x=element_text(angle=90,hjust=1))+
+ ylab("Betweenness centrality")+
+ xlab("Bacterial sequence")+
+ ggtitle("2015 & 2016 (pre-Irma)")
+
+btwn.17.tosort <- data.frame(betweenness(net_17))
+btwn.17.tosort$vector <- btwn.17.tosort$betweenness.net_17.
+btwn.17.sorted <- btwn.17.tosort[order(-btwn.17.tosort$vector),]
+btwn.17.sorted$sqs <- rownames(btwn.17.sorted)
+btwn.17.top10 <- btwn.17.sorted[1:10,]
+sqs_order <- btwn.17.top10$sqs
+btwn.17.top10$sqs <- factor(btwn.17.top10$sqs,levels=sqs_order)
+gg.btwn.17 <- ggplot(btwn.17.top10,aes(x=sqs,y=vector))+
+ geom_bar(stat="identity")+
+ theme(axis.text.x=element_text(angle=90,hjust=1))+
+ ylab("Betweenness centrality")+
+ xlab("Bacterial sequence")+
+ ggtitle("2017 (Irma)")
+
+btwn.18.tosort <- data.frame(betweenness(net_18))
+btwn.18.tosort$vector <- btwn.18.tosort$betweenness.net_18.
+btwn.18.sorted <- btwn.18.tosort[order(-btwn.18.tosort$vector),]
+btwn.18.sorted$sqs <- rownames(btwn.18.sorted)
+btwn.18.top10 <- btwn.18.sorted[1:10,]
+sqs_order <- btwn.18.top10$sqs
+btwn.18.top10$sqs <- factor(btwn.18.top10$sqs,levels=sqs_order)
+gg.btwn.18 <- ggplot(btwn.18.top10,aes(x=sqs,y=vector))+
+ geom_bar(stat="identity")+
+ theme(axis.text.x=element_text(angle=90,hjust=1))+
+ ylab("Betweenness centrality")+
+ xlab("Bacterial sequence")+
+ ggtitle("2018 (post-Irma)")
+
+ggarrange(gg.btwn.1516,gg.btwn.17,gg.btwn.18)
Similar pattern as in the eigenvector centrality plots
+Left to do: - More metric analysis & all the statistics!
+
