The default multiple testing correction applied in the 'results' function is Benjamini-Hochberg. However, OTUs are not independent. There is information available about the specific structure of their relations in the form of a phylogenetic tree.
The package structSSI outlines a hierarchical procedure for multiple testing correction that takes this specific structure of the hypotheses into account.
I'd like to use this hierarchical FDR correction for DESeq results. I understand that running DESeq involves 3 functions.
dds <- estimateSizeFactors(dds) dds <- estimateDispersions(dds) dds <- nbinomWaldTest(dds)
The first two can be run once for the entire dataset (at the OTU level). I then want to be able to fit negative binomial models for the combined abundance of all of the children of each node in the phylogenetic tree somehow using the information from the dispersion fit.
Is there a way to extract the dispersion fit calculated on all the species for use in the model fitting, given that this will be on a new DESeq dataset containing the combined abundances of lots of OTUs?