zenith_gsa-methods: Perform gene set analysis using zenith

In GabrielHoffman/zenith: Gene set analysis following differential expression using linear (mixed) modeling with dream

Perform gene set analysis using zenith

Description

Perform a competitive gene set analysis accounting for correlation between genes.

Usage

zenith_gsa(
  fit,
  geneSets,
  coefs,
  use.ranks = FALSE,
  n_genes_min = 10,
  inter.gene.cor = 0.01,
  progressbar = TRUE,
  ...
)

## S4 method for signature 'MArrayLM,GeneSetCollection'
zenith_gsa(
  fit,
  geneSets,
  coefs,
  use.ranks = FALSE,
  n_genes_min = 10,
  inter.gene.cor = 0.01,
  progressbar = TRUE,
  ...
)

Arguments

fit	results from dream()
geneSets	GeneSetCollection
coefs	list of coefficients to test using topTable(fit, coef=coefs[[i]])
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
inter.gene.cor	if NA, estimate correlation from data. Otherwise, use specified value
progressbar	if TRUE, show progress bar
...	other arguments

Details

This code adapts the widely used camera() analysis \insertCitewu2012camerazenith in the limma package \insertCiteritchie2015limmazenith to the case of linear (mixed) models used by variancePartition::dream().

Value

data.frame of results for each gene set and cell type

References

\insertAllCited

See Also

limma::camera

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")
   
# Run zenith analysis
res.gsa = zenith_gsa(fit, gs, 'gendermale', progressbar=FALSE )

# Show top gene sets
head(res.gsa, 2)

# for each cell type select 3 genesets with largest t-statistic
# and 1 geneset with the lowest
# Grey boxes indicate the gene set could not be evaluted because
#    to few genes were represented
plotZenithResults(res.gsa)

GabrielHoffman/zenith documentation built on March 10, 2024, 11:43 p.m.