Skip to content. | Skip to navigation

Personal tools

Sections

R ecodistances fig6

### 16s - ecological distances - figure 6
# applied to rarefaction_script output
# paper version

rm(list=ls())
setwd()
options(stringsAsFactors = T)
freq=read.table("filename_taxotofreq_rarefied.txt", header=T, sep="\t")

library(vegan)
Methode="bray"
DVEC=c()
GROUP=c()

for(j in levels(freq$sample.type)){
for(k in levels(freq$level)){
name=paste(k,j,sep="_")
sub=freq[freq$sample.type==j & freq$level==k, -(1:4)]
if(nrow(sub)!=0){
d=vegdist(sub, method = Methode)
}else{
d=NULL
}

if(length(d)!=0){
d=as.matrix(d)
n=nrow(d)
for(i in 1:(ncol(d)-1)){
dvec=d[(i+1):n,i]
DVEC=c(DVEC,dvec)
group=rep(name,length(dvec))
GROUP=c(GROUP,group)
}
}
}
}

DIST=data.frame(Distance=DVEC, Condition=GROUP)

boxplot(Distance~Condition,DIST, ylab=paste(Methode, "distances"))
stripchart(Distance~Condition,DIST,
vertical = TRUE, method="stack",
pch = 21, cex=1.2, col = "black", bg="darkgoldenrod1",
add = TRUE)
title(paste("Intra-condition", Methode, "distances"))
pairwise.t.test(DIST$Distance,DIST$Condition,p.adjust.method = "none")


## clustering based on ecological distances
rownames(freq)=freq$sample.ID
groupCodes <- paste(freq$sample.type, freq$level)
groups=unique(groupCodes)
colorCodes <- data.frame(group=groups,colors=c("mediumblue", "tomato", "mediumvioletred", "goldenrod", "dodgerblue"))
d=vegdist(freq[,-(1:4)], method=Methode)
h=hclust(d, "ward.D2")
plot(h)
dend <- as.dendrogram(h)
library(dendextend)
labels_colors(dend) <- as.character(colorCodes$colors[match(groupCodes,colorCodes$group)][order.dendrogram(dend)])
plot(dend, ylab="Height")
legend("topright", legend=colorCodes$group, col=as.character(colorCodes$colors), pch=16)