Rarefaction curve

[TuomasBorman]

To visualize the effect of rarefaction, one can visualize the species richness as a function of size of random sample.

This can be done with vegan, however, we might consider adding ggplot implementation to miaViz.

library(vegan)
data("dune")
rarecurve(dune, step = 1, sample = min(rowSums(dune)))

[vsuksi]

iNEXT package has ggplot2 functionality which is useful for thinking about this and considering overlap.

It could be useful to allow coverage-based rarefaction in mia and the option to label the sample-based curves at specific Good-Turing coverage in miaViz rarefaction curves. As per Chao & Jost (2014), the slope == 1 - coverage. Coverage-based rarefaction lets us compare samples at same completeness level instead of the technical variable of library size. This is useful especially when sample-based curve shapes have varying slopes at a set library size. In the following example perhaps one would not do rarefaction to begin with but assuming that one wants to do rarefaction one would like to do rarefaction to the same slope/coverage deficit.

library(vegan)
library(iNEXT)
data(dune)

# Get min and max observed richness samples (for clear example)
max_rich <- which(rowSums(dune > 0) %in% max(rowSums(dune > 0)))
min_rich <- which(rowSums(dune > 0) %in% min(rowSums(dune > 0))) # also min coverage
dune <- dune[c(min_rich, max_rich), ]

# Get library size of higher coverage sample at coverage of lower coverage sample
coverage <- iNEXT::iNEXT(t(dune), datatype = "abundance", q = 0)$iNextEst$coverage_based
min_sc <- min(coverage[coverage$Method == "Observed", ]$SC)
n5 <- max(coverage[coverage$Assemblage == 5 & coverage$SC < min_sc, ]$m)

vegan::rarecurve(dune, sample = n5)

[TuomasBorman]

Thank you, seems useful. Can you check if there is ggplot function in ggvegan?

Your script seems rather straightforwards, would you be able to propose a miaViz function for it?

[vsuksi]

No rarefaction curve in ggvegan. I'm still looking into overlap but will most likely make a draft PR on basic rarefaction curve for miaViz to make this more concrete.

My interest (for research use case) is primarily visualizing a set slope (Good-Turing coverage) in coverage-based rarefaction. So additionally it would be useful to have:

1. Good-Turing coverage, perhaps as an option in mia::addAlpha()
2. coverage-based rarefaction option in the rarefaction infrastructure in mia

Any thoughts on these in the meanwhile?

[TuomasBorman]

Sounds good!

1. Have you checked vegan::specpool? Does that do similar things?
2. Sounds also good, but @antagomir might have better view on this.

[antagomir]

I think (1) is prmary, and (2) is interesting depending on availability of time to implement (compared to other ideas to implement).

[vsuksi]

Implemented basic rarefaction curve based on vegan::rarecurve #226

vegan::specpool estimates extrapolated species richness accumulated across samples, so not suitable for the basic individual-based rarefaction curve commonly used for QC. Do we also want to support species accumulation curves?

Might be odd to include Good-Turing coverage in addAlpha in that it is not an alpha diversity metric.

I now hesitate to include coverage-based rarefaction in the wider rarefaction code because it would likely be motivated to consider Chao et al. (2016) singleton estimator in the presence of sequencing errors but that is a lot of thinking and work and doesn't seem to be available except partially in supplement.

[antagomir]

1. species accumulation curves: can be useful for other purposes; can be a separate issue to be done later ifit cannot be implemented on the same go (?)
2. Good-Turing coverage: imo addAlpha could have any scalar summaries of the community (e.g. global measures, one number per sample)
3. coverage-based rarefaction: ok to skip

[vsuksi]

If the simple version wrapping vegan::rarecurve() finds support, I can do a second pass, create tests and finish up.

1. If we want a more general wrapper which also does species accumulation curves, I'm tempted to try to adapt iNEXT package which provides a rather general solution already. It has only a few dependencies. The major challenge if I recall was performance problems due to not being able to turn off steps, for example extrapolation.

2. I can add Good-Turing coverage to addAlpha. An alternative if we want to avoid generalising addAlpha to scalar summaries is to save coverage and other central parameters in colData when using plotRarefaction.

3. Just a note on general approach to rarefaction (for example, coverage-based) in mia framework, one could:

• Use miaViz rarefaction curve (iNEXT wrapper) plot for inspection and save rarefaction and coverage information to colData. A useful feature here would be to mark rarefaction curves at a set coverage.
• Change rarefaction functionality in mia such that sample can be a vector, for example library size at same coverage.
• Looking into separating saving of coverage etc. into a mia function so miaViz could be clean plotting interface (no luck). Alternative is to demonstrate general rarefaction approach in miaViz rarefaction curve examples where user would use iNEXT package to save values (such as library size at set coverage) which one would pass to other.fields for visualization and sample in rarefaction functionality.