Select what conditionsto graph in Report tools with Deseq containing multiple conditions
Entering edit mode
Last seen 6.1 years ago


I have been trying to setup a HTML report of DESEQ2 analysis using ReportingTools and have encountered some snags.

My DESEQ2 analysis has 5 conditions and I want to create a report for each contrast. After reading some posts on this site I was able to get that to happen. Basically I can get logfold change of contrast I want but the graphs created by ReportingTools include all conditions. 

If possible I would like to show only  box plots of the two conditions that are in my contrast and  in order C(my control sample), then CE (treatment sample):

## annotate using biomaRt
add.anns <-  function(df, mart, ...) {
  nm <- rownames(df)
  anns <- getBM( attributes = c("ensembl_gene_id", "mgi_symbol", "description"), filters = "ensembl_gene_id", values = nm, mart = mart)
  anns <- anns[match(nm, anns[, 1]), ]
  colnames(anns) <- c("Ensembl", "Symbol", "Gene Description")
  df <- cbind(anns, df[, 1:ncol(df)]) 
  rownames(df) <- nm

# to add links to ensembl
ensemblLinks <-
  function(df, ...){
    naind <-$Ensembl)
    df$Ensembl <- hwrite(as.character(df$Ensembl), link = paste0("",
                                                                 as.character(df$Ensembl)), table = FALSE)
    df$Ensembl[naind] <- ""

# to add links to Wikigenes
wikiLinks <-
  function(df, ...){
    naind <-$Symbol)
    df$Symbol <- hwrite(as.character(df$Symbol), link = paste0("",
                                                               as.character(df$Symbol)), table = FALSE)
    df$Symbol[naind] <- ""
#More colorful graphs
theme <- reporting.theme.alternate()
lattice.options(default.theme = theme)
#Publish report
des2Report <- HTMLReport(shortName = "C vs CE RNASeq Analysis",title = "C vs CE RNA-Ssq Analysis of differential expression using DESeq2", reportDirectory = "./reports")
publish(dds.dse, des2Report,  .modifyDF=list(add.anns, modifyReportDF, ensemblLinks, wikiLinks), mart = m, pvalueCutoff=0.05, contrast = c("condition","CE","C"), resultsName = "CvsCE", factor = colData(dds.dse)$condition, reportDir="./reports")


I think it has something to do with:

factor = colData(dds.dse)$condition

But nothing I have tried works.

Any help/hints would be greatly appreciated.

Also is it possible to alter the plots to be made with box plot of ggplot2 instead?

I want to make something like this (Without the other conditions as stated above):

I think ReportingTools calls the following function  from "plotting_functions.R", found in ReportingTools source, to make a graph:

.make.gene.plots <- function(df, expression.dat, factor,,
    ylab.type = "Expression Value", scales = list(), par.settings = list(),
    xlab = NULL, ...){
    scales <- c(scales, list(x = list(rot = 45)))
    if(is(expression.dat, "CountDataSet")){
        ## Get the normalized counts, but add a pseudocount to all of the
        ## entries, because we're going to plot on a log-scale
        expression.dat <- counts(expression.dat, normalized=TRUE)+1
        scales <- c(scales, list(y = list(log = 10)))
    } else if(inherits(expression.dat, "eSet")){
        expression.dat <- exprs(expression.dat)
    } else if(is(expression.dat, "DGEList")){
        ## If we get a DGEList, get the normalized cpm value. Again, we add a
        ## pseudocount to avoid problems with zero count data
        expression.dat <- cpm(expression.dat$counts) + 1
    } else if(is(expression.dat, "data.frame")){
        ## If it's a data.frame, try to coerce it to a matrix,
        ## but this might not work.
        expression.dat <- as.matrix(expression.dat)
    } else if(is(expression.dat, "DESeqDataSet")){
        expression.dat <- counts(expression.dat, normalized = TRUE) + 1
        scales <- c(scales, list(y = list(log = 10)))
    } else if(is.null(expression.dat)) {
        stop("No expression data was provided. Nothing to plot.")
    if(any(!rownames(df) %in% rownames(expression.dat))){
        stop("Can't find expression data for some features\n")
    for(probe in rownames(df)){
        if("Symbol" %in% colnames(df)){
            ylab <- paste(df[probe, 'Symbol'], ylab.type)
        } else {
            ylab <- paste(probe, ylab.type)
        bigplot <- stripplot(expression.dat[ probe, ] ~ factor,
            panel = panel.boxandstrip, groups = factor, ylab = ylab,
            scales = scales, par.settings = par.settings, xlab = xlab)
        miniplot <- .remove.axis.and.padding(bigplot)
        minipng.filename <- paste("mini", probe ,"png", sep='.')
        minipng.file <- file.path(, minipng.filename)
        png(minipng.file, height=40, width=200)
        pushViewport(viewport(angle = 270, height = unit(220, 'points'), 
            width = unit(44, 'points'), name = "VP"))
        print(miniplot, newpage = FALSE)
        pdf.filename <- paste("boxplot", probe, "pdf", sep=".")
        pdf.file <- file.path(, pdf.filename)
        pdf(pdf.file, height=4.5, width=4.5)


I tried changing 

 bigplot <- stripplot(expression.dat[ probe, ] ~ factor,
            panel = panel.boxandstrip, groups = factor, ylab = ylab,
            scales = scales, par.settings = par.settings, xlab = xlab)


      geneCounts <- plotCounts(expression.dat, gene=probe, intgroup=c("condition"), returnData=TRUE)
      geneCounts$count<- geneCounts$count+1
      bigplot<-ggplot(geneCounts, aes(x=condition, y=count, fill=condition)) +
        scale_y_log10() + 
        stat_boxplot(geom ='errorbar') + 
        geom_point( position=position_dodge(width=0.75), size=3)+
        ggtitle(gene) +
        theme(text = element_text(color="#000000", size=25))+
        theme(plot.title = element_text(color="#000000", face="bold", size=32, hjust = 0.5)) +
        theme(axis.title = element_text(color="#000000", face="bold", size=28)) 

This did not work for me.  If anyone has any ideas on how to do this or if it's impossible/not worth it please let me know...


Thank you for your time!

deseq2 ggplot2 reportingtools • 1.2k views

Login before adding your answer.

Traffic: 283 users visited in the last hour
Help About
Access RSS

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

Powered by the version 2.3.6