1. RNA-seq差异表达分析入门指南
第一次接触RNA-seq数据分析的朋友可能会被各种专业术语吓到，但别担心，我刚开始做生信分析时也是这样。简单来说，RNA-seq差异表达分析就是找出不同实验条件下表达量有显著差异的基因。比如比较癌症组织和正常组织，找出哪些基因在癌症中异常活跃或沉默。
limma、Glimma和edgeR这三个R包是处理RNA-seq数据的黄金组合。limma擅长线性建模，edgeR精于计数数据建模，而Glimma则让结果可视化变得生动直观。这三个工具配合使用，就像厨房里的刀、锅和铲子，各司其职又能完美配合。
在实际项目中，我经常遇到这样的问题：测序数据拿到了，但不知道从何下手。这时候一个清晰的流程就特别重要。完整的分析通常包括：数据预处理（质量控制、标准化）→ 差异表达分析 → 结果可视化 → 功能注释。今天我们重点讲解最核心的差异表达分析环节。
2. 数据准备与预处理
2.1 原始数据质控
拿到原始测序数据后，我习惯先用FastQC做个快速检查。这个工具会生成一份详细的质检报告，包括测序质量分布、GC含量、重复序列比例等关键指标。有一次分析肿瘤样本时，就发现某个样本的GC含量异常高，后来证实是样本污染。
处理RNA-seq数据时，edgeR的DGEList对象是我们的起点。这个对象不仅存储原始计数矩阵，还能记录样本分组信息和基因注释。创建方法很简单：
library(edgeR) counts <- read.delim("raw_counts.txt", row.names=1) groups <- factor(c("Control","Control","Treat","Treat")) dge <- DGEList(counts=counts, group=groups)2.2 表达量标准化
RNA-seq数据有个特点：不同样本的测序深度可能差异很大。就像比较两个水杯里的水量，必须先确认杯子大小是否相同。edgeR的calcNormFactors函数采用TMM方法进行标准化，这种方法对表达量差异大的基因不敏感，特别适合真实生物数据。
dge <- calcNormFactors(dge)标准化后的数据通常转换为log2-CPM值（每百万计数的log2值）。这个转换使数据更接近正态分布，也方便后续比较。我常用下面这段代码：
lcpm <- cpm(dge, log=TRUE, prior.count=2)3. limma差异表达分析实战
3.1 设计矩阵构建
limma分析的第一步是创建设计矩阵（design matrix），这相当于给实验"画图纸"。比如我们比较三种细胞类型（Basal、LP、ML）的基因表达差异，同时考虑测序批次效应：
design <- model.matrix(~0 + group + lane) colnames(design) <- gsub("group", "", colnames(design))这里~0表示不设截距项，这样后续的对比矩阵设置会更直观。有一次我忘记加0，结果对比时各种系数对不上，调试了好久才发现问题。
3.2 均值-方差关系校正
RNA-seq数据的方差通常与均值相关，这违反了线性模型的同方差假设。voom函数就是解决这个问题的神器，它会为每个基因计算权重：
v <- voom(dge, design, plot=TRUE)voom生成的图特别重要，它能告诉我们数据质量如何。好的数据应该呈现典型的"喇叭形"趋势线。如果左端突然下降，说明低表达基因过滤不够；如果整体太平，可能样本间生物学差异太小。
3.3 线性模型拟合
有了权重后，就可以用limma的经典三部曲进行差异分析：
fit <- lmFit(v, design) cont.matrix <- makeContrasts(BasalvsLP=Basal-LP, levels=colnames(design)) fit2 <- contrasts.fit(fit, cont.matrix) efit <- eBayes(fit2)eBayes这一步采用了经验贝叶斯方法，借用所有基因的信息来更准确估计单个基因的差异。这就像班级考试后老师根据全班表现调整评分标准，避免对个别同学打分过于极端。
4. 结果解读与可视化
4.1 差异基因筛选
使用topTreat函数可以查看差异最显著的基因。我一般会同时关注p值和logFC（差异倍数对数）：
topGenes <- topTreat(efit, coef=1, n=Inf) sigGenes <- topGenes[topGenes$adj.P.Val < 0.05 & abs(topGenes$logFC) > 1, ]在实际项目中，我发现单纯依赖p值可能会漏掉一些生物学意义重要的基因。比如某个基因p值=0.051，但logFC很大，可能值得进一步验证。
4.2 交互式可视化
静态图已经不能满足现代研究需求了，这就是Glimma的用武之地。它生成的交互式MD图可以在网页中自由探索：
library(Glimma) glMDPlot(efit, coef=1, status=dt, main=colnames(efit)[1], counts=lcpm, groups=group, launch=TRUE)最近一次项目汇报中，我用Glimma做的图让合作生物学家特别兴奋，因为他们可以直接点击感兴趣的基因查看表达模式，不用再问我"那个红色点是什么基因"了。
4.3 热图展示
热图是展示差异基因表达模式的经典方式。我常用pheatmap包，因为它自动优化的聚类和标注效果很棒：
library(pheatmap) select <- order(abs(efit$coefficients[,1]), decreasing=TRUE)[1:50] pheatmap(lcpm[select,], scale="row", show_rownames=TRUE, cluster_cols=TRUE, annotation_col=annot)有个小技巧：热图前对行进行缩放（scale="row"），这样不同基因的表达量可以放在同一尺度比较，红色代表高表达，蓝色代表低表达。
5. 高级分析与注意事项
5.1 基因集富集分析
差异基因列表出来后，下一步就是理解它们的生物学意义。camera基因集检验是个不错的选择：
library(limma) idx <- ids2indices(Mm.c2, id=rownames(v)) cam <- camera(v, idx, design, contrast=cont.matrix[,1]) head(cam, 5)这个方法考虑了基因间的相关性，比传统的超几何检验更可靠。我最近分析阿尔茨海默症数据时，就通过它发现了突触功能相关通路显著富集。
5.2 批次效应处理
当发现样本明显按批次而非实验条件聚类时，就需要考虑批次校正。ComBat是个常用选择，但要注意不要过度校正：
library(sva) batch <- factor(c(1,1,2,2)) modcombat <- model.matrix(~1, data=pheno) combat_edata <- ComBat(dat=lcpm, batch=batch, mod=modcombat)曾经有个项目，校正后生物学差异反而消失了，后来发现是因为批次和实验条件完全混杂。这种情况只能重新设计实验。
5.3 常见问题排查
新手最容易犯的错误包括：
- 忘记过滤低表达基因（建议保留CPM>1的基因）
- 设计矩阵设置错误（一定要检查contrast矩阵）
- 忽略voom权重直接建模
- 多重检验校正方法选择不当（RNA-seq推荐用BH方法）
我习惯在分析脚本中加入检查点，比如：
stopifnot(ncol(design)==ncol(cont.matrix)) stopifnot(all(colnames(design)==rownames(cont.matrix)))这些小技巧能避免很多低级错误。
