[BioC] limma: adding weights to mroast

[Gordon K Smyth]

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}}