Question

ReportingTools publish(): no image when .modifyDF argument is included

0

Entering edit mode

psutton • 0

@psutton-20110

Last seen 5.7 years ago

I am learning about using ReportingTools for DESeqDataSet. Using the example in the user guide:

library(DESeq2)
library(ReportingTools)
data(mockRnaSeqData)

conditions <- c(rep("case",3), rep("control", 3))

mockRna.dse <- DESeqDataSetFromMatrix(countData = mockRnaSeqData,
                                        colData = as.data.frame(conditions), 
                                        design = ~ conditions)

colData(mockRna.dse)$conditions <- factor(colData(mockRna.dse)$conditions, 
                                        levels=c("control", "case"))

# Get a DESeqDataSet object
mockRna.dse <- DESeq(mockRna.dse)

# Now the results can be written to a report using the DESeqDataSet object.
des2Report <- HTMLReport(shortName = 'RNAseq_analysis_with_DESeq2',
                           title = 'RNA-seq analysis of differential expression using DESeq2',
                           reportDirectory = "./reports")


# COMMENT 1: I added n=100 so it doesn't take so long
# COMMENT 2: .modifyDF = list(myFunc1, myFunc2) or .modifyDF = list() 
#                         makes the Images disappear from the HTML
publish(mockRna.dse,des2Report, pvalueCutoff=0.05, n = 100, 
          annotation.db="org.Mm.eg.db", factor = colData(mockRna.dse)$conditions,
          reportDir="./reports")  

url <- finish(des2Report)
browseURL(url)

The HTML output has an Image column with stripplots of expression counts for each gene.

It is possible to write functions which modify the data frame, and publish can call these functions if they are passed as a list via the .modifyDF argument, e.g. .modifyDF = list(myFunc1, myFunc2)

However, when I do this, and in particular, even if I pass in an empty list (.modifyDF = list()), the Image column in the resulting HTML is now empty. How can I pass in .modifyDF and also preserve the generation of the stripplot images?

ReportingTools • 1.0k views

ADD COMMENT • link updated 5.7 years ago by James W. MacDonald 67k • written 5.7 years ago by psutton • 0

score 0 · Answer 1 · 2019-03-15

0

Entering edit mode

James W. MacDonald 67k

@james-w-macdonald-5106

Last seen 12 hours ago

United States

In the vignette for RNA-Seq, in section 3, there is this code:

> desReport <- HTMLReport(shortName ='RNAseq_analysis_with_DESeq',
 title ='RNA-seq analysis of differential expression using DESeq', reportDirectory = "./reports")
> publish(res,desReport,name="df",countTable=mockRnaSeqData, pvalueCutoff=0.05,    conditions=conditions,annotation.db="org.Mm.eg.db",   
expName="deseq",reportDir="./reports", .modifyDF=makeDESeqDF)
> finish(desReport)

From which (and the preceding paragraph) you can surmise that the 'default' .modifyDf argument is makeDESeqDF.

If you pass in an argument for .modifyDF, it overwrites the existing function, so if you want to add functions rather than overwriting, you need to pass a list that includes the original function that you like, something like

publish(mockRna.dse,des2Report, pvalueCutoff=0.05, n = 100, 
          annotation.db="org.Mm.eg.db", factor = colData(mockRna.dse)$conditions,
          reportDir="./reports", .modifyDF = list(makeDESeqDF, myFunc1, myFunc2))

ADD COMMENT • link 5.7 years ago James W. MacDonald 67k

0

Entering edit mode

Thank for you this helpful information! The default argument of .modifyDF makes sense.

However, in this particular case, it appears that makeDESeqDF works on data.frames created by DESeq (not DESeq2) after running nbinomTest.

1) The DESeq example from the vignette works perfectly:

library(ReportingTools)
data(mockRnaSeqData)
library(DESeq)                        ### NOT DESeq2

conditions <- c(rep("case",3), rep("control", 3))
cds<-newCountDataSet(mockRnaSeqData, conditions)
cds<-estimateSizeFactors(cds)
cds<-estimateDispersions(cds)
res<-nbinomTest(cds,"control", "case" )

desReport <- HTMLReport(shortName = 'RNAseq_analysis_with_DESeq',
    title = 'RNA-seq analysis of differential expression using DESeq',
    reportDirectory = "./reports")

publish(res,desReport,name="df",countTable=mockRnaSeqData, pvalueCutoff=0.05,
    conditions=conditions,annotation.db="org.Mm.eg.db",
    expName="deseq",reportDir="./reports", .modifyDF=makeDESeqDF)

url <- finish(desReport)

2) In contrast, trying to pass a DESeqDataSet object (DESeq2) to publish works, but if I put in .modifyDF=makeDESeqDF, it fails:

library(ReportingTools)
data(mockRnaSeqData)
library(DESeq2)         ### not DESeq

conditions <- c(rep("case",3), rep("control", 3))

mockRna.dse <- DESeqDataSetFromMatrix(countData = mockRnaSeqData,
                                      colData = as.data.frame(conditions), 
                                      design = ~ conditions)

colData(mockRna.dse)$conditions <- factor(colData(mockRna.dse)$conditions, 
                                          levels=c("control", "case"))

mockRna.dse <- DESeq(mockRna.dse)

des2Report <- HTMLReport(shortName = 'RNAseq_analysis_with_DESeq2',
    title = 'RNA-seq analysis of differential expression using DESeq2',
    reportDirectory = "./reports")

# This works (no .modifyDF argument)
publish(mockRna.dse, des2Report, pvalueCutoff=0.05,
    annotation.db="org.Mm.eg.db", factor = colData(mockRna.dse)$conditions,
    reportDir="./reports")

# This does not work
publish(mockRna.dse, des2Report, pvalueCutoff=0.05,
    annotation.db="org.Mm.eg.db", factor = colData(mockRna.dse)$conditions,
    reportDir="./reports", .modifyDF=makeDESeqDF)

Error in f(df, report, object = object, ...) :
  No genes meet the selection criteria. Try changing the p-value cutoff.

because I think makeDESeqDF is expecting a data.frame from DESeq::nbinomTest and will not work with a DESeqDataSet object.

So, what is the default argument for .modifyDF when a DESeqDataSet is passed to publish? This way I can pass a list with my custom functions and also include the default argument.

ADD REPLY • link 5.7 years ago psutton • 0

1

Entering edit mode

Ah, you're right. For DESeqDataSet objects, there are a couple of functions that coerce the DESeqDataSet to a DESeqResults object, and then get the output DataFrame. So there isn't a 'stock' .modifyDF function for this class.

You probably need to ensure that you set make.plots to TRUE and then set up your functions to modify the resulting data.frame based on what you should expect to get. In other words, here is the function that will be called:

> showMethods(toReportDF, class="DESeqDataSet", includeDefs=TRUE)
Function: toReportDF (package ReportingTools)
object="DESeqDataSet"
function (object, ...) 
{
    .local <- function (object, htmlReport, contrast = NULL, 
        resultName = NULL, annotation.db = NULL, pvalueCutoff = 0.01, 
        lfc = 0, n = 500, sort.by = "pvalue", make.plots = FALSE, 
        ...) 
    {
        if (!"results" %in% mcols(mcols(object))$type) {
            stop("No results found in DESeqDataSet, please run DESeq first.")
        }
        if (is.null(contrast) && !is.null(resultName)) {
            resTab <- results(object, name = resultName)
        }
        else if (!is.null(contrast) && is.null(resultName)) {
            resTab <- results(object, contrast = contrast)
        }
        else {
            resTab <- results(object)
        }
        resTab <- as.data.frame(resTab[which(resTab$padj < pvalueCutoff & 
            abs(resTab$log2FoldChange) > abs(lfc)), ])
        if (!is.nullsort.by) & sort.by %in% colnames(resTab)) 
            resTab <- resTab[order(resTab[, sort.by]), ]
        if (n < nrow(resTab)) 
            resTab <- resTab[1:n, ]
        ann.map.available <- tryCatch(getAnnMap("SYMBOL", annotation.db), 
            error = function(e) {
                return(FALSE)
            })
        if (inherits(ann.map.available, "AnnDbBimap")) {
            check.eg.ids(rownames(resTab), annotation.db)
            fdata <- annotate.genes(rownames(resTab), annotation.db, 
                keytype = "ENTREZID", columns = list(EntrezId = "ENTREZID", 
                  Symbol = "SYMBOL", GeneName = "GENENAME"))
        }
        else {
            fdata <- data.frame(ID = rownames(resTab), stringsAsFactors = FALSE)
        }
        if (make.plots) {
            ret <- data.frame(fdata, Image = rep("", nrow(fdata)), 
                logFC = resTab$log2FoldChange, `p-Value` = resTab$pvalue, 
                `Adjusted p-Value` = resTab$padj, stringsAsFactors = FALSE, 
                check.names = FALSE, row.names = rownames(resTab))
        }
        else {
            ret <- data.frame(fdata, logFC = resTab$log2FoldChange, 
                `p-Value` = resTab$pvalue, `Adjusted p-Value` = resTab$padj, 
                stringsAsFactors = FALSE, check.names = FALSE, 
                row.names = rownames(resTab))
        }
        ret
    }
    .local(object, ...)
}

So your function(s) should assume that they will be getting a data.frame that looks like the output from this function, and then do whatever you want based on that assumption.