本文关键词：geo芯片数据r语言处理代码

搞微阵芯片数据的同行估计都懂那种痛，原始数据一堆乱麻，limma或者sva跑起来报错，或者就是结果跟预期对不上。这篇不整虚的，直接告诉你怎么用 R 语言搞定 GEO 芯片数据，尤其是那些坑爹的 QC 步骤。

别信什么包安装一下万事大吉。

真的，别。

我之前处理一个 GEO 上的 Affymetrix 数据集，RMA 方法直接报内存溢出。

后来才发现是原始文件解压时路径不对，或者字符编码出问题了。

R 语言对路径和编码很敏感，尤其是你在 Windows 上跑，那些反斜杠 \ 能把你搞死。

一定要用 / 或者 file.path() 函数拼接路径，这是基础中的基础。

先说读取。

大多数 GEO 芯片数据是 CEL 文件。

用 affy 包里的 read.affybatch() 是最常见的。

但这里有个大坑，很多人直接读进去就预处理了。

错大特错。

你得先看原始表达谱图，rma() 之前先做个 simpleaffy::rawPlot()。

看看那些 spot 的正态分布。

如果明显偏态，你的预处理结果绝对歪。

我见过太多人直接出热图，结果全是假信号。

RMA 预处理是标配，但 backgroundCorrect 和 normalize 参数别乱动。

默认值通常是最稳的。

除非你有非常特定的理由，否则别自己去调 method="rma" 里的细节参数，容易翻车。

清洗完背景噪声，接下来是过滤。

这里有个很真实的经验，别把所有探针都留下。

低表达的基因噪声太大，会影响下游的显著性分析。

我建议用 filterfun = function(x) var(x) > 0.5 * mean(x) 或者类似的变异系数过滤掉一半的低信噪比探针。

具体阈值得看你数据的整体质量，一般 GEO 上的公开数据，质量参差不齐，过滤比例大点没事。

然后就是标准化。

RMA 其实自带了归一化，但如果你的数据集是混合批次，或者不同平台混在一起，那就得用 sva 或者 comBat 去批次效应了。

limma 里的 duplicateCorrelation() 处理生物重复时，记得要传对 block 信息，不然协方差矩阵算错，后面 limma 的 t-test 全是废的。

我特别想吐槽一下 annotate 包。

很多新手装 biomaRt 去下载注释，然后卡在镜像服务器上。

国内网络环境，biomaRt 真的慢得令人发指，有时候直接超时。

建议直接用 NCBI 下载的 gene symbol 映射表，手动 merge 进去。

虽然土，但快，还稳。

还有，关于 P 值校正。

很多人跑完 row.ttests() 直接看 adj P < 0.05。

但 GEO 数据里，很多批次效应没去干净时，P 值分布是偏的。

这时候看火山图，会发现两边不对称。

这时候别急着发文，先查查是不是批次没归一好。

我写了一套通用的 geo_chip_pipeline.R，里面包含了从下载、读取、RMA、过滤到 limma 分析的全流程。

核心代码大概长这样：

`r

library(affy)

library(limma)

library(sva)

读取 CEL 文件

eset <- read.affybatch("path/to/cels/", sample.pheno.csv)

预处理

eset_norm <- rma(eset)

表达矩阵

exprs_mat <- exprs(eset_norm)

过滤低变探针

index <- apply(exprs_mat, 1, var) > 1e5

exprs_mat <- exprs_mat[index, ]

设计矩阵和拟合模型

假设分组在 phenotype 的 group 列

design <- model.matrix(~ 0 + group, data = phenodata(eset_norm))

colnames(design) <- make.names(colnames(design), unique=TRUE)

fit <- lmFit(exprs_mat, design)

fit2 <- contrasts.fit(fit, makeContrasts("Treatment - Control", levels=design))

fit2 <- eBayes(fit2)

获取显著差异基因

topTable(fit2, coef = 1, ntop = Inf, adjust.method = "BH")

`

注意，makeContrasts 那里的 levels 一定要对得上 design 矩阵里的列名，格式转换经常出错。

如果报错 Error in solve(...) system is computationally singular，那就是你的样本太少，或者完全共线性了，检查一下 pheno 数据里有没有重复样本或者缺失值。

最后输出结果时，记得把 probe ID 换成 gene symbol，别给人家一堆 1000_f.at 这种 ID，没法看。

用 mapIds() 批量转换，处理一下多对多的情况，取第一个或平均值都行，保持一致就好。

这套流程跑下来，大概几十分钟。

别指望 R 语言有魔法，它只是工具。

数据本身质量不行，代码再精妙也没用。

GEO 上的数据，下载完先花半小时看看 raw plot，这是省心的关键。

别偷懒，不然后面返工更痛苦。

希望这些踩过的坑能帮到你。

代码只是骨架，理解数据才是灵魂。

有问题评论区聊，别藏着。