Results from the univariate regressions performed using can be combined in a post-processing step to perform multivariate hypothesis testing. In this example, we fit on transcript-level counts and then perform multivariate hypothesis testing by combining transcripts at the gene-level. This is done with the function.

Import transcript-level counts

Read in transcript counts from the package.

NOTE: this works with version <= 1.40.0. It fails as of Nov 2, 2026.

library(readr)
library(tximport)
library(tximportData)

# specify directory
path <- system.file("extdata", package = "tximportData")

# read sample meta-data
samples <- read.table(file.path(path, "samples.txt"), header = TRUE)
samples.ext <- read.table(file.path(path, "samples_extended.txt"), header = TRUE, sep = "\t")

# read assignment of transcripts to genes
# remove genes on the PAR, since these are present twice
tx2gene <- read_csv(file.path(path, "tx2gene.gencode.v27.csv"))
tx2gene <- tx2gene[grep("PAR_Y", tx2gene$GENEID, invert = TRUE), ]

# read transcript-level quatifictions
files <- file.path(path, "salmon", samples$run, "quant.sf.gz")
txi <- tximport(files, type = "salmon", txOut = TRUE)

# Create metadata simulating two conditions
sampleTable <- data.frame(condition = factor(rep(c("A", "B"), each = 3)))
rownames(sampleTable) <- paste0("Sample", 1:6)

Standard dream analysis

Perform standard analysis at the transcript-level

library(variancePartition)
library(edgeR)

# Prepare transcript-level reads
dge <- DGEList(txi$counts)
design <- model.matrix(~condition, data = sampleTable)
isexpr <- filterByExpr(dge, design)
dge <- dge[isexpr, ]
dge <- calcNormFactors(dge)

# Estimate precision weights
vobj <- voomWithDreamWeights(dge, ~condition, sampleTable)

# Fit regression model one transcript at a time
fit <- dream(vobj, ~condition, sampleTable)
fit <- eBayes(fit)

Multivariate analysis

Combine the transcript-level results at the gene-level. The mapping between transcript and gene is stored in as a list.

# Prepare transcript to gene mapping
# keep only transcripts present in vobj
# then convert to list with key GENEID and values TXNAMEs
keep <- tx2gene$TXNAME %in% rownames(vobj)
tx2gene.lst <- unstack(tx2gene[keep, ])

# Run multivariate test on entries in each feature set
# Default method is "FE.empirical", but use "FE" here to reduce runtime
res <- mvTest(fit, vobj, tx2gene.lst, coef = "conditionB", method = "FE")

# truncate gene names since they have version numbers
# ENST00000498289.5 -> ENST00000498289
res$ID.short <- gsub("\\..+", "", res$ID)

Gene set analysis

Perform gene set analysis using on the gene-level test statistics.

# must have zenith > v1.0.2
library(zenith)
library(GSEABase)

gs <- get_MSigDB("C1", to = "ENSEMBL")

df_gsa <- zenithPR_gsa(res$stat, res$ID.short, gs, inter.gene.cor = .05)

head(df_gsa)

Session info

## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.5 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB             
##  [4] LC_COLLATE=C               LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
## [10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tximportData_1.41.1      tximport_1.41.1          readr_2.2.0             
##  [4] edgeR_4.99.12            pander_0.6.6             variancePartition_1.43.2
##  [7] BiocParallel_1.47.0      limma_3.99.4             ggplot2_4.0.3           
## [10] knitr_1.52              
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1     dplyr_1.2.1          farver_2.1.2         S7_0.2.2            
##  [5] bitops_1.1-0         fastmap_1.2.0        digest_0.6.39        lifecycle_1.0.5     
##  [9] statmod_1.5.2        magrittr_2.0.5       compiler_4.6.1       rlang_1.3.0         
## [13] sass_0.4.10          tools_4.6.1          yaml_2.3.12          labeling_0.4.3      
## [17] bit_4.6.0            plyr_1.8.9           RColorBrewer_1.1-3   KernSmooth_2.23-27  
## [21] withr_3.0.3          purrr_1.2.2          numDeriv_2016.8-1.1  BiocGenerics_0.59.12
## [25] grid_4.6.1           aod_1.3.3            caTools_1.18.4       scales_1.4.0        
## [29] gtools_3.9.5         iterators_1.0.14     MASS_7.3-66          dichromat_2.0-1     
## [33] cli_3.6.6            mvtnorm_1.4-2        crayon_1.5.3         rmarkdown_2.32      
## [37] reformulas_0.4.4     generics_0.1.4       otel_0.2.0           tzdb_0.5.0          
## [41] reshape2_1.4.5       minqa_1.2.8          cachem_1.1.0         stringr_1.6.0       
## [45] splines_4.6.1        parallel_4.6.1       matrixStats_1.5.0    vctrs_0.7.3         
## [49] boot_1.3-32          Matrix_1.7-6         jsonlite_2.0.0       hms_1.1.4           
## [53] bit64_4.8.6          pbkrtest_0.5.5       jquerylib_0.1.4      tidyr_1.3.2         
## [57] snow_0.4-4           glue_1.8.1           nloptr_2.2.1         codetools_0.2-20    
## [61] stringi_1.8.9        gtable_0.3.6         EnvStats_3.1.0       lme4_2.1-0          
## [65] lmerTest_3.2-1       tibble_3.3.1         remaCor_0.0.20       pillar_1.11.1       
## [69] htmltools_0.5.9      gplots_3.3.0         R6_2.6.1             Rdpack_2.6.6        
## [73] vroom_1.7.1          evaluate_1.0.5       lattice_0.23-1       Biobase_2.73.2      
## [77] rbibutils_2.4.1      backports_1.5.1      RhpcBLASctl_0.23-42  broom_1.0.13        
## [81] fANCOVA_0.6-1        corpcor_1.6.10       bslib_0.12.0         Rcpp_1.1.2          
## [85] nlme_3.1-172         xfun_0.62            pkgconfig_2.0.3

<>

References