[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