Skip to content. | Skip to navigation

Personal tools

Sections

R rarefaction

### 16s - rarefaction
# applied to taxotofreq_script output
# paper version


rm(list=ls())
setwd()

library(vegan)

freq=read.table("filename_taxotofreq.txt", header=T, sep="\t")

library(data.table, quietly = T)
for (j in 5:ncol(freq)){
set(freq,which(freq[[j]]<0),j,0)
}

freq$unclassified=NULL


S <- specnumber(freq[,-(1:4)]) # observed number of taxa
raremax <- min(rowSums(freq[,-(1:4)])) # minimum frequency

out=rarecurve(freq[,-(1:4)], step = 1, sample = raremax, col = "blue", cex = 0.6)

Nmax <- sapply(out, function(x) max(attr(x, "Subsample")))
Smax <- sapply(out, max)
library(RColorBrewer)
cols=colorRampPalette(colors = c("black","red","purple","blue","green","gold","orange","brown"))(32)

plot(c(1, max(Nmax)), c(1, max(Smax)), xlab = "Sample Size",
ylab = "orders", type = "n")
abline(v = raremax)
for (i in seq_along(out)) {
N <- attr(out[[i]], "Subsample")
lines(N, out[[i]], col = cols[i], lwd=0.5)
}
legend("topright", legend=freq$sample.ID, col=cols, pch=16, cex=0.8)
title("Family rarefaction curve")

R=rarefy(freq[,-(1:4)],raremax)
freqR=rrarefy(freq[,-(1:4)],raremax)

freqR=data.frame(freq[,1:4], freqR)
write.table(freqR, "filename_taxotofreq_rarefied.txt", sep="\t", row.names=F, quote=F)