[alisrpp]

You are using DESeq2 v1.2 and reading the documentation off the Bioc website for version v1.0.

This is why I recommended you to look up the man page for ?results.

And it's always safer to use:

browseVignettes("DESeq2")

to ensure you are reading the documentation for the version you are using.

Ups! Thanks!

By using:

browseVignettes("DESeq2")

i am obtaining

/library/DESeq2/doc/DESeq2.pdf not found

Any idea about why this is happening?

[Michael Love]

You are using DESeq2 v1.2 and reading the documentation off the Bioc website for version v1.0.

This is why I recommended you to look up the man page for ?results.

And it's always safer to use:

browseVignettes("DESeq2")

to ensure you are reading the documentation for the version you are using.

[alisrpp]

The local fit helps you out by dipping down at the high counts, giving you lower estimates of dispersion there than by taking the mean overall. The large circles are genes whose dispersion estimates are not shrunk because they deviate too far from the fitted line. The local fit is useable as well I would say.

I presume -- because you haven't posted the output of sessionInfo() -- that you could be using a version of DESeq2 which performs automatic independent filtering to optimize the number of rejections. Check the Details section of ?results and the argument 'independentFiltering'.

Here is my sessionInfo()

> 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] grid      parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] gplots_2.11.3           MASS_7.3-29             KernSmooth_2.23-10      caTools_1.14           
 [5] gdata_2.13.2            gtools_3.1.0            RColorBrewer_1.0-5      DESeq2_1.2.0           
 [9] RcppArmadillo_0.3.920.1 Rcpp_0.10.5             GenomicRanges_1.14.1    XVector_0.2.0          
[13] IRanges_1.20.0          BiocGenerics_0.8.0     

loaded via a namespace (and not attached):
 [1] annotate_1.40.0      AnnotationDbi_1.24.0 Biobase_2.22.0       bitops_1.0-6        
 [5] DBI_0.2-7            genefilter_1.44.0    lattice_0.20-24      locfit_1.5-9.1      
 [9] RSQLite_0.11.4       splines_3.0.2        stats4_3.0.2         survival_2.37-4     
[13] tools_3.0.2          XML_3.95-0.2         xtable_1.7-1

I'm still having some questions about the independent filtering. I am understanding that I have to perform it right after the:

dds<- estimateSizeFactors(dds)
ddsLocal <- estimateDispersions(dds, fitType="local")
ddsLocal <- nbinomWaldTest(ddsLocal)
res <- results(ddsLocal, name= "condition_MR_vs_BR") 
res <- res[order(res$padj),]

With the code:

use <- res$baseMean >= 10 & !is.na(res$pvalue)
resFilt <- res [use, ]
resFilt$padj <- p.adjust(resFilt$pvalue, method="BH")

And then continue with (for example):

#filter for upregulated and downregulated genes
resSig <- resFilt[ resFilt$padj < .1, ]
write.table(resSig[ order( resSig$log2FoldChange, -resSig$baseMean ), ]
            ,"~\\DESeq2\\BRvsMR_old_DownRegulated.txt")
write.table(resSig[ order( -resSig$log2FoldChange, -resSig$baseMean ),
                    ],"~\\DESeq2\\BRvsMR_old_UpRegulated.txt")

But i also found this code:

filterThreshold <- 2.0
keep <- rowMeans ( counts(ddsLocal, normalized=TRUE)) > filterThreshold

And now i'm really confused. What is this for? (the last code) When i have to use it? Am I doing something wrong?

Thank you again for the help.

[sindrle]

Ah!
Always a good answer from you. Thanks!

[dpryan]

If you're trying to do a comparison, the edgeR equivalent would be the "RLE" method in calcNormFactors(). estimateSizeFactors() doesn't take an option in this regard (see help(estimateSizeFactors) for what options it does have).

[sindrle]

Ok, but typing in "estimateSizeFactors()" to look for options gives no info just (object, ...)

But what is the default? Im curious because Im running edgeR with "upper quartile".

Thanks!

[Michael Love]

estimateSizeFactors() is the same function used for normalization as in DESeq, so the method is the same median ratio method as in DESeq. I can try to make this more clear in the man page for DESeq().

[sindrle]

Hi!
Quick question; what is the default normalisation method in DESeq2?
Cant find anything in the manual.

Thanks!

[Michael Love]

The local fit helps you out by dipping down at the high counts, giving you lower estimates of dispersion there than by taking the mean overall. The large circles are genes whose dispersion estimates are not shrunk because they deviate too far from the fitted line. The local fit is useable as well I would say.

I presume -- because you haven't posted the output of sessionInfo() -- that you could be using a version of DESeq2 which performs automatic independent filtering to optimize the number of rejections. Check the Details section of ?results and the argument 'independentFiltering'.

[alisrpp]

There is a lot of shrinkage of dispersion estimates (the blue points are close to the red line) presumably because you have few degrees of freedom (number of samples minus number of parameters to estimate). Yes, I would use fitType="mean" for this dataset as there is not a clear trend.

Hi Michael,

First of all, thank you for all the help.

I tried fitType="mean" with the code:

#Count matrix input
Cele_SPvsLR = read.csv (file.choose(), header=TRUE, row.names=1)
CeleDesign <- data.frame(
  row.names = colnames(Cele_SPvsLR),
  condition = factor(c("SP", "SP", "LR", "LR")))
dds <- DESeqDataSetFromMatrix(countData = Cele_SPvsLR,
                              colData = CeleDesign,
                              design = ~ condition)
dds


#Est size factor = normalize for library size
dds<- estimateSizeFactors(dds)
ddsMean <- estimateDispersions(dds, fitType="mean")
ddsMean <- nbinomWaldTest(ddsMean)
#For a parametric dispersion
dds <- DESeq(dds)
#plot dispersion
plotDispEsts(ddsMean)
plotDispEsts(dds)

#Differential expression analysis
resultsNames(ddsMean)
res <- results(ddsMean, name= "condition_SP_vs_LR") 
res <- res[order(res$padj),]
head(res)
plotMA(ddsMean)
sum(res$padj < .1, na.rm=TRUE)

This resulted in:

> sum(res$padj < .1, na.rm=TRUE)
[1] 0

But when I tried fitType="local" I get:

> sum(res$padj < .1, na.rm=TRUE)
[1] 337

(*see attached plots)

The other "weird" (at least for me) thing is that by filtering by overall count to my fitType="local" results I get less genes:

> plot(res$baseMean, pmin(-log10(res$pvalue), 50),
+ log="x", xlab="mean of normalized counts",
+ ylab=expression(-log[10](pvalue)))
Warning message:
In xy.coords(x, y, xlabel, ylabel, log) :
  19407 x values <= 0 omitted from logarithmic plot
> abline(v=10, col="red", lwd=1)
> use <- res$baseMean >= 10 & !is.na(res$pvalue)
> table(use)
use
 FALSE   TRUE 
 93445 105376 
> resFilt <- res [use, ]
> resFilt$padj <- p.adjust(resFilt$pvalue, method="BH")
> sum(resFilt$padj < .1, na.rm=TRUE)
[1] 152

Just as a refresher, I'm working with 4 samples (2 replicates/reproductive stages).

Do you think that it is a bad idea to use fitType="local" given that it is the fitType that gives me results?

Can you help me understand why by filtering my "local" results I am obtaining less genes?

Thanks, thanks, and thanks.

[Michael Love]

There is a lot of shrinkage of dispersion estimates (the blue points are close to the red line) presumably because you have few degrees of freedom (number of samples minus number of parameters to estimate). Yes, I would use fitType="mean" for this dataset as there is not a clear trend.

[alisrpp]

As there is not a clear dependence of the dispersion estimates on the mean of normalized counts, I would recommend to use fitType="mean".

Hi,

By doing what you suggested:

dds<- estimateSizeFactors(dds)
ddsMean <- estimateDispersions(dds, fitType="mean")
ddsMean <- nbinomWaldTest(ddsMean)
plotDispEsts(ddsMean)

I'm obtaining the subsequent plot.
It's really weird, right?
Should i continue with fitType="mean"?

Thanks.

[Michael Love]

Then these might just be difficult rows for some reason.

You can examine them with:

head(counts(dds[!mcols(dds)$betaConv,]))

You can exclude them from the results table with:

res <- results(dds[mcols(dds)$betaConv,])

[sindrle]

I upgraded all as you said, but still 159 not fitted at maxit 100.

[Michael Love]

This dispersion plot looks good/typical.

re:maxit: Unlike with linear regression, there is not a closed form solution for generalized linear models (the model used by DESeq2). Instead there is an iterative algorithm to find the optimal betas (the log fold changes). 'maxit' is the maximum number of iterations to run, and if the convergence criteria has not been met then the software prints a warnings and tells you how to find these genes for which this was the case.

However, I believe most datasets would have all genes converge if you used DESeq2 v1.2 which was just released.

This involves upgrading R to 3.0.2 and Bioconductor to v2.13. Information here:

http://www.bioconductor.org/install/