1. 差异基因分析前传从原始测序数据到表达矩阵的必经之路很多刚接触转录组的人会有一个误解拿到测序公司的报告里面有现成的基因表达量表格直接拿来做差异分析就行了。实际上那份报告的生成过程比你想象的要复杂得多而筛选出差异基因表达矩阵这件事误差恰恰容易在最不起眼的环节悄悄埋下。先说整条链路下机得到的FASTQ原始数据要先做质量评估然后比对到参考基因组或转录组再按基因或者转录本级别统计read数最后才得到我们熟悉的表达矩阵count matrix。这个矩阵的行是基因列是样本每个格子里面是一个非负整数代表这个基因在这个样本里检测到的reads数量。后续所有差异分析、富集分析、可视化全部建立在这个矩阵之上。换句话说表达矩阵是原料原料出了问题后面做出来的菜再好看也是废的。所以第一个关键决策点就是用什么工具比对和定量。目前主流的方案有两种一种是STAR比对到参考基因组再用featureCounts统计基因水平的read数另一种是直接用Salmon或者kallisto这类伪比对工具跳到转录组层面定量拿到转录本水平的TPM或者read count后再用tximport汇总到基因水平。我自己用得最多的是STAR featureCounts这条路线原因是它对后续差异分析的兼容性最好尤其是做基因水平分析的时候。STAR的映射速度快内存占用虽然在比对工具里算偏高的一个人类基因组索引大概需要30GB左右内存但换来的是极高的比对准确性。featureCounts可以非常方便地从BAM文件里统计出每个基因的read计数同时支持链特异性文库的设定这一步很容易被忽略。这里就踩过一个大坑链特异性参数设错了也不行。如果文库是链特异性的比如dUTP法建库这是目前最主流的建库方式之一但你用featureCounts统计的时候没有加-s 2表示read1和基因反义链对应那统计出来的计数会有相当一部分比对到错误的链上表现为反义链的read数显著高于正常值最终导致基因为一个基因的正反链计数几乎一样后续差异分析的结果自然会出现大面积失真。还有一个容易忽略的参数是多比对reads的处理方式。featureCounts默认只统计quality score足够高的唯一比对reads这没什么问题。但如果你的物种是多倍体或者基因组里有高度相似的重复区域那大量reads会被丢弃这时候就需要评估一下unique mapped reads的比例。我习惯在建完索引之后先跑一个样本做快速检查看一眼STAR日志里的Uniquely mapped reads percentage如果低于75%就要考虑是否是参考基因组版本不合适或者是转录组里重复序列比例太高的问题。表达矩阵拿到了接下来还有一个逃不掉的问题要不要做基因注释过滤。参考基因组的GTF注释文件里有大量的非编码RNA、假基因、线粒体基因这些基因在标准的mRNA-seq差异分析里通常是要根据研究目的过滤掉的。比如线粒体基因占比过高往往代表细胞状态有问题rRNA残留会导致某些核糖体蛋白基因的计数异常偏高。我会在拿到count矩阵之后根据GTF的gene_biotype字段只保留protein_coding和部分lncRNA具体保留哪些取决于文章的研究对象——如果课题本身就是研究lncRNA的那就另当别论。另外表达矩阵的行名问题看似不起眼却经常在后续分析里造成不可挽回的错误。不同来源的GTF文件基因ID有可能是Ensembl ID如ENSG00000141510、也有可能是Symbol如TP53还有可能是NCBI的Gene ID。一旦某个环节混用了ID类型后续的注释、富集分析都会对不上号。我的做法是从头到尾统一用Ensembl ID作为唯一标识所有转换都集中在最后做注释时才进行这样即使SQL脚本、R脚本来回切换也不容易出错。走到这一步原料基本就准备好了。但真正的差异基因分析是从这里才开始的。2. 差异分析工具选型DESeq2、edgeR、limma-voom到底怎么选很多人总爱问哪个差异分析工具最好其实这个问题本身就不严谨。转录组数据是一种典型的计数数据分布高度偏离正态而且不同基因之间的表达量基数差异极大——有的基因平均表达量上万有的只有几十甚至几个read。直接拿普通的t检验去处理这样的数据假阳性会高得离谱。所以差异分析工具的核心任务是把这种计数数据的分布特征合理地建模从而在有限样本量下给出尽量可靠的统计推断。目前最主流的三个工具分别是DESeq2、edgeR和limma-voom。它们在统计思路上有本质的区别。DESeq2用的是负二项分布模型Negative Binomial它假设每个基因的read count符合均值与方差相关的负二项分布然后用shrinkage收缩估计的方式估计基因的离散度参数。这个设计的好处是当你样本量较小、基因表达极差大的时候它会借助所有基因的离散度信息来借用信息把方差估计拉向一个合理的整体水平从而减少假阳性。坏处是它对每个基因都要拟合两层广义线性模型计算速度相对较慢但目前数据量下通常几分钟到十几分钟就能跑完不算瓶颈。edgeR的底层模型也是负二项分布但它的标准化方式用的是TMMTrimmed Mean of M-values整体思路比DESeq2更轻量。它更适合做样本数量较多的差异分析也适合基因数量较大时的快速筛选。如果你的项目里样本数特别多超过20个或者你只想做一个快速的初步筛选edgeR是个好选择。limma-voom走的是另一条路。limma本身是为芯片数据设计的线性模型工具voom算法则把RNA-seq的counts转换成log2-CPM值同时估计每个观测值的均值-方差关系并把这个精度权重融入到线性模型里。这套组合拳的优势是极其稳定尤其是在做复杂实验设计多个因素交叉的时候limma的线性模型框架处理起来比DESeq2和edgeR都更自然。此外limma速度极快几千个样本跑起来也不吃力。我自己的选型经验是这样的场景推荐工具理由经典两组比较对照 vs 处理每组3-5个重复DESeq2收缩估计让少量样本的离散度估计更稳定样本量大每组10个以上想快速出结果edgeRTMM标准化速度快结果稳定有多个因素如时间、剂量、组织等的复杂设计limma-voom线性模型框架处理多因素最舒服交互项也好做单细胞RNA-seq的差异分析pseudobulkDESeq2对低表达基因的收缩处理更好这些工具有一点是共通的它们的输出里都包含log2FoldChange对数倍数变化、pvalue和padj校正后的p值而差异基因表达矩阵正是从这些结果里筛选出来的。不过在真正运行之前还有一个绕不开的预处理步骤标准化。原始count矩阵的每个样本测序深度不一样总reads数可能差两三倍直接比较的话测序深度大的样本看起来所有基因都上调了这显然是错误结论。DESeq2内部通过median-of-ratios方法计算了size factor来自动处理测序深度差异edgeR用TMMlimma-voom则用CPM标准化。所以在跑差异分析的时候不需要也不应该在跑之前手动把count值转换成TPM——因为这些工具本身就对count数据做了标准化处理你手动转了一道反而会破坏负二项模型的计数分布假设。选择工具之后还有一个细节值得单独说说因子水平的设定。R语言里因子水平的顺序会影响结果的解读方向。比如你的处理组叫Treat对照组叫Ctrl如果R默认把Ctrl作为第一个水平那么DESeq2默认输出的log2FoldChange就是Treat相对Ctrl的变化但如果因子的水平顺序反过来log2FoldChange的符号也会反过来火山图里上调下调的方向就完全颠倒。这个问题我在给人debug时遇到过好几次明明代码每一步都对结果一查发现是对照组的水平排到后面去了。把工具确定下来接下来才真正进入到筛选差异基因的核心环节。3. 筛选差异基因的核心流程阈值设定、结果解读、可视化三板斧我见过不少人在拿到差异分析结果表之后随手设一个padj 0.05、|log2FC| 1就完事了。这个流程本身没错但如果你不理解这两个阈值背后的含义后续就会被各种伪差异基因折磨。先说padj这个概念是怎么来的。对一万多个基因分别做假设检验每个基因都会得到一个p值。p值的含义是在原假设为真基因无差异的情况下观察到当前这么大的差异甚至更大差异的概率。假设你设置的显著性水平是0.05那么即使所有基因都没差异每一万个基因里也大约有500个会因为你运气不好而出现p0.05的结果。这就是多重假设检验问题。为了解决这个问题Benjamini和Hochberg提出了一个叫BH校正的方法通过控制FDRFalse Discovery Rate错误发现率把所有p值按大小排序后逐步调整得到padj。padj的含义不再是单个基因出错概率而是在所有你判定为显著的基因里预期错误的比例。比如你筛出了1000个padj0.05的基因预期其中50个其实是假阳性。这是一种全局的风险控制比单个p值更可信。|log2FC| 1的意思则是倍数变化fold change不低于2倍。因为log2FC 1代表表达量相差2倍log2FC 2代表相差4倍。设这个阈值是为了过滤掉那些统计学显著但生物学效应微小的基因。如果一个基因的padj只有0.001但log2FC只有0.1说明它只是测序深度高、统计功效强真实效应可能没有实际意义。但是请注意这两个阈值并不是铁律。在某些场景下需要更严格比如你做全基因组筛选研究的又是一种表达变化幅度普遍很小的调控因子那么用|log2FC| 1可能一个基因都筛不出来。这时候可以放宽到|log2FC| 0.5甚至结合具体情况单独看padj。反过来如果你的研究是临床标志物筛选要求结果可重复性极高那可以把padj设成0.01甚至更严格。筛选阈值应该服务于你的生物学问题而不是反过来。跑完差异分析之后R里面通常是这样操作library(DESeq2) # 读取count矩阵 countData - read.table(gene_counts.txt, header TRUE, row.names 1) # 样本信息保证列名与count矩阵的列名对应 colData - data.frame( condition factor(c(Ctrl,Ctrl,Ctrl,Treat,Treat,Treat), levels c(Ctrl,Treat)) ) dds - DESeqDataSetFromMatrix( countData countData, colData colData, design ~ condition ) # 过滤保留至少在3个样本里总read数大于10的基因 keep - rowSums(counts(dds) 10) 3 dds - dds[keep, ] dds - DESeq(dds) res - results(dds, contrast c(condition, Treat, Ctrl)) res_df - as.data.frame(res) # 添加基因ID列 res_df$gene_id - rownames(res_df) # 筛选差异基因 sig_genes - subset(res_df, padj 0.05 abs(log2FoldChange) 1)筛选完之后最重要的一步是检查结果是否合理。我通常会做三个检查第一差异基因的数量和上下调比例。如果上调基因数量是下调基因的十倍先不要激动排除一下是不是某组样本质量有问题比如处理组细胞污染、样本标签搞反等等。第二已知的阳性对照基因是否出现在名单里。比如你做某种药物处理如果这个药物公认会上调某个关键通路基因而这个基因没在差异基因列表里那大概率是实验环节出问题了这时候无论统计结果多么漂亮都要回去查数据来源。第三主成分分析PCA的可视化结果。在DESeq2里用vst或者rlog变换后的数据做PCA如果不同组的样本能明显分开说明组间差异显著大于组内差异差异分析的结果就相对可信如果PCA图上两组样本混在一起那即使跑出一堆差异基因也很可能是批次效应而非真实生物学差异。拿到差异基因列表之后接下里就是差异基因表达矩阵的真正落地环节了。很多人会忽略这一环纯粹输出一份筛选结果表还不算完还要把这组基因在样本里的表达情况整理成方便可视化和下游分析的矩阵格式。这一步的操作逻辑是这样的从原始的vst标准化矩阵或者TPM矩阵中用差异基因的行名去子集得到一个新的小矩阵行是显著差异基因列是样本。这个矩阵才是后续画热图、做聚类分析、跑富集分析的输入文件。vsd - vst(dds, blind FALSE) vsd_mat - assay(vsd) # 用差异基因子集表达矩阵 deg_exp - vsd_mat[rownames(sig_genes), ] # 将矩阵排序并存储 deg_exp_sorted - deg_exp[order(sig_genes$log2FoldChange, decreasing TRUE), ] write.csv(deg_exp_sorted, DEG_expression_matrix.csv)这个deg_exp_sorted就是我们要的东西。它不但是热图绘制的直接输入也是后续GSEA、GO/KEGG富集分析、蛋白互作网络构建等一系列下游分析的起点。4. 常见坑点与排查思路从批次效应到多重检验的完整链路我在前面提到PCA可以用来检查样本分组但PCA还有一个更重要的作用发现批次效应。批次效应是转录组项目里最常见、也最致命的问题之一。简单来说如果你的3个对照组样本在一个批次跑完、3个处理组样本在另一批次跑完那么样本间的差异可能更多来源于不同批次的技术误差而不是真实的生物学差异。排查批次效应最直观的方法是看PCA图的前两个主成分。如果PC1主成分1能完美分开对照组和处理组同时PC2或PC3上还出现了一批按上机日期或者送样时间聚类的样本那就要高度怀疑批次效应的存在。还有一种情况更隐蔽两组样本本身生物学差异很大比如野生型和基因敲除型这时候第一批次的整体效应会叠加在真实差异之上PCA图看不出来但不代表批次效应不存在。针对批次效应的处理有几种常规手段。最理想的情况是实验设计时就把不同条件随机打散到不同批次俗称随机化但这在现实中经常做不到。如果已经做了能补救的标准方案是用limma包的removeBatchEffect函数传入已知的批次变量来校正数据或者用RUVSeq、sva这类工具先估计隐藏的因子再在差异分析模型里把这个因子作为协变量放进去。我实际跑过的最典型的坑是在DESeq2的设计公式里漏掉了批次变量。比如设计公式只写了~ condition但样本实际上分了几个批次结果差异分析把批次效应和生物学效应混在一起导致大量假阳性基因——尤其是那些对批次敏感的低表达基因。后来把设计公式改成~ batch condition差异基因数量从7000多降到了1000多富集结果也跟文献报道更吻合了。这件事给我的教训是任何与生物学处理无关的已知技术变量都应该尽可能纳入设计公式。再说一个容易被忽视的低表达基因过滤问题。很多人一上来就直接跑DESeq2不过滤任何基因。这时候大量在所有样本里读数都是0或者只有1~2的基因也会参与离散度估计和多重检验校正。这些基因不仅会增加多重检验矫正的负担padj会变得更保守还会在计算离散度时引入噪音导致真实表达、稳定表达的基因反而被误伤了。那到底怎么过滤通常的做法是保留至少在最小样本数比如每组最少样本数里count数达到一定阈值比如10或者1的基因。注意这个过滤用的count数一定要用原始的整数count而不是TPM或者CPM。原因是在模型假设里count服从负二项分布低计数的基因信息量本来就少过滤掉它们不至于损失太多真实的生物学信号。如果用TPM过滤大量本身转录本长度很短的基因会被错误地丢掉——因为TPM对转录本长度做了归一化短基因的TPM天然偏低但这不代表它不表达。接下来是多重检验校正顺序的问题。R里results()出来的padj是已经用BH方法校正过的如果你自作聪明地对pvalue再调一次p.adjust会把结果弄得更严格或者更宽松——这取决于你用的方法。我见过有同学在拿到res之后发现padj列全是NA然后自己重新用p.adjust(res$pvalue, method bonferroni)跑了一遍结果把差异基因筛得一个不剩最后跑来问我为什么结果为空。一问才知道他把已经被独立筛选过一次的p值又做了一次非常保守的Bonferroni校正。这不是工具问题而是对统计流程的理解问题。还有一组需要单独拎出来说的基因那些padj为NA的表达量极低基因。在DESeq2里如果一个基因在所有样本里的count都很低它可能没有足够的数据来估算离散度对应的padj就会显示为NA。在筛差异基因时千万不要把含NA的阈值判断直接拿去比较——R里NA 1的结果是NA这会让你的子集结果莫名其妙少掉大量的行。安全的做法是先过滤掉padj是NA的行再做筛选。低样本量下的p值分布偏移问题也值得一提。当每组只有2~3个样本时单基因的p值分布往往不服从均匀分布而是偏向0或偏向1——这是因为极度缺乏自由度导致的检验功效不足。这种情况下DSSeq2的收缩估计虽然能在一定程度上缓解但仍要警惕如果你在每组只有2个样本、甚至没有生物学重复的情况下硬跑差异分析得到的p值分布会非常糟糕很多p值极小的基因其实只是在两组样本间恰好差异大而已没有可重复性可言。这也是为什么很多期刊现在强制要求至少3个生物学重复作为RNA-seq实验的最低门槛。最后要说的排查手段是把自己钉在结果的逻辑合理性上。跑完差异分析后脱离开统计术语回到你的实验本身去问三个问题这个处理在生物学上到底应该影响哪些通路这些通路里的代表性基因在结果里出现了没有结果里有没有哪些基因跟已知机制完全相悖这些问题的答案往往比任何统计指标更能帮你发现深层的数据问题。有一次我跑完一个药物浓度梯度的转录组数据差异基因列表里竟然出现了大量的糖酵解相关基因上调而原课题做的是免疫调控怎么也想不通。后来重新检查才发现处理组的细胞培养时间比对照组长了整整一天细胞密度过高导致培养基缺氧糖酵解通路是被环境逼出来的。这个例子说明差异基因分析筛出来的东西很多时候反映的是你不知道的实验偏差这时候回到实验记录找原因而不是急着写结论。5. 差异基因表达矩阵的落地应用热图、富集、实验验证一条龙差异基因表达矩阵做出来之后实际工作中最常见的下游使用方式有三种可视化热图和火山图、通路富集分析、以及实验验证前的候选基因清单整理。每个环节的侧重点各不相同但都可以围绕这个矩阵展开。热图是展示差异基因表达模式最强有力的工具。绘制的时候选基因数量如果太多比如超过1000个整张热图会变成一团密集的马赛克看不出任何模式。这时候我通常会过滤掉那些在各个样本间表达量波动特别小的基因或者先按方差排序取top50~100画图。重要的一点是热图里的颜色映射要基于标准化后的表达量vst或者z-score不能直接拿原始count画因为不同样本的测序深度差异会让所有基因的颜色都偏暗或偏亮掩盖真实的差异模式。library(pheatmap) # 对每一行基因做z-score标准化 heatmap_mat - t(scale(t(deg_exp_sorted))) pheatmap( heatmap_mat, cluster_cols TRUE, cluster_rows TRUE, show_rownames FALSE, annotation_col data.frame(condition colData$condition, row.names colnames(heatmap_mat)) )热图只是看差异功能富集分析才是帮助你把差异基因列表转化成一串能被审稿人理解的生物学结论的关键步骤。常规做法是使用clusterProfiler做GO和KEGG通路富集分析。GO富集分析的输入通常是一组差异基因ID加上背景基因ID也就是你测到的所有基因而不是全基因组所有基因。这一点非常关键如果你不指定universe背景集富集分析会拿全部基因作为背景计算富集显著性测序技术天生只能检测到一部分基因这样算出来的富集p值会被严重高估出现大量本来不显著的伪富集结果。library(clusterProfiler) library(org.Hs.eg.db) # 转换ID sig_entrez - bitr(sig_genes$gene_id, fromType ENSEMBL, toType ENTREZID, OrgDb org.Hs.eg.db) # GO富集分析 ego - enrichGO( gene sig_entrez$ENTREZID, universe background_entrez$ENTREZID, OrgDb org.Hs.eg.db, keyType ENTREZID, ont BP, pAdjustMethod BH, qvalueCutoff 0.05 )GO富集分析虽然通用但有个老毛病术语体系太庞大结果动辄几十上百条显著性通路读起来像天书。这时候我的经验是先用simpleEnrichment或者手动的冗余去除reduce把相似的GO term合并再做一层GSEA分析来辅助判断哪些通路里的基因是一致上调或一致下调的。GSEA和普通的GO富集不同它不要求你先设定一个差异基因列表而是直接拿所有基因的表达变化排序和foldchange信息计算通路级别的得分。这样能抓住那些单个基因变化不显著、但整个通路协同上调或者下调的生物学信号非常适合处理那种效应积少成多的通路。矩阵的下一个落地场合是实验验证这个环节往往是整个分析链条里最容易被低估的。差异分析出来一堆基因热图富集都很漂亮但真要落到下一步wet lab实验你不能把1000个基因都拿去做qPCR或者Western验证。这时候矩阵的价值就体现出来了你要从这个矩阵里挑出那些表达变化幅度大、组内一致性高、且在已有文献里有明确机制支撑的代表性基因而不是单纯挑log2FC最大的那几个。什么叫组内一致性高举个例子某个基因在三个处理样本的log2FC分别是2.0、2.1和0.9平均值1.67看起来显著上调但实际上第三个样本的数据明显偏低可能是样本处理不均匀或者操作误差导致的。如果这个基因恰好属于一个跟课题方向高度相关的通路你拿它去做qPCR验证十有八九会验证失败。反过来如果某个基因log2FC只有1.2但三个样本之间几乎无波动它在qPCR里重复出来的概率反而更高。在矩阵里检查组内一致性我通常会看每个基因在不同样本里的原始标准化表达值绘制箱线图或者直接计算每组内变异系数。还有一个实用技巧把差异表达矩阵直接作为验证实验设计的参考。你在设计qPCR验证引物或者挑选做免疫组化的抗体时先确认备选基因在你的矩阵里确实有足够的表达量比如vst标准化后的均值大于某个阈值而不是只盯着log2FC大小。表达量太低的基因在qPCR里CT值会超过35循环验证结果会因为技术噪音而变得完全不可信。说到这其实差异基因表达矩阵的下游应用是很灵活的。有位同事做过一个有意思的尝试他把筛选出的差异基因表达矩阵和转录因子结合位点数据库结合通过共表达网络分析找到了一批与核心转录因子共调控的模块基因然后把这些模块基因和临床患者的预后数据做了生存分析最后硬是从一个普通的细胞系转录组项目里挖出了三个新的预后标志物候选基因。这个例子很好地说明差异基因表达矩阵不仅仅是分析链路的终点它更像一个蓄水池后续只要你愿意随时可以从池子里引出不同的分析支流。6. 实操总结我踩过的坑希望你不用再踩一遍最后聊几句个人经验层面的东西不展开写大理论只说那些我最初不理解、后来被现实教育过才明白的细节。第一永远不要直接拿公司给的基因表达量表去做差异分析。很多公司给的表格里面对应的可能是FPKM、TPM或者已经做了某种标准化的值。FPKM/TPM作为跨样本比较的表达量没问题但它们不应该被输入到DESeq2或edgeR里做差异检验。这两个工具天然假设输入是整数count你把小数输入进去要么报错要么结果畸变。如果不幸你手上只有FPKM矩阵可以用tximport配合转录本定量结果重建count数据或者换用limma-voom流程来处理voom内部会对log-counts做拟合对这类数据的容忍度更高一些。第二log2FoldChange和Log2CPM是两个完全不同的东西不要混用。差异结果表里的log2FoldChange是指两个组之间标准化后表达量的对数比率正负号代表上下调方向而Log2CPM是各个样本里基因的表达量归一化值只用于展示和比较表达水平。曾经有人问我为什么某个基因log2FoldChange 2.5但热图上颜色不深仔细一看热图用的是原始count矩阵样本内其他基因的表达量太高把颜色分布整个压扁了。这个问题说大不大说小不小但在审稿人眼里这就是数据处理的硬伤。第三运行DESeq2之前仔细检查样本列名是否和colData的行名完全一致。R里面DESeqDataSetFromMatrix要求列名与colData的DataFrame行名精确匹配哪怕多一个空格、一个下划线函数都会报错。这个错误信息有时候并不直观最常见的报错是colData must have the same samples as colnames of countData。如果你看到这行报错先检查一下是不是样本名没对上而不是急着改代码。第四不要忽略基因注释版本的核对。同一个物种不同版本的GTF文件基因ID可能发生改变gene symbol改了Ensembl ID也可能有版本后缀比如ENSG00000141510.17和ENSG00000141510.10是不同的东西。在多个GTF版本之间混用时差异基因筛选结果的可靠性会大打折扣。最稳妥的操作是在项目开始时固定一个参考基因组版本和一个GTF版本后续所有比对、计数、注释都基于它完成绝不中途更改。写到这里差异基因分析从原始数据到差异基因表达矩阵的链路差不多就完整了。这个流程里没有一步是可以用全自动工具跳过的黑箱每一步都藏着值得仔细推敲的统计假设和生物学前提。我自己从最初只会跑DESeq2默认参数的小白到现在每次拿到结果先做三遍逆向检查中间踩过的坑远比这篇文章写的要多。如果你正在跑自己的第一批转录组数据我的建议是别急着追求一键汇总先把表达矩阵的质量检查、工具选型的理由、阈值设定的依据想清楚后面再复杂的结果也只是水到渠成的事。