DESeq2 without biol replicates

[bioBob]

Is there an implicit normalization across the samples in this?

[Simon Anders]

You are comparing DESeq(1) and DESeq2. DESeq(1) does not have and never had any shrinkage on fold changes.

[Michael Love]

Yes, shrunken log fold changes were incorporated for the Wald test since the beginning of DESeq2.

You can always check the methods of the current version with ?DESeq, vignette("DESeq2"), or for more specific updates: news(package="DESeq2").

[visserm]

Hi Michael,
Is the shrunken log fold change also applied in DESeq2 v 1.0.17?

[Michael Love]

hi,

I made a change to rlogTransformation in the development branch of DESeq2 (v 1.1.21), which might be relevant for your analysis.

Previously, the rlog transformation by default would use a blind design formula and use the MAP dispersion estimates (the blue points in a plot produced by plotDispEsts). True, large log fold changes would have high MAP dispersions and be shrunk toward 0 (especially with no replicates). Instead we now (in the devel branch) use the fitted dispersion estimates, which shrinks log fold changes based on the general mean-dispersion trend line, and preserves true large log fold changes.

It is described here:

http://www.bioconductor.org/packages/2.13/bioc/news/DESeq2/NEWS

[sisterdot]

hey,

could you please help again, i am a bit as it seems that the DESeq and DESeq2 regularized LogFC calculations differ strongly if i stick to the recommended procedure

looked for a potentially differentially expressed gene from pasilla (FBgn0039155 counts are 38 to 831), wanted to see what the DESeq2 regularized LogFC would be and tried to compared it to the DESeq rLogFC as depicted below.
DESeq says rLogFC is 3.095034, DESeq2 says 0.1595547 (and the later doesn't change much if i try out blind rlogTransformation, estimateSizeFactors, estimateDispersions). somehow i believe the DESeq2 procedure i used isn't doing what i want.

thanks!

----------------- deseq2

library(DESeq2)
library(pasilla)
data("pasillaGenes")
countData <- counts(pasillaGenes)
countData<-countData[,c("treated1fb","untreated1fb")]
head(countData)
colData <- pData(pasillaGenes)[c("treated1fb","untreated1fb"),c("condition","type")]
dds <- DESeqDataSetFromMatrix(
       countData = countData,
       colData = colData,
       design = ~ condition)
rld <- rlogTransformation( dds)#,blind = TRUE)
res <- data.frame(
   assay(rld), 
   avgLogExpr = ( assay(rld)[,2] + assay(rld)[,1] ) / 2,
   rLogFC = assay(rld)[,2] - assay(rld)[,1] )
res[c("FBgn0039155"),]

treated1fb untreated1fb avgLogExpr rLogFC
FBgn0039155 8.850547 9.010102 8.930324 0.1595547

----------------- deseq

library(DESeq)
library(pasilla)
data("pasillaGenes")
countData <- counts(pasillaGenes)
countData<-countData[,c("treated1fb","untreated1fb")]
condition <- factor( c( "treated", "untreated") )
cds <- newCountDataSet( countData, condition )
cds <- estimateSizeFactors( cds )
cds <- estimateDispersions( cds, method="blind", sharingMode="fit-only" )
vsd <- varianceStabilizingTransformation( cds )
res <- nbinomTest( cds, "treated", "untreated") 
rlfc <- (exprs(vsd)[, c("untreated1fb")]  -exprs(vsd)[, c("treated1fb")] )
res$rLogFC<-rlfc
res[res$id=="FBgn0039155",]

id baseMean baseMeanA baseMeanB foldChange log2FoldChange
9696 FBgn0039155 513.8001 31.70826 995.8919 31.40796 4.973058
pval padj rLogFC
9696 9.669445e-23 1.062188e-18 3.095034

[sisterdot]

thanks for your help Simon- you are the best!

> So, if you promise to not use any such comparisons for actual science
hihihi...
you are totally right of course!

[Simon Anders]

To be honest, we couldn't yet be bothered to explain how to analyse such data in DESeq2. It's tricky to write up because too many people will misinterpret whatever I write as if it were actually possible to conduct a meaningful statistical analysis when comparing just two samples.

So, if you promise to not use any such comparisons for actual science, here is how you do it:

Start as above:

library(DESeq2)
library(pasilla)
data("pasillaGenes")
countData <- counts(pasillaGenes)
countData<-countData[,c("treated1fb","untreated1fb")]
colData <- pData(pasillaGenes)[c("treated1fb","untreated1fb"),c("condition","type")]
dds <- DESeqDataSetFromMatrix(
       countData = countData,
       colData = colData,
       design = ~ condition)

Now use DESeq2's new "rlog transformation". This replaced the VST we had before. It transforms the average of the genes across samples to a log2 scale but "pulls in" those genes for which the evidence for strong fold changes is weak due to low counts.

rld <- rlogTransformation( dds )

As this is a logarithm-like scale, the differences between the samples can be considered as a "regularized log2 fold change". Let's make a result data frame:

res <- data.frame(
   assay(rld), 
   avgLogExpr = ( assay(rld)[,2] + assay(rld)[,1] ) / 2,
   rLogFC = assay(rld)[,2] - assay(rld)[,1] )

And now we have a ranking of genes by regularized fold changes:

> head( res[ order(res$rLogFC), ] )
            treated1fb untreated1fb avgLogExpr    rLogFC
FBgn0011260   7.830359     6.627326   7.228842 -1.203033
FBgn0001226  10.128636     8.929985   9.529311 -1.198652
FBgn0034718   8.503006     7.315640   7.909323 -1.187366
FBgn0003501   7.927864     6.743974   7.335919 -1.183889
FBgn0033635  11.126300     9.973979  10.550139 -1.152321
FBgn0033367  13.411814    12.269436  12.840625 -1.142378

This ranking put ones genes on top which are strongly downregulated in the second sample compared to the first one. If you do this with normal log expressions, the weak genes will appear at the top because they are noisiest and hence tend to have exaggerated fold changes.

The advantage of this procedure is that it does not produce any p values (which would be misleading anyway).

[sisterdot]

Hey all,

i wanted to try DESeq2 in a pairwise comparison setting without replicates.

unfortunately using the code below on a pasilla-test-case i can only find differentially expressed genes using DESeq and not DESeq2.
(DESeq2 says "same number of samples and coefficients to fit, estimating dispersion by treating samples as replicates" so should be getting the task)

i know it must be me, but where am I going wrong?

thanks a lot!!


#-------------------------- deseq2

library(DESeq2)
library(pasilla)
data("pasillaGenes")
countData <- counts(pasillaGenes)
countData<-countData[,c("treated1fb","untreated1fb")]
colData <- pData(pasillaGenes)[c("treated1fb","untreated1fb"),c("condition","type")]
dds <- DESeqDataSetFromMatrix(
       countData = countData,
       colData = colData,
       design = ~ condition)
colData(dds)$condition <- factor(colData(dds)$condition,levels=c("treated","untreated"))
dds
dds <- DESeq(dds)
res <- results(dds)
res <- res[order(res$padj),]
head(res)

DataFrame with 6 rows and 5 columns
baseMean log2FoldChange lfcSE pvalue padj
<numeric> <numeric> <numeric> <numeric> <numeric>
FBgn0000008 60.7056967 -0.09993564 0.36063169 0.7816935 0.9973207
FBgn0000014 0.8344279 -0.03920709 0.08405437 0.6408940 0.9973207
FBgn0000015 0.4172140 -0.02034815 0.06095374 0.7385084 0.9973207

#-------------------------- deseq

library(DESeq)
library(pasilla)
data("pasillaGenes")
countData <- counts(pasillaGenes)
countData<-countData[,c("treated1fb","untreated1fb")]
condition <- factor( c( "treated", "untreated") )
cds <- newCountDataSet( countData, condition )
cds <- estimateSizeFactors( cds )
cds <- estimateDispersions( cds, method="blind", sharingMode="fit-only" )
res <- nbinomTest( cds, "treated", "untreated") 
res <- res[order(res$padj),]
head(res)

id baseMean baseMeanA baseMeanB foldChange log2FoldChange
9696 FBgn0039155 513.8001 31.70826 995.89187 31.40796187 4.973058
2605 FBgn0026562 20645.8226 6151.40275 35140.24243 5.71255758 2.514137
29 FBgn0000071 241.5499 444.75009 38.34963 0.08622736 -3.535710
pval padj
9696 9.669445e-23 1.062188e-18
2605 1.477349e-12 8.114339e-09
29 3.097019e-10 1.134025e-06