setwd("/SET PATH TO DIRECTORY/Sandwave")
source("/SET PATH TO DIRECTORY/ErrorBar.R")
source("/SET PATH TO DIRECTORY/EigenBar.R")
require(ade4)
require(dunn.test)

#Data
data <- read.csv("ALL_MACROFAUNA_June2017.csv", h = T, check.names = F, as.is = F)
colnames(data)
tab <- data[,17:ncol(data)]
#Remove rare and inaccurate taxa
tab <- tab[,colnames(tab) %in% c("Annelida", "Bivalvia sp", "Brachyura", "Brachyura megalopa", "Mollusca",
                                 "Decapoda", "Gastropoda sp", "Phyllodocida sp",
                                 "Spio sp", "Tellinidae") == F]
#Taxonomic merging
colnames(tab)
tx <- colnames(tab)
tx[tx %in% c("Abra alba", "Abra prismatica")] <- "Abra sp"
tx[tx %in% c("Amphiura filiformis", "Ophiura albida", "Ophiura ophiura")] <- "Ophiuroidea"
tx[tx %in% c("Bathyporeia elegans")] <- "Bathyporeia sp"
tx[tx %in% c("Brachyura megalopa")] <- "Brachyura sp"
tx[tx %in% c("Diastylis bradyi", "Monopseudocuma gilsoni",
             "Pseudocuma longicornis", "Pseudocuma similis")] <- "Cumacea"
tx[tx %in% c("Chaetozone sp", "Tharyx sp")] <- "Cirratulidae"
tx[tx %in% c("Cumacea sp")] <- "Cumacea"
tx[tx %in% c("Dosinia exoleta")] <- "Dosinia sp"
tx[tx %in% c("Echinocardium cordatum")] <- "Echinoidea"
tx[tx %in% c("Echinodermata")] <- "Echinoidea"
tx[tx %in% c("Echinoida sp")] <- "Echinoidea"
tx[tx %in% c("Echinoidea sp")] <- "Echinoidea"
tx[tx %in% c("Ensis ensis", "Ensis leei", "Ensis silliqua", "Ensis sp", "Pharidae sp")] <- "Pharidae"
tx[tx %in% c("Magelona johnstoni", "Magelona mirabilis")] <- "Magelona sp"
tx[tx %in% c("Mysella bidentata")] <- "Kurtiella bidentata"
tx[tx %in% c("Nematoda sp")] <- "Nematoda"
tx[tx %in% c("Nemertea sp")] <- "Nemertea"
tx[tx %in% c("Nephtys caeca", "Nephtys cirrosa", "Nephtys hombergii",
             "Nephtys longosetosa","Nephtyidae sp")] <- "Nephtys sp"
tx[tx %in% c("Oligocheata")] <- "Haplotaxida"
tx[tx %in% c("Ophelia borealis", "Opheliinae sp")] <- "Ophelia sp"
tx[tx %in% c("Ophiuroidea sp")] <- "Ophiuroidea"
tx[tx %in% c("Pestarella tyrrhena")] <- "Gilyossius tyrrhenus"
tx[tx %in% c("Phaxas sp")] <- "Pharidae"
tx[tx %in% c("Phoronida sp")] <- "Phoronida"
tx[tx %in% c("Spio decorata")] <- "Spio decoratus"
tx[tx %in% c("Spisula subtruncata")] <- "Spisula sp"
tx[tx %in% c("Lanice conchilega")] <- "Terebellidae"
tx[tx %in% c("Venerupis corrugata")] <- "Venerupis sp"
#Sum taxa per higher taxonomic level
tab <- apply(tab, 1, function(x) tapply(x, factor(tx), sum))
tab <- data.frame(t(tab), check.names = F)
tab <- tab[,order(colnames(tab))]
colnames(tab)[colnames(tab) != "Gastrosaccus spinifer"] <- gsub(" sp", " sp.",
  colnames(tab)[colnames(tab) != "Gastrosaccus spinifer"], fixed = T)
colnames(tab)

#Samples x taxa matrices
w <- factor(data$Name)[data$Type2 == "counts" & data$Type1 != "Fragment"]
tab.ind <- tab[data$Type2 == "counts" & data$Type1 != "Fragment",]
tab.ind <- apply(tab.ind, 2, function(x) tapply(x, w, sum))
tab.biom <- tab[data$Type2 == "biomass" & data$Type1 != "Fragment",]
tab.biom <- apply(tab.biom, 2, function(x) tapply(x, w, sum))

#Habitat
hab <- unique(data[,c("MU3", "Name")])
hab <- hab[order(hab[,2]),][,1]
levels(hab) <- c("C", "S", "G", "T")

#Basic community descriptors
x <- hab
levels(x) <- c(2, 3, 1, 4)
x <- factor(as.character(x))
tx <- apply(ifelse(tab.ind == 0, 0, 1), 1, sum)
ind <- apply(tab.ind, 1, sum) / (pi * 0.16^2)
biom <- apply(tab.biom, 1, sum) / (pi * 0.16^2)
df <- data.frame(tx, ind, biom)
d <- read.csv("reldist.csv", h = T)
d <- d[order(d$Name),]
unique(d$Name == rownames(df))
d <- d[order(d$Name),][,2]

#Tests
for(i in 1:3){
  w1 <- kruskal.test(x = df[,i], g = x)
  w2 <- dunn.test(x = df[,i], g = x, method = "bonferroni")
  print(c("Taxon richness", "Individual density", "Biomass density")[i])
  print(w1)
  print(w2)
}

pdf("Fig2.pdf", width = 7, height = 12)
col <- x
levels(col) <- c("gold3", "gold", "gold4", 1)
col <- as.character(col)
lab = c("Number of taxa",
        expression(paste("Individual density (ind./m"^"2", ")", sep = "")),
        expression(paste("Fresh biomass density (g/m"^"2", ")", sep = "")))
m <- matrix(1:8, nr = 4)
layout(m, width = c(0.17, 1), height = c(1, 1, 1, 0.2))
par(mar = rep(0, 4))
for(i in 1:3){
  plot.new()
  mtext(side = 2, at = 0.5, lab[i], cex = 1.4, line = -2.7)
}
plot.new()
par(mar = c(1.5, 1.5, 1, 1.5), mgp = c(3, 2, 0), cex.axis = 1.8)
for(i in 1:3){
  plot(d, df[,i], bty = "n", pch = 21, col = "white", bg = col, cex = 2.4, xlab = "", ylab = "",
       ylim = list(c(0, 30), c(0, 3000), c(0, 2500))[[i]], las = 1, xaxt = "n")
  axis(side = 1, at = seq(-0.8, 0.4, 0.2), labels = c(F, F, T)[i])
  lines(lowess(df[,i] ~ d, f = 0.5))
  abline(v = 0, lty = 5)
  title(main = letters[i], font.main = 1, adj = 0.98, cex.main = 3, line = -1)
  legend(-0.8, c(30, 2750, 2250)[i], pch = 21, pt.cex = 2.5, col = "white",
         pt.bg = c("gold", "gold3", "gold4", 1),
         legend = c("Crest", "Gentle slope", "Steep slope", "Trough"),
         bty = "n", cex = 1.8)
  legend(-0.45, c(30, 2750, 2250)[i], legend = c("a", "a", "b", "ab"), bty = "n", cex = 1.8)
}
plot.new()
mtext(side = 1, "Relative distance", cex = 1.4, line = -0.7)
dev.off()

#Taxon distributions
df1 <- tab.ind[,apply(tab.ind, 2, sum) > 0]
fac <- x
df1 <- apply(df1 / (pi * 0.16^2), 2, function(x) tapply(x, fac, mean))
df1 <- 100 * sweep(df1, 2, apply(df1, 2, sum), "/")
df1 <- data.frame(t(df1))
for(i in 1:ncol(df1)){
  w <- df1[,i]
  b <- c(0, 25, 50, 75, 101)
  w[w > 0] <- cut(w[w > 0], breaks = b, labels = 1:4)
  df1[,i] <- as.numeric(w)
}

dudi <- dudi.pca(df1, scale = F, scan = F)

df2 <- tab.ind[,apply(tab.ind, 2, sum) > 0] / (pi * 0.16^2)
df2 <- apply(df2, 2, function(x) tapply(x, fac, mean))
df2 <- data.frame(t(df2))
for(i in 1:ncol(df2)){
  w <- df2[,i]
  b <- c(0, 1, 10, 50, 100, 1000)
  w[w > 0] <- cut(w[w > 0], breaks = b, labels = 1:5)
  df2[,i] <- as.numeric(w)
}

pdf("Fig3.pdf", width = 4, height = 12)
table.value(df2[order(dudi$li[,1]),]^2, cleg = 0, csi = 2,
            col.labels = c("G", "C", "S", "T"), grid = T,
            clabel.row = 1, clabel.col = 1.2)
legend(0.5, 64.0, pch = rep(15, 3), pt.cex = c(0.5, 0.8, 1.3, 1.8, 2.2, 2.6),
       legend = c("  <1", "  1-10", "  10-50", "  50-100", "  >100"),
       bty = "n", y.intersp = 1.5, cex = 1)
dev.off()

#Richness per habitat
w <- apply(ifelse(df == 0, 0, 1), 2, sum)
w

#Specific taxa
w <- apply(ifelse(df == 0, 0, 1), 1, sum)
length(df[,3][df[,3] > 0 & w == 1])

#Absolute densities
tab <- tab.ind[,apply(tab.ind, 2, sum) > 0]
fac <- x
df <- apply(tab, 2, function(x) tapply(x, fac, mean))
round(t(df)[order(df[1,]),], 4)

#Between-habitat analysis
#trace(scatterutil.eti.circ, edit = T)
tab <- tab.ind[,apply(tab.ind, 2, sum) > 0]
dudi <- dudi.pca(log(tab + 1), scale = F, scan = F, nf = 3)
b <- bca(dudi, hab, scan = F, nf = 3)
randtest(b, 999)

pdf("Fig4.pdf", width = 11, height = 10)
li <- b$ls
li[,1] <- -li[,1]
li[,2] <- -li[,2]
co <- b$c1
co[,1] <- -co[,1]
co[,2] <- -co[,2]
par(mfcol = c(2, 2))
for(i in 1:2){
  par(font = 1)
  s.label(co, clab = 0, cgrid = 1.8, cpoint = 0,
          xlim = c(-0.45, 1.0), yax = c(2, 3)[i],
          sub = c("a", "c")[i], csub = 2.8, possub = "topleft")
  w <- co
  w11 <- w[abs(w[,1]) < 0.1 | abs(w[,2]) < 0.1,]
  w12 <- w[abs(w[,1]) > 0.1 | abs(w[,2]) > 0.1,]
  w21 <- w[abs(w[,1]) < 0.1 | abs(w[,3]) < 0.1,]
  w22 <- w[abs(w[,1]) > 0.1 | abs(w[,3]) > 0.1,]
  s.arrow(list(w11, w21)[[i]], clab = 0, yax = c(2, 3)[i], add.p = T)
  s.arrow(list(w12, w22)[[i]], clab = 1.1, yax = c(2, 3)[i], add.p = T)
  col = list(c(1, 1, "white"), c(1, "white", 1))[[i]]
  eigen.bar(b$eig, ratio = 0.3, col = col)
}
for(i in 1:2){
  par(font = 1)
  s.class(li, hab, cpoint = 2, clab = 0,
          pch = as.character(hab), cstar = 0, axesell = F,
          cgrid = 1.8, xlim = c(-4, 6), yax = c(2, 3)[i],
          sub = c("b", "d")[i], csub = 2.8, possub = "topleft")
}
dev.off()

#Interpretation of axes
b$c1[order(b$c1[,1]),]

#Within-habitat analysis
w <- wca(dudi, hab, scan = F)
x <- hab
levels(x) <- c(2, 3, 1, 4)
li <- w$li
li[,1] <- -li[,1]
c1 <- w$c1
c1[,1] <- -c1[,1]
pdf("w.pdf", width = 8, height = 16, family = "Times")
m <- matrix(c(1, 1, 2, 4, 1, 1, 3, 5), nr = 4)
layout(m)
s.arrow(c1, clab= 2, xlim = c(-0.5, 0.8), cgrid = 3)
eig(w$eig, col = c(1, 1, rep("white", length(w$eig) - 2)), ratio = 0.5)
par(mar = rep(0.1, 4))
for(i in 1:4){
  s.label(li, clab = 0, cpoint = 0, cgrid = 3)
  text(li[x == i,], as.character(hab[x == i]), cex = 3)
  s.class(li[x == i,], x[x == i], clab = 0, cpoint = 0, cstar = 0, axesell = F, add.p = T)
}
dev.off()

#Biological traits
tab <- tab.ind[,apply(tab.ind, 2, sum) > 0]
tr <- read.csv("traits.csv", h = T)
colnames(tr)
tab.tr <- as.matrix(tr[,9:ncol(tr)])
rownames(tab.tr) <- tr$tx
colnames(tab.tr)
tab.tr <- tab.tr[,-c(6:10, 22)]
blo <- rep(1:5, c(5, 4, 4, 3, 5))
tab.tr <- tab.tr[rownames(tab.tr) %in% colnames(tab),]
tab <- tab[,colnames(tab) %in% rownames(tab.tr)]
tab <- as.matrix(tab) %*% tab.tr
tab <- prep.fuzzy.var(data.frame(tab), as.vector(table(blo)))

mod <- c(
  "Very small", "Small", "Interm.", "Large", "Very large",
  "Sessile", "Tubicolous", "Crawler", "Swimmer",
  "0", "0-5", "5-15", ">15",
  "Depos.", "Suspens.", "Carn.-Scav.",
  "Surf. mod.", "Diffus.", "D. conv.", "U. conv.", "Regen.")
x <- hab
levels(x) <- c(2, 3, 1, 4)
x <- factor(as.character(x))
title <- c("Body mass", "Motility", "Burrowing depth (cm)",
           "Feeding type", "Sediment mixing type")

df <- data.frame()
for(i in 1:max(blo)){
  x2 <- rep(x, table(blo)[i])
  x1 <- rep(1:table(blo)[i], rep(nrow(tab), table(blo)[i]))
  y <- tab[,blo == i]
  y <- 100 * matrix(as.matrix(y))[,1]
  for(j in 1:max(x1)){
    w <- dunn.test(x = y[x1 == j], g = x2[x1 == j], method = "bonferroni", list = T)
    w <- cbind(w$comparisons, round(w$P.adjusted, 3))
    name.tr <- rep(title[i], nrow(w))
    name.mo <- rep(mod[blo == i][j], nrow(w))
    w <- cbind(name.tr, name.mo, w)
    df <- rbind(df, w)
  }
}

sig <- c(
  rep("a", 4),
  "a", "ab", "ab", "b",
  "a", "b", "b", "b",
  rep("a", 4),
  rep("a", 4),
  rep("a", 4),
  rep("a", 4),
  rep("a", 4),
  "a", "ab", "b", "ab",
  "a", "ab", "b", "b",
  "a", "ab", "b", "ab",
  "a", "ab", "b", "b",
  "a", "a", "bc", "abc",
  rep("a", 4),
  "a", "ab", "b", "ab",
  "a", "ab", "b", "ab",
  "a", "ab", "b", "ab",
  rep("a", 4),
  "a", "ab", "b", "ab",
  rep("a", 4),
  rep("a", 4)
)
blo.sig <- sort(rep(blo, 4))

pdf("Fig5.pdf", width = 12, height = 8, useDingbats = F)
par(mfrow = c(2, 3), mar = c(3, 5, 3, 1.5))
for(i in 1:max(blo)){
  x2 <- rep(x, table(blo)[i])
  x1 <- rep(1:table(blo)[i], rep(nrow(tab), table(blo)[i]))
  y <- tab[,blo == i]
  y <- 100 * matrix(as.matrix(y))[,1]
  mod.hab = rep(c("G", "C", "S", "T"), table(blo)[i])
  w <- eb(data.frame(x1, x2), y, horiz = T, depth = T, m = median, bar = "pc",
     pch.m = 21, col.m = "white", bg = 1, cex.m = 2, cex.axis = 1.2,
     pch.p = "|", col.p = "grey40", mod = mod.hab, cex.mod = 0.7,
     title = paste(letters[i], ") ", title[i], sep = ""),
     line.title = 1, cex.title = 1.7, return.x = T)
  from <- round(par("usr")[3])
  to <- round(par("usr")[4])
  at <- rep(2.5, table(blo)[i]) + (5 * 0:(table(blo)[i] - 1))
  mtext(side = 2, text = mod[blo == i], at = -at, line = 3, cex = 0.8)
  text(rep(101, length(unique(w))), unique(w), sig[blo.sig == i],
       las = 1, cex = 1.1, col = "grey50")
}
dev.off()
