Cacher is the code snippet organizer for pro developers

plot_map_bc <- function(physeq = physeq){

library(vegan)
library(tidyverse)
library(broom)
library(Hmisc)
library(viridis)

# calculate bray-curtis

physeq_data <-  sample_data(physeq_hel)
physeq_data <- as(physeq_data, "data.frame")

physeq_data <- physeq_data %>% 
  filter(Metagenomes %in% sample_names(physeq)) %>% 
  select(Metagenomes, Longitude, Latitude, Depth_Id)

srf <- physeq_data %>% filter(Depth_Id == "SRF")

physeq_hel <- prune_samples(srf$Metagenomes, physeq_hel)

physeq_bc <- vegdist(otu_table(physeq), "bray")
physeq_bc_graph <- tidy(physeq_bc) %>%
  mutate(weight = 1 - distance) %>%
  graph_from_data_frame() %>%
  simplify()

binary_search_filter_1 <- function(g, low, high) {
  x<-g
  if ( high < low ) {
    
    return(NULL)
    
  }else {
    mid <- floor((low + high) / 2)
    x<-igraph::delete.edges(g, which(E(g)$weight <= weights[mid])) 
    cat("low:", low, " mid:", mid, " high:", high, " mid_weight:", weights[mid], " components:", igraph::count_components(x), " edges:", ecount(x), "\n")
    
    if ((mid - low) == 0){
      x<-igraph::delete.edges(g, which(E(g)$weight <= weights[mid-1])) 
      cat("Final: low:", low - 1, " mid:", mid - 1, " high:", high - 1, " mid_weight:", weights[mid - 1], " components:",igraph::count_components(x), " edges:", ecount(x), "\n")
      return(list(weight=mid - 1, graph=x))
      exit
    }
    
    if (((high - low) == 0)){
      x<-igraph::delete.edges(g, which(E(g)$weight <= weights[mid+1])) 
      cat("Final: low:", low, " mid:", mid, " high:", high, " mid_weight:", weights[mid], " components:",igraph::count_components(x), " edges:", ecount(x), "\n")
      return(list(weight=mid , graph=x))
      exit
    }
    
    
    if (!(igraph::is.connected(x))){
      binary_search_filter_1(g, low, mid - 1)
    }else if (igraph::is.connected(x)){
      binary_search_filter_1(g, mid + 1, high)
    }
  }
}

weights <- unique(sort(E(physeq_bc_graph)$weight, decreasing = FALSE))
physeq_bc_graph_filt <- binary_search_filter_1(g = physeq_bc_graph, low = 1, high = length(weights))$graph

layoutCoordinates <- physeq_data
adjacencyList <- igraph::as_data_frame(physeq_bc_graph_filt) 

edgeMaker <- function(whichRow, len = 100, curved = TRUE){
  # Get coordinates for each metagenome
  fromC <- as.numeric(subset(layoutCoordinates, Metagenomes == as.character(adjacencyList[whichRow,]$from))[2:3])  # Origin
  toC <- as.numeric(subset(layoutCoordinates, Metagenomes == as.character(adjacencyList[whichRow,]$to))[2:3])  # Terminus
  
  # Add curve:
  graphCenter <- colMeans(layoutCoordinates[2:3])  # Center of the overall graph
  bezierMid <- c(fromC[1], toC[2])  # A midpoint, for bended edges
  distance1 <- sum((graphCenter - bezierMid)^2)
  if(distance1 < sum((graphCenter - c(toC[1], fromC[2]))^2)){
    bezierMid <- c(toC[1], fromC[2])
  }  # To select the best Bezier midpoint
  bezierMid <- (fromC + toC + bezierMid) / 3  # Moderate the Bezier midpoint
  if(curved == FALSE){bezierMid <- (fromC + toC) / 2}  # Remove the curve
  
  edge <- data.frame(bezier(c(fromC[1], bezierMid[1], toC[1]),  # Generate
                            c(fromC[2], bezierMid[2], toC[2]),  # X & y
                            evaluation = len))  # Bezier path coordinates
  edge$Sequence <- 1:len  # For size and colour weighting in plot
  edge$Group <- paste(adjacencyList[whichRow,]$from,adjacencyList[whichRow,]$to, sep = "#")
  edge$weight <- adjacencyList[whichRow,]$weight
  return(edge)
}

# Generate a (curved) edge path for each pair of connected nodes
allEdges <- lapply(1:nrow(adjacencyList), edgeMaker, len = 500, curved = TRUE)
allEdges <- do.call(rbind, allEdges)  # a fine-grained path ^, with bend ^

# Plot 
mapWorld <- borders("world", colour="#0a0c2a", fill="#0a0c2a") # create a layer of borders

zp1 <- ggplot(allEdges) + mapWorld  # Pretty simple plot code
zp1 <- zp1 + geom_path(aes(x = x, y = y, group = Group,  # Edges with gradient
                           colour = weight, size = -Sequence))  # and taper
zp1 <- zp1 + geom_point(data = data.frame(physeq_data),  # Add nodes
                        aes(x = Longitude, y = Latitude), size = 2, pch = 21,
                        colour = "black", fill = "gray") # Customize gradient v
zp1 <- zp1 + scale_colour_gradientn(colours = rev(viridis(100, option = "C", direction = -1)), labels = scales::percent, name="Similarity")
zp1 <- zp1 + scale_size(range = c(1/10, 1), guide = "none")  # Customize taper
zp1 <- zp1 + theme_bw() + coord_equal(ratio=1) + theme(legend.position="top") +
  xlab("") + ylab("") + scale_x_continuous(breaks=c(-100,0,100))# Clean up plot
zp1
}