Question: ReportingTools publish(): no image when .modifyDF argument is included
0
gravatar for psutton
9 weeks ago by
psutton0
psutton0 wrote:

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 • 84 views
ADD COMMENTlink modified 9 weeks ago by James W. MacDonald50k • written 9 weeks ago by psutton0
Answer: ReportingTools publish(): no image when .modifyDF argument is included
0
gravatar for James W. MacDonald
9 weeks ago by
United States
James W. MacDonald50k wrote:

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 COMMENTlink written 9 weeks ago by James W. MacDonald50k

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 REPLYlink written 9 weeks ago by psutton0
1

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.

ADD REPLYlink written 9 weeks ago by James W. MacDonald50k

Thanks! I'll look into that.

ADD REPLYlink written 9 weeks ago by psutton0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 202 users visited in the last hour