Skip to contents

Gene set analysis using pre-computed test statistic

Source: R/zenithPR_gsa.R
zenithPR_gsa.Rd

Perform gene set analysis on the result of a pre-computed test statistic. Test whether statistics in a gene set are larger/smaller than statistics not in the set.

Usage

zenithPR_gsa(
  statistics,
  ids,
  geneSets,
  use.ranks = FALSE,
  n_genes_min = 10,
  progressbar = TRUE,
  inter.gene.cor = 0.01,
  coef.name = "zenithPR"
)

Arguments

statistics

pre-computed test statistics

ids

name of gene for each entry in statistics

geneSets

GeneSetCollection

use.ranks

do a rank-based test TRUE or a parametric test FALSE? default: FALSE

n_genes_min

minumum number of genes in a geneset

progressbar

if TRUE, show progress bar

inter.gene.cor

correlation of test statistics with in gene set

coef.name

name of column to store test statistic

Value

  • NGenes: number of genes in this set

  • Correlation: mean correlation between expression of genes in this set

  • delta: difference in mean t-statistic for genes in this set compared to genes not in this set

  • se: standard error of delta

  • p.less: p-value for hypothesis test of H0: delta < 0

  • p.greater: p-value for hypothesis test of H0: delta > 0

  • PValue: p-value for hypothesis test H0: delta != 0

  • Direction: direction of effect based on sign(delta)

  • FDR: false discovery rate based on Benjamini-Hochberg method in p.adjust

  • coef.name: name for pre-computed test statistics. Default: zenithPR

Details

This is the same as zenith_gsa(), but uses pre-computed test statistics. Note that zenithPR_gsa() may give slightly different results for small samples sizes, if zenithPR_gsa() is fed t-statistics instead of z-statistics.

See also

zenith_gsa(), limma::cameraPR()

Examples

# Load packages
library(edgeR)
library(variancePartition)
library(tweeDEseqCountData)

# Load RNA-seq data from LCL's
data(pickrell)
geneCounts = exprs(pickrell.eset)
df_metadata = pData(pickrell.eset)

# Filter genes
# Note this is low coverage data, so just use as code example
dsgn = model.matrix(~ gender, df_metadata)
keep = filterByExpr(geneCounts, dsgn, min.count=5)

# Compute library size normalization
dge = DGEList(counts = geneCounts[keep,])
dge = calcNormFactors(dge)

# Estimate precision weights using voom
vobj = voomWithDreamWeights(dge, ~ gender, df_metadata)

# Apply dream analysis
fit = dream(vobj, ~ gender, df_metadata)
fit = eBayes(fit)

# Load Hallmark genes from MSigDB
# use gene 'SYMBOL', or 'ENSEMBL' id
# use get_GeneOntology() to load Gene Ontology
gs = get_MSigDB("H", to="ENSEMBL")
#> Using cached version from 2024-03-08 18:44:35
   
# Run zenithPR analysis with a test statistic for each gene
tab = topTable(fit, coef='gendermale', number=Inf)
  
res.gsa = zenithPR_gsa(tab$t, rownames(tab), gs)