######################################################################
# Load libraries, some of these are new and you will need to install
######################################################################
library(vegan)
library(BiodiversityR)
library(labdsv)
library(gclus)

# Site by species matrix of sample data
# Data is stored in the variable "community"
community<-read.csv("sample.csv", header=TRUE, row.names="site")
community<-as.matrix(community)
sample_factors<-read.csv("factors.csv",header=T)

# Rank abundance curves
rank_data<-rankabundance(community)
rankabunplot(rank_data)

# Compare rank abundance among creeks or years
# THIS FUNCTION WILL WAIT FOR YOU TO CLICK ON THE GRAPH TO PLACE THE LEGEND
rank_data<-rankabuncomp(community,y=sample_factors,factor="creek")

# Compare rank abundance among creeks or years
# THIS FUNCTION WILL WAIT FOR YOU TO CLICK ON THE GRAPH TO PLACE THE LEGEND
rank_data<-rankabuncomp(community,y=sample_factors,factor="year")

# Disimilarity indices
braycurt<-vegdist(community,method="bray")
# note the structure of the resulting matix
braycurt

euclidean<-vegdist(community,method="euclidean")
euclidean

# Cody(1975) measure of beta diversity
cody_1975<-designdist(community,"(a+c)/2", abcd=TRUE)
cody_1975

# use of dsvdis function in labdsv package
steinhaus<-dsvdis(community,index="steinhaus")
steinhaus

# genetic distances
library(adegenet)
#load a microsat dataset that comes with the package
data(microsatt) 
# 9 loci, 18 populations, allele frequency
microsatt
# table of allele frequencies (18 populations x 112 allels)
microsatt$tab
dim(microsatt$tab)
# euclidean distance
msat<-euclidean<-dist(microsatt$tab)
# Nei's distance
nei<-dist.genpop(as.genpop(microsatt$tab), method=1)
plot(msat,nei)

# heatmap of genetic distances
heatmap(as.matrix(nei),Rowv=T,Colv=T)

# heatmap of raw data
heatmap(as.matrix(community),Rowv=NA,Colv=NA)

# heatmap of similarity matrix
# note that this is a triangular matrix that is made square by copying
# the top values
heatmap(x=as.matrix(braycurt),Rowv=NA,Colv=NA)

# reorder values so that similar samples next to each other and dendrograms are plotted in margins
heatmap(x=as.matrix(braycurt), Rowv=T,Colv=T)

# Colors are not great...try RColorBrewer package
library(RColorBrewer)

display.brewer.all(type="seq")

heatmap(x=as.matrix(braycurt),col = brewer.pal(9,"PuRd"), Rowv=T,Colv=T)