[BioC] limma: adding weights to mroast

[Anthoula Gaigneaux]

Dear Gordon,

Thanks a lot for your answer and for the update.
I don't really need to use a list of vectors as it was a workaround.

Best wishes,
Anthoula.

Anthoula Gaigneaux, PhD
Laboratoire de Biologie Moléculaire et Cellulaire du Cancer (LBMCC)
Hôpital Kirchberg
9, rue Steichen
L-2540 Luxembourg
tel: 00352/ 2468 4046
anthoula.gaigneaux at lbmcc.lu

> 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 inte...{{dropped:6}}