hi Kristina,


On Thu, Jan 16, 2014 at 1:33 PM, Kristina M Fontanez <fontanez@mit.edu>wrote:

> Hi Simon-
>
> Thank you for your detailed reply. I agree that it is the case often in
> marine microbiology that the null hypothesis of no changes between
> treatments, often does not really apply. When you frame the discussion this
> way, it makes sense that I am seeing large numbers of differentially
> abundant taxa.
>
> Your discussion about how to narrow down the large number of
> differentially abundant is quite interesting. If I am understanding you
> correctly you are suggesting that a reasonable cutoff on adjust p-values
> (0.1 or 10% false discovery rate) is simply too lenient for my data. I
> agree and that is why I initially dropped to 1e-5 however that number is
> really quite abritrary.
>
> Prior to reading your e-mail I had been considering using the base Mean of
> normalized counts across samples as a primary cutoff. I believe that this
> is exactly what independent Filtering does, in that it attempts to find a
> specific number of normalized mean counts across samples which maximizes
> the number of taxa with a FDR below alpha (0.1). However, this optimization
> isn’t very effective for my data because in one comparison (Dead vs None)
> 655/1596 OTUs are below the 0.1 cutoff and yet have base mean of normalized
> counts less than 10. It filters the data somewhat but not to a very helpful
> degree.


​I think the point is that this is not even the kind of filter that you
would like. It is filtering based on the tests that the log fold changes
are zero, when what you would like is ​to filter based on a test that the
absolute log fold changes are larger than some threshold greater than zero.
We have in fact implemented this in the development version of DESeq2 (v1.3
paired with Bioc 2.14), as arguments to the results() function, which will
be released in a few months. I should resist pasting some code in this
email thread, which might become out of date and is not under source
control. You can try to download the Package Source, a .tar.gz file here,
http://www.bioconductor.org/packages/devel/bioc/html/DESeq2.html, and then
install this by filling in 'xx' with the appropriate version number you
download:

install.packages("DESeq2_1.3.xx.tar.gz",repos=NULL)



> In order to use it effectively I would need to lower alpha to 1e-5 or
> something similar, which is again - an arbitrary cutoff.
>
> In the end there are really two questions I need to answer:
> 1) How many taxa are “truly" differentially abundant?
>
> 2) Of those truly differentially abundant taxa, what log2fold changes are
> biologically meaningful?
>
> In regards to (1), it seems that based on your email and what I’ve read in
> the DESeq manual that anything with a FDR of 10% ought to be regarded as
> “truly” differentially abundant. Is that correct? (Notably, if you turn on
> independent filtering there is a threshold selected for the base mean of
> normalized cutoffs that optimizes the total number of taxa below the 0.1
> (alpha) - so in that case the mean of normalized counts is indeed taken
> into account and well as the FDR). Cooks cutoff for outliers doesn’t apply
> for my samples because hugely differentially abundant taxa in my treatment
> are likely true biological events (microbes rapidly growing in response to
> influx of nutrients) rather than technical or experimental artifacts.
>
> In regards to (2), I think this is where the log2fold changes but not
> standard errors become useful. In my data, small standard errors (high
> precision) (< 0.3) are correlated with small log2fold changes (< 3) and low
> mean of normalized counts across samples (< 1000). If I choose to only keep
> those taxa with small standard errors then I will be throwing away some of
> the most differentially abundant taxa in my dataset (log2fold changes > 3,
> mean counts > 1000 and standard errors > 0.3). That seems backwards.
> Alternatively, if I choose to keep those taxa with high log2fold changes
> (>3), considering those to be biologically meaningful, then I reduce the
> dataset considerably.
>
> That leads me to the final and most difficult question: what level of
> log2fold change in a microbial community is biologically relevant on the
> time scales and in the ecological location in which these samples were
> collected?
>
> I can think of one possible ways to assess this:
>
> We have metagenomic samples of microbial communities taken at this
> location, at particular depths, going back several years (albeit sequenced
> with different genomic technologies). Previous research has shown that the
> microbial community in the seawater at any given depth doesn’t change very
> much from year to year. If I then make the assumption that any changes in
> the seawater community at a particular depth from one year to the next are
> the “baseline” level of variation that can be expected, I should be able to
> calculate the maximum log2fold change in abundance for any given taxa from
> one year to the next. I could then use that log2fold change as my threshold
> for biologically meaningful change when analyzing my treatments.
>
> Any insights you may have to these ideas would be greatly appreciated.
>
> Thank you,
> Kristina
> ------------------------------------------------------------------
> Kristina Fontanez, Postdoctoral Fellow
> fontanez@mit.edu<mailto:fontanez@mit.edu>
> Massachusetts Institute of Technology
> Department of Civil and Environmental Engineering
> 48-120E
> 15 Vassar Street
> Cambridge, MA 02139
>
>
>
> On Jan 16, 2014, at 10:02 AM, Simon Anders <anders@embl.de<mailto:
> anders@embl.de>> wrote:
>
> Hi Kristina
>
> I'll write a bit a longer reply, to advertise a kind of paradigm shift
> that we are now exploring in DESeq2.
>
> First a note to other readers of this thread: Kristina has sent me a
> private mail explaining the experiment in more detail. To sum this up
> without disclosing anything, the three "treatments" labelled "none",
> "life" and "dead" are not so much different treatments applied to one
> kind of sample but rather three different (actually: very different)
> manners in which the samples were collected, which are hence expected to
> capture quite different subsets of the species that are present at the
> sample collection spots.
>
> Therefore, it is no surprise that observed OTU abundances differ so
> greatly between the "treatments". Remember that the purpose of a test
> for differential expression in its usual sense is to find gene for which
> you can reject with confidence the null hypothesis that the gene is
> _not_at_all_ affected by the treatment, i.e., that its abundance stays
> exactly the same.
>
> It is quite common that this specific null hypothesis turns out to be a
> bit silly. All genes are connected by the cell's complicated regulatory
> network, and hence, there typically might not be a single gene which is
> not at least very slightly affected by the treatment. In essence, what a
> significant p value hence really means is that the effect was big enough
> that we can say with confidence, in which _direction_ the gene has
> changed. Quite often, the effect strength required to establish the
> direction of change is larger than the effect strength required to make
> the change of biological interest, and therefore, one looks at p values.
> If it is the other way round, one needs to also have a cut-off on fold
> change or do a banded test (see below).
>
> In your case, it is also clear a priori that no OTU will have the same
> abundance in two "treatments", given that the different sample
> collection approaches favour quite different species. Hence, it is no
> surprise that at a reasonable cut-off on the adjusted p values (say,
> 0.1), nearly all OTUs are significant.
>
> This means that you are in the favourable situation that you can take
> the estimated differences between treatments pretty much at face value.
> However, you want to know _how_precise they actually are, and for this,
> DESeq2 now reports standard errors alongside all estimated log fold
> changes. This data, and not the p values, is likely the more useful
> output for you.
>
>  Simon
>
>
>         [[alternative HTML version deleted]]
>
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>

	[[alternative HTML version deleted]]

