[TomHarrop]

Hi again,

I'm using the LRT with DESeq2 as follows.

> design(dds)
~batch + stage
> dds < - DESeq(dds, test = "LRT", reduced = ~ batch)

I see that I can access the calculated deviance for each gene:

> mcols(mcols(dds), use.names=TRUE)[c(20,25),]
DataFrame with 2 rows and 2 columns
                    type                                   description
             <character>                                   <character>
LRTStatistic     results LRT statistic: '~ batch + stage' vs '~ batch'
deviance         results                    deviance of the full model

Just so I can illustrate to some colleagues what I'm doing using some actual data, I'm wondering, for a given gene, if it's possible to retrieve the expression values predicted by the full and reduced models from dds after running the LRT?

Thanks,

Tom

[apredeus]

I see. I'll just add that call to my script then, should be enough. Thank you!

[Michael Love]

We cannot print it during model fitting, because we fit a model without knowing what comparison will be made. Will it be B vs A? Or C vs B? Or all the interaction terms together as a likelihood ratio test? This is why it is printed when you print the results object to the console, and the information is also stored as metadata columns:

mcols(res)

Check the DESeq2 vignette section on factor levels, which addresses your question.

R's factor function chooses a level ordering alphabetically unless you specify a meaningful order.

[apredeus]

How exactly have you called the "results" function? In any case, the first line of the output of the results function indicates the exact comparison that is reported. Can you check whether this line was correct in your case?

Oh, I see what the problem was. I have an R function that does the following:

res <- results(ddsMatrix)
resOrdered <- res[order(res$padj),]

Thus I did not see the output you are talking about. Still though, maybe it's a good idea to print it along with messages "estimating size factors" etc?

At any rate, thank you for your reply!

[Simon Anders]

How exactly have you called the "results" function?

In any case, the first line of the output of the results function indicates the exact comparison that is reported. Can you check whether this line was correct in your case?

[apredeus]

Here's a bit random and fascinating observation I've made when using DESeq2.

I've had three groups of experiments - with conditions labelled as "B", "12", and "24". I've ran pairwise comparison with donor-sensitive design on all three groups, getting differential expression for them in the same way:

B vs 12 (cond1 = B, cond2 = 12)
B vs 24 (cond1 = B, cond2 = 24)
12 vs 24 (cond1 = 12, cond2 = 24)

so, in first two cases, log2FC was reported for cond1/cond2 (i.e. positive numbers were for genes down-regulated at 12 or 24). In the latter case, log2FC was reported for cond2/cond1.

The differences were subtle enough that it took me a long long while to catch this error. I figured it's probably due to some sorting - either internal or in the condition file. Maybe it would make sense to write explicitly in the logs which one is cond1 and which one is cond2?

Cheers

[IsBeth]

Hi Michael,

You were absolutely right I had filtered the object returned from the results() function with a minimum log2foldChange of 1.5 and an adjusted p-value of 0.01. It was definitively too much. The lfc of the result data is less than 1 in most cases, and the adj p-value is superior than one...

With no filtering, I get the results table with the statistics I need, but no DE genes are detected and so none are plotted (with red points). Is there some way tu increase the statistical power (maybe in the results() function)?Or to change some parameters in order to mark genes as differentially expressed (even if their values are not magnific)?

Thank you for your time =)

[Michael Love]

DESeq() in DESeq2 also estimates p-values, and adjusted p-values (for those genes which pass independent filtering step).

I would guess that you are filtering with a p-value or adjusted p-value threshold which is less than all the p-values or adjusted p-values in the object returned by results(), and generating a 0-row dataframe.

[IsBeth]

The without-replicates analysis is not expected to produce small p-values or adjusted p-values. We have overestimated the variance by treating the samples from different conditions as replicates for the purpose of estimating variance

Ok, so when estimating the variance this way, it is not sufficient to detect DE genes? I know that to work with no replicates is not appropiate, but in these cases I need to provide some information of the data for the end user . Isn't there some way to obtain a DE genes list (even if the results are to be taken with care)?
(I think that it was possible with DESeq)

[Michael Love]

I guess the point I want to emphasize is the quote from the original DESeq paper, regarding the without-replicates analysis:

"Some overestimation of the variance may be expected, which will make that approach conservative."
Furthermore, "while one may not want to draw strong conclusions from such an analysis, it
may still be useful for exploration and hypothesis generation."

The analysis without replicates is very conservative, and only recommended for exploration (for example, you can look at which log fold changes are largest in the MA plot).

So to answer your question from before:

"it doesn't detect the DE genes. The result table is empty, and my MAplot doesn't show differential expression either"

The without-replicates analysis is not expected to produce small p-values or adjusted p-values. We have overestimated the variance by treating the samples from different conditions as replicates for the purpose of estimating variance.

[IsBeth]

Yes, I did. And I'm sorry, since I've seen that this problem is treated in other posts (such as "DESeq2 without biol replicates"). I tried the following code, but it doesn't work by now (I need to obtain a result table with statistics like the adjusted p-values and lfc):

design_vector <- c("C", "T")

coldat=DataFrame(grp=factor(design_vector), each=1)
dds <- DESeqDataSetFromMatrix(raw_filter, colData=coldat, design =~grp)
dds <- DESeq(dds)

rld <- rlogTransformation( dds )
NOTE: fitType='parametric', but the dispersion trend was not well captured by the
function: y = a/x + b, and a local regression fit was automatically substituted.
specify fitType='local' or 'mean' to avoid this message next time.

deseq2.res <- results(dds)


What do I have to change? Cause I'm sure I'm doing it completely wrong

[Michael Love]

When you ran DESeq(), you should have seen the following message printed to console:

estimating dispersion by treating samples as replicates.
read the ?DESeq section on 'Experiments without replicates'

Did you read this section?

[IsBeth]

Hello!

I'm trying DESeq2 with synthetic RNA-seq count data without replicates (40% of the synthetic genes are differentially expressed (DE), with a up - and downregulation ratio of 50%). Although DESeq2 detects these ratios correctly, it doesn't detect the DE genes. The result table is empty, and my MAplot doesn't show differential expression either. Is there any solution to this? Do I have to change any parameter?

[Gonza]

Hi all,

I am new to DEseq2. Maybe this is silly question, but I do not understand why meanSdPlot fails.... any ideas to solve this problem?


library("vsn")
par(mfrow=c(1,3))
notAllZero <- (rowSums(counts(dds))>0)
meanSdPlot(log2(counts(dds,normalized=TRUE)[notAllZero,] + 1))
meanSdPlot(assay(rld[notAllZero,]))
meanSdPlot(assay(vsd[notAllZero,]))

Error: could not find function "meanSdPlot"


Thanks

[Neuromancer]

Dear Michael and Devon,

thank you for taking the time to think about this and write an answer.
These are importnat consideradtions and we will take them into account now.

The other thing that has changed between the two analyses is of course the DESeq2 version (v1.4.1 vs 1.6.3), which might also influence results slightly, I guess.

Removing the new IDs is not really an option because it not only new ones, but also substitutes and other changes that are updated with new ensembl releases.

We'll try the option beta-prior=false to see what effect this has on our data set.

Thanks again,

Roman