genomewalker
6/13/2017 - 12:51 PM

plot_distance_decay.r

plot_distance_decay <- funtion(physeq){
  
  physeq_data <-  sample_data(physeq)
  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 <- prune_samples(srf$Metagenomes, physeq)
  
  physeq_bc <- vegdist(otu_table(physeq), "bray")
  
  physeq_data <- physeq_data %>% 
    select(Metagenomes, Longitude, Latitude)
  
  rownames(physeq_data) <- physeq_data$Metagenomes
  physeq_data$Metagenomes <- NULL
  physeq_data <- as.matrix(physeq_data)
  physeq_data <- physeq_data[base::labels(physeq_bc),]
  geo_dm <- geosphere::distm(physeq_data, fun=geosphere::distHaversine)
  
  rownames(geo_dm) <- rownames(physeq_data)
  colnames(geo_dm) <- rownames(physeq_data)
  geo_dm_df <- broom::tidy(as.dist(geo_dm))
  
  physeq_bc <- broom::tidy(physeq_bc)
  
  dis <- geo_dm_df %>% left_join(physeq_bc, by = c("item1", "item2"))
  colnames(dis) <- c("loc1","loc2","value_geo","value_dis")
  
  ggplot(dis, aes(value_geo/1000, value_dis)) + 
    geom_point(alpha = 0.4, size = 1) + 
    geom_smooth(se = FALSE, method = "lm") +
    #scale_color_manual(values = c("#2095B8", "#EF2719", "#E0AE27")) +
    xlab("Spatial distance (Km)") + ylab("Community dissimilarity (Bray-Curtis)") + 
    theme_bw() +
    theme(plot.title = element_text(size=12),
          legend.position = "top") +
    ggtitle("Similarity-distance decay relationship")
}