Defines functions aggregatePvals

Documented in aggregatePvals

##This function takes as input a matrix of p-values for example obtained 
##from a GSEA on multiple phenotypes, with a row for each gene set and a 
##column for each phenotype and aggregates the p-values by row (i.e. one 
##aggregated p-value for each gene set) according to Fisher or Stouffer's 

aggregatePvals <- function(pvalMatrix, method = "fishers", pAdjustMethod 
	= "BH", order=TRUE) {
	##check 'pvalMatrix'
	if(!is(pvalMatrix, "matrix")) 
		stop("The argument pvalMatrix should be a matrix")
		stop("The argument pvalMatrix should be a matrix with rownames")
	##check 'method'
	if(!(method %in% c("stouffers", "fishers") && length(method)==1 )) 
		stop('The argument method should be one of "stouffers" or 
	##check 'pAdjustMethod'
	paraCheck("pAdjustMethod", pAdjustMethod)
	##do aggregation
	aggr.pval <- rep(1, nrow(pvalMatrix))
	names(aggr.pval) <- rownames(pvalMatrix)
	if (method == "fishers") {
		pvalMatrixLogged <- log(pvalMatrix)
		aggr.pval <- (-2) * apply(pvalMatrixLogged, 1, sum)		
		aggr.pval <- pchisq(aggr.pval, df = 2*ncol(pvalMatrix), 
			lower.tail = FALSE)
	if (method == "stouffers") {
		aggr.pval <- apply(pvalMatrix , 1, 
				function(x) {
					newPval <- pnorm(sum(qnorm(x)) / sqrt(length(x)))
	adj.aggr.pval <- p.adjust(aggr.pval, method = pAdjustMethod)	
	aggr.pval.table <- cbind(aggr.pval, adj.aggr.pval)
		aggr.pval.table <- aggr.pval.table[
			order(aggr.pval.table[, "Aggregated.p.value"]), ]

Try the HTSanalyzeR package in your browser

Any scripts or data that you put into this service are public.

HTSanalyzeR documentation built on Oct. 31, 2019, 7:10 a.m.