Simon-

I agree that for the treatment effect, combining the dispersions across depths is the way to proceed with these data. I think your reasoning for depth as an ordered factor is sound and likely applies to the depth profile I am working with. I attempted to follow your suggestion (snipped below) to create a design that allows me to use depth as an ordered factor. However, I found it confusing that the design is ~treatment + depth. With depth as the second variable, doesn’t that mean that the dispersions will be combined across treatments? I thought the goal was the combine dispersions across depths so that the treatments could be compared.

To do this with DESeq2, specify depth in you sample table not as a
factor, but as numerical vector, e.g. by putting the depth in metres.
The estimated log2 fold change can then be interpreted as reciprocal
halving/doubling distances: A log2 fold change of +0.2, for example,
would mean that the taxon's abundance doubles every 5 metres of depth
(1/0.2=5), -0.4 would mean that it's abundance halves every 2.5 metres.
If this assumption of exponential change is not too far from reality
(and at least for light, it is true, and hence maybe also for
light-dependent microbes), you will get valid p values.

Your design for this should then be: ~ treatment + depth

I should caution that we haven't tested regression on continuous
variables much, so if it does not work, ask again. It should work but
may need some tweaking, especially you may need to switch off
coefficient shrinkage.
/snip

In the code below I used a ~TREATMENT+DEPTH design with depth as numerical vectors and the “None” (seawater) treatment as the base level of the treatment factor. My main concern is that the output results in hundreds of differentially abundant taxa (344 OTUs for the Dead vs Seawater comparison with p-adjusted values <1e-5). These results include both independent Filtering and the cooksCutoff (by default). Given the really high numbers I am inclined to think that a large number of these are false positives and we have modeled the dispersion incorrectly. What do you think of this assessment?
> jam
class: DESeqDataSet
dim: 2151 12
exptData(0):
assays(1): counts
rownames(2151): OTU1 OTU2 ... OTU2150 OTU2151
rowData metadata column names(0):
colnames(12): HD5_150L HD5_150D ... HD9_SW_300_22 HD9_SW_500_22
colData names(5): DEPTH TREATMENT Cmg.m2.day Nmg.m2.day Pmg.m2.day
> design(jam)
~TREATMENT + DEPTH
> colData(jam)$TREATMENT<-relevel(colData(jam)$TREATMENT,"None")
> colData(jam)$DEPTH<-relevel(colData(jam)$DEPTH,"150")
Error in relevel.default(colData(jam)$DEPTH, "150") :
  'relevel' only for factors
> jam<-estimateSizeFactors(jam)
> jam<-estimateDispersions(jam,fitType="local")
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
> plotDispEsts(jam)
> jam<-nbinomWaldTest(jam)
> resultsNames(jam)
[1] "Intercept"              "TREATMENT_Dead_vs_None" "TREATMENT_Live_vs_None"
[4] "DEPTH"
> DvN_clean<-results(jam,"TREATMENT_Dead_vs_None",independentFiltering=TRUE)
> DvN_clean<-DvN_clean[order(DvN_clean$padj,na.last=NA),]
> write.csv(DvN_clean,file="Dead_vs_Seawater.TreatDepth.iF_true.cook_true.csv")

Thanks,
Kristina

> sessionInfo()
R version 3.0.2 (2013-09-25)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods
[8] base

other attached packages:
 [1] DESeq2_1.2.5            RcppArmadillo_0.3.920.1 Rcpp_0.10.6
 [4] snow_0.3-13             baySeq_1.16.0           GenomicRanges_1.14.3
 [7] XVector_0.2.0           IRanges_1.20.5          BiocGenerics_0.8.0
[10] ggplot2_0.9.3.1         scales_0.2.3            phyloseq_1.7.9

loaded via a namespace (and not attached):
 [1] ade4_1.5-2           annotate_1.40.0      AnnotationDbi_1.24.0
 [4] ape_3.0-11           Biobase_2.22.0       biom_0.3.10
 [7] Biostrings_2.30.1    cluster_1.14.4       codetools_0.2-8
[10] colorspace_1.2-4     DBI_0.2-7            dichromat_2.0-0
[13] digest_0.6.3         foreach_1.4.1        genefilter_1.44.0
[16] grid_3.0.2           gtable_0.1.2         igraph_0.6.6
[19] iterators_1.0.6      labeling_0.2         lattice_0.20-24
[22] locfit_1.5-9.1       MASS_7.3-29          Matrix_1.1-0
[25] multtest_2.18.0      munsell_0.4.2        nlme_3.1-113
[28] permute_0.7-0        plyr_1.8             proto_0.3-10
[31] RColorBrewer_1.0-5   reshape2_1.2.2       RJSONIO_1.0-3
[34] RSQLite_0.11.4       splines_3.0.2        stats4_3.0.2
[37] stringr_0.6.2        survival_2.37-4      tools_3.0.2
[40] vegan_2.0-9          XML_3.95-0.2         xtable_1.7-1

------------------------------------------------------------------
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 13, 2014, at 5:54 PM, Simon Anders <anders@embl.de<mailto:anders@embl.de>> wrote:

Hi Kristina

On 13/01/14 22:08, Kristina M Fontanez wrote:
Thank you for your insights. Replicates aren’t an option often in marine
metagenomics because many of us aren’t working at depths where dropping
a bucket is not a feasible option. Notably, we are working with
medium-term deployments of in situ incubations where the number of
sampling spots was extremely limited. It’s a first pass at a difficult
issue in marine microbiology but until we have the physical
infrastructure (sampling mechanisms) to get a large number of samples,
we often sacrifice replicates for better coverage across depths. Looking
forward, replicates is definitely a goal.

First, you may (or may not) have a valid substitute for replication, see
below, and if so, there might be little to complain about your design.

Of course, doing all four depths twice might be overkill for a first
try. If I had money for four samples, I might have opted for three
depths and taken, say, the middle depth twice, to get some idea how much
variation there can be. Precisely because this is a pilot study and a
first try for a new field, it seems important to me to establish that
the measurements for a given depth are stable and reproducible, _before_
rolling out a bigger experiment, and this is why I'm always so puzzled
that people consider replicates as important for follow-ups but not for
pilot studies. I would say it is the other way round.

The treatments are in situ incubations representing live microbes,
poisoned microbes (dead) and surrounding seawater (not an incubation,
can be thought of as a control).

Your point about guessing the dispersions is a difficult one for me to
agree with. First, the differences between the treatments are likely to
far exceed the differences within the depths. Preliminary analyses in
bayseq where samples are grouped by treatment bear this point out. So, I
think it’s reasonable to combine the dispersions across depths so that I
can compare live and dead treatments.

Actually, I fully agree.

If treatment causes much stronger differences than depth, then using the
samples from different depths as replicates to assess the effect of
treatment is entirely reasonable and valid. I somehow thought you wanted
to to it the other way round (pooling treatments to find differences
between depths -- though this is not what you wrote), and this won't
work, because nothing will be significant.

Second, as you point out, guessing
the appropriate dispersion makes the data difficult to publish given I
won’t have a good reason to argue for any value. These particular marine
metagenomes would be the first of their kind so there really isn’t a
good reference point.

For the treatment effect, you are fine because you can argue that the
dispersion estimated across different depths is certainly a conservative
estimate for what you would have gotten from replicates -- so you are
fine there.

Problem is, you cannot reason about the effect of depth this way, but
you would not want to throw away your knowledge about depths.

However, I neglected in my first post that depth is an _ordered_ factor,
and this might make quite a difference.

I am not a marine biologist, but I imagine that most taxa will either
keep increasing or keep decreasing in abundance when you go deeper and
deeper, at least if your four depths are not that extreme. I assume that
only few of your taxa will be happiest or unhappiest at the medium
depths, i.e. have strongest or weakest abundance not at the most shallow
or deepest sample but in one of the middle ones.

If so (and only if so), you have a good substitute for replication.
Strong differences between depths would then nearly always be
consistent: the direction of change for all three steps (from depth 1 to
depth 2, from depth 2 to depth 3 and from depth 3 to depth 4) is the
same. Taxa without depth preference will fluctuate randomly: the
differences between adjacent depth are not only weak but also
inconsistent in their directions. Hence, the taxa with the strongest
differences between depths will nearly always be consistent while those
with the weakest differences will have random directions of change, and
between the two, the line for significance can be drawn. This is why an
ordered factor can substitute for replication.

To do this with DESeq2, specify depth in you sample table not as a
factor, but as numerical vector, e.g. by putting the depth in metres.
The estimated log2 fold change can then be interpreted as reciprocal
halving/doubling distances: A log2 fold change of +0.2, for example,
would mean that the taxon's abundance doubles every 5 metres of depth
(1/0.2=5), -0.4 would mean that it's abundance halves every 2.5 metres.
If this assumption of exponential change is not too far from reality
(and at least for light, it is true, and hence maybe also for
light-dependent microbes), you will get valid p values.

Your design for this should then be: ~ treatment + depth

I should caution that we haven't tested regression on continuous
variables much, so if it does not work, ask again. It should work but
may need some tweaking, especially you may need to switch off
coefficient shrinkage.

One remaining question is whether I need the design ~ 1 option. I’m still
not clear on that.

When you create a DESeqDataSet, you have to pass a design formula.
However, the 'estimateSizeFactor' function's output does not depend on
this. The function 'rlogTransformation' takes an option, 'blind'. Its
default, TRUE, means that the function should _ignore_ the design
formula and dispersion values stored in the object and recompute them
using the design formula '~ 1'. With '~ 1', the function does not know
anything about the samples' assignment to groups, which is why we called
it 'blind'. If you use 'rlogData', you have to specify '~ 1' manually.

 Simon


	[[alternative HTML version deleted]]

