0

Simplified version of my problem:

I have an OTU table with with 60 sample points (each with abundances for various OTUs) across three sites: A, B and C. Each site has 20 samples.

I want to plot rarefaction curves for each of the sites: A, B and C. I want to look to see if the curves for the sites plateau- to check if 'enough' sequences were sampled from each site.

At the moment I'm doing this:

raremax <- min(rowSums(otu.table))
rarecurve(otu.table, sample = raremax)

Currently every single sample point is being plotted as individual curves- this is not what I want- I want to plot the overall curves for each sampling site.

Frodo Baggings
  • 57
  • 3
  • 6
  • Do you mean the site average of rarefactions for individual sampling units, or the rarefaction of the site sums of sampling units? These are different. (This comment contains a hint how to achieve these results.) – Jari Oksanen Aug 03 '17 at 06:18
  • To be completely honest, I'm not sure- but I would hazard a guess that it's the the second option; I need to estimate whether the amount of sequenced samples for each site was sufficient to cover all the OTUs and so to accurately estimate species richness. Thanks! – Frodo Baggings Aug 03 '17 at 09:14

1 Answers1

0

For the second option, you must sum the species counts by sites. This will ignore all differences among sampling units by sites (that some people call beta diversity), and treat sampling units as simple replicates. You may need to justify your approach to reviewers if you want to publish your results.

With standard R you can use aggregate to get the sums. The following example uses vegan data sets mite and mite.env:

rarecurve(aggregate(mite, list(mite.env$Shrub), sum)[,-1])

The aggregate will add one column that you must remove, and you will not have labels for the curves. A little bit more winding way gives you also the names:

rarecurve(sapply(mite, function(x) sapply(split(x, mite.env$Shrub), sum)))
rarecurve(t(sapply(split(mite, mite.env$Shrub), colSums))) # alternatively

For this your data must be a data.frame--- matrix will fail.

Jari Oksanen
  • 3,287
  • 1
  • 11
  • 15