limma: adding weights to mroast
1
0
Entering edit mode
@gordon-smyth
Last seen 14 hours ago
WEHI, Melbourne, Australia
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}}
Cancer limma Cancer limma • 1.8k views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 14 hours ago
WEHI, Melbourne, Australia
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}}
ADD COMMENT
0
Entering edit mode
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}}
ADD REPLY

Login before adding your answer.

Traffic: 677 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6