Directory
R diversity table2
### 16s - diversity - table 2
# applied to rarefaction_script output
# paper version
rm(list=ls())
setwd()
freq=read.table("filename_taxotofreq_rarefied.txt", header=T, sep="\t")
library(data.table, quietly = T)
colnames(freq)[4]="exp"
for (j in 5:ncol(freq)){
set(freq,which(freq[[j]]<0),j,0)
}
# removal of unclassified + Wolbachia sequences
freq$unclassified=NULL
freq$Rickettsiales = NULL
freq$Anaplasmataceae = NULL
# nb of taxa in each group
test=freq[,c(5:ncol(freq))]
summary(test)
for(i in 1:nrow(test)){
test[i,]=round(100*test[i,]/sum(test[i,]),2)
}
prop=cbind(freq[,c(1:4)], test)
# threshold
s=1
library(data.table, quietly = T)
for (j in 5:ncol(prop)){
set(prop,which(prop[[j]]<s),j,0)
}
ALPHA=prop[,1:4]
alpha=apply(prop[,5:ncol(prop)],1,function(x){return(sum(x!=0))})
ALPHA$alpha=alpha
mod=lm(alpha~exp*sample.type*level, ALPHA)
anova(mod)
ANOVA <- drop1(mod, .~., test="F") # type III ANOVA
write.table(ANOVA, "filename_taxotofreq_rarefied__ANOVA.txt", sep="\t",dec="," ,row.names=T, quote=F)