笔记写到26.3.2了其实这个编号是我自己的一套生信学习笔记体系——第26章是转录组分析专题3.2节正好是差异表达分析的核心实操部分。前面刚把STAR比对、featureCounts定量这些上游步骤跑完这期笔记就接着记录从count矩阵往下走的过程如何用DESeq2找出组间显著差异表达的基因以及怎么用PCA图、MA图、火山图验证结果到底靠不靠谱。这套流程是每一个做RNA-seq相关的生信学习者都需要掌握的基础链路也是从“会跑命令”过渡到“会判断结果是否可信”的关键一步。这篇笔记同样假设你已经对Linux基本操作、R语言基础有大致了解不需要太深但至少要能跑起来。如果你的情况和我差不多都是零基础转生信、一个人折腾那这篇笔记应该能帮你把差异分析这条路走通顺带避开我踩过的几个坑。1. 先把整条链路对齐count矩阵到底是怎么产生、怎么被消费的1.1 从FASTQ到count矩阵上游这几步各干各的事标准的RNA-seq分析流程大概是FASTQ原始数据开始经过质量控制、比对或定量最后得到一张以基因为行、样本为列的表达矩阵。这张矩阵就是绝大多数下游分析的起点差异表达分析、聚类热图、富集分析全都是围绕它展开。我自己的流程是拿到测序下机的FASTQ文件后先用fastp做质控和接头去除这一步会把低质量碱基trim掉避免后续比对的时候出现大量没比上的read。然后做STAR比对把read定位到参考基因组上得到每个样本的BAM文件。最后用featureCounts基于GTF注释文件统计每个基因覆盖了多少条read输出count矩阵。这里有一个容易被新手忽略的点GTF注释文件的版本必须和参考基因组配套。比如你用了GRCh38.p13的参考基因组注释文件却用了GRCm39的featureCounts跑起来大概率会报染色体命名对不上甚至什么错都不报只是最终的count数明显异常。这类问题排查起来很费时间最好在建立参考基因组索引时就把注释文件一起整理清楚。1.2 两种定量路线的取舍STARfeatureCounts还是Salmon做定量这一步社区里现在主流有两个路线第一个路线是先比对到参考基因组再做基因计数STARfeatureCounts的组合是最常见的第二个路线是绕过基因组直接把read比对到转录本上定量代表性工具是Salmon和kallisto。两条路线从差异分析的角度说都能用但实际场景不太一样我整理了一个对比对比维度STAR featureCountsSalmon tximport是否生成BAM文件生成可后续查看IGV、排查污染不生成只输出定量结果依赖参考基因组版本是基因组和GTF必须配套否主要依赖转录本序列耗时与内存高尤其STAR建索引和比对阶段低一般几十分钟就能完成适合场景需要核查比对情况、看junction、做变异检测大样本量、服务器资源紧张下游接口featureCounts直接产出count矩阵通过tximport转为基因水平count矩阵我个人的习惯是样本量不大、需要仔细核查数据质量时用STAR这条线因为BAM文件能做的事情太多了——IGV里看某个基因的reads分布、检查是否存在样本污染、确认junction位点这些都是Salmon给不了的。但如果是几百个样本的表达谱项目跑STAR会相当吃资源Salmon的轻量优势就体现出来了。这里没有谁替代谁的问题纯粹是场景匹配。我这个项目样本不多就继续用STARfeatureCounts后面的分析全部围绕它产生的count矩阵展开。2. 样本信息表和数据合并一个字段顺序就能毁掉整个分析2.1 featureCounts的输出长什么样怎么批量合并featureCounts对每个样本输出一个表格文件默认列包括Geneid、Chr、Start、End、Strand、Length以及该样本的count列。如果一批样本跑完你会得到十几个这样的txt文件接下来要做的就是把它们合并成一张以基因为行、样本为列的矩阵。合并的做法不复杂但手动复制粘贴肯定是不行的。我习惯在R里批量读取只保留Geneid和count列然后按Geneid作为行名合并library(dplyr) sample_ids - c(s1,s2,s3,s4,s5,s6) count_list - list() for (id in sample_ids) { file_path - paste0(counts/, id, _counts.txt) tmp - read.table(file_path, header TRUE, skip 2, sep \t, check.names FALSE) tmp - tmp[, c(Geneid, id)] count_list[[id]] - tmp } merged - count_list[[1]] for (id in sample_ids[-1]) { merged - full_join(merged, count_list[[id]], by Geneid) } rownames(merged) - merged$Geneid merged - merged[, sample_ids, drop TRUE]注意我用了skip 2因为featureCounts输出文件前面有两行以井号开头的注释信息不跳过的话表头会错位。合并之后建议检查一下行数是不是和原来一致如果某个样本因为基因注释版本不同导致基因数量不一样full_join会把行数变多这时候就要回头查数据。2.2 样本信息表colData的构建细节有了count矩阵下一步就是构建DESeq2要求的样本信息表colData。这个东西看起来只是一个简单的data.frame但构建的时候有几个细节值得多说两句。第一行名必须和count矩阵的列名完全一致顺序也要一致。DESeq2在内部是按列位置来对应样本信息的只要顺序错一个后面所有结果全错而且不报错。第二除了实验分组强烈建议把批次信息也放进来。这里说的批次包括但不限于文库制备批次、测序上机流通池编号、RNA提取时间。哪怕你当时觉得“这些样本是同一天做的”只要不是同一个流池上机的就有批次效应的可能。后面第5章我会说一个我自己遇到的真实案例批次效应绝对值得一开始就放进colData里。2.3 因子水平的顺序一个容易忽略却致命的细节colData里分组那一列不能直接用字符串要转成factor而且要明确基准水平是谁。DESeq2比较方向的默认逻辑是“非基准水平 vs 基准水平”比如condition这一列有ctrl和treat两个水平如果你设置的levels是c(treat,ctrl)那么DESeq2默认比较的就是ctrl相对于treatlog2FoldChange的正负号会完全反过来。colData - data.frame( condition factor(c(ctrl,ctrl,ctrl,treat,treat,treat), levels c(ctrl,treat)), batch factor(c(B1,B1,B2,B1,B1,B2), levels c(B1,B2)) )我刚开始跑的时候就是没注意这个把levels写反了结果差异基因出来一大半是反的看起来还完全自洽。后来拿RT-qPCR验证的数据对照才发现符号不对。这个坑太容易被忽略了所以每次构建dds之前我都会先table(colData$condition)看一眼因子顺序。3. DESeq2跑差异分析核心代码与统计逻辑3.1 为什么差异分析普遍选择DESeq2现在做转录组差异表达分析最常用的工具不外乎DESeq2、edgeR和limma-voom。三者都能输出差异基因列表但统计模型不同。DESeq2用的是负二项分布模型它能天然处理生物学重复之间的过离散现象不需要人为去估算一个经验方差阈值。简单说真实生物学重复之间的波动往往比泊松分布预期的更大负二项分布多一个离散参数这个参数会在DESeq2内部自动估计对用户来说基本透明。工具统计模型标准化方式适用场景DESeq2负二项分布几何均值中位数比常规RNA-seq、重复数适中edgeR负二项分布TMM样本量大、需要精细的文库校正limma-voom均值-方差加权线性模型基于log2CPM数据量大、计算速度快对于大多数常规RNA-seq项目我的建议是优先DESeq2它是目前社区里最稳的选择。DESeq2把文库大小归一化、离散度估计、统计检验、多重检验校正全包在一个流程里默认参数在绝大多数情况下都合理不太需要手动干预。3.2 核心流程代码dds构建、design公式与results提取核心流程其实就几行R代码library(DESeq2) count_matrix - read.csv(merged_counts.csv, row.names 1) colData - data.frame( condition factor(c(ctrl,ctrl,ctrl,treat,treat,treat), levels c(ctrl,treat)), batch factor(c(B1,B1,B2,B1,B1,B2), levels c(B1,B2)) ) dds - DESeqDataSetFromMatrix(countData count_matrix, colData colData, design ~ batch condition) dds - DESeq(dds) res - results(dds, contrast c(condition, treat, ctrl)) summary(res)design公式这里要解释一下。~ batch condition的意思是在估计离散度和拟合模型时把批次效应作为一个协变量纳入再把condition作为主要的组间因素。batch放前面不意味着它更重要只是告诉模型先吸收掉这一部分技术变异再来评估分组差异。如果batch和condition是完全混杂的——比如所有处理组样本都在B2批次、对照组都在B1批次——那么设计公式里同时放batch和condition会出问题因为模型无法区分变异到底来自哪个因素。这种情况应该在实验设计阶段就避免如果已经踩了这种混杂的局面就要靠纯统计手段去尽力弥补但远没有一开始设计好来得好。跑完DESeq之后results()函数提取的就是组间比较的结果默认是control vs treatment注意看summary输出里的方向。结果表里每一行是一个基因每一列的意义下面单独说。3.3 结果表字段与差异基因筛选标准DESeq2的results表包含了多个字段刚开始看会有点懵但实际上每个字段都有明确的统计学含义字段含义实际用途baseMean该基因在所有样本中的平均归一化表达量低表达基因往往不可靠log2FoldChange组间表达差异的log2倍变化差异幅度正负号代表方向lfcSElog2FoldChange的标准误估计的稳健程度statWald检验统计量内部使用pvalue未校正的P值单基因统计显著性padjBH法校正后的P值多重检验校正后用于筛选筛选差异基因我习惯用两个阈值同时卡padj 0.05且|log2FoldChange| 1。padj是经过Benjamini-Hochberg多重假设检验校正后的P值你一次比较几万个基因如果不做校正会有大量假阳性出现所以一定不要只拿原始pvalue筛选。3.4 为什么推荐用lfcShrink收缩效应量有一件事我一开始总跳过后面才意识到很重要直接用results()拿到的log2FoldChange在小样本、低表达基因上往往被高估而且方差很大。社区里推荐的方案是再用lfcShrink()做一步效应量收缩。res_shrink - lfcShrink(dds, coef condition_treat_vs_ctrl, type apeglm)收缩的本质是把那些低信息量基因的log2FoldChange向0拉近信息量高、差异大的基因受影响很小。这样处理之后画MA图会好看很多排名的稳定性也更好。做差异基因筛选时我一般用收缩后的log2FoldChange做MA图展示用原始结果的padj做显著性筛选两者各取所长。4. 验证结果PCA、MA图和火山图的操作思路4.1 样本关系检查PCA图必须用归一化变换后的数据差异分析得到结果之后第一步不是急着看差异基因列表而是先确认样本整体关系是否符合预期。最常用的工具是PCA图。这里有个细节不能直接拿count矩阵做PCA因为原始count数据的方差跟表达量相关低表达基因的高太多噪音。需要先对count数据做variance stabilizing transformation或regularized log transformation让数据的方差更稳定再做PCA。vsd - vst(dds, blind TRUE) plotPCA(vsd, intgroup c(condition, batch))PCA图怎么看最核心一点是看生物学重复是否聚在一起看不同分组是否能在PC1或PC2方向明显分开再看图上的分组模式是否能解释。如果对照组和处理组分得清清楚楚那说明处理效应非常强如果样本不是按实验分组聚而是按批次分成两团那就是批次效应在捣乱。PCA图虽然只是无监督降维但它能把很多隐藏问题暴露出来是结果验证的第一道关口。4.2 MA图先看分布再看显著性MA图把每个基因的平均表达量baseMean放在x轴把log2FoldChange放在y轴每个点是一个基因。这个图的理想状态是整体呈现一个对称的“箭头”形状低表达区域的点在0附近也有一定散布高表达区域的点更加紧密。如果你用收缩后的log2FoldChange画MA图会发现低表达基因的y轴被明显收拢整体更干净。plotMA(res_shrink, ylim c(-4, 4))当你看到MA图上有一群高表达基因且log2FoldChange很大时要警惕是否真的生物学差异还是某个样本出了问题——比如一个样本在另一个流池发生污染这种异常会在高表达基因上体现得特别明显。4.3 火山图与差异基因标记火山图本质上是把P值和效应量画在一张图上x轴是log2FoldChangey轴是-log10(padj)。x轴越靠两边说明差异越大y轴越向上说明越显著。我一般用ggplot2自己画library(ggplot2) res_df - as.data.frame(res_shrink) res_df$gene - rownames(res_df) res_df$sig - ifelse(res_df$padj 0.05 abs(res_df$log2FoldChange) 1, sig, ns) ggplot(res_df, aes(x log2FoldChange, y -log10(padj), color sig)) geom_point(size 0.5) scale_color_manual(values c(ns grey70, sig firebrick)) theme_minimal()火山图的意义不只是好看而是让你快速判断差异基因的分布是否合理。正常情况下显著基因应该呈V字形分布在两侧而且最好比较对称。如果所有显著基因全在一侧另一侧几乎没有要思考是生物学本质还是某个处理方向定义反了。另外如果显著基因数量远超预期比如上万个基本可以判断数据有问题常见原因包括样本污染、批次未去除干净、设计公式错误导致过度拟合等。PCA、MA、火山图这三个图一起看基本能把绝大多数数据质量问题暴露出来。5. 实证分享三个月来踩过的坑和排查链路5.1 坑一某个样本比对率从85%掉到31%这个坑花了几乎一整天排查记录下完整链路供参考。现象六个样本里有一个样本在STAR比对时uniquely mapped比例只有31%其他样本都在82%-86%之间。如果用这个样本硬往下走差异分析的结果会很奇怪但表面上不会报错。排查第一步我重新跑了一遍fastp看这个样本的质控报告。果然它的接头残留率异常高reads末尾大量adaptor序列没有去除干净。重启了trim把参数调得更严格重新跑STAR比对率恢复到了80%以上。第二步是检查参考基因组版本发现这个样本本身并没有问题问题就是当时建库后原始数据里接头比例本来就偏高质控这一步没有仔细看报告就放过了。这个坑给的经验是比对率不仅仅是一个质检数字它直接决定下游表达量估计的可信度。每次拿到新样本我都会先看一眼fastp的html报告重点看adapter残留、Q20/Q30比例和duplication rate再决定要不要跑下游。5.2 坑二PCA图里样本按测序批次聚类而不是按实验分组另一个项目里PCA图一画出来明显能看出样本不是按处理分组聚在一起的而是整齐地分成两团。对照colData才发现这两团恰好对应两次测序上机。这就是典型的批次效应。解决的思路有三个层次。最简单的做法就是在DESeq2的design公式里加入batch变量即写成~ batch condition让统计模型把批次效应吸收掉。第二种做法是limma框架里的removeBatchEffect适合你后面要做的分析不是差异检验而是聚类、热图这类需要先去除批次影响的场景。第三种是ComBat-seq它对RNA-seq计数数据的批次效应处理更强适合批次效应非常严重的情况。但这里必须提醒一声design加batch的前提是batch与condition没有完全混杂。如果某个处理组全都跑在同一个批次你加batch进去会让condition估计变得极其不稳结果可能padj大量变大甚至没有显著基因。遇到这种情况就别指望统计校正兜底了应该回头梳理实验设计补充样本或者调换批次。5.3 坑三padj整列全是1或接近1这个情况我在一个样本量很小的测试数据上碰到过只有两个重复差异结果出来后padj几乎全部是1summary显示差异基因数量为0。当时第一反应是DESeq2的离散度估计出了问题。样本量太少时每个基因的离散度无法被准确估计模型会给出极大的方差自然一个都不显著。另一个可能原因是DESeq2的independent filtering机制把太多低count的基因过滤掉了剩余的基因检验能力又不足。解决办法主要有两个方向一是增加生物学重复数量这是最根本的二是在构建DESeqDataSetFromMatrix之前做预过滤把低表达基因先删掉一部分减少多重检验的次数给真正有信号的基因更多机会。keep - rowSums(counts(dds)) 10 dds - dds[keep, ]5.4 坑四count矩阵里大量基因是0该不该删拿到count矩阵先别急着跑数一数有多少基因在所有样本里都是0。如果这种基因占了很大的比例一方面会增加多重检验的负担另一方面会让离群样本对结果的影响更大。我会先按至少在最小样本组中有部分样本count大于一定阈值这个思路做预过滤比如保留至少在3个样本里count大于10的基因。这不会让你的结果失真反而会让差异检验更稳健。现象可能原因处理建议单个样本比对率骤降接头残留、建库污染重跑trim查看fastp报告如有污染考虑剔除PCA按批次聚类技术性批次效应design加batch、removeBatchEffect或ComBat-seqpadj大面积接近1重复数太少、过滤过度增加重复数、预过滤降低多重检验负担count矩阵大量0低表达基因、注释不全预过滤保留有表达信号的基因6. 从差异表到下一步富集分析前必须做的一件事6.1 导出规范化的差异基因表差异分析跑完筛选出显著基因后许多人会直接拿着基因Symbol列表去做富集分析。这里我要强调一个容易出问题的点先把你用的基因ID类型搞清楚。featureCounts输出的是GTF里gene_id那一列你要是之前用的是Ensembl的GTF那这一列基本就是ENSG开头的Ensembl ID。很多注释工具、数据库查询接口默认认的是基因Symbol或Entrez ID不转换就扔进去结果会有大量基因匹配不上。导出结果表的时候我会把Ensembl ID保留同时并行做一次ID转换用AnnotationHub或biomaRt拉一下对应关系write.csv(as.data.frame(res_shrink), DESeq2_results.csv, quote FALSE) library(AnnotationDbi) library(org.Hs.eg.db) ens2symbol - AnnotationDbi::select(org.Hs.eg.db, keys rownames(res_df), keytype ENSEMBL, columns c(SYMBOL,ENTREZID))转换之后检查一下匹配率低于90%就要想想是不是参考基因组版本跟注释库的温度对错了。6.2 进入富集分析前的两个准备动作第一个动作是把背景基因集准备好。很多初学者直接拿全部差异基因去做超几何检验这在逻辑上是错的。富集分析的基础是比较你的差异基因列表里某个通路的基因比例与背景基因集里该通路的比例是否有显著差异。背景基因集应该是所有在检测范围内、且被有效表达检测到的基因而不是全基因组所有基因。用DESeq2跑完背景集可以用所有baseMean大于0的基因来定义这样富集结果的可靠度会高很多。第二个动作是明确你做的是GO还是KEGG。GO富集用的是基因本体注释分生物学过程、细胞组分、分子功能三大类结果常出现冗余需要REVIGO之类的工具做去冗余KEGG富集侧重代谢和信号通路结果更直观但注释覆盖范围比GO窄。具体用哪个取决于你的生物学问题不要两个都订阅完直接跑结果里会出现大量和你研究无关的冗余条目。我还习惯在富集分析之前用差异基因画一张简单热图选top差异基因或者通路代表性基因直观看看它们在样本间的表达模式是否和分组一致。这一步不复杂但能在进入富集之前最后再筛掉一批假阳性基因——如果某个“显著差异基因”在热图上的表达模式完全是乱的那它在富集分析里贡献的分数也不可信。下一步我准备把基因ID转成entrez格式之后直接跑clusterProfiler到时候再接着记26.3.3的笔记。如果你也在跑这条流程建议先把本章里的样本信息表和因子水平两步检查一遍这两处看起来不起眼但在结果里放大得最明显。