高通量测序数据与差异表达分析概述
好，咱们正式开始。这一章我把高通量测序数据——特别是RNA-seq数据——的差异表达分析流程梳理一遍。你想想看，做这行的人几乎天天跟这东西打交道。我从头说，基础打牢，后面学edgeR和limma-voom才不吃力。
高通量测序数据长什么样？
说实话，我第一次接触原始测序数据时有点懵。一堆FASTQ文件，里面全是碱基和质量值。但我不建议大家直接碰FASTQ，除非你想从头搭流程。通常我们拿到的，是经过比对和定量之后的表达矩阵。长这样——行是基因，列是样本，格子里的数字是 counts。
说白了，就是一张计数表。我习惯叫它“原始计数矩阵”。
表达矩阵的结构
给你一个最简单的例子，人工造了点数据，你感受一下结构。
Sample1_Control Sample2_Treated Sample3_Control Sample4_Treated
GeneA 1523 842 1789 912
GeneB 234 1987 189 2103
GeneC 4567 3210 4987 3456
GeneD 89 567 78 612看见没？每行一个基因，每列一个样本。数字是比对到这个基因上的reads数目。嗯，这些数字看起来挺整齐的，但背后的问题可不少。后面我会慢慢讲。
差异表达分析的核心逻辑
我自己的理解是这样的：差异表达分析本质上是个统计推断问题。给定了两组（或多组）样本的表达量，你问基因在这两组间有没有显著变化。很多人以为直接算个倍数变化，大于2就是上调，小于0.5就是下调。错，完全不是那么回事。
实际流程通常包含这几步：
- 数据预处理——过滤低表达基因、标准化
- 统计建模——用负二项分布或准似然模型拟合离散度
- 假设检验——计算p值和校正多重检验
- 结果注释与可视化——火山图、热图、富集分析
关键点：差异表达不是只看倍数变化，还得看统计显著性。我曾经有个项目，某个基因倍数变化高达10倍，但p值0.6，后来发现是样本量太少加技术重复导致的假象。从那以后我再也不单看fold change了。
为什么不能直接做t检验？
刚入行的人经常问这个问题。常规t检验假设数据服从正态分布，并且方差齐性。RNA-seq的count数据呢？它是离散的，低表达基因方差很大，高表达基因方差又相对稳定。尤其重要的是——RNA-seq数据有明显的“均值-方差”关系，说白了你不能直接用t检验那套东西。
我做个表，把几种常见模型比一下，你一看就明白。
| 模型 | 适用数据类型 | 核心假设 | 典型工具 |
|---|---|---|---|
| 泊松分布 | 计数数据 | 均值=方差 | 早期方法 (已淘汰) |
| 负二项分布 | 过离散计数数据 | 均值<方差（有离散度参数） | edgeR, DESeq2 |
| 准似然 | 灵活，对分布不敏感 | 仅要求一阶矩正确 | limma-voom, edgeR QL |
看到最后一行没？limma-voom就是走准似然路线。说白了它比纯负二项更灵活一些，尤其当你的样本量不大或者有批次效应的时候，voom转换完再做经验贝叶斯，效果挺好。
批次效应是什么？为什么需要消除？
这个问题我得好好说一下。我在项目中遇到过很多次——样本明明处理方式一样，只是分了批次上机测序，结果差异基因全跟着批次走，跟处理条件没半毛钱关系。
注意：批次效应不是异常值，它是系统性偏差。比如两个批次的文库制备用了不同试剂盒，或者不同天跑测序仪，都会引入批次效应。不消除它，你的差异分析结果可能就是垃圾。
我个人的经验是：实验设计阶段考虑批次，比后期用什么算法都重要。能随机化就随机化，不要把全部对照组放在Batch1，全部处理组放在Batch2——那叫混淆，神仙也救不了。
消除批次效应的常见策略
做法有好几种，我大致列一下：
- ComBat：基于贝叶斯框架，调整表达值
- limma中的removeBatchEffect：在线性模型中去掉批次变量
- RUVseq (Remove Unwanted Variation)：用内部对照基因估算“不想管的变异”
- limma-voom + 批次作为协变量：直接在模型中包含批次因子
咱们这个课程重点讲的就是最后一种——在limma-voom框架里把批次当作协变量放进去。为什么选它？因为voom本身把count数据转换成了近似正态的log2-CPM值，再加个线性模型里的批次项，干净利落，结果容易解释。
小技巧：如果你做探索性数据可视化（比如PCA、聚类图），可以用removeBatchEffect先去除批次对表达量的影响再画图。但在差异检验阶段，我个人习惯还是把批次放在模型里，而不是预先调整表达值——后者可能会导致统计检验的p值偏差。
本章知识框架
下面这张SVG图把本章的核心概念串了一遍。你看一遍，后面动手做的时候脑子里会清楚很多。
图里的逻辑很清楚：从原始FASTQ开始，拿到定量矩阵，预处理以后进入差异表达分析核心。三种主流路线，咱们选的是limma-voom加批次效应消除这条。最终输出差异基因列表——但要记住，这只是第一步，后面还要做功能注释和生物学验证。
怎么评估差异表达结果的好坏？
好问题。我在项目里常用的几项指标：
- p值分布：如果p值直方图在0附近有个尖峰，其余区域平坦，说明模型合理。如果p值整体往1偏移，可能模型有问题。
- MA图：差异基因应该大致对称分布在横轴两侧。如果全偏向一边，可能是标准化没做好。
- 主成分分析：先看主成分图，检查批次效应明显不明显。如果batch1和batch2把样本分开了，赶紧想办法。
经验之谈：我习惯拿到表达矩阵后，第一件事不是跑差异分析，而是先做PCA+热图。花10分钟看看数据长得怎么样。如果对照和处理组混在一起，先别急着找原因，检查样本标签有没有弄反——我吃过这个亏。
从RNA-seq到差异分析，我建议你记住这几条
嗯，总结一下吧，不啰嗦：
- counts数据不等于表达量——必须经过标准化（如CPM、TPM、RPKM）才能做样本间比较。
- 批次效应是真实存在的——不是你分析能力的问题，是实验设计的问题。接纳它，处理它。
- 没有完美的模型——edgeR、DESeq2、limma-voom各有适用场景。选工具前先搞清楚你的数据特点：样本量、离散度、批次数量。
- 差异分析只是过程，不是终点——最后你还得回答生物学问题。