# limma ----
library(limma)
## Step 1: log2
if (max(data) > 100) {
expr <- log2(data + 1)
}
## Step 2: between-array normalization
expr.norm <- normalizeBetweenArrays(expr, method = "quantile")
## Step 3: limma DE
design <- model.matrix(~ colData$condition)
fit <- lmFit(expr.norm, design)
fit <- eBayes(fit)
res <- topTable(fit, number = Inf)
# GSVA ----
log_norm <- as.matrix(log2(expr.norm + 1))
gene_set <- list(GENE_123 = c('GENE_1', 'GENE_2', 'GENE_3'),
GENE_456 = c('GENE_4', 'GENE_5', 'GENE_6'))
param <- gsvaParam(
exprData = log_norm,
geneSets = gene_set,
kcdf = "Gaussian"
)
gsva_mat <- gsva(param)
res <- data.frame(score = gsva_mat[GENE_123,], group = colData$condition)
res
with(res, wilcox.test(score ~ group))
# plot ----
p <- res |> ggplot(aes(group, score)) + geom_boxplot() + theme_bw()
p基因表达芯片
基因表达芯片（Gene Expression Microarray）目前较少使用，但是公共数据库中仍有很多值得探索的数据集。
GSE67501 的数据文件标题是 non-normalized data，但其数据不是整数。GEO 网站注明 Experiment type 为 Expression profiling by array，平台为 Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip 也表明这是基因芯片数据，不是测序数据。使用 limma 包计算：先取 log 再 normalize，如果要计算 logTPM 则需再 log 一次。和测序得到的 raw counts 先 normalize 再 log 不一样。
此处展示芯片数据的整理并衔接 GSVA。