[BioC] limma: adding weights to mroast

Gordon K Smyth smyth at wehi.EDU.AU
Wed Jan 11 00:19:03 CET 2012


Dear Anthoula,

Sorry, my advice applied to roast() rather than to mroast().  You are 
quite correct that there is a problem with passing gene.weights to 
mroast(), and your suggested fix should work ok.  Thanks for bringing it 
to my attention.

I have today committed updates to roast() and mroast() so that mroast() 
will now handle a gene.weights vector such as yours automatically. 
gene.weights can simply be a vector of the same length as nrow(y) like 
gene.weights=fit$Amean.  I've actually made the changes to roast() rather 
than to mroast().

mroast() still doesn't allow the user to provide a list of gene.weights 
vectors, but that isn't necessary for the example you give.  Let me know 
if you do need this full generality.

Best wishes
Gordon

---------------------------------------------
Professor Gordon K Smyth,
Bioinformatics Division,
Walter and Eliza Hall Institute of Medical Research,
1G Royal Parade, Parkville, Vic 3052, Australia.
http://www.wehi.edu.au
http://www.statsci.org/smyth

On Mon, 2 Jan 2012, Gordon K Smyth wrote:

> Dear Anthoula,
>
> You are getting an error message from mroast() because the vector of gene 
> weights must be of the same length as the number of the genes in the set, not 
> of length equal to all the genes in your dataset.  This is what the error 
> message is telling you.  Hence you need
>
>  gene.weights=fitA.genes$Amean[inds]
>
> rather than
>
>  gene.weights=fitA.genes$Amean
>
> There is no need to change the code of the function itself.
>
> Best wishes
> Gordon
>
>> Date: Thu, 22 Dec 2011 14:18:11 +0100
>> From: Anthoula Gaigneaux <anthoula.gaigneaux at lbmcc.lu>
>> To: bioconductor at r-project.org
>> Subject: [BioC] limma: adding weights to mroast
>> 
>> Dear List,
>> 
>> First thank you for this very nice package.
>> 
>> I'm currently using mroast in limma for differential geneset testing on
>> the basis of a MaList. Adding expression-related weights like explained
>> in the paper (Wu 2010) generated errors.
>> 
>> When adding weights as a vector, dimensions did not fit (see below).
>> Then I tried to build a weight list similar to iset object, but the same
>> error prompted. I've searched in several documentations, including
>> mailing list and changelog, but have seen nothing related.
>> 
>> Looking at the code, I've fixed the issue by adding a subscript in the
>> mroast function (gene.weights[[i]]). Is this fix is ok, or if there's
>> another way to define weights for mroast.
>> 
>> Here is the code:
>> ##### test
>> identical(ma.genes$genes$SYMBOL, fitA.genes$genes$SYMBOL )
>> gmt <- getGmt("c2.cp.v3.0.symbols.gmt" )
>> universe = ma.genes$genes$SYMBOL
>> gmt.list <- geneIds(gmt)
>> inds <- symbols2indices(gmt.list, universe)
>> design <- fitA.genes$design
>> contrast <- fitA.genes$contrasts
>> test <- mroast(iset=inds, y= ma.genes, design, contrast=
>> contrast[,4],gene.weights=fitA.genes$Amean, nrot=10000)
>>     Error in roast(iset = iset[[i]], y = y, design = design, contrast =
>>     contrast,  :
>>       length of gene.weights disagrees with size of set
>> 
>> weights <- lapply(inds, function(x) w <- fitA.genes$Amean[x] )
>> test <- mroast(iset=inds, y= ma.genes, design, contrast= contrast[,4],
>> gene.weights= weights, nrot=10000)
>>     Error in roast(iset = iset[[i]], y = y, design = design, contrast =
>>     contrast,  :
>>       length of gene.weights disagrees with size of set
>> 
>> ##### adapt function
>> mroast2 <- function (iset = NULL, y, design, contrast = ncol(design),
>> set.statistic = "mean",
>>     gene.weights = NULL, array.weights = NULL, block = NULL,
>>     correlation, var.prior = NULL, df.prior = NULL, trend.var = FALSE,
>>     nrot = 999, adjust.method = "BH")
>> {
>>     if (is.null(iset))
>>         iset <- rep(TRUE, nrow(y))
>>     if (!is.list(iset))
>>         iset <- list(set = iset)
>>     nsets <- length(iset)
>>     if (is.null(names(iset)))
>>         names(iset) <- paste("set", 1:nsets, sep = "")
>>     fit <- lmFit(y, design = design, weights = array.weights,
>>         block = block, correlation = correlation)
>>     covariate <- NULL
>>     if (trend.var) {
>>         covariate <- fit$Amean
>>         if (is.null(covariate))
>>             covariate <- rowMeans(as.matrix(y))
>>     }
>>     sv <- squeezeVar(fit$sigma^2, df = fit$df.residual, covariate =
>> covariate)
>>     var.prior <- sv$var.prior
>>     df.prior <- sv$df.prior
>>     pv <- adjpv <- active <- array(0, c(nsets, 3), dimnames =
>> list(names(iset),
>>         c("Mixed", "Up", "Down")))
>>     if (nsets < 1)
>>         return(pv)
>>     for (i in 1:nsets) {
>>         out <- roast(iset = iset[[i]], y = y, design = design,
>>             contrast = contrast, set.statistic = set.statistic,
>>             gene.weights = gene.weights[[i]], array.weights =
>> array.weights,
>>             block = block, correlation = correlation, var.prior =
>> var.prior,
>>             df.prior = df.prior, nrot = nrot)[[1]]
>>         pv[i, ] <- out$P.Value
>>         active[i, ] <- out$Active.Prop
>>     }
>>     adjpv[, "Mixed"] <- p.adjust(pv[, "Mixed"], method = adjust.method)
>>     adjpv[, "Up"] <- p.adjust(pv[, "Up"], method = adjust.method)
>>     adjpv[, "Down"] <- p.adjust(pv[, "Down"], method = adjust.method)
>>     list(P.Value = pv, Adj.P.Value = adjpv, Active.Proportion = active)
>> }
>> 
>> test <- mroast2(iset=inds, y= ma.genes, design, contrast= contrast[,4],
>> gene.weights= weights, nrot=10000) #ok
>> 
>> ##### sessionInfo()
>> R version 2.14.0 (2011-10-31)
>> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>> 
>> locale:
>> [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
>> 
>> attached base packages:
>> [1] stats     graphics  grDevices utils     datasets  methods   base
>> 
>> other attached packages:
>> [1] GSEABase_1.16.0      graph_1.32.0         annotate_1.32.0
>> [4] AnnotationDbi_1.16.0 Biobase_2.14.0       limma_3.10.0
>> 
>> loaded via a namespace (and not attached):
>> [1] DBI_0.2-5      IRanges_1.12.1 RSQLite_0.10.0 tools_2.14.0
>> [5] XML_3.4-3      xtable_1.6-0
>> 
>> Thanks for helping and merry Christmas,
>> 
>> Anthoula Gaigneaux, PhD
>> Bioinformatician
>> Laboratoire de Biologie Moléculaire et Cellulaire du Cancer (LBMCC)
>> Hôpital Kirchberg
>> 9, rue Steichen
>> L-2540 Luxembourg
>> anthoula.gaigneaux at lbmcc.lu

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:5}}


More information about the Bioconductor mailing list