SEQanswers

Go Back   SEQanswers > Applications Forums > RNA Sequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
infinite fold change RNAseq biofreak Bioinformatics 6 07-02-2012 10:02 AM
DESeq DiffExpress Vs. Fold Change KellerMac Bioinformatics 3 06-10-2011 10:49 AM
fold change value-cuffdiff madsaan Bioinformatics 4 02-10-2011 06:51 AM
fold change for genes with 0 reads rnaseq RNA Sequencing 15 11-05-2010 01:23 AM
fold change for genes with 0 reads rnaseq General 2 11-04-2010 07:44 PM

Reply
 
Thread Tools
Old 08-20-2013, 05:15 AM   #1
Alex234
Member
 
Location: UK

Join Date: Aug 2013
Posts: 31
Default Contradictory log-fold-change figures from DESeq?

Hello, I have made a couple of figures in DESeq, the first should plot genes with adjusted p-value <0.05 in red with >0.05 genes in grey, and was made with:

plotMA(res,
+ col = ifelse(res$padj>0.05, "gray69", "red3"),
+ linecol = "#ff000080",
+ xlab = "mean of normalized counts", ylab = expression(log[2]~fold~change),
+ log = "x", cex=0.45)
> dev.off()


The second should just plot genes that have a p-value >0.05 (i.e. all the red ones from the first figure):

resSig = res[res$padj <0.05, ]
> plotMA(resSig)


And yet there are red points in the first (padj<0.05) that are not present in the second! Can anybody explain this to me? I have attached the two figures.

And on a related note, is it better to use the p-value or adjusted p-value that DESeq gives you?

Many thanks

Alex
Attached Files
File Type: pdf Fold change padj less than 0.05 in red.pdf (134.4 KB, 33 views)
File Type: pdf padj less than 0.05 only.pdf (5.8 KB, 17 views)
Alex234 is offline   Reply With Quote
Old 08-20-2013, 06:04 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Always use the adjusted p-value.

The differences look to be simply due to the different scales on the graphs and perhaps the fact that you're using a slightly different coloring/inclusion threshold in each.

If you instead did something like:
Code:
plotDE <- function(res)
    plot(res$baseMean, res$log2FoldChange, col = ifelse(res$padj>=0.05, "gray69", "red3"), log="x", pch=6, xlim=c(0.1,1e5), ylim=c(-4,4)))
plotDE(res)
resSig <- res[which(res$padj<0.05),]
plotDE(resSig)
then the two graphs should be more or less identical. The plotMA function is just a nicer version of the "plotDE" function I pasted above (in fact, if you just type "plotMA", without quotes, you'll see exactly what the function does).
dpryan is offline   Reply With Quote
Old 08-20-2013, 06:36 AM   #3
Alex234
Member
 
Location: UK

Join Date: Aug 2013
Posts: 31
Default

Thanks, I'm still not sure what was wrong with plotMA, but the plotDE line you gave me works fine!


Thanks

Alex
Alex234 is offline   Reply With Quote
Old 08-21-2013, 12:08 AM   #4
Alex234
Member
 
Location: UK

Join Date: Aug 2013
Posts: 31
Default

Hi, got a related problem now where the scatterplot of ordinary vs moderated log ratios in DESeq seems to be wrong - the points seem to be upside down (picture attached and lines of code below)

Would be grateful for any advice! (and do you know how I can retrieve the variance stabilised data into an excel/txt/csv file?)

Thanks

Alex

> mod_lfc = (rowMeans( exprs(vsd)[, conditions(cds)=="WT", drop=FALSE] ) -
+ rowMeans( exprs(vsd)[, conditions(cds)=="KO", drop=FALSE] ))

> lfc = res$log2FoldChange
> table(lfc[!is.finite(lfc)], useNA="always")

> logdecade = 1 + round( log10( 1+rowMeans(counts(cdsBlind, normalized=TRUE)) ) )
> lfccol = colorRampPalette( c( "gray", "blue" ) )(6)[logdecade]

> ymax = 4.5
> pdf("Ordinary vs Moderated Log ratios.pdf")
> plot( pmax(-ymax, pmin(ymax, lfc)), mod_lfc,
+ xlab = "ordinary log-ratio", ylab = "moderated log-ratio",
+ cex=0.45, asp=1, col = lfccol,
+ pch = ifelse(lfc<(-ymax), 60, ifelse(lfc>ymax, 62, 16)))
> abline( a=0, b=1, col="red3")
Attached Files
File Type: pdf Ordinary vs Moderated Log ratios.pdf (180.0 KB, 14 views)
Alex234 is offline   Reply With Quote
Old 08-21-2013, 12:36 AM   #5
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Regarding the plot, I would check that res$log2FoldChange is log2(WT/KO) rather than log2(KO/WT), as it appears that there's just a swapped axis somewhere. Regarding writing the vsd data, the exprs() function returns a matrix, so you can just write.csv or write.delim to write to a file.
dpryan is offline   Reply With Quote
Old 08-21-2013, 01:28 AM   #6
Alex234
Member
 
Location: UK

Join Date: Aug 2013
Posts: 31
Default

Great, both worked - re: the vsd data (the numbers are counts?): can/should this be used in the same way as the original counts data to look at DEGs (since the problem of 0 counts seems to have been dealt with), or is it just a form that the data needs to be in to produce PCAs and sample-to-sample difference heatmaps?


Thanks, I really appreciate your help!

Alex
Alex234 is offline   Reply With Quote
Old 08-21-2013, 04:03 AM   #7
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Well, variance stabilized numbers are transformed counts. If you're interested in doing the DEG analysis in DESeq (or edgeR or other similar packages), then you'll find these sorts of transformations most useful for QC, such as PCA and heatmaps. That doesn't mean that you can't use variance stabilized data for DEG analysis. Voom, which is used to transform RNAseq data for use in analysis by limma, is effectively doing something like this, i.e., it tries to model the mean-variance relationship and then transforms the data such that it should have a more normal distribution such that the regular microarray tools work. In practice, both approaches tend to give pretty similar results, at least in my hands.
dpryan is offline   Reply With Quote
Old 08-22-2013, 12:37 AM   #8
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

Hi,

I've seen this thread only yesterday, and now had a close look at the problem described in the original post. In short, this is a bug: The 'col' argument of 'plotMA' is broken and messes up the assignment of colours to points. (Internally, we subset the 'res' object to get rid of genes with all-zero counts and we forgot to subset 'col', too.) I'll fix this later today.

Simon
Simon Anders is offline   Reply With Quote
Old 08-22-2013, 02:56 AM   #9
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

I fixed the issue with the 'col' argument in DESeq 1.12.1. In DESeq2, it was already correct.
Simon Anders is offline   Reply With Quote
Old 08-22-2013, 03:31 AM   #10
Alex234
Member
 
Location: UK

Join Date: Aug 2013
Posts: 31
Default

Thanks - if I re-load DeSEQ 1 from bioconductor, will it update the copy on my system?

Is there any advantage of plotMA over plotDE anyway?


Thanks

Alex
Alex234 is offline   Reply With Quote
Old 08-22-2013, 03:34 AM   #11
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

You'll have to wait a day or two for the change to get propagated to the web server (or you download from the SVN server, if you know how to do this).

What is plotDE?
Simon Anders is offline   Reply With Quote
Old 08-22-2013, 03:36 AM   #12
Alex234
Member
 
Location: UK

Join Date: Aug 2013
Posts: 31
Default

I thought it was part of DESeq - dpryan introduced it to me:

Quote:
Originally Posted by dpryan View Post
Always use the adjusted p-value.

The differences look to be simply due to the different scales on the graphs and perhaps the fact that you're using a slightly different coloring/inclusion threshold in each.

If you instead did something like:
Code:
plotDE <- function(res)
    plot(res$baseMean, res$log2FoldChange, col = ifelse(res$padj>=0.05, "gray69", "red3"), log="x", pch=6, xlim=c(0.1,1e5), ylim=c(-4,4)))
plotDE(res)
resSig <- res[which(res$padj<0.05),]
plotDE(resSig)
then the two graphs should be more or less identical. The plotMA function is just a nicer version of the "plotDE" function I pasted above (in fact, if you just type "plotMA", without quotes, you'll see exactly what the function does).
Alex234 is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 12:18 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO