1 安装

GSEAlens 可从 Bioconductor 安装。使用以下命令安装稳定版本：

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GSEAlens")

开发版本可从 GitHub 安装：

if (!requireNamespace("pak", quietly = TRUE))
    install.packages("pak")
pak::pkg_install("DDL095/GSEAlens")

2 简介

本 vignette 描述如何准备 GSEAlens 所需的输入对象。GSEAlens 本身不执行差异表达分析（DEG）；它接收来自 limma 的已拟合 MArrayLM 对象或来自 DESeq2 的 DESeqDataSet 对象。以下步骤基于这些包的标准工作流程，此处提供是为了完整性。完整的 DEG 工作流程文档请参阅：

limma 和 edgeR 的 vignette（limma-voom 流程）
DESeq2 的 vignette
airway 包提供示例数据集

3 示例数据

我们使用 airway 数据集（4 vs 4 样本，dex 处理）。

library(airway)
data(airway)
expression_data <- airway

4 limma-voom 流程

GSEAlens 要求无截距设计（~0+group），使列名直接对应组别名称，从而精确构建对比。

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
)
# 过滤低表达基因，保留蛋白质编码 RNA，去重复 symbol
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, ]

进行 limma-voom 拟合，获取 fit 对象。

VoomOutPut <- voom(gsea_limma_voom_data, design)
fit <- lmFit(object = VoomOutPut, design = design) |>
  contrasts.fit(contrasts = contrast_matrix) |>
  eBayes()

5 DESeq2 流程（SummarizedExperiment 输入）

library("DESeq2")
dds_se <- DESeqDataSet(expression_data, design = ~ cell + dex)
# 仅保留蛋白编码基因
gene_biotypes <- SummarizedExperiment::rowData(dds_se)$gene_biotype
keep_protein_coding <- gene_biotypes == "protein_coding"
dds_se <- dds_se[keep_protein_coding, ]
# 去除低表达基因（要求 >=10 reads 且至少出现在 >=3 个样本中）
smallestGroupSize <- 3
keep <- rowSums(DESeq2::counts(dds_se) >= 10) >= smallestGroupSize
dds_se <- dds_se[keep, ]
# 基因符号去重，保留表达量最高的那一行
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, ]
# 运行 DESeq2 拟合
dds_se <- DESeq(dds_se)

6 DESeq2 流程（Count Matrix 输入）

DDS_rawdata <- expression_data
# 仅保留蛋白编码基因
gene_biotypes <- SummarizedExperiment::rowData(DDS_rawdata)$gene_biotype
keep_protein_coding <- gene_biotypes == "protein_coding"
DDS_rawdata <- DDS_rawdata[keep_protein_coding, ]
# 去除低表达基因
smallestGroupSize <- 3
keep_epd <- rowSums(SummarizedExperiment::assay(DDS_rawdata, "counts") >= 10) >= smallestGroupSize
DDS_rawdata <- DDS_rawdata[keep_epd, ]
# 基因符号去重
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, ]
# 抽取 counts 矩阵和样本元数据
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)
# 从矩阵构建 DESeqDataSet 并运行 DESeq2
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ dex)
dds <- DESeq(dds)

7 下一步

fit、dds_se 和 dds 准备就绪后，返回主 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

GSEAlens：从 limma-voom 和 DESeq2 流程准备输入数据

2026-08-09