What the heck is going on with this DESeq2 MA Plot?
1
2
Entering edit mode
5.2 years ago

Hello,

I'm getting a very strange MA plot, despite following the DESeq2 instructions. Effectively, my MA plot looks exactly like a line, with x-values ranging from 0 to 50,000, and y ranging from 0 to 90,000 (so, not the y-axis range of -2 to 2 the code predicts). Furthermore, when I look at the data itself, it doesn't look so odd - in fact, log2FoldChange (this is what is plotted on y axis right?) seems to have an even number of points below and above 0, which is not the case in the MA Plot.

What is going on? Any suggestions would be appreciated!

Image here: https://ibb.co/dhwUcv

   dds <- DESeq(myData)
   deseqOutput <- results(dds)
   deseqOutput_p05 <- results(dds, alpha=0.05)
   summary(deseqOutput_p05)

    out of 23719 with nonzero total read count
    adjusted p-value < 0.05
    LFC > 0 (up)     : 339, 1.4% 
    LFC < 0 (down)   : 381, 1.6% 
    outliers [1]     : 0, 0% 
    low counts [2]   : 4143, 17% 
    (mean count < 3)

    sum(deseqOutput$log2FoldChange < 0)
        [1] 13241

    sum(deseqOutput$log2FoldChange > 0)
        [1] 10478

plotMA(deseqOutput , ylim=c(-2,2))
plotMA(deseqOutput, ylim=c(-2,100000), xlim=c(0,100000))

EDIT: I found a clue! If I do a remove objects AND clear plots AND type dev.off() into the console in R, I still get this error. However, the error disappears if I completely close RStudio and open it up again/re-run code...UNLESS I include a piece of code that plots a PCA using ggplot2. Any ideas what this means? Thanks for your help!

EDIT 2: Here is the code from the PCA code that seems to cause this problem:

quickPCA <- function(rawCounts, sampleNames, condition, title){ #sometimes have to change function names to match what you are passing in for some reason
  print('starting')

  # load ggplot2
  library(edgeR)
  library(DESeq2)
  library(ggplot2)
  library(ggrepel)

  sizeFactors <- calcNormFactors(rawCounts, method="RLE")
  rawCounts_normd <- asinh(rawCounts*sizeFactors)
  colnames(rawCounts_normd) <- sampleNames

  # PCA
  cds.pca <- prcomp(t(rawCounts_normd),
                           center = TRUE)

  # create data frame with scores
  scores <- as.data.frame(cds.pca$x)
  print(rownames(scores))

  # plot of first two principal components
  p= ggplot(data = as.data.frame(cds.pca$x), aes(x = PC1, y = PC2, label = rownames( as.data.frame(cds.pca$x) ) )) +
    geom_hline(yintercept = 0, colour = "gray65") +
    geom_vline(xintercept = 0, colour = "gray65") +
    geom_point() +
    geom_text_repel(size = 5, aes(colour = factor(condition), label = rownames( as.data.frame(cds.pca$x)), fontface='bold')) +
    scale_colour_discrete(l=40) +
    ggtitle(title)
  print(p)
}
R MA Plot DESeq2 RNA-Seq • 4.6k views
ADD COMMENT
0
Entering edit mode

UNLESS I include a piece of code that plots a PCA using ggplot2.

This seems crucial. If you add that code to your post this might also provide us with clues about what goes wrong.

ADD REPLY
0
Entering edit mode

Good point! I've included my PCA-plotting code under Edit2 in the original post.

ADD REPLY
9
Entering edit mode
2.5 years ago
Anna ▴ 130

I had the same problem, but not necessarily linked to a previous PCA plot. In my case, the plotMA function was being called from a different package. Be aware that there are many packages that have a plotMA and not all of them take the same class of object.

I solved it by using DESeq2::plotMA(res)

ADD COMMENT

Login before adding your answer.

Traffic: 1194 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6