Thanks Christopher, I've added a photo of my rarefaction, is it still however showing uneven sequencing depth though?

Hi Christopher, thank you. So I first need to group all sequences passing 97% filter and then rarefy, using normalized number of reads?

I'd also like to draw attention to this paper by McMurdoe (http://arxiv.org/abs/1310.0424) called 'Waste not want not: Why rarefying microbiome data is inadmissable" as it brings up some interesting points. I have (as you can see) different sequence cluster abundances in my samples. If I am just characterizing the communities, seeing which microbes are in each sample, is it possible to do something else as the paper suggests so as not to get rid of other sequences, or to just leave the curves as they are?

If you ever want to compare any beta-diversity index, you have to rarefy/subsample your samples to equal sequencing depth.

If you have biom-format OTU table, for instance, you can run a qiime script called "single_rarefaction"http://qiime.org/scripts/single_rarefaction.html

Hi Fang, what if I'm just visualising what otu's are in each environment or anatomical site? Is rarefaction needed then?

Hi Alex, I mean that by bringing the no of sequences down to that of the sample with the lowest number of reads, I lose the reads from the samples with the greater numbers 

It all depends on what question you want to answer.

If you want to compare your samples to each other, then you must rarefy, in order to make a fair comparison. Otherwise the difference will mainly appear due to library sizes.

If you want to make a catalogue of OTUs found in your samples, without comparing them, then by all means use all of your data.

Dear Philip,

you pointed yourself to the PLOS CB paper by McMurdie & Holmes, 2014 which I believe gives some very good theoretical and pragmatic reasons why it is almost never necessary (or indeed a good idea) to rarefy your data. However, there is clearly some disagreement about this among authors, as you have pointed out in your OP. For example, the default QIIME SOPs include data rarefying steps, and you may indeed run into discussions with reviewers if you do not rarefy.

I have found the solution proposed in the McMurdie & Holmes paper quite helpful in many applications: they have implemented data transformation procedures from existing methods developed for gene expression analysis (edgeR & DESeq; basically, this is a transformation of the OTU table). Also, the implementation of these techniques in phyloseq works like a charm (see attached link for a tutorial).

However, this will not provide a solution to all your problems and neither "fix" all your un-rarefied analyses. In particular, some diversity estimators can be rather sensitive to varying sample sizes, while others prove to be more robust. Personally, I do usually not rarefy e.g. prior to alpha diversity estimation, but rather kick out too small samples entirely. Similarly, with weighted UniFrac (preferably even on unclustered trees, i.e. trees of "100% OTUs") you should also be on the safe side. The easiest way to test this is to simply run your analysis twice: once on rarefied and once on unrarefied data (the latter potentially corrected using the above-mentioned procedure by McMurdie & Holmes) and compare your results. Since you are dealing with OTU tables, the computational overhead of this should not be too large. In my experience, for sufficiently deep NGS data (and the rarefaction curves you posted look deep enough to me) you have a good chance of getting very similar results using both approaches for several standard ecological analyses.

Hope this helps,

Sebastian

http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003531

http://joey711.github.io/phyloseq-extensions/edgeR.html

Hi Philip, if you just want to visualize the OTU distribution across different samples, you might want to use heatmap. But the data need to be transformed to relative abundance（%），rather than the original reads number.

Thank you all very much for your very helpful answers. I really appreciate it

I'm looking at the rarefaction debate at the moment, and my natural leaning is towards 'don't rarefy', but it seems the main reason for rarefying is for unbiased diversity estimation.  Rather than sticking with standard pipelines, are there other diversity estimates which do account for sampling depth?

I have spent some time looking at this debate. I do not know of any diversity estimators that account for varying sampling depth. Typically these metrics were designed for sampling plant communities where you would sample the same size plots across the landscape. The issue has arisen with next-gen sequencing data and I doubt anyone has developed those methods yet. Would love to hear though if someone knows of one.

The best way I plan to proceed with my 16s amplicon data set is remove samples that had low sequences -- now this is an arbitrary distinction. I suppose you could look at the distribution of sequences per sample and try to come up with a standard place to make the cut-off point. My samples vary from 250 to 10,000 sequences per sample. With only 2-3 samples below 500 and the majority around 1,000-3,000 sequences. I removed samples that had below 500 sequences.

Then proceed with alpha and beta diversity metrics. I typically do not use any abundance weighted metrics here. If you do such as Shannon for alpha diversity or Bray-Curtis for beta diversity you should normalize your abundance counts beforehand. I have been using MetagenomeSeq to normalize the counts. Here's a post from github https://github.com/nosson/metagenomeSeq/issues/16 that I've been trying to figure out how this normalization works, but my explanations are not correct in the post. I am having a hard time understanding the mathematical jargon but Joseph Paulson is supposed to provide an example soon to help understand how cumulative sum scaling method works. If anyone can explain it I would love to hear.

Then I have been looking into how to compare relative abundances. You would first filter your dataset to only contain the 'common' OTUs to help reduce the amount of comparisons you are doing. I think some of the code may do this for you though....You can use MetagenomeSeq and their zero-inflated gaussian model to look at differential abundances. Also, there are some tutorials with Phyloseq in R. http://joey711.github.io/phyloseq-extensions/DESeq2.html. I haven't completed my reading on these statistical analyses yet though. This is as far as I can comment. 

Carly, thanks for those useful comments and links.

I also came across Simpson's Ds or Hurlbert's (1971) Probability of an Interspecific Encounter (PIE). PIE seems to be independent of sampling depth, so seems like it may be a good option to pursue in absence of rarefaction. It may have other downsides (doesn't everything?) that I'm not aware of yet.

I also read into the zero-inflated Gaussian models for differential abundance testing, and I agree these look like a very promising approach.

Christopher, I looked again at your paper today, and you also evaluated Simpson's metric and found it robust to sample size.  I couldn't figure out the scale quite from the image in the RG copy of the paper, what minimum library size did you find you needed for Simpson's Ds to stabilise?

Given that rarefaction is controversial, I'm increasingly liking the idea of not rarefying, and using Simpson's Ds for diversity estimates, but if anyone has any contraindications, I'd very much welcome hearing them.

New paper giving some guidelines on the best way to approach this type of data and associated statistical issues -- https://peerj.com/preprints/1157/

Effects of library size variance, sparsity, and compositionality on the analysis of microbiome data

Thanks Giorgio Marozzi for sharing recently below. This looks like the peer-reviewed version of the original paper I posted:

https://www.researchgate.net/publication/314194282_Normalization_and_microbial_differential_abundance_strategies_depend_upon_data_characteristics

Article Normalization and microbial differential abundance strategie...

 Thanks Carly, that's a really useful paper!  I guess it may change a little as it goes through peer-review, but I am told the authors are highly respected in the field.

I want to suggest this paper on data normalization. It is extremely interesting and it helps giving some guidelines on the best way to approach a 16s dataset:

Normalization and microbial differential abundance strategies depend upon data characteristics

http://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-017-0237-y

Hello there,

Rarefaction is a way to normalize your samples, so that you are not comparing 150,000 samples to 35,000. This is important if you are comparing the absolute abundances of OTUs in your samples. I ONLY compare relative abundances and I choose to rarefy when quantifying Alpha-diversity, but with beta-diversity I use non-rarefied samples. I'm not convinced that either one is better than the other. Just make sure that you are CLEAR in your methods, so that the readers understand how you analyzed your data.

Good luck!

I am going to discard 5 or 6 samples if I have to rarefied as the sequences is below 1000, If I discard what remain, suppose If i have three biological replicates or four replicates , in those 5, two are from one treatment and three of another treatment, It means I am losing three replicates of one sample and one replicates of other sample. Now when to calculate their mean will be low than other samples, this is biased, I will never rarefy the samples, if reviewer reject, I will tell them my final sample has 4000 sequences and last one has 200, therefore I donot want to loss those, so I did not rarefy,