Last active
January 16, 2017 03:57
-
-
Save wilkox/a3faa160a4c52a8cc756a6b1d27d77a4 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # Libraries | |
| library(readr) | |
| library(SpiecEasi) | |
| library(phyloseq) | |
| library(tibble) | |
| library(igraph) | |
| library(dplyr) | |
| library(stringr) | |
| library(tidyr) | |
| library(ggplot2) | |
| library(Matrix) | |
| library(orca) | |
| library(ggnetwork) | |
| library(forcats) | |
| library(intergraph) | |
| library(network) | |
| # 1 Prepare phyloseq OTU table | |
| OTUTableADMa3zW <- "ADMa3z_Winter.txt" %>% | |
| read_tsv %>% | |
| select(OTU, Sample, Count) %>% | |
| mutate(OTU = word(OTU, 1, sep = ";")) %>% | |
| # 2 Select only most abundant OTUs | |
| # Delete the marked lines below to include all OTUs | |
| group_by(OTU) %>% # <- Delete | |
| mutate(TotalCount = sum(Count)) %>% # <- Delete | |
| ungroup %>% # <- Delete | |
| filter(TotalCount >= 100) %>% # <- Delete | |
| select(-TotalCount) %>% # <- Delete | |
| spread(OTU, Count, fill = 0) %>% | |
| as.data.frame %>% | |
| remove_rownames %>% | |
| column_to_rownames("Sample") %>% | |
| as.matrix %>% | |
| otu_table(taxa_are_rows = FALSE) | |
| # 3 Prepare phyloseq taxa table | |
| TaxaTableADMa3zW <- "ADMa3z_Winter.txt" %>% | |
| read_tsv %>% | |
| select(OTU, Domain = Kingdom, Phylum, Class, Order, Family, Genus, Species) %>% | |
| mutate(OTU = word(OTU, 1, sep = ";")) %>% | |
| unique %>% | |
| as.data.frame %>% | |
| remove_rownames %>% | |
| column_to_rownames("OTU") %>% | |
| as.matrix %>% | |
| tax_table | |
| # 4 Combine OTU and taxonomy into single phyloseq object | |
| PhyseqADMa3zW <- phyloseq(OTUTableADMa3zW,TaxaTableADMa3zW) | |
| # 4.1 Extract OTU names for later use | |
| OTUNames <- PhyseqADMa3zW %>% otu_table %>% colnames | |
| # 5 Perform SPIEC_EASI | |
| se.ADMa3zW <- PhyseqADMa3zW %>% | |
| spiec.easi( | |
| method = 'mb', | |
| lambda.min.ratio = 1e-2, | |
| nlambda = 20, | |
| icov.select.params = list(rep.num = 20, ncores = 2) | |
| ) | |
| # 6 Convert to igraph and network objects | |
| ig.ADMa3zW <- se.ADMa3zW %>% | |
| .$refit %>% | |
| adj2igraph(vertex.attr = list(name = taxa_names(PhyseqADMa3zW))) | |
| Network <- asNetwork(ig.ADMa3zW) | |
| # 7 Determine correlation weights and direction (positive or negative | |
| # correlation), and save | |
| elist.weighted.ADMa3zW <- se.ADMa3zW %>% | |
| getOptBeta %>% | |
| symBeta(mode = 'maxabs') %>% | |
| summary %>% | |
| # Assuming OTU numbers in i and j correlate to the order of OTUs in | |
| # se.ADMa3zw, restore OTU names | |
| left_join(data.frame(OTU1 = OTUNames, i = 1:length(OTUNames))) %>% | |
| left_join(data.frame(OTU2 = OTUNames, j = 1:length(OTUNames))) | |
| elist.weighted.ADMa3zW %>% | |
| write_tsv("elist_ADMa3zW_weighted.tsv") | |
| # 8 Identify top seven genera by number of edges (i.e. sum of degrees across all | |
| # OTUs in that genus) | |
| TopDegreeGenera <- ig.ADMa3zW %>% | |
| igraph::degree() %>% | |
| data.frame(OTU = names(.), Degree = .) %>% | |
| remove_rownames %>% | |
| as.tbl %>% | |
| left_join( | |
| "ADMa3z_Winter.txt" %>% | |
| read_tsv %>% | |
| select(OTU, Domain = Kingdom, Phylum, Class, Order, Family, Genus, Species) %>% | |
| unique %>% | |
| mutate(OTU = word(OTU, 1, sep = ";")) | |
| ) %>% | |
| group_by(Genus) %>% | |
| summarise(SumDegree = sum(Degree)) %>% | |
| arrange(desc(SumDegree)) %>% | |
| na.omit %>% | |
| .[1:9, ] %>% | |
| .$Genus | |
| # Prepare vertex data to add to network | |
| VertexData <- | |
| # Get OTU names for each vertex, as they are ordered in the network | |
| Network %>% | |
| network::get.vertex.attribute("vertex.names") %>% | |
| data.frame(OTU = .) %>% | |
| as.tbl %>% | |
| # Add genus | |
| left_join( | |
| "ADMa3z_Winter.txt" %>% | |
| read_tsv %>% | |
| select(OTU, Genus) %>% | |
| unique | |
| ) %>% | |
| # Collapse genera to top by number of edges | |
| mutate(Genus = ifelse(Genus %in% TopDegreeGenera, Genus, "Low-degree genera")) %>% | |
| mutate(Genus = as.character(Genus)) | |
| # Add vertex data to network | |
| Network %v% "Genus" <- VertexData$Genus | |
| # Prepare weights to add to network: weights < 0.15 set to 0 | |
| elist.weighted.ADMa3zW <- elist.weighted.ADMa3zW %>% | |
| mutate(x = ifelse(x < 0.15, NA, x)) | |
| Network <- Network %>% | |
| set.edge.attribute("Correlation", elist.weighted.ADMa3zW$x) | |
| # Flatten network to dataframe with default layout | |
| NetworkDF <- ggnetwork(Network) | |
| # 10 Draw plot | |
| Plot <- NetworkDF %>% | |
| ggplot(aes(x = x, y = y, xend = xend, yend = yend)) + | |
| geom_edges(aes(alpha = Correlation)) + | |
| geom_nodes(aes(colour = Genus)) + | |
| theme_blank() | |
| ggsave("spiec_easi_genus_ADMa3z_winter.png", plot = Plot, width = 4.74, height = 5.33, units = "in") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment