Skip to content. | Skip to navigation

Personal tools

Sections

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)