简介：面向生物信息学入门与进阶研究者，这份资料聚焦GEO数据库芯片数据的差异表达分析，系统讲解R语言limma包从数据下载到结果可视化的完整流程。内容涵盖GEOquery获取数据、affy/oligo预处理、实验设计矩阵构建、lmFit线性建模、eBayes经验贝叶斯修正，以及topTable筛选差异基因和火山图、热图展示等核心环节，并延伸至clusterProfiler等功能富集分析，帮助读者快速掌握高通量数据挖掘的实用技能。包体共12个文件，以R脚本、mp4操作录屏和辅助代码文件为主，另有可执行程序与示例数据压缩包，整体约378.5MB。配套视频按数据下载、数据转换、差异分析、热图绘制逐步演示，R脚本可直接修改运行，适合边看边练。目前已有2362人学习下载，资源中特意收录了R语言软件与万能代码，对零基础用户较为友好，可有效降低环境配置门槛，让使用者将精力集中于分析思路与生物学解读。 做GEO数据库的差异分析，limma包差不多是绕不开的第一选择。很多刚接触生信的同学下载完GEO数据，手里拿到一个表达矩阵，却不知道下一步怎么走；还有些人被各种教程带着跑了一遍，结果换一个数据集就报错。我前阵子正好整理了一套完整的差异分析流程，从数据下载到limma输出差异基因列表，每一步都踩了一遍坑。这篇文章就把这套流程掰开揉碎讲清楚，把那些教程里没写透的细节补上。
这套东西适合谁？主要给两类人：一是刚开始接触芯片数据、转录组数据，想用公共数据库做分析的学生和科研人员；二是已经在用GEO但总在数据预处理和limma参数上传出问题，想搞清楚底层逻辑的同学。读完你至少能独立跑通一个GSE数据集的差异分析，并且明白每一步为什么要这么做。
1. 差异分析的整体思路拆解
1.1 为什么选limma包而不是其他工具
limma（Linear Models for Microarray and RNA-Seq Data）最初是专门为芯片数据设计的线性模型分析工具，后来扩展到RNA-seq数据。它在GEO差异分析里占据统治地位，核心原因是它对小样本量的处理极其稳健。医学和生物学实验的重复样本通常只有3到6个，这种规模下很多统计方法容易出问题，而limma通过经验贝叶斯方法从全局基因中借用信息来稳定方差估计，即使样本量很小也能给出可靠的统计推断。
我在实际使用中对比过edgeR、DESeq2和limma，edgeR和DESeq2虽然也常用，但它们更偏RNA-seq的count数据，输入格式和处理思路都不同。GEO数据库里大量的Affymetrix、Agilent芯片数据，原始格式本身就不是count，limma天然适配这些平台。另外limma的设计矩阵和对比矩阵逻辑非常灵活，多样本分组、配对设计、时间序列都能处理，一套代码框架可以应对绝大多数差异分析场景。
1.2 分析前必须搞清楚的三件事
拿到一个GEO数据集，别急着写代码，先花十分钟搞清楚三件事：表达矩阵、分组信息、平台注释。这三件事是limma分析的输入前提，缺一不可，而且每件都有坑。
表达矩阵就是基因或探针在样本中的表达量，行是基因/探针，列是样本。GEO下载的数据格式五花八门，有原始CEL文件、有处理过的series matrix文件、还有补充文件里的各种表格。做差异分析建议直接用series matrix文件里提取的表达矩阵，省去自己读CEL文件的麻烦。要注意拿到手的值是否已log2转换，这一步后面会重点讲。
分组信息是差异分析的关键。你要清楚比较的双方是谁和谁，比如疾病组vs正常组、处理组vs对照组。分组信息通常可以在GEO页面的sample信息里找到，但也常常需要从论文的补充材料里找，因为GEO有时候存储的分组信息并不完整。
平台注释是探针到基因名的映射关系，这是GEO分析里最大的坑之一。芯片的每个探针对应哪个基因，需要平台的注释文件（GPL文件）来翻译。不同芯片平台的注释文件格式不一样，注释的好坏直接影响后续分析的质量。
2. 数据下载和表达矩阵整理
2.1 GEO数据怎么选、怎么判断质量
在GEO数据库（NCBI GEO）搜索框输入疾病或组织相关关键词，限定数据集类型，会出来一堆结果。选数据集时有几个关键标准：样本量不要太少（每组至少3个以上）、明确的组别设计、平台信息完整、样本注释不混乱。还有一点很多人忽略，就是看看这个数据集是否已经做过分层或基质校正，也就是normalization，因为有些数据是官方已经处理过的，有些则需要自己处理。
还有一个很实用的小技巧，先点开数据集的GEO2R看看它内置的分组信息是否和论文一致。GEO2R本身就是用GEOquery和limma做的在线分析工具，它给出的分组标签能帮你快速判断这个数据集的分组信息是否明确。有时候论文里的分组是A组和B组，GEO页面上却写着“disease: 1/0”，这种就需要自己根据样本属性重新构建分组。
2.2 用GEOquery下载并提取表达矩阵
下载GEO数据，R语言里最方便的是GEOquery包。基础用法是getGEO("GSE编号", destdir = "数据存放目录")，返回的对象里有表达矩阵和样本信息。但这里有个容易困惑的点：返回对象的结构因数据类型而异，有的是ExpressionSet对象，有的只是列表。
从ExpressionSet里提取表达矩阵用exprs()函数，提取样本信息用pData()。这里有个常见问题，exprs()出来的矩阵有时候行名是探针ID，有些平台则是基因名，你需要先看清楚。另外，getGEO默认下载的Series Matrix文件里，表达值大多是已经做过RMA或MAS5标准化的，可以直接用，但有些数据集的表达值范围很奇怪，比如出现负数或很大值，就需要检查是否该取log。
library(GEOquery) # 下载并读取GSE数据 gse <- getGEO("GSE1000", GSEMatrix = TRUE, getGPL = FALSE) expr_matrix <- exprs(gse[[1]]) sample_info <- pData(gse[[1]]) # 查看表达矩阵前几行，判断数据范围和格式 head(expr_matrix[, 1:5])2.3 探针注释与ID转换的完整流程
探针注释这个环节，很多人会在这里卡住或者用错方法。最标准的做法是下载对应GPL平台的注释文件，用getGEO("GPL编号")获取，如果网络不好也可以直接从GEO的FTP站点下载。注释文件里包含探针ID到基因Symbol的对应关系。
拿到注释后要做的三件事：一是把探针ID和表达矩阵的行名对齐，二是合并注释得到基因级别的表达矩阵，三是处理多个探针对应同一基因的情况，通常是取最大值或平均值。我用apply(x, 2, max)取最大值，这样能保留表达较强的探针信号，但也有人偏好取平均，各有千秋，关键是保持一致并在论文里说明。
# 获取平台注释，以GPL570为例 gpl <- getGEO("GPL570", getGPL = TRUE) # 提取探针到基因的映射表格 gpl_table <- Table(gpl) probe2gene <- gpl_table[, c("ID", "Gene Symbol")] # 合并探针表达矩阵和基因注释，然后去掉无基因注释的探针 expr_with_gene <- merge(probe2gene, expr_matrix, by.x = "ID", by.y = "row.names")这里有个比较隐蔽的坑：merge之后，原来是数值矩阵的表达数据可能变成数据框，行名也丢了，需要重新处理。此外，很多平台注释文件里“Gene Symbol”列会有---表示未知基因，注意把它过滤掉。
3. limma差异分析完整实操
3.1 构建分组、设计矩阵与对比矩阵
构建设计矩阵是limma分析的核心环节，这一步错了后面全错。用model.matrix()创建设计矩阵，公式写法取决于你的实验设计。最常见的是两分组比较，公式是~ 0 + group，其中group是因子。为什么要用0 +这种写法而不直接写成~ group？区别在于前者建模的是每组自己的均值，对比时可以直接提取两组的差值；后者建模的是基线组和差值，对比时反而麻烦。
分组时还有一个关键点：group因子的level顺序决定了对比的方向。比如你想比较“疾病组相对于正常组上调”，就要把正常组设为第一个level，用levels(group) <- c("normal", "disease")，确保后续对比是disease减normal。方向搞反了，整个差异列表的上下调就会颠倒。
library(limma) # 构建分组因子，normal在前，disease在后 group <- factor(c("disease", "normal", "disease", "normal"), levels = c("normal", "disease")) design <- model.matrix(~ 0 + group) colnames(design) <- levels(group) # 构建对比矩阵：disease组相比normal组 contrast_matrix <- makeContrasts(disease - normal, levels = design)3.2 lmFit、eBayes、topTable三步走的核心逻辑
limma的主程序就是三行代码：lmFit()拟合线性模型，eBayes()做经验贝叶斯统计检验，topTable()输出结果。每步都有值得留意的细节。
lmFit()输入表达矩阵和设计矩阵，表达矩阵必须是数值矩阵，行是基因，列是样本。这里有个易错点：表达矩阵必须是完整矩阵，不能有NA值。有时候GEO数据自带缺失值，需要提前用impute包补全，或者删掉有缺失的行，否则lmFit会报错。
eBayes()做的是对每个基因的统计量进行经验贝叶斯修正，这一步输出的就是最终用于筛选的p-value和调整后的p-value（adj.P.Val）。默认用BH方法做多重检验校正，控制FDR。对于差异基因筛选，我看大家常用adjusted p-value < 0.05和|logFC| > 1这两个标准，但也要根据数据实际情况调节，比如样本太少时p值很难小于0.05，可以适当放宽。
topTable()是输出结果的函数，可以指定排序方式、筛选阈值、输出数量。coef参数要指定你想查看的对比，默认是第一个。这个参数经常被忽略，导致多对比时输出错了列。
# 拟合线性模型 fit <- lmFit(expr_matrix, design) # 应用对比矩阵 fit2 <- contrasts.fit(fit, contrast_matrix) # 经验贝叶斯调整 fit2 <- eBayes(fit2) # 输出全部基因的结果，按调整后p值排序 results <- topTable(fit2, coef = 1, number = Inf, sort.by = "P") # 筛选差异基因 deg_list <- results[results$adj.P.Val < 0.05 & abs(results$logFC) > 1, ]3.3 差异基因结果的解读与输出
拿到差异基因列表后，先别急着画火山图，先自己检查一下结果是否合理。我习惯先看上下调基因的数量比例。比如刺激实验，通常会上调基因和下调基因数量不会差太多；如果是敲除实验，可能一个方向的基因会特别多。如果结果显示几千个上调、十几个下调，先怀疑一下分组方向是否搞反了，或者数据归一化是否有问题。
输出方面，除了差异基因的表格，通常还要保存完整的结果表，方便后续做注释、富集分析、或者补充材料。推荐输出CSV或TSV格式，保留所有列信息：logFC、AveExpr、t、P.Value、adj.P.Val、B。基因Symbol最好也绑到结果表里，方便后续Human Mine或者Enrichr直接使用。
# 绑定基因名 final_result <- data.frame(gene = rownames(results), results) # 写出完整结果 write.csv(final_result, "all_results.csv", row.names = FALSE) # 写出差异基因列表 write.csv(deg_list, "deg_filtered.csv", row.names = FALSE)4. 常见问题与排查技巧实录
4.1 表达数据是否该log2转换
这是新手最容易困惑的问题之一。判断方法其实很简单：看表达值的中位数或最大值。芯片数据经过RMA标准化后通常是log2值，范围大致在2到15之间，中位数在6到8；如果是RNA-seq的count数据，数值会很大，动辄上千上万。limma针对芯片数据直接处理log2值，如果是count数据，需要先做log2转换或者用limma的voom功能。
有次我用一个GEO数据集，表达矩阵里出现了负数，一开始以为是数据坏了，后来才发现是平台本身用log2转换后又有中心化处理，负值代表低于平均水平。这种数据直接拿去做差异分析没问题，但画热图和展示时要用原始值做细节调整。关键是别拿到一个矩阵就开始跑，先看看数值范围，心里有数。
4.2 探针注释失败、多个探针对应同一基因怎么办
探针注释失败的情况有几种：一个是芯片平台太老，注释文件里大量探针没有基因映射；另一个是我用merge合并表达矩阵后，发现行数比原来少了非常多，许多探针被过滤掉了。遇到这种情况，先检查注释文件和表达矩阵的行名格式是否完全一致，有些平台注释里探针ID有后缀，比如“_at”结尾，而表达矩阵里没有，就需要预处理对齐。
多对一基因的处理，除了之前说的取最大值，还可以用aggregate()函数按基因名聚合。我实际用下来，aggregate()配max或mean效果都稳定。取最大值的理由是这个基因在该样本的真实表达水平应该由其活性最高的转录本代表，这个逻辑在绝大多数情况下成立。
# 按基因名聚合，每基因保留最大表达值 expr_by_gene <- aggregate(expr_with_gene[, -1], by = list(gene = expr_with_gene$`Gene Symbol`), FUN = max) # 将第一列基因名设置为行名，删除原列 rownames(expr_by_gene) <- expr_by_gene$gene expr_by_gene$gene <- NULL4.3 批次效应与样本量少等坑
批次效应是公共数据做差异分析无法回避的问题。不同时间、不同实验室、不同芯片批次做出来的数据存在系统性偏差，如果分组刚好和批次完全重合，差异分析结果会严重失真。通常在下载数据时就要留意样本信息的“characteristics”列，看看是否有“batch”“date”“lab”等信息。如果发现批次混杂，可以用sva包的ComBat函数做批次效应校正，但校正前要保证设计方程正确，否则可能误删真实的生物学差异。
样本量少的情况也很常见，两个组各自只有3个样本，统计效力本来就弱，差异基因数量少是正常现象。前阵子我拿到一个数据集，每组就3个样本，用默认阈值筛完只剩下十几个基因。后来我把logFC阈值放宽到0.5，p值阈值放宽到0.1，再结合表达量绝对值做筛选，总算得到一批可用的候选基因。这种做法适合做初步探索，发论文投稿时还是要按严格标准来。
# 批次效应校正示例，假设batch信息已知道 library(sva) expr_combat <- ComBat(dat = expr_matrix, batch = batch_info, mod = design)5. 一些实操中的心得体会
用limma做GEO差异分析，跑通流程只是第一步，最花时间的往往不是代码而是数据清洗。我常跟朋友说，差异分析的成败在现场看不到，全在数据准备阶段。GEO数据库里的数据质量参差不齐，有些上传者把分组信息写得乱七八糟，有些芯片平台注释信息缺失严重，这些情况都需要在分析前发现并处理好，而不是硬着头皮往下跑。
我个人习惯是每次拿到新数据集，先建一个文件夹，把原始下载的series matrix文件、GPL注释文件、代码脚本、输出结果全部归类保存。别小看这个习惯，生信分析反复迭代很正常，有完整的数据目录结构，复盘和修改时能省下一大堆时间。另外，所有分析脚本开头我都会加上日期和数据集编号的注释，方便几个月后再回来时快速回忆当初的思路。
最后说一个很多人都忽略的细节：GEO数据集的引用问题。使用别人的数据做二次分析，记得在文章里引用原始数据的文献和GEO编号，这不仅是学术规范，也是数据共享的伦理要求。分析时顺手把数据集的基本信息，比如平台、样本量、分组、下载时间记录下来，写文章时直接就能用。
本文还有配套的精品资源，点击获取
