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