In tschauer/HelpersforDESeq2: This is a collection of helper functions for DESeq2

knitr::opts_chunk$set(echo = TRUE)

Load Example Results

library(HelpersforDESeq2)

data_dir <- system.file("extdata/", package = "HelpersforDESeq2")

res1 <- read.delim(paste0(data_dir, "res.HighMLEpATP-HighMLEmATP.txt"), row.names = 1, stringsAsFactors = F)
res2 <- read.delim(paste0(data_dir, "res.LowMLEpATP-LowMLEmATP.txt"), row.names = 1, stringsAsFactors = F)

res1$log10baseMean <- log10(res1$baseMean+1)
res2$log10baseMean <- log10(res2$baseMean+1)

head(res1)
head(res2)

Goal

par(mfrow=c(1,2), mar = c(4,4,2,2), oma = c(1,1,1,1), mgp = c(2,1,0))

plot(res1$log10baseMean,
     res1$log2FoldChange,
     main = "ugly?")

plottingMA(res =  res1,
           main_title = "better?",
           selection_ids = c("roX1","roX2"),
           selection_id_type = "gene_symbol",
           selection_point_size = 1,
           selection_text_label = TRUE,
           selection_shadow = FALSE,
           xlims = c(0, 6),
           ylims = c(-10,10),
           x_axis_by = 2,
           padj_cutoff = 0.01,
           show_legend = TRUE)

Build Plot Step-by-Step

par(mfrow=c(1,1), mar = c(4,4,2,2), oma = c(4,4,4,4), mgp = c(2,1,0))

plot(res1$log10baseMean,
     res1$log2FoldChange)

par(mfrow=c(1,1), mar = c(4,4,2,2), oma = c(4,4,4,4), mgp = c(2,1,0))

plot(res1$log10baseMean,
     res1$log2FoldChange,
     xlab = "log10 mean counts",
     ylab = "log2 fold change",
     xlim = c(0,6),
     ylim = c(-10,10),
     col = rgb(0.7,0.7,0.7,0.5), pch=19, cex = 0.25)

abline(h=0, col="grey32")

par(mfrow=c(1,1), mar = c(4,4,2,2), oma = c(4,4,4,4), mgp = c(2,1,0))

plot(res1$log10baseMean,
     res1$log2FoldChange,
     xlab = "log10 mean counts",
     ylab = "log2 fold change",
     xlim = c(0,6),
     ylim = c(-10,10),
     col = rgb(0.7,0.7,0.7,0.5), pch=19, cex = 0.25)

abline(h=0, col="grey32")

selection_ids <- c("roX1","roX2")
selection_vector <- res1$gene_symbol %in% selection_ids
selection_color <- rgb(0.9,0.6,0,1)

points(res1$log10baseMean[selection_vector],
       res1$log2FoldChange[selection_vector],
       col = selection_color, pch=19, cex = 1.0)

legend("topright", legend =  "labeled", col = selection_color,
       bg = "white", border = NA, bty = "n", cex = 0.8, pch = 19)

par(mfrow=c(1,1), mar = c(4,4,2,2), oma = c(4,4,4,4), mgp = c(2,1,0))

plot(res1$log10baseMean,
     res1$log2FoldChange,
     xlab = "log10 mean counts",
     ylab = "log2 fold change",
     xlim = c(0,6),
     ylim = c(-10,10),
     col = rgb(0.7,0.7,0.7,0.5), pch=19, cex = 0.25)

abline(h=0, col="grey32")

selection_ids <- c("roX1","roX2")
selection_vector <- res1$gene_symbol %in% selection_ids
selection_color <- rgb(0.8,0,0,1)

points(res1$log10baseMean[selection_vector],
       res1$log2FoldChange[selection_vector],
       col = selection_color, pch=19, cex = 1.0)

text(res1$log10baseMean[selection_vector],
     res1$log2FoldChange[selection_vector], 
     res1$gene_symbol[selection_vector],
     col = selection_color, adj = c(0,-0.5))

legend("topright", legend =  "labeled", col = selection_color,
       bg = "white", border = NA, bty = "n", cex = 0.8, pch = 19)

par(mfrow=c(1,1), mar = c(4,4,2,2), oma = c(4,4,4,4), mgp = c(2,1,0))

plot(res1$log10baseMean,
     res1$log2FoldChange,
     xlab = "log10 mean counts",
     ylab = "log2 fold change",
     xlim = c(0,6),
     ylim = c(-10,10),
     col = rgb(0.7,0.7,0.7,0.5), pch=19, cex = 0.25)

abline(h=0, col="grey32")

selection_vector <- res1$padj < 0.01

points(res1$log10baseMean[selection_vector],
       res1$log2FoldChange[selection_vector],
       col = selection_color, pch=19, cex = 0.25)

legend("topright", legend =  "significant", col = selection_color,
       bg = "white", border = NA, bty = "n", cex = 0.8, pch = 19)

Write a function

par(mfrow=c(1,2), mar = c(4,4,2,2), oma = c(1,1,1,1), mgp = c(2,1,0))

plottingMAbasics <- function(res,
                             padj_cutoff = 0.01){

        plot(res$log10baseMean,
             res$log2FoldChange,
             xlab = "log10 mean counts",
             ylab = "log2 fold change",
             xlim = c(0,6),
             ylim = c(-10,10),
             col = rgb(0.7,0.7,0.7,0.5), pch=19, cex = 0.25)

        abline(h=0, col="grey32")

        selection_vector <- res$padj < padj_cutoff

        points(res$log10baseMean[selection_vector],
               res$log2FoldChange[selection_vector],
               col = selection_color, pch=19, cex = 0.25)

        legend("topright", legend =  "significant", col = selection_color,
               bg = "white", border = NA, bty = "n", cex = 0.8, pch = 19) 
}

plottingMAbasics(res1)
plottingMAbasics(res2)

Final MA plot Function

par(mfrow=c(1,1), mar = c(4,4,2,2), oma = c(4,4,4,4), mgp = c(2,1,0))

res_file_name <- "res.HighMLEpATP-HighMLEmATP.txt"
res_tmp <- read.delim(paste0(data_dir, res_file_name), row.names = 1, stringsAsFactors = F)

res_name <- gsub(".txt","", res_file_name)
assign(res_name, res_tmp)

plottingMA(res =  get(res_name),
           main_title = gsub("res.","",res_name),
           selection_ids = c("roX1","roX2"),
           selection_id_type = "gene_symbol",
           selection_point_size = 1,
           selection_text_label = TRUE,
           selection_shadow = FALSE,
           xlims = c(0, 6),
           ylims = c(-10,10),
           x_axis_by = 2,
           padj_cutoff = 0.01,
           show_legend = TRUE)

Log2FC - log2FC plot

par(mfrow=c(1,1), mar = c(4,4,2,2), oma = c(4,4,4,4), mgp = c(2,1,0))

res_merged <- merge(res1, res2, by = "row.names")

plot(res_merged$log2FoldChange.x,
     res_merged$log2FoldChange.y,
     xlab = "log2FC res1",
     ylab = "log2FC res2",
     xlim = c(-6,6),
     ylim = c(-6,6),
     col = rgb(0.7,0.7,0.7,0.5), pch=19, cex = 0.25)

abline(h=0, v=0, col="grey32")
abline(coef = c(0,1), col="grey32", lty=2)

selection_ids = c("roX1","roX2")
selection_vector <- res_merged$gene_symbol.x %in% selection_ids
selection_color <- rgb(0.8,0,0,1)

points(res_merged$log2FoldChange.x[selection_vector],
       res_merged$log2FoldChange.y[selection_vector],
       col = selection_color, pch=19, cex = 1.0)

text(res_merged$log2FoldChange.x[selection_vector],
     res_merged$log2FoldChange.y[selection_vector], 
     res_merged$gene_symbol.x[selection_vector],
     col = selection_color, adj = c(0,-0.5))

legend("topleft", legend =  "labeled", col = selection_color,
       bg = "white", cex = 1, pch = 19)

Final log2FC plot Function

par(mfrow=c(1,1), mar = c(5,5,1,1), oma = c(4,4,4,4), mgp = c(3,1,0))

res_file_names <- c("res.HighMLEpATP-HighMLEmATP.txt", "res.LowMLEpATP-LowMLEmATP.txt")

for(res_file_name in res_file_names){

        res_tmp <- read.delim(paste0(data_dir, res_file_name), row.names = 1, stringsAsFactors = F)

        res_name <- gsub(".txt","", res_file_name)
        assign(res_name, res_tmp)

}

res_names <- ls(pattern = "^res\\.")
res_names

plotLog2FC(res1 = get(res_names[1]),
           res2 = get(res_names[2]),
           main_title = "",
           x_label = paste("log2FC \n", gsub("res.","",res_names[1])),
           y_label = paste("log2FC \n", gsub("res.","",res_names[2])),
           lims = c(-6,6),
           point_size = 0.25,
           selection_ids = c("roX1","roX2"),
           selection_id_type = "gene_symbol",
           selection_point_size = 1,
           selection_legend = "labeled",
           selection_text_label = TRUE)