微阵列数据分析

2025年12月19日 · Tutorials

摘要： 使用 R 语言中的 Bioconductor 分析 Affymetrix 微阵列数据。使用 affy/oligo 加载 CEL 文件，使用 RMA 进行归一化，使用 limma 进行差异表达分析。尽管 RNA-seq 已成为主流，但对于现有数据集和 GEO 数据的再分析，该方法仍然相关。

为什么在 2026 年还要使用微阵列？

虽然 RNA-seq 在很大程度上已取代微阵列用于新研究，但在以下情况下你仍会遇到微阵列数据：

✅ 重新分析已发表的数据集（GEO 拥有数百万个样本）
✅ 跨研究的荟萃分析（标准化平台）
✅ 临床诊断（FDA 批准的平台）
✅ 预算限制（对于大型队列研究仍然更便宜）

微阵列 vs RNA-seq

特性	微阵列	RNA-seq
成本	$50-100/样本	$200-500/样本
新转录本	❌ 否	✅ 是
动态范围	有限（~3 个数量级）	宽（>5 个数量级）
定量	相对	绝对
归一化	需要跨样本	可单样本进行
数据大小	~10 MB/样本	~2-10 GB/样本
最佳适用	已知基因，大型队列	发现，异构体，新基因

工作流程概述

CEL 文件 → RMA 归一化 → 差异表达 → 功能分析

微阵列分析流程

第一部分：数据结构

ExpressionSet 对象

Bioconductor 中用于表达数据的标准容器：

library(Biobase)

# 结构
ExpressionSet (storageMode: lockedEnvironment)
  assayData: 12625 个特征，22 个样本      # 表达矩阵
  phenoData:                                 # 样本元数据
    sampleNames: CLL11.CEL, CLL12.CEL, ...
    varLabels: SampleID, Disease, Age
  featureData:                               # 探针注释
  annotation: hgu95av2                       # 芯片类型

提取组件

# 表达矩阵
expr_matrix <- exprs(eset)
# 行 = 探针/基因，列 = 样本

# 样本元数据
metadata <- pData(eset)
# 列 = 临床变量（疾病、年龄等）

# 探针注释
features <- fData(eset)
# 将探针映射到基因符号

第二部分：加载数据

从 GEO 下载（推荐）

library(GEOquery)

# 下载处理后的数据
gse <- getGEO("GSE1009", GSEMatrix = TRUE)
eset <- gse[[1]]  # ExpressionSet 对象

# 检查平台
annotation(eset)

# 查看样本信息
pData(eset)[1:5, 1:5]

从原始 CEL 文件加载（Affymetrix）

选项 1：affy 包（较旧芯片）

library(affy)

# 从目录读取 CEL 文件
raw_data <- ReadAffy(celfile.path = "CEL_files/")

# 检查数据
sampleNames(raw_data)
dim(exprs(raw_data))  # 探针水平

# RMA 归一化（背景校正 + 分位数归一化 + 汇总）
eset <- rma(raw_data)

选项 2：oligo 包（较新芯片）

library(oligo)

# 列出 CEL 文件
cel_files <- list.celfiles("CEL_files/", full.names = TRUE)

# 读取
raw_data <- read.celfiles(cel_files)

# RMA 归一化
eset <- rma(raw_data)

从矩阵创建 ExpressionSet

library(Biobase)

# 如果你有独立的表达矩阵和元数据
expr_matrix <- read.csv("expression.csv", row.names = 1)
metadata <- read.csv("metadata.csv", row.names = 1)

# 创建带注释的元数据
pheno <- new("AnnotatedDataFrame", data = metadata)

# 创建 ExpressionSet
eset <- ExpressionSet(
  assayData = as.matrix(expr_matrix),
  phenoData = pheno,
  annotation = "hgu95av2"  # 芯片类型
)

第三部分：归一化

RMA（稳健多阵列平均）

Affymetrix 数据的行业标准。 结合了：

背景校正（使用 PM/MM 探针）
分位数归一化（跨阵列）
汇总（探针集 → 基因）

# RMA 已包含分位数归一化
eset <- rma(raw_data)

# 验证归一化
boxplot(exprs(eset), main = "RMA 后", las = 2)

重要： 在 RMA 之后不要再运行 normalizeBetweenArrays()（已经完成）。

替代方法

# MAS5（旧方法）
eset_mas5 <- mas5(raw_data)

# GCRMA（考虑 GC 含量）
library(gcrma)
eset_gcrma <- gcrma(raw_data)

# VSN（方差稳定化）
library(vsn)
eset_vsn <- vsnrma(raw_data)

建议： 除非有特定原因，否则坚持使用 RMA。

第四部分：质量控制

归一化前

library(affyPLM)

# 拟合探针水平模型
plm_fit <- fitPLM(raw_data)

# RNA 降解图
RNAdeg <- AffyRNAdeg(raw_data)
plotAffyRNAdeg(RNAdeg)

# 伪图像（空间伪影）
image(raw_data[, 1])  # 第一个阵列

归一化后

# 箱线图（应已对齐）
boxplot(exprs(eset), main = "表达分布")

# PCA
library(ggplot2)
pca <- prcomp(t(exprs(eset)), scale. = TRUE)
pca_df <- data.frame(PC1 = pca$x[,1], PC2 = pca$x[,2], 
                     Group = pData(eset)$Disease)
ggplot(pca_df, aes(PC1, PC2, color = Group)) + 
  geom_point(size = 3)

# 样本距离热图
library(pheatmap)
dists <- dist(t(exprs(eset)))
pheatmap(as.matrix(dists))

第五部分：探针到基因的映射

注释探针

library(annotate)
library(hgu95av2.db)  # 更改为你的芯片

# 获取基因符号
symbols <- getSYMBOL(featureNames(eset), "hgu95av2.db")
fData(eset)$GeneSymbol <- symbols

# 获取 Entrez ID
entrez <- getEG(featureNames(eset), "hgu95av2.db")
fData(eset)$EntrezID <- entrez

# 查看映射
head(fData(eset))

处理每个基因的多个探针

library(limma)

# 当多个探针映射到同一基因时，取平均表达
averaged <- avereps(exprs(eset), ID = fData(eset)$GeneSymbol)

第六部分：使用 limma 进行差异表达分析

基本的两组比较

library(limma)

# 设计矩阵
design <- model.matrix(~ 0 + factor(eset$Disease))
colnames(design) <- c("Control", "Disease")

# 拟合线性模型
fit <- lmFit(eset, design)

# 定义对比
contrast_matrix <- makeContrasts(
  Disease_vs_Control = Disease - Control,
  levels = design
)

# 拟合对比
fit2 <- contrasts.fit(fit, contrast_matrix)
fit2 <- eBayes(fit2)

# 提取结果
results <- topTable(fit2, coef = "Disease_vs_Control", number = Inf, adjust.method = "BH")
head(results)

多组分析

# 多组设计
design <- model.matrix(~ 0 + eset$CellType)
colnames(design) <- levels(factor(eset$CellType))

fit <- lmFit(eset, design)

# 多重对比
contrast_matrix <- makeContrasts(
  TypeA_vs_TypeB = TypeA - TypeB,
  TypeA_vs_TypeC = TypeA - TypeC,
  TypeB_vs_TypeC = TypeB - TypeC,
  levels = design
)

fit2 <- contrasts.fit(fit, contrast_matrix)
fit2 <- eBayes(fit2)

# 查看所有比较
summary(decideTests(fit2))

包含协变量

# 调整批次效应或其他协变量
design <- model.matrix(~ Disease + Batch + Age, data = pData(eset))

fit <- lmFit(eset, design)
fit <- eBayes(fit)

# 疾病效应的结果（已针对批次和年龄调整）
results <- topTable(fit, coef = "DiseaseCase", number = Inf)

第七部分：可视化

火山图

library(ggplot2)

# 准备数据
results$significant <- results$adj.P.Val < 0.05 & abs(results$logFC) > 1

ggplot(results, aes(x = logFC, y = -log10(P.Value), color = significant)) +
  geom_point(alpha = 0.5) +
  scale_color_manual(values = c("gray", "red")) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  labs(title = "火山图", x = "Log2 倍数变化", y = "-Log10 P 值")

顶部基因热图

# 前 50 个差异表达基因
top_genes <- rownames(results)[1:50]
heatmap_data <- exprs(eset)[top_genes, ]

pheatmap(heatmap_data,
         scale = "row",
         annotation_col = pData(eset)[, c("Disease"), drop = FALSE],
         show_rownames = TRUE,
         fontsize_row = 6)

MA 图

limma::plotMA(fit2, coef = 1)

第八部分：功能富集

GO 富集

library(clusterProfiler)
library(org.Hs.eg.db)

# 获取显著基因
sig_genes <- results$EntrezID[results$adj.P.Val < 0.05]
sig_genes <- na.omit(sig_genes)

# GO 富集
go_results <- enrichGO(
  gene = sig_genes,
  OrgDb = org.Hs.eg.db,
  ont = "BP",  # 生物过程
  pAdjustMethod = "BH",
  qvalueCutoff = 0.05
)

# 可视化
dotplot(go_results, showCategory = 20)

KEGG 通路

kegg_results <- enrichKEGG(
  gene = sig_genes,
  organism = "hsa",  # 人类
  pAdjustMethod = "BH",
  qvalueCutoff = 0.05
)

dotplot(kegg_results)

完整分析脚本

# 加载库
library(GEOquery)
library(limma)
library(pheatmap)

# 1. 下载数据
gse <- getGEO("GSE1009", GSEMatrix = TRUE)
eset <- gse[[1]]

# 2. 过滤低表达探针
eset <- eset[rowMeans(exprs(eset)) > 5, ]

# 3. 检查归一化
boxplot(exprs(eset), las = 2)

# 4. 设计矩阵
design <- model.matrix(~ 0 + Disease, data = pData(eset))
colnames(design) <- c("Control", "Disease")

# 5. 差异表达
fit <- lmFit(eset, design)
contrast <- makeContrasts(Disease - Control, levels = design)
fit2 <- contrasts.fit(fit, contrast)
fit2 <- eBayes(fit2)

# 6. 结果
results <- topTable(fit2, number = Inf, adjust.method = "BH")

# 7. 导出
write.csv(results, "DEG_results.csv")

# 8. 前 50 个基因热图
top_genes <- rownames(results)[1:50]
pheatmap(exprs(eset)[top_genes, ], scale = "row", 
         annotation_col = pData(eset)[, "Disease", drop = FALSE])

故障排除

"未找到注释包"

修复： 安装芯片特定注释：

# 列出可用包
BiocManager::available("hgu")

# 安装你的芯片包
BiocManager::install("hgu95av2.db")

RMA 因内存错误失败

修复： 分批处理：

# 减少内存占用
rm(raw_data)
gc()

PCA 中可见批次效应

修复： 使用 ComBat：

library(sva)
combat_expr <- ComBat(exprs(eset), batch = pData(eset)$Batch)
exprs(eset) <- combat_expr

工具与软件

R 包

核心： affy, oligo, limma, Biobase
注释： annotate, org.Hs.eg.db, hgu95av2.db
可视化： pheatmap, ggplot2
富集： clusterProfiler, DOSE

GUI 替代方案

Chipster： chipster.csc.fi
用于微阵列分析的点选式界面（无需编码）。

资源

Bioconductor 工作流程： bioconductor.org/packages/workflows
limma 用户指南： bioconductor.org/packages/limma
GEO 数据库： ncbi.nlm.nih.gov/geo

另请参阅：

Microarray