GSEAlens is available on Bioconductor. Install the stable release with:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("GSEAlens")The development version can be installed from GitHub:
if (!requireNamespace("pak", quietly = TRUE))
install.packages("pak")
pak::pkg_install("DDL095/GSEAlens")This vignette describes how to prepare the input objects that GSEAlens
consumes. GSEAlens itself does not perform differential expression analysis (DEG);
it accepts already-fitted MArrayLM objects from limma or DESeqDataSet
objects from DESeq2. The following steps are based on the standard
workflows of those packages and are provided here for completeness.
For full DEG workflow documentation, please consult:
We use the airway dataset (4 vs 4 samples, dex treatment).
library(airway)
data(airway)
expression_data <- airwayGSEAlens requires a no-intercept design (~0+group) so that column names
directly correspond to group names, enabling precise contrast construction.
library(edgeR)
library(limma)
group_level <- expression_data$dex
design <- model.matrix(~0+group_level)
colnames(design) <- levels(group_level)
compare_end <- combn(levels(group_level), 2, simplify = FALSE)
contrast_strings <- sapply(compare_end, function(x) paste(x[2], x[1], sep = " - "))
contrast_matrix <- makeContrasts(contrasts = contrast_strings, levels = design)
genes_df <- data.frame(
gene_id = SummarizedExperiment::rowData(expression_data)$gene_id,
symbol = SummarizedExperiment::rowData(expression_data)$symbol,
gene_biotype = SummarizedExperiment::rowData(expression_data)$gene_biotype
)
genes_df$Length <- SummarizedExperiment::rowData(expression_data)$gene_seq_end -
SummarizedExperiment::rowData(expression_data)$gene_seq_start + 1
gsea_limma_voom_data <- edgeR::DGEList(
counts = SummarizedExperiment::assay(expression_data, "counts"),
genes = genes_df,
norm.factors = NULL,
group = group_level,
remove.zeros = TRUE
)
# Filter low-expression genes, keep protein-coding RNAs, deduplicate symbols
total_counts <- edgeR::cpm(gsea_limma_voom_data) |> rowSums()
dup_symbols <- gsea_limma_voom_data$genes$symbol[duplicated(gsea_limma_voom_data$genes$symbol)]
keep <- rep(TRUE, nrow(gsea_limma_voom_data))
for (gene in dup_symbols) {
idx <- which(gsea_limma_voom_data$genes$symbol == gene)
best_idx <- idx[which.max(total_counts[idx])]
remove_idx <- idx[idx != best_idx]
keep[remove_idx] <- FALSE
}
gsea_limma_voom_data <- gsea_limma_voom_data[keep, ]
rownames(gsea_limma_voom_data) <- gsea_limma_voom_data$genes$symbol
keep_biotype <- gsea_limma_voom_data$genes$gene_biotype == "protein_coding"
gsea_limma_voom_data <- gsea_limma_voom_data[keep_biotype, ]
gsea_limma_voom_data <- edgeR::normLibSizes(gsea_limma_voom_data, method = "TMM")
isexpr <- rowSums(edgeR::cpm(gsea_limma_voom_data) > 1) >= 3
gsea_limma_voom_data <- gsea_limma_voom_data[isexpr, ]Perform limma-voom fitting to obtain the fit object.
VoomOutPut <- voom(gsea_limma_voom_data, design)
fit <- lmFit(object = VoomOutPut, design = design) |>
contrasts.fit(contrasts = contrast_matrix) |>
eBayes()library("DESeq2")
dds_se <- DESeqDataSet(expression_data, design = ~ cell + dex)
# Keep protein-coding genes only
gene_biotypes <- SummarizedExperiment::rowData(dds_se)$gene_biotype
keep_protein_coding <- gene_biotypes == "protein_coding"
dds_se <- dds_se[keep_protein_coding, ]
# Remove low-expression genes (require >=10 reads in >=3 samples)
smallestGroupSize <- 3
keep <- rowSums(DESeq2::counts(dds_se) >= 10) >= smallestGroupSize
dds_se <- dds_se[keep, ]
# Deduplicate gene symbols, keeping the highest-count row
rownames(dds_se) <- SummarizedExperiment::rowData(dds_se)$gene_name
total_counts <- rowSums(SummarizedExperiment::assay(dds_se))
dup_genes <- rownames(dds_se)[duplicated(rownames(dds_se))]
keep <- rep(TRUE, nrow(dds_se))
for (gene in dup_genes) {
idx <- which(rownames(dds_se) == gene)
best_idx <- idx[which.max(total_counts[idx])]
remove_idx <- idx[idx != best_idx]
keep[remove_idx] <- FALSE
}
dds_se <- dds_se[keep, ]
# Run DESeq2 fitting
dds_se <- DESeq(dds_se)DDS_rawdata <- expression_data
# Keep protein-coding genes only
gene_biotypes <- SummarizedExperiment::rowData(DDS_rawdata)$gene_biotype
keep_protein_coding <- gene_biotypes == "protein_coding"
DDS_rawdata <- DDS_rawdata[keep_protein_coding, ]
# Remove low-expression genes
smallestGroupSize <- 3
keep_epd <- rowSums(SummarizedExperiment::assay(DDS_rawdata, "counts") >= 10) >= smallestGroupSize
DDS_rawdata <- DDS_rawdata[keep_epd, ]
# Deduplicate gene symbols
rownames(DDS_rawdata) <- SummarizedExperiment::rowData(DDS_rawdata)$gene_name
total_counts <- rowSums(SummarizedExperiment::assay(DDS_rawdata))
keep_name <- rep(TRUE, nrow(DDS_rawdata))
dup_genes <- rownames(DDS_rawdata)[duplicated(rownames(DDS_rawdata))]
for (gene in dup_genes) {
idx <- which(rownames(DDS_rawdata) == gene)
best_idx <- idx[which.max(total_counts[idx])]
remove_idx <- idx[idx != best_idx]
keep_name[remove_idx] <- FALSE
}
DDS_rawdata <- DDS_rawdata[keep_name, ]
# Extract counts and sample metadata into plain matrices
cts <- SummarizedExperiment::assay(DDS_rawdata, "counts")
coldata <- as.data.frame(SummarizedExperiment::colData(DDS_rawdata))
coldata <- coldata[, c("cell", "dex")]
coldata$cell <- factor(coldata$cell)
coldata$dex <- factor(coldata$dex)
# Build DESeqDataSet from matrices and run DESeq2
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ dex)
dds <- DESeq(dds)Once fit, dds_se, and dds are prepared, return to the main vignette:
vignette("GSEAlens")
sessionInfo()## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 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
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] DESeq2_1.53.2 edgeR_4.11.6
## [3] limma_3.69.2 airway_1.33.2
## [5] SummarizedExperiment_1.43.0 Biobase_2.73.2
## [7] GenomicRanges_1.65.1 Seqinfo_1.3.0
## [9] IRanges_2.47.2 S4Vectors_0.51.6
## [11] BiocGenerics_0.59.11 generics_0.1.4
## [13] MatrixGenerics_1.25.0 matrixStats_1.5.0
## [15] BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 SparseArray_1.13.2 lattice_0.22-9
## [4] magrittr_2.0.5 digest_0.6.39 RColorBrewer_1.1-3
## [7] evaluate_1.0.5 grid_4.6.1 bookdown_0.47
## [10] fastmap_1.2.0 jsonlite_2.0.0 Matrix_1.7-6
## [13] BiocManager_1.30.27 scales_1.4.0 codetools_0.2-20
## [16] jquerylib_0.1.4 abind_1.4-8 cli_3.6.6
## [19] rlang_1.3.0 XVector_0.53.0 cachem_1.1.0
## [22] DelayedArray_0.39.4 yaml_2.3.12 otel_0.2.0
## [25] S4Arrays_1.13.0 tools_4.6.1 parallel_4.6.1
## [28] BiocParallel_1.47.0 dplyr_1.2.1 ggplot2_4.0.3
## [31] locfit_1.5-9.12 vctrs_0.7.3 R6_2.6.1
## [34] lifecycle_1.0.5 pkgconfig_2.0.3 pillar_1.11.1
## [37] bslib_0.12.0 gtable_0.3.6 glue_1.8.1
## [40] Rcpp_1.1.2 statmod_1.5.2 tidyselect_1.2.1
## [43] tibble_3.3.1 xfun_0.60 dichromat_2.0-1
## [46] knitr_1.51 farver_2.1.2 htmltools_0.5.9
## [49] rmarkdown_2.31 compiler_4.6.1 S7_0.2.2