hi Noa,

When I look up the column data for the dataset you sent, I don't see the
samples grouped as I was expecting from your email.

I show below how it is possible to email the column data, even if it is not
for public viewing, by "scrubbing" the sample and factor level names. (Note
though that this does not change the results names stored in the
DESeqDataSet though, so you have to use the original factor level names
used when DESeq() was run, in order to extract results).

> library(DESeq2)
> load("~/Downloads/example.RData")
> dim(ddsSub)
[1]  1 45
> nlevels(colData(ddsSub)$condition)
[1] 26
> colnames(ddsSub) <- paste0("sample",seq_len(ncol(ddsSub)))
> colnames(ddsSub)
 [1] "sample1"  "sample2"  "sample3"  "sample4"  "sample5"  "sample6"
 [7] "sample7"  "sample8"  "sample9"  "sample10" "sample11" "sample12"
[13] "sample13" "sample14" "sample15" "sample16" "sample17" "sample18"
[19] "sample19" "sample20" "sample21" "sample22" "sample23" "sample24"
[25] "sample25" "sample26" "sample27" "sample28" "sample29" "sample30"
[31] "sample31" "sample32" "sample33" "sample34" "sample35" "sample36"
[37] "sample37" "sample38" "sample39" "sample40" "sample41" "sample42"
[43] "sample43" "sample44" "sample45"
> levels(colData(ddsSub)$condition) <-
paste0("lvl",seq_len(nlevels(ddsSub$condition)))
> levels(colData(ddsSub)$condition)
 [1] "lvl1"  "lvl2"  "lvl3"  "lvl4"  "lvl5"  "lvl6"  "lvl7"  "lvl8"  "lvl9"
[10] "lvl10" "lvl11" "lvl12" "lvl13" "lvl14" "lvl15" "lvl16" "lvl17" "lvl18"
[19] "lvl19" "lvl20" "lvl21" "lvl22" "lvl23" "lvl24" "lvl25" "lvl26"

The normalized counts and conditions listed in your email do not seem to
correspond in the condition column. You had:

T0: 29852 and 39519; T1: 37352 and 44385; T3:52539 and 58,556.
> The log2 fold change of T0/T1 was 3.28 (p-adjusted was 0.04), while the
> log2 fold change of T2/T0 was 2.96 (p-adjusted was 0.105)


The normalized counts alongside the condition:

> nc <- round(counts(ddsSub,normalized=TRUE)[1,])
> data.frame(condition=colData(ddsSub)$condition, normcounts=nc)
         condition normcounts
sample1       lvl1       4267
sample2       lvl2       4285
sample3       lvl3       5235
sample4       lvl4       2406
sample5       lvl5       3312
sample6       lvl4       7117
sample7       lvl5       3356
sample8       lvl1       2375
sample9       lvl2       5355
sample10      lvl3       3746
sample11      lvl6        284
sample12      lvl7       8375
sample13      lvl8        313
sample14      lvl6       3865
sample15      lvl7        400
sample16      lvl9       9866
sample17     lvl10        223
sample18     lvl11        370
sample19      lvl9        365
sample20     lvl10        431
sample21     lvl11        452
sample22     lvl12      25262
sample23     lvl13      18322
sample24     lvl14      23999
sample25     lvl12      20305
sample26     lvl13      35235
sample27     lvl14      31593
sample28     lvl15      32076
sample29     lvl16      32923
sample30     lvl17      40614
sample31     lvl18      42732
sample32     lvl19      36969
sample33     lvl20      45901
sample34     lvl21      29852
sample35     lvl22      39520
sample36     lvl23      37353
sample37     lvl24      44386
sample38     lvl25      52540
sample39     lvl26      58556
sample40     lvl18      16779
sample41     lvl19      20490
sample42     lvl20      18864
sample43     lvl15      17213
sample44     lvl16      16181
sample45     lvl17      16639

The normalized counts for 29852 and 39519 are in different levels: lvl21
and lvl22

37352 and 44385 are in lvl23 and lvl24

52539 and 58,556 are in lvl25 and lvl26

If I extract the contrast for, say, lvl23 vs lvl21 (using but scrubbing
your original level names), the log fold change (having run DESeq with
betaPrior=FALSE) is equal to the log ratio of the condition means:

> results(ddsSub, contrast=c("condition","lvl23","lvl21"))
DataFrame with 1 row and 6 columns
           baseMean log2FoldChange     lfcSE      stat    pvalue      padj
          <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
####  18237.81      0.3233761  1.492571  0.216657 0.8284756 0.8284756

> log2(mean(nc[colData(ddsSub)$condition == "lvl23"]) /
mean(nc[colData(ddsSub)$condition == "lvl21"]))
[1] 0.3233965


It might be easier for you to subset the DESeqDataSet to smaller groups
before running DESeq(), and to examine
as.data.frame(colData(dds)$condition) before running to make sure that the
sample and conditions line up as you expect. Maybe I should address the
question: what's the advantage/disadvantage to running all the samples in
one go, vs subsetting and running separately? One advantage of 'all in one'
is convenience, in that any contrast can be extracted after the model has
run. Another advantage is that if there are replicates for other
experimental conditions, and we expect that the dispersion within these
experimental conditions is similar to the dispersion within the conditions
we are interested in, then estimating the dispersion altogether should
provide better estimates. However, this assumes that the dispersion will be
similar across these groups, and that you have replicates for these
"out-group" conditions. If that is not the case, I recommend splitting the
analysis into subgroups of interest. A disadvantage of 'all in one' is that
the running time is longer, because for each gene we have to find the
maximum of the likelihood or posterior in a high dimensional space. It's
much faster to find this maximal point in a lower dimensional space (when
there are only a few variables, and a few levels per variable). Another
disadvantage of 'all in one' is that it becomes a bit more difficult to
keep track of samples and conditions.

Regarding the email from your coworker (sorry to have missed it, I couldn't
find this email in searching the Bioc mailing list)

I also wanted to raise an additional question that my colleague sent few
> weeks ago to the bioC list and was skipped. It's related to version 1.2.5
> but still we would like to better understand the following issue:
> A Dataset contains 4 different conditions with 3 replicates each.
> In one of the genes examined, the raw counts of all 3 replicates of two
> conditions (I will call them A and B) are zero (while other conditions
> received non-zero values).
> When viewing the results given by the DESeq2 package (after building the
> data set from a matrix, estimating size factors and dispersions, using the
> Wald test and the function results with contrast) I see the following:
> In the comparison : A vs. B, even though all samples received zero counts:
> log2FoldChagne = -1.08976
> p-value = 0.3044
> padj = 0.36323
> The threshold set by the independent filtering for that comparison was
> around 1. If all replicates of conditions A and B received zero counts,
> shouldn’t the independent filtering filter this gene from this comparison?
> (meaning the padj value should be NA) or is the condition checked by the
> independent filtering refers to all samples in the data (and not only those
> of the conditions being compared), even though the threshold is set for
> each comparison separately?


That the log2FoldChange is not equal to zero for this contrast, can be
answered by the thread I originally posted. (
http://permalink.gmane.org/gmane.science.biology.informatics.conductor/51749
​
 ).  In DESeq2 version 1.2.x, when a DESeqDataSet is constructed with
factors in the design with 3 or more levels, we recommend setting
betaPrior=FALSE. (this is printed as a message in the construction of the
DESeqDataSet in the updated release version since November). Then there
might be a nonzero log2FoldChange, but the standard error will be very
large, as the profile likelihood for this coefficient is sloping toward
negative infinity, so there will not be anything near significance. In
DESeq2 version 1.3.x (devel version which will be released in Spring 2014)
we have a comprehensive solution for using a beta prior when factors are
present in the design formula with 3 or more levels.

Yes, the independent filtering operates on the mean of the normalized
counts across *all* samples in the DESeqDataSet. The particular threshold
is set by optimizing the number of rejections based on the adjusted
p-values from the test currently being generated. However, such an
operation as you describe can be easily performed by supplied any kind of
independent filtering statistic to the 'filter' argument of the results
function. For filtering based on the mean of only the two conditions
currently being investigated, this would look like:

> dds <- makeExampleDESeqDataSet()
> colData(dds)$condition <- factor(rep(c("a","b","c"),each=4))
> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> bcMean <- rowMeans(counts(dds,normalized=TRUE)[,colData(dds)$condition
%in% c("b","c")])
> results(dds, contrast=c("condition","c","b"), filter=bcMean)
DataFrame with 1000 rows and 6 columns
             baseMean log2FoldChange     lfcSE       stat    pvalue
 padj
            <numeric>      <numeric> <numeric>  <numeric> <numeric>
<numeric>
feature1    64.518521    -0.28198624 0.3521963 -0.8006507 0.4233339
 NA
feature2     3.954962     0.45020401 0.6358071  0.7080827 0.4788939
 NA
feature3    50.301443     0.04948108 0.4502842  0.1098886 0.9124978
 NA
feature4    77.491191    -0.34507847 0.3806662 -0.9065119 0.3646649
 0.917873
feature5    76.342876    -0.24041099 0.3612686 -0.6654633 0.5057542
 0.917873
...               ...            ...       ...        ...       ...
...
feature996   2.614749     -0.3276690 0.6030464 -0.5433561 0.5868846
 NA
feature997   9.856852     -0.4837603 0.5661241 -0.8545128 0.3928209
 NA
feature998   5.403219      0.4152869 0.6299118  0.6592778 0.5097174
 NA
feature999   9.738524     -0.1803672 0.5557705 -0.3245354 0.7455327
 NA
feature1000 40.187737      0.6003838 0.4249903  1.4127001 0.1577439
 NA

Mike




On Sun, Jan 12, 2014 at 11:39 AM, Michael Love
<michaelisaiahlove@gmail.com>wrote:

> hi Noa,
>
> Could you email me a small reproducible example? e.g.:
>
> ...
>  dds <- nbinomWaldTest(dds, betaPrior = FALSE)
> idx <- 1234 # the row you are referencing
> ddsSub <- dds[idx, ]
> save(ddsSub, "small_example.RData")
>
> best,
>
> Mike
>
>
>
> On Sun, Jan 12, 2014 at 9:03 AM, Noa Henig <ntzunz@gmail.com> wrote:
>
>> Hi Simon,
>> I forgot to mention that the counts are the *normalized* counts, I guess
>> this is the reason you asked for the size factors.
>> I ran DESeq2 again using the code below (this time with beta prior = FALSE
>> to handle the issue of zeros and low counts I originally had). The fold
>> change in this run was a little higher than before:
>> log2FC (T1/T0) = 3.97;  log2FC (T3/T0) = 3.64
>> How can this be explained considering the values of the normalized counts
>> for that gene? should I do some modification to adjust also for the high
>> counts?
>> Thanks !
>> Noa
>>
>>
>> This is the code:
>> countData = read.table("merged_counts.txt", header=TRUE, row.names=1)
>> colData = data.frame(row.names=colnames(countData), condition =
>> c(param[,2])) ## it's a long list of samples names
>> colData$condition = factor( colData$condition, levels =
>> c(param[1:numOfConditions,7]))
>>
>> dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData,
>> design = ~condition)
>> dds <- estimateSizeFactors(dds)
>> dds <- estimateDispersions(dds, maxit = max_iterations) ## a user-defined
>> parameter
>> dds <- nbinomWaldTest(dds, betaPrior = FALSE)
>> tmp_counts = counts(dds, normalized=TRUE)
>> normalized_counts = data.frame(tmp_counts)
>>
>> final_stats = NULL
>> filtering_thresholds = matrix(ncol = 2, nrow = as.numeric(param[1,3]),
>> byrow = TRUE)
>> for(i in 1:numOfComparisons)
>> {
>> curr_res = results(dds, alpha = 0.05, contrast = c("condition",
>> param[i,4],
>> param[i,5]))
>>          ..writing the data...
>> }
>>
>>
>>
>>
>> On Wed, Jan 8, 2014 at 3:21 PM, Simon Anders <anders@embl.de> wrote:
>>
>> > Hi Noa
>> >
>> > On 08/01/14 14:11, Noa Henig wrote:
>> > > I would expect, for example:
>> > > genes having these counts T0: 29852 and 39519; T1: 37352 and 44385;
>> > > T3:52539 and 58,556.
>> > > The log2 fold change of T0/T1 was 3.28 (p-adjusted was 0.04), while
>> the
>> > > log2 fold change of T2/T0 was 2.96 (p-adjusted was 0.105)
>> >
>> > Please tell us the values of the size factors (run "sizeFactors(dds)")
>> > to see how far off this is.
>> >
>> > And please always post your code, not just the output.
>> >
>> >   Simon
>> >
>> > _______________________________________________
>> > 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]]
>>
>> _______________________________________________
>> 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]]

