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