Change log2FoldChange range - plotMA

Question

[Johan Largo]

Hi all

I am trying to make a graph with the plotMA tool and it throws me the graph below, now I would like from the script in RStudio to color those that have a log2FoldChange greater than 2.

Try to use

plotMA(res,
       main="MA plot",
       cex=0.5, lwd=0.5,colNonSig = "#595959",colSig = "#00BFC4",colLine = "#F8766D",col=ifelse(res$log2FoldChange>2),
       ylim=c(-16,16),
       alpha=0.001)

But it does not work. Thanks :D

[ATpoint]

Just set pvalues for all genes to 1 that have logFCs below your cutoff. That way they do not show up as significant. For anything else you probably need some custom code such as ggplot.

[Johan Largo]

mmm

You're right about "Just set pvalues" and my coconut lit up.

I think that what I want to do can basically be done by another graph, for example a VolcanoPlot. I've been reading here http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html and basically "In DESeq2, the function plotMA shows the log2 fold changes attributable to a given variable over the mean of normalized counts for all the samples in the DESeqDataSet. Points will be colored red if the adjusted p value is less than 0.1"

As I understand it, plotMA only reports which genes have "significant" p-values and gives me an idea of the log2FoldChange and normalized counts from my assay. But if you wanted to "colorize" for two conditions, it would be better to use a VolcanoPlot, right? correct me if I'm crazy?

Thanks for your response. :D

[ATpoint]

These are two different types of plots. Use what you find appropriate. MA is baseMean vs logFC and Volcano is logFC vs -log10(p)

[Johan Largo]

Oh, okay I understand, ok thanks :D

Answer 1

[seidel]

You can use base R graphics to make these plots. The data is sitting there in columns of the res object, so you can filter it directly, and use boolean vectors to pick out the things you need:

# make sure there are no NA values
sum(is.na(res$log2FoldChange))

# choose some significance threshold
# (i use the name iv to remind myself it's an index vector)
iv.sig <- res$padj < 0.001 

# genes up > 1
iv.up <- res$log2FoldChange > 1 & iv.sig

# genes down  < -1
iv.dn <- res$log2FoldChange < -1 & iv.sig

# make an MA plot
plot(log2(res$baseMean + 1), res$log2FoldChange, pch=".", col="grey",
     main="MA plot", xlab="log2(baseMean)", ylab="Log2FC")
points(log2(res$baseMean + 1)[iv.up], res$log2FoldChange[iv.up], col="red", pch=20)
points(log2(res$baseMean + 1)[iv.dn], res$log2FoldChange[iv.dn], col="green", pch=20)

# or make a Volcano plot
plot(res$log2FoldChange, abs(res$stat), pch=".", col="grey",
     main="Volcano Plot", xlab="Log2FC", ylab="Sig Stat")
points(res$log2FoldChange[iv.up], abs(res$stat)[iv.up], col="red", pch=20)
points(res$log2FoldChange[iv.dn], abs(res$stat)[iv.dn], col="green", pch=20)

# get the genes of interest, up or down
goi <- rownames(res)[iv.up | iv.dn]