I described my most abundant species within groups by looking at relative abundance, and
in QIIME 1 I used “group_significance.py” (Kruskall-Wallis corrected by the Benjamini-Hochberg FDR) to look for differences between groups for taxa found in >25% of samples. The reviewer of my paper said this is not valid due to uneven group sizes (e.g. n = 13, 26 and 53).
When I asked a respected bioinformatician he felt the stats I had used had been reasonable, but unfortunately didn’t point me towards something I could say to the reviewer to justify keeping the results in my paper. Can anybody give me a good reason why I can (my preference) or can’t keep my analyses?
I am happy to do further analyses (for example DeSeq2) but I don’t want to have to change the whole paper if I don’t have to.
With regards to other methods, when I used DeSeq2 I felt it was picking up far too many differences based on really rare taxa. I’m interested in ones that are very prevalent and/or abundant as in general these tend to be more clinically relevant. Is there a good guide to what I can and can’t do with DeSeq2 in terms of filtering or reporting results?
Is there something else, easy to understand, that will help me look for taxa that are differentially expressed between groups of different sizes? I am working with fungal ITS data so can only use non-phylogenetic statistical approaches.
I am afraid that evil reviewer #3 is correct in this case: you do need to re-run your differential abundance tests to use something other than kruskal wallis (I have suggestions below)… a lot has changed since qiime1 was started over a decade ago, and differential abundance testing in particular is a hot area of research because microbiome datasets present several challenges that make many conventional statistical approaches prone to error… uneven group sizes is just the tip of the iceberg, sparsity and compositionality are much larger concerns.
Here are some good review papers on the topic:
Both of the papers I listed above recommend against using DESeq due to high false discovery rates, which is quite consistent with your observations of “too many differences” — a lot has changed with DESeq2 but this method may still not be appropriate for compositional microbiome data (though others have recommended it in the 6-yr-old literature… again, this is a rapidly evolving area of research!)
In general, filtering low-abundance features is okay (and is even recommended for some methods, like ANCOM, to reduce runtime if you are not interested in rare features). You can also sort your outputs based on abundance or prevalence to focus on biologically interesting features.
So what can you use instead?
QIIME 2 has a few differential abundance methods worth testing out here. The best for your case (as a compositionally aware replacement for kruskal-wallis) would be to use ANCOM (part of the q2-composition plugin that comes installed with the current QIIME 2 2020.2 release) or q2-aldex2 (installation instructions and details can be found in the library: https://library.qiime2.org/).
Thank you so much for your response and the papers, which I shall read in full.
A somewhat cheeky question, given this is the QIIME2 forum. I am trying to transition over to QIIME2, and will certainly be using it for future datasets, but the dataset in question all the work to date was done in QIIME 1 version 1.9.1. I am on a tight deadline to respond to reviewers questions to resubmit my article. I want to make as few changes as possible. Is there a way to use ANCOM with QIIME 1.9.1?
Although I’ve started to learn QIIME 2 I feel much more comfortable with QIIME 1 and its been a while since I started learning QIIME 2 so I’d almost be starting from scratch again. I am unable to install software myself on the University’s HPC system but it does have qiime2/2019.4.
You say “In general, filtering low-abundance features is okay” - do you have a level of filtering you would recommend for ANCOM or DESeq2 (as I could use QIIME 1.9.1 differential_abundance.py for DESeq2).
There is a way to use ANCOM without touching QIIME 2: use its R implementation.
If you are comfortable working with R, that would be the way to go. However as you are already familiar with QIIME 1, then importing your biom table and running this in QIIME 2 would be easier than putting together some R commands, since importing a biom table should be fairly seamless. Here’s what to do:
Two steps, a total of 4 commands. (and if you want to use q2-aldex2 instead of ancom, follow the first 2 steps then follow the ancom tutorial for the other steps)
That works! Since that version is a year old, some of the commands may have changed slightly, but if you run into trouble just open a new topic on this forum and we can help.
Wow that’s a big question — it all depends on what level of rarity you are interested in, and on the number of reads per sample you have in your dataset. So there is no one-size-fits-all, but in general rare OTUs tend to be spurious and should probably be filtered anyway — since you have OTU clustered data from QIIME 1 you could use that paper to choose a filtering depth and see how that impacts your data, then go from there.
Which reminds me: DESeq and kruskal-wallis are both prone to high false-positive error rates with microbiome data, so do not be alarmed if aldex2 or ANCOM yield far fewer significant features.